- Timestamp:
- Nov 25, 2008 4:35:13 PM (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/MultScatter_MonteCarlo_2D.ipf
r455 r456 10 10 // 11 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 be fore 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? 12 16 // - Most importantly, this needs to be checked for correctness of the MC simulation 13 17 // X how can I get the "data" on absolute scale? This would be a great comparison vs. the ideal model calculation 14 18 // X why does my integrated tau not match up with John's analytical calculations? where are the assumptions? 15 19 // - get rid of all small angle assumptions - to make sure that the calculation is correct at all angles 16 20 // - my simulated transmission is larger than what is measured, even after correcting for the quartz cell. 21 // Why? Do I need to include absorption? Just inherent problems with incoherent cross sections? 22 // 17 23 // X at the larger angles, is the "flat" detector being properly accounted for - in terms of 18 24 // the solid angle and how many counts fall in that pixel. Am I implicitly defining a spherical detector … … 268 274 Variable isOn=0,testQ,testPhi,xPixel,yPixel 269 275 Variable NSingleIncoherent,NSingleCoherent,NScatterEvents,incoherentEvent,coherentEvent 270 Variable NDoubleCoherent,NMultipleScatter 276 Variable NDoubleCoherent,NMultipleScatter,countIt,detEfficiency 277 278 detEfficiency = 1.0 //70% counting efficiency = 0.7 271 279 272 280 imon = inputWave[0] … … 291 299 num_bins = 200 //number of 1-D bins (not really used) 292 300 293 // my additions - calculate the random edeviate function as needed294 // and calculate the scattering power from the model function 301 // my additions - calculate the random deviate function as needed 302 // and calculate the scattering power from the model function (passed in as a wave) 295 303 // 296 304 Variable left = leftx(ran_dev) … … 301 309 zpow = sig_sas*thick //since I now calculate the sig_sas from the model 302 310 SIG_ABS = SIGABS_0 * WAVElength 311 sig_abs = 0.0 //cm-1 303 312 SIG_TOTAL =SIG_ABS + SIG_SAS + sig_incoh 304 313 // Print "The TOTAL XSECTION. (CM-1) is ",sig_total … … 306 315 // results[0] = sig_total //assign these after everything's done 307 316 // results[1] = sig_sas 308 // RATIO = SIG_ABS / SIG_TOTAL 317 // variable ratio1,ratio2 318 // ratio1 = sig_abs/sig_total 319 // ratio2 = sig_incoh/sig_total 320 // // 0->ratio1 = abs 321 // // ratio1 -> ratio2 = incoh 322 // // > ratio2 = coherent 309 323 RATIO = sig_incoh / SIG_TOTAL 310 324 … … 373 387 IF (((zz > 0.0) && (zz < THICK)) && (rr < r2)) 374 388 //NEUTRON INTERACTED. 389 ran = abs(enoise(1)) //[0,1] 390 391 // if(ran<ratio1) 392 // //absorption event 393 // n3 +=1 394 // done=1 395 // else 396 375 397 INDEX += 1 //Increment counter of scattering events. 376 398 IF(INDEX == 1) 377 399 N2 += 1 //Increment # of scat. neutrons 378 400 Endif 379 ran = abs(enoise(1)) //[0,1]380 401 //Split neutron interactions into scattering and absorption events 381 IF(ran > ratio ) //C NEUTRON SCATTERED coherently 402 // IF(ran > (ratio1+ratio2) ) //C NEUTRON SCATTERED coherently 403 IF(ran > ratio) //C NEUTRON SCATTERED coherently 382 404 coherentEvent = 1 383 405 FIND_THETA = 0 //false … … 414 436 PHI = 2.0*PI*Ran //Chooses azimuthal scattering angle. 415 437 ENDIF //(ran > ratio) 438 // endif // event was absorption 416 439 ELSE 417 440 //NEUTRON ESCAPES FROM SAMPLE -- bin it somewhere 418 441 DONE = 1 //done = true, will exit from loop 442 443 // countIt = 1 444 // if(abs(enoise(1)) > detEfficiency) //efficiency of 70% wired @top 445 // countIt = 0 //detector does not register 446 // endif 447 419 448 //Increment #scattering events array 420 449 If (index <= N_Index) … … 430 459 // since the scattering is isotropic, I can safely pick a new, random value 431 460 // this would not be true if simulating anisotropic scattering. 432 testPhi = abs(enoise(1))*2*Pi 461 //testPhi = abs(enoise(1))*2*Pi 462 testPhi = FindPhi(Vx,Vy) //use the exiting phi value as defined by Vx and Vy 463 433 464 // is it on the detector? 434 465 FindPixel(testQ,testPhi,wavelength,sdd,pixSize,xCtr,yCtr,xPixel,yPixel) … … 471 502 endif 472 503 //Print "n1,index (x,y) = ",n1,index, xpixel,ypixel 504 473 505 else // if neutron escaped without interacting 506 474 507 // then it must be a transmitted neutron 475 508 // don't need to calculate, just increment the proper counters … … 478 511 nt[0] += 1 479 512 480 // (Don't do this - it's a waste of time) 481 // re-initialize and go back and count this neutron again 482 // go back to xy position 483 // xx -= ll*vx 484 // yy -= ll*vy 485 // zz -= ll*vz 486 // RR = sqrt(xx*xx+yy*yy) //radial position of scattering event. 487 // 488 // Vx = 0.0 // Initialize incident direction vector. 489 // Vy = 0.0 490 // Vz = 1.0 491 // zz = 0 492 // theta = 0 493 // phi = 0 494 endif 495 513 endif //if interacted 496 514 ENDIF 497 515 while (!done) … … 507 525 results[6] = NMultipleScatter //# of multiple scattering events 508 526 509 527 // Print "# absorbed = ",n3 528 510 529 // trans_th = exp(-sig_total*thick) 511 530 // TRANS_exp = (N1-N2) / N1 // Transmission … … 581 600 return(0) 582 601 End 602 603 //phi is defined from +x axis, proceeding CCW around [0,2Pi] 604 Threadsafe Function FindPhi(vx,vy) 605 variable vx,vy 606 607 variable phi 608 609 phi = atan(vy/vx) //returns a value from -pi/2 to pi/2 610 611 // special cases 612 if(vx==0 && vy > 0) 613 return(pi/2) 614 endif 615 if(vx==0 && vy < 0) 616 return(3*pi/2) 617 endif 618 if(vx >= 0 && vy == 0) 619 return(0) 620 endif 621 if(vx < 0 && vy == 0) 622 return(pi) 623 endif 624 625 626 627 if(vx > 0 && vy > 0) 628 return(phi) 629 endif 630 if(vx < 0 && vy > 0) 631 return(abs(phi) + pi/2) 632 endif 633 if(vx < 0 && vy < 0) 634 return(phi + pi) 635 endif 636 if( vx > 0 && vy < 0) 637 return(abs(phi) + 3*pi/2) 638 endif 639 640 return(phi) 641 end 642 583 643 584 644 ThreadSafe Function FindPixel(testQ,testPhi,lam,sdd,pixSize,xCtr,yCtr,xPixel,yPixel)
Note: See TracChangeset
for help on using the changeset viewer.
