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