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

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

minor changes and comments

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;"
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,60,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.