1 | #pragma rtGlobals=1 // Use modern global access method. |
---|
2 | |
---|
3 | // |
---|
4 | // Monte Carlo simulator for SASCALC |
---|
5 | // October 2008 SRK |
---|
6 | // |
---|
7 | // This code simulates the scattering for a selected model based on the instrument configuration |
---|
8 | // This code is based directly on John Barker's code, which includes multiple scattering effects. |
---|
9 | // A lot of the setup, wave creation, and post-calculations are done in SASCALC->ReCalculateInten() |
---|
10 | // |
---|
11 | // |
---|
12 | // - Why am I off by a factor of 2.7 - 3.7 (MC too high) relative to real data? |
---|
13 | // I need to include efficiency (70%?) - do I knock these off before the simulation or do I |
---|
14 | // really simulate that some fraction of neutrons on the detector don't actually get counted? |
---|
15 | // Is the flux estimate up-to-date? !! Flux estimates at NG3 are out-of-date.... |
---|
16 | // - my simulated transmission is larger than what is measured, even after correcting for the quartz cell. |
---|
17 | // Why? Do I need to include absorption? Just inherent problems with incoherent cross sections? |
---|
18 | |
---|
19 | // - Most importantly, this needs to be checked for correctness of the MC simulation |
---|
20 | // X how can I get the "data" on absolute scale? This would be a great comparison vs. the ideal model calculation |
---|
21 | // X why does my integrated tau not match up with John's analytical calculations? where are the assumptions? |
---|
22 | // - get rid of all small angle assumptions - to make sure that the calculation is correct at all angles |
---|
23 | |
---|
24 | // |
---|
25 | // X at the larger angles, is the "flat" detector being properly accounted for - in terms of |
---|
26 | // the solid angle and how many counts fall in that pixel. Am I implicitly defining a spherical detector |
---|
27 | // so that what I see is already "corrected"? |
---|
28 | // X the MC will, of course benefit greatly from being XOPized. Maybe think about parallel implementation |
---|
29 | // by allowing the data arrays to accumulate. First pass at the XOP is done. Not pretty, not the speediest (5.8x) |
---|
30 | // but it is functional. Add spinCursor() for long calculations. See WaveAccess XOP example. |
---|
31 | // X the background parameter for the model MUST be zero, or the integration for scattering |
---|
32 | // 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 |
---|
33 | // X fully use the SASCALC input, most importantly, flux on sample. |
---|
34 | // X if no MC desired, still use the selected model |
---|
35 | // X better display of MC results on panel |
---|
36 | // X settings for "count for X seconds" or "how long to 1E6 cts on detector" (but 1E6 is typically too many counts...) |
---|
37 | // X warn of projected simulation time |
---|
38 | // - add quartz window scattering to the simulation somehow |
---|
39 | // -?- do smeared models make any sense?? Yes, John agrees that they do, and may be used in a more realistic simulation |
---|
40 | // -?- but the random deviate can't be properly calculated... |
---|
41 | // - make sure that the ratio of scattering coherent/incoherent is properly adjusted for the sample composition |
---|
42 | // or the volume fraction of solvent. |
---|
43 | // |
---|
44 | // X add to the results the fraction of coherently scattered neutrons that are singly scattered, different than |
---|
45 | // the overall fraction of singly scattered, and maybe more important to know. |
---|
46 | // |
---|
47 | // X change the fraction reaching the detector to exclude those that don't interact. These transmitted neutrons |
---|
48 | // aren't counted. Is the # that interact a better number? |
---|
49 | // |
---|
50 | // - do we want to NOT offset the data by a multiplicative factor as it is "frozen" , so that the |
---|
51 | // effects on the absolute scale can be seen? |
---|
52 | // |
---|
53 | // X why is "pure" incoherent scattering giving me a q^-1 slope, even with the detector all the way back? |
---|
54 | // -NO- can I speed up by assuming everything interacts? This would compromise the ability to calculate multiple scattering |
---|
55 | // X ask John how to verify what is going on |
---|
56 | // - a number of models are now found to be ill-behaved when q=1e-10. Then the random deviate calculation blows up. |
---|
57 | // a warning has been added - but some models need a proper limiting value, and some (power-law) are simply unuseable |
---|
58 | // unless something else can be done. |
---|
59 | // - if the MC gags on a simulation, it often gets "stuck" and can't do the normal calculation from the model, which it |
---|
60 | // should always default to... |
---|
61 | // |
---|
62 | // |
---|
63 | |
---|
64 | // threaded call to the main function, adds up the individual runs, and returns what is to be displayed |
---|
65 | // results is calculated and sent back for display |
---|
66 | Function Monte_SANS_Threaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) |
---|
67 | WAVE inputWave,ran_dev,nt,j1,j2,nn,linear_data,results |
---|
68 | |
---|
69 | //initialize ran1 in the XOP by passing a negative integer |
---|
70 | // does nothing in the Igor code |
---|
71 | Duplicate/O results retWave |
---|
72 | //results[0] = -1*(datetime) |
---|
73 | |
---|
74 | Variable NNeutron=inputWave[0] |
---|
75 | Variable i,nthreads= ThreadProcessorCount |
---|
76 | if(nthreads>2) //only support 2 processors until I can figure out how to properly thread the XOP and to loop it |
---|
77 | nthreads=2 |
---|
78 | endif |
---|
79 | |
---|
80 | // nthreads = 1 |
---|
81 | |
---|
82 | variable mt= ThreadGroupCreate(nthreads) |
---|
83 | |
---|
84 | inputWave[0] = NNeutron/nthreads //split up the number of neutrons |
---|
85 | |
---|
86 | for(i=0;i<nthreads;i+=1) |
---|
87 | Duplicate/O nt $("nt"+num2istr(i)) //new instance for each thread |
---|
88 | Duplicate/O j1 $("j1"+num2istr(i)) |
---|
89 | Duplicate/O j2 $("j2"+num2istr(i)) |
---|
90 | Duplicate/O nn $("nn"+num2istr(i)) |
---|
91 | Duplicate/O linear_data $("linear_data"+num2istr(i)) |
---|
92 | Duplicate/O retWave $("retWave"+num2istr(i)) |
---|
93 | Duplicate/O inputWave $("inputWave"+num2istr(i)) |
---|
94 | Duplicate/O ran_dev $("ran_dev"+num2istr(i)) |
---|
95 | |
---|
96 | // ?? I need explicit wave references? |
---|
97 | // maybe I need to have everything in separate data folders - bu tI haven't tried that. seems like a reach. |
---|
98 | // more likely there is something bad going on in the XOP code. |
---|
99 | if(i==0) |
---|
100 | WAVE inputWave0,ran_dev0,nt0,j10,j20,nn0,linear_data0,retWave0 |
---|
101 | retWave0[0] = -1*datetime //to initialize ran3 |
---|
102 | ThreadStart mt,i,Monte_SANS_W1(inputWave0,ran_dev0,nt0,j10,j20,nn0,linear_data0,retWave0) |
---|
103 | Print "started thread 0" |
---|
104 | endif |
---|
105 | if(i==1) |
---|
106 | WAVE inputWave1,ran_dev1,nt1,j11,j21,nn1,linear_data1,retWave1 |
---|
107 | //retWave1[0] = -1*datetime //to initialize ran3 |
---|
108 | ThreadStart mt,i,Monte_SANS_W2(inputWave1,ran_dev1,nt1,j11,j21,nn1,linear_data1,retWave1) |
---|
109 | Print "started thread 1" |
---|
110 | endif |
---|
111 | // if(i==2) |
---|
112 | // WAVE inputWave2,ran_dev2,nt2,j12,j22,nn2,linear_data2,retWave2 |
---|
113 | // retWave2[0] = -1*datetime //to initialize ran3 |
---|
114 | // ThreadStart mt,i,Monte_SANS_W(inputWave2,ran_dev2,nt2,j12,j22,nn2,linear_data2,retWave2) |
---|
115 | // endif |
---|
116 | // if(i==3) |
---|
117 | // WAVE inputWave3,ran_dev3,nt3,j13,j23,nn3,linear_data3,retWave3 |
---|
118 | // retWave3[0] = -1*datetime //to initialize ran3 |
---|
119 | // ThreadStart mt,i,Monte_SANS_W(inputWave3,ran_dev3,nt3,j13,j23,nn3,linear_data3,retWave3) |
---|
120 | // endif |
---|
121 | endfor |
---|
122 | |
---|
123 | // wait until done |
---|
124 | do |
---|
125 | variable tgs= ThreadGroupWait(mt,100) |
---|
126 | while( tgs != 0 ) |
---|
127 | variable dummy= ThreadGroupRelease(mt) |
---|
128 | Print "done with all threads" |
---|
129 | |
---|
130 | // calculate all of the bits for the results |
---|
131 | if(nthreads == 1) |
---|
132 | nt = nt0 // add up each instance |
---|
133 | j1 = j10 |
---|
134 | j2 = j20 |
---|
135 | nn = nn0 |
---|
136 | linear_data = linear_data0 |
---|
137 | retWave = retWave0 |
---|
138 | endif |
---|
139 | if(nthreads == 2) |
---|
140 | nt = nt0+nt1 // add up each instance |
---|
141 | j1 = j10+j11 |
---|
142 | j2 = j20+j21 |
---|
143 | nn = nn0+nn1 |
---|
144 | linear_data = linear_data0+linear_data1 |
---|
145 | retWave = retWave0+retWave1 |
---|
146 | endif |
---|
147 | // if(nthreads == 3) |
---|
148 | // nt = nt0+nt1+nt2 // add up each instance |
---|
149 | // j1 = j10+j11+j12 |
---|
150 | // j2 = j20+j21+j22 |
---|
151 | // nn = nn0+nn1+nn2 |
---|
152 | // linear_data = linear_data0+linear_data1+linear_data2 |
---|
153 | // retWave = retWave0+retWave1+retWave2 |
---|
154 | // endif |
---|
155 | // if(nthreads == 4) |
---|
156 | // nt = nt0+nt1+nt2+nt3 // add up each instance |
---|
157 | // j1 = j10+j11+j12+j13 |
---|
158 | // j2 = j20+j21+j22+j23 |
---|
159 | // nn = nn0+nn1+nn2+nn3 |
---|
160 | // linear_data = linear_data0+linear_data1+linear_data2+linear_data3 |
---|
161 | // retWave = retWave0+retWave1+retWave2+retWave3 |
---|
162 | // endif |
---|
163 | |
---|
164 | // fill up the results wave |
---|
165 | Variable xc,yc |
---|
166 | xc=inputWave[3] |
---|
167 | yc=inputWave[4] |
---|
168 | results[0] = inputWave[9]+inputWave[10] //total XS |
---|
169 | results[1] = inputWave[10] //SAS XS |
---|
170 | results[2] = retWave[1] //number that interact n2 |
---|
171 | results[3] = retWave[2] - linear_data[xc][yc] //# reaching detector minus Q(0) |
---|
172 | results[4] = retWave[3]/retWave[1] //avg# times scattered |
---|
173 | results[5] = retWave[4]/retWave[1] //single coherent fraction |
---|
174 | results[6] = retWave[5]/retWave[1] //double coherent fraction |
---|
175 | results[7] = retWave[6]/retWave[1] //multiple scatter fraction |
---|
176 | results[8] = (retWave[0]-retWave[1])/retWave[0] //transmitted fraction |
---|
177 | |
---|
178 | return(0) |
---|
179 | End |
---|
180 | |
---|
181 | // worker function for threads, does nothing except switch between XOP and Igor versions |
---|
182 | ThreadSafe Function Monte_SANS_W1(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) |
---|
183 | WAVE inputWave,ran_dev,nt,j1,j2,nn,linear_data,results |
---|
184 | |
---|
185 | #if exists("Monte_SANSX") |
---|
186 | Monte_SANSX(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) |
---|
187 | #else |
---|
188 | Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) |
---|
189 | #endif |
---|
190 | |
---|
191 | return (0) |
---|
192 | End |
---|
193 | // worker function for threads, does nothing except switch between XOP and Igor versions |
---|
194 | ThreadSafe Function Monte_SANS_W2(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) |
---|
195 | WAVE inputWave,ran_dev,nt,j1,j2,nn,linear_data,results |
---|
196 | |
---|
197 | #if exists("Monte_SANSX2") |
---|
198 | Monte_SANSX2(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) |
---|
199 | #else |
---|
200 | Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) |
---|
201 | #endif |
---|
202 | |
---|
203 | return (0) |
---|
204 | End |
---|
205 | |
---|
206 | // NON-threaded call to the main function returns what is to be displayed |
---|
207 | // results is calculated and sent back for display |
---|
208 | Function Monte_SANS_NotThreaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) |
---|
209 | WAVE inputWave,ran_dev,nt,j1,j2,nn,linear_data,results |
---|
210 | |
---|
211 | //initialize ran1 in the XOP by passing a negative integer |
---|
212 | // does nothing in the Igor code, enoise is already initialized |
---|
213 | Duplicate/O results retWave |
---|
214 | WAVE retWave |
---|
215 | retWave[0] = -1*abs(trunc(100000*enoise(1))) |
---|
216 | |
---|
217 | #if exists("Monte_SANSX") |
---|
218 | Monte_SANSX(inputWave,ran_dev,nt,j1,j2,nn,linear_data,retWave) |
---|
219 | #else |
---|
220 | Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,retWave) |
---|
221 | #endif |
---|
222 | |
---|
223 | // fill up the results wave |
---|
224 | Variable xc,yc |
---|
225 | xc=inputWave[3] |
---|
226 | yc=inputWave[4] |
---|
227 | results[0] = inputWave[9]+inputWave[10] //total XS |
---|
228 | results[1] = inputWave[10] //SAS XS |
---|
229 | results[2] = retWave[1] //number that interact n2 |
---|
230 | results[3] = retWave[2] - linear_data[xc][yc] //# reaching detector minus Q(0) |
---|
231 | results[4] = retWave[3]/retWave[1] //avg# times scattered |
---|
232 | results[5] = retWave[4]/retWave[1] //single coherent fraction |
---|
233 | results[6] = retWave[5]/retWave[1] //double coherent fraction |
---|
234 | results[7] = retWave[6]/retWave[1] //multiple scatter fraction |
---|
235 | results[8] = (retWave[0]-retWave[1])/retWave[0] //transmitted fraction |
---|
236 | |
---|
237 | return(0) |
---|
238 | End |
---|
239 | |
---|
240 | |
---|
241 | |
---|
242 | ////////// |
---|
243 | // PROGRAM Monte_SANS |
---|
244 | // PROGRAM simulates multiple SANS. |
---|
245 | // revised 2/12/99 JGB |
---|
246 | // added calculation of random deviate, and 2D 10/2008 SRK |
---|
247 | |
---|
248 | // N1 = NUMBER OF INCIDENT NEUTRONS. |
---|
249 | // N2 = NUMBER INTERACTED IN THE SAMPLE. |
---|
250 | // N3 = NUMBER ABSORBED. |
---|
251 | // THETA = SCATTERING ANGLE. |
---|
252 | |
---|
253 | // IMON = 'Enter number of neutrons to use in simulation.' |
---|
254 | // NUM_BINS = 'Enter number of THETA BINS TO use. (<5000).' |
---|
255 | // R1 = 'Enter beam radius. (cm)' |
---|
256 | // R2 = 'Enter sample radius. (cm)' |
---|
257 | // thick = 'Enter sample thickness. (cm)' |
---|
258 | // wavelength = 'Enter neutron wavelength. (A)' |
---|
259 | // R0 = 'Enter sphere radius. (A)' |
---|
260 | // |
---|
261 | |
---|
262 | ThreadSafe Function Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,MC_linear_data,results) |
---|
263 | WAVE inputWave,ran_dev,nt,j1,j2,nn,MC_linear_data,results |
---|
264 | |
---|
265 | Variable imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,sig_sas |
---|
266 | Variable NUM_BINS,N_INDEX |
---|
267 | Variable RHO,SIGSAS,SIGABS_0 |
---|
268 | Variable ii,jj,IND,idum,INDEX,IR,NQ |
---|
269 | Variable qmax,theta_max,R_DAB,R0,BOUND,I0,q0,zpow |
---|
270 | Variable N1,N2,N3,DTH,zz,tt,SIG_SINGLE,xx,yy,PHI,UU,SIG |
---|
271 | Variable THETA,Ran,ll,D_OMEGA,RR,Tabs,Ttot,I1_sumI |
---|
272 | Variable G0,E_NT,E_NN,TRANS_th,Trans_exp,rat |
---|
273 | Variable GG,GG_ED,dS_dW,ds_dw_double,ds_dw_single |
---|
274 | Variable DONE,FIND_THETA,err //used as logicals |
---|
275 | |
---|
276 | Variable Vx,Vy,Vz,Theta_z,qq |
---|
277 | Variable Sig_scat,Sig_abs,Ratio,Sig_total |
---|
278 | Variable isOn=0,testQ,testPhi,xPixel,yPixel |
---|
279 | Variable NSingleIncoherent,NSingleCoherent,NScatterEvents,incoherentEvent,coherentEvent |
---|
280 | Variable NDoubleCoherent,NMultipleScatter,countIt,detEfficiency |
---|
281 | |
---|
282 | detEfficiency = 1.0 //70% counting efficiency = 0.7 |
---|
283 | |
---|
284 | imon = inputWave[0] |
---|
285 | r1 = inputWave[1] |
---|
286 | r2 = inputWave[2] |
---|
287 | xCtr = inputWave[3] |
---|
288 | yCtr = inputWave[4] |
---|
289 | sdd = inputWave[5] |
---|
290 | pixSize = inputWave[6] |
---|
291 | thick = inputWave[7] |
---|
292 | wavelength = inputWave[8] |
---|
293 | sig_incoh = inputWave[9] |
---|
294 | sig_sas = inputWave[10] |
---|
295 | |
---|
296 | // SetRandomSeed 0.1 //to get a reproduceable sequence |
---|
297 | |
---|
298 | //scattering power and maximum qvalue to bin |
---|
299 | // zpow = .1 //scattering power, calculated below |
---|
300 | qmax = 4*pi/wavelength //maximum Q to bin 1D data. (A-1) (not really used, so set to a big value) |
---|
301 | sigabs_0 = 0.0 // ignore absorption cross section/wavelength [1/(cm A)] |
---|
302 | N_INDEX = 50 // maximum number of scattering events per neutron |
---|
303 | num_bins = 200 //number of 1-D bins (not really used) |
---|
304 | |
---|
305 | // my additions - calculate the random deviate function as needed |
---|
306 | // and calculate the scattering power from the model function (passed in as a wave) |
---|
307 | // |
---|
308 | Variable left = leftx(ran_dev) |
---|
309 | Variable delta = deltax(ran_dev) |
---|
310 | |
---|
311 | //c total SAS cross-section |
---|
312 | // SIG_SAS = zpow/thick |
---|
313 | zpow = sig_sas*thick //since I now calculate the sig_sas from the model |
---|
314 | SIG_ABS = SIGABS_0 * WAVElength |
---|
315 | sig_abs = 0.0 //cm-1 |
---|
316 | SIG_TOTAL =SIG_ABS + SIG_SAS + sig_incoh |
---|
317 | // Print "The TOTAL XSECTION. (CM-1) is ",sig_total |
---|
318 | // Print "The TOTAL SAS XSECTION. (CM-1) is ",sig_sas |
---|
319 | // results[0] = sig_total //assign these after everything's done |
---|
320 | // results[1] = sig_sas |
---|
321 | // variable ratio1,ratio2 |
---|
322 | // ratio1 = sig_abs/sig_total |
---|
323 | // ratio2 = sig_incoh/sig_total |
---|
324 | // // 0->ratio1 = abs |
---|
325 | // // ratio1 -> ratio2 = incoh |
---|
326 | // // > ratio2 = coherent |
---|
327 | RATIO = sig_incoh / SIG_TOTAL |
---|
328 | |
---|
329 | //c assuming theta = sin(theta)...OK |
---|
330 | theta_max = wavelength*qmax/(2*pi) |
---|
331 | //C SET Theta-STEP SIZE. |
---|
332 | DTH = Theta_max/NUM_BINS |
---|
333 | // Print "theta bin size = dth = ",dth |
---|
334 | |
---|
335 | //C INITIALIZE COUNTERS. |
---|
336 | N1 = 0 |
---|
337 | N2 = 0 |
---|
338 | N3 = 0 |
---|
339 | NSingleIncoherent = 0 |
---|
340 | NSingleCoherent = 0 |
---|
341 | NDoubleCoherent = 0 |
---|
342 | NMultipleScatter = 0 |
---|
343 | NScatterEvents = 0 |
---|
344 | |
---|
345 | //C INITIALIZE ARRAYS. |
---|
346 | j1 = 0 |
---|
347 | j2 = 0 |
---|
348 | nt = 0 |
---|
349 | nn=0 |
---|
350 | |
---|
351 | //C MONITOR LOOP - looping over the number of incedent neutrons |
---|
352 | //note that zz, is the z-position in the sample - NOT the scattering power |
---|
353 | |
---|
354 | // NOW, start the loop, throwing neutrons at the sample. |
---|
355 | do |
---|
356 | Vx = 0.0 // Initialize direction vector. |
---|
357 | Vy = 0.0 |
---|
358 | Vz = 1.0 |
---|
359 | |
---|
360 | Theta = 0.0 // Initialize scattering angle. |
---|
361 | Phi = 0.0 // Intialize azimuthal angle. |
---|
362 | N1 += 1 // Increment total number neutrons counter. |
---|
363 | DONE = 0 // True when neutron is scattered out of the sample. |
---|
364 | INDEX = 0 // Set counter for number of scattering events. |
---|
365 | zz = 0.0 // Set entering dimension of sample. |
---|
366 | incoherentEvent = 0 |
---|
367 | coherentEvent = 0 |
---|
368 | |
---|
369 | do // Makes sure position is within circle. |
---|
370 | ran = abs(enoise(1)) //[0,1] |
---|
371 | xx = 2.0*R1*(Ran-0.5) //X beam position of neutron entering sample. |
---|
372 | ran = abs(enoise(1)) //[0,1] |
---|
373 | yy = 2.0*R1*(Ran-0.5) //Y beam position ... |
---|
374 | RR = SQRT(xx*xx+yy*yy) //Radial position of neutron in incident beam. |
---|
375 | while(rr>r1) |
---|
376 | |
---|
377 | do //Scattering Loop, will exit when "done" == 1 |
---|
378 | // keep scattering multiple times until the neutron exits the sample |
---|
379 | ran = abs(enoise(1)) //[0,1] RANDOM NUMBER FOR DETERMINING PATH LENGTH |
---|
380 | ll = PATH_len(ran,Sig_total) |
---|
381 | //Determine new scattering direction vector. |
---|
382 | err = NewDirection(vx,vy,vz,Theta,Phi) //vx,vy,vz is updated, theta, phi unchanged by function |
---|
383 | |
---|
384 | //X,Y,Z-POSITION OF SCATTERING EVENT. |
---|
385 | xx += ll*vx |
---|
386 | yy += ll*vy |
---|
387 | zz += ll*vz |
---|
388 | RR = sqrt(xx*xx+yy*yy) //radial position of scattering event. |
---|
389 | |
---|
390 | //Check whether interaction occurred within sample volume. |
---|
391 | IF (((zz > 0.0) && (zz < THICK)) && (rr < r2)) |
---|
392 | //NEUTRON INTERACTED. |
---|
393 | ran = abs(enoise(1)) //[0,1] |
---|
394 | |
---|
395 | // if(ran<ratio1) |
---|
396 | // //absorption event |
---|
397 | // n3 +=1 |
---|
398 | // done=1 |
---|
399 | // else |
---|
400 | |
---|
401 | INDEX += 1 //Increment counter of scattering events. |
---|
402 | IF(INDEX == 1) |
---|
403 | N2 += 1 //Increment # of scat. neutrons |
---|
404 | Endif |
---|
405 | //Split neutron interactions into scattering and absorption events |
---|
406 | // IF(ran > (ratio1+ratio2) ) //C NEUTRON SCATTERED coherently |
---|
407 | IF(ran > ratio) //C NEUTRON SCATTERED coherently |
---|
408 | coherentEvent = 1 |
---|
409 | FIND_THETA = 0 //false |
---|
410 | DO |
---|
411 | //ran = abs(enoise(1)) //[0,1] |
---|
412 | //theta = Scat_angle(Ran,R_DAB,wavelength) // CHOOSE DAB ANGLE -- this is 2Theta |
---|
413 | //Q0 = 2*PI*THETA/WAVElength // John chose theta, calculated Q |
---|
414 | |
---|
415 | // pick a q-value from the deviate function |
---|
416 | // pnt2x truncates the point to an integer before returning the x |
---|
417 | // so get it from the wave scaling instead |
---|
418 | Q0 =left + binarysearchinterp(ran_dev,abs(enoise(1)))*delta |
---|
419 | theta = Q0/2/Pi*wavelength //SAS approximation |
---|
420 | |
---|
421 | //Print "q0, theta = ",q0,theta |
---|
422 | |
---|
423 | FIND_THETA = 1 //always accept |
---|
424 | |
---|
425 | while(!find_theta) |
---|
426 | ran = abs(enoise(1)) //[0,1] |
---|
427 | PHI = 2.0*PI*Ran //Chooses azimuthal scattering angle. |
---|
428 | ELSE |
---|
429 | //NEUTRON scattered incoherently |
---|
430 | // N3 += 1 |
---|
431 | // DONE = 1 |
---|
432 | // phi and theta are random over the entire sphere of scattering |
---|
433 | // !can't just choose random theta and phi, won't be random over sphere solid angle |
---|
434 | incoherentEvent = 1 |
---|
435 | |
---|
436 | ran = abs(enoise(1)) //[0,1] |
---|
437 | theta = acos(2*ran-1) |
---|
438 | |
---|
439 | ran = abs(enoise(1)) //[0,1] |
---|
440 | PHI = 2.0*PI*Ran //Chooses azimuthal scattering angle. |
---|
441 | ENDIF //(ran > ratio) |
---|
442 | // endif // event was absorption |
---|
443 | ELSE |
---|
444 | //NEUTRON ESCAPES FROM SAMPLE -- bin it somewhere |
---|
445 | DONE = 1 //done = true, will exit from loop |
---|
446 | |
---|
447 | // countIt = 1 |
---|
448 | // if(abs(enoise(1)) > detEfficiency) //efficiency of 70% wired @top |
---|
449 | // countIt = 0 //detector does not register |
---|
450 | // endif |
---|
451 | |
---|
452 | //Increment #scattering events array |
---|
453 | If (index <= N_Index) |
---|
454 | NN[INDEX] += 1 |
---|
455 | Endif |
---|
456 | |
---|
457 | if(index != 0) //the neutron interacted at least once, figure out where it ends up |
---|
458 | |
---|
459 | Theta_z = acos(Vz) // Angle WITH respect to z axis. |
---|
460 | testQ = 2*pi*sin(theta_z)/wavelength |
---|
461 | |
---|
462 | // pick a random phi angle, and see if it lands on the detector |
---|
463 | // since the scattering is isotropic, I can safely pick a new, random value |
---|
464 | // this would not be true if simulating anisotropic scattering. |
---|
465 | //testPhi = abs(enoise(1))*2*Pi |
---|
466 | testPhi = FindPhi(Vx,Vy) //use the exiting phi value as defined by Vx and Vy |
---|
467 | |
---|
468 | // is it on the detector? |
---|
469 | FindPixel(testQ,testPhi,wavelength,sdd,pixSize,xCtr,yCtr,xPixel,yPixel) |
---|
470 | |
---|
471 | if(xPixel != -1 && yPixel != -1) |
---|
472 | //if(index==1) // only the single scattering events |
---|
473 | MC_linear_data[xPixel][yPixel] += 1 //this is the total scattering, including multiple scattering |
---|
474 | //endif |
---|
475 | isOn += 1 // neutron that lands on detector |
---|
476 | endif |
---|
477 | |
---|
478 | If(theta_z < theta_max) |
---|
479 | //Choose index for scattering angle array. |
---|
480 | //IND = NINT(THETA_z/DTH + 0.4999999) |
---|
481 | ind = round(THETA_z/DTH + 0.4999999) //round is eqivalent to nint() |
---|
482 | NT[ind] += 1 //Increment bin for angle. |
---|
483 | //Increment angle array for single scattering events. |
---|
484 | IF(INDEX == 1) |
---|
485 | j1[ind] += 1 |
---|
486 | Endif |
---|
487 | //Increment angle array for double scattering events. |
---|
488 | IF (INDEX == 2) |
---|
489 | j2[ind] += 1 |
---|
490 | Endif |
---|
491 | EndIf |
---|
492 | |
---|
493 | // increment all of the counters now since done==1 here and I'm sure to exit and get another neutron |
---|
494 | NScatterEvents += index //total number of scattering events |
---|
495 | if(index == 1 && incoherentEvent == 1) |
---|
496 | NSingleIncoherent += 1 |
---|
497 | endif |
---|
498 | if(index == 1 && coherentEvent == 1) |
---|
499 | NSingleCoherent += 1 |
---|
500 | endif |
---|
501 | if(index == 2 && coherentEvent == 1 && incoherentEvent == 0) |
---|
502 | NDoubleCoherent += 1 |
---|
503 | endif |
---|
504 | if(index > 1) |
---|
505 | NMultipleScatter += 1 |
---|
506 | endif |
---|
507 | //Print "n1,index (x,y) = ",n1,index, xpixel,ypixel |
---|
508 | |
---|
509 | else // if neutron escaped without interacting |
---|
510 | |
---|
511 | // then it must be a transmitted neutron |
---|
512 | // don't need to calculate, just increment the proper counters |
---|
513 | MC_linear_data[xCtr][yCtr] += 1 |
---|
514 | isOn += 1 |
---|
515 | nt[0] += 1 |
---|
516 | |
---|
517 | endif //if interacted |
---|
518 | ENDIF |
---|
519 | while (!done) |
---|
520 | while(n1 < imon) |
---|
521 | |
---|
522 | // Print "Monte Carlo Done" |
---|
523 | results[0] = n1 |
---|
524 | results[1] = n2 |
---|
525 | results[2] = isOn |
---|
526 | results[3] = NScatterEvents //sum of # of times that neutrons scattered |
---|
527 | results[4] = NSingleCoherent //# of events that are single, coherent |
---|
528 | results[5] = NDoubleCoherent |
---|
529 | results[6] = NMultipleScatter //# of multiple scattering events |
---|
530 | |
---|
531 | // Print "# absorbed = ",n3 |
---|
532 | |
---|
533 | // trans_th = exp(-sig_total*thick) |
---|
534 | // TRANS_exp = (N1-N2) / N1 // Transmission |
---|
535 | // dsigma/domega assuming isotropic scattering, with no absorption. |
---|
536 | // Print "trans_exp = ",trans_exp |
---|
537 | // Print "total # of neutrons reaching 2D detector",isOn |
---|
538 | // Print "fraction of incident neutrons reaching detector ",isOn/iMon |
---|
539 | |
---|
540 | // Print "Total number of neutrons = ",N1 |
---|
541 | // Print "Total number of neutrons that interact = ",N2 |
---|
542 | // Print "Fraction of singly scattered neutrons = ",sum(j1,-inf,inf)/N2 |
---|
543 | // results[2] = N2 //number that scatter |
---|
544 | // results[3] = isOn - MC_linear_data[xCtr][yCtr] //# scattered reaching detector minus zero angle |
---|
545 | |
---|
546 | |
---|
547 | // Tabs = (N1-N3)/N1 |
---|
548 | // Ttot = (N1-N2)/N1 |
---|
549 | // I1_sumI = NN[0]/(N2-N3) |
---|
550 | // Print "Tabs = ",Tabs |
---|
551 | // Print "Transmitted neutrons = ",Ttot |
---|
552 | // results[8] = Ttot |
---|
553 | // Print "I1 / all I1 = ", I1_sumI |
---|
554 | |
---|
555 | End |
---|
556 | //////// end of main function for calculating multiple scattering |
---|
557 | |
---|
558 | |
---|
559 | // returns the random deviate as a wave |
---|
560 | // and the total SAS cross-section [1/cm] sig_sas |
---|
561 | Function CalculateRandomDeviate(func,coef,lam,outWave,SASxs) |
---|
562 | FUNCREF SANSModelAAO_MCproto func |
---|
563 | WAVE coef |
---|
564 | Variable lam |
---|
565 | String outWave |
---|
566 | Variable &SASxs |
---|
567 | |
---|
568 | Variable nPts_ran=10000,qu |
---|
569 | qu = 4*pi/lam |
---|
570 | |
---|
571 | // Make/O/N=(nPts_ran)/D root:Packages:NIST:SAS:Gq,root:Packages:NIST:SAS:xw // if these waves are 1000 pts, the results are "pixelated" |
---|
572 | // WAVE Gq = root:Packages:NIST:SAS:gQ |
---|
573 | // WAVE xw = root:Packages:NIST:SAS:xw |
---|
574 | |
---|
575 | // hard-wired into the Simulation directory rather than the SAS folder. |
---|
576 | // plotting resolution-smeared models won't work any other way |
---|
577 | Make/O/N=(nPts_ran)/D root:Simulation:Gq,root:Simulation:xw // if these waves are 1000 pts, the results are "pixelated" |
---|
578 | WAVE Gq = root:Simulation:gQ |
---|
579 | WAVE xw = root:Simulation:xw |
---|
580 | SetScale/I x (0+1e-6),qu*(1-1e-10),"", Gq,xw //don't start at zero or run up all the way to qu to avoid numerical errors |
---|
581 | |
---|
582 | /// |
---|
583 | /// if all of the coefficients are well-behaved, then the last point is the background |
---|
584 | // and I can set it to zero here (only for the calculation) |
---|
585 | Duplicate/O coef,tmp_coef |
---|
586 | Variable num=numpnts(coef) |
---|
587 | tmp_coef[num-1] = 0 |
---|
588 | |
---|
589 | xw=x //for the AAO |
---|
590 | func(tmp_coef,Gq,xw) //call as AAO |
---|
591 | |
---|
592 | // Gq = x*Gq // SAS approximation |
---|
593 | Gq = Gq*sin(2*asin(x/qu))/sqrt(1-(x/qu)) // exact |
---|
594 | // |
---|
595 | Integrate/METH=1 Gq/D=Gq_INT |
---|
596 | |
---|
597 | // SASxs = lam*lam/2/pi*Gq_INT[nPts_ran-1] //if the approximation is used |
---|
598 | SASxs = lam*Gq_INT[nPts_ran-1] |
---|
599 | |
---|
600 | Gq_INT /= Gq_INT[nPts_ran-1] |
---|
601 | |
---|
602 | Duplicate/O Gq_INT $outWave |
---|
603 | |
---|
604 | return(0) |
---|
605 | End |
---|
606 | |
---|
607 | |
---|
608 | |
---|
609 | ThreadSafe Function FindPixel(testQ,testPhi,lam,sdd,pixSize,xCtr,yCtr,xPixel,yPixel) |
---|
610 | Variable testQ,testPhi,lam,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel |
---|
611 | |
---|
612 | Variable theta,dy,dx,qx,qy |
---|
613 | //decompose to qx,qy |
---|
614 | qx = testQ*cos(testPhi) |
---|
615 | qy = testQ*sin(testPhi) |
---|
616 | |
---|
617 | //convert qx,qy to pixel locations relative to # of pixels x, y from center |
---|
618 | theta = 2*asin(qy*lam/4/pi) |
---|
619 | dy = sdd*tan(theta) |
---|
620 | yPixel = round(yCtr + dy/pixSize) |
---|
621 | |
---|
622 | theta = 2*asin(qx*lam/4/pi) |
---|
623 | dx = sdd*tan(theta) |
---|
624 | xPixel = round(xCtr + dx/pixSize) |
---|
625 | |
---|
626 | //if on detector, return xPix and yPix values, otherwise -1 |
---|
627 | if(yPixel > 127 || yPixel < 0) |
---|
628 | yPixel = -1 |
---|
629 | endif |
---|
630 | if(xPixel > 127 || xPixel < 0) |
---|
631 | xPixel = -1 |
---|
632 | endif |
---|
633 | |
---|
634 | return(0) |
---|
635 | End |
---|
636 | |
---|
637 | Function MC_CheckFunctionAndCoef(funcStr,coefStr) |
---|
638 | String funcStr,coefStr |
---|
639 | |
---|
640 | SVAR/Z listStr=root:Packages:NIST:coefKWStr |
---|
641 | if(SVAR_Exists(listStr) == 1) |
---|
642 | String properCoefStr = StringByKey(funcStr, listStr ,"=",";",0) |
---|
643 | if(cmpstr("",properCoefStr)==0) |
---|
644 | return(0) //false, no match found, so properCoefStr is returned null |
---|
645 | endif |
---|
646 | if(cmpstr(coefStr,properCoefStr)==0) |
---|
647 | return(1) //true, the coef is the correct match |
---|
648 | endif |
---|
649 | endif |
---|
650 | return(0) //false, wrong coef |
---|
651 | End |
---|
652 | |
---|
653 | Function/S MC_getFunctionCoef(funcStr) |
---|
654 | String funcStr |
---|
655 | |
---|
656 | SVAR/Z listStr=root:Packages:NIST:coefKWStr |
---|
657 | String coefStr="" |
---|
658 | if(SVAR_Exists(listStr) == 1) |
---|
659 | coefStr = StringByKey(funcStr, listStr ,"=",";",0) |
---|
660 | endif |
---|
661 | return(coefStr) |
---|
662 | End |
---|
663 | |
---|
664 | Function SANSModelAAO_MCproto(w,yw,xw) |
---|
665 | Wave w,yw,xw |
---|
666 | |
---|
667 | Print "in SANSModelAAO_MCproto function" |
---|
668 | return(1) |
---|
669 | end |
---|
670 | |
---|
671 | Function/S MC_FunctionPopupList() |
---|
672 | String list,tmp |
---|
673 | list = FunctionList("*",";","KIND:10") //get every user defined curve fit function |
---|
674 | |
---|
675 | //now start to remove everything the user doesn't need to see... |
---|
676 | |
---|
677 | tmp = FunctionList("*_proto",";","KIND:10") //prototypes |
---|
678 | list = RemoveFromList(tmp, list ,";") |
---|
679 | //prototypes that show up if GF is loaded |
---|
680 | list = RemoveFromList("GFFitFuncTemplate", list) |
---|
681 | list = RemoveFromList("GFFitAllAtOnceTemplate", list) |
---|
682 | list = RemoveFromList("NewGlblFitFunc", list) |
---|
683 | list = RemoveFromList("NewGlblFitFuncAllAtOnce", list) |
---|
684 | list = RemoveFromList("GlobalFitFunc", list) |
---|
685 | list = RemoveFromList("GlobalFitAllAtOnce", list) |
---|
686 | list = RemoveFromList("GFFitAAOStructTemplate", list) |
---|
687 | list = RemoveFromList("NewGF_SetXWaveInList", list) |
---|
688 | list = RemoveFromList("NewGlblFitFuncAAOStruct", list) |
---|
689 | |
---|
690 | // more to remove as a result of 2D/Gizmo |
---|
691 | list = RemoveFromList("A_WMRunLessThanDelta", list) |
---|
692 | list = RemoveFromList("WMFindNaNValue", list) |
---|
693 | list = RemoveFromList("WM_Make3DBarChartParametricWave", list) |
---|
694 | list = RemoveFromList("UpdateQxQy2Mat", list) |
---|
695 | list = RemoveFromList("MakeBSMask", list) |
---|
696 | |
---|
697 | // MOTOFIT/GenFit bits |
---|
698 | tmp = "GEN_allatoncefitfunc;GEN_fitfunc;GetCheckBoxesState;MOTO_GFFitAllAtOnceTemplate;MOTO_GFFitFuncTemplate;MOTO_NewGF_SetXWaveInList;MOTO_NewGlblFitFunc;MOTO_NewGlblFitFuncAllAtOnce;" |
---|
699 | list = RemoveFromList(tmp, list ,";") |
---|
700 | |
---|
701 | // SANS Reduction bits |
---|
702 | tmp = "ASStandardFunction;Ann_1D_Graph;Avg_1D_Graph;BStandardFunction;CStandardFunction;Draw_Plot1D;MyMat2XYZ;NewDirection;SANSModelAAO_MCproto;Monte_SANS_Threaded;Monte_SANS_NotThreaded;Monte_SANS_W1;Monte_SANS_W2;" |
---|
703 | list = RemoveFromList(tmp, list ,";") |
---|
704 | list = RemoveFromList("Monte_SANS", list) |
---|
705 | |
---|
706 | tmp = FunctionList("f*",";","NPARAMS:2") //point calculations |
---|
707 | list = RemoveFromList(tmp, list ,";") |
---|
708 | |
---|
709 | tmp = FunctionList("fSmear*",";","NPARAMS:3") //smeared dependency calculations |
---|
710 | list = RemoveFromList(tmp, list ,";") |
---|
711 | |
---|
712 | //non-fit functions that I can't seem to filter out |
---|
713 | list = RemoveFromList("BinaryHS_PSF11;BinaryHS_PSF12;BinaryHS_PSF22;EllipCyl_Integrand;PP_Inner;PP_Outer;Phi_EC;TaE_Inner;TaE_Outer;",list,";") |
---|
714 | //////////////// |
---|
715 | |
---|
716 | //more functions from analysis models (2008) |
---|
717 | tmp = "Barbell_Inner;Barbell_Outer;Barbell_integrand;BCC_Integrand;Integrand_BCC_Inner;Integrand_BCC_Outer;" |
---|
718 | list = RemoveFromList(tmp, list ,";") |
---|
719 | tmp = "CapCyl;CapCyl_Inner;CapCyl_Outer;ConvLens;ConvLens_Inner;ConvLens_Outer;" |
---|
720 | list = RemoveFromList(tmp, list ,";") |
---|
721 | tmp = "Dumb;Dumb_Inner;Dumb_Outer;FCC_Integrand;Integrand_FCC_Inner;Integrand_FCC_Outer;" |
---|
722 | list = RemoveFromList(tmp, list ,";") |
---|
723 | tmp = "Integrand_SC_Inner;Integrand_SC_Outer;SC_Integrand;SphCyl;SphCyl_Inner;SphCyl_Outer;" |
---|
724 | list = RemoveFromList(tmp, list ,";") |
---|
725 | |
---|
726 | //simplify the display, forcing smeared calculations behind the scenes |
---|
727 | tmp = FunctionList("Smear*",";","NPARAMS:1") //smeared dependency calculations |
---|
728 | list = RemoveFromList(tmp, list ,";") |
---|
729 | |
---|
730 | if(strlen(list)==0) |
---|
731 | list = "No functions plotted" |
---|
732 | endif |
---|
733 | |
---|
734 | list = SortList(list) |
---|
735 | |
---|
736 | list = "default;"+list |
---|
737 | return(list) |
---|
738 | End |
---|
739 | |
---|
740 | |
---|
741 | //Function Scat_Angle(Ran,R_DAB,wavelength) |
---|
742 | // Variable Ran,r_dab,wavelength |
---|
743 | // |
---|
744 | // Variable qq,arg,theta |
---|
745 | // qq = 1. / ( R_DAB*sqrt(1.0/Ran - 1.0) ) |
---|
746 | // arg = qq*wavelength/(4*pi) |
---|
747 | // If (arg < 1.0) |
---|
748 | // theta = 2.*asin(arg) |
---|
749 | // else |
---|
750 | // theta = pi |
---|
751 | // endif |
---|
752 | // Return (theta) |
---|
753 | //End |
---|
754 | |
---|
755 | //calculates new direction (xyz) from an old direction |
---|
756 | //theta and phi don't change |
---|
757 | ThreadSafe Function NewDirection(vx,vy,vz,theta,phi) |
---|
758 | Variable &vx,&vy,&vz |
---|
759 | Variable theta,phi |
---|
760 | |
---|
761 | Variable err=0,vx0,vy0,vz0 |
---|
762 | Variable nx,ny,mag_xy,tx,ty,tz |
---|
763 | |
---|
764 | //store old direction vector |
---|
765 | vx0 = vx |
---|
766 | vy0 = vy |
---|
767 | vz0 = vz |
---|
768 | |
---|
769 | mag_xy = sqrt(vx0*vx0 + vy0*vy0) |
---|
770 | if(mag_xy < 1e-12) |
---|
771 | //old vector lies along beam direction |
---|
772 | nx = 0 |
---|
773 | ny = 1 |
---|
774 | tx = 1 |
---|
775 | ty = 0 |
---|
776 | tz = 0 |
---|
777 | else |
---|
778 | Nx = -Vy0 / Mag_XY |
---|
779 | Ny = Vx0 / Mag_XY |
---|
780 | Tx = -Vz0*Vx0 / Mag_XY |
---|
781 | Ty = -Vz0*Vy0 / Mag_XY |
---|
782 | Tz = Mag_XY |
---|
783 | endif |
---|
784 | |
---|
785 | //new scattered direction vector |
---|
786 | Vx = cos(phi)*sin(theta)*Tx + sin(phi)*sin(theta)*Nx + cos(theta)*Vx0 |
---|
787 | Vy = cos(phi)*sin(theta)*Ty + sin(phi)*sin(theta)*Ny + cos(theta)*Vy0 |
---|
788 | Vz = cos(phi)*sin(theta)*Tz + cos(theta)*Vz0 |
---|
789 | |
---|
790 | Return(err) |
---|
791 | End |
---|
792 | |
---|
793 | ThreadSafe Function path_len(aval,sig_tot) |
---|
794 | Variable aval,sig_tot |
---|
795 | |
---|
796 | Variable retval |
---|
797 | |
---|
798 | retval = -1*ln(1-aval)/sig_tot |
---|
799 | |
---|
800 | return(retval) |
---|
801 | End |
---|
802 | |
---|
803 | // globals are initialized in SASCALC.ipf |
---|
804 | Window MC_SASCALC() : Panel |
---|
805 | PauseUpdate; Silent 1 // building window... |
---|
806 | NewPanel /W=(92,556,390,1028)/K=1 as "SANS Simulator" |
---|
807 | SetVariable MC_setvar0,pos={28,73},size={144,15},bodyWidth=80,title="# of neutrons" |
---|
808 | SetVariable MC_setvar0,format="%5.4g" |
---|
809 | SetVariable MC_setvar0,limits={-inf,inf,100},value= root:Packages:NIST:SAS:gImon |
---|
810 | SetVariable MC_setvar0_1,pos={28,119},size={131,15},bodyWidth=60,title="Thickness (cm)" |
---|
811 | SetVariable MC_setvar0_1,limits={-inf,inf,0.1},value= root:Packages:NIST:SAS:gThick |
---|
812 | SetVariable MC_setvar0_2,pos={28,96},size={149,15},bodyWidth=60,title="Incoherent XS (cm)" |
---|
813 | SetVariable MC_setvar0_2,limits={-inf,inf,0.1},value= root:Packages:NIST:SAS:gSig_incoh |
---|
814 | SetVariable MC_setvar0_3,pos={28,142},size={150,15},bodyWidth=60,title="Sample Radius (cm)" |
---|
815 | SetVariable MC_setvar0_3,limits={-inf,inf,0.1},value= root:Packages:NIST:SAS:gR2 |
---|
816 | PopupMenu MC_popup0,pos={13,13},size={165,20},proc=MC_ModelPopMenuProc,title="Model Function" |
---|
817 | PopupMenu MC_popup0,mode=1,value= #"MC_FunctionPopupList()" |
---|
818 | Button MC_button0,pos={17,181},size={130,20},proc=MC_DoItButtonProc,title="Do MC Simulation" |
---|
819 | Button MC_button1,pos={181,181},size={80,20},proc=MC_Display2DButtonProc,title="Show 2D" |
---|
820 | SetVariable setvar0_3,pos={105,484},size={50,20},disable=1 |
---|
821 | GroupBox group0,pos={15,42},size={267,130},title="Monte Carlo" |
---|
822 | SetVariable cntVar,pos={190,73},size={80,15},proc=CountTimeSetVarProc,title="time(s)" |
---|
823 | SetVariable cntVar,format="%d" |
---|
824 | SetVariable cntVar,limits={1,10,1},value= root:Packages:NIST:SAS:gCntTime |
---|
825 | |
---|
826 | String fldrSav0= GetDataFolder(1) |
---|
827 | SetDataFolder root:Packages:NIST:SAS: |
---|
828 | Edit/W=(13,217,283,450)/HOST=# results_desc,results |
---|
829 | ModifyTable format(Point)=1,width(Point)=0,width(results_desc)=150 |
---|
830 | SetDataFolder fldrSav0 |
---|
831 | RenameWindow #,T_results |
---|
832 | SetActiveSubwindow ## |
---|
833 | EndMacro |
---|
834 | |
---|
835 | |
---|
836 | Function CountTimeSetVarProc(sva) : SetVariableControl |
---|
837 | STRUCT WMSetVariableAction &sva |
---|
838 | |
---|
839 | switch( sva.eventCode ) |
---|
840 | case 1: // mouse up |
---|
841 | case 2: // Enter key |
---|
842 | case 3: // Live update |
---|
843 | Variable dval = sva.dval |
---|
844 | |
---|
845 | // get the neutron flux, multiply, and reset the global for # neutrons |
---|
846 | NVAR imon=root:Packages:NIST:SAS:gImon |
---|
847 | imon = dval*beamIntensity() |
---|
848 | |
---|
849 | break |
---|
850 | endswitch |
---|
851 | |
---|
852 | return 0 |
---|
853 | End |
---|
854 | |
---|
855 | |
---|
856 | Function MC_ModelPopMenuProc(pa) : PopupMenuControl |
---|
857 | STRUCT WMPopupAction &pa |
---|
858 | |
---|
859 | switch( pa.eventCode ) |
---|
860 | case 2: // mouse up |
---|
861 | Variable popNum = pa.popNum |
---|
862 | String popStr = pa.popStr |
---|
863 | SVAR gStr = root:Packages:NIST:SAS:gFuncStr |
---|
864 | gStr = popStr |
---|
865 | |
---|
866 | break |
---|
867 | endswitch |
---|
868 | |
---|
869 | return 0 |
---|
870 | End |
---|
871 | |
---|
872 | Function MC_DoItButtonProc(ba) : ButtonControl |
---|
873 | STRUCT WMButtonAction &ba |
---|
874 | |
---|
875 | switch( ba.eventCode ) |
---|
876 | case 2: // mouse up |
---|
877 | // click code here |
---|
878 | NVAR doMC = root:Packages:NIST:SAS:gDoMonteCarlo |
---|
879 | doMC = 1 |
---|
880 | ReCalculateInten(1) |
---|
881 | doMC = 0 //so the next time won't be MC |
---|
882 | break |
---|
883 | endswitch |
---|
884 | |
---|
885 | return 0 |
---|
886 | End |
---|
887 | |
---|
888 | |
---|
889 | Function MC_Display2DButtonProc(ba) : ButtonControl |
---|
890 | STRUCT WMButtonAction &ba |
---|
891 | |
---|
892 | switch( ba.eventCode ) |
---|
893 | case 2: // mouse up |
---|
894 | // click code here |
---|
895 | Execute "ChangeDisplay(\"SAS\")" |
---|
896 | break |
---|
897 | endswitch |
---|
898 | |
---|
899 | return 0 |
---|
900 | End |
---|
901 | |
---|
902 | // after a 2d data image is averaged in the usual way, take the waves and generate a "fake" folder of the 1d |
---|
903 | // data, to appear as if it was loaded from a real data file. |
---|
904 | // |
---|
905 | // currently only works with SANS data, but can later be expanded to generate fake USANS data sets |
---|
906 | // |
---|
907 | Function Fake1DDataFolder(qval,aveint,sigave,sigmaQ,qbar,fSubs,dataFolder) |
---|
908 | WAVE qval,aveint,sigave,sigmaQ,qbar,fSubs |
---|
909 | String dataFolder |
---|
910 | |
---|
911 | String baseStr=dataFolder |
---|
912 | if(DataFolderExists("root:"+baseStr)) |
---|
913 | SetDataFolder $("root:"+baseStr) |
---|
914 | else |
---|
915 | NewDataFolder/S $("root:"+baseStr) |
---|
916 | endif |
---|
917 | |
---|
918 | ////overwrite the existing data, if it exists |
---|
919 | Duplicate/O qval, $(baseStr+"_q") |
---|
920 | Duplicate/O aveint, $(baseStr+"_i") |
---|
921 | Duplicate/O sigave, $(baseStr+"_s") |
---|
922 | // Duplicate/O sigmaQ, $(baseStr+"sq") |
---|
923 | // Duplicate/O qbar, $(baseStr+"qb") |
---|
924 | // Duplicate/O fSubS, $(baseStr+"fs") |
---|
925 | |
---|
926 | // need to switch based on SANS/USANS |
---|
927 | if (isSANSResolution(sigave[0])) //checks to see if the first point of the wave is <0] |
---|
928 | // make a resolution matrix for SANS data |
---|
929 | Variable np=numpnts(qval) |
---|
930 | Make/D/O/N=(np,4) $(baseStr+"_res") |
---|
931 | Wave res=$(baseStr+"_res") |
---|
932 | |
---|
933 | res[][0] = sigmaQ[p] //sigQ |
---|
934 | res[][1] = qBar[p] //qBar |
---|
935 | res[][2] = fSubS[p] //fShad |
---|
936 | res[][3] = qval[p] //Qvalues |
---|
937 | |
---|
938 | // keep a copy of everything in SAS too... the smearing wrapper function looks for |
---|
939 | // data in folders based on waves it is passed - an I lose control of that |
---|
940 | Duplicate/O res, $("root:Packages:NIST:SAS:"+baseStr+"_res") |
---|
941 | Duplicate/O qval, $("root:Packages:NIST:SAS:"+baseStr+"_q") |
---|
942 | Duplicate/O aveint, $("root:Packages:NIST:SAS:"+baseStr+"_i") |
---|
943 | Duplicate/O sigave, $("root:Packages:NIST:SAS:"+baseStr+"_s") |
---|
944 | else |
---|
945 | //the data is USANS data |
---|
946 | // nothing done here yet |
---|
947 | // dQv = -$w3[0] |
---|
948 | |
---|
949 | // USANS_CalcWeights(baseStr,dQv) |
---|
950 | |
---|
951 | endif |
---|
952 | |
---|
953 | //clean up |
---|
954 | SetDataFolder root: |
---|
955 | |
---|
956 | End |
---|
957 | |
---|
958 | |
---|
959 | |
---|
960 | /////UNUSED, testing routines that have not been updated to work with SASCALC |
---|
961 | // |
---|
962 | //Macro Simulate2D_MonteCarlo(imon,r1,r2,xCtr,yCtr,sdd,thick,wavelength,sig_incoh,funcStr) |
---|
963 | // Variable imon=100000,r1=0.6,r2=0.8,xCtr=100,yCtr=64,sdd=400,thick=0.1,wavelength=6,sig_incoh=0.1 |
---|
964 | // String funcStr |
---|
965 | // Prompt funcStr, "Pick the model function", popup, MC_FunctionPopupList() |
---|
966 | // |
---|
967 | // String coefStr = MC_getFunctionCoef(funcStr) |
---|
968 | // Variable pixSize = 0.5 // can't have 11 parameters in macro! |
---|
969 | // |
---|
970 | // if(!MC_CheckFunctionAndCoef(funcStr,coefStr)) |
---|
971 | // Abort "The coefficients and function type do not match. Please correct the selections in the popup menus." |
---|
972 | // endif |
---|
973 | // |
---|
974 | // Make/O/D/N=10 root:Packages:NIST:SAS:results |
---|
975 | // Make/O/T/N=10 root:Packages:NIST:SAS:results_desc |
---|
976 | // |
---|
977 | // RunMonte(imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,funcStr,coefStr,results) |
---|
978 | // |
---|
979 | //End |
---|
980 | // |
---|
981 | //Function RunMonte(imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,funcStr,coefStr) |
---|
982 | // Variable imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh |
---|
983 | // String funcStr,coefStr |
---|
984 | // WAVE results |
---|
985 | // |
---|
986 | // FUNCREF SANSModelAAO_MCproto func=$funcStr |
---|
987 | // |
---|
988 | // Monte_SANS(imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,sig_sas,ran_dev,linear_data,results) |
---|
989 | //End |
---|
990 | // |
---|
991 | ////// END UNUSED BLOCK |
---|