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

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

added a button on 2D Sim input panel to link to the suggestions for incoherent cross section.

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