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

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

Added a calculation of the total fraction of multiple coherent scattering. This value is now displayed on the 1D sim panel and on the UCALC panel. If the total fraction of multiple coherent scattering (relative to the fraction coherently scattered) is greater than 0.10, then the field turns red as a warning.

Changed the numbering and suffix scheme for saving simulated VAX raw data files.

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