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

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

Added a "signature" to the simulated 2D data set, setting a few pixels in the lower-left corner to known values. These points are behind the default mask.

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