Context Navigation

MultScatter_MonteCarlo_2D.ipf @ 622

Last change on this file since 622 was 622, checked in by srkline, 13 years ago
Updated MonteCarlo? to use a maximum of 4 processors. If XOP is not present (with 4 functions), the Igor code reverts to a single thread, even reverting to ipf code if necessary.
File size: 63.0 KB

Line
1	#pragma rtGlobals=1 // Use modern global access method.
2	#pragma IgorVersion=6.1
3
4	//
5	// Monte Carlo simulator for SASCALC
6	// October 2008 SRK
7	//
8	// This code simulates the scattering for a selected model based on the instrument configuration
9	// This code is based directly on John Barker's code, which includes multiple scattering effects.
10	// A lot of the setup, wave creation, and post-calculations are done in SASCALC->ReCalculateInten()
11	//
12	//
13	// at the end of the procedure fils is a very simple example of scripting for unattended simulations
14	// - not for the casual user at all.
15
16
17
18	// the RNG issue is really not worth the effort. multiple copies with different RNG is as good as I need. Plus,
19	// whatever XOP crashing was happining during threading is really unlikely to be from the RNG
20	//
21	// *** look into erand48() as the (pseudo) random number generator (it's a standard c-lib function, at least on unix)
22	// and is apparantly thread safe. drand48() returns values [0.0,1.0)
23	//http://qnxcs.unomaha.edu/help/product/neutrino/lib_ref/e/erand48.html
24	//
25
26
27	// - Why am I off by a factor of 2.7 - 3.7 (MC too high) relative to real data?
28	// I need to include efficiency (70%?) - do I knock these off before the simulation or do I
29	// really simulate that some fraction of neutrons on the detector don't actually get counted?
30	// Is the flux estimate up-to-date? !! Flux estimates at NG3 are out-of-date....
31	// - my simulated transmission is larger than what is measured, even after correcting for the quartz cell.
32	// Why? Do I need to include absorption? Just inherent problems with incoherent cross sections?
33
34	// - Most importantly, this needs to be checked for correctness of the MC simulation
35	// X how can I get the "data" on absolute scale? This would be a great comparison vs. the ideal model calculation
36	// X why does my integrated tau not match up with John's analytical calculations? where are the assumptions?
37	// X get rid of all small angle assumptions - to make sure that the calculation is correct at all angles
38
39	//
40	// X at the larger angles, is the "flat" detector being properly accounted for - in terms of
41	// the solid angle and how many counts fall in that pixel. Am I implicitly defining a spherical detector
42	// so that what I see is already "corrected"?
43	// X the MC will, of course benefit greatly from being XOPized. Maybe think about parallel implementation
44	// by allowing the data arrays to accumulate. First pass at the XOP is done. Not pretty, not the speediest (5.8x)
45	// but it is functional. Add spinCursor() for long calculations. See WaveAccess XOP example.
46	// X the background parameter for the model MUST be zero, or the integration for scattering
47	// power will be incorrect. (now the LAST point in a copy of the coef wave is set to zero, only for the rad_dev calculation
48	// X fully use the SASCALC input, most importantly, flux on sample.
49	// X if no MC desired, still use the selected model
50	// X better display of MC results on panel
51	// X settings for "count for X seconds" or "how long to 1E6 cts on detector" (but 1E6 is typically too many counts...)
52	// X warn of projected simulation time
53	// - add quartz window scattering to the simulation somehow
54	// -?- do smeared models make any sense?? Yes, John agrees that they do, and may be used in a more realistic simulation
55	// -?- but the random deviate can't be properly calculated...
56	// - make sure that the ratio of scattering coherent/incoherent is properly adjusted for the sample composition
57	// or the volume fraction of solvent.
58	//
59	// X add to the results the fraction of coherently scattered neutrons that are singly scattered, different than
60	// the overall fraction of singly scattered, and maybe more important to know.
61	//
62	// X change the fraction reaching the detector to exclude those that don't interact. These transmitted neutrons
63	// aren't counted. Is the # that interact a better number?
64	//
65	// - do we want to NOT offset the data by a multiplicative factor as it is "frozen" , so that the
66	// effects on the absolute scale can be seen?
67	//
68	// X why is "pure" incoherent scattering giving me a q^-1 slope, even with the detector all the way back?
69	// -NO- can I speed up by assuming everything interacts? This would compromise the ability to calculate multiple scattering
70	// X ask John how to verify what is going on
71	// - a number of models are now found to be ill-behaved when q=1e-10. Then the random deviate calculation blows up.
72	// a warning has been added - but some models need a proper limiting value, and some (power-law) are simply unuseable
73	// unless something else can be done. Using a log-spacing of points doesn't seem to help, and it introduces a lot of
74	// other problems. Not the way to go.
75	// - if the MC gags on a simulation, it often gets "stuck" and can't do the normal calculation from the model, which it
76	// should always default to...
77	//
78	//
79
80
81
82
83
84	// setting the flag to 1 == 2D simulation
85	// any other value for flag == 1D simulation
86	//
87	// must remember to close/reopen the simulation control panel to get the correct panel
88	//
89	Function Set_2DMonteCarlo_Flag(value)
90	Variable value
91
92	NVAR flag=root:Packages:NIST:SAS:gDoMonteCarlo
93	flag=value
94	return(0)
95	end
96
97	// threaded call to the main function, adds up the individual runs, and returns what is to be displayed
98	// results is calculated and sent back for display
99	Function Monte_SANS_Threaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
100	WAVE inputWave,ran_dev,nt,j1,j2,nn,linear_data,results
101
102	//initialize ran1 in the XOP by passing a negative integer
103	// does nothing in the Igor code
104	Duplicate/O results retWave
105
106	Variable NNeutron=inputWave[0]
107	Variable i,nthreads= ThreadProcessorCount
108
109	// make sure that the XOP exists if we are going to thread
110	#if exists("Monte_SANSX4")
111	//OK
112	if(nthreads>4) //only support 4 processors until I can figure out how to properly thread the XOP and to loop it
113	nthreads=4
114	endif
115	#else
116	nthreads = 1
117	#endif
118
119	// nthreads = 1
120	NVAR mt=root:myGlobals:gThreadGroupID
121	mt = ThreadGroupCreate(nthreads)
122	NVAR gInitTime = root:Packages:NIST:SAS:gRanDateTime //time that SASCALC was started
123	// Print "thread group ID = ",mt
124
125	inputWave[0] = NNeutron/nthreads //split up the number of neutrons
126
127	for(i=0;i<nthreads;i+=1)
128	Duplicate/O nt $("nt"+num2istr(i)) //new instance for each thread
129	Duplicate/O j1 $("j1"+num2istr(i))
130	Duplicate/O j2 $("j2"+num2istr(i))
131	Duplicate/O nn $("nn"+num2istr(i))
132	Duplicate/O linear_data $("linear_data"+num2istr(i))
133	Duplicate/O retWave $("retWave"+num2istr(i))
134	Duplicate/O inputWave $("inputWave"+num2istr(i))
135	Duplicate/O ran_dev $("ran_dev"+num2istr(i))
136
137	// ?? I need explicit wave references?
138	// maybe I need to have everything in separate data folders - bu tI haven't tried that. seems like a reach.
139	// more likely there is something bad going on in the XOP code.
140	if(i==0)
141	WAVE inputWave0,ran_dev0,nt0,j10,j20,nn0,linear_data0,retWave0
142	retWave0 = 0 //clear the return wave
143	retWave0[0] = -1*trunc(datetime-gInitTime) //to initialize ran3
144	ThreadStart mt,i,Monte_SANS_W1(inputWave0,ran_dev0,nt0,j10,j20,nn0,linear_data0,retWave0)
145	Print "started thread 0"
146	endif
147	if(i==1)
148	WAVE inputWave1,ran_dev1,nt1,j11,j21,nn1,linear_data1,retWave1
149	retWave1 = 0 //clear the return wave
150	retWave1[0] = -1*trunc(datetime-gInitTime-2) //to initialize ran1
151	ThreadStart mt,i,Monte_SANS_W2(inputWave1,ran_dev1,nt1,j11,j21,nn1,linear_data1,retWave1)
152	Print "started thread 1"
153	endif
154	if(i==2)
155	WAVE inputWave2,ran_dev2,nt2,j12,j22,nn2,linear_data2,retWave2
156	retWave2[0] = -1*trunc(datetime-gInitTime-3) //to initialize ran3a
157	ThreadStart mt,i,Monte_SANS_W3(inputWave2,ran_dev2,nt2,j12,j22,nn2,linear_data2,retWave2)
158	endif
159	if(i==3)
160	WAVE inputWave3,ran_dev3,nt3,j13,j23,nn3,linear_data3,retWave3
161	retWave3[0] = -1*trunc(datetime-gInitTime-4) //to initialize ran1a
162	ThreadStart mt,i,Monte_SANS_W4(inputWave3,ran_dev3,nt3,j13,j23,nn3,linear_data3,retWave3)
163	endif
164	endfor
165
166	// wait until done
167	do
168	variable tgs= ThreadGroupWait(mt,100)
169	while( tgs != 0 )
170	variable dummy= ThreadGroupRelease(mt)
171	mt=0
172	Print "done with all threads"
173
174	// calculate all of the bits for the results
175	if(nthreads == 1)
176	nt = nt0 // add up each instance
177	j1 = j10
178	j2 = j20
179	nn = nn0
180	linear_data = linear_data0
181	retWave = retWave0
182	endif
183	if(nthreads == 2)
184	nt = nt0+nt1 // add up each instance
185	j1 = j10+j11
186	j2 = j20+j21
187	nn = nn0+nn1
188	linear_data = linear_data0+linear_data1
189	retWave = retWave0+retWave1
190	endif
191	if(nthreads == 3)
192	nt = nt0+nt1+nt2 // add up each instance
193	j1 = j10+j11+j12
194	j2 = j20+j21+j22
195	nn = nn0+nn1+nn2
196	linear_data = linear_data0+linear_data1+linear_data2
197	retWave = retWave0+retWave1+retWave2
198	endif
199	if(nthreads == 4)
200	nt = nt0+nt1+nt2+nt3 // add up each instance
201	j1 = j10+j11+j12+j13
202	j2 = j20+j21+j22+j23
203	nn = nn0+nn1+nn2+nn3
204	linear_data = linear_data0+linear_data1+linear_data2+linear_data3
205	retWave = retWave0+retWave1+retWave2+retWave3
206	endif
207
208	// fill up the results wave
209	Variable xc,yc
210	xc=inputWave[3]
211	yc=inputWave[4]
212	results[0] = inputWave[9]+inputWave[10] //total XS
213	results[1] = inputWave[10] //SAS XS
214	results[2] = retWave[1] //number that interact n2
215	results[3] = retWave[2] - linear_data[xc][yc] //# reaching detector minus Q(0)
216	results[4] = retWave[3]/retWave[1] //avg# times scattered
217	results[5] = retWave[4]/retWave[7] //single coherent fraction
218	results[6] = retWave[5]/retWave[7] //multiple coherent fraction
219	results[7] = retWave[6]/retWave[1] //multiple scatter fraction
220	results[8] = (retWave[0]-retWave[1])/retWave[0] //transmitted fraction
221
222	return(0)
223	End
224
225	// worker function for threads, does nothing except switch between XOP and Igor versions
226	//
227	// uses ran3
228	ThreadSafe Function Monte_SANS_W1(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
229	WAVE inputWave,ran_dev,nt,j1,j2,nn,linear_data,results
230
231	#if exists("Monte_SANSX")
232	Monte_SANSX(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
233	#else
234	Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
235	#endif
236
237	return (0)
238	End
239
240	// worker function for threads, does nothing except switch between XOP and Igor versions
241	//
242	// uses ran1
243	ThreadSafe Function Monte_SANS_W2(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
244	WAVE inputWave,ran_dev,nt,j1,j2,nn,linear_data,results
245
246	#if exists("Monte_SANSX2")
247	Monte_SANSX2(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
248	#else
249	// Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
250	#endif
251
252	return (0)
253	End
254
255	// uses ran3a
256	ThreadSafe Function Monte_SANS_W3(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
257	WAVE inputWave,ran_dev,nt,j1,j2,nn,linear_data,results
258
259	#if exists("Monte_SANSX3")
260	Monte_SANSX3(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
261	#else
262	// Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
263	#endif
264
265	return (0)
266	End
267
268	// uses ran1a
269	ThreadSafe Function Monte_SANS_W4(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
270	WAVE inputWave,ran_dev,nt,j1,j2,nn,linear_data,results
271
272	#if exists("Monte_SANSX4")
273	Monte_SANSX4(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
274	#else
275	// Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
276	#endif
277
278	return (0)
279	End
280
281
282
283	// NON-threaded call to the main function returns what is to be displayed
284	// results is calculated and sent back for display
285	Function Monte_SANS_NotThreaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
286	WAVE inputWave,ran_dev,nt,j1,j2,nn,linear_data,results
287
288	//initialize ran1 in the XOP by passing a negative integer
289	// does nothing in the Igor code, enoise is already initialized
290	Duplicate/O results retWave
291	WAVE retWave
292	retWave[0] = -1abs(trunc(100000enoise(1)))
293
294	#if exists("Monte_SANSX")
295	Monte_SANSX(inputWave,ran_dev,nt,j1,j2,nn,linear_data,retWave)
296	#else
297	Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,retWave)
298	#endif
299
300	// fill up the results wave
301	Variable xc,yc
302	xc=inputWave[3]
303	yc=inputWave[4]
304	results[0] = inputWave[9]+inputWave[10] //total XS
305	results[1] = inputWave[10] //SAS XS
306	results[2] = retWave[1] //number that interact n2
307	results[3] = retWave[2] - linear_data[xc][yc] //# reaching detector minus Q(0)
308	results[4] = retWave[3]/retWave[1] //avg# times scattered
309	results[5] = retWave[4]/retWave[7] //single coherent fraction
310	results[6] = retWave[5]/retWave[7] //double coherent fraction
311	results[7] = retWave[6]/retWave[1] //multiple scatter fraction
312	results[8] = (retWave[0]-retWave[1])/retWave[0] //transmitted fraction
313
314	return(0)
315	End
316
317
318
319	//////////
320	// PROGRAM Monte_SANS
321	// PROGRAM simulates multiple SANS.
322	// revised 2/12/99 JGB
323	// added calculation of random deviate, and 2D 10/2008 SRK
324
325	// N1 = NUMBER OF INCIDENT NEUTRONS.
326	// N2 = NUMBER INTERACTED IN THE SAMPLE.
327	// N3 = NUMBER ABSORBED.
328	// THETA = SCATTERING ANGLE.
329
330	// IMON = 'Enter number of neutrons to use in simulation.'
331	// NUM_BINS = 'Enter number of THETA BINS TO use. (<5000).'
332	// R1 = 'Enter beam radius. (cm)'
333	// R2 = 'Enter sample radius. (cm)'
334	// thick = 'Enter sample thickness. (cm)'
335	// wavelength = 'Enter neutron wavelength. (A)'
336	// R0 = 'Enter sphere radius. (A)'
337	//
338
339	ThreadSafe Function Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,MC_linear_data,results)
340	WAVE inputWave,ran_dev,nt,j1,j2,nn,MC_linear_data,results
341
342	Variable imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,sig_sas
343	Variable NUM_BINS,N_INDEX
344	Variable RHO,SIGSAS,SIGABS_0
345	Variable ii,jj,IND,idum,INDEX,IR,NQ
346	Variable qmax,theta_max,R_DAB,R0,BOUND,I0,q0,zpow
347	Variable N1,N2,N3,DTH,zz,tt,SIG_SINGLE,xx,yy,PHI,UU,SIG
348	Variable THETA,Ran,ll,D_OMEGA,RR,Tabs,Ttot,I1_sumI
349	Variable G0,E_NT,E_NN,TRANS_th,Trans_exp,rat
350	Variable GG,GG_ED,dS_dW,ds_dw_double,ds_dw_single
351	Variable DONE,FIND_THETA,err //used as logicals
352
353	Variable Vx,Vy,Vz,Theta_z,qq
354	Variable Sig_scat,Sig_abs,Ratio,Sig_total
355	Variable isOn=0,testQ,testPhi,xPixel,yPixel
356	Variable NSingleIncoherent,NSingleCoherent,NScatterEvents,incoherentEvent,coherentEvent
357	Variable NDoubleCoherent,NMultipleScatter,countIt,detEfficiency
358	Variable NMultipleCoherent,NCoherentEvents
359
360	detEfficiency = 1.0 //70% counting efficiency = 0.7
361
362	imon = inputWave[0]
363	r1 = inputWave[1]
364	r2 = inputWave[2]
365	xCtr = inputWave[3]
366	yCtr = inputWave[4]
367	sdd = inputWave[5]
368	pixSize = inputWave[6]
369	thick = inputWave[7]
370	wavelength = inputWave[8]
371	sig_incoh = inputWave[9]
372	sig_sas = inputWave[10]
373
374	// SetRandomSeed 0.1 //to get a reproduceable sequence
375
376	//scattering power and maximum qvalue to bin
377	// zpow = .1 //scattering power, calculated below
378	qmax = 4*pi/wavelength //maximum Q to bin 1D data. (A-1) (not really used, so set to a big value)
379	sigabs_0 = 0.0 // ignore absorption cross section/wavelength [1/(cm A)]
380	N_INDEX = 50 // maximum number of scattering events per neutron
381	num_bins = 200 //number of 1-D bins (not really used)
382
383	// my additions - calculate the random deviate function as needed
384	// and calculate the scattering power from the model function (passed in as a wave)
385	//
386	Variable left = leftx(ran_dev)
387	Variable delta = deltax(ran_dev)
388
389	//c total SAS cross-section
390	// SIG_SAS = zpow/thick
391	zpow = sig_sas*thick //since I now calculate the sig_sas from the model
392	SIG_ABS = SIGABS_0 * WAVElength
393	sig_abs = 0.0 //cm-1
394	SIG_TOTAL =SIG_ABS + SIG_SAS + sig_incoh
395	// Print "The TOTAL XSECTION. (CM-1) is ",sig_total
396	// Print "The TOTAL SAS XSECTION. (CM-1) is ",sig_sas
397	// results[0] = sig_total //assign these after everything's done
398	// results[1] = sig_sas
399	// variable ratio1,ratio2
400	// ratio1 = sig_abs/sig_total
401	// ratio2 = sig_incoh/sig_total
402	// // 0->ratio1 = abs
403	// // ratio1 -> ratio2 = incoh
404	// // > ratio2 = coherent
405	RATIO = sig_incoh / SIG_TOTAL
406
407	//c assuming theta = sin(theta)...OK
408	theta_max = wavelengthqmax/(2pi)
409	//C SET Theta-STEP SIZE.
410	DTH = Theta_max/NUM_BINS
411	// Print "theta bin size = dth = ",dth
412
413	//C INITIALIZE COUNTERS.
414	N1 = 0
415	N2 = 0
416	N3 = 0
417	NSingleIncoherent = 0
418	NSingleCoherent = 0
419	NDoubleCoherent = 0
420	NMultipleScatter = 0
421	NScatterEvents = 0
422	NMultipleCoherent = 0
423	NCoherentEvents = 0
424
425	//C INITIALIZE ARRAYS.
426	j1 = 0
427	j2 = 0
428	nt = 0
429	nn=0
430
431	//C MONITOR LOOP - looping over the number of incedent neutrons
432	//note that zz, is the z-position in the sample - NOT the scattering power
433
434	// NOW, start the loop, throwing neutrons at the sample.
435	do
436	Vx = 0.0 // Initialize direction vector.
437	Vy = 0.0
438	Vz = 1.0
439
440	Theta = 0.0 // Initialize scattering angle.
441	Phi = 0.0 // Intialize azimuthal angle.
442	N1 += 1 // Increment total number neutrons counter.
443	DONE = 0 // True when neutron is scattered out of the sample.
444	INDEX = 0 // Set counter for number of scattering events.
445	zz = 0.0 // Set entering dimension of sample.
446	incoherentEvent = 0
447	coherentEvent = 0
448
449	do // Makes sure position is within circle.
450	ran = abs(enoise(1)) //[0,1]
451	xx = 2.0R1(Ran-0.5) //X beam position of neutron entering sample.
452	ran = abs(enoise(1)) //[0,1]
453	yy = 2.0R1(Ran-0.5) //Y beam position ...
454	RR = SQRT(xxxx+yyyy) //Radial position of neutron in incident beam.
455	while(rr>r1)
456
457	do //Scattering Loop, will exit when "done" == 1
458	// keep scattering multiple times until the neutron exits the sample
459	ran = abs(enoise(1)) //[0,1] RANDOM NUMBER FOR DETERMINING PATH LENGTH
460	ll = PATH_len(ran,Sig_total)
461	//Determine new scattering direction vector.
462	err = NewDirection(vx,vy,vz,Theta,Phi) //vx,vy,vz is updated, theta, phi unchanged by function
463
464	//X,Y,Z-POSITION OF SCATTERING EVENT.
465	xx += ll*vx
466	yy += ll*vy
467	zz += ll*vz
468	RR = sqrt(xxxx+yyyy) //radial position of scattering event.
469
470	//Check whether interaction occurred within sample volume.
471	IF (((zz > 0.0) && (zz < THICK)) && (rr < r2))
472	//NEUTRON INTERACTED.
473	ran = abs(enoise(1)) //[0,1]
474
475	// if(ran<ratio1)
476	// //absorption event
477	// n3 +=1
478	// done=1
479	// else
480
481	INDEX += 1 //Increment counter of scattering events.
482	IF(INDEX == 1)
483	N2 += 1 //Increment # of scat. neutrons
484	Endif
485	//Split neutron interactions into scattering and absorption events
486	// IF(ran > (ratio1+ratio2) ) //C NEUTRON SCATTERED coherently
487	IF(ran > ratio) //C NEUTRON SCATTERED coherently
488	coherentEvent += 1
489	FIND_THETA = 0 //false
490	DO
491	//ran = abs(enoise(1)) //[0,1]
492	//theta = Scat_angle(Ran,R_DAB,wavelength) // CHOOSE DAB ANGLE -- this is 2Theta
493	//Q0 = 2PITHETA/WAVElength // John chose theta, calculated Q
494
495	// pick a q-value from the deviate function
496	// pnt2x truncates the point to an integer before returning the x
497	// so get it from the wave scaling instead
498	Q0 =left + binarysearchinterp(ran_dev,abs(enoise(1)))*delta
499	theta = Q0/2/Pi*wavelength //SAS approximation. 1% error at theta=30 deg (theta/2=15deg)
500
501	//Print "q0, theta = ",q0,theta
502
503	FIND_THETA = 1 //always accept
504
505	while(!find_theta)
506	ran = abs(enoise(1)) //[0,1]
507	PHI = 2.0PIRan //Chooses azimuthal scattering angle.
508	ELSE
509	//NEUTRON scattered incoherently
510	// N3 += 1
511	// DONE = 1
512	// phi and theta are random over the entire sphere of scattering
513	// !can't just choose random theta and phi, won't be random over sphere solid angle
514	incoherentEvent += 1
515
516	ran = abs(enoise(1)) //[0,1]
517	theta = acos(2*ran-1)
518
519	ran = abs(enoise(1)) //[0,1]
520	PHI = 2.0PIRan //Chooses azimuthal scattering angle.
521	ENDIF //(ran > ratio)
522	// endif // event was absorption
523	ELSE
524	//NEUTRON ESCAPES FROM SAMPLE -- bin it somewhere
525	DONE = 1 //done = true, will exit from loop
526
527	// countIt = 1
528	// if(abs(enoise(1)) > detEfficiency) //efficiency of 70% wired @top
529	// countIt = 0 //detector does not register
530	// endif
531
532	//Increment #scattering events array
533	If (index <= N_Index)
534	NN[INDEX] += 1
535	Endif
536
537	if(index != 0) //the neutron interacted at least once, figure out where it ends up
538
539	Theta_z = acos(Vz) // Angle WITH respect to z axis.
540	testQ = 2pisin(theta_z)/wavelength
541
542	// pick a random phi angle, and see if it lands on the detector
543	// since the scattering is isotropic, I can safely pick a new, random value
544	// this would not be true if simulating anisotropic scattering.
545	//testPhi = abs(enoise(1))2Pi
546	testPhi = MC_FindPhi(Vx,Vy) //use the exiting phi value as defined by Vx and Vy
547
548	// is it on the detector?
549	FindPixel(testQ,testPhi,wavelength,sdd,pixSize,xCtr,yCtr,xPixel,yPixel)
550
551	if(xPixel != -1 && yPixel != -1)
552	//if(index==1) // only the single scattering events
553	MC_linear_data[xPixel][yPixel] += 1 //this is the total scattering, including multiple scattering
554	//endif
555	isOn += 1 // neutron that lands on detector
556	endif
557
558	If(theta_z < theta_max)
559	//Choose index for scattering angle array.
560	//IND = NINT(THETA_z/DTH + 0.4999999)
561	ind = round(THETA_z/DTH + 0.4999999) //round is eqivalent to nint()
562	NT[ind] += 1 //Increment bin for angle.
563	//Increment angle array for single scattering events.
564	IF(INDEX == 1)
565	j1[ind] += 1
566	Endif
567	//Increment angle array for double scattering events.
568	IF (INDEX == 2)
569	j2[ind] += 1
570	Endif
571	EndIf
572
573	// increment all of the counters now since done==1 here and I'm sure to exit and get another neutron
574	NScatterEvents += index //total number of scattering events
575	if(index == 1 && incoherentEvent == 1)
576	NSingleIncoherent += 1
577	endif
578	if(index == 1 && coherentEvent == 1)
579	NSingleCoherent += 1
580	endif
581	if(index == 2 && coherentEvent == 1 && incoherentEvent == 0)
582	NDoubleCoherent += 1
583	endif
584	if(index > 1)
585	NMultipleScatter += 1
586	endif
587	if(coherentEvent >= 1 && incoherentEvent == 0)
588	NCoherentEvents += 1
589	endif
590	if(coherentEvent > 1 && incoherentEvent == 0)
591	NMultipleCoherent += 1
592	endif
593
594
595	//Print "n1,index (x,y) = ",n1,index, xpixel,ypixel
596
597	else // if neutron escaped without interacting
598
599	// then it must be a transmitted neutron
600	// don't need to calculate, just increment the proper counters
601
602	MC_linear_data[xCtr+xx/pixsize][yCtr+yy/pixsize] += 1
603	isOn += 1
604	nt[0] += 1
605
606	endif //if interacted
607	ENDIF
608	while (!done)
609	while(n1 < imon)
610
611	// Print "Monte Carlo Done"
612	results[0] = n1
613	results[1] = n2
614	results[2] = isOn
615	results[3] = NScatterEvents //sum of # of times that neutrons scattered (coh+incoh)
616	results[4] = NSingleCoherent //# of events that are single, coherent
617	results[5] = NMultipleCoherent //# of scattered neutrons that are coherently scattered more than once
618	results[6] = NMultipleScatter //# of scattered neutrons that are scattered more than once (coh and/or incoh)
619	results[7] = NCoherentEvents //# of scattered neutrons that are scattered coherently one or more times
620
621	// Print "# absorbed = ",n3
622
623	// trans_th = exp(-sig_total*thick)
624	// TRANS_exp = (N1-N2) / N1 // Transmission
625	// dsigma/domega assuming isotropic scattering, with no absorption.
626	// Print "trans_exp = ",trans_exp
627	// Print "total # of neutrons reaching 2D detector",isOn
628	// Print "fraction of incident neutrons reaching detector ",isOn/iMon
629
630	// Print "Total number of neutrons = ",N1
631	// Print "Total number of neutrons that interact = ",N2
632	// Print "Fraction of singly scattered neutrons = ",sum(j1,-inf,inf)/N2
633	// results[2] = N2 //number that scatter
634	// results[3] = isOn - MC_linear_data[xCtr][yCtr] //# scattered reaching detector minus zero angle
635
636
637	// Tabs = (N1-N3)/N1
638	// Ttot = (N1-N2)/N1
639	// I1_sumI = NN[0]/(N2-N3)
640	// Print "Tabs = ",Tabs
641	// Print "Transmitted neutrons = ",Ttot
642	// results[8] = Ttot
643	// Print "I1 / all I1 = ", I1_sumI
644
645	End
646	//////// end of main function for calculating multiple scattering
647
648
649	// returns the random deviate as a wave
650	// and the total SAS cross-section [1/cm] sig_sas
651	Function CalculateRandomDeviate(func,coef,lam,outWave,SASxs)
652	FUNCREF SANSModelAAO_MCproto func
653	WAVE coef
654	Variable lam
655	String outWave
656	Variable &SASxs
657
658	Variable nPts_ran=10000,qu
659	qu = 4*pi/lam
660
661	// hard-wired into the Simulation directory rather than the SAS folder.
662	// plotting resolution-smeared models won't work any other way
663	Make/O/N=(nPts_ran)/D root:Simulation:Gq,root:Simulation:xw // if these waves are 1000 pts, the results are "pixelated"
664	WAVE Gq = root:Simulation:gQ
665	WAVE xw = root:Simulation:xw
666	SetScale/I x (0+1e-4),qu*(1-1e-10),"", Gq,xw //don't start at zero or run up all the way to qu to avoid numerical errors
667
668	/// if all of the coefficients are well-behaved, then the last point is the background
669	// and I can set it to zero here (only for the calculation)
670	Duplicate/O coef,tmp_coef
671	Variable num=numpnts(coef)
672	tmp_coef[num-1] = 0
673
674	xw=x //for the AAO
675	func(tmp_coef,Gq,xw) //call as AAO
676
677	// Gq = x*Gq // SAS approximation
678	Gq = Gqsin(2asin(x/qu))/sqrt(1-(x/qu)) // exact
679	//
680	//
681	Integrate/METH=1 Gq/D=Gq_INT
682
683	// SASxs = lamlam/2/piGq_INT[nPts_ran-1] //if the approximation is used
684	SASxs = lam*Gq_INT[nPts_ran-1]
685
686	Gq_INT /= Gq_INT[nPts_ran-1]
687
688	Duplicate/O Gq_INT $outWave
689
690	return(0)
691	End
692
693	// returns the random deviate as a wave
694	// and the total SAS cross-section [1/cm] sig_sas
695	//
696	// uses a log spacing of x for better coverage
697	// downside is that it doesn't use built-in integrate, which is automatically cumulative
698	//
699	// --- Currently does not work - since the usage of the random deviate in the MC routine is based on the
700	// wave properties of ran_dev, that is it must have the proper scaling and be equally spaced.
701	//
702	// -- not really sure that this will solve any of the problems with some functions (notably those with power-laws)
703	// giving unreasonably large SAS cross sections. (>>10)
704	//
705	Function CalculateRandomDeviate_log(func,coef,lam,outWave,SASxs)
706	FUNCREF SANSModelAAO_MCproto func
707	WAVE coef
708	Variable lam
709	String outWave
710	Variable &SASxs
711
712	Variable nPts_ran=1000,qu,qmin,ii
713	qmin=1e-5
714	qu = 4*pi/lam
715
716	// hard-wired into the Simulation directory rather than the SAS folder.
717	// plotting resolution-smeared models won't work any other way
718	Make/O/N=(nPts_ran)/D root:Simulation:Gq,root:Simulation:xw // if these waves are 1000 pts, the results are "pixelated"
719	WAVE Gq = root:Simulation:gQ
720	WAVE xw = root:Simulation:xw
721	// SetScale/I x (0+1e-4),qu*(1-1e-10),"", Gq,xw //don't start at zero or run up all the way to qu to avoid numerical errors
722	xw = alog(log(qmin) + x*((log(qu)-log(qmin))/nPts_ran))
723
724	/// if all of the coefficients are well-behaved, then the last point is the background
725	// and I can set it to zero here (only for the calculation)
726	Duplicate/O coef,tmp_coef
727	Variable num=numpnts(coef)
728	tmp_coef[num-1] = 0
729
730	func(tmp_coef,Gq,xw) //call as AAO
731	Gq = Gqsin(2asin(xw/qu))/sqrt(1-(xw/qu)) // exact
732
733
734	Duplicate/O Gq Gq_INT
735	Gq_INT = 0
736	for(ii=0;ii<nPts_ran;ii+=1)
737	Gq_INT[ii] = AreaXY(xw,Gq,qmin,xw[ii])
738	endfor
739
740	SASxs = lam*Gq_INT[nPts_ran-1]
741
742	Gq_INT /= Gq_INT[nPts_ran-1]
743
744	Duplicate/O Gq_INT $outWave
745
746	return(0)
747	End
748
749	ThreadSafe Function FindPixel(testQ,testPhi,lam,sdd,pixSize,xCtr,yCtr,xPixel,yPixel)
750	Variable testQ,testPhi,lam,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel
751
752	Variable theta,dy,dx,qx,qy
753	//decompose to qx,qy
754	qx = testQ*cos(testPhi)
755	qy = testQ*sin(testPhi)
756
757	//convert qx,qy to pixel locations relative to # of pixels x, y from center
758	theta = 2asin(qylam/4/pi)
759	dy = sdd*tan(theta)
760	yPixel = round(yCtr + dy/pixSize)
761
762	theta = 2asin(qxlam/4/pi)
763	dx = sdd*tan(theta)
764	xPixel = round(xCtr + dx/pixSize)
765
766	NVAR pixelsX = root:myGlobals:gNPixelsX
767	NVAR pixelsY = root:myGlobals:gNPixelsY
768
769	//if on detector, return xPix and yPix values, otherwise -1
770	if(yPixel > pixelsY \|\| yPixel < 0)
771	yPixel = -1
772	endif
773	if(xPixel > pixelsX \|\| xPixel < 0)
774	xPixel = -1
775	endif
776
777	return(0)
778	End
779
780	Function MC_CheckFunctionAndCoef(funcStr,coefStr)
781	String funcStr,coefStr
782
783	SVAR/Z listStr=root:Packages:NIST:coefKWStr
784	if(SVAR_Exists(listStr) == 1)
785	String properCoefStr = StringByKey(funcStr, listStr ,"=",";",0)
786	if(cmpstr("",properCoefStr)==0)
787	return(0) //false, no match found, so properCoefStr is returned null
788	endif
789	if(cmpstr(coefStr,properCoefStr)==0)
790	return(1) //true, the coef is the correct match
791	endif
792	endif
793	return(0) //false, wrong coef
794	End
795
796	Function/S MC_getFunctionCoef(funcStr)
797	String funcStr
798
799	SVAR/Z listStr=root:Packages:NIST:coefKWStr
800	String coefStr=""
801	if(SVAR_Exists(listStr) == 1)
802	coefStr = StringByKey(funcStr, listStr ,"=",";",0)
803	endif
804	return(coefStr)
805	End
806
807	Function SANSModelAAO_MCproto(w,yw,xw)
808	Wave w,yw,xw
809
810	Print "in SANSModelAAO_MCproto function"
811	return(1)
812	end
813
814	Function/S MC_FunctionPopupList()
815	String list,tmp
816	list = User_FunctionPopupList()
817
818	//simplify the display, forcing smeared calculations behind the scenes
819	tmp = FunctionList("Smear*",";","NPARAMS:1") //smeared dependency calculations
820	list = RemoveFromList(tmp, list,";")
821
822
823	if(strlen(list)==0)
824	list = "No functions plotted"
825	endif
826
827	list = SortList(list)
828
829	return(list)
830	End
831
832
833	//Function Scat_Angle(Ran,R_DAB,wavelength)
834	// Variable Ran,r_dab,wavelength
835	//
836	// Variable qq,arg,theta
837	// qq = 1. / ( R_DAB*sqrt(1.0/Ran - 1.0) )
838	// arg = qqwavelength/(4pi)
839	// If (arg < 1.0)
840	// theta = 2.*asin(arg)
841	// else
842	// theta = pi
843	// endif
844	// Return (theta)
845	//End
846
847	//calculates new direction (xyz) from an old direction
848	//theta and phi don't change
849	ThreadSafe Function NewDirection(vx,vy,vz,theta,phi)
850	Variable &vx,&vy,&vz
851	Variable theta,phi
852
853	Variable err=0,vx0,vy0,vz0
854	Variable nx,ny,mag_xy,tx,ty,tz
855
856	//store old direction vector
857	vx0 = vx
858	vy0 = vy
859	vz0 = vz
860
861	mag_xy = sqrt(vx0vx0 + vy0vy0)
862	if(mag_xy < 1e-12)
863	//old vector lies along beam direction
864	nx = 0
865	ny = 1
866	tx = 1
867	ty = 0
868	tz = 0
869	else
870	Nx = -Vy0 / Mag_XY
871	Ny = Vx0 / Mag_XY
872	Tx = -Vz0*Vx0 / Mag_XY
873	Ty = -Vz0*Vy0 / Mag_XY
874	Tz = Mag_XY
875	endif
876
877	//new scattered direction vector
878	Vx = cos(phi)sin(theta)Tx + sin(phi)sin(theta)Nx + cos(theta)*Vx0
879	Vy = cos(phi)sin(theta)Ty + sin(phi)sin(theta)Ny + cos(theta)*Vy0
880	Vz = cos(phi)sin(theta)Tz + cos(theta)*Vz0
881
882	Return(err)
883	End
884
885	ThreadSafe Function path_len(aval,sig_tot)
886	Variable aval,sig_tot
887
888	Variable retval
889
890	retval = -1*ln(1-aval)/sig_tot
891
892	return(retval)
893	End
894
895	// globals are initialized in SASCALC.ipf
896	// coordinates if I make this part of the panel - but this breaks other things...
897	//
898	//Proc MC_SASCALC()
899	//// PauseUpdate; Silent 1 // building window...
900	//
901	//// NewPanel /W=(92,556,390,1028)/K=1 as "SANS Simulator"
902	// SetVariable MC_setvar0,pos={491,73},size={144,15},bodyWidth=80,title="# of neutrons"
903	// SetVariable MC_setvar0,format="%5.4g"
904	// SetVariable MC_setvar0,limits={-inf,inf,100},value= root:Packages:NIST:SAS:gImon
905	// SetVariable MC_setvar0_1,pos={491,119},size={131,15},bodyWidth=60,title="Thickness (cm)"
906	// SetVariable MC_setvar0_1,limits={-inf,inf,0.1},value= root:Packages:NIST:SAS:gThick
907	// SetVariable MC_setvar0_2,pos={491,96},size={149,15},bodyWidth=60,title="Incoherent XS (cm)"
908	// SetVariable MC_setvar0_2,limits={-inf,inf,0.1},value= root:Packages:NIST:SAS:gSig_incoh
909	// SetVariable MC_setvar0_3,pos={491,142},size={150,15},bodyWidth=60,title="Sample Radius (cm)"
910	// SetVariable MC_setvar0_3,limits={-inf,inf,0.1},value= root:Packages:NIST:SAS:gR2
911	// PopupMenu MC_popup0,pos={476,13},size={165,20},proc=MC_ModelPopMenuProc,title="Model Function"
912	// PopupMenu MC_popup0,mode=1,value= #"MC_FunctionPopupList()"
913	// Button MC_button0,pos={480,181},size={130,20},proc=MC_DoItButtonProc,title="Do MC Simulation"
914	// Button MC_button1,pos={644,181},size={80,20},proc=MC_Display2DButtonProc,title="Show 2D"
915	// SetVariable setvar0_3,pos={568,484},size={50,20},disable=1
916	// GroupBox group0,pos={478,42},size={267,130},title="Monte Carlo"
917	// SetVariable cntVar,pos={653,73},size={80,15},proc=CountTimeSetVarProc,title="time(s)"
918	// SetVariable cntVar,format="%d"
919	// SetVariable cntVar,limits={1,10,1},value= root:Packages:NIST:SAS:gCntTime
920	//
921	// String fldrSav0= GetDataFolder(1)
922	// SetDataFolder root:Packages:NIST:SAS:
923	// Edit/W=(476,217,746,450)/HOST=# results_desc,results
924	// ModifyTable format(Point)=1,width(Point)=0,width(results_desc)=150
925	// SetDataFolder fldrSav0
926	// RenameWindow #,T_results
927	// SetActiveSubwindow ##
928	//EndMacro
929
930	// as a stand-alone panel, extra control bar (right) and subwindow implementations don't work right
931	// for various reasons...
932	Proc MC_SASCALC()
933
934	// when opening the panel, set the raw counts check to 1
935	root:Packages:NIST:SAS:gRawCounts = 1
936
937	PauseUpdate; Silent 1 // building window...
938	NewPanel /W=(92,556,713,818)/K=1 as "SANS Simulator"
939	DoWindow/C MC_SASCALC
940
941	SetVariable MC_setvar0,pos={28,73},size={144,15},bodyWidth=80,title="# of neutrons"
942	SetVariable MC_setvar0,format="%5.4g"
943	SetVariable MC_setvar0,limits={0,inf,100},value= root:Packages:NIST:SAS:gImon
944	SetVariable MC_setvar0_1,pos={28,119},size={131,15},bodyWidth=60,title="Thickness (cm)"
945	SetVariable MC_setvar0_1,limits={0,inf,0.1},value= root:Packages:NIST:SAS:gThick
946	SetVariable MC_setvar0_2,pos={28,96},size={149,15},bodyWidth=60,title="Incoherent XS (1/cm)"
947	SetVariable MC_setvar0_2,limits={0,inf,0.1},value= root:Packages:NIST:SAS:gSig_incoh
948	SetVariable MC_setvar0_3,pos={28,142},size={150,15},bodyWidth=60,title="Sample Radius (cm)"
949	SetVariable MC_setvar0_3,limits={-inf,inf,0.1},value= root:Packages:NIST:SAS:gR2
950	PopupMenu MC_popup0,pos={13,13},size={165,20},proc=MC_ModelPopMenuProc,title="Model Function"
951	PopupMenu MC_popup0,mode=1,value= #"MC_FunctionPopupList()"
952	Button MC_button0,pos={17,181},size={130,20},proc=MC_DoItButtonProc,title="Do MC Simulation"
953	Button MC_button0,fColor=(3,52428,1)
954	Button MC_button1,pos={17,208},size={80,20},proc=MC_Display2DButtonProc,title="Show 2D"
955	Button MC_button3,pos={182,94},size={25,20},proc=showIncohXSHelp,title="?"
956	SetVariable setvar0_3,pos={105,484},size={50,20},disable=1
957	GroupBox group0,pos={15,42},size={267,130},title="Monte Carlo"
958	SetVariable cntVar,pos={185,73},size={90,15},proc=CountTimeSetVarProc,title="time(s)"
959	SetVariable cntVar,format="%d"
960	SetVariable cntVar,limits={1,3600,1},value= root:Packages:NIST:SAS:gCntTime
961	Button MC_button2,pos={17,234},size={100,20},proc=SaveAsVAXButtonProc,title="Save 2D VAX"
962	CheckBox check0,pos={216,180},size={68,14},title="Raw counts",variable = root:Packages:NIST:SAS:gRawCounts
963	CheckBox check0_1,pos={216,199},size={60,14},title="Yes Offset",variable= root:Packages:NIST:SAS:gDoTraceOffset
964	CheckBox check0_2,pos={216,199+19},size={60,14},title="Beam Stop in",variable= root:Packages:NIST:SAS:gBeamStopIn
965	CheckBox check0_3,pos={216,199+2*19},size={60,14},title="use XOP",variable= root:Packages:NIST:SAS:gUse_MC_XOP
966
967	String fldrSav0= GetDataFolder(1)
968	SetDataFolder root:Packages:NIST:SAS:
969	Edit/W=(344,23,606,248)/HOST=# results_desc,results
970	ModifyTable format(Point)=1,width(Point)=0,width(results_desc)=150
971	SetDataFolder fldrSav0
972	RenameWindow #,T_results
973	SetActiveSubwindow ##
974	EndMacro
975
976
977	//
978	Proc showIncohXSHelp(ctrlName): ButtonControl
979	String ctrlName
980	DisplayHelpTopic/K=1/Z "Approximate Incoherent Cross Section"
981	if(V_flag !=0)
982	DoAlert 0,"The SANS Simulation Help file could not be found"
983	endif
984	end
985
986
987	Function CountTimeSetVarProc(sva) : SetVariableControl
988	STRUCT WMSetVariableAction &sva
989
990	switch( sva.eventCode )
991	case 1: // mouse up
992	case 2: // Enter key
993	case 3: // Live update
994	Variable dval = sva.dval
995
996	// get the neutron flux, multiply, and reset the global for # neutrons
997	NVAR imon=root:Packages:NIST:SAS:gImon
998	imon = dval*beamIntensity()
999
1000	break
1001	endswitch
1002
1003	return 0
1004	End
1005
1006
1007	Function MC_ModelPopMenuProc(pa) : PopupMenuControl
1008	STRUCT WMPopupAction &pa
1009
1010	switch( pa.eventCode )
1011	case 2: // mouse up
1012	Variable popNum = pa.popNum
1013	String popStr = pa.popStr
1014	SVAR gStr = root:Packages:NIST:SAS:gFuncStr
1015	gStr = popStr
1016
1017	break
1018	endswitch
1019
1020	return 0
1021	End
1022
1023	Function MC_DoItButtonProc(ba) : ButtonControl
1024	STRUCT WMButtonAction &ba
1025
1026	switch( ba.eventCode )
1027	case 2: // mouse up
1028	// click code here
1029	NVAR doMC = root:Packages:NIST:SAS:gDoMonteCarlo
1030	doMC = 1
1031	ReCalculateInten(1)
1032	doMC = 0 //so the next time won't be MC
1033	break
1034	endswitch
1035
1036	return 0
1037	End
1038
1039
1040	Function MC_Display2DButtonProc(ba) : ButtonControl
1041	STRUCT WMButtonAction &ba
1042
1043	switch( ba.eventCode )
1044	case 2: // mouse up
1045	// click code here
1046	Execute "ChangeDisplay(\"SAS\")"
1047	break
1048	endswitch
1049
1050	return 0
1051	End
1052
1053	// after a 2d data image is averaged in the usual way, take the waves and generate a "fake" folder of the 1d
1054	// data, to appear as if it was loaded from a real data file.
1055	//
1056	// ---- use FakeUSANSDataFolder() if you want to fake a 1D USANS data set ----
1057	//
1058	Function Fake1DDataFolder(qval,aveint,sigave,sigmaQ,qbar,fSubs,dataFolder)
1059	WAVE qval,aveint,sigave,sigmaQ,qbar,fSubs
1060	String dataFolder
1061
1062	String baseStr=dataFolder
1063	if(DataFolderExists("root:"+baseStr))
1064	SetDataFolder $("root:"+baseStr)
1065	else
1066	NewDataFolder/S $("root:"+baseStr)
1067	endif
1068
1069	////overwrite the existing data, if it exists
1070	Duplicate/O qval, $(baseStr+"_q")
1071	Duplicate/O aveint, $(baseStr+"_i")
1072	Duplicate/O sigave, $(baseStr+"_s")
1073
1074
1075	// make a resolution matrix for SANS data
1076	Variable np=numpnts(qval)
1077	Make/D/O/N=(np,4) $(baseStr+"_res")
1078	Wave res=$(baseStr+"_res")
1079
1080	res[][0] = sigmaQ[p] //sigQ
1081	res[][1] = qBar[p] //qBar
1082	res[][2] = fSubS[p] //fShad
1083	res[][3] = qval[p] //Qvalues
1084
1085	// keep a copy of everything in SAS too... the smearing wrapper function looks for
1086	// data in folders based on waves it is passed - an I lose control of that
1087	Duplicate/O res, $("root:Packages:NIST:SAS:"+baseStr+"_res")
1088	Duplicate/O qval, $("root:Packages:NIST:SAS:"+baseStr+"_q")
1089	Duplicate/O aveint, $("root:Packages:NIST:SAS:"+baseStr+"_i")
1090	Duplicate/O sigave, $("root:Packages:NIST:SAS:"+baseStr+"_s")
1091
1092	//clean up
1093	SetDataFolder root:
1094
1095	End
1096
1097	// writes out a VAX binary data file
1098	// automatically generates a name
1099	// will prompt for the sample label
1100	//
1101	// currently hard-wired for SAS data folder
1102	//
1103	Function SaveAsVAXButtonProc(ctrlName,[runIndex,simLabel])
1104	String ctrlName
1105	Variable runIndex
1106	String simLabel
1107
1108
1109	// if default parameters were passed in, use them
1110	// if not, set them to "bad" values so that the user will be prompted later
1111	NVAR autoSaveIndex = root:Packages:NIST:SAS:gAutoSaveIndex
1112	SVAR autoSaveLabel = root:Packages:NIST:SAS:gAutoSaveLabel
1113
1114	// Determine if the optional parameters were supplied
1115	if( ParamIsDefault(runIndex)) //==1 if parameter was NOT specified
1116	print "runIndex not specified"
1117	autoSaveIndex=0 // 0 == bad value, test for this later
1118	else
1119	autoSaveIndex=runIndex
1120	endif
1121
1122	if( ParamIsDefault(simLabel)) //==1 if parameter was NOT specified
1123	print "simLabel not specified"
1124	autoSaveLabel="" // "" == bad value, test for this later
1125	else
1126	autoSaveLabel=simLabel
1127	endif
1128
1129	String fullpath="",destStr=""
1130	Variable refnum
1131
1132	fullpath = Write_RawData_File("SAS","",0)
1133
1134	// write out the results into a text file
1135	destStr = "root:Packages:NIST:SAS:"
1136	SetDataFolder $destStr
1137
1138	WAVE results=results
1139	WAVE/T results_desc=results_desc
1140
1141	//check each wave
1142	If(!(WaveExists(results)))
1143	Abort "results DNExist WriteVAXData()"
1144	Endif
1145	If(!(WaveExists(results_desc)))
1146	Abort "results_desc DNExist WriteVAXData()"
1147	Endif
1148
1149	NVAR actSimTime = root:Packages:NIST:SAS:g_actSimTime
1150	String str = ""
1151	sprintf str,"%30s\t\t%g seconds\r","MonteCarlo Simulation time = ",actSimTime
1152
1153	Open refNum as fullpath+".txt"
1154	wfprintf refNum, "%30s\t\t%g\r",results_desc,results
1155	FBinWrite refNum,str
1156	FStatus refNum
1157	FSetPos refNum,V_logEOF
1158	Close refNum
1159
1160	///////////////////////////////
1161
1162	// could also automatically do the average here, but probably not worth the the effort...
1163
1164	SetDataFolder root:
1165
1166	return(0)
1167	End
1168
1169	// calculates the fraction of the scattering that reaches the detector, given the random deviate function
1170	// and qmin and qmax
1171	//
1172	//
1173	// still some question of the corners and number of pixels per q-bin
1174	Function FractionReachingDetector(ran_dev,Qmin,Qmax)
1175	wave ran_dev
1176	Variable Qmin,Qmax
1177
1178	Variable r1,r2,frac
1179	r1=x2pnt(ran_dev,Qmin)
1180	r2=x2pnt(ran_dev,Qmax)
1181
1182	// no normalization needed - the full q-range is defined as [0,1]
1183	frac = ran_dev[r2] - ran_dev[r1]
1184
1185	return frac
1186	End
1187
1188
1189	/// called in SASCALC:ReCalculateInten()
1190	Function Simulate_2D_MC(funcStr,aveint,qval,sigave,sigmaq,qbar,fsubs)
1191	String funcStr
1192	WAVE aveint,qval,sigave,sigmaq,qbar,fsubs
1193
1194	NVAR doMonteCarlo = root:Packages:NIST:SAS:gDoMonteCarlo // == 1 if 2D MonteCarlo set by hidden flag
1195	WAVE rw=root:Packages:NIST:SAS:realsRead
1196
1197	// Try to nicely exit from a threading error, if possible
1198	Variable err=0
1199	if(!exists("root:myGlobals:gThreadGroupID"))
1200	Variable/G root:myGlobals:gThreadGroupID=0
1201	endif
1202	NVAR mt=root:myGlobals:gThreadGroupID
1203
1204	if(mt!=0) //there was an error with the stopping of the threads, possibly user abort
1205	err = ThreadGroupRelease(mt)
1206	Print "threading err = ",err
1207	if(err == 0)
1208	// all should be OK
1209	else
1210	return(0)
1211	endif
1212	endif
1213
1214	NVAR imon = root:Packages:NIST:SAS:gImon
1215	NVAR thick = root:Packages:NIST:SAS:gThick
1216	NVAR sig_incoh = root:Packages:NIST:SAS:gSig_incoh
1217	NVAR r2 = root:Packages:NIST:SAS:gR2
1218
1219	// do the simulation here, or not
1220	Variable r1,xCtr,yCtr,sdd,pixSize,wavelength
1221	String coefStr,abortStr,str
1222
1223	r1 = rw[24]/2/10 // sample diameter convert diam in [mm] to radius in cm
1224	xCtr = rw[16]
1225	yCtr = rw[17]
1226	sdd = rw[18]*100 //conver header of [m] to [cm]
1227	pixSize = rw[10]/10 // convert pix size in mm to cm
1228	wavelength = rw[26]
1229	coefStr = MC_getFunctionCoef(funcStr)
1230
1231	if(!MC_CheckFunctionAndCoef(funcStr,coefStr))
1232	doMonteCarlo = 0 //we're getting out now, reset the flag so we don't continually end up here
1233	Abort "The coefficients and function type do not match. Please correct the selections in the popup menus."
1234	endif
1235
1236	Variable sig_sas
1237	// FUNCREF SANSModelAAO_MCproto func=$("fSmeared"+funcStr) //a wrapper for the structure version
1238	FUNCREF SANSModelAAO_MCproto func=$(funcStr) //unsmeared
1239	WAVE results = root:Packages:NIST:SAS:results
1240	WAVE linear_data = root:Packages:NIST:SAS:linear_data
1241	WAVE data = root:Packages:NIST:SAS:data
1242
1243	results = 0
1244	linear_data = 0
1245
1246	CalculateRandomDeviate(func,$coefStr,wavelength,"root:Packages:NIST:SAS:ran_dev",SIG_SAS)
1247	if(sig_sas > 100)
1248	sprintf abortStr,"sig_sas = %g. Please check that the model coefficients have a zero background, or the low q is well-behaved.",sig_sas
1249	Abort abortStr
1250	endif
1251
1252	WAVE ran_dev=$"root:Packages:NIST:SAS:ran_dev"
1253
1254	Make/O/D/N=5000 root:Packages:NIST:SAS:nt=0,root:Packages:NIST:SAS:j1=0,root:Packages:NIST:SAS:j2=0
1255	Make/O/D/N=100 root:Packages:NIST:SAS:nn=0
1256	Make/O/D/N=11 root:Packages:NIST:SAS:inputWave=0
1257
1258	WAVE nt = root:Packages:NIST:SAS:nt
1259	WAVE j1 = root:Packages:NIST:SAS:j1
1260	WAVE j2 = root:Packages:NIST:SAS:j2
1261	WAVE nn = root:Packages:NIST:SAS:nn
1262	WAVE inputWave = root:Packages:NIST:SAS:inputWave
1263
1264	inputWave[0] = imon
1265	inputWave[1] = r1
1266	inputWave[2] = r2
1267	inputWave[3] = xCtr
1268	inputWave[4] = yCtr
1269	inputWave[5] = sdd
1270	inputWave[6] = pixSize
1271	inputWave[7] = thick
1272	inputWave[8] = wavelength
1273	inputWave[9] = sig_incoh
1274	inputWave[10] = sig_sas
1275
1276	linear_data = 0 //initialize
1277
1278	Variable t0,trans
1279
1280	// get a time estimate, and give the user a chance to exit if they're unsure.
1281	t0 = stopMStimer(-2)
1282	inputWave[0] = 1000
1283	NVAR useXOP = root:Packages:NIST:SAS:gUse_MC_XOP //if zero, will use non-threaded Igor code
1284
1285	if(useXOP)
1286	//use a single thread, otherwise time is dominated by overhead
1287	Monte_SANS_NotThreaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
1288	else
1289	Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
1290	endif
1291
1292	t0 = (stopMSTimer(-2) - t0)*1e-6
1293	t0 *= imon/1000/ThreadProcessorCount //projected time, in seconds (using threads for the calculation)
1294
1295	Print "Estimated Simulation time (s) = ",t0
1296
1297	// to correct for detector efficiency, send only the fraction of neutrons that are actually counted
1298	NVAR detectorEff = root:Packages:NIST:SAS:g_detectorEff
1299	NVAR actSimTime = root:Packages:NIST:SAS:g_actSimTime
1300	NVAR SimTimeWarn = root:Packages:NIST:SAS:g_SimTimeWarn
1301
1302	inputWave[0] = imon * detectorEff //reset number of input neutrons before full simulation
1303
1304	if(t0>SimTimeWarn)
1305	sprintf str,"The simulation will take approximately %d seconds.\r- Proceed?",t0
1306	DoAlert 1,str
1307	if(V_flag == 2)
1308	doMonteCarlo = 0
1309	reCalculateInten(1) //come back in and do the smeared calculation
1310	return(0)
1311	endif
1312	endif
1313
1314	linear_data = 0 //initialize
1315	// threading crashes!! - there must be some operation in the XOP that is not threadSafe. What, I don't know...
1316	// I think it's the ran() calls, being "non-reentrant". So the XOP now defines two separate functions, that each
1317	// use a different rng. This works. 1.75x speedup.
1318	t0 = stopMStimer(-2)
1319
1320	if(useXOP)
1321	Monte_SANS_Threaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
1322	else
1323	Monte_SANS_NotThreaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
1324	endif
1325
1326	t0 = (stopMSTimer(-2) - t0)*1e-6
1327	Printf "MC sim time = %g seconds\r",t0
1328	actSimTime = t0
1329
1330	trans = results[8] //(n1-n2)/n1
1331	if(trans == 0)
1332	trans = 1
1333	endif
1334
1335	Print "counts on detector, including transmitted = ",sum(linear_data,-inf,inf)
1336
1337	// linear_data[xCtr][yCtr] = 0 //snip out the transmitted spike
1338	// Print "counts on detector not transmitted = ",sum(linear_data,-inf,inf)
1339
1340	// or simulate a beamstop
1341	NVAR MC_BS_in = root:Packages:NIST:SAS:gBeamStopIn //if zero, beam stop is "out", as in a transmission measurement
1342
1343	Variable rad=beamstopDiam()/2 //beamstop radius in cm
1344	if(MC_BS_in)
1345	rad /= 0.5 //convert cm to pixels
1346	rad += 0. // (no - it cuts off the low Q artificially) add an extra pixel to each side to account for edge
1347	Duplicate/O linear_data,root:Packages:NIST:SAS:tmp_mask//,root:Packages:NIST:SAS:MC_linear_data
1348	WAVE tmp_mask = root:Packages:NIST:SAS:tmp_mask
1349	tmp_mask = (sqrt((p-xCtr)^2+(q-yCtr)^2) < rad) ? 0 : 1 //behind beamstop = 0, away = 1
1350
1351	linear_data *= tmp_mask
1352	endif
1353
1354	results[9] = sum(linear_data,-inf,inf)
1355	// Print "counts on detector not behind beamstop = ",results[9]
1356
1357	// convert to absolute scale
1358	Variable kappa //,beaminten = beamIntensity()
1359	// kappa = beamIntenpir1r1thick(pixSize/sdd)^2trans*(iMon/beaminten)
1360	kappa = thick(pixSize/sdd)^2trans*iMon
1361
1362	//use kappa to get back to counts => linear_data = round(linear_data*kappa)
1363	Note/K linear_data ,"KAPPA="+num2str(kappa)+";"
1364
1365	NVAR rawCts = root:Packages:NIST:SAS:gRawCounts
1366	if(!rawCts) //go ahead and do the linear scaling
1367	linear_data = linear_data / kappa
1368	linear_data /= detectorEff
1369	endif
1370	data = linear_data
1371
1372	// re-average the 2D data
1373	S_CircularAverageTo1D("SAS")
1374
1375	// put the new result into the simulation folder
1376	Fake1DDataFolder(qval,aveint,sigave,sigmaQ,qbar,fSubs,"Simulation")
1377
1378
1379	return(0)
1380	end
1381
1382	//phi is defined from +x axis, proceeding CCW around [0,2Pi]
1383	ThreadSafe Function MC_FindPhi(vx,vy)
1384	variable vx,vy
1385
1386	variable phi
1387
1388	phi = atan(vy/vx) //returns a value from -pi/2 to pi/2
1389
1390	// special cases
1391	if(vx==0 && vy > 0)
1392	return(pi/2)
1393	endif
1394	if(vx==0 && vy < 0)
1395	return(3*pi/2)
1396	endif
1397	if(vx >= 0 && vy == 0)
1398	return(0)
1399	endif
1400	if(vx < 0 && vy == 0)
1401	return(pi)
1402	endif
1403
1404
1405	if(vx > 0 && vy > 0)
1406	return(phi)
1407	endif
1408	if(vx < 0 && vy > 0)
1409	return(phi + pi)
1410	endif
1411	if(vx < 0 && vy < 0)
1412	return(phi + pi)
1413	endif
1414	if( vx > 0 && vy < 0)
1415	return(phi + 2*pi)
1416	endif
1417
1418	return(phi)
1419	end
1420
1421
1422
1423
1424
1425	Proc Sim_1D_Panel()
1426	PauseUpdate; Silent 1 // building window...
1427	NewPanel /W=(92,556,713,818)/K=1 as "1D SANS Simulator"
1428	DoWindow/C Sim_1D_Panel
1429
1430	SetVariable cntVar,pos={26,68},size={160,15},title="Counting time(s)",format="%d"
1431	SetVariable cntVar,limits={1,36000,10},value= root:Packages:NIST:SAS:gCntTime
1432	SetVariable cntVar, proc=Sim_1D_CountTimeSetVarProc
1433	SetVariable MC_setvar0_1,pos={26,91},size={160,15},title="Thickness (cm)"
1434	SetVariable MC_setvar0_1,limits={0,inf,0.1},value= root:Packages:NIST:SAS:gThick
1435	SetVariable MC_setvar0_1, proc=Sim_1D_SamThickSetVarProc
1436
1437	SetVariable MC_setvar0_3,pos={26,114},size={160,15},title="Sample Transmission"
1438	SetVariable MC_setvar0_3,limits={0,1,0.01},value= root:Packages:NIST:SAS:gSamTrans
1439	SetVariable MC_setvar0_3, proc=Sim_1D_SamTransSetVarProc
1440
1441	PopupMenu MC_popup0,pos={13,13},size={165,20},proc=Sim_1D_ModelPopMenuProc,title="Model Function"
1442	PopupMenu MC_popup0,mode=1,value= #"MC_FunctionPopupList()"
1443	Button MC_button0,pos={17,181},size={130,20},proc=Sim_1D_DoItButtonProc,title="Do 1D Simulation"
1444	Button MC_button0,fColor=(3,52428,1)
1445	Button MC_button1,pos={17,211},size={150,20},proc=Save_1DSimData,title="Save Simulated Data"
1446	GroupBox group0,pos={15,42},size={280,130},title="Sample Setup"
1447	CheckBox check0_1,pos={216,179},size={60,14},title="Yes Offset",variable= root:Packages:NIST:SAS:gDoTraceOffset
1448	CheckBox check0_2,pos={216,199},size={60,14},title="Abs scale?",variable= root:Packages:NIST:SAS:g_1D_DoABS
1449	CheckBox check0_3,pos={216,219},size={60,14},title="Noise?",variable= root:Packages:NIST:SAS:g_1D_AddNoise
1450
1451	// a box for the results
1452	GroupBox group1,pos={314,23},size={277,163},title="Simulation Results"
1453	ValDisplay valdisp0,pos={326,48},size={220,13},title="Total detector counts"
1454	ValDisplay valdisp0,limits={0,0,0},barmisc={0,1000},value= root:Packages:NIST:SAS:g_1DTotCts
1455	ValDisplay valdisp0_1,pos={326,72},size={220,13},title="Estimated count rate (1/s)"
1456	ValDisplay valdisp0_1,limits={0,0,0},barmisc={0,1000},value=root:Packages:NIST:SAS:g_1DEstDetCR
1457	ValDisplay valdisp0_2,pos={326,96},size={220,13},title="Fraction of beam scattered"
1458	ValDisplay valdisp0_2,limits={0,0,0},barmisc={0,1000},value= root:Packages:NIST:SAS:g_1DFracScatt
1459	ValDisplay valdisp0_3,pos={326,121},size={220,13},title="Estimated transmission"
1460	ValDisplay valdisp0_3,limits={0,0,0},barmisc={0,1000},value=root:Packages:NIST:SAS:g_1DEstTrans
1461	ValDisplay valdisp0_4,pos={326,145},size={220,13},title="Multiple Coherent Scattering"
1462	ValDisplay valdisp0_4,limits={0,0,0},barmisc={0,1000},value=root:Packages:NIST:SAS:g_MultScattFraction
1463	// set the flags here -- do the simulation, but not 2D
1464
1465	root:Packages:NIST:SAS:doSimulation = 1 // == 1 if 1D simulated data, 0 if other from the checkbox
1466	root:Packages:NIST:SAS:gDoMonteCarlo = 0 // == 1 if 2D MonteCarlo set by hidden flag
1467
1468
1469	EndMacro
1470
1471	Function Sim_1D_CountTimeSetVarProc(sva) : SetVariableControl
1472	STRUCT WMSetVariableAction &sva
1473
1474	switch( sva.eventCode )
1475	case 1: // mouse up
1476	case 2: // Enter key
1477	case 3: // Live update
1478	Variable dval = sva.dval
1479
1480	ReCalculateInten(1)
1481
1482	break
1483	endswitch
1484
1485	return 0
1486	End
1487
1488	Function Sim_1D_SamThickSetVarProc(sva) : SetVariableControl
1489	STRUCT WMSetVariableAction &sva
1490
1491	switch( sva.eventCode )
1492	case 1: // mouse up
1493	case 2: // Enter key
1494	case 3: // Live update
1495	Variable dval = sva.dval
1496
1497	ReCalculateInten(1)
1498
1499	break
1500	endswitch
1501
1502	return 0
1503	End
1504
1505	Function Sim_1D_SamTransSetVarProc(sva) : SetVariableControl
1506	STRUCT WMSetVariableAction &sva
1507
1508	switch( sva.eventCode )
1509	case 1: // mouse up
1510	case 2: // Enter key
1511	case 3: // Live update
1512	Variable dval = sva.dval
1513
1514	ReCalculateInten(1)
1515
1516	break
1517	endswitch
1518
1519	return 0
1520	End
1521
1522
1523	Function Sim_1D_ModelPopMenuProc(pa) : PopupMenuControl
1524	STRUCT WMPopupAction &pa
1525
1526	switch( pa.eventCode )
1527	case 2: // mouse up
1528	Variable popNum = pa.popNum
1529	String popStr = pa.popStr
1530	SVAR gStr = root:Packages:NIST:SAS:gFuncStr
1531	gStr = popStr
1532
1533	break
1534	endswitch
1535
1536	return 0
1537	End
1538
1539
1540	Function Sim_1D_DoItButtonProc(ba) : ButtonControl
1541	STRUCT WMButtonAction &ba
1542
1543	switch( ba.eventCode )
1544	case 2: // mouse up
1545
1546	ReCalculateInten(1)
1547
1548	break
1549	endswitch
1550
1551	return 0
1552	End
1553
1554
1555	//
1556	//
1557	//
1558	Function Save_1DSimData(ctrlName) : ButtonControl
1559	String ctrlName
1560
1561	String type="SAS",fullpath=""
1562	Variable dialog=1 //=1 will present dialog for name
1563
1564	String destStr=""
1565	destStr = "root:Packages:NIST:"+type
1566
1567	Variable refNum
1568	String formatStr = "%15.4g %15.4g %15.4g %15.4g %15.4g %15.4g\r\n"
1569	String fname,ave="C",hdrStr1="",hdrStr2=""
1570	Variable step=1
1571
1572	If(1)
1573	//setup a "fake protocol" wave, sice I have no idea of the current state of the data
1574	Make/O/T/N=8 root:myGlobals:Protocols:SIMProtocol
1575	Wave/T SIMProtocol = $"root:myGlobals:Protocols:SIMProtocol"
1576	String junk="**SIMULATED DATA**"
1577	//stick in the fake protocol...
1578	NVAR ctTime = root:Packages:NIST:SAS:gCntTime
1579	NVAR totalCts = root:Packages:NIST:SAS:g_1DTotCts //summed counts (simulated)
1580	NVAR detCR = root:Packages:NIST:SAS:g_1DEstDetCR // estimated detector count rate
1581	NVAR fractScat = root:Packages:NIST:SAS:g_1DFracScatt
1582	NVAR mScat = root:Packages:NIST:SAS:g_MultScattFraction
1583
1584	SIMProtocol[0] = junk
1585	SIMProtocol[1] = "\tCounting time (s) = "+num2str(ctTime)
1586	SIMProtocol[2] = "\tTotal detector counts = "+num2str(totalCts)
1587	SIMProtocol[3] = "\tDetector countrate (1/s) = "+num2str(detCR)
1588	SIMProtocol[4] = "\tFraction of beam scattered coherently = "+num2str(fractScat)
1589	SIMProtocol[5] = "\tFraction of multiple coherent scattering = "+num2str(mScat)
1590	SIMProtocol[6] = ""
1591	SIMProtocol[7] = ""
1592	//set the global
1593	String/G root:myGlobals:Protocols:gProtoStr = "SIMProtocol"
1594	Endif
1595
1596
1597	//***these waves MUST EXIST, or IGOR Pro will crash, with a type 2 error**
1598	WAVE intw=$(destStr + ":integersRead")
1599	WAVE rw=$(destStr + ":realsRead")
1600	WAVE/T textw=$(destStr + ":textRead")
1601	WAVE qvals =$(destStr + ":qval")
1602	WAVE inten=$(destStr + ":aveint")
1603	WAVE sig=$(destStr + ":sigave")
1604	WAVE qbar = $(destStr + ":QBar")
1605	WAVE sigmaq = $(destStr + ":SigmaQ")
1606	WAVE fsubs = $(destStr + ":fSubS")
1607
1608	SVAR gProtoStr = root:myGlobals:Protocols:gProtoStr
1609	Wave/T proto=$("root:myGlobals:Protocols:"+gProtoStr)
1610
1611	//check each wave
1612	If(!(WaveExists(intw)))
1613	Abort "intw DNExist Save_1DSimData()"
1614	Endif
1615	If(!(WaveExists(rw)))
1616	Abort "rw DNExist Save_1DSimData()"
1617	Endif
1618	If(!(WaveExists(textw)))
1619	Abort "textw DNExist Save_1DSimData()"
1620	Endif
1621	If(!(WaveExists(qvals)))
1622	Abort "qvals DNExist Save_1DSimData()"
1623	Endif
1624	If(!(WaveExists(inten)))
1625	Abort "inten DNExist Save_1DSimData()"
1626	Endif
1627	If(!(WaveExists(sig)))
1628	Abort "sig DNExist Save_1DSimData()"
1629	Endif
1630	If(!(WaveExists(qbar)))
1631	Abort "qbar DNExist Save_1DSimData()"
1632	Endif
1633	If(!(WaveExists(sigmaq)))
1634	Abort "sigmaq DNExist Save_1DSimData()"
1635	Endif
1636	If(!(WaveExists(fsubs)))
1637	Abort "fsubs DNExist Save_1DSimData()"
1638	Endif
1639	If(!(WaveExists(proto)))
1640	Abort "current protocol wave DNExist Save_1DSimData()"
1641	Endif
1642
1643	//strings can be too long to print-- must trim to 255 chars
1644	Variable ii,num=8
1645	Make/O/T/N=(num) tempShortProto
1646	for(ii=0;ii<num;ii+=1)
1647	tempShortProto[ii] = (proto[ii])[0,240]
1648	endfor
1649
1650	if(dialog)
1651	PathInfo/S catPathName
1652	fullPath = DoSaveFileDialog("Save data as")
1653	If(cmpstr(fullPath,"")==0)
1654	//user cancel, don't write out a file
1655	Close/A
1656	Abort "no data file was written"
1657	Endif
1658	//Print "dialog fullpath = ",fullpath
1659	Endif
1660
1661	NVAR monCt = root:Packages:NIST:SAS:gImon
1662	NVAR thick = root:Packages:NIST:SAS:gThick
1663	NVAR trans = root:Packages:NIST:SAS:gSamTrans //for 1D, default value
1664
1665
1666
1667	hdrStr1 = num2str(monCt)+" "+num2str(rw[26])+" "+num2str(rw[19])+" "+num2str(rw[18])
1668	hdrStr1 += " "+num2str(trans)+" "+num2str(thick) + ave +" "+num2str(step) + "\r\n"
1669
1670	hdrStr2 = num2str(rw[16])+" "+num2str(rw[17])+" "+num2str(rw[23])+" "+num2str(rw[24])+" "
1671	hdrStr2 += num2str(rw[25])+" "+num2str(rw[27])+" "+num2str(rw[21])+" "+"ORNL " + "\r\n"
1672
1673	//actually open the file here
1674	Open refNum as fullpath
1675
1676	//write out the standard header information
1677	fprintf refnum,"FILE: %s\t\t CREATED: %s\r\n","SIMULATED DATA",(date() +" "+ time())
1678	fprintf refnum,"LABEL: %s\r\n","SIMULATED DATA"
1679	fprintf refnum,"MON CNT LAMBDA DET ANG DET DIST TRANS THICK AVE STEP\r\n"
1680	fprintf refnum,hdrStr1
1681	fprintf refnum,"BCENT(X,Y) A1(mm) A2(mm) A1A2DIST(m) DL/L BSTOP(mm) DET_TYP \r\n"
1682	fprintf refnum,hdrStr2
1683	// fprintf refnum,headerFormat,rw[0],rw[26],rw[19],rw[18],rw[4],rw[5],ave,step
1684
1685	//insert protocol information here
1686	//-1 list of sample files
1687	//0 - bkg
1688	//1 - emp
1689	//2 - div
1690	//3 - mask
1691	//4 - abs params c2-c5
1692	//5 - average params
1693	fprintf refnum, "SAM: %s\r\n",tempShortProto[0]
1694	fprintf refnum, "BGD: %s\r\n",tempShortProto[1]
1695	fprintf refnum, "EMP: %s\r\n",tempShortProto[2]
1696	fprintf refnum, "DIV: %s\r\n",tempShortProto[3]
1697	fprintf refnum, "MASK: %s\r\n",tempShortProto[4]
1698	fprintf refnum, "ABS: %s\r\n",tempShortProto[5]
1699	fprintf refnum, "Average Choices: %s\r\n",tempShortProto[6]
1700
1701	//write out the data columns
1702	fprintf refnum,"The 6 columns are \| Q (1/A) \| I(Q) (1/cm) \| std. dev. I(Q) (1/cm) \| sigmaQ \| meanQ \| ShadowFactor\|\r\n"
1703	wfprintf refnum, formatStr, qvals,inten,sig,sigmaq,qbar,fsubs
1704
1705	Close refnum
1706
1707	SetDataFolder root: //(redundant)
1708
1709	//write confirmation of write operation to history area
1710	Print "Averaged File written: ", GetFileNameFromPathNoSemi(fullPath)
1711	KillWaves/Z tempShortProto
1712
1713	//clear the stuff that was created for case of saving files
1714	If(1)
1715	Killwaves/Z root:myGlobals:Protocols:SIMProtocol
1716	String/G root:myGlobals:Protocols:gProtoStr = ""
1717	Endif
1718
1719
1720	return(0)
1721
1722	End
1723
1724
1725	/// called in SASCALC:ReCalculateInten()
1726	Function Simulate_1D(funcStr,aveint,qval,sigave,sigmaq,qbar,fsubs)
1727	String funcStr
1728	WAVE aveint,qval,sigave,sigmaq,qbar,fsubs
1729
1730	Variable r1,xCtr,yCtr,sdd,pixSize,wavelength
1731	String coefStr,abortStr,str
1732
1733	FUNCREF SANSModelAAO_MCproto func=$("fSmeared"+funcStr) //a wrapper for the structure version
1734	FUNCREF SANSModelAAO_MCproto funcUnsmeared=$(funcStr) //unsmeared
1735	coefStr = MC_getFunctionCoef(funcStr)
1736
1737	if(!MC_CheckFunctionAndCoef(funcStr,coefStr))
1738	Abort "Function and coefficients do not match. You must plot the unsmeared function before simulation."
1739	endif
1740
1741	Wave inten=$"root:Simulation:Simulation_i" // this will exist and send the smeared calculation to the corect DF
1742
1743	// the resolution-smeared intensity is calculated, including the incoherent background
1744	func($coefStr,inten,qval)
1745
1746	NVAR imon = root:Packages:NIST:SAS:gImon
1747	NVAR ctTime = root:Packages:NIST:SAS:gCntTime
1748	NVAR thick = root:Packages:NIST:SAS:gThick
1749	NVAR trans = root:Packages:NIST:SAS:gSamTrans
1750	NVAR SimDetCts = root:Packages:NIST:SAS:g_1DTotCts //summed counts (simulated)
1751	NVAR estDetCR = root:Packages:NIST:SAS:g_1DEstDetCR // estimated detector count rate
1752	NVAR fracScat = root:Packages:NIST:SAS:g_1DFracScatt // fraction of beam captured on detector
1753	NVAR estTrans = root:Packages:NIST:SAS:g_1DEstTrans // estimated transmission of sample
1754	NVAR mScat = root:Packages:NIST:SAS:g_MultScattFraction
1755	NVAR detectorEff = root:Packages:NIST:SAS:g_detectorEff
1756
1757	WAVE rw=root:Packages:NIST:SAS:realsRead
1758	WAVE nCells=root:Packages:NIST:SAS:nCells
1759
1760	pixSize = rw[10]/10 // convert pix size in mm to cm
1761	sdd = rw[18]*100 //convert header of [m] to [cm]
1762	wavelength = rw[26] // in 1/A
1763
1764	imon = beamIntensity()*ctTime
1765
1766	// calculate the scattering cross section simply to be able to estimate the transmission
1767	Variable sig_sas=0
1768
1769	// remember that the random deviate is the coherent portion ONLY - the incoherent background is
1770	// subtracted before the calculation.
1771	CalculateRandomDeviate(funcUnsmeared,$coefStr,wavelength,"root:Packages:NIST:SAS:ran_dev",sig_sas)
1772
1773	// if(sig_sas > 100)
1774	// sprintf abortStr,"sig_sas = %g. Please check that the model coefficients have a zero background, or the low q is well-behaved.",sig_sas
1775	// endif
1776
1777	// calculate the multiple scattering fraction for display (10/2009)
1778	Variable ii,nMax=10,tau
1779	mScat=0
1780	tau = thick*sig_sas
1781	// this sums the normalized scattering P', so the result is the fraction of multiply coherently scattered
1782	// neutrons out of those that were scattered
1783	for(ii=2;ii<nMax;ii+=1)
1784	mScat += tau^(ii)/factorial(ii)
1785	// print tau^(ii)/factorial(ii)
1786	endfor
1787	estTrans = exp(-1thicksig_sas) //thickness and sigma both in units of cm
1788	mscat *= (estTrans)/(1-estTrans)
1789
1790	if(mScat > 0.1) // /Z to supress error if this was a 1D calc with the 2D panel open
1791	ValDisplay/Z valdisp0_4 win=Sim_1D_Panel,labelBack=(65535,32768,32768)
1792	else
1793	ValDisplay/Z valdisp0_4 win=Sim_1D_Panel,labelBack=0
1794	endif
1795
1796
1797	Print "Sig_sas = ",sig_sas
1798
1799	Duplicate/O qval prob_i,countsInAnnulus
1800
1801	// not needed - nCells takes care of this when the error is correctly calculated
1802	// Duplicate/O qval circle_fraction,rval,nCells_expected
1803	// rval = sddtan(2asin(qval*wavelength/4/pi)) //radial distance in cm
1804	// nCells_expected = 2pirval/pixSize //does this need to be an integer?
1805	// circle_fraction = nCells / nCells_expected
1806
1807
1808	// prob_i = transthicknCells(pixSize/sdd)^2inten //probability of a neutron in q-bin(i) that has nCells
1809	prob_i = transthick(pixSize/sdd)^2*inten //probability of a neutron in q-bin(i)
1810
1811	Variable P_on = sum(prob_i,-inf,inf)
1812	Print "P_on = ",P_on
1813
1814	// fracScat = P_on
1815	fracScat = 1-estTrans
1816
1817	// aveint = (Imon)*prob_i / circle_fraction / nCells_expected
1818	// added correction for detector efficiency, since SASCALC is flux on sample
1819	aveint = (Imon)prob_idetectorEff
1820
1821	countsInAnnulus = aveint*nCells
1822	SimDetCts = sum(countsInAnnulus,-inf,inf)
1823	estDetCR = SimDetCts/ctTime
1824
1825
1826	NVAR doABS = root:Packages:NIST:SAS:g_1D_DoABS
1827	NVAR addNoise = root:Packages:NIST:SAS:g_1D_AddNoise
1828
1829	// this is where the number of cells comes in - the calculation of the error bars
1830	// sigma[i] = SUM(sigma[ij]^2) / nCells^2
1831	// and since in the simulation, SUM(sigma[ij]^2) = nCellssigma[ij]^2 = nCellsInten
1832	// then...
1833	sigave = sqrt(aveint/nCells) // corrected based on John's memo, from 8/9/99
1834
1835	// add in random error in aveint based on the sigave
1836	if(addNoise)
1837	aveint += gnoise(sigave)
1838	endif
1839
1840	// signature in the standard deviation, do this after the noise is added
1841	// start at 10 to be out of the beamstop (makes for nicer plotting)
1842	// end at 50 to leave the natural statistics at the end of the set (may have a total of 80+ points if no offset)
1843	sigave[10,50;10] = 10*sigave[p]
1844
1845	// convert to absolute scale, remembering to un-correct for the detector efficiency
1846	if(doABS)
1847	Variable kappa = thick(pixSize/sdd)^2trans*iMon
1848	aveint /= kappa
1849	sigave /= kappa
1850	aveint /= detectorEff
1851	sigave /= detectorEff
1852	endif
1853
1854
1855	return(0)
1856	End
1857
1858	/// for testing only
1859	Function testProbability(sas,thick)
1860	Variable sas,thick
1861
1862	Variable tau,trans,p2p,p3p,p4p
1863
1864	tau = sas*thick
1865	trans = exp(-tau)
1866
1867	Print "tau = ",tau
1868	Print "trans = ",trans
1869
1870	p2p = tau^2/factorial(2)*trans/(1-trans)
1871	p3p = tau^3/factorial(3)*trans/(1-trans)
1872	p4p = tau^4/factorial(4)*trans/(1-trans)
1873
1874	Print "double scattering = ",p2p
1875	Print "triple scattering = ",p3p
1876	Print "quadruple scattering = ",p4p
1877
1878	Variable ii,nMax=10,mScat=0
1879	for(ii=2;ii<nMax;ii+=1)
1880	mScat += tau^(ii)/factorial(ii)
1881	// print tau^(ii)/factorial(ii)
1882	endfor
1883	mscat *= (Trans)/(1-Trans)
1884
1885	Print "Total fraction of multiple scattering = ",mScat
1886
1887	return(mScat)
1888	End
1889
1890
1891	//// this is a very simple example of how to script the MC simulation to run unattended
1892	//
1893	// you need to supply for each "run": the run index (you increment manually)
1894	// the sample label (as a string)
1895	//
1896	// changing the various configuration paramters will have to be done on a case-by-case basis
1897	// looking into SASCALC to see what is really changed,
1898	// or the configuration parameters of the MC_SASCALC panel
1899	//
1900	//
1901	//Function Script_2DMC()
1902	//
1903	//
1904	// NVAR SimTimeWarn = root:Packages:NIST:SAS:g_SimTimeWarn
1905	// SimTimeWarn = 36000 //sets the threshold for the warning dialog to 10 hours
1906	// STRUCT WMButtonAction ba
1907	// ba.eventCode = 2 //fake mouse click on button
1908	//
1909	// NVAR detDist = root:Packages:NIST:SAS:gDetDist
1910	//
1911	// detDist = 200 //set directly in cm
1912	// MC_DoItButtonProc(ba)
1913	// SaveAsVAXButtonProc("",runIndex=105,simLabel="this is run 105, SDD = 200")
1914	//
1915	// detDist = 300 //set directly in cm
1916	// MC_DoItButtonProc(ba)
1917	// SaveAsVAXButtonProc("",runIndex=106,simLabel="this is run 106, SDD = 300")
1918	//
1919	// detDist = 400 //set directly in cm
1920	// MC_DoItButtonProc(ba)
1921	// SaveAsVAXButtonProc("",runIndex=107,simLabel="this is run 107, SDD = 400")
1922	//
1923	//
1924	// SimTimeWarn = 10 //back to 10 seconds for manual operation
1925	// return(0)
1926	//end

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format