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