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

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

Added UCALC window to list of windows is package loader

Box Sum (from the marquee) asks for normalized (SAM) or RAW data

SASCALC uses real beamstop size rather than projected sizewhen doing simulation

Fake1DDataFolder() in MultiScatter_MonteCarlo ipf bug fixed where a bogus beamstop diameter (not 1"-4") would cause oddities in the resolution wave, and an (obviously) incorrect smearing of the simulation.

Fixed range 7 defaults in UCALC panel.

File size: 51.4 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// ---- use FakeUSANSDataFolder() if you want to fake a 1D USANS data set ----
985//
986Function        Fake1DDataFolder(qval,aveint,sigave,sigmaQ,qbar,fSubs,dataFolder)
987        WAVE qval,aveint,sigave,sigmaQ,qbar,fSubs
988        String dataFolder
989
990        String baseStr=dataFolder
991        if(DataFolderExists("root:"+baseStr))
992                SetDataFolder $("root:"+baseStr)
993        else
994                NewDataFolder/S $("root:"+baseStr)
995        endif
996
997        ////overwrite the existing data, if it exists
998        Duplicate/O qval, $(baseStr+"_q")
999        Duplicate/O aveint, $(baseStr+"_i")
1000        Duplicate/O sigave, $(baseStr+"_s")
1001
1002
1003        // make a resolution matrix for SANS data
1004        Variable np=numpnts(qval)
1005        Make/D/O/N=(np,4) $(baseStr+"_res")
1006        Wave res=$(baseStr+"_res")
1007       
1008        res[][0] = sigmaQ[p]            //sigQ
1009        res[][1] = qBar[p]              //qBar
1010        res[][2] = fSubS[p]             //fShad
1011        res[][3] = qval[p]              //Qvalues
1012       
1013        // keep a copy of everything in SAS too... the smearing wrapper function looks for
1014        // data in folders based on waves it is passed - an I lose control of that
1015        Duplicate/O res, $("root:Packages:NIST:SAS:"+baseStr+"_res")
1016        Duplicate/O qval,  $("root:Packages:NIST:SAS:"+baseStr+"_q")
1017        Duplicate/O aveint,  $("root:Packages:NIST:SAS:"+baseStr+"_i")
1018        Duplicate/O sigave,  $("root:Packages:NIST:SAS:"+baseStr+"_s")
1019       
1020        //clean up             
1021        SetDataFolder root:
1022       
1023End
1024
1025// writes out a VAX binary data file
1026// automatically generates a name
1027// will prompt for the sample label
1028//
1029Function SaveAsVAXButtonProc(ctrlName) : ButtonControl
1030        String ctrlName
1031
1032        WriteVAXData("SAS","",0)
1033End
1034
1035// calculates the fraction of the scattering that reaches the detector, given the random deviate function
1036// and qmin and qmax
1037//
1038//
1039// still some question of the corners and number of pixels per q-bin
1040Function FractionReachingDetector(ran_dev,Qmin,Qmax)
1041        wave ran_dev
1042        Variable Qmin,Qmax
1043       
1044        Variable r1,r2,frac
1045        r1=x2pnt(ran_dev,Qmin)
1046        r2=x2pnt(ran_dev,Qmax)
1047       
1048        // no normalization needed - the full q-range is defined as [0,1]
1049        frac = ran_dev[r2] - ran_dev[r1]
1050       
1051        return frac
1052End
1053
1054
1055/// called in SASCALC:ReCalculateInten()
1056Function        Simulate_2D_MC(funcStr,aveint,qval,sigave,sigmaq,qbar,fsubs)
1057        String funcStr
1058        WAVE aveint,qval,sigave,sigmaq,qbar,fsubs
1059
1060        NVAR doMonteCarlo = root:Packages:NIST:SAS:gDoMonteCarlo                // == 1 if 2D MonteCarlo set by hidden flag
1061        WAVE rw=root:Packages:NIST:SAS:realsRead
1062
1063        NVAR imon = root:Packages:NIST:SAS:gImon
1064        NVAR thick = root:Packages:NIST:SAS:gThick
1065        NVAR sig_incoh = root:Packages:NIST:SAS:gSig_incoh
1066        NVAR r2 = root:Packages:NIST:SAS:gR2
1067
1068        // do the simulation here, or not
1069        Variable r1,xCtr,yCtr,sdd,pixSize,wavelength
1070        String coefStr,abortStr,str
1071
1072        r1 = rw[24]/2/10                // sample diameter convert diam in [mm] to radius in cm
1073        xCtr = rw[16]
1074        yCtr = rw[17]
1075        sdd = rw[18]*100                //conver header of [m] to [cm]
1076        pixSize = rw[10]/10             // convert pix size in mm to cm
1077        wavelength = rw[26]
1078        coefStr = MC_getFunctionCoef(funcStr)
1079       
1080        if(!MC_CheckFunctionAndCoef(funcStr,coefStr))
1081                doMonteCarlo = 0                //we're getting out now, reset the flag so we don't continually end up here
1082                Abort "The coefficients and function type do not match. Please correct the selections in the popup menus."
1083        endif
1084       
1085        Variable sig_sas
1086//              FUNCREF SANSModelAAO_MCproto func=$("fSmeared"+funcStr)         //a wrapper for the structure version
1087        FUNCREF SANSModelAAO_MCproto func=$(funcStr)            //unsmeared
1088        WAVE results = root:Packages:NIST:SAS:results
1089        WAVE linear_data = root:Packages:NIST:SAS:linear_data
1090        WAVE data = root:Packages:NIST:SAS:data
1091
1092        results = 0
1093        linear_data = 0
1094       
1095        CalculateRandomDeviate(func,$coefStr,wavelength,"root:Packages:NIST:SAS:ran_dev",SIG_SAS)
1096        if(sig_sas > 100)
1097                sprintf abortStr,"sig_sas = %g. Please check that the model coefficients have a zero background, or the low q is well-behaved.",sig_sas
1098                Abort abortStr
1099        endif
1100       
1101        WAVE ran_dev=$"root:Packages:NIST:SAS:ran_dev"
1102       
1103        Make/O/D/N=5000 root:Packages:NIST:SAS:nt=0,root:Packages:NIST:SAS:j1=0,root:Packages:NIST:SAS:j2=0
1104        Make/O/D/N=100 root:Packages:NIST:SAS:nn=0
1105        Make/O/D/N=11 root:Packages:NIST:SAS:inputWave=0
1106       
1107        WAVE nt = root:Packages:NIST:SAS:nt
1108        WAVE j1 = root:Packages:NIST:SAS:j1
1109        WAVE j2 = root:Packages:NIST:SAS:j2
1110        WAVE nn = root:Packages:NIST:SAS:nn
1111        WAVE inputWave = root:Packages:NIST:SAS:inputWave
1112
1113        inputWave[0] = imon
1114        inputWave[1] = r1
1115        inputWave[2] = r2
1116        inputWave[3] = xCtr
1117        inputWave[4] = yCtr
1118        inputWave[5] = sdd
1119        inputWave[6] = pixSize
1120        inputWave[7] = thick
1121        inputWave[8] = wavelength
1122        inputWave[9] = sig_incoh
1123        inputWave[10] = sig_sas
1124
1125        linear_data = 0         //initialize
1126
1127        Variable t0,trans
1128       
1129        // get a time estimate, and give the user a chance to exit if they're unsure.
1130        t0 = stopMStimer(-2)
1131        inputWave[0] = 1000
1132        NVAR useXOP = root:Packages:NIST:SAS:gUse_MC_XOP                //if zero, will use non-threaded Igor code
1133       
1134        if(useXOP)
1135                //use a single thread, otherwise time is dominated by overhead
1136                Monte_SANS_NotThreaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
1137        else
1138                Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
1139        endif
1140       
1141        t0 = (stopMSTimer(-2) - t0)*1e-6
1142        t0 *= imon/1000/ThreadProcessorCount                    //projected time, in seconds (using threads for the calculation)
1143        inputWave[0] = imon             //reset
1144       
1145        if(t0>10)
1146                sprintf str,"The simulation will take approximately %d seconds.\r- Proceed?",t0
1147                DoAlert 1,str
1148                if(V_flag == 2)
1149                        doMonteCarlo = 0
1150                        reCalculateInten(1)             //come back in and do the smeared calculation
1151                        return(0)
1152                endif
1153        endif
1154       
1155        linear_data = 0         //initialize
1156// threading crashes!! - there must be some operation in the XOP that is not threadSafe. What, I don't know...
1157// I think it's the ran() calls, being "non-reentrant". So the XOP now defines two separate functions, that each
1158// use a different rng. This works. 1.75x speedup.     
1159        t0 = stopMStimer(-2)
1160
1161        if(useXOP)
1162                Monte_SANS_Threaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
1163        else
1164                Monte_SANS_NotThreaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
1165        endif
1166       
1167        t0 = (stopMSTimer(-2) - t0)*1e-6
1168        Printf  "MC sim time = %g seconds\r",t0
1169       
1170        trans = results[8]                      //(n1-n2)/n1
1171        if(trans == 0)
1172                trans = 1
1173        endif
1174
1175        Print "counts on detector, including transmitted = ",sum(linear_data,-inf,inf)
1176       
1177//              linear_data[xCtr][yCtr] = 0                     //snip out the transmitted spike
1178//              Print "counts on detector not transmitted = ",sum(linear_data,-inf,inf)
1179
1180        // or simulate a beamstop
1181        NVAR MC_BS_in = root:Packages:NIST:SAS:gBeamStopIn              //if zero, beam stop is "out", as in a transmission measurement
1182       
1183        Variable rad=beamstopDiam()/2           //beamstop radius in cm
1184        if(MC_BS_in)
1185                rad /= 0.5                              //convert cm to pixels
1186                rad += 0.                                       // (no - it cuts off the low Q artificially) add an extra pixel to each side to account for edge
1187                Duplicate/O linear_data,root:Packages:NIST:SAS:tmp_mask//,root:Packages:NIST:SAS:MC_linear_data
1188                WAVE tmp_mask = root:Packages:NIST:SAS:tmp_mask
1189                tmp_mask = (sqrt((p-xCtr)^2+(q-yCtr)^2) < rad) ? 0 : 1          //behind beamstop = 0, away = 1
1190               
1191                linear_data *= tmp_mask
1192        endif
1193       
1194        results[9] = sum(linear_data,-inf,inf)
1195        //              Print "counts on detector not behind beamstop = ",results[9]
1196       
1197        // convert to absolute scale
1198        Variable kappa          //,beaminten = beamIntensity()
1199//              kappa = beamInten*pi*r1*r1*thick*(pixSize/sdd)^2*trans*(iMon/beaminten)
1200        kappa = thick*(pixSize/sdd)^2*trans*iMon
1201       
1202        //use kappa to get back to counts => linear_data = round(linear_data*kappa)
1203        Note/K linear_data ,"KAPPA="+num2str(kappa)+";"
1204       
1205        NVAR rawCts = root:Packages:NIST:SAS:gRawCounts
1206        if(!rawCts)                     //go ahead and do the linear scaling
1207                linear_data = linear_data / kappa
1208        endif           
1209        data = linear_data
1210       
1211        // re-average the 2D data
1212        S_CircularAverageTo1D("SAS")
1213       
1214        // put the new result into the simulation folder
1215        Fake1DDataFolder(qval,aveint,sigave,sigmaQ,qbar,fSubs,"Simulation")     
1216                               
1217
1218        return(0)
1219end
1220
1221//phi is defined from +x axis, proceeding CCW around [0,2Pi]
1222ThreadSafe Function MC_FindPhi(vx,vy)
1223        variable vx,vy
1224       
1225        variable phi
1226       
1227        phi = atan(vy/vx)               //returns a value from -pi/2 to pi/2
1228       
1229        // special cases
1230        if(vx==0 && vy > 0)
1231                return(pi/2)
1232        endif
1233        if(vx==0 && vy < 0)
1234                return(3*pi/2)
1235        endif
1236        if(vx >= 0 && vy == 0)
1237                return(0)
1238        endif
1239        if(vx < 0 && vy == 0)
1240                return(pi)
1241        endif
1242       
1243       
1244        if(vx > 0 && vy > 0)
1245                return(phi)
1246        endif
1247        if(vx < 0 && vy > 0)
1248                return(phi + pi)
1249        endif
1250        if(vx < 0 && vy < 0)
1251                return(phi + pi)
1252        endif
1253        if( vx > 0 && vy < 0)
1254                return(phi + 2*pi)
1255        endif
1256       
1257        return(phi)
1258end
1259
1260
1261
1262
1263
1264Window Sim_1D_Panel() : Panel
1265        PauseUpdate; Silent 1           // building window...
1266        NewPanel /W=(92,556,713,818)/K=1 as "1D SANS Simulator"
1267        SetVariable cntVar,pos={26,68},size={160,15},title="Counting time(s)",format="%d"
1268        SetVariable cntVar,limits={1,36000,10},value= root:Packages:NIST:SAS:gCntTime
1269        SetVariable cntVar, proc=Sim_1D_CountTimeSetVarProc
1270        SetVariable MC_setvar0_1,pos={26,91},size={160,15},title="Thickness (cm)"
1271        SetVariable MC_setvar0_1,limits={0,inf,0.1},value= root:Packages:NIST:SAS:gThick
1272        SetVariable MC_setvar0_1, proc=Sim_1D_SamThickSetVarProc
1273
1274        SetVariable MC_setvar0_3,pos={26,114},size={160,15},title="Sample Transmission"
1275        SetVariable MC_setvar0_3,limits={0,1,0.01},value= root:Packages:NIST:SAS:gSamTrans
1276        SetVariable MC_setvar0_3, proc=Sim_1D_SamTransSetVarProc
1277
1278        PopupMenu MC_popup0,pos={13,13},size={165,20},proc=Sim_1D_ModelPopMenuProc,title="Model Function"
1279        PopupMenu MC_popup0,mode=1,value= #"MC_FunctionPopupList()"
1280        Button MC_button0,pos={17,181},size={130,20},proc=Sim_1D_DoItButtonProc,title="Do 1D Simulation"
1281        Button MC_button0,fColor=(3,52428,1)
1282        Button MC_button1,pos={17,211},size={150,20},proc=Save_1DSimData,title="Save Simulated Data"
1283        GroupBox group0,pos={15,42},size={280,130},title="Sample Setup"
1284        CheckBox check0_1,pos={216,179},size={60,14},title="Yes Offset",variable= root:Packages:NIST:SAS:gDoTraceOffset
1285        CheckBox check0_2,pos={216,199},size={60,14},title="Abs scale?",variable= root:Packages:NIST:SAS:g_1D_DoABS
1286        CheckBox check0_3,pos={216,219},size={60,14},title="Noise?",variable= root:Packages:NIST:SAS:g_1D_AddNoise
1287       
1288// a box for the results
1289        GroupBox group1,pos={314,23},size={277,163},title="Simulation Results"
1290        ValDisplay valdisp0,pos={326,48},size={220,13},title="Total detector counts"
1291        ValDisplay valdisp0,limits={0,0,0},barmisc={0,1000},value= root:Packages:NIST:SAS:g_1DTotCts
1292        ValDisplay valdisp0_1,pos={326,72},size={220,13},title="Estimated count rate (1/s)"
1293        ValDisplay valdisp0_1,limits={0,0,0},barmisc={0,1000},value=root:Packages:NIST:SAS:g_1DEstDetCR
1294        ValDisplay valdisp0_2,pos={326,96},size={220,13},title="Fraction of beam scattered"
1295        ValDisplay valdisp0_2,limits={0,0,0},barmisc={0,1000},value= root:Packages:NIST:SAS:g_1DFracScatt
1296        ValDisplay valdisp0_3,pos={326,121},size={220,13},title="Estimated transmission"
1297        ValDisplay valdisp0_3,limits={0,0,0},barmisc={0,1000},value=root:Packages:NIST:SAS:g_1DEstTrans
1298        // set the flags here -- do the simulation, but not 2D
1299       
1300        root:Packages:NIST:SAS:doSimulation     = 1     // == 1 if 1D simulated data, 0 if other from the checkbox
1301        root:Packages:NIST:SAS:gDoMonteCarlo     = 0  // == 1 if 2D MonteCarlo set by hidden flag
1302
1303       
1304EndMacro
1305
1306Function Sim_1D_CountTimeSetVarProc(sva) : SetVariableControl
1307        STRUCT WMSetVariableAction &sva
1308
1309        switch( sva.eventCode )
1310                case 1: // mouse up
1311                case 2: // Enter key
1312                case 3: // Live update
1313                        Variable dval = sva.dval
1314
1315                        ReCalculateInten(1)
1316                       
1317                        break
1318        endswitch
1319
1320        return 0
1321End
1322
1323Function Sim_1D_SamThickSetVarProc(sva) : SetVariableControl
1324        STRUCT WMSetVariableAction &sva
1325
1326        switch( sva.eventCode )
1327                case 1: // mouse up
1328                case 2: // Enter key
1329                case 3: // Live update
1330                        Variable dval = sva.dval
1331
1332                        ReCalculateInten(1)
1333                       
1334                        break
1335        endswitch
1336
1337        return 0
1338End
1339
1340Function Sim_1D_SamTransSetVarProc(sva) : SetVariableControl
1341        STRUCT WMSetVariableAction &sva
1342
1343        switch( sva.eventCode )
1344                case 1: // mouse up
1345                case 2: // Enter key
1346                case 3: // Live update
1347                        Variable dval = sva.dval
1348
1349                        ReCalculateInten(1)
1350                       
1351                        break
1352        endswitch
1353
1354        return 0
1355End
1356
1357
1358Function Sim_1D_ModelPopMenuProc(pa) : PopupMenuControl
1359        STRUCT WMPopupAction &pa
1360
1361        switch( pa.eventCode )
1362                case 2: // mouse up
1363                        Variable popNum = pa.popNum
1364                        String popStr = pa.popStr
1365                        SVAR gStr = root:Packages:NIST:SAS:gFuncStr
1366                        gStr = popStr
1367                       
1368                        break
1369        endswitch
1370
1371        return 0
1372End
1373
1374
1375Function Sim_1D_DoItButtonProc(ba) : ButtonControl
1376        STRUCT WMButtonAction &ba
1377
1378        switch( ba.eventCode )
1379                case 2: // mouse up
1380               
1381                        ReCalculateInten(1)
1382                       
1383                        break
1384        endswitch
1385
1386        return 0
1387End
1388
1389
1390//
1391//
1392//
1393Function Save_1DSimData(ctrlName) : ButtonControl
1394        String ctrlName
1395
1396        String type="SAS",fullpath=""
1397        Variable dialog=1               //=1 will present dialog for name
1398       
1399        String destStr=""
1400        destStr = "root:Packages:NIST:"+type
1401       
1402        Variable refNum
1403        String formatStr = "%15.4g %15.4g %15.4g %15.4g %15.4g %15.4g\r\n"
1404        String fname,ave="C",hdrStr1="",hdrStr2=""
1405        Variable step=1
1406       
1407        If(1)
1408                //setup a "fake protocol" wave, sice I have no idea of the current state of the data
1409                Make/O/T/N=8 root:myGlobals:Protocols:SIMProtocol
1410                Wave/T SIMProtocol = $"root:myGlobals:Protocols:SIMProtocol"
1411                String junk="****SIMULATED DATA****"
1412                //stick in the fake protocol...
1413                NVAR ctTime = root:Packages:NIST:SAS:gCntTime
1414                NVAR totalCts = root:Packages:NIST:SAS:g_1DTotCts                       //summed counts (simulated)
1415                NVAR detCR = root:Packages:NIST:SAS:g_1DEstDetCR                // estimated detector count rate
1416                NVAR fractScat = root:Packages:NIST:SAS:g_1DFracScatt
1417       
1418                SIMProtocol[0] = junk
1419                SIMProtocol[1] = "\tCounting time (s) = "+num2str(ctTime)
1420                SIMProtocol[2] = "\tTotal detector counts = "+num2str(totalCts)
1421                SIMProtocol[3] = "\tDetector countrate (1/s) = "+num2str(detCR)
1422                SIMProtocol[4] = "\tFraction of beam scattered coherently = "+num2str(fractScat)
1423                SIMProtocol[5] = junk
1424                SIMProtocol[6] = ""
1425                SIMProtocol[7] = ""
1426                //set the global
1427                String/G root:myGlobals:Protocols:gProtoStr = "SIMProtocol"
1428        Endif
1429       
1430       
1431        //*****these waves MUST EXIST, or IGOR Pro will crash, with a type 2 error****
1432        WAVE intw=$(destStr + ":integersRead")
1433        WAVE rw=$(destStr + ":realsRead")
1434        WAVE/T textw=$(destStr + ":textRead")
1435        WAVE qvals =$(destStr + ":qval")
1436        WAVE inten=$(destStr + ":aveint")
1437        WAVE sig=$(destStr + ":sigave")
1438        WAVE qbar = $(destStr + ":QBar")
1439        WAVE sigmaq = $(destStr + ":SigmaQ")
1440        WAVE fsubs = $(destStr + ":fSubS")
1441
1442        SVAR gProtoStr = root:myGlobals:Protocols:gProtoStr
1443        Wave/T proto=$("root:myGlobals:Protocols:"+gProtoStr)
1444       
1445        //check each wave
1446        If(!(WaveExists(intw)))
1447                Abort "intw DNExist Save_1DSimData()"
1448        Endif
1449        If(!(WaveExists(rw)))
1450                Abort "rw DNExist Save_1DSimData()"
1451        Endif
1452        If(!(WaveExists(textw)))
1453                Abort "textw DNExist Save_1DSimData()"
1454        Endif
1455        If(!(WaveExists(qvals)))
1456                Abort "qvals DNExist Save_1DSimData()"
1457        Endif
1458        If(!(WaveExists(inten)))
1459                Abort "inten DNExist Save_1DSimData()"
1460        Endif
1461        If(!(WaveExists(sig)))
1462                Abort "sig DNExist Save_1DSimData()"
1463        Endif
1464        If(!(WaveExists(qbar)))
1465                Abort "qbar DNExist Save_1DSimData()"
1466        Endif
1467        If(!(WaveExists(sigmaq)))
1468                Abort "sigmaq DNExist Save_1DSimData()"
1469        Endif
1470        If(!(WaveExists(fsubs)))
1471                Abort "fsubs DNExist Save_1DSimData()"
1472        Endif
1473        If(!(WaveExists(proto)))
1474                Abort "current protocol wave DNExist Save_1DSimData()"
1475        Endif
1476
1477        //strings can be too long to print-- must trim to 255 chars
1478        Variable ii,num=8
1479        Make/O/T/N=(num) tempShortProto
1480        for(ii=0;ii<num;ii+=1)
1481                tempShortProto[ii] = (proto[ii])[0,240]
1482        endfor
1483       
1484        if(dialog)
1485                PathInfo/S catPathName
1486                fullPath = DoSaveFileDialog("Save data as")
1487                If(cmpstr(fullPath,"")==0)
1488                        //user cancel, don't write out a file
1489                        Close/A
1490                        Abort "no data file was written"
1491                Endif
1492                //Print "dialog fullpath = ",fullpath
1493        Endif
1494       
1495        NVAR monCt = root:Packages:NIST:SAS:gImon
1496        NVAR thick = root:Packages:NIST:SAS:gThick
1497        NVAR trans = root:Packages:NIST:SAS:gSamTrans                   //for 1D, default value
1498       
1499
1500       
1501        hdrStr1 = num2str(monCt)+"  "+num2str(rw[26])+"       "+num2str(rw[19])+"     "+num2str(rw[18])
1502        hdrStr1 += "     "+num2str(trans)+"     "+num2str(thick) + ave +"   "+num2str(step) + "\r\n"
1503
1504        hdrStr2 = num2str(rw[16])+"  "+num2str(rw[17])+"  "+num2str(rw[23])+"    "+num2str(rw[24])+"    "
1505        hdrStr2 += num2str(rw[25])+"    "+num2str(rw[27])+"    "+num2str(rw[21])+"    "+"ORNL  " + "\r\n"
1506       
1507        //actually open the file here
1508        Open refNum as fullpath
1509       
1510        //write out the standard header information
1511        fprintf refnum,"FILE: %s\t\t CREATED: %s\r\n","SIMULATED DATA",(date() +"  "+ time())
1512        fprintf refnum,"LABEL: %s\r\n","SIMULATED DATA"
1513        fprintf refnum,"MON CNT   LAMBDA   DET ANG   DET DIST   TRANS   THICK   AVE   STEP\r\n"
1514        fprintf refnum,hdrStr1
1515        fprintf refnum,"BCENT(X,Y)   A1(mm)   A2(mm)   A1A2DIST(m)   DL/L   BSTOP(mm)   DET_TYP \r\n"
1516        fprintf refnum,hdrStr2
1517//      fprintf refnum,headerFormat,rw[0],rw[26],rw[19],rw[18],rw[4],rw[5],ave,step
1518
1519        //insert protocol information here
1520        //-1 list of sample files
1521        //0 - bkg
1522        //1 - emp
1523        //2 - div
1524        //3 - mask
1525        //4 - abs params c2-c5
1526        //5 - average params
1527        fprintf refnum, "SAM: %s\r\n",tempShortProto[0]
1528        fprintf refnum, "BGD: %s\r\n",tempShortProto[1]
1529        fprintf refnum, "EMP: %s\r\n",tempShortProto[2]
1530        fprintf refnum, "DIV: %s\r\n",tempShortProto[3]
1531        fprintf refnum, "MASK: %s\r\n",tempShortProto[4]
1532        fprintf refnum, "ABS: %s\r\n",tempShortProto[5]
1533        fprintf refnum, "Average Choices: %s\r\n",tempShortProto[6]
1534       
1535        //write out the data columns
1536        fprintf refnum,"The 6 columns are | Q (1/A) | I(Q) (1/cm) | std. dev. I(Q) (1/cm) | sigmaQ | meanQ | ShadowFactor|\r\n"
1537        wfprintf refnum, formatStr, qvals,inten,sig,sigmaq,qbar,fsubs
1538       
1539        Close refnum
1540       
1541        SetDataFolder root:             //(redundant)
1542       
1543        //write confirmation of write operation to history area
1544        Print "Averaged File written: ", GetFileNameFromPathNoSemi(fullPath)
1545        KillWaves/Z tempShortProto
1546
1547        //clear the stuff that was created for case of saving files
1548        If(1)
1549                Killwaves/Z root:myGlobals:Protocols:SIMProtocol
1550                String/G root:myGlobals:Protocols:gProtoStr = ""
1551        Endif
1552       
1553       
1554        return(0)
1555       
1556End
Note: See TracBrowser for help on using the repository browser.