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

Last change on this file since 544 was 544, checked in by srkline, 14 years ago

Lots of changes to add the first version of a USANS simulator, like SASCALC

(!) first issue is with the entanglement of dependencies - need to load SANS macros first!

otherwise, worth a first test for interface and accuracy. behavior is similar to the 1D SASCALC

File size: 45.6 KB
Line 
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. Using a log-spacing of points doesn't seem to help, and it introduces a lot of
59//   other problems. Not the way to go.
60// - if the MC gags on a simulation, it often gets "stuck" and can't do the normal calculation from the model, which it
61//   should always default to...
62//
63//
64
65// threaded call to the main function, adds up the individual runs, and returns what is to be displayed
66// results is calculated and sent back for display
67Function Monte_SANS_Threaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
68        WAVE inputWave,ran_dev,nt,j1,j2,nn,linear_data,results
69
70        //initialize ran1 in the XOP by passing a negative integer
71        // does nothing in the Igor code
72        Duplicate/O results retWave
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        NVAR gInitTime = root:Packages:NIST:SAS:gRanDateTime            //time that SASCALC was started
84
85       
86        inputWave[0] = NNeutron/nthreads                //split up the number of neutrons
87       
88        for(i=0;i<nthreads;i+=1)
89                Duplicate/O nt $("nt"+num2istr(i))              //new instance for each thread
90                Duplicate/O j1 $("j1"+num2istr(i))
91                Duplicate/O j2 $("j2"+num2istr(i))
92                Duplicate/O nn $("nn"+num2istr(i))
93                Duplicate/O linear_data $("linear_data"+num2istr(i))
94                Duplicate/O retWave $("retWave"+num2istr(i))
95                Duplicate/O inputWave $("inputWave"+num2istr(i))
96                Duplicate/O ran_dev $("ran_dev"+num2istr(i))
97                               
98                // ?? I need explicit wave references?
99                // maybe I need to have everything in separate data folders - bu tI haven't tried that. seems like a reach.
100                // more likely there is something bad going on in the XOP code.
101                if(i==0)
102                        WAVE inputWave0,ran_dev0,nt0,j10,j20,nn0,linear_data0,retWave0
103                        retWave0 = 0            //clear the return wave
104                        retWave0[0] = -1*(datetime-gInitTime)           //to initialize ran3
105                        ThreadStart mt,i,Monte_SANS_W1(inputWave0,ran_dev0,nt0,j10,j20,nn0,linear_data0,retWave0)
106                        Print "started thread 0"
107                endif
108                if(i==1)
109                        WAVE inputWave1,ran_dev1,nt1,j11,j21,nn1,linear_data1,retWave1
110                        retWave1 = 0                    //clear the return wave
111                        retWave1[0] = -1*(datetime-gInitTime)           //to initialize ran1
112                        ThreadStart mt,i,Monte_SANS_W2(inputWave1,ran_dev1,nt1,j11,j21,nn1,linear_data1,retWave1)
113                        Print "started thread 1"
114                endif
115//              if(i==2)
116//                      WAVE inputWave2,ran_dev2,nt2,j12,j22,nn2,linear_data2,retWave2
117//                      retWave2[0] = -1*datetime               //to initialize ran3
118//                      ThreadStart mt,i,Monte_SANS_W(inputWave2,ran_dev2,nt2,j12,j22,nn2,linear_data2,retWave2)
119//              endif
120//              if(i==3)
121//                      WAVE inputWave3,ran_dev3,nt3,j13,j23,nn3,linear_data3,retWave3
122//                      retWave3[0] = -1*datetime               //to initialize ran3
123//                      ThreadStart mt,i,Monte_SANS_W(inputWave3,ran_dev3,nt3,j13,j23,nn3,linear_data3,retWave3)
124//              endif
125        endfor
126
127// wait until done
128        do
129                variable tgs= ThreadGroupWait(mt,100)
130        while( tgs != 0 )
131        variable dummy= ThreadGroupRelease(mt)
132        Print "done with all threads"
133
134        // calculate all of the bits for the results
135        if(nthreads == 1)
136                nt = nt0                // add up each instance
137                j1 = j10
138                j2 = j20
139                nn = nn0
140                linear_data = linear_data0
141                retWave = retWave0
142        endif
143        if(nthreads == 2)
144                nt = nt0+nt1            // add up each instance
145                j1 = j10+j11
146                j2 = j20+j21
147                nn = nn0+nn1
148                linear_data = linear_data0+linear_data1
149                retWave = retWave0+retWave1
150        endif
151//      if(nthreads == 3)
152//              nt = nt0+nt1+nt2                // add up each instance
153//              j1 = j10+j11+j12
154//              j2 = j20+j21+j22
155//              nn = nn0+nn1+nn2
156//              linear_data = linear_data0+linear_data1+linear_data2
157//              retWave = retWave0+retWave1+retWave2
158//      endif
159//      if(nthreads == 4)
160//              nt = nt0+nt1+nt2+nt3            // add up each instance
161//              j1 = j10+j11+j12+j13
162//              j2 = j20+j21+j22+j23
163//              nn = nn0+nn1+nn2+nn3
164//              linear_data = linear_data0+linear_data1+linear_data2+linear_data3
165//              retWave = retWave0+retWave1+retWave2+retWave3
166//      endif
167       
168        // fill up the results wave
169        Variable xc,yc
170        xc=inputWave[3]
171        yc=inputWave[4]
172        results[0] = inputWave[9]+inputWave[10]         //total XS
173        results[1] = inputWave[10]                                              //SAS XS
174        results[2] = retWave[1]                                                 //number that interact n2
175        results[3] = retWave[2] - linear_data[xc][yc]                           //# reaching detector minus Q(0)
176        results[4] = retWave[3]/retWave[1]                              //avg# times scattered
177        results[5] = retWave[4]/retWave[1]                                              //single coherent fraction
178        results[6] = retWave[5]/retWave[1]                              //double coherent fraction
179        results[7] = retWave[6]/retWave[1]                              //multiple scatter fraction
180        results[8] = (retWave[0]-retWave[1])/retWave[0]                 //transmitted fraction
181       
182        return(0)
183End
184
185// worker function for threads, does nothing except switch between XOP and Igor versions
186//
187// uses ran3
188ThreadSafe Function Monte_SANS_W1(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
189        WAVE inputWave,ran_dev,nt,j1,j2,nn,linear_data,results
190       
191#if exists("Monte_SANSX")
192        Monte_SANSX(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
193#else
194        Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
195#endif
196
197        return (0)
198End
199
200// worker function for threads, does nothing except switch between XOP and Igor versions
201//
202// uses ran1
203ThreadSafe Function Monte_SANS_W2(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
204        WAVE inputWave,ran_dev,nt,j1,j2,nn,linear_data,results
205       
206#if exists("Monte_SANSX2")
207        Monte_SANSX2(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
208#else
209        Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
210#endif
211
212        return (0)
213End
214
215// NON-threaded call to the main function returns what is to be displayed
216// results is calculated and sent back for display
217Function Monte_SANS_NotThreaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
218        WAVE inputWave,ran_dev,nt,j1,j2,nn,linear_data,results
219
220        //initialize ran1 in the XOP by passing a negative integer
221        // does nothing in the Igor code, enoise is already initialized
222        Duplicate/O results retWave
223        WAVE retWave
224        retWave[0] = -1*abs(trunc(100000*enoise(1)))
225       
226#if exists("Monte_SANSX")
227        Monte_SANSX(inputWave,ran_dev,nt,j1,j2,nn,linear_data,retWave)
228#else
229        Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,retWave)
230#endif
231
232        // fill up the results wave
233        Variable xc,yc
234        xc=inputWave[3]
235        yc=inputWave[4]
236        results[0] = inputWave[9]+inputWave[10]         //total XS
237        results[1] = inputWave[10]                                              //SAS XS
238        results[2] = retWave[1]                                                 //number that interact n2
239        results[3] = retWave[2] - linear_data[xc][yc]                           //# reaching detector minus Q(0)
240        results[4] = retWave[3]/retWave[1]                              //avg# times scattered
241        results[5] = retWave[4]/retWave[1]                                              //single coherent fraction
242        results[6] = retWave[5]/retWave[1]                              //double coherent fraction
243        results[7] = retWave[6]/retWave[1]                              //multiple scatter fraction
244        results[8] = (retWave[0]-retWave[1])/retWave[0]                 //transmitted fraction
245       
246        return(0)
247End
248
249
250
251//////////
252//    PROGRAM Monte_SANS
253//    PROGRAM simulates multiple SANS.
254//       revised 2/12/99  JGB
255//            added calculation of random deviate, and 2D 10/2008 SRK
256
257//    N1 = NUMBER OF INCIDENT NEUTRONS.
258//    N2 = NUMBER INTERACTED IN THE SAMPLE.
259//    N3 = NUMBER ABSORBED.
260//    THETA = SCATTERING ANGLE.
261
262//        IMON = 'Enter number of neutrons to use in simulation.'
263//        NUM_BINS = 'Enter number of THETA BINS TO use. (<5000).'
264//        R1 = 'Enter beam radius. (cm)'
265//        R2 = 'Enter sample radius. (cm)'
266//        thick = 'Enter sample thickness. (cm)'
267//        wavelength = 'Enter neutron wavelength. (A)'
268//        R0 = 'Enter sphere radius. (A)'
269//
270
271ThreadSafe Function Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,MC_linear_data,results)
272        WAVE inputWave,ran_dev,nt,j1,j2,nn,MC_linear_data,results
273
274        Variable imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,sig_sas
275        Variable NUM_BINS,N_INDEX
276        Variable RHO,SIGSAS,SIGABS_0
277        Variable ii,jj,IND,idum,INDEX,IR,NQ
278        Variable qmax,theta_max,R_DAB,R0,BOUND,I0,q0,zpow
279        Variable N1,N2,N3,DTH,zz,tt,SIG_SINGLE,xx,yy,PHI,UU,SIG
280        Variable THETA,Ran,ll,D_OMEGA,RR,Tabs,Ttot,I1_sumI
281        Variable G0,E_NT,E_NN,TRANS_th,Trans_exp,rat
282        Variable GG,GG_ED,dS_dW,ds_dw_double,ds_dw_single
283        Variable DONE,FIND_THETA,err            //used as logicals
284
285        Variable Vx,Vy,Vz,Theta_z,qq
286        Variable Sig_scat,Sig_abs,Ratio,Sig_total
287        Variable isOn=0,testQ,testPhi,xPixel,yPixel
288        Variable NSingleIncoherent,NSingleCoherent,NScatterEvents,incoherentEvent,coherentEvent
289        Variable NDoubleCoherent,NMultipleScatter,countIt,detEfficiency
290       
291        detEfficiency = 1.0             //70% counting efficiency = 0.7
292       
293        imon = inputWave[0]
294        r1 = inputWave[1]
295        r2 = inputWave[2]
296        xCtr = inputWave[3]
297        yCtr = inputWave[4]
298        sdd = inputWave[5]
299        pixSize = inputWave[6]
300        thick = inputWave[7]
301        wavelength = inputWave[8]
302        sig_incoh = inputWave[9]
303        sig_sas = inputWave[10]
304       
305//      SetRandomSeed 0.1               //to get a reproduceable sequence
306
307//scattering power and maximum qvalue to bin
308//      zpow = .1               //scattering power, calculated below
309        qmax = 4*pi/wavelength          //maximum Q to bin 1D data. (A-1) (not really used, so set to a big value)
310        sigabs_0 = 0.0          // ignore absorption cross section/wavelength [1/(cm A)]
311        N_INDEX = 50            // maximum number of scattering events per neutron
312        num_bins = 200          //number of 1-D bins (not really used)
313       
314// my additions - calculate the random deviate function as needed
315// and calculate the scattering power from the model function (passed in as a wave)
316//
317        Variable left = leftx(ran_dev)
318        Variable delta = deltax(ran_dev)
319       
320//c       total SAS cross-section
321//      SIG_SAS = zpow/thick
322        zpow = sig_sas*thick                    //since I now calculate the sig_sas from the model
323        SIG_ABS = SIGABS_0 * WAVElength
324        sig_abs = 0.0           //cm-1
325        SIG_TOTAL =SIG_ABS + SIG_SAS + sig_incoh
326//      Print "The TOTAL XSECTION. (CM-1) is ",sig_total
327//      Print "The TOTAL SAS XSECTION. (CM-1) is ",sig_sas
328//      results[0] = sig_total          //assign these after everything's done
329//      results[1] = sig_sas
330//      variable ratio1,ratio2
331//      ratio1 = sig_abs/sig_total
332//      ratio2 = sig_incoh/sig_total
333//      // 0->ratio1 = abs
334//      // ratio1 -> ratio2 = incoh
335//      // > ratio2 = coherent
336        RATIO = sig_incoh / SIG_TOTAL
337       
338//c       assuming theta = sin(theta)...OK
339        theta_max = wavelength*qmax/(2*pi)
340//C     SET Theta-STEP SIZE.
341        DTH = Theta_max/NUM_BINS
342//      Print "theta bin size = dth = ",dth
343
344//C     INITIALIZE COUNTERS.
345        N1 = 0
346        N2 = 0
347        N3 = 0
348        NSingleIncoherent = 0
349        NSingleCoherent = 0
350        NDoubleCoherent = 0
351        NMultipleScatter = 0
352        NScatterEvents = 0
353
354//C     INITIALIZE ARRAYS.
355        j1 = 0
356        j2 = 0
357        nt = 0
358        nn=0
359       
360//C     MONITOR LOOP - looping over the number of incedent neutrons
361//note that zz, is the z-position in the sample - NOT the scattering power
362
363// NOW, start the loop, throwing neutrons at the sample.
364        do
365                Vx = 0.0                        // Initialize direction vector.
366                Vy = 0.0
367                Vz = 1.0
368               
369                Theta = 0.0             //      Initialize scattering angle.
370                Phi = 0.0                       //      Intialize azimuthal angle.
371                N1 += 1                 //      Increment total number neutrons counter.
372                DONE = 0                        //      True when neutron is scattered out of the sample.
373                INDEX = 0                       //      Set counter for number of scattering events.
374                zz = 0.0                        //      Set entering dimension of sample.
375                incoherentEvent = 0
376                coherentEvent = 0
377               
378                do                                      //      Makes sure position is within circle.
379                        ran = abs(enoise(1))            //[0,1]
380                        xx = 2.0*R1*(Ran-0.5)           //X beam position of neutron entering sample.
381                        ran = abs(enoise(1))            //[0,1]
382                        yy = 2.0*R1*(Ran-0.5)           //Y beam position ...
383                        RR = SQRT(xx*xx+yy*yy)          //Radial position of neutron in incident beam.
384                while(rr>r1)
385
386                do    //Scattering Loop, will exit when "done" == 1
387                                // keep scattering multiple times until the neutron exits the sample
388                        ran = abs(enoise(1))            //[0,1]  RANDOM NUMBER FOR DETERMINING PATH LENGTH
389                        ll = PATH_len(ran,Sig_total)
390                        //Determine new scattering direction vector.
391                        err = NewDirection(vx,vy,vz,Theta,Phi)          //vx,vy,vz is updated, theta, phi unchanged by function
392
393                        //X,Y,Z-POSITION OF SCATTERING EVENT.
394                        xx += ll*vx
395                        yy += ll*vy
396                        zz += ll*vz
397                        RR = sqrt(xx*xx+yy*yy)          //radial position of scattering event.
398
399                        //Check whether interaction occurred within sample volume.
400                        IF (((zz > 0.0) && (zz < THICK)) && (rr < r2))
401                                //NEUTRON INTERACTED.
402                                ran = abs(enoise(1))            //[0,1]
403                               
404//                              if(ran<ratio1)
405//                                      //absorption event
406//                                      n3 +=1
407//                                      done=1
408//                              else
409
410                                INDEX += 1                      //Increment counter of scattering events.
411                                IF(INDEX == 1)
412                                        N2 += 1                 //Increment # of scat. neutrons
413                                Endif
414                                //Split neutron interactions into scattering and absorption events
415//                              IF(ran > (ratio1+ratio2) )              //C             NEUTRON SCATTERED coherently
416                                IF(ran > ratio)         //C             NEUTRON SCATTERED coherently
417                                        coherentEvent = 1
418                                        FIND_THETA = 0                  //false
419                                        DO
420                                                //ran = abs(enoise(1))          //[0,1]
421                                                //theta = Scat_angle(Ran,R_DAB,wavelength)      // CHOOSE DAB ANGLE -- this is 2Theta
422                                                //Q0 = 2*PI*THETA/WAVElength                                    // John chose theta, calculated Q
423
424                                                // pick a q-value from the deviate function
425                                                // pnt2x truncates the point to an integer before returning the x
426                                                // so get it from the wave scaling instead
427                                                Q0 =left + binarysearchinterp(ran_dev,abs(enoise(1)))*delta
428                                                theta = Q0/2/Pi*wavelength              //SAS approximation. 1% error at theta=30 deg (theta/2=15deg)
429                                               
430                                                //Print "q0, theta = ",q0,theta
431                                               
432                                                FIND_THETA = 1          //always accept
433
434                                        while(!find_theta)
435                                        ran = abs(enoise(1))            //[0,1]
436                                        PHI = 2.0*PI*Ran                        //Chooses azimuthal scattering angle.
437                                ELSE
438                                        //NEUTRON scattered incoherently
439                   // N3 += 1
440                  // DONE = 1
441                  // phi and theta are random over the entire sphere of scattering
442                  // !can't just choose random theta and phi, won't be random over sphere solid angle
443                        incoherentEvent = 1
444                       
445                        ran = abs(enoise(1))            //[0,1]
446                                        theta = acos(2*ran-1)           
447                       
448                        ran = abs(enoise(1))            //[0,1]
449                                        PHI = 2.0*PI*Ran                        //Chooses azimuthal scattering angle.
450                                ENDIF           //(ran > ratio)
451//                              endif           // event was absorption
452                        ELSE
453                                //NEUTRON ESCAPES FROM SAMPLE -- bin it somewhere
454                                DONE = 1                //done = true, will exit from loop
455                               
456//                              countIt = 1
457//                              if(abs(enoise(1)) > detEfficiency)              //efficiency of 70% wired @top
458//                                      countIt = 0                                     //detector does not register
459//                              endif
460                               
461                                //Increment #scattering events array
462                                If (index <= N_Index)
463                                        NN[INDEX] += 1
464                                Endif
465                               
466                                if(index != 0)          //the neutron interacted at least once, figure out where it ends up
467
468                                        Theta_z = acos(Vz)              // Angle WITH respect to z axis.
469                                        testQ = 2*pi*sin(theta_z)/wavelength
470                                       
471                                        // pick a random phi angle, and see if it lands on the detector
472                                        // since the scattering is isotropic, I can safely pick a new, random value
473                                        // this would not be true if simulating anisotropic scattering.
474                                        //testPhi = abs(enoise(1))*2*Pi
475                                        testPhi = MC_FindPhi(Vx,Vy)             //use the exiting phi value as defined by Vx and Vy
476                                       
477                                        // is it on the detector?       
478                                        FindPixel(testQ,testPhi,wavelength,sdd,pixSize,xCtr,yCtr,xPixel,yPixel)
479                                       
480                                        if(xPixel != -1 && yPixel != -1)
481                                                //if(index==1)  // only the single scattering events
482                                                        MC_linear_data[xPixel][yPixel] += 1             //this is the total scattering, including multiple scattering
483                                                //endif
484                                                        isOn += 1               // neutron that lands on detector
485                                        endif
486       
487                                        If(theta_z < theta_max)
488                                                //Choose index for scattering angle array.
489                                                //IND = NINT(THETA_z/DTH + 0.4999999)
490                                                ind = round(THETA_z/DTH + 0.4999999)            //round is eqivalent to nint()
491                                                NT[ind] += 1                    //Increment bin for angle.
492                                                //Increment angle array for single scattering events.
493                                                IF(INDEX == 1)
494                                                        j1[ind] += 1
495                                                Endif
496                                                //Increment angle array for double scattering events.
497                                                IF (INDEX == 2)
498                                                        j2[ind] += 1
499                                                Endif
500                                        EndIf
501                                       
502                                        // increment all of the counters now since done==1 here and I'm sure to exit and get another neutron
503                                        NScatterEvents += index         //total number of scattering events
504                                        if(index == 1 && incoherentEvent == 1)
505                                                NSingleIncoherent += 1
506                                        endif
507                                        if(index == 1 && coherentEvent == 1)
508                                                NSingleCoherent += 1
509                                        endif
510                                        if(index == 2 && coherentEvent == 1 && incoherentEvent == 0)
511                                                NDoubleCoherent += 1
512                                        endif
513                                        if(index > 1)
514                                                NMultipleScatter += 1
515                                        endif
516                                        //Print "n1,index (x,y) = ",n1,index, xpixel,ypixel
517                                       
518                                else    // if neutron escaped without interacting
519                               
520                                        // then it must be a transmitted neutron
521                                        // don't need to calculate, just increment the proper counters
522                                       
523                                        MC_linear_data[xCtr+xx/pixsize][yCtr+yy/pixsize] += 1
524                                        isOn += 1
525                                        nt[0] += 1
526                                       
527                                endif           //if interacted
528                        ENDIF
529                while (!done)
530        while(n1 < imon)
531
532//      Print "Monte Carlo Done"
533        results[0] = n1
534        results[1] = n2
535        results[2] = isOn
536        results[3] = NScatterEvents             //sum of # of times that neutrons scattered
537        results[4] = NSingleCoherent            //# of events that are single, coherent
538        results[5] = NDoubleCoherent
539        results[6] = NMultipleScatter           //# of multiple scattering events
540       
541//      Print "# absorbed = ",n3
542
543//      trans_th = exp(-sig_total*thick)
544//      TRANS_exp = (N1-N2) / N1                        // Transmission
545        // dsigma/domega assuming isotropic scattering, with no absorption.
546//      Print "trans_exp = ",trans_exp
547//      Print "total # of neutrons reaching 2D detector",isOn
548//      Print "fraction of incident neutrons reaching detector ",isOn/iMon
549       
550//      Print "Total number of neutrons = ",N1
551//      Print "Total number of neutrons that interact = ",N2
552//      Print "Fraction of singly scattered neutrons = ",sum(j1,-inf,inf)/N2
553//      results[2] = N2                                         //number that scatter
554//      results[3] = isOn - MC_linear_data[xCtr][yCtr]                  //# scattered reaching detector minus zero angle
555
556       
557//      Tabs = (N1-N3)/N1
558//      Ttot = (N1-N2)/N1
559//      I1_sumI = NN[0]/(N2-N3)
560//      Print "Tabs = ",Tabs
561//      Print "Transmitted neutrons = ",Ttot
562//      results[8] = Ttot
563//      Print "I1 / all I1 = ", I1_sumI
564
565End
566////////        end of main function for calculating multiple scattering
567
568
569// returns the random deviate as a wave
570// and the total SAS cross-section [1/cm] sig_sas
571Function CalculateRandomDeviate(func,coef,lam,outWave,SASxs)
572        FUNCREF SANSModelAAO_MCproto func
573        WAVE coef
574        Variable lam
575        String outWave
576        Variable &SASxs
577
578        Variable nPts_ran=10000,qu
579        qu = 4*pi/lam           
580       
581// hard-wired into the Simulation directory rather than the SAS folder.
582// plotting resolution-smeared models won't work any other way
583        Make/O/N=(nPts_ran)/D root:Simulation:Gq,root:Simulation:xw             // if these waves are 1000 pts, the results are "pixelated"
584        WAVE Gq = root:Simulation:gQ
585        WAVE xw = root:Simulation:xw
586        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
587
588/// if all of the coefficients are well-behaved, then the last point is the background
589// and I can set it to zero here (only for the calculation)
590        Duplicate/O coef,tmp_coef
591        Variable num=numpnts(coef)
592        tmp_coef[num-1] = 0
593       
594        xw=x                                                                                            //for the AAO
595        func(tmp_coef,Gq,xw)                                                                    //call as AAO
596
597//      Gq = x*Gq                                                                                                       // SAS approximation
598        Gq = Gq*sin(2*asin(x/qu))/sqrt(1-(x/qu))                        // exact
599        //
600        //
601        Integrate/METH=1 Gq/D=Gq_INT
602       
603//      SASxs = lam*lam/2/pi*Gq_INT[nPts_ran-1]                 //if the approximation is used
604        SASxs = lam*Gq_INT[nPts_ran-1]
605       
606        Gq_INT /= Gq_INT[nPts_ran-1]
607       
608        Duplicate/O Gq_INT $outWave
609
610        return(0)
611End
612
613// returns the random deviate as a wave
614// and the total SAS cross-section [1/cm] sig_sas
615//
616// uses a log spacing of x for better coverage
617// downside is that it doesn't use built-in integrate, which is automatically cumulative
618//
619// --- Currently does not work - since the usage of the random deviate in the MC routine is based on the
620// wave properties of ran_dev, that is it must have the proper scaling and be equally spaced.
621//
622// -- not really sure that this will solve any of the problems with some functions (notably those with power-laws)
623// giving unreasonably large SAS cross sections. (>>10)
624//
625Function CalculateRandomDeviate_log(func,coef,lam,outWave,SASxs)
626        FUNCREF SANSModelAAO_MCproto func
627        WAVE coef
628        Variable lam
629        String outWave
630        Variable &SASxs
631
632        Variable nPts_ran=1000,qu,qmin,ii
633        qmin=1e-5
634        qu = 4*pi/lam           
635
636// hard-wired into the Simulation directory rather than the SAS folder.
637// plotting resolution-smeared models won't work any other way
638        Make/O/N=(nPts_ran)/D root:Simulation:Gq,root:Simulation:xw             // if these waves are 1000 pts, the results are "pixelated"
639        WAVE Gq = root:Simulation:gQ
640        WAVE xw = root:Simulation:xw
641//      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
642        xw =  alog(log(qmin) + x*((log(qu)-log(qmin))/nPts_ran))
643
644/// if all of the coefficients are well-behaved, then the last point is the background
645// and I can set it to zero here (only for the calculation)
646        Duplicate/O coef,tmp_coef
647        Variable num=numpnts(coef)
648        tmp_coef[num-1] = 0
649       
650        func(tmp_coef,Gq,xw)                                                                    //call as AAO
651        Gq = Gq*sin(2*asin(xw/qu))/sqrt(1-(xw/qu))                      // exact
652
653       
654        Duplicate/O Gq Gq_INT
655        Gq_INT = 0
656        for(ii=0;ii<nPts_ran;ii+=1)
657                Gq_INT[ii] = AreaXY(xw,Gq,qmin,xw[ii])
658        endfor
659       
660        SASxs = lam*Gq_INT[nPts_ran-1]
661       
662        Gq_INT /= Gq_INT[nPts_ran-1]
663       
664        Duplicate/O Gq_INT $outWave
665
666        return(0)
667End
668
669ThreadSafe Function FindPixel(testQ,testPhi,lam,sdd,pixSize,xCtr,yCtr,xPixel,yPixel)
670        Variable testQ,testPhi,lam,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel
671
672        Variable theta,dy,dx,qx,qy
673        //decompose to qx,qy
674        qx = testQ*cos(testPhi)
675        qy = testQ*sin(testPhi)
676
677        //convert qx,qy to pixel locations relative to # of pixels x, y from center
678        theta = 2*asin(qy*lam/4/pi)
679        dy = sdd*tan(theta)
680        yPixel = round(yCtr + dy/pixSize)
681       
682        theta = 2*asin(qx*lam/4/pi)
683        dx = sdd*tan(theta)
684        xPixel = round(xCtr + dx/pixSize)
685
686        //if on detector, return xPix and yPix values, otherwise -1
687        if(yPixel > 127 || yPixel < 0)
688                yPixel = -1
689        endif
690        if(xPixel > 127 || xPixel < 0)
691                xPixel = -1
692        endif
693       
694        return(0)
695End
696
697Function MC_CheckFunctionAndCoef(funcStr,coefStr)
698        String funcStr,coefStr
699       
700        SVAR/Z listStr=root:Packages:NIST:coefKWStr
701        if(SVAR_Exists(listStr) == 1)
702                String properCoefStr = StringByKey(funcStr, listStr  ,"=",";",0)
703                if(cmpstr("",properCoefStr)==0)
704                        return(0)               //false, no match found, so properCoefStr is returned null
705                endif
706                if(cmpstr(coefStr,properCoefStr)==0)
707                        return(1)               //true, the coef is the correct match
708                endif
709        endif
710        return(0)                       //false, wrong coef
711End
712
713Function/S MC_getFunctionCoef(funcStr)
714        String funcStr
715
716        SVAR/Z listStr=root:Packages:NIST:coefKWStr
717        String coefStr=""
718        if(SVAR_Exists(listStr) == 1)
719                coefStr = StringByKey(funcStr, listStr  ,"=",";",0)
720        endif
721        return(coefStr)
722End
723
724Function SANSModelAAO_MCproto(w,yw,xw)
725        Wave w,yw,xw
726       
727        Print "in SANSModelAAO_MCproto function"
728        return(1)
729end
730
731Function/S MC_FunctionPopupList()
732        String list,tmp
733        list = User_FunctionPopupList()
734       
735        //simplify the display, forcing smeared calculations behind the scenes
736        tmp = FunctionList("Smear*",";","NPARAMS:1")            //smeared dependency calculations
737        list = RemoveFromList(tmp, list,";")
738
739
740        if(strlen(list)==0)
741                list = "No functions plotted"
742        endif
743       
744        list = SortList(list)
745       
746        return(list)
747End             
748
749
750//Function Scat_Angle(Ran,R_DAB,wavelength)
751//      Variable Ran,r_dab,wavelength
752//
753//      Variable qq,arg,theta
754//      qq = 1. / ( R_DAB*sqrt(1.0/Ran - 1.0) )
755//      arg = qq*wavelength/(4*pi)
756//      If (arg < 1.0)
757//              theta = 2.*asin(arg)
758//      else
759//              theta = pi
760//      endif
761//      Return (theta)
762//End
763
764//calculates new direction (xyz) from an old direction
765//theta and phi don't change
766ThreadSafe Function NewDirection(vx,vy,vz,theta,phi)
767        Variable &vx,&vy,&vz
768        Variable theta,phi
769       
770        Variable err=0,vx0,vy0,vz0
771        Variable nx,ny,mag_xy,tx,ty,tz
772       
773        //store old direction vector
774        vx0 = vx
775        vy0 = vy
776        vz0 = vz
777       
778        mag_xy = sqrt(vx0*vx0 + vy0*vy0)
779        if(mag_xy < 1e-12)
780                //old vector lies along beam direction
781                nx = 0
782                ny = 1
783                tx = 1
784                ty = 0
785                tz = 0
786        else
787                Nx = -Vy0 / Mag_XY
788                Ny = Vx0 / Mag_XY
789                Tx = -Vz0*Vx0 / Mag_XY
790                Ty = -Vz0*Vy0 / Mag_XY
791                Tz = Mag_XY
792        endif
793       
794        //new scattered direction vector
795        Vx = cos(phi)*sin(theta)*Tx + sin(phi)*sin(theta)*Nx + cos(theta)*Vx0
796        Vy = cos(phi)*sin(theta)*Ty + sin(phi)*sin(theta)*Ny + cos(theta)*Vy0
797        Vz = cos(phi)*sin(theta)*Tz + cos(theta)*Vz0
798       
799        Return(err)
800End
801
802ThreadSafe Function path_len(aval,sig_tot)
803        Variable aval,sig_tot
804       
805        Variable retval
806       
807        retval = -1*ln(1-aval)/sig_tot
808       
809        return(retval)
810End
811
812// globals are initialized in SASCALC.ipf
813// coordinates if I make this part of the panel - but this breaks other things...
814//
815//Proc MC_SASCALC()
816////    PauseUpdate; Silent 1           // building window...
817//
818////    NewPanel /W=(92,556,390,1028)/K=1 as "SANS Simulator"
819//      SetVariable MC_setvar0,pos={491,73},size={144,15},bodyWidth=80,title="# of neutrons"
820//      SetVariable MC_setvar0,format="%5.4g"
821//      SetVariable MC_setvar0,limits={-inf,inf,100},value= root:Packages:NIST:SAS:gImon
822//      SetVariable MC_setvar0_1,pos={491,119},size={131,15},bodyWidth=60,title="Thickness (cm)"
823//      SetVariable MC_setvar0_1,limits={-inf,inf,0.1},value= root:Packages:NIST:SAS:gThick
824//      SetVariable MC_setvar0_2,pos={491,96},size={149,15},bodyWidth=60,title="Incoherent XS (cm)"
825//      SetVariable MC_setvar0_2,limits={-inf,inf,0.1},value= root:Packages:NIST:SAS:gSig_incoh
826//      SetVariable MC_setvar0_3,pos={491,142},size={150,15},bodyWidth=60,title="Sample Radius (cm)"
827//      SetVariable MC_setvar0_3,limits={-inf,inf,0.1},value= root:Packages:NIST:SAS:gR2
828//      PopupMenu MC_popup0,pos={476,13},size={165,20},proc=MC_ModelPopMenuProc,title="Model Function"
829//      PopupMenu MC_popup0,mode=1,value= #"MC_FunctionPopupList()"
830//      Button MC_button0,pos={480,181},size={130,20},proc=MC_DoItButtonProc,title="Do MC Simulation"
831//      Button MC_button1,pos={644,181},size={80,20},proc=MC_Display2DButtonProc,title="Show 2D"
832//      SetVariable setvar0_3,pos={568,484},size={50,20},disable=1
833//      GroupBox group0,pos={478,42},size={267,130},title="Monte Carlo"
834//      SetVariable cntVar,pos={653,73},size={80,15},proc=CountTimeSetVarProc,title="time(s)"
835//      SetVariable cntVar,format="%d"
836//      SetVariable cntVar,limits={1,10,1},value= root:Packages:NIST:SAS:gCntTime
837//     
838//      String fldrSav0= GetDataFolder(1)
839//      SetDataFolder root:Packages:NIST:SAS:
840//      Edit/W=(476,217,746,450)/HOST=#  results_desc,results
841//      ModifyTable format(Point)=1,width(Point)=0,width(results_desc)=150
842//      SetDataFolder fldrSav0
843//      RenameWindow #,T_results
844//      SetActiveSubwindow ##
845EndMacro
846
847// as a stand-alone panel, extra control bar  (right) and subwindow implementations don't work right
848// for various reasons...
849Window MC_SASCALC() : Panel
850        PauseUpdate; Silent 1           // building window...
851        NewPanel /W=(92,556,713,818)/K=1 as "SANS Simulator"
852        SetVariable MC_setvar0,pos={28,73},size={144,15},bodyWidth=80,title="# of neutrons"
853        SetVariable MC_setvar0,format="%5.4g"
854        SetVariable MC_setvar0,limits={0,inf,100},value= root:Packages:NIST:SAS:gImon
855        SetVariable MC_setvar0_1,pos={28,119},size={131,15},bodyWidth=60,title="Thickness (cm)"
856        SetVariable MC_setvar0_1,limits={0,inf,0.1},value= root:Packages:NIST:SAS:gThick
857        SetVariable MC_setvar0_2,pos={28,96},size={149,15},bodyWidth=60,title="Incoherent XS (cm)"
858        SetVariable MC_setvar0_2,limits={0,inf,0.1},value= root:Packages:NIST:SAS:gSig_incoh
859        SetVariable MC_setvar0_3,pos={28,142},size={150,15},bodyWidth=60,title="Sample Radius (cm)"
860        SetVariable MC_setvar0_3,limits={-inf,inf,0.1},value= root:Packages:NIST:SAS:gR2
861        PopupMenu MC_popup0,pos={13,13},size={165,20},proc=MC_ModelPopMenuProc,title="Model Function"
862        PopupMenu MC_popup0,mode=1,value= #"MC_FunctionPopupList()"
863        Button MC_button0,pos={17,181},size={130,20},proc=MC_DoItButtonProc,title="Do MC Simulation"
864        Button MC_button0,fColor=(3,52428,1)
865        Button MC_button1,pos={17,208},size={80,20},proc=MC_Display2DButtonProc,title="Show 2D"
866        SetVariable setvar0_3,pos={105,484},size={50,20},disable=1
867        GroupBox group0,pos={15,42},size={267,130},title="Monte Carlo"
868        SetVariable cntVar,pos={190,73},size={80,15},proc=CountTimeSetVarProc,title="time(s)"
869        SetVariable cntVar,format="%d"
870        SetVariable cntVar,limits={1,600,1},value= root:Packages:NIST:SAS:gCntTime
871        Button MC_button2,pos={17,234},size={100,20},proc=SaveAsVAXButtonProc,title="Save 2D VAX"
872        CheckBox check0,pos={216,180},size={68,14},title="Raw counts",variable = root:Packages:NIST:SAS:gRawCounts
873        CheckBox check0_1,pos={216,199},size={60,14},title="Yes Offset",variable= root:Packages:NIST:SAS:gDoTraceOffset
874        CheckBox check0_2,pos={216,199+19},size={60,14},title="Beam Stop in",variable= root:Packages:NIST:SAS:gBeamStopIn
875        CheckBox check0_3,pos={216,199+2*19},size={60,14},title="use XOP",variable= root:Packages:NIST:SAS:gUse_MC_XOP
876       
877        String fldrSav0= GetDataFolder(1)
878        SetDataFolder root:Packages:NIST:SAS:
879        Edit/W=(344,23,606,248)/HOST=#  results_desc,results
880        ModifyTable format(Point)=1,width(Point)=0,width(results_desc)=150
881        SetDataFolder fldrSav0
882        RenameWindow #,T_results
883        SetActiveSubwindow ##
884EndMacro
885
886Function CountTimeSetVarProc(sva) : SetVariableControl
887        STRUCT WMSetVariableAction &sva
888
889        switch( sva.eventCode )
890                case 1: // mouse up
891                case 2: // Enter key
892                case 3: // Live update
893                        Variable dval = sva.dval
894
895                        // get the neutron flux, multiply, and reset the global for # neutrons
896                        NVAR imon=root:Packages:NIST:SAS:gImon
897                        imon = dval*beamIntensity()
898                       
899                        break
900        endswitch
901
902        return 0
903End
904
905
906Function MC_ModelPopMenuProc(pa) : PopupMenuControl
907        STRUCT WMPopupAction &pa
908
909        switch( pa.eventCode )
910                case 2: // mouse up
911                        Variable popNum = pa.popNum
912                        String popStr = pa.popStr
913                        SVAR gStr = root:Packages:NIST:SAS:gFuncStr
914                        gStr = popStr
915                       
916                        break
917        endswitch
918
919        return 0
920End
921
922Function MC_DoItButtonProc(ba) : ButtonControl
923        STRUCT WMButtonAction &ba
924
925        switch( ba.eventCode )
926                case 2: // mouse up
927                        // click code here
928                        NVAR doMC = root:Packages:NIST:SAS:gDoMonteCarlo
929                        doMC = 1
930                        ReCalculateInten(1)
931                        doMC = 0                //so the next time won't be MC
932                        break
933        endswitch
934
935        return 0
936End
937
938
939Function MC_Display2DButtonProc(ba) : ButtonControl
940        STRUCT WMButtonAction &ba
941
942        switch( ba.eventCode )
943                case 2: // mouse up
944                        // click code here
945                        Execute "ChangeDisplay(\"SAS\")"
946                        break
947        endswitch
948
949        return 0
950End
951
952// after a 2d data image is averaged in the usual way, take the waves and generate a "fake" folder of the 1d
953// data, to appear as if it was loaded from a real data file.
954//
955// currently only works with SANS data, but can later be expanded to generate fake USANS data sets
956// ---- use FakeUSANSDataFolder() instead----
957//
958Function        Fake1DDataFolder(qval,aveint,sigave,sigmaQ,qbar,fSubs,dataFolder)
959        WAVE qval,aveint,sigave,sigmaQ,qbar,fSubs
960        String dataFolder
961
962        String baseStr=dataFolder
963        if(DataFolderExists("root:"+baseStr))
964                SetDataFolder $("root:"+baseStr)
965        else
966                NewDataFolder/S $("root:"+baseStr)
967        endif
968
969        ////overwrite the existing data, if it exists
970        Duplicate/O qval, $(baseStr+"_q")
971        Duplicate/O aveint, $(baseStr+"_i")
972        Duplicate/O sigave, $(baseStr+"_s")
973//      Duplicate/O sigmaQ, $(baseStr+"sq")
974//      Duplicate/O qbar, $(baseStr+"qb")
975//      Duplicate/O fSubS, $(baseStr+"fs")
976
977        // need to switch based on SANS/USANS
978        if (isSANSResolution(sigave[0]))                //checks to see if the first point of the wave is <0]
979                // make a resolution matrix for SANS data
980                Variable np=numpnts(qval)
981                Make/D/O/N=(np,4) $(baseStr+"_res")
982                Wave res=$(baseStr+"_res")
983               
984                res[][0] = sigmaQ[p]            //sigQ
985                res[][1] = qBar[p]              //qBar
986                res[][2] = fSubS[p]             //fShad
987                res[][3] = qval[p]              //Qvalues
988               
989                // keep a copy of everything in SAS too... the smearing wrapper function looks for
990                // data in folders based on waves it is passed - an I lose control of that
991                Duplicate/O res, $("root:Packages:NIST:SAS:"+baseStr+"_res")
992                Duplicate/O qval,  $("root:Packages:NIST:SAS:"+baseStr+"_q")
993                Duplicate/O aveint,  $("root:Packages:NIST:SAS:"+baseStr+"_i")
994                Duplicate/O sigave,  $("root:Packages:NIST:SAS:"+baseStr+"_s")
995        else
996                //the data is USANS data
997                // nothing done here yet
998//              dQv = -$w3[0]
999               
1000//              USANS_CalcWeights(baseStr,dQv)
1001               
1002        endif
1003
1004        //clean up             
1005        SetDataFolder root:
1006       
1007End
1008
1009// writes out a VAX binary data file
1010// automatically generates a name
1011// will prompt for the sample label
1012//
1013Function SaveAsVAXButtonProc(ctrlName) : ButtonControl
1014        String ctrlName
1015
1016        WriteVAXData("SAS","",0)
1017End
1018
1019// calculates the fraction of the scattering that reaches the detector, given the random deviate function
1020// and qmin and qmax
1021//
1022//
1023// still some question of the corners and number of pixels per q-bin
1024Function FractionReachingDetector(ran_dev,Qmin,Qmax)
1025        wave ran_dev
1026        Variable Qmin,Qmax
1027       
1028        Variable r1,r2,frac
1029        r1=x2pnt(ran_dev,Qmin)
1030        r2=x2pnt(ran_dev,Qmax)
1031       
1032        // no normalization needed - the full q-range is defined as [0,1]
1033        frac = ran_dev[r2] - ran_dev[r1]
1034       
1035        return frac
1036End
1037
1038
1039/// called in SASCALC:ReCalculateInten()
1040Function        Simulate_2D_MC(funcStr,aveint,qval,sigave,sigmaq,qbar,fsubs)
1041        String funcStr
1042        WAVE aveint,qval,sigave,sigmaq,qbar,fsubs
1043
1044        NVAR doMonteCarlo = root:Packages:NIST:SAS:gDoMonteCarlo                // == 1 if 2D MonteCarlo set by hidden flag
1045        WAVE rw=root:Packages:NIST:SAS:realsRead
1046
1047        NVAR imon = root:Packages:NIST:SAS:gImon
1048        NVAR thick = root:Packages:NIST:SAS:gThick
1049        NVAR sig_incoh = root:Packages:NIST:SAS:gSig_incoh
1050        NVAR r2 = root:Packages:NIST:SAS:gR2
1051
1052        // do the simulation here, or not
1053        Variable r1,xCtr,yCtr,sdd,pixSize,wavelength
1054        String coefStr,abortStr,str
1055
1056        r1 = rw[24]/2/10                // sample diameter convert diam in [mm] to radius in cm
1057        xCtr = rw[16]
1058        yCtr = rw[17]
1059        sdd = rw[18]*100                //conver header of [m] to [cm]
1060        pixSize = rw[10]/10             // convert pix size in mm to cm
1061        wavelength = rw[26]
1062        coefStr = MC_getFunctionCoef(funcStr)
1063       
1064        if(!MC_CheckFunctionAndCoef(funcStr,coefStr))
1065                doMonteCarlo = 0                //we're getting out now, reset the flag so we don't continually end up here
1066                Abort "The coefficients and function type do not match. Please correct the selections in the popup menus."
1067        endif
1068       
1069        Variable sig_sas
1070//              FUNCREF SANSModelAAO_MCproto func=$("fSmeared"+funcStr)         //a wrapper for the structure version
1071        FUNCREF SANSModelAAO_MCproto func=$(funcStr)            //unsmeared
1072        WAVE results = root:Packages:NIST:SAS:results
1073        WAVE linear_data = root:Packages:NIST:SAS:linear_data
1074        WAVE data = root:Packages:NIST:SAS:data
1075
1076        results = 0
1077        linear_data = 0
1078       
1079        CalculateRandomDeviate(func,$coefStr,wavelength,"root:Packages:NIST:SAS:ran_dev",SIG_SAS)
1080        if(sig_sas > 100)
1081                sprintf abortStr,"sig_sas = %g. Please check that the model coefficients have a zero background, or the low q is well-behaved.",sig_sas
1082                Abort abortStr
1083        endif
1084       
1085        WAVE ran_dev=$"root:Packages:NIST:SAS:ran_dev"
1086       
1087        Make/O/D/N=5000 root:Packages:NIST:SAS:nt=0,root:Packages:NIST:SAS:j1=0,root:Packages:NIST:SAS:j2=0
1088        Make/O/D/N=100 root:Packages:NIST:SAS:nn=0
1089        Make/O/D/N=11 root:Packages:NIST:SAS:inputWave=0
1090       
1091        WAVE nt = root:Packages:NIST:SAS:nt
1092        WAVE j1 = root:Packages:NIST:SAS:j1
1093        WAVE j2 = root:Packages:NIST:SAS:j2
1094        WAVE nn = root:Packages:NIST:SAS:nn
1095        WAVE inputWave = root:Packages:NIST:SAS:inputWave
1096
1097        inputWave[0] = imon
1098        inputWave[1] = r1
1099        inputWave[2] = r2
1100        inputWave[3] = xCtr
1101        inputWave[4] = yCtr
1102        inputWave[5] = sdd
1103        inputWave[6] = pixSize
1104        inputWave[7] = thick
1105        inputWave[8] = wavelength
1106        inputWave[9] = sig_incoh
1107        inputWave[10] = sig_sas
1108
1109        linear_data = 0         //initialize
1110
1111        Variable t0,trans
1112       
1113        // get a time estimate, and give the user a chance to exit if they're unsure.
1114        t0 = stopMStimer(-2)
1115        inputWave[0] = 1000
1116        NVAR useXOP = root:Packages:NIST:SAS:gUse_MC_XOP                //if zero, will use non-threaded Igor code
1117       
1118        if(useXOP)
1119                //use a single thread, otherwise time is dominated by overhead
1120                Monte_SANS_NotThreaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
1121        else
1122                Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
1123        endif
1124       
1125        t0 = (stopMSTimer(-2) - t0)*1e-6
1126        t0 *= imon/1000/ThreadProcessorCount                    //projected time, in seconds (using threads for the calculation)
1127        inputWave[0] = imon             //reset
1128       
1129        if(t0>10)
1130                sprintf str,"The simulation will take approximately %d seconds.\r- Proceed?",t0
1131                DoAlert 1,str
1132                if(V_flag == 2)
1133                        doMonteCarlo = 0
1134                        reCalculateInten(1)             //come back in and do the smeared calculation
1135                        return(0)
1136                endif
1137        endif
1138       
1139        linear_data = 0         //initialize
1140// threading crashes!! - there must be some operation in the XOP that is not threadSafe. What, I don't know...
1141// I think it's the ran() calls, being "non-reentrant". So the XOP now defines two separate functions, that each
1142// use a different rng. This works. 1.75x speedup.     
1143        t0 = stopMStimer(-2)
1144
1145        if(useXOP)
1146                Monte_SANS_Threaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
1147        else
1148                Monte_SANS_NotThreaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
1149        endif
1150       
1151        t0 = (stopMSTimer(-2) - t0)*1e-6
1152        Printf  "MC sim time = %g seconds\r",t0
1153       
1154        trans = results[8]                      //(n1-n2)/n1
1155        if(trans == 0)
1156                trans = 1
1157        endif
1158
1159        Print "counts on detector, including transmitted = ",sum(linear_data,-inf,inf)
1160       
1161//              linear_data[xCtr][yCtr] = 0                     //snip out the transmitted spike
1162//              Print "counts on detector not transmitted = ",sum(linear_data,-inf,inf)
1163
1164        // or simulate a beamstop
1165        NVAR MC_BS_in = root:Packages:NIST:SAS:gBeamStopIn              //if zero, beam stop is "out", as in a transmission measurement
1166       
1167        Variable rad=beamstopDiam()/2           //beamstop radius in cm
1168        if(MC_BS_in)
1169                rad /= 0.5                              //convert cm to pixels
1170                rad += 0.                                       // (no - it cuts off the low Q artificially) add an extra pixel to each side to account for edge
1171                Duplicate/O linear_data,root:Packages:NIST:SAS:tmp_mask//,root:Packages:NIST:SAS:MC_linear_data
1172                WAVE tmp_mask = root:Packages:NIST:SAS:tmp_mask
1173                tmp_mask = (sqrt((p-xCtr)^2+(q-yCtr)^2) < rad) ? 0 : 1          //behind beamstop = 0, away = 1
1174               
1175                linear_data *= tmp_mask
1176        endif
1177       
1178        results[9] = sum(linear_data,-inf,inf)
1179        //              Print "counts on detector not behind beamstop = ",results[9]
1180       
1181        // convert to absolute scale
1182        Variable kappa          //,beaminten = beamIntensity()
1183//              kappa = beamInten*pi*r1*r1*thick*(pixSize/sdd)^2*trans*(iMon/beaminten)
1184        kappa = thick*(pixSize/sdd)^2*trans*iMon
1185       
1186        //use kappa to get back to counts => linear_data = round(linear_data*kappa)
1187        Note/K linear_data ,"KAPPA="+num2str(kappa)+";"
1188       
1189        NVAR rawCts = root:Packages:NIST:SAS:gRawCounts
1190        if(!rawCts)                     //go ahead and do the linear scaling
1191                linear_data = linear_data / kappa
1192        endif           
1193        data = linear_data
1194       
1195        // re-average the 2D data
1196        S_CircularAverageTo1D("SAS")
1197       
1198        // put the new result into the simulation folder
1199        Fake1DDataFolder(qval,aveint,sigave,sigmaQ,qbar,fSubs,"Simulation")     
1200                               
1201
1202        return(0)
1203end
1204
1205//phi is defined from +x axis, proceeding CCW around [0,2Pi]
1206ThreadSafe Function MC_FindPhi(vx,vy)
1207        variable vx,vy
1208       
1209        variable phi
1210       
1211        phi = atan(vy/vx)               //returns a value from -pi/2 to pi/2
1212       
1213        // special cases
1214        if(vx==0 && vy > 0)
1215                return(pi/2)
1216        endif
1217        if(vx==0 && vy < 0)
1218                return(3*pi/2)
1219        endif
1220        if(vx >= 0 && vy == 0)
1221                return(0)
1222        endif
1223        if(vx < 0 && vy == 0)
1224                return(pi)
1225        endif
1226       
1227       
1228        if(vx > 0 && vy > 0)
1229                return(phi)
1230        endif
1231        if(vx < 0 && vy > 0)
1232                return(phi + pi)
1233        endif
1234        if(vx < 0 && vy < 0)
1235                return(phi + pi)
1236        endif
1237        if( vx > 0 && vy < 0)
1238                return(phi + 2*pi)
1239        endif
1240       
1241        return(phi)
1242end
1243
1244
1245
1246
1247
1248Window Sim_1D_Panel() : Panel
1249        PauseUpdate; Silent 1           // building window...
1250        NewPanel /W=(92,556,713,818)/K=1 as "1D SANS Simulator"
1251        SetVariable cntVar,pos={26,68},size={160,15},title="Counting time(s)",format="%d"
1252        SetVariable cntVar,limits={1,36000,10},value= root:Packages:NIST:SAS:gCntTime
1253        SetVariable cntVar, proc=Sim_1D_CountTimeSetVarProc
1254        SetVariable MC_setvar0_1,pos={26,91},size={160,15},title="Thickness (cm)"
1255        SetVariable MC_setvar0_1,limits={0,inf,0.1},value= root:Packages:NIST:SAS:gThick
1256        SetVariable MC_setvar0_1, proc=Sim_1D_SamThickSetVarProc
1257
1258        SetVariable MC_setvar0_3,pos={26,114},size={160,15},title="Sample Transmission"
1259        SetVariable MC_setvar0_3,limits={0,1,0.01},value= root:Packages:NIST:SAS:gSamTrans
1260        SetVariable MC_setvar0_3, proc=Sim_1D_SamTransSetVarProc
1261
1262        PopupMenu MC_popup0,pos={13,13},size={165,20},proc=Sim_1D_ModelPopMenuProc,title="Model Function"
1263        PopupMenu MC_popup0,mode=1,value= #"MC_FunctionPopupList()"
1264        Button MC_button0,pos={17,181},size={130,20},proc=Sim_1D_DoItButtonProc,title="Do 1D Simulation"
1265        Button MC_button0,fColor=(3,52428,1)
1266        GroupBox group0,pos={15,42},size={280,130},title="Sample Setup"
1267        CheckBox check0_1,pos={216,179},size={60,14},title="Yes Offset",variable= root:Packages:NIST:SAS:gDoTraceOffset
1268        CheckBox check0_2,pos={216,199},size={60,14},title="Abs scale?",variable= root:Packages:NIST:SAS:g_1D_DoABS
1269        CheckBox check0_3,pos={216,219},size={60,14},title="Noise?",variable= root:Packages:NIST:SAS:g_1D_AddNoise
1270       
1271// a box for the results
1272        GroupBox group1,pos={314,23},size={277,163},title="Simulation Results"
1273        ValDisplay valdisp0,pos={326,48},size={220,13},title="Total detector counts"
1274        ValDisplay valdisp0,limits={0,0,0},barmisc={0,1000},value= root:Packages:NIST:SAS:g_1DTotCts
1275        ValDisplay valdisp0_1,pos={326,72},size={220,13},title="Estimated count rate (1/s)"
1276        ValDisplay valdisp0_1,limits={0,0,0},barmisc={0,1000},value=root:Packages:NIST:SAS:g_1DEstDetCR
1277        ValDisplay valdisp0_2,pos={326,96},size={220,13},title="Fraction of beam scattered"
1278        ValDisplay valdisp0_2,limits={0,0,0},barmisc={0,1000},value= root:Packages:NIST:SAS:g_1DFracScatt
1279        ValDisplay valdisp0_3,pos={326,121},size={220,13},title="Estimated transmission"
1280        ValDisplay valdisp0_3,limits={0,0,0},barmisc={0,1000},value=root:Packages:NIST:SAS:g_1DEstTrans
1281        // set the flags here -- do the simulation, but not 2D
1282       
1283        root:Packages:NIST:SAS:doSimulation     = 1     // == 1 if 1D simulated data, 0 if other from the checkbox
1284        root:Packages:NIST:SAS:gDoMonteCarlo     = 0  // == 1 if 2D MonteCarlo set by hidden flag
1285
1286       
1287EndMacro
1288
1289Function Sim_1D_CountTimeSetVarProc(sva) : SetVariableControl
1290        STRUCT WMSetVariableAction &sva
1291
1292        switch( sva.eventCode )
1293                case 1: // mouse up
1294                case 2: // Enter key
1295                case 3: // Live update
1296                        Variable dval = sva.dval
1297
1298                        ReCalculateInten(1)
1299                       
1300                        break
1301        endswitch
1302
1303        return 0
1304End
1305
1306Function Sim_1D_SamThickSetVarProc(sva) : SetVariableControl
1307        STRUCT WMSetVariableAction &sva
1308
1309        switch( sva.eventCode )
1310                case 1: // mouse up
1311                case 2: // Enter key
1312                case 3: // Live update
1313                        Variable dval = sva.dval
1314
1315                        ReCalculateInten(1)
1316                       
1317                        break
1318        endswitch
1319
1320        return 0
1321End
1322
1323Function Sim_1D_SamTransSetVarProc(sva) : SetVariableControl
1324        STRUCT WMSetVariableAction &sva
1325
1326        switch( sva.eventCode )
1327                case 1: // mouse up
1328                case 2: // Enter key
1329                case 3: // Live update
1330                        Variable dval = sva.dval
1331
1332                        ReCalculateInten(1)
1333                       
1334                        break
1335        endswitch
1336
1337        return 0
1338End
1339
1340
1341Function Sim_1D_ModelPopMenuProc(pa) : PopupMenuControl
1342        STRUCT WMPopupAction &pa
1343
1344        switch( pa.eventCode )
1345                case 2: // mouse up
1346                        Variable popNum = pa.popNum
1347                        String popStr = pa.popStr
1348                        SVAR gStr = root:Packages:NIST:SAS:gFuncStr
1349                        gStr = popStr
1350                       
1351                        break
1352        endswitch
1353
1354        return 0
1355End
1356
1357
1358Function Sim_1D_DoItButtonProc(ba) : ButtonControl
1359        STRUCT WMButtonAction &ba
1360
1361        switch( ba.eventCode )
1362                case 2: // mouse up
1363               
1364                        ReCalculateInten(1)
1365                       
1366                        break
1367        endswitch
1368
1369        return 0
1370End
Note: See TracBrowser for help on using the repository browser.