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

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

Changes to how Patch interprets the match string. 3 cases now:
(1) - interpret as a run number list, just like all other run number lists and ranges.
(2) grep for text in the raw VAX binary. can be slow. did not allow generic regular expressions to be entered since of the three people on earth that know how to use them, none are neutron scatterers.
(3) filter on SDD, just like in the MRED beta ops. Radio buttons control how it's interpreted. Removed the redundant "Show Header" button, as the header is displayed for the currently popped item without user intervention.

This closes Cedric's ticket #212

Plus a few tweaks of MC and Wrapper.

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