source: sans/XOP_Dev/MonteCarlo/MonteCarlo2.c @ 789

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

Added wavelength distribution as part of the simulation

Fixed bug where pixel locations of (-1,0) row or column were incorrectly rounded towards zero, and ended up on the detector in row 0 or column 0, effectively doubling the count there.

Separated the 4 functions into 4 separate files, for ease in checking that they are all the same. Now the only difference in each file should be the RNG used. And they had better be different - I'm almost positive now that cross-calling of the same RNG by different threads is responsible for the crashing.
File size: 15.9 KB
Line 
1/*
2 *  MonteCarlo2.c
3 *  SANSMonteCarlo
4 *
5 *  Created by Steve Kline on 7/1/10.
6 *  Copyright 2010 NCNR. All rights reserved.
7 *
8 */
9
10
11#include "XOPStandardHeaders.h"                 // Include ANSI headers, Mac headers, IgorXOP.h, XOP.h and XOPSupport.h
12#include "MonteCarlo.h"
13
14//static int gCallSpinProcess = 1;              // Set to 1 to all user abort (cmd dot) and background processing.
15
16// these versions are DIRECT COPIES of the main version in MonteCarlo.c
17// make changes there and copy them here. All that changes here is that the random
18// number calls are different.
19//
20// version X uses ran3
21// version X2 uses ran1
22// version X3 uses ran3a
23// version X4 usus ran1a
24
25int
26Monte_SANSX2(MC_ParamsPtr p) {
27        double *inputWave;                              /* pointer to double precision wave data */
28        double *ran_dev;                                /* pointer to double precision wave data */
29        double *nt;                             /* pointer to double precision wave data */
30        double *j1;                             /* pointer to double precision wave data */
31        double *j2;                             /* pointer to double precision wave data */
32        double *nn;                             /* pointer to double precision wave data */
33//      double *MC_linear_data;                         /* pointer to double precision wave data */
34        double *results;                                /* pointer to double precision wave data */
35        double retVal;                          //return value
36
37        long imon;
38        double r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,sig_sas;
39        long ind,index,n_index;
40        double qmax,theta_max,q0,zpow;
41        long n1,n2,n3;
42        double dth,zz,xx,yy,phi;
43        double theta,ran,ll,rr;
44        long done,find_theta,err;               //used as logicals
45        long xPixel,yPixel;
46        double vx,vy,vz,theta_z;
47        double sig_abs,ratio,sig_total;
48        double testQ,testPhi,left,delta,dummy,pi;
49        double sigabs_0,num_bins;
50        long NSingleIncoherent,NSingleCoherent,NScatterEvents,incoherentEvent,coherentEvent;
51        long NDoubleCoherent,NMultipleScatter,isOn,xCtr_long,yCtr_long;
52        long NMultipleCoherent,NCoherentEvents;
53        double deltaLam,v1,v2,currWavelength,rsq,fac;           //for simulating wavelength distribution
54
55
56        // for accessing the 2D wave data, direct method (see the WaveAccess example XOP)
57        waveHndl wavH;
58//      int waveType,hState;
59        long numDimensions;
60        long dimensionSizes[MAX_DIMENSIONS+1];
61//      char* dataStartPtr;
62//      long dataOffset;
63//      long numRows, numColumns;
64        long numRows_ran_dev;
65//      double *dp0, *dp;
66        double value[2];                                // Pointers used for double data.
67        long seed;
68        long indices[MAX_DIMENSIONS];
69       
70//      char buf[256];
71               
72        /* check that wave handles are all valid */
73        if (p->inputWaveH == NIL) {
74                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
75                return(NON_EXISTENT_WAVE);
76        }
77        if (p->ran_devH == NIL) {
78                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
79                return(NON_EXISTENT_WAVE);
80        }       
81        if (p->ntH == NIL) {
82                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
83                return(NON_EXISTENT_WAVE);
84        }
85        if (p->j1H == NIL) {
86                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
87                return(NON_EXISTENT_WAVE);
88        }
89        if (p->j2H == NIL) {
90                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
91                return(NON_EXISTENT_WAVE);
92        }
93        if (p->nnH == NIL) {
94                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
95                return(NON_EXISTENT_WAVE);
96        }
97        if (p->MC_linear_dataH == NIL) {
98                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
99                return(NON_EXISTENT_WAVE);
100        }
101        if (p->resultsH == NIL) {
102                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
103                return(NON_EXISTENT_WAVE);
104        }
105       
106        p->retVal = 0;
107       
108// trusting that all inputs are DOUBLE PRECISION WAVES!!!
109        inputWave = WaveData(p->inputWaveH);
110        ran_dev = WaveData(p->ran_devH);
111        nt = WaveData(p->ntH);
112        j1 = WaveData(p->j1H);
113        j2 = WaveData(p->j2H);
114        nn = WaveData(p->nnH);
115//      MC_linear_data = WaveData(p->MC_linear_dataH);
116        results = WaveData(p->resultsH);
117       
118        seed = (long)results[0];
119       
120//      sprintf(buf, "input seed = %ld\r", seed);
121//      XOPNotice(buf);
122       
123        if(seed >= 0) {
124                seed = -1234509876;
125        }
126
127        dummy = ran1(&seed);            //initialize the random sequence by passing in a negative value
128        seed = 12348765;                //non-negative after that does nothing
129
130        imon = (int)inputWave[0];
131        r1 = inputWave[1];
132        r2 = inputWave[2];
133        xCtr = inputWave[3];
134        yCtr = inputWave[4];
135        sdd = inputWave[5];
136        pixSize = inputWave[6];
137        thick = inputWave[7];
138        wavelength = inputWave[8];
139        sig_incoh = inputWave[9];
140        sig_sas = inputWave[10];
141        deltaLam = inputWave[11];
142        xCtr_long = (long)(xCtr+0.5);
143        yCtr_long = (long)(yCtr+0.5);
144       
145        dummy = MDGetWaveScaling(p->ran_devH, 0, &delta, &left);                //0 is the rows
146        if (retVal = MDGetWaveDimensions(p->ran_devH, &numDimensions, dimensionSizes))
147                return retVal;
148        numRows_ran_dev = dimensionSizes[0];
149       
150        pi = 4.0*atan(1.0);     
151       
152        // access the 2D wave data for writing using the direct method
153        wavH = p->MC_linear_dataH;
154        if (wavH == NIL)
155                return NOWAV;
156
157//      waveType = WaveType(wavH);
158//      if (waveType & NT_CMPLX)
159//              return NO_COMPLEX_WAVE;
160//      if (waveType==TEXT_WAVE_TYPE)
161//              return NUMERIC_ACCESS_ON_TEXT_WAVE;
162//      if (retVal = MDGetWaveDimensions(wavH, &numDimensions, dimensionSizes))
163//              return retVal;
164//      numRows = dimensionSizes[0];
165//      numColumns = dimensionSizes[1];
166       
167//      if (retVal = MDAccessNumericWaveData(wavH, kMDWaveAccessMode0, &dataOffset))
168//              return retVal;
169               
170//      hState = MoveLockHandle(wavH);          // So wave data can't move. Remember to call HSetState when done.
171//      dataStartPtr = (char*)(*wavH) + dataOffset;
172//      dp0 = (double*)dataStartPtr;                    // Pointer to the start of the 2D wave data.
173       
174//scattering power and maximum qvalue to bin
175//      zpow = .1               //scattering power, calculated below
176        qmax = 4.0*pi/wavelength;       //maximum Q to bin 1D data. (A-1) (not really used)
177        sigabs_0 = 0.0;         // ignore absorption cross section/wavelength [1/(cm A)]
178        n_index = 50;   // maximum number of scattering events per neutron
179        num_bins = 200;         //number of 1-D bins (not really used)
180       
181//c       total SAS cross-section
182//
183        zpow = sig_sas*thick;                   //since I now calculate the sig_sas from the model
184        sig_abs = sigabs_0 * wavelength;
185        sig_total = sig_abs + sig_sas + sig_incoh;
186//      Print "The TOTAL XSECTION. (CM-1) is ",sig_total
187//      Print "The TOTAL SAS XSECTION. (CM-1) is ",sig_sas
188//      results[0] = sig_total;
189//      results[1] = sig_sas;
190//      RATIO = SIG_ABS / SIG_TOTAL
191        ratio = sig_incoh / sig_total;
192       
193        theta_max = wavelength*qmax/(2*pi);
194//C     SET Theta-STEP SIZE.
195        dth = theta_max/num_bins;
196//      Print "theta bin size = dth = ",dth
197
198//C     INITIALIZE COUNTERS.
199        n1 = 0;
200        n2 = 0;
201        n3 = 0;
202        NSingleIncoherent = 0;
203        NSingleCoherent = 0;
204        NDoubleCoherent = 0;
205        NMultipleScatter = 0;
206        NScatterEvents = 0;
207        NMultipleCoherent = 0;
208        NCoherentEvents = 0;
209       
210        isOn = 0;
211       
212//C     MONITOR LOOP - looping over the number of incedent neutrons
213//note that zz, is the z-position in the sample - NOT the scattering power
214// NOW, start the loop, throwing neutrons at the sample.
215        do {
216                ////SpinProcess() IS A CALLBACK, and not good for Threading!
217//              if ((n1 % 1000 == 0) && gCallSpinProcess && SpinProcess()) {            // Spins cursor and allows background processing.
218//                              retVal = -1;                                                            // User aborted.
219//                              break;
220//              }
221       
222                vx = 0.0;                       // Initialize direction vector.
223                vy = 0.0;
224                vz = 1.0;
225               
226                theta = 0.0;            //      Initialize scattering angle.
227                phi = 0.0;                      //      Intialize azimuthal angle.
228                n1 += 1;                        //      Increment total number neutrons counter.
229                done = 0;                       //      True when neutron is absorbed or when  scattered out of the sample.
230                index = 0;                      //      Set counter for number of scattering events.
231                zz = 0.0;                       //      Set entering dimension of sample.
232                incoherentEvent = 0;
233                coherentEvent = 0;
234               
235                do      {                               //      Makes sure position is within circle.
236                        ran = ran1(&seed);              //[0,1]
237                        xx = 2.0*r1*(ran-0.5);          //X beam position of neutron entering sample.
238                        ran = ran1(&seed);              //[0,1]
239                        yy = 2.0*r1*(ran-0.5);          //Y beam position ...
240                        rr = sqrt(xx*xx+yy*yy);         //Radial position of neutron in incident beam.
241                } while(rr>r1);
242
243                //pick the wavelength out of the wavelength spread, approximate as a gaussian
244                // from NR - pg 288. Needs random # from [0,1]. del is deltaLam/lam (as FWHM) and the
245                // 2.35 converts to a gaussian std dev.
246                do {
247                        v1=2.0*ran1(&seed)-1.0;
248                        v2=2.0*ran1(&seed)-1.0;
249                        rsq=v1*v1+v2*v2;
250                } while (rsq >= 1.0 || rsq == 0.0);
251                fac=sqrt(-2.0*log(rsq)/rsq);
252               
253                //              gset=v1*fac             //technically, I'm throwing away one of the two values
254                currWavelength = (v2*fac)*deltaLam*wavelength/2.35 + wavelength;
255               
256                do {   //Scattering Loop, will exit when "done" == 1
257                                // keep scattering multiple times until the neutron exits the sample
258                        ran = ran1(&seed);              //[0,1]  RANDOM NUMBER FOR DETERMINING PATH LENGTH
259                        ll = path_len(ran,sig_total);
260                        //Determine new scattering direction vector.
261                        err = NewDirection(&vx,&vy,&vz,theta,phi);              //vx,vy,vz updated, theta, phi unchanged by function
262                                                                       
263                        //X,Y,Z-POSITION OF SCATTERING EVENT.
264                        xx += ll*vx;
265                        yy += ll*vy;
266                        zz += ll*vz;
267                        rr = sqrt(xx*xx+yy*yy);         //radial position of scattering event.
268
269                        //sprintf(buf, "xx,yy,zz,vx,vy,vz,ll = %g %g %g %g %g %g %g\r",xx,yy,zz,vx,vy,vz,ll);
270                        //XOPNotice(buf);
271                                               
272                        //Check whether interaction occurred within sample volume.
273                        if (((zz > 0.0) && (zz < thick)) && (rr < r2)) {
274                                //NEUTRON INTERACTED.
275                                //sprintf(buf,"neutron interacted\r");
276                                //XOPNotice(buf);
277                               
278                                index += 1;                     //Increment counter of scattering events.
279                                if (index == 1) {
280                                        n2 += 1;                //Increment # of scat. neutrons
281                                }
282                                ran = ran1(&seed);              //[0,1]
283                                //Split neutron interactions into scattering and absorption events
284                                if (ran > ratio ) {             //C             NEUTRON SCATTERED coherently
285                                        //sprintf(buf,"neutron scatters coherently\r");
286                                        //XOPNotice(buf);
287                                        coherentEvent += 1;
288                                        find_theta = 0;                 //false
289                                        do {
290                                                // pick a q-value from the deviate function
291                                                // pnt2x truncates the point to an integer before returning the x
292                                                // so get it from the wave scaling instead
293//                                              q0 =left + binarysearchinterp(ran_dev,ran1(seed))*delta;
294                                               
295                                                q0 =left + locate_interp(ran_dev,numRows_ran_dev,ran1(&seed))*delta;
296                                                theta = q0/2/pi*currWavelength;         //SAS approximation
297                                               
298                                                find_theta = 1;         //always accept
299
300                                                //sprintf(buf, "after locate_interp call q0 = %g, theta = %g,left = %g,delta = %g\r",q0,theta,left,delta);
301                                                //XOPNotice(buf);
302
303                                        } while(!find_theta);
304                                       
305                                        ran = ran1(&seed);              //[0,1]
306                                        phi = 2.0*pi*ran;                       //Chooses azimuthal scattering angle.
307                                } else {
308                                        //NEUTRON scattered incoherently
309                                        //sprintf(buf,"neutron scatters incoherent\r");
310                                        //XOPNotice(buf);
311                                        incoherentEvent += 1;
312                                  // phi and theta are random over the entire sphere of scattering
313                                        // !can't just choose random theta and phi, won't be random over sphere solid angle
314
315                                        ran = ran1(&seed);              //[0,1]
316                                        theta = acos(2.0*ran-1);
317                       
318                                        ran = ran1(&seed);              //[0,1]
319                                        phi = 2.0*pi*ran;                       //Chooses azimuthal scattering angle.
320                                }               //(ran > ratio)
321                        } else {
322                                //NEUTRON ESCAPES FROM SAMPLE -- bin it somewhere                                                               
323                                done = 1;               //done = true, will exit from loop
324                                //Increment #scattering events array
325                                MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
326                                indices[0] =index;                      //this sets access to nn[index]
327                                if (index <= n_index) {
328                                        if (retVal = MDGetNumericWavePointValue(p->nnH, indices, value))
329                                                return retVal;
330                                        value[0] += 1; // add one to the value
331                                        if (retVal = MDSetNumericWavePointValue(p->nnH, indices, value))
332                                                return retVal;
333                                //      nn[index] += 1;
334                                }
335                                                                                               
336                                if( index != 0) {               //neutron was scattered, figure out where it went
337                                        theta_z = acos(vz);             // Angle (= 2theta) WITH respect to z axis.
338                                        testQ = 2*pi*sin(theta_z)/currWavelength;
339                                       
340                                        // pick a random phi angle, and see if it lands on the detector
341                                        // since the scattering is isotropic, I can safely pick a new, random value
342                                        // this would not be true if simulating anisotropic scattering.
343                                        testPhi = ran1(&seed)*2*pi;
344                                       
345                                        // is it on the detector?       
346                                        FindPixel(testQ,testPhi,currWavelength,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel);
347                                                                                                       
348                                        if(xPixel != -1 && yPixel != -1) {
349                                                isOn += 1;
350                                                MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
351                                                indices[0] = xPixel;
352                                                indices[1] = yPixel;
353                                                if (retVal = MDGetNumericWavePointValue(wavH, indices, value))
354                                                        return retVal;
355                                                value[0] += 1; // Real part
356                                                if (retVal = MDSetNumericWavePointValue(wavH, indices, value))
357                                                        return retVal;
358                                                //if(index==1)  // only the single scattering events
359                                                        //dp = dp0 + xPixel + yPixel*numColumns;                //offset the pointer to the exact memory location
360                                                        //*dp += 1;             //increment the value there
361                                                //endif
362                                        }
363                                       
364
365        /*              is this causing me a problem since I'm not locking these? Probably not, since it crashes even if I comment these out... */
366                                        if(theta_z < theta_max) {
367                                                //Choose index for scattering angle array.
368                                                //IND = NINT(THETA_z/DTH + 0.4999999)
369                                                ind = (long)(theta_z/dth + 0.4999999);          //round is eqivalent to nint()
370                                                nt[ind] += 1;                   //Increment bin for angle.
371                                                //Increment angle array for single scattering events.
372                                                if (index == 1) {
373                                                        j1[ind] += 1;
374                                                }
375                                                //Increment angle array for double scattering events.
376                                                if (index == 2) {
377                                                        j2[ind] += 1;
378                                                }
379                                        }
380        /**/
381                                       
382                                        // increment all of the counters now since done==1 here and I'm sure to exit and get another neutron
383                                        NScatterEvents += index;                //total number of scattering events
384                                        if(index == 1 && incoherentEvent == 1) {
385                                                NSingleIncoherent += 1;
386                                        }
387                                        if(index == 1 && coherentEvent == 1) {
388                                                NSingleCoherent += 1;
389                                        }
390                                        if(index == 2 && coherentEvent == 1 && incoherentEvent == 0) {
391                                                NDoubleCoherent += 1;
392                                        }
393                                        if(index > 1) {
394                                                NMultipleScatter += 1;
395                                        }
396                                        if(coherentEvent >= 1 && incoherentEvent == 0) {
397                                                NCoherentEvents += 1;
398                                        }
399                                        if(coherentEvent > 1 && incoherentEvent == 0) {
400                                                NMultipleCoherent += 1;
401                                        }
402
403                                } else {        // index was zero, neutron must be transmitted, so just increment the proper counters and data
404                                                isOn += 1;
405                                                nt[0] += 1;
406                                                MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
407                                                //indices[0] = xCtr_long;               //don't put everything in one pixel
408                                                //indices[1] = yCtr_long;
409                                                indices[0] = (long)(xCtr+xx/pixSize+0.5);
410                                                indices[1] = (long)(yCtr+yy/pixSize+0.5);
411                                                // check for valid indices - got an XOP error, probably from here
412                                                if(indices[0] > 127) indices[0] = 127;
413                                                if(indices[0] < 0) indices[0] = 0;
414                                                if(indices[1] > 127) indices[1] = 127;
415                                                if(indices[1] < 0) indices[1] = 0;
416                                               
417                                                if (retVal = MDGetNumericWavePointValue(wavH, indices, value))
418                                                        return retVal;
419                                                value[0] += 1; // Real part
420                                                if (retVal = MDSetNumericWavePointValue(wavH, indices, value))
421                                                        return retVal;
422                                        }       
423                        }
424                 } while (!done);
425        } while(n1 < imon);
426       
427// assign the results to the wave
428
429        MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
430        value[0] = (double)n1;
431        indices[0] = 0;
432        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
433                return retVal;
434        value[0] = (double)n2;
435        indices[0] = 1;
436        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
437                return retVal;
438        value[0] = (double)isOn;
439        indices[0] = 2;
440        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
441                return retVal;
442        value[0] = (double)NScatterEvents;
443        indices[0] = 3;
444        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
445                return retVal;
446        value[0] = (double)NSingleCoherent;
447        indices[0] = 4;
448        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
449                return retVal;
450        value[0] = (double)NMultipleCoherent;
451        indices[0] = 5;
452        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
453                return retVal;
454        value[0] = (double)NMultipleScatter;
455        indices[0] = 6;
456        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
457                return retVal; 
458        value[0] = (double)NCoherentEvents;
459        indices[0] = 7;
460        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
461                return retVal;
462
463//      HSetState((Handle)wavH, hState);                //release the handle of the 2D data wave
464//      WaveHandleModified(wavH);                       // Inform Igor that we have changed the wave. (CALLBACK! needed, but not allowed in Threading)
465       
466        return(0);
467}
468
Note: See TracBrowser for help on using the repository browser.