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

Last change on this file since 456 was 456, checked in by srkline, 14 years ago

minor changes to MonteCarlo? procedures to add in a realistic beamstop, plus some tests of including absorption and detector efficiency to the simulation. These have been commented out as is seems that they do little to improve the agreement betwen Mc countrate and real data countrate. MC simulation is 2.7 or 3.7x higher CR on detector (only two data points). Need to verify flux estimates first.

File size: 34.5 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2
3//
4// Monte Carlo simulator for SASCALC
5// October 2008 SRK
6//
7// This code simulates the scattering for a selected model based on the instrument configuration
8// This code is based directly on John Barker's code, which includes multiple scattering effects.
9// A lot of the setup, wave creation, and post-calculations are done in SASCALC->ReCalculateInten()
10//
11//
12// - Why am I off by a factor of 2.7 - 3.7 (MC too high) relative to real data?
13//   I need to include efficiency (70%?) - do I knock these off be fore the simulation or do I
14//    really simulate that some fraction of neutrons on the detector don't actually get counted?
15//   Is the flux estimate up-to-date?
16// - Most importantly, this needs to be checked for correctness of the MC simulation
17// X how can I get the "data" on absolute scale? This would be a great comparison vs. the ideal model calculation
18// X why does my integrated tau not match up with John's analytical calculations? where are the assumptions?
19// - get rid of all small angle assumptions - to make sure that the calculation is correct at all angles
20// - my simulated transmission is larger than what is measured, even after correcting for the quartz cell.
21//   Why? Do I need to include absorption? Just inherent problems with incoherent cross sections?
22//
23// X at the larger angles, is the "flat" detector being properly accounted for - in terms of
24//   the solid angle and how many counts fall in that pixel. Am I implicitly defining a spherical detector
25//   so that what I see is already "corrected"?
26// X the MC will, of course benefit greatly from being XOPized. Maybe think about parallel implementation
27//   by allowing the data arrays to accumulate. First pass at the XOP is done. Not pretty, not the speediest (5.8x)
28//   but it is functional. Add spinCursor() for long calculations. See WaveAccess XOP example.
29// X the background parameter for the model MUST be zero, or the integration for scattering
30//    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
31// X fully use the SASCALC input, most importantly, flux on sample.
32// X if no MC desired, still use the selected model
33// X better display of MC results on panel
34// X settings for "count for X seconds" or "how long to 1E6 cts on detector" (but 1E6 is typically too many counts...)
35// X warn of projected simulation time
36// - add quartz window scattering to the simulation somehow
37// -?- do smeared models make any sense?? Yes, John agrees that they do, and may be used in a more realistic simulation
38//   -?- but the random deviate can't be properly calculated...
39// - make sure that the ratio of scattering coherent/incoherent is properly adjusted for the sample composition
40//   or the volume fraction of solvent.
41//
42// X add to the results the fraction of coherently scattered neutrons that are singly scattered, different than
43//   the overall fraction of singly scattered, and maybe more important to know.
44//
45// X change the fraction reaching the detector to exclude those that don't interact. These transmitted neutrons
46//   aren't counted. Is the # that interact a better number?
47//
48// - do we want to NOT offset the data by a multiplicative factor as it is "frozen" , so that the
49//   effects on the absolute scale can be seen?
50//
51// X why is "pure" incoherent scattering giving me a q^-1 slope, even with the detector all the way back?
52// - can I speed up by assuming everything interacts? This would compromise the ability to calculate multiple scattering
53// X ask John how to verify what is going on
54// - a number of models are now found to be ill-behaved when q=1e-10. Then the random deviate calculation blows up.
55//   a warning has been added - but some models need a proper limiting value, and some (power-law) are simply unuseable
56//   unless something else can be done.
57//
58//
59
60// threaded call to the main function, adds up the individual runs, and returns what is to be displayed
61// results is calculated and sent back for display
62Function Monte_SANS_Threaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
63        WAVE inputWave,ran_dev,nt,j1,j2,nn,linear_data,results
64
65        //initialize ran1 in the XOP by passing a negative integer
66        // does nothing in the Igor code
67        Duplicate/O results retWave
68        //results[0] = -1*(datetime)
69
70        Variable NNeutron=inputWave[0]
71        Variable i,nthreads= ThreadProcessorCount
72        if(nthreads>2)          //only support 2 processors until I can figure out how to properly thread the XOP and to loop it
73                nthreads=2
74        endif
75       
76//      nthreads = 1
77       
78        variable mt= ThreadGroupCreate(nthreads)
79       
80        inputWave[0] = NNeutron/nthreads                //split up the number of neutrons
81       
82        for(i=0;i<nthreads;i+=1)
83                Duplicate/O nt $("nt"+num2istr(i))              //new instance for each thread
84                Duplicate/O j1 $("j1"+num2istr(i))
85                Duplicate/O j2 $("j2"+num2istr(i))
86                Duplicate/O nn $("nn"+num2istr(i))
87                Duplicate/O linear_data $("linear_data"+num2istr(i))
88                Duplicate/O retWave $("retWave"+num2istr(i))
89                Duplicate/O inputWave $("inputWave"+num2istr(i))
90                Duplicate/O ran_dev $("ran_dev"+num2istr(i))
91               
92                // ?? I need explicit wave references?
93                // maybe I need to have everything in separate data folders - bu tI haven't tried that. seems like a reach.
94                // more likely there is something bad going on in the XOP code.
95                if(i==0)
96                        WAVE inputWave0,ran_dev0,nt0,j10,j20,nn0,linear_data0,retWave0
97                        retWave0[0] = -1*datetime               //to initialize ran3
98                        ThreadStart mt,i,Monte_SANS_W1(inputWave0,ran_dev0,nt0,j10,j20,nn0,linear_data0,retWave0)
99                        Print "started thread 0"
100                endif
101                if(i==1)
102                        WAVE inputWave1,ran_dev1,nt1,j11,j21,nn1,linear_data1,retWave1
103                        //retWave1[0] = -1*datetime             //to initialize ran3
104                        ThreadStart mt,i,Monte_SANS_W1(inputWave1,ran_dev1,nt1,j11,j21,nn1,linear_data1,retWave1)
105                        Print "started thread 1"
106                endif
107//              if(i==2)
108//                      WAVE inputWave2,ran_dev2,nt2,j12,j22,nn2,linear_data2,retWave2
109//                      retWave2[0] = -1*datetime               //to initialize ran3
110//                      ThreadStart mt,i,Monte_SANS_W(inputWave2,ran_dev2,nt2,j12,j22,nn2,linear_data2,retWave2)
111//              endif
112//              if(i==3)
113//                      WAVE inputWave3,ran_dev3,nt3,j13,j23,nn3,linear_data3,retWave3
114//                      retWave3[0] = -1*datetime               //to initialize ran3
115//                      ThreadStart mt,i,Monte_SANS_W(inputWave3,ran_dev3,nt3,j13,j23,nn3,linear_data3,retWave3)
116//              endif
117        endfor
118
119// wait until done
120        do
121                variable tgs= ThreadGroupWait(mt,100)
122        while( tgs != 0 )
123        variable dummy= ThreadGroupRelease(mt)
124        Print "done with all threads"
125
126        // calculate all of the bits for the results
127        if(nthreads == 1)
128                nt = nt0                // add up each instance
129                j1 = j10
130                j2 = j20
131                nn = nn0
132                linear_data = linear_data0
133                retWave = retWave0
134        endif
135        if(nthreads == 2)
136                nt = nt0+nt1            // add up each instance
137                j1 = j10+j11
138                j2 = j20+j21
139                nn = nn0+nn1
140                linear_data = linear_data0+linear_data1
141                retWave = retWave0+retWave1
142        endif
143//      if(nthreads == 3)
144//              nt = nt0+nt1+nt2                // add up each instance
145//              j1 = j10+j11+j12
146//              j2 = j20+j21+j22
147//              nn = nn0+nn1+nn2
148//              linear_data = linear_data0+linear_data1+linear_data2
149//              retWave = retWave0+retWave1+retWave2
150//      endif
151//      if(nthreads == 4)
152//              nt = nt0+nt1+nt2+nt3            // add up each instance
153//              j1 = j10+j11+j12+j13
154//              j2 = j20+j21+j22+j23
155//              nn = nn0+nn1+nn2+nn3
156//              linear_data = linear_data0+linear_data1+linear_data2+linear_data3
157//              retWave = retWave0+retWave1+retWave2+retWave3
158//      endif
159       
160        // fill up the results wave
161        Variable xc,yc
162        xc=inputWave[3]
163        yc=inputWave[4]
164        results[0] = inputWave[9]+inputWave[10]         //total XS
165        results[1] = inputWave[10]                                              //SAS XS
166        results[2] = retWave[1]                                                 //number that interact n2
167        results[3] = retWave[2] - linear_data[xc][yc]                           //# reaching detector minus Q(0)
168        results[4] = retWave[3]/retWave[1]                              //avg# times scattered
169        results[5] = retWave[4]/retWave[1]                                              //single coherent fraction
170        results[6] = retWave[5]/retWave[1]                              //double coherent fraction
171        results[7] = retWave[6]/retWave[1]                              //multiple scatter fraction
172        results[8] = (retWave[0]-retWave[1])/retWave[0]                 //transmitted fraction
173       
174        return(0)
175End
176
177// worker function for threads, does nothing except switch between XOP and Igor versions
178ThreadSafe Function Monte_SANS_W1(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
179        WAVE inputWave,ran_dev,nt,j1,j2,nn,linear_data,results
180       
181#if exists("Monte_SANSX")
182        Monte_SANSX(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
183#else
184        Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
185#endif
186
187        return (0)
188End
189// worker function for threads, does nothing except switch between XOP and Igor versions
190ThreadSafe Function Monte_SANS_W2(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
191        WAVE inputWave,ran_dev,nt,j1,j2,nn,linear_data,results
192       
193#if exists("xxxxMonte_SANSX")
194        Monte_SANSX(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
195#else
196        Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
197#endif
198
199        return (0)
200End
201
202// NON-threaded call to the main function returns what is to be displayed
203// results is calculated and sent back for display
204Function Monte_SANS_NotThreaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
205        WAVE inputWave,ran_dev,nt,j1,j2,nn,linear_data,results
206
207        //initialize ran1 in the XOP by passing a negative integer
208        // does nothing in the Igor code, enoise is already initialized
209        Duplicate/O results retWave
210        WAVE retWave
211        retWave[0] = -1*abs(trunc(100000*enoise(1)))
212       
213#if exists("Monte_SANSX")
214        Monte_SANSX(inputWave,ran_dev,nt,j1,j2,nn,linear_data,retWave)
215#else
216        Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,retWave)
217#endif
218
219        // fill up the results wave
220        Variable xc,yc
221        xc=inputWave[3]
222        yc=inputWave[4]
223        results[0] = inputWave[9]+inputWave[10]         //total XS
224        results[1] = inputWave[10]                                              //SAS XS
225        results[2] = retWave[1]                                                 //number that interact n2
226        results[3] = retWave[2] - linear_data[xc][yc]                           //# reaching detector minus Q(0)
227        results[4] = retWave[3]/retWave[1]                              //avg# times scattered
228        results[5] = retWave[4]/retWave[1]                                              //single coherent fraction
229        results[6] = retWave[5]/retWave[1]                              //double coherent fraction
230        results[7] = retWave[6]/retWave[1]                              //multiple scatter fraction
231        results[8] = (retWave[0]-retWave[1])/retWave[0]                 //transmitted fraction
232       
233        return(0)
234End
235
236
237
238//////////
239//    PROGRAM Monte_SANS
240//    PROGRAM simulates multiple SANS.
241//       revised 2/12/99  JGB
242//            added calculation of random deviate, and 2D 10/2008 SRK
243
244//    N1 = NUMBER OF INCIDENT NEUTRONS.
245//    N2 = NUMBER INTERACTED IN THE SAMPLE.
246//    N3 = NUMBER ABSORBED.
247//    THETA = SCATTERING ANGLE.
248
249//        IMON = 'Enter number of neutrons to use in simulation.'
250//        NUM_BINS = 'Enter number of THETA BINS TO use. (<5000).'
251//        R1 = 'Enter beam radius. (cm)'
252//        R2 = 'Enter sample radius. (cm)'
253//        thick = 'Enter sample thickness. (cm)'
254//        wavelength = 'Enter neutron wavelength. (A)'
255//        R0 = 'Enter sphere radius. (A)'
256//
257
258ThreadSafe Function Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,MC_linear_data,results)
259        WAVE inputWave,ran_dev,nt,j1,j2,nn,MC_linear_data,results
260
261        Variable imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,sig_sas
262        Variable NUM_BINS,N_INDEX
263        Variable RHO,SIGSAS,SIGABS_0
264        Variable ii,jj,IND,idum,INDEX,IR,NQ
265        Variable qmax,theta_max,R_DAB,R0,BOUND,I0,q0,zpow
266        Variable N1,N2,N3,DTH,zz,tt,SIG_SINGLE,xx,yy,PHI,UU,SIG
267        Variable THETA,Ran,ll,D_OMEGA,RR,Tabs,Ttot,I1_sumI
268        Variable G0,E_NT,E_NN,TRANS_th,Trans_exp,rat
269        Variable GG,GG_ED,dS_dW,ds_dw_double,ds_dw_single
270        Variable DONE,FIND_THETA,err            //used as logicals
271
272        Variable Vx,Vy,Vz,Theta_z,qq
273        Variable Sig_scat,Sig_abs,Ratio,Sig_total
274        Variable isOn=0,testQ,testPhi,xPixel,yPixel
275        Variable NSingleIncoherent,NSingleCoherent,NScatterEvents,incoherentEvent,coherentEvent
276        Variable NDoubleCoherent,NMultipleScatter,countIt,detEfficiency
277       
278        detEfficiency = 1.0             //70% counting efficiency = 0.7
279       
280        imon = inputWave[0]
281        r1 = inputWave[1]
282        r2 = inputWave[2]
283        xCtr = inputWave[3]
284        yCtr = inputWave[4]
285        sdd = inputWave[5]
286        pixSize = inputWave[6]
287        thick = inputWave[7]
288        wavelength = inputWave[8]
289        sig_incoh = inputWave[9]
290        sig_sas = inputWave[10]
291       
292//      SetRandomSeed 0.1               //to get a reproduceable sequence
293
294//scattering power and maximum qvalue to bin
295//      zpow = .1               //scattering power, calculated below
296        qmax = 4*pi/wavelength          //maximum Q to bin 1D data. (A-1) (not really used, so set to a big value)
297        sigabs_0 = 0.0          // ignore absorption cross section/wavelength [1/(cm A)]
298        N_INDEX = 50            // maximum number of scattering events per neutron
299        num_bins = 200          //number of 1-D bins (not really used)
300       
301// my additions - calculate the random deviate function as needed
302// and calculate the scattering power from the model function (passed in as a wave)
303//
304        Variable left = leftx(ran_dev)
305        Variable delta = deltax(ran_dev)
306       
307//c       total SAS cross-section
308//      SIG_SAS = zpow/thick
309        zpow = sig_sas*thick                    //since I now calculate the sig_sas from the model
310        SIG_ABS = SIGABS_0 * WAVElength
311        sig_abs = 0.0           //cm-1
312        SIG_TOTAL =SIG_ABS + SIG_SAS + sig_incoh
313//      Print "The TOTAL XSECTION. (CM-1) is ",sig_total
314//      Print "The TOTAL SAS XSECTION. (CM-1) is ",sig_sas
315//      results[0] = sig_total          //assign these after everything's done
316//      results[1] = sig_sas
317//      variable ratio1,ratio2
318//      ratio1 = sig_abs/sig_total
319//      ratio2 = sig_incoh/sig_total
320//      // 0->ratio1 = abs
321//      // ratio1 -> ratio2 = incoh
322//      // > ratio2 = coherent
323        RATIO = sig_incoh / SIG_TOTAL
324       
325//c       assuming theta = sin(theta)...OK
326        theta_max = wavelength*qmax/(2*pi)
327//C     SET Theta-STEP SIZE.
328        DTH = Theta_max/NUM_BINS
329//      Print "theta bin size = dth = ",dth
330
331//C     INITIALIZE COUNTERS.
332        N1 = 0
333        N2 = 0
334        N3 = 0
335        NSingleIncoherent = 0
336        NSingleCoherent = 0
337        NDoubleCoherent = 0
338        NMultipleScatter = 0
339        NScatterEvents = 0
340
341//C     INITIALIZE ARRAYS.
342        j1 = 0
343        j2 = 0
344        nt = 0
345        nn=0
346       
347//C     MONITOR LOOP - looping over the number of incedent neutrons
348//note that zz, is the z-position in the sample - NOT the scattering power
349
350// NOW, start the loop, throwing neutrons at the sample.
351        do
352                Vx = 0.0                        // Initialize direction vector.
353                Vy = 0.0
354                Vz = 1.0
355               
356                Theta = 0.0             //      Initialize scattering angle.
357                Phi = 0.0                       //      Intialize azimuthal angle.
358                N1 += 1                 //      Increment total number neutrons counter.
359                DONE = 0                        //      True when neutron is scattered out of the sample.
360                INDEX = 0                       //      Set counter for number of scattering events.
361                zz = 0.0                        //      Set entering dimension of sample.
362                incoherentEvent = 0
363                coherentEvent = 0
364               
365                do                                      //      Makes sure position is within circle.
366                        ran = abs(enoise(1))            //[0,1]
367                        xx = 2.0*R1*(Ran-0.5)           //X beam position of neutron entering sample.
368                        ran = abs(enoise(1))            //[0,1]
369                        yy = 2.0*R1*(Ran-0.5)           //Y beam position ...
370                        RR = SQRT(xx*xx+yy*yy)          //Radial position of neutron in incident beam.
371                while(rr>r1)
372
373                do    //Scattering Loop, will exit when "done" == 1
374                                // keep scattering multiple times until the neutron exits the sample
375                        ran = abs(enoise(1))            //[0,1]  RANDOM NUMBER FOR DETERMINING PATH LENGTH
376                        ll = PATH_len(ran,Sig_total)
377                        //Determine new scattering direction vector.
378                        err = NewDirection(vx,vy,vz,Theta,Phi)          //vx,vy,vz is updated, theta, phi unchanged by function
379
380                        //X,Y,Z-POSITION OF SCATTERING EVENT.
381                        xx += ll*vx
382                        yy += ll*vy
383                        zz += ll*vz
384                        RR = sqrt(xx*xx+yy*yy)          //radial position of scattering event.
385
386                        //Check whether interaction occurred within sample volume.
387                        IF (((zz > 0.0) && (zz < THICK)) && (rr < r2))
388                                //NEUTRON INTERACTED.
389                                ran = abs(enoise(1))            //[0,1]
390                               
391//                              if(ran<ratio1)
392//                                      //absorption event
393//                                      n3 +=1
394//                                      done=1
395//                              else
396
397                                INDEX += 1                      //Increment counter of scattering events.
398                                IF(INDEX == 1)
399                                        N2 += 1                 //Increment # of scat. neutrons
400                                Endif
401                                //Split neutron interactions into scattering and absorption events
402//                              IF(ran > (ratio1+ratio2) )              //C             NEUTRON SCATTERED coherently
403                                IF(ran > ratio)         //C             NEUTRON SCATTERED coherently
404                                        coherentEvent = 1
405                                        FIND_THETA = 0                  //false
406                                        DO
407                                                //ran = abs(enoise(1))          //[0,1]
408                                                //theta = Scat_angle(Ran,R_DAB,wavelength)      // CHOOSE DAB ANGLE -- this is 2Theta
409                                                //Q0 = 2*PI*THETA/WAVElength                                    // John chose theta, calculated Q
410
411                                                // pick a q-value from the deviate function
412                                                // pnt2x truncates the point to an integer before returning the x
413                                                // so get it from the wave scaling instead
414                                                Q0 =left + binarysearchinterp(ran_dev,abs(enoise(1)))*delta
415                                                theta = Q0/2/Pi*wavelength              //SAS approximation
416                                               
417                                                //Print "q0, theta = ",q0,theta
418                                               
419                                                FIND_THETA = 1          //always accept
420
421                                        while(!find_theta)
422                                        ran = abs(enoise(1))            //[0,1]
423                                        PHI = 2.0*PI*Ran                        //Chooses azimuthal scattering angle.
424                                ELSE
425                                        //NEUTRON scattered incoherently
426                   // N3 += 1
427                  // DONE = 1
428                  // phi and theta are random over the entire sphere of scattering
429                  // !can't just choose random theta and phi, won't be random over sphere solid angle
430                        incoherentEvent = 1
431                       
432                        ran = abs(enoise(1))            //[0,1]
433                                        theta = acos(2*ran-1)           
434                       
435                        ran = abs(enoise(1))            //[0,1]
436                                        PHI = 2.0*PI*Ran                        //Chooses azimuthal scattering angle.
437                                ENDIF           //(ran > ratio)
438//                              endif           // event was absorption
439                        ELSE
440                                //NEUTRON ESCAPES FROM SAMPLE -- bin it somewhere
441                                DONE = 1                //done = true, will exit from loop
442                               
443//                              countIt = 1
444//                              if(abs(enoise(1)) > detEfficiency)              //efficiency of 70% wired @top
445//                                      countIt = 0                                     //detector does not register
446//                              endif
447                               
448                                //Increment #scattering events array
449                                If (index <= N_Index)
450                                        NN[INDEX] += 1
451                                Endif
452                               
453                                if(index != 0)          //the neutron interacted at least once, figure out where it ends up
454
455                                        Theta_z = acos(Vz)              // Angle WITH respect to z axis.
456                                        testQ = 2*pi*sin(theta_z)/wavelength
457                                       
458                                        // pick a random phi angle, and see if it lands on the detector
459                                        // since the scattering is isotropic, I can safely pick a new, random value
460                                        // this would not be true if simulating anisotropic scattering.
461                                        //testPhi = abs(enoise(1))*2*Pi
462                                        testPhi = FindPhi(Vx,Vy)                //use the exiting phi value as defined by Vx and Vy
463                                       
464                                        // is it on the detector?       
465                                        FindPixel(testQ,testPhi,wavelength,sdd,pixSize,xCtr,yCtr,xPixel,yPixel)
466                                       
467                                        if(xPixel != -1 && yPixel != -1)
468                                                //if(index==1)  // only the single scattering events
469                                                        MC_linear_data[xPixel][yPixel] += 1             //this is the total scattering, including multiple scattering
470                                                //endif
471                                                        isOn += 1               // neutron that lands on detector
472                                        endif
473       
474                                        If(theta_z < theta_max)
475                                                //Choose index for scattering angle array.
476                                                //IND = NINT(THETA_z/DTH + 0.4999999)
477                                                ind = round(THETA_z/DTH + 0.4999999)            //round is eqivalent to nint()
478                                                NT[ind] += 1                    //Increment bin for angle.
479                                                //Increment angle array for single scattering events.
480                                                IF(INDEX == 1)
481                                                        j1[ind] += 1
482                                                Endif
483                                                //Increment angle array for double scattering events.
484                                                IF (INDEX == 2)
485                                                        j2[ind] += 1
486                                                Endif
487                                        EndIf
488                                       
489                                        // increment all of the counters now since done==1 here and I'm sure to exit and get another neutron
490                                        NScatterEvents += index         //total number of scattering events
491                                        if(index == 1 && incoherentEvent == 1)
492                                                NSingleIncoherent += 1
493                                        endif
494                                        if(index == 1 && coherentEvent == 1)
495                                                NSingleCoherent += 1
496                                        endif
497                                        if(index == 2 && coherentEvent == 1 && incoherentEvent == 0)
498                                                NDoubleCoherent += 1
499                                        endif
500                                        if(index > 1)
501                                                NMultipleScatter += 1
502                                        endif
503                                        //Print "n1,index (x,y) = ",n1,index, xpixel,ypixel
504                                       
505                                else    // if neutron escaped without interacting
506                               
507                                        // then it must be a transmitted neutron
508                                        // don't need to calculate, just increment the proper counters
509                                        MC_linear_data[xCtr][yCtr] += 1
510                                        isOn += 1
511                                        nt[0] += 1
512                                       
513                                endif           //if interacted
514                        ENDIF
515                while (!done)
516        while(n1 < imon)
517
518//      Print "Monte Carlo Done"
519        results[0] = n1
520        results[1] = n2
521        results[2] = isOn
522        results[3] = NScatterEvents             //sum of # of times that neutrons scattered
523        results[4] = NSingleCoherent            //# of events that are single, coherent
524        results[5] = NDoubleCoherent
525        results[6] = NMultipleScatter           //# of multiple scattering events
526       
527//      Print "# absorbed = ",n3
528
529//      trans_th = exp(-sig_total*thick)
530//      TRANS_exp = (N1-N2) / N1                        // Transmission
531        // dsigma/domega assuming isotropic scattering, with no absorption.
532//      Print "trans_exp = ",trans_exp
533//      Print "total # of neutrons reaching 2D detector",isOn
534//      Print "fraction of incident neutrons reaching detector ",isOn/iMon
535       
536//      Print "Total number of neutrons = ",N1
537//      Print "Total number of neutrons that interact = ",N2
538//      Print "Fraction of singly scattered neutrons = ",sum(j1,-inf,inf)/N2
539//      results[2] = N2                                         //number that scatter
540//      results[3] = isOn - MC_linear_data[xCtr][yCtr]                  //# scattered reaching detector minus zero angle
541
542       
543//      Tabs = (N1-N3)/N1
544//      Ttot = (N1-N2)/N1
545//      I1_sumI = NN[0]/(N2-N3)
546//      Print "Tabs = ",Tabs
547//      Print "Transmitted neutrons = ",Ttot
548//      results[8] = Ttot
549//      Print "I1 / all I1 = ", I1_sumI
550
551End
552////////        end of main function for calculating multiple scattering
553
554
555// returns the random deviate as a wave
556// and the total SAS cross-section [1/cm] sig_sas
557Function CalculateRandomDeviate(func,coef,lam,outWave,SASxs)
558        FUNCREF SANSModelAAO_MCproto func
559        WAVE coef
560        Variable lam
561        String outWave
562        Variable &SASxs
563
564        Variable nPts_ran=10000,qu
565        qu = 4*pi/lam           
566       
567//      Make/O/N=(nPts_ran)/D root:Packages:NIST:SAS:Gq,root:Packages:NIST:SAS:xw               // if these waves are 1000 pts, the results are "pixelated"
568//      WAVE Gq = root:Packages:NIST:SAS:gQ
569//      WAVE xw = root:Packages:NIST:SAS:xw
570
571// hard-wired into the Simulation directory rather than the SAS folder.
572// plotting resolution-smeared models won't work any other way
573        Make/O/N=(nPts_ran)/D root:Simulation:Gq,root:Simulation:xw             // if these waves are 1000 pts, the results are "pixelated"
574        WAVE Gq = root:Simulation:gQ
575        WAVE xw = root:Simulation:xw
576        SetScale/I x (0+1e-6),qu*(1-1e-10),"", Gq,xw                    //don't start at zero or run up all the way to qu to avoid numerical errors
577
578///
579/// if all of the coefficients are well-behaved, then the last point is the background
580// and I can set it to zero here (only for the calculation)
581        Duplicate/O coef,tmp_coef
582        Variable num=numpnts(coef)
583        tmp_coef[num-1] = 0
584       
585        xw=x                                                                                            //for the AAO
586        func(tmp_coef,Gq,xw)                                                                    //call as AAO
587
588//      Gq = x*Gq                                                                                                       // SAS approximation
589        Gq = Gq*sin(2*asin(x/qu))/sqrt(1-(x/qu))                        // exact
590        //
591        Integrate/METH=1 Gq/D=Gq_INT
592       
593//      SASxs = lam*lam/2/pi*Gq_INT[nPts_ran-1]                 //if the approximation is used
594        SASxs = lam*Gq_INT[nPts_ran-1]
595       
596        Gq_INT /= Gq_INT[nPts_ran-1]
597       
598        Duplicate/O Gq_INT $outWave
599
600        return(0)
601End
602
603//phi is defined from +x axis, proceeding CCW around [0,2Pi]
604Threadsafe Function FindPhi(vx,vy)
605        variable vx,vy
606       
607        variable phi
608       
609        phi = atan(vy/vx)               //returns a value from -pi/2 to pi/2
610       
611        // special cases
612        if(vx==0 && vy > 0)
613                return(pi/2)
614        endif
615        if(vx==0 && vy < 0)
616                return(3*pi/2)
617        endif
618        if(vx >= 0 && vy == 0)
619                return(0)
620        endif
621        if(vx < 0 && vy == 0)
622                return(pi)
623        endif
624       
625       
626       
627        if(vx > 0 && vy > 0)
628                return(phi)
629        endif
630        if(vx < 0 && vy > 0)
631                return(abs(phi) + pi/2)
632        endif
633        if(vx < 0 && vy < 0)
634                return(phi + pi)
635        endif
636        if( vx > 0 && vy < 0)
637                return(abs(phi) + 3*pi/2)
638        endif
639       
640        return(phi)
641end
642
643
644ThreadSafe Function FindPixel(testQ,testPhi,lam,sdd,pixSize,xCtr,yCtr,xPixel,yPixel)
645        Variable testQ,testPhi,lam,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel
646
647        Variable theta,dy,dx,qx,qy
648        //decompose to qx,qy
649        qx = testQ*cos(testPhi)
650        qy = testQ*sin(testPhi)
651
652        //convert qx,qy to pixel locations relative to # of pixels x, y from center
653        theta = 2*asin(qy*lam/4/pi)
654        dy = sdd*tan(theta)
655        yPixel = round(yCtr + dy/pixSize)
656       
657        theta = 2*asin(qx*lam/4/pi)
658        dx = sdd*tan(theta)
659        xPixel = round(xCtr + dx/pixSize)
660
661        //if on detector, return xPix and yPix values, otherwise -1
662        if(yPixel > 127 || yPixel < 0)
663                yPixel = -1
664        endif
665        if(xPixel > 127 || xPixel < 0)
666                xPixel = -1
667        endif
668       
669        return(0)
670End
671
672Function MC_CheckFunctionAndCoef(funcStr,coefStr)
673        String funcStr,coefStr
674       
675        SVAR/Z listStr=root:Packages:NIST:coefKWStr
676        if(SVAR_Exists(listStr) == 1)
677                String properCoefStr = StringByKey(funcStr, listStr  ,"=",";",0)
678                if(cmpstr("",properCoefStr)==0)
679                        return(0)               //false, no match found, so properCoefStr is returned null
680                endif
681                if(cmpstr(coefStr,properCoefStr)==0)
682                        return(1)               //true, the coef is the correct match
683                endif
684        endif
685        return(0)                       //false, wrong coef
686End
687
688Function/S MC_getFunctionCoef(funcStr)
689        String funcStr
690
691        SVAR/Z listStr=root:Packages:NIST:coefKWStr
692        String coefStr=""
693        if(SVAR_Exists(listStr) == 1)
694                coefStr = StringByKey(funcStr, listStr  ,"=",";",0)
695        endif
696        return(coefStr)
697End
698
699Function SANSModelAAO_MCproto(w,yw,xw)
700        Wave w,yw,xw
701       
702        Print "in SANSModelAAO_MCproto function"
703        return(1)
704end
705
706Function/S MC_FunctionPopupList()
707        String list,tmp
708        list = FunctionList("*",";","KIND:10")          //get every user defined curve fit function
709
710        //now start to remove everything the user doesn't need to see...
711               
712        tmp = FunctionList("*_proto",";","KIND:10")             //prototypes
713        list = RemoveFromList(tmp, list  ,";")
714        //prototypes that show up if GF is loaded
715        list = RemoveFromList("GFFitFuncTemplate", list)
716        list = RemoveFromList("GFFitAllAtOnceTemplate", list)
717        list = RemoveFromList("NewGlblFitFunc", list)
718        list = RemoveFromList("NewGlblFitFuncAllAtOnce", list)
719        list = RemoveFromList("GlobalFitFunc", list)
720        list = RemoveFromList("GlobalFitAllAtOnce", list)
721        list = RemoveFromList("GFFitAAOStructTemplate", list)
722        list = RemoveFromList("NewGF_SetXWaveInList", list)
723        list = RemoveFromList("NewGlblFitFuncAAOStruct", list)
724       
725        // more to remove as a result of 2D/Gizmo
726        list = RemoveFromList("A_WMRunLessThanDelta", list)
727        list = RemoveFromList("WMFindNaNValue", list)
728        list = RemoveFromList("WM_Make3DBarChartParametricWave", list)
729        list = RemoveFromList("UpdateQxQy2Mat", list)
730        list = RemoveFromList("MakeBSMask", list)
731       
732        // MOTOFIT/GenFit bits
733        tmp = "GEN_allatoncefitfunc;GEN_fitfunc;GetCheckBoxesState;MOTO_GFFitAllAtOnceTemplate;MOTO_GFFitFuncTemplate;MOTO_NewGF_SetXWaveInList;MOTO_NewGlblFitFunc;MOTO_NewGlblFitFuncAllAtOnce;"
734        list = RemoveFromList(tmp, list  ,";")
735
736        // SANS Reduction bits
737        tmp = "ASStandardFunction;Ann_1D_Graph;Avg_1D_Graph;BStandardFunction;CStandardFunction;Draw_Plot1D;MyMat2XYZ;NewDirection;SANSModelAAO_MCproto;Monte_SANS_Threaded;Monte_SANS_NotThreaded;Monte_SANS_W1;Monte_SANS_W2;"
738        list = RemoveFromList(tmp, list  ,";")
739        list = RemoveFromList("Monte_SANS", list)
740
741        tmp = FunctionList("f*",";","NPARAMS:2")                //point calculations
742        list = RemoveFromList(tmp, list  ,";")
743       
744        tmp = FunctionList("fSmear*",";","NPARAMS:3")           //smeared dependency calculations
745        list = RemoveFromList(tmp, list  ,";")
746       
747        //non-fit functions that I can't seem to filter out
748        list = RemoveFromList("BinaryHS_PSF11;BinaryHS_PSF12;BinaryHS_PSF22;EllipCyl_Integrand;PP_Inner;PP_Outer;Phi_EC;TaE_Inner;TaE_Outer;",list,";")
749////////////////
750
751        //more functions from analysis models (2008)
752        tmp = "Barbell_Inner;Barbell_Outer;Barbell_integrand;BCC_Integrand;Integrand_BCC_Inner;Integrand_BCC_Outer;"
753        list = RemoveFromList(tmp, list  ,";")
754        tmp = "CapCyl;CapCyl_Inner;CapCyl_Outer;ConvLens;ConvLens_Inner;ConvLens_Outer;"
755        list = RemoveFromList(tmp, list  ,";")
756        tmp = "Dumb;Dumb_Inner;Dumb_Outer;FCC_Integrand;Integrand_FCC_Inner;Integrand_FCC_Outer;"
757        list = RemoveFromList(tmp, list  ,";")
758        tmp = "Integrand_SC_Inner;Integrand_SC_Outer;SC_Integrand;SphCyl;SphCyl_Inner;SphCyl_Outer;"
759        list = RemoveFromList(tmp, list  ,";")
760
761        //simplify the display, forcing smeared calculations behind the scenes
762        tmp = FunctionList("Smear*",";","NPARAMS:1")            //smeared dependency calculations
763        list = RemoveFromList(tmp, list  ,";")
764
765        if(strlen(list)==0)
766                list = "No functions plotted"
767        endif
768       
769        list = SortList(list)
770       
771        list = "default;"+list
772        return(list)
773End             
774
775
776//Function Scat_Angle(Ran,R_DAB,wavelength)
777//      Variable Ran,r_dab,wavelength
778//
779//      Variable qq,arg,theta
780//      qq = 1. / ( R_DAB*sqrt(1.0/Ran - 1.0) )
781//      arg = qq*wavelength/(4*pi)
782//      If (arg < 1.0)
783//              theta = 2.*asin(arg)
784//      else
785//              theta = pi
786//      endif
787//      Return (theta)
788//End
789
790//calculates new direction (xyz) from an old direction
791//theta and phi don't change
792ThreadSafe Function NewDirection(vx,vy,vz,theta,phi)
793        Variable &vx,&vy,&vz
794        Variable theta,phi
795       
796        Variable err=0,vx0,vy0,vz0
797        Variable nx,ny,mag_xy,tx,ty,tz
798       
799        //store old direction vector
800        vx0 = vx
801        vy0 = vy
802        vz0 = vz
803       
804        mag_xy = sqrt(vx0*vx0 + vy0*vy0)
805        if(mag_xy < 1e-12)
806                //old vector lies along beam direction
807                nx = 0
808                ny = 1
809                tx = 1
810                ty = 0
811                tz = 0
812        else
813                Nx = -Vy0 / Mag_XY
814                Ny = Vx0 / Mag_XY
815                Tx = -Vz0*Vx0 / Mag_XY
816                Ty = -Vz0*Vy0 / Mag_XY
817                Tz = Mag_XY
818        endif
819       
820        //new scattered direction vector
821        Vx = cos(phi)*sin(theta)*Tx + sin(phi)*sin(theta)*Nx + cos(theta)*Vx0
822        Vy = cos(phi)*sin(theta)*Ty + sin(phi)*sin(theta)*Ny + cos(theta)*Vy0
823        Vz = cos(phi)*sin(theta)*Tz + cos(theta)*Vz0
824       
825        Return(err)
826End
827
828ThreadSafe Function path_len(aval,sig_tot)
829        Variable aval,sig_tot
830       
831        Variable retval
832       
833        retval = -1*ln(1-aval)/sig_tot
834       
835        return(retval)
836End
837
838// globals are initialized in SASCALC.ipf
839Window MC_SASCALC() : Panel
840        PauseUpdate; Silent 1           // building window...
841        NewPanel /W=(92,556,390,1028)/K=1 as "SANS Simulator"
842        SetVariable MC_setvar0,pos={28,73},size={144,15},bodyWidth=80,title="# of neutrons"
843        SetVariable MC_setvar0,format="%5.4g"
844        SetVariable MC_setvar0,limits={-inf,inf,100},value= root:Packages:NIST:SAS:gImon
845        SetVariable MC_setvar0_1,pos={28,119},size={131,15},bodyWidth=60,title="Thickness (cm)"
846        SetVariable MC_setvar0_1,limits={-inf,inf,0.1},value= root:Packages:NIST:SAS:gThick
847        SetVariable MC_setvar0_2,pos={28,96},size={149,15},bodyWidth=60,title="Incoherent XS (cm)"
848        SetVariable MC_setvar0_2,limits={-inf,inf,0.1},value= root:Packages:NIST:SAS:gSig_incoh
849        SetVariable MC_setvar0_3,pos={28,142},size={150,15},bodyWidth=60,title="Sample Radius (cm)"
850        SetVariable MC_setvar0_3,limits={-inf,inf,0.1},value= root:Packages:NIST:SAS:gR2
851        PopupMenu MC_popup0,pos={13,13},size={165,20},proc=MC_ModelPopMenuProc,title="Model Function"
852        PopupMenu MC_popup0,mode=1,value= #"MC_FunctionPopupList()"
853        Button MC_button0,pos={17,181},size={130,20},proc=MC_DoItButtonProc,title="Do MC Simulation"
854        Button MC_button1,pos={181,181},size={80,20},proc=MC_Display2DButtonProc,title="Show 2D"
855        SetVariable setvar0_3,pos={105,484},size={50,20},disable=1
856        GroupBox group0,pos={15,42},size={267,130},title="Monte Carlo"
857        SetVariable cntVar,pos={190,73},size={80,15},proc=CountTimeSetVarProc,title="time(s)"
858        SetVariable cntVar,format="%d"
859        SetVariable cntVar,limits={1,10,1},value= root:Packages:NIST:SAS:gCntTime
860       
861        String fldrSav0= GetDataFolder(1)
862        SetDataFolder root:Packages:NIST:SAS:
863        Edit/W=(13,217,283,450)/HOST=#  results_desc,results
864        ModifyTable format(Point)=1,width(Point)=0,width(results_desc)=150
865        SetDataFolder fldrSav0
866        RenameWindow #,T_results
867        SetActiveSubwindow ##
868EndMacro
869
870
871Function CountTimeSetVarProc(sva) : SetVariableControl
872        STRUCT WMSetVariableAction &sva
873
874        switch( sva.eventCode )
875                case 1: // mouse up
876                case 2: // Enter key
877                case 3: // Live update
878                        Variable dval = sva.dval
879
880                        // get the neutron flux, multiply, and reset the global for # neutrons
881                        NVAR imon=root:Packages:NIST:SAS:gImon
882                        imon = dval*beamIntensity()
883                       
884                        break
885        endswitch
886
887        return 0
888End
889
890
891Function MC_ModelPopMenuProc(pa) : PopupMenuControl
892        STRUCT WMPopupAction &pa
893
894        switch( pa.eventCode )
895                case 2: // mouse up
896                        Variable popNum = pa.popNum
897                        String popStr = pa.popStr
898                        SVAR gStr = root:Packages:NIST:SAS:gFuncStr
899                        gStr = popStr
900                       
901                        break
902        endswitch
903
904        return 0
905End
906
907Function MC_DoItButtonProc(ba) : ButtonControl
908        STRUCT WMButtonAction &ba
909
910        switch( ba.eventCode )
911                case 2: // mouse up
912                        // click code here
913                        NVAR doMC = root:Packages:NIST:SAS:gDoMonteCarlo
914                        doMC = 1
915                        ReCalculateInten(1)
916                        doMC = 0                //so the next time won't be MC
917                        break
918        endswitch
919
920        return 0
921End
922
923
924Function MC_Display2DButtonProc(ba) : ButtonControl
925        STRUCT WMButtonAction &ba
926
927        switch( ba.eventCode )
928                case 2: // mouse up
929                        // click code here
930                        Execute "ChangeDisplay(\"SAS\")"
931                        break
932        endswitch
933
934        return 0
935End
936
937// after a 2d data image is averaged in the usual way, take the waves and generate a "fake" folder of the 1d
938// data, to appear as if it was loaded from a real data file.
939//
940// currently only works with SANS data, but can later be expanded to generate fake USANS data sets
941//
942Function        Fake1DDataFolder(qval,aveint,sigave,sigmaQ,qbar,fSubs,dataFolder)
943        WAVE qval,aveint,sigave,sigmaQ,qbar,fSubs
944        String dataFolder
945
946        String baseStr=dataFolder
947        if(DataFolderExists("root:"+baseStr))
948                SetDataFolder $("root:"+baseStr)
949        else
950                NewDataFolder/S $("root:"+baseStr)
951        endif
952
953        ////overwrite the existing data, if it exists
954        Duplicate/O qval, $(baseStr+"_q")
955        Duplicate/O aveint, $(baseStr+"_i")
956        Duplicate/O sigave, $(baseStr+"_s")
957//      Duplicate/O sigmaQ, $(baseStr+"sq")
958//      Duplicate/O qbar, $(baseStr+"qb")
959//      Duplicate/O fSubS, $(baseStr+"fs")
960
961        // need to switch based on SANS/USANS
962        if (isSANSResolution(sigave[0]))                //checks to see if the first point of the wave is <0]
963                // make a resolution matrix for SANS data
964                Variable np=numpnts(qval)
965                Make/D/O/N=(np,4) $(baseStr+"_res")
966                Wave res=$(baseStr+"_res")
967               
968                res[][0] = sigmaQ[p]            //sigQ
969                res[][1] = qBar[p]              //qBar
970                res[][2] = fSubS[p]             //fShad
971                res[][3] = qval[p]              //Qvalues
972               
973                // keep a copy of everything in SAS too... the smearing wrapper function looks for
974                // data in folders based on waves it is passed - an I lose control of that
975                Duplicate/O res, $("root:Packages:NIST:SAS:"+baseStr+"_res")
976                Duplicate/O qval,  $("root:Packages:NIST:SAS:"+baseStr+"_q")
977                Duplicate/O aveint,  $("root:Packages:NIST:SAS:"+baseStr+"_i")
978                Duplicate/O sigave,  $("root:Packages:NIST:SAS:"+baseStr+"_s")
979        else
980                //the data is USANS data
981                // nothing done here yet
982//              dQv = -$w3[0]
983               
984//              USANS_CalcWeights(baseStr,dQv)
985               
986        endif
987
988        //clean up             
989        SetDataFolder root:
990       
991End
992
993
994
995/////UNUSED, testing routines that have not been updated to work with SASCALC
996//
997//Macro Simulate2D_MonteCarlo(imon,r1,r2,xCtr,yCtr,sdd,thick,wavelength,sig_incoh,funcStr)
998//      Variable imon=100000,r1=0.6,r2=0.8,xCtr=100,yCtr=64,sdd=400,thick=0.1,wavelength=6,sig_incoh=0.1
999//      String funcStr
1000//      Prompt funcStr, "Pick the model function", popup,       MC_FunctionPopupList()
1001//     
1002//      String coefStr = MC_getFunctionCoef(funcStr)
1003//      Variable pixSize = 0.5          // can't have 11 parameters in macro!
1004//     
1005//      if(!MC_CheckFunctionAndCoef(funcStr,coefStr))
1006//              Abort "The coefficients and function type do not match. Please correct the selections in the popup menus."
1007//      endif
1008//     
1009//      Make/O/D/N=10 root:Packages:NIST:SAS:results
1010//      Make/O/T/N=10 root:Packages:NIST:SAS:results_desc
1011//     
1012//      RunMonte(imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,funcStr,coefStr,results)
1013//     
1014//End
1015//
1016//Function RunMonte(imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,funcStr,coefStr)
1017//      Variable imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh
1018//      String funcStr,coefStr
1019//      WAVE results
1020//     
1021//      FUNCREF SANSModelAAO_MCproto func=$funcStr
1022//     
1023//      Monte_SANS(imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,sig_sas,ran_dev,linear_data,results)
1024//End
1025//
1026////// END UNUSED BLOCK
Note: See TracBrowser for help on using the repository browser.