source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/MultScatter_MonteCarlo_2D.ipf @ 622

Last change on this file since 622 was 622, checked in by srkline, 13 years ago

Updated MonteCarlo? to use a maximum of 4 processors. If XOP is not present (with 4 functions), the Igor code reverts to a single thread, even reverting to ipf code if necessary.

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