source: sans/XOP_Dev/MonteCarlo/MonteCarlo3.c @ 711

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