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

Last change on this file since 890 was 890, checked in by srkline, 10 years ago

Added a procedure file to allow "scripting" of the MonteCarlo? simulation: MC_SimulationScripting.ipf. See the beginning of the file for instructions, and then the examples for numbered instructions of what needs to be edited to write your own script. Not casual-user-friendly, but makes it easire for more advanced users to simulate an entire experiment.

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