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

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

Updated the SingleModelTemplate? to be compatible with the new release

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