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

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

A variety of changes...

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