Changeset 437 for sans/XOP_Dev/MonteCarlo
 Timestamp:
 Nov 3, 2008 5:33:38 PM (14 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

sans/XOP_Dev/MonteCarlo/MonteCarlo.c
r435 r437 63 63 double sigabs_0,num_bins; 64 64 long NSingleIncoherent,NSingleCoherent,NScatterEvents,incoherentEvent,coherentEvent; 65 long NDoubleCoherent,NMultipleScatter,isOn ;65 long NDoubleCoherent,NMultipleScatter,isOn,xCtr_long,yCtr_long; 66 66 67 67 // for accessing the 2D wave data, direct method (see the WaveAccess example XOP) … … 148 148 sig_incoh = inputWave[9]; 149 149 sig_sas = inputWave[10]; 150 150 xCtr_long = round(xCtr); 151 yCtr_long = round(yCtr); 152 151 153 dummy = MDGetWaveScaling(p>ran_devH, 0, &delta, &left); //0 is the rows 152 154 if (result = MDGetWaveDimensions(p>ran_devH, &numDimensions, dimensionSizes)) … … 324 326 } 325 327 326 //IF (VZ > 1.0) // FIX INVALID ARGUMENT 327 //VZ = 1.0  1.2e7 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 366 354 } 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 } 370 371 } 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 } 388 400 } 389 401 } while (!done);
Note: See TracChangeset
for help on using the changeset viewer.