Changeset 437


Ignore:
Timestamp:
Nov 3, 2008 5:33:38 PM (14 years ago)
Author:
srkline
Message:

Change to calculation to avoid unnecessary calculations for transmitted neutrons. Marginal speedup. Done to be consistent with non-XOP version.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • sans/XOP_Dev/MonteCarlo/MonteCarlo.c

    r435 r437  
    6363        double sigabs_0,num_bins; 
    6464        long NSingleIncoherent,NSingleCoherent,NScatterEvents,incoherentEvent,coherentEvent; 
    65         long NDoubleCoherent,NMultipleScatter,isOn; 
     65        long NDoubleCoherent,NMultipleScatter,isOn,xCtr_long,yCtr_long; 
    6666 
    6767        // for accessing the 2D wave data, direct method (see the WaveAccess example XOP) 
     
    148148        sig_incoh = inputWave[9]; 
    149149        sig_sas = inputWave[10]; 
    150                  
     150        xCtr_long = round(xCtr); 
     151        yCtr_long = round(yCtr); 
     152         
    151153        dummy = MDGetWaveScaling(p->ran_devH, 0, &delta, &left);                //0 is the rows 
    152154        if (result = MDGetWaveDimensions(p->ran_devH, &numDimensions, dimensionSizes)) 
     
    324326                                } 
    325327                                                                                                 
    326                                 //IF (VZ > 1.0)         // FIX INVALID ARGUMENT 
    327                                         //VZ = 1.0 - 1.2e-7 
    328                                 //ENDIF 
    329                                 theta_z = acos(vz);             // Angle (= 2theta) WITH respect to z axis. 
    330                                 testQ = 2*pi*sin(theta_z)/wavelength; 
    331                                  
    332                                 // pick a random phi angle, and see if it lands on the detector 
    333                                 // since the scattering is isotropic, I can safely pick a new, random value 
    334                                 // this would not be true if simulating anisotropic scattering. 
    335                                 testPhi = ran3(&seed)*2*pi; 
    336                                  
    337                                 // is it on the detector?        
    338                                 FindPixel(testQ,testPhi,wavelength,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel); 
    339                                                                                                  
    340                                 if(xPixel != -1 && yPixel != -1) { 
    341                                         isOn += 1; 
    342                                         MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 
    343                                         indices[0] = xPixel; 
    344                                         indices[1] = yPixel; 
    345                                         if (result = MDGetNumericWavePointValue(wavH, indices, value)) 
    346                                                 return result; 
    347                                         value[0] += 1; // Real part 
    348                                         if (result = MDSetNumericWavePointValue(wavH, indices, value)) 
    349                                                 return result; 
    350                                         //if(index==1)  // only the single scattering events 
    351                                                 //dp = dp0 + xPixel + yPixel*numColumns;                //offset the pointer to the exact memory location 
    352                                                 //*dp += 1;             //increment the value there 
    353                                         //endif 
    354                                 } 
    355                                  
    356  
    357 /*              is this causing me a problem since I'm not locking these? Probably not, since it crashes even if I comment these out... */ 
    358                                 if(theta_z < theta_max) { 
    359                                         //Choose index for scattering angle array. 
    360                                         //IND = NINT(THETA_z/DTH + 0.4999999) 
    361                                         ind = round(theta_z/dth + 0.4999999);           //round is eqivalent to nint() 
    362                                         nt[ind] += 1;                   //Increment bin for angle. 
    363                                         //Increment angle array for single scattering events. 
    364                                         if (index == 1) { 
    365                                                 j1[ind] += 1; 
     328                                if( index != 0) {               //neutron was scattered, figure out where it went 
     329                                        theta_z = acos(vz);             // Angle (= 2theta) WITH respect to z axis. 
     330                                        testQ = 2*pi*sin(theta_z)/wavelength; 
     331                                         
     332                                        // pick a random phi angle, and see if it lands on the detector 
     333                                        // since the scattering is isotropic, I can safely pick a new, random value 
     334                                        // this would not be true if simulating anisotropic scattering. 
     335                                        testPhi = ran3(&seed)*2*pi; 
     336                                         
     337                                        // is it on the detector?        
     338                                        FindPixel(testQ,testPhi,wavelength,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel); 
     339                                                                                                         
     340                                        if(xPixel != -1 && yPixel != -1) { 
     341                                                isOn += 1; 
     342                                                MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 
     343                                                indices[0] = xPixel; 
     344                                                indices[1] = yPixel; 
     345                                                if (result = MDGetNumericWavePointValue(wavH, indices, value)) 
     346                                                        return result; 
     347                                                value[0] += 1; // Real part 
     348                                                if (result = MDSetNumericWavePointValue(wavH, indices, value)) 
     349                                                        return result; 
     350                                                //if(index==1)  // only the single scattering events 
     351                                                        //dp = dp0 + xPixel + yPixel*numColumns;                //offset the pointer to the exact memory location 
     352                                                        //*dp += 1;             //increment the value there 
     353                                                //endif 
    366354                                        } 
    367                                         //Increment angle array for double scattering events. 
    368                                         if (index == 2) { 
    369                                                 j2[ind] += 1; 
     355                                         
     356 
     357        /*              is this causing me a problem since I'm not locking these? Probably not, since it crashes even if I comment these out... */ 
     358                                        if(theta_z < theta_max) { 
     359                                                //Choose index for scattering angle array. 
     360                                                //IND = NINT(THETA_z/DTH + 0.4999999) 
     361                                                ind = round(theta_z/dth + 0.4999999);           //round is eqivalent to nint() 
     362                                                nt[ind] += 1;                   //Increment bin for angle. 
     363                                                //Increment angle array for single scattering events. 
     364                                                if (index == 1) { 
     365                                                        j1[ind] += 1; 
     366                                                } 
     367                                                //Increment angle array for double scattering events. 
     368                                                if (index == 2) { 
     369                                                        j2[ind] += 1; 
     370                                                } 
    370371                                        } 
    371                                 } 
    372 /**/ 
    373                                  
    374                                 // increment all of the counters now since done==1 here and I'm sure to exit and get another neutron 
    375                                 NScatterEvents += index;                //total number of scattering events 
    376                                 if(index == 1 && incoherentEvent == 1) { 
    377                                         NSingleIncoherent += 1; 
    378                                 } 
    379                                 if(index == 1 && coherentEvent == 1) { 
    380                                         NSingleCoherent += 1; 
    381                                 } 
    382                                 if(index == 2 && coherentEvent == 1 && incoherentEvent == 0) { 
    383                                         NDoubleCoherent += 1; 
    384                                 } 
    385                                 if(index > 1) { 
    386                                         NMultipleScatter += 1; 
    387                                 } 
     372        /**/ 
     373                                         
     374                                        // increment all of the counters now since done==1 here and I'm sure to exit and get another neutron 
     375                                        NScatterEvents += index;                //total number of scattering events 
     376                                        if(index == 1 && incoherentEvent == 1) { 
     377                                                NSingleIncoherent += 1; 
     378                                        } 
     379                                        if(index == 1 && coherentEvent == 1) { 
     380                                                NSingleCoherent += 1; 
     381                                        } 
     382                                        if(index == 2 && coherentEvent == 1 && incoherentEvent == 0) { 
     383                                                NDoubleCoherent += 1; 
     384                                        } 
     385                                        if(index > 1) { 
     386                                                NMultipleScatter += 1; 
     387                                        } 
     388                                } else {        // index was zero, neutron must be transmitted, so just increment the proper counters and data 
     389                                                isOn += 1; 
     390                                                nt[0] += 1; 
     391                                                MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 
     392                                                indices[0] = xCtr_long; 
     393                                                indices[1] = yCtr_long; 
     394                                                if (result = MDGetNumericWavePointValue(wavH, indices, value)) 
     395                                                        return result; 
     396                                                value[0] += 1; // Real part 
     397                                                if (result = MDSetNumericWavePointValue(wavH, indices, value)) 
     398                                                        return result; 
     399                                        }        
    388400                        } 
    389401                 } while (!done); 
Note: See TracChangeset for help on using the changeset viewer.