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

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

Added conditional compile instructions to some reduction procedure files that have only minor facility-specific changes. These changes are not significant enough to merit a separate facility file that must be maintained with essentially duplicate functions.

#define SYMBOL was attempted, but did not work since the symbols weren't actually defined until after the compile... and I couldn't figure out how to define - then compile. in additon, the table is static until Igor is quit - so multiple symbols could be defined, and compiling would fail.

So... the method now that appears to work is to put a dummy function for each facility in its facility specific "Includes" file. Then the conditional compilation checks for exists("function").

File size: 58.6 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        //if on detector, return xPix and yPix values, otherwise -1
723        if(yPixel > 127 || yPixel < 0)
724                yPixel = -1
725        endif
726        if(xPixel > 127 || xPixel < 0)
727                xPixel = -1
728        endif
729       
730        return(0)
731End
732
733Function MC_CheckFunctionAndCoef(funcStr,coefStr)
734        String funcStr,coefStr
735       
736        SVAR/Z listStr=root:Packages:NIST:coefKWStr
737        if(SVAR_Exists(listStr) == 1)
738                String properCoefStr = StringByKey(funcStr, listStr  ,"=",";",0)
739                if(cmpstr("",properCoefStr)==0)
740                        return(0)               //false, no match found, so properCoefStr is returned null
741                endif
742                if(cmpstr(coefStr,properCoefStr)==0)
743                        return(1)               //true, the coef is the correct match
744                endif
745        endif
746        return(0)                       //false, wrong coef
747End
748
749Function/S MC_getFunctionCoef(funcStr)
750        String funcStr
751
752        SVAR/Z listStr=root:Packages:NIST:coefKWStr
753        String coefStr=""
754        if(SVAR_Exists(listStr) == 1)
755                coefStr = StringByKey(funcStr, listStr  ,"=",";",0)
756        endif
757        return(coefStr)
758End
759
760Function SANSModelAAO_MCproto(w,yw,xw)
761        Wave w,yw,xw
762       
763        Print "in SANSModelAAO_MCproto function"
764        return(1)
765end
766
767Function/S MC_FunctionPopupList()
768        String list,tmp
769        list = User_FunctionPopupList()
770       
771        //simplify the display, forcing smeared calculations behind the scenes
772        tmp = FunctionList("Smear*",";","NPARAMS:1")            //smeared dependency calculations
773        list = RemoveFromList(tmp, list,";")
774
775
776        if(strlen(list)==0)
777                list = "No functions plotted"
778        endif
779       
780        list = SortList(list)
781       
782        return(list)
783End             
784
785
786//Function Scat_Angle(Ran,R_DAB,wavelength)
787//      Variable Ran,r_dab,wavelength
788//
789//      Variable qq,arg,theta
790//      qq = 1. / ( R_DAB*sqrt(1.0/Ran - 1.0) )
791//      arg = qq*wavelength/(4*pi)
792//      If (arg < 1.0)
793//              theta = 2.*asin(arg)
794//      else
795//              theta = pi
796//      endif
797//      Return (theta)
798//End
799
800//calculates new direction (xyz) from an old direction
801//theta and phi don't change
802ThreadSafe Function NewDirection(vx,vy,vz,theta,phi)
803        Variable &vx,&vy,&vz
804        Variable theta,phi
805       
806        Variable err=0,vx0,vy0,vz0
807        Variable nx,ny,mag_xy,tx,ty,tz
808       
809        //store old direction vector
810        vx0 = vx
811        vy0 = vy
812        vz0 = vz
813       
814        mag_xy = sqrt(vx0*vx0 + vy0*vy0)
815        if(mag_xy < 1e-12)
816                //old vector lies along beam direction
817                nx = 0
818                ny = 1
819                tx = 1
820                ty = 0
821                tz = 0
822        else
823                Nx = -Vy0 / Mag_XY
824                Ny = Vx0 / Mag_XY
825                Tx = -Vz0*Vx0 / Mag_XY
826                Ty = -Vz0*Vy0 / Mag_XY
827                Tz = Mag_XY
828        endif
829       
830        //new scattered direction vector
831        Vx = cos(phi)*sin(theta)*Tx + sin(phi)*sin(theta)*Nx + cos(theta)*Vx0
832        Vy = cos(phi)*sin(theta)*Ty + sin(phi)*sin(theta)*Ny + cos(theta)*Vy0
833        Vz = cos(phi)*sin(theta)*Tz + cos(theta)*Vz0
834       
835        Return(err)
836End
837
838ThreadSafe Function path_len(aval,sig_tot)
839        Variable aval,sig_tot
840       
841        Variable retval
842       
843        retval = -1*ln(1-aval)/sig_tot
844       
845        return(retval)
846End
847
848// globals are initialized in SASCALC.ipf
849// coordinates if I make this part of the panel - but this breaks other things...
850//
851//Proc MC_SASCALC()
852////    PauseUpdate; Silent 1           // building window...
853//
854////    NewPanel /W=(92,556,390,1028)/K=1 as "SANS Simulator"
855//      SetVariable MC_setvar0,pos={491,73},size={144,15},bodyWidth=80,title="# of neutrons"
856//      SetVariable MC_setvar0,format="%5.4g"
857//      SetVariable MC_setvar0,limits={-inf,inf,100},value= root:Packages:NIST:SAS:gImon
858//      SetVariable MC_setvar0_1,pos={491,119},size={131,15},bodyWidth=60,title="Thickness (cm)"
859//      SetVariable MC_setvar0_1,limits={-inf,inf,0.1},value= root:Packages:NIST:SAS:gThick
860//      SetVariable MC_setvar0_2,pos={491,96},size={149,15},bodyWidth=60,title="Incoherent XS (cm)"
861//      SetVariable MC_setvar0_2,limits={-inf,inf,0.1},value= root:Packages:NIST:SAS:gSig_incoh
862//      SetVariable MC_setvar0_3,pos={491,142},size={150,15},bodyWidth=60,title="Sample Radius (cm)"
863//      SetVariable MC_setvar0_3,limits={-inf,inf,0.1},value= root:Packages:NIST:SAS:gR2
864//      PopupMenu MC_popup0,pos={476,13},size={165,20},proc=MC_ModelPopMenuProc,title="Model Function"
865//      PopupMenu MC_popup0,mode=1,value= #"MC_FunctionPopupList()"
866//      Button MC_button0,pos={480,181},size={130,20},proc=MC_DoItButtonProc,title="Do MC Simulation"
867//      Button MC_button1,pos={644,181},size={80,20},proc=MC_Display2DButtonProc,title="Show 2D"
868//      SetVariable setvar0_3,pos={568,484},size={50,20},disable=1
869//      GroupBox group0,pos={478,42},size={267,130},title="Monte Carlo"
870//      SetVariable cntVar,pos={653,73},size={80,15},proc=CountTimeSetVarProc,title="time(s)"
871//      SetVariable cntVar,format="%d"
872//      SetVariable cntVar,limits={1,10,1},value= root:Packages:NIST:SAS:gCntTime
873//     
874//      String fldrSav0= GetDataFolder(1)
875//      SetDataFolder root:Packages:NIST:SAS:
876//      Edit/W=(476,217,746,450)/HOST=#  results_desc,results
877//      ModifyTable format(Point)=1,width(Point)=0,width(results_desc)=150
878//      SetDataFolder fldrSav0
879//      RenameWindow #,T_results
880//      SetActiveSubwindow ##
881//EndMacro
882
883// as a stand-alone panel, extra control bar  (right) and subwindow implementations don't work right
884// for various reasons...
885Window MC_SASCALC() : Panel
886
887        // when opening the panel, set the raw counts check to 1
888        root:Packages:NIST:SAS:gRawCounts = 1
889       
890        PauseUpdate; Silent 1           // building window...
891        NewPanel /W=(92,556,713,818)/K=1 as "SANS Simulator"
892        SetVariable MC_setvar0,pos={28,73},size={144,15},bodyWidth=80,title="# of neutrons"
893        SetVariable MC_setvar0,format="%5.4g"
894        SetVariable MC_setvar0,limits={0,inf,100},value= root:Packages:NIST:SAS:gImon
895        SetVariable MC_setvar0_1,pos={28,119},size={131,15},bodyWidth=60,title="Thickness (cm)"
896        SetVariable MC_setvar0_1,limits={0,inf,0.1},value= root:Packages:NIST:SAS:gThick
897        SetVariable MC_setvar0_2,pos={28,96},size={149,15},bodyWidth=60,title="Incoherent XS (cm)"
898        SetVariable MC_setvar0_2,limits={0,inf,0.1},value= root:Packages:NIST:SAS:gSig_incoh
899        SetVariable MC_setvar0_3,pos={28,142},size={150,15},bodyWidth=60,title="Sample Radius (cm)"
900        SetVariable MC_setvar0_3,limits={-inf,inf,0.1},value= root:Packages:NIST:SAS:gR2
901        PopupMenu MC_popup0,pos={13,13},size={165,20},proc=MC_ModelPopMenuProc,title="Model Function"
902        PopupMenu MC_popup0,mode=1,value= #"MC_FunctionPopupList()"
903        Button MC_button0,pos={17,181},size={130,20},proc=MC_DoItButtonProc,title="Do MC Simulation"
904        Button MC_button0,fColor=(3,52428,1)
905        Button MC_button1,pos={17,208},size={80,20},proc=MC_Display2DButtonProc,title="Show 2D"
906        SetVariable setvar0_3,pos={105,484},size={50,20},disable=1
907        GroupBox group0,pos={15,42},size={267,130},title="Monte Carlo"
908        SetVariable cntVar,pos={185,73},size={90,15},proc=CountTimeSetVarProc,title="time(s)"
909        SetVariable cntVar,format="%d"
910        SetVariable cntVar,limits={1,3600,1},value= root:Packages:NIST:SAS:gCntTime
911        Button MC_button2,pos={17,234},size={100,20},proc=SaveAsVAXButtonProc,title="Save 2D VAX"
912        CheckBox check0,pos={216,180},size={68,14},title="Raw counts",variable = root:Packages:NIST:SAS:gRawCounts
913        CheckBox check0_1,pos={216,199},size={60,14},title="Yes Offset",variable= root:Packages:NIST:SAS:gDoTraceOffset
914        CheckBox check0_2,pos={216,199+19},size={60,14},title="Beam Stop in",variable= root:Packages:NIST:SAS:gBeamStopIn
915        CheckBox check0_3,pos={216,199+2*19},size={60,14},title="use XOP",variable= root:Packages:NIST:SAS:gUse_MC_XOP
916       
917        String fldrSav0= GetDataFolder(1)
918        SetDataFolder root:Packages:NIST:SAS:
919        Edit/W=(344,23,606,248)/HOST=#  results_desc,results
920        ModifyTable format(Point)=1,width(Point)=0,width(results_desc)=150
921        SetDataFolder fldrSav0
922        RenameWindow #,T_results
923        SetActiveSubwindow ##
924EndMacro
925
926Function CountTimeSetVarProc(sva) : SetVariableControl
927        STRUCT WMSetVariableAction &sva
928
929        switch( sva.eventCode )
930                case 1: // mouse up
931                case 2: // Enter key
932                case 3: // Live update
933                        Variable dval = sva.dval
934
935                        // get the neutron flux, multiply, and reset the global for # neutrons
936                        NVAR imon=root:Packages:NIST:SAS:gImon
937                        imon = dval*beamIntensity()
938                       
939                        break
940        endswitch
941
942        return 0
943End
944
945
946Function MC_ModelPopMenuProc(pa) : PopupMenuControl
947        STRUCT WMPopupAction &pa
948
949        switch( pa.eventCode )
950                case 2: // mouse up
951                        Variable popNum = pa.popNum
952                        String popStr = pa.popStr
953                        SVAR gStr = root:Packages:NIST:SAS:gFuncStr
954                        gStr = popStr
955                       
956                        break
957        endswitch
958
959        return 0
960End
961
962Function MC_DoItButtonProc(ba) : ButtonControl
963        STRUCT WMButtonAction &ba
964
965        switch( ba.eventCode )
966                case 2: // mouse up
967                        // click code here
968                        NVAR doMC = root:Packages:NIST:SAS:gDoMonteCarlo
969                        doMC = 1
970                        ReCalculateInten(1)
971                        doMC = 0                //so the next time won't be MC
972                        break
973        endswitch
974
975        return 0
976End
977
978
979Function MC_Display2DButtonProc(ba) : ButtonControl
980        STRUCT WMButtonAction &ba
981
982        switch( ba.eventCode )
983                case 2: // mouse up
984                        // click code here
985                        Execute "ChangeDisplay(\"SAS\")"
986                        break
987        endswitch
988
989        return 0
990End
991
992// after a 2d data image is averaged in the usual way, take the waves and generate a "fake" folder of the 1d
993// data, to appear as if it was loaded from a real data file.
994//
995// ---- use FakeUSANSDataFolder() if you want to fake a 1D USANS data set ----
996//
997Function        Fake1DDataFolder(qval,aveint,sigave,sigmaQ,qbar,fSubs,dataFolder)
998        WAVE qval,aveint,sigave,sigmaQ,qbar,fSubs
999        String dataFolder
1000
1001        String baseStr=dataFolder
1002        if(DataFolderExists("root:"+baseStr))
1003                SetDataFolder $("root:"+baseStr)
1004        else
1005                NewDataFolder/S $("root:"+baseStr)
1006        endif
1007
1008        ////overwrite the existing data, if it exists
1009        Duplicate/O qval, $(baseStr+"_q")
1010        Duplicate/O aveint, $(baseStr+"_i")
1011        Duplicate/O sigave, $(baseStr+"_s")
1012
1013
1014        // make a resolution matrix for SANS data
1015        Variable np=numpnts(qval)
1016        Make/D/O/N=(np,4) $(baseStr+"_res")
1017        Wave res=$(baseStr+"_res")
1018       
1019        res[][0] = sigmaQ[p]            //sigQ
1020        res[][1] = qBar[p]              //qBar
1021        res[][2] = fSubS[p]             //fShad
1022        res[][3] = qval[p]              //Qvalues
1023       
1024        // keep a copy of everything in SAS too... the smearing wrapper function looks for
1025        // data in folders based on waves it is passed - an I lose control of that
1026        Duplicate/O res, $("root:Packages:NIST:SAS:"+baseStr+"_res")
1027        Duplicate/O qval,  $("root:Packages:NIST:SAS:"+baseStr+"_q")
1028        Duplicate/O aveint,  $("root:Packages:NIST:SAS:"+baseStr+"_i")
1029        Duplicate/O sigave,  $("root:Packages:NIST:SAS:"+baseStr+"_s")
1030       
1031        //clean up             
1032        SetDataFolder root:
1033       
1034End
1035
1036// writes out a VAX binary data file
1037// automatically generates a name
1038// will prompt for the sample label
1039//
1040// currently hard-wired for SAS data folder
1041//
1042Function SaveAsVAXButtonProc(ctrlName) : ButtonControl
1043        String ctrlName
1044
1045        String fullpath="",destStr=""
1046        Variable refnum
1047       
1048        fullpath = Write_RawData_File("SAS","",0)
1049       
1050        // write out the results into a text file
1051        destStr = "root:Packages:NIST:SAS:"
1052        SetDataFolder $destStr
1053
1054        WAVE results=results
1055        WAVE/T results_desc=results_desc
1056       
1057        //check each wave
1058        If(!(WaveExists(results)))
1059                Abort "results DNExist WriteVAXData()"
1060        Endif
1061        If(!(WaveExists(results_desc)))
1062                Abort "results_desc DNExist WriteVAXData()"
1063        Endif
1064       
1065        Open refNum as fullpath+".txt"
1066                wfprintf refNum, "%30s\t\t%g\r",results_desc,results
1067                FStatus refNum
1068                FSetPos refNum,V_logEOF
1069        Close refNum
1070       
1071        ///////////////////////////////
1072       
1073        // could also automatically do the average here, but probably not worth the the effort...
1074       
1075        SetDataFolder root:
1076       
1077        return(0)
1078End
1079
1080// calculates the fraction of the scattering that reaches the detector, given the random deviate function
1081// and qmin and qmax
1082//
1083//
1084// still some question of the corners and number of pixels per q-bin
1085Function FractionReachingDetector(ran_dev,Qmin,Qmax)
1086        wave ran_dev
1087        Variable Qmin,Qmax
1088       
1089        Variable r1,r2,frac
1090        r1=x2pnt(ran_dev,Qmin)
1091        r2=x2pnt(ran_dev,Qmax)
1092       
1093        // no normalization needed - the full q-range is defined as [0,1]
1094        frac = ran_dev[r2] - ran_dev[r1]
1095       
1096        return frac
1097End
1098
1099
1100/// called in SASCALC:ReCalculateInten()
1101Function        Simulate_2D_MC(funcStr,aveint,qval,sigave,sigmaq,qbar,fsubs)
1102        String funcStr
1103        WAVE aveint,qval,sigave,sigmaq,qbar,fsubs
1104
1105        NVAR doMonteCarlo = root:Packages:NIST:SAS:gDoMonteCarlo                // == 1 if 2D MonteCarlo set by hidden flag
1106        WAVE rw=root:Packages:NIST:SAS:realsRead
1107       
1108// Try to nicely exit from a threading error, if possible
1109        Variable err=0
1110        if(!exists("root:myGlobals:gThreadGroupID"))
1111                Variable/G root:myGlobals:gThreadGroupID=0
1112        endif
1113        NVAR mt=root:myGlobals:gThreadGroupID
1114
1115        if(mt!=0)       //there was an error with the stopping of the threads, possibly user abort
1116                err = ThreadGroupRelease(mt)
1117                Print "threading err = ",err
1118                if(err == 0)
1119                        // all *should* be OK
1120                else
1121                        return(0)
1122                endif
1123        endif
1124
1125        NVAR imon = root:Packages:NIST:SAS:gImon
1126        NVAR thick = root:Packages:NIST:SAS:gThick
1127        NVAR sig_incoh = root:Packages:NIST:SAS:gSig_incoh
1128        NVAR r2 = root:Packages:NIST:SAS:gR2
1129
1130        // do the simulation here, or not
1131        Variable r1,xCtr,yCtr,sdd,pixSize,wavelength
1132        String coefStr,abortStr,str
1133
1134        r1 = rw[24]/2/10                // sample diameter convert diam in [mm] to radius in cm
1135        xCtr = rw[16]
1136        yCtr = rw[17]
1137        sdd = rw[18]*100                //conver header of [m] to [cm]
1138        pixSize = rw[10]/10             // convert pix size in mm to cm
1139        wavelength = rw[26]
1140        coefStr = MC_getFunctionCoef(funcStr)
1141       
1142        if(!MC_CheckFunctionAndCoef(funcStr,coefStr))
1143                doMonteCarlo = 0                //we're getting out now, reset the flag so we don't continually end up here
1144                Abort "The coefficients and function type do not match. Please correct the selections in the popup menus."
1145        endif
1146       
1147        Variable sig_sas
1148//              FUNCREF SANSModelAAO_MCproto func=$("fSmeared"+funcStr)         //a wrapper for the structure version
1149        FUNCREF SANSModelAAO_MCproto func=$(funcStr)            //unsmeared
1150        WAVE results = root:Packages:NIST:SAS:results
1151        WAVE linear_data = root:Packages:NIST:SAS:linear_data
1152        WAVE data = root:Packages:NIST:SAS:data
1153
1154        results = 0
1155        linear_data = 0
1156       
1157        CalculateRandomDeviate(func,$coefStr,wavelength,"root:Packages:NIST:SAS:ran_dev",SIG_SAS)
1158        if(sig_sas > 100)
1159                sprintf abortStr,"sig_sas = %g. Please check that the model coefficients have a zero background, or the low q is well-behaved.",sig_sas
1160                Abort abortStr
1161        endif
1162       
1163        WAVE ran_dev=$"root:Packages:NIST:SAS:ran_dev"
1164       
1165        Make/O/D/N=5000 root:Packages:NIST:SAS:nt=0,root:Packages:NIST:SAS:j1=0,root:Packages:NIST:SAS:j2=0
1166        Make/O/D/N=100 root:Packages:NIST:SAS:nn=0
1167        Make/O/D/N=11 root:Packages:NIST:SAS:inputWave=0
1168       
1169        WAVE nt = root:Packages:NIST:SAS:nt
1170        WAVE j1 = root:Packages:NIST:SAS:j1
1171        WAVE j2 = root:Packages:NIST:SAS:j2
1172        WAVE nn = root:Packages:NIST:SAS:nn
1173        WAVE inputWave = root:Packages:NIST:SAS:inputWave
1174
1175        inputWave[0] = imon
1176        inputWave[1] = r1
1177        inputWave[2] = r2
1178        inputWave[3] = xCtr
1179        inputWave[4] = yCtr
1180        inputWave[5] = sdd
1181        inputWave[6] = pixSize
1182        inputWave[7] = thick
1183        inputWave[8] = wavelength
1184        inputWave[9] = sig_incoh
1185        inputWave[10] = sig_sas
1186
1187        linear_data = 0         //initialize
1188
1189        Variable t0,trans
1190       
1191        // get a time estimate, and give the user a chance to exit if they're unsure.
1192        t0 = stopMStimer(-2)
1193        inputWave[0] = 1000
1194        NVAR useXOP = root:Packages:NIST:SAS:gUse_MC_XOP                //if zero, will use non-threaded Igor code
1195       
1196        if(useXOP)
1197                //use a single thread, otherwise time is dominated by overhead
1198                Monte_SANS_NotThreaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
1199        else
1200                Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
1201        endif
1202       
1203        t0 = (stopMSTimer(-2) - t0)*1e-6
1204        t0 *= imon/1000/ThreadProcessorCount                    //projected time, in seconds (using threads for the calculation)
1205        inputWave[0] = imon             //reset
1206       
1207        if(t0>10)
1208                sprintf str,"The simulation will take approximately %d seconds.\r- Proceed?",t0
1209                DoAlert 1,str
1210                if(V_flag == 2)
1211                        doMonteCarlo = 0
1212                        reCalculateInten(1)             //come back in and do the smeared calculation
1213                        return(0)
1214                endif
1215        endif
1216       
1217        linear_data = 0         //initialize
1218// threading crashes!! - there must be some operation in the XOP that is not threadSafe. What, I don't know...
1219// I think it's the ran() calls, being "non-reentrant". So the XOP now defines two separate functions, that each
1220// use a different rng. This works. 1.75x speedup.     
1221        t0 = stopMStimer(-2)
1222
1223        if(useXOP)
1224                Monte_SANS_Threaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
1225        else
1226                Monte_SANS_NotThreaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
1227        endif
1228       
1229        t0 = (stopMSTimer(-2) - t0)*1e-6
1230        Printf  "MC sim time = %g seconds\r",t0
1231       
1232        trans = results[8]                      //(n1-n2)/n1
1233        if(trans == 0)
1234                trans = 1
1235        endif
1236
1237        Print "counts on detector, including transmitted = ",sum(linear_data,-inf,inf)
1238       
1239//              linear_data[xCtr][yCtr] = 0                     //snip out the transmitted spike
1240//              Print "counts on detector not transmitted = ",sum(linear_data,-inf,inf)
1241
1242        // or simulate a beamstop
1243        NVAR MC_BS_in = root:Packages:NIST:SAS:gBeamStopIn              //if zero, beam stop is "out", as in a transmission measurement
1244       
1245        Variable rad=beamstopDiam()/2           //beamstop radius in cm
1246        if(MC_BS_in)
1247                rad /= 0.5                              //convert cm to pixels
1248                rad += 0.                                       // (no - it cuts off the low Q artificially) add an extra pixel to each side to account for edge
1249                Duplicate/O linear_data,root:Packages:NIST:SAS:tmp_mask//,root:Packages:NIST:SAS:MC_linear_data
1250                WAVE tmp_mask = root:Packages:NIST:SAS:tmp_mask
1251                tmp_mask = (sqrt((p-xCtr)^2+(q-yCtr)^2) < rad) ? 0 : 1          //behind beamstop = 0, away = 1
1252               
1253                linear_data *= tmp_mask
1254        endif
1255       
1256        results[9] = sum(linear_data,-inf,inf)
1257        //              Print "counts on detector not behind beamstop = ",results[9]
1258       
1259        // convert to absolute scale
1260        Variable kappa          //,beaminten = beamIntensity()
1261//              kappa = beamInten*pi*r1*r1*thick*(pixSize/sdd)^2*trans*(iMon/beaminten)
1262        kappa = thick*(pixSize/sdd)^2*trans*iMon
1263       
1264        //use kappa to get back to counts => linear_data = round(linear_data*kappa)
1265        Note/K linear_data ,"KAPPA="+num2str(kappa)+";"
1266       
1267        NVAR rawCts = root:Packages:NIST:SAS:gRawCounts
1268        if(!rawCts)                     //go ahead and do the linear scaling
1269                linear_data = linear_data / kappa
1270        endif           
1271        data = linear_data
1272       
1273        // re-average the 2D data
1274        S_CircularAverageTo1D("SAS")
1275       
1276        // put the new result into the simulation folder
1277        Fake1DDataFolder(qval,aveint,sigave,sigmaQ,qbar,fSubs,"Simulation")     
1278                               
1279
1280        return(0)
1281end
1282
1283//phi is defined from +x axis, proceeding CCW around [0,2Pi]
1284ThreadSafe Function MC_FindPhi(vx,vy)
1285        variable vx,vy
1286       
1287        variable phi
1288       
1289        phi = atan(vy/vx)               //returns a value from -pi/2 to pi/2
1290       
1291        // special cases
1292        if(vx==0 && vy > 0)
1293                return(pi/2)
1294        endif
1295        if(vx==0 && vy < 0)
1296                return(3*pi/2)
1297        endif
1298        if(vx >= 0 && vy == 0)
1299                return(0)
1300        endif
1301        if(vx < 0 && vy == 0)
1302                return(pi)
1303        endif
1304       
1305       
1306        if(vx > 0 && vy > 0)
1307                return(phi)
1308        endif
1309        if(vx < 0 && vy > 0)
1310                return(phi + pi)
1311        endif
1312        if(vx < 0 && vy < 0)
1313                return(phi + pi)
1314        endif
1315        if( vx > 0 && vy < 0)
1316                return(phi + 2*pi)
1317        endif
1318       
1319        return(phi)
1320end
1321
1322
1323
1324
1325
1326Window Sim_1D_Panel() : Panel
1327        PauseUpdate; Silent 1           // building window...
1328        NewPanel /W=(92,556,713,818)/K=1 as "1D SANS Simulator"
1329        SetVariable cntVar,pos={26,68},size={160,15},title="Counting time(s)",format="%d"
1330        SetVariable cntVar,limits={1,36000,10},value= root:Packages:NIST:SAS:gCntTime
1331        SetVariable cntVar, proc=Sim_1D_CountTimeSetVarProc
1332        SetVariable MC_setvar0_1,pos={26,91},size={160,15},title="Thickness (cm)"
1333        SetVariable MC_setvar0_1,limits={0,inf,0.1},value= root:Packages:NIST:SAS:gThick
1334        SetVariable MC_setvar0_1, proc=Sim_1D_SamThickSetVarProc
1335
1336        SetVariable MC_setvar0_3,pos={26,114},size={160,15},title="Sample Transmission"
1337        SetVariable MC_setvar0_3,limits={0,1,0.01},value= root:Packages:NIST:SAS:gSamTrans
1338        SetVariable MC_setvar0_3, proc=Sim_1D_SamTransSetVarProc
1339
1340        PopupMenu MC_popup0,pos={13,13},size={165,20},proc=Sim_1D_ModelPopMenuProc,title="Model Function"
1341        PopupMenu MC_popup0,mode=1,value= #"MC_FunctionPopupList()"
1342        Button MC_button0,pos={17,181},size={130,20},proc=Sim_1D_DoItButtonProc,title="Do 1D Simulation"
1343        Button MC_button0,fColor=(3,52428,1)
1344        Button MC_button1,pos={17,211},size={150,20},proc=Save_1DSimData,title="Save Simulated Data"
1345        GroupBox group0,pos={15,42},size={280,130},title="Sample Setup"
1346        CheckBox check0_1,pos={216,179},size={60,14},title="Yes Offset",variable= root:Packages:NIST:SAS:gDoTraceOffset
1347        CheckBox check0_2,pos={216,199},size={60,14},title="Abs scale?",variable= root:Packages:NIST:SAS:g_1D_DoABS
1348        CheckBox check0_3,pos={216,219},size={60,14},title="Noise?",variable= root:Packages:NIST:SAS:g_1D_AddNoise
1349       
1350// a box for the results
1351        GroupBox group1,pos={314,23},size={277,163},title="Simulation Results"
1352        ValDisplay valdisp0,pos={326,48},size={220,13},title="Total detector counts"
1353        ValDisplay valdisp0,limits={0,0,0},barmisc={0,1000},value= root:Packages:NIST:SAS:g_1DTotCts
1354        ValDisplay valdisp0_1,pos={326,72},size={220,13},title="Estimated count rate (1/s)"
1355        ValDisplay valdisp0_1,limits={0,0,0},barmisc={0,1000},value=root:Packages:NIST:SAS:g_1DEstDetCR
1356        ValDisplay valdisp0_2,pos={326,96},size={220,13},title="Fraction of beam scattered"
1357        ValDisplay valdisp0_2,limits={0,0,0},barmisc={0,1000},value= root:Packages:NIST:SAS:g_1DFracScatt
1358        ValDisplay valdisp0_3,pos={326,121},size={220,13},title="Estimated transmission"
1359        ValDisplay valdisp0_3,limits={0,0,0},barmisc={0,1000},value=root:Packages:NIST:SAS:g_1DEstTrans
1360        ValDisplay valdisp0_4,pos={326,145},size={220,13},title="Multiple Coherent Scattering"
1361        ValDisplay valdisp0_4,limits={0,0,0},barmisc={0,1000},value=root:Packages:NIST:SAS:g_MultScattFraction
1362        // set the flags here -- do the simulation, but not 2D
1363       
1364        root:Packages:NIST:SAS:doSimulation     = 1     // == 1 if 1D simulated data, 0 if other from the checkbox
1365        root:Packages:NIST:SAS:gDoMonteCarlo     = 0  // == 1 if 2D MonteCarlo set by hidden flag
1366
1367       
1368EndMacro
1369
1370Function Sim_1D_CountTimeSetVarProc(sva) : SetVariableControl
1371        STRUCT WMSetVariableAction &sva
1372
1373        switch( sva.eventCode )
1374                case 1: // mouse up
1375                case 2: // Enter key
1376                case 3: // Live update
1377                        Variable dval = sva.dval
1378
1379                        ReCalculateInten(1)
1380                       
1381                        break
1382        endswitch
1383
1384        return 0
1385End
1386
1387Function Sim_1D_SamThickSetVarProc(sva) : SetVariableControl
1388        STRUCT WMSetVariableAction &sva
1389
1390        switch( sva.eventCode )
1391                case 1: // mouse up
1392                case 2: // Enter key
1393                case 3: // Live update
1394                        Variable dval = sva.dval
1395
1396                        ReCalculateInten(1)
1397                       
1398                        break
1399        endswitch
1400
1401        return 0
1402End
1403
1404Function Sim_1D_SamTransSetVarProc(sva) : SetVariableControl
1405        STRUCT WMSetVariableAction &sva
1406
1407        switch( sva.eventCode )
1408                case 1: // mouse up
1409                case 2: // Enter key
1410                case 3: // Live update
1411                        Variable dval = sva.dval
1412
1413                        ReCalculateInten(1)
1414                       
1415                        break
1416        endswitch
1417
1418        return 0
1419End
1420
1421
1422Function Sim_1D_ModelPopMenuProc(pa) : PopupMenuControl
1423        STRUCT WMPopupAction &pa
1424
1425        switch( pa.eventCode )
1426                case 2: // mouse up
1427                        Variable popNum = pa.popNum
1428                        String popStr = pa.popStr
1429                        SVAR gStr = root:Packages:NIST:SAS:gFuncStr
1430                        gStr = popStr
1431                       
1432                        break
1433        endswitch
1434
1435        return 0
1436End
1437
1438
1439Function Sim_1D_DoItButtonProc(ba) : ButtonControl
1440        STRUCT WMButtonAction &ba
1441
1442        switch( ba.eventCode )
1443                case 2: // mouse up
1444               
1445                        ReCalculateInten(1)
1446                       
1447                        break
1448        endswitch
1449
1450        return 0
1451End
1452
1453
1454//
1455//
1456//
1457Function Save_1DSimData(ctrlName) : ButtonControl
1458        String ctrlName
1459
1460        String type="SAS",fullpath=""
1461        Variable dialog=1               //=1 will present dialog for name
1462       
1463        String destStr=""
1464        destStr = "root:Packages:NIST:"+type
1465       
1466        Variable refNum
1467        String formatStr = "%15.4g %15.4g %15.4g %15.4g %15.4g %15.4g\r\n"
1468        String fname,ave="C",hdrStr1="",hdrStr2=""
1469        Variable step=1
1470       
1471        If(1)
1472                //setup a "fake protocol" wave, sice I have no idea of the current state of the data
1473                Make/O/T/N=8 root:myGlobals:Protocols:SIMProtocol
1474                Wave/T SIMProtocol = $"root:myGlobals:Protocols:SIMProtocol"
1475                String junk="****SIMULATED DATA****"
1476                //stick in the fake protocol...
1477                NVAR ctTime = root:Packages:NIST:SAS:gCntTime
1478                NVAR totalCts = root:Packages:NIST:SAS:g_1DTotCts                       //summed counts (simulated)
1479                NVAR detCR = root:Packages:NIST:SAS:g_1DEstDetCR                // estimated detector count rate
1480                NVAR fractScat = root:Packages:NIST:SAS:g_1DFracScatt
1481                NVAR mScat = root:Packages:NIST:SAS:g_MultScattFraction
1482       
1483                SIMProtocol[0] = junk
1484                SIMProtocol[1] = "\tCounting time (s) = "+num2str(ctTime)
1485                SIMProtocol[2] = "\tTotal detector counts = "+num2str(totalCts)
1486                SIMProtocol[3] = "\tDetector countrate (1/s) = "+num2str(detCR)
1487                SIMProtocol[4] = "\tFraction of beam scattered coherently = "+num2str(fractScat)
1488                SIMProtocol[5] = "\tFraction of multiple coherent scattering = "+num2str(mScat)
1489                SIMProtocol[6] = ""
1490                SIMProtocol[7] = ""
1491                //set the global
1492                String/G root:myGlobals:Protocols:gProtoStr = "SIMProtocol"
1493        Endif
1494       
1495       
1496        //*****these waves MUST EXIST, or IGOR Pro will crash, with a type 2 error****
1497        WAVE intw=$(destStr + ":integersRead")
1498        WAVE rw=$(destStr + ":realsRead")
1499        WAVE/T textw=$(destStr + ":textRead")
1500        WAVE qvals =$(destStr + ":qval")
1501        WAVE inten=$(destStr + ":aveint")
1502        WAVE sig=$(destStr + ":sigave")
1503        WAVE qbar = $(destStr + ":QBar")
1504        WAVE sigmaq = $(destStr + ":SigmaQ")
1505        WAVE fsubs = $(destStr + ":fSubS")
1506
1507        SVAR gProtoStr = root:myGlobals:Protocols:gProtoStr
1508        Wave/T proto=$("root:myGlobals:Protocols:"+gProtoStr)
1509       
1510        //check each wave
1511        If(!(WaveExists(intw)))
1512                Abort "intw DNExist Save_1DSimData()"
1513        Endif
1514        If(!(WaveExists(rw)))
1515                Abort "rw DNExist Save_1DSimData()"
1516        Endif
1517        If(!(WaveExists(textw)))
1518                Abort "textw DNExist Save_1DSimData()"
1519        Endif
1520        If(!(WaveExists(qvals)))
1521                Abort "qvals DNExist Save_1DSimData()"
1522        Endif
1523        If(!(WaveExists(inten)))
1524                Abort "inten DNExist Save_1DSimData()"
1525        Endif
1526        If(!(WaveExists(sig)))
1527                Abort "sig DNExist Save_1DSimData()"
1528        Endif
1529        If(!(WaveExists(qbar)))
1530                Abort "qbar DNExist Save_1DSimData()"
1531        Endif
1532        If(!(WaveExists(sigmaq)))
1533                Abort "sigmaq DNExist Save_1DSimData()"
1534        Endif
1535        If(!(WaveExists(fsubs)))
1536                Abort "fsubs DNExist Save_1DSimData()"
1537        Endif
1538        If(!(WaveExists(proto)))
1539                Abort "current protocol wave DNExist Save_1DSimData()"
1540        Endif
1541
1542        //strings can be too long to print-- must trim to 255 chars
1543        Variable ii,num=8
1544        Make/O/T/N=(num) tempShortProto
1545        for(ii=0;ii<num;ii+=1)
1546                tempShortProto[ii] = (proto[ii])[0,240]
1547        endfor
1548       
1549        if(dialog)
1550                PathInfo/S catPathName
1551                fullPath = DoSaveFileDialog("Save data as")
1552                If(cmpstr(fullPath,"")==0)
1553                        //user cancel, don't write out a file
1554                        Close/A
1555                        Abort "no data file was written"
1556                Endif
1557                //Print "dialog fullpath = ",fullpath
1558        Endif
1559       
1560        NVAR monCt = root:Packages:NIST:SAS:gImon
1561        NVAR thick = root:Packages:NIST:SAS:gThick
1562        NVAR trans = root:Packages:NIST:SAS:gSamTrans                   //for 1D, default value
1563       
1564
1565       
1566        hdrStr1 = num2str(monCt)+"  "+num2str(rw[26])+"       "+num2str(rw[19])+"     "+num2str(rw[18])
1567        hdrStr1 += "     "+num2str(trans)+"     "+num2str(thick) + ave +"   "+num2str(step) + "\r\n"
1568
1569        hdrStr2 = num2str(rw[16])+"  "+num2str(rw[17])+"  "+num2str(rw[23])+"    "+num2str(rw[24])+"    "
1570        hdrStr2 += num2str(rw[25])+"    "+num2str(rw[27])+"    "+num2str(rw[21])+"    "+"ORNL  " + "\r\n"
1571       
1572        //actually open the file here
1573        Open refNum as fullpath
1574       
1575        //write out the standard header information
1576        fprintf refnum,"FILE: %s\t\t CREATED: %s\r\n","SIMULATED DATA",(date() +"  "+ time())
1577        fprintf refnum,"LABEL: %s\r\n","SIMULATED DATA"
1578        fprintf refnum,"MON CNT   LAMBDA   DET ANG   DET DIST   TRANS   THICK   AVE   STEP\r\n"
1579        fprintf refnum,hdrStr1
1580        fprintf refnum,"BCENT(X,Y)   A1(mm)   A2(mm)   A1A2DIST(m)   DL/L   BSTOP(mm)   DET_TYP \r\n"
1581        fprintf refnum,hdrStr2
1582//      fprintf refnum,headerFormat,rw[0],rw[26],rw[19],rw[18],rw[4],rw[5],ave,step
1583
1584        //insert protocol information here
1585        //-1 list of sample files
1586        //0 - bkg
1587        //1 - emp
1588        //2 - div
1589        //3 - mask
1590        //4 - abs params c2-c5
1591        //5 - average params
1592        fprintf refnum, "SAM: %s\r\n",tempShortProto[0]
1593        fprintf refnum, "BGD: %s\r\n",tempShortProto[1]
1594        fprintf refnum, "EMP: %s\r\n",tempShortProto[2]
1595        fprintf refnum, "DIV: %s\r\n",tempShortProto[3]
1596        fprintf refnum, "MASK: %s\r\n",tempShortProto[4]
1597        fprintf refnum, "ABS: %s\r\n",tempShortProto[5]
1598        fprintf refnum, "Average Choices: %s\r\n",tempShortProto[6]
1599       
1600        //write out the data columns
1601        fprintf refnum,"The 6 columns are | Q (1/A) | I(Q) (1/cm) | std. dev. I(Q) (1/cm) | sigmaQ | meanQ | ShadowFactor|\r\n"
1602        wfprintf refnum, formatStr, qvals,inten,sig,sigmaq,qbar,fsubs
1603       
1604        Close refnum
1605       
1606        SetDataFolder root:             //(redundant)
1607       
1608        //write confirmation of write operation to history area
1609        Print "Averaged File written: ", GetFileNameFromPathNoSemi(fullPath)
1610        KillWaves/Z tempShortProto
1611
1612        //clear the stuff that was created for case of saving files
1613        If(1)
1614                Killwaves/Z root:myGlobals:Protocols:SIMProtocol
1615                String/G root:myGlobals:Protocols:gProtoStr = ""
1616        Endif
1617       
1618       
1619        return(0)
1620       
1621End
1622
1623
1624/// called in SASCALC:ReCalculateInten()
1625Function Simulate_1D(funcStr,aveint,qval,sigave,sigmaq,qbar,fsubs)
1626        String funcStr
1627        WAVE aveint,qval,sigave,sigmaq,qbar,fsubs
1628
1629        Variable r1,xCtr,yCtr,sdd,pixSize,wavelength
1630        String coefStr,abortStr,str     
1631
1632        FUNCREF SANSModelAAO_MCproto func=$("fSmeared"+funcStr)                 //a wrapper for the structure version
1633        FUNCREF SANSModelAAO_MCproto funcUnsmeared=$(funcStr)           //unsmeared
1634        coefStr = MC_getFunctionCoef(funcStr)
1635       
1636        if(!MC_CheckFunctionAndCoef(funcStr,coefStr))
1637                Abort "Function and coefficients do not match. You must plot the unsmeared function before simulation."
1638        endif
1639       
1640        Wave inten=$"root:Simulation:Simulation_i"              // this will exist and send the smeared calculation to the corect DF
1641       
1642        // the resolution-smeared intensity is calculated, including the incoherent background
1643        func($coefStr,inten,qval)
1644
1645        NVAR imon = root:Packages:NIST:SAS:gImon
1646        NVAR ctTime = root:Packages:NIST:SAS:gCntTime
1647        NVAR thick = root:Packages:NIST:SAS:gThick
1648        NVAR trans = root:Packages:NIST:SAS:gSamTrans
1649        NVAR SimDetCts = root:Packages:NIST:SAS:g_1DTotCts                      //summed counts (simulated)
1650        NVAR estDetCR = root:Packages:NIST:SAS:g_1DEstDetCR                     // estimated detector count rate
1651        NVAR fracScat = root:Packages:NIST:SAS:g_1DFracScatt            // fraction of beam captured on detector
1652        NVAR estTrans = root:Packages:NIST:SAS:g_1DEstTrans             // estimated transmission of sample
1653        NVAR mScat = root:Packages:NIST:SAS:g_MultScattFraction
1654       
1655        WAVE rw=root:Packages:NIST:SAS:realsRead
1656        WAVE nCells=root:Packages:NIST:SAS:nCells                               
1657                                       
1658        pixSize = rw[10]/10             // convert pix size in mm to cm
1659        sdd = rw[18]*100                //convert header of [m] to [cm]
1660        wavelength = rw[26]             // in 1/A
1661       
1662        imon = beamIntensity()*ctTime
1663       
1664        // calculate the scattering cross section simply to be able to estimate the transmission
1665        Variable sig_sas=0
1666       
1667        // remember that the random deviate is the coherent portion ONLY - the incoherent background is
1668        // subtracted before the calculation.
1669        CalculateRandomDeviate(funcUnsmeared,$coefStr,wavelength,"root:Packages:NIST:SAS:ran_dev",sig_sas)
1670       
1671//                              if(sig_sas > 100)
1672//                                      sprintf abortStr,"sig_sas = %g. Please check that the model coefficients have a zero background, or the low q is well-behaved.",sig_sas
1673//                              endif
1674
1675        // calculate the multiple scattering fraction for display (10/2009)
1676        Variable ii,nMax=10,tau
1677        mScat=0
1678        tau = thick*sig_sas
1679        // this sums the normalized scattering P', so the result is the fraction of multiply coherently scattered
1680        // neutrons out of those that were scattered
1681        for(ii=2;ii<nMax;ii+=1)
1682                mScat += tau^(ii)/factorial(ii)
1683//              print tau^(ii)/factorial(ii)
1684        endfor
1685        estTrans = exp(-1*thick*sig_sas)                //thickness and sigma both in units of cm
1686        mscat *= (estTrans)/(1-estTrans)
1687
1688        if(mScat > 0.1)         //  /Z to supress error if this was a 1D calc with the 2D panel open
1689                ValDisplay/Z valdisp0_4 win=Sim_1D_Panel,labelBack=(65535,32768,32768)
1690        else
1691                ValDisplay/Z valdisp0_4 win=Sim_1D_Panel,labelBack=0
1692        endif
1693
1694
1695        Print "Sig_sas = ",sig_sas
1696       
1697        Duplicate/O qval prob_i,countsInAnnulus
1698       
1699        // not needed - nCells takes care of this when the error is correctly calculated
1700//                              Duplicate/O qval circle_fraction,rval,nCells_expected
1701//                              rval = sdd*tan(2*asin(qval*wavelength/4/pi))            //radial distance in cm
1702//                              nCells_expected = 2*pi*rval/pixSize                                     //does this need to be an integer?
1703//                              circle_fraction = nCells / nCells_expected
1704       
1705                               
1706//                              prob_i = trans*thick*nCells*(pixSize/sdd)^2*inten                       //probability of a neutron in q-bin(i) that has nCells
1707        prob_i = trans*thick*(pixSize/sdd)^2*inten                      //probability of a neutron in q-bin(i)
1708       
1709        Variable P_on = sum(prob_i,-inf,inf)
1710        Print "P_on = ",P_on
1711       
1712//                              fracScat = P_on
1713        fracScat = 1-estTrans
1714       
1715//                              aveint = (Imon)*prob_i / circle_fraction / nCells_expected
1716        aveint = (Imon)*prob_i
1717
1718        countsInAnnulus = aveint*nCells
1719        SimDetCts = sum(countsInAnnulus,-inf,inf)
1720        estDetCR = SimDetCts/ctTime
1721       
1722       
1723        NVAR doABS = root:Packages:NIST:SAS:g_1D_DoABS
1724        NVAR addNoise = root:Packages:NIST:SAS:g_1D_AddNoise
1725       
1726        // this is where the number of cells comes in - the calculation of the error bars
1727        // sigma[i] = SUM(sigma[ij]^2) / nCells^2
1728        // and since in the simulation, SUM(sigma[ij]^2) = nCells*sigma[ij]^2 = nCells*Inten
1729        // then...
1730        sigave = sqrt(aveint/nCells)            // corrected based on John's memo, from 8/9/99
1731       
1732        // add in random error in aveint based on the sigave
1733        if(addNoise)
1734                aveint += gnoise(sigave)
1735        endif
1736
1737        // signature in the standard deviation, do this after the noise is added
1738        // start at 10 to be out of the beamstop (makes for nicer plotting)
1739        // end at 50 to leave the natural statistics at the end of the set (may have a total of 80+ points if no offset)
1740        sigave[10,50;10] = 10*sigave[p]
1741
1742        // convert to absolute scale
1743        if(doABS)
1744                Variable kappa = thick*(pixSize/sdd)^2*trans*iMon
1745                aveint /= kappa
1746                sigave /= kappa
1747        endif
1748                               
1749                               
1750        return(0)
1751End
1752
1753/// for testing only
1754Function testProbability(sas,thick)
1755        Variable sas,thick
1756       
1757        Variable tau,trans,p2p,p3p,p4p
1758       
1759        tau = sas*thick
1760        trans = exp(-tau)
1761       
1762        Print "tau = ",tau
1763        Print "trans = ",trans
1764       
1765        p2p = tau^2/factorial(2)*trans/(1-trans)
1766        p3p = tau^3/factorial(3)*trans/(1-trans)
1767        p4p = tau^4/factorial(4)*trans/(1-trans)
1768       
1769        Print "double scattering = ",p2p
1770        Print "triple scattering = ",p3p
1771        Print "quadruple scattering = ",p4p
1772       
1773        Variable ii,nMax=10,mScat=0
1774        for(ii=2;ii<nMax;ii+=1)
1775                mScat += tau^(ii)/factorial(ii)
1776//              print tau^(ii)/factorial(ii)
1777        endfor
1778        mscat *= (Trans)/(1-Trans)
1779
1780        Print "Total fraction of multiple scattering = ",mScat
1781
1782        return(mScat)
1783End
Note: See TracBrowser for help on using the repository browser.