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

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

Bug fixes in the Monte Carlo routines for SASCALC, and some improvements for functionality, including the calculation on smeared model functions. For various reasons, the MC simulation must use the unsmeared model.

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