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

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

in SASCALC - changed some of the logic to allow 2D simulation

in MultiScatter_MC:

  • changed results wave to report coherent multiple scattering counters correctly (also changed XOP)
  • added save of 1D sim data

in U_CALC:

  • adjusted display for optimal view on windows
  • 7 angle ranges now (added defaults, propagated where all controls are modified)
  • added the angle axis correctly now as a free axis, with a hook linked to the bottom axis
  • removed XS and trans calculations, as integration of ran_dev isn't well

in USANS_EmptyWaves:

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