 Timestamp:
 Oct 17, 2008 3:12:18 PM (14 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/MultScatter_MonteCarlo_2D.ipf
r429 r430 7 7 // This code simulates the scattering for a selected model based on the instrument configuration 8 8 // This code is based directly on John Barker's code, which includes multiple scattering effects. 9 // A lot of the setup, wave creation, and postcalculations are done in SASCALC>ReCalculateInten() 9 10 // 10 11 // 11 12 //  Most importantly, this needs to be checked for correctness of the MC simulation 12 // how can I get the "data" on absolute scale? This would be a great comparison vs. the ideal model calculation13 // X how can I get the "data" on absolute scale? This would be a great comparison vs. the ideal model calculation 13 14 //  why does my integrated tau not match up with John's analytical calculations? where are the assumptions? 14 15 //  get rid of all small angle assumptions  to make sure that the calculation is correct at all angles 15 16 //  what is magical about Qu? Is this an assumpution? 16 //  why am I off a magical factor of 2 in my calculation of kappa to get to absolute scale (in ReCalculateInten())17 17 18 18 //  at the larger angles, is the "flat" detector being properly accounted for  in terms of 19 19 // the solid angle and how many counts fall in that pixel. Am I implicitly defining a spherical detector 20 20 // so that what I see is already "corrected"? 21 //  the MC will, of course benefit greatly from being XOPized. Maybe think about parallel implementation 22 // by allowing the data arrays to accumulate. 21 // X the MC will, of course benefit greatly from being XOPized. Maybe think about parallel implementation 22 // by allowing the data arrays to accumulate. First pass at the XOP is done. Not pretty, not the speediest (5.8x) 23 // but it is functional. Add spinCursor() for long calculations. See WaveAccess XOP example. 23 24 //  the background parameter for the model MUST be zero, or the integration for scattering 24 25 // power will be incorrect. 25 26 //  fully use the SASCALC input, most importantly, flux on sample. 26 // if no MC desired, still use the selected model27 // better display of MC results on panel27 // X if no MC desired, still use the selected model 28 // X better display of MC results on panel 28 29 //  settings for "count for X seconds" or "how long to 1E6 cts on detector" (run short sim, then multiply) 29 30 //  add quartz window scattering to the simulation somehow 30 31 //  do smeared models make any sense?? 31 // 32 33 Macro Simulate2D_MonteCarlo(imon,r1,r2,xCtr,yCtr,sdd,thick,wavelength,sig_incoh,funcStr) 34 Variable imon=100000,r1=0.6,r2=0.8,xCtr=100,yCtr=64,sdd=400,thick=0.1,wavelength=6,sig_incoh=0.1 35 String funcStr 36 Prompt funcStr, "Pick the model function", popup, MC_FunctionPopupList() 37 38 String coefStr = MC_getFunctionCoef(funcStr) 39 Variable pixSize = 0.5 // can't have 11 parameters in macro! 40 41 if(!MC_CheckFunctionAndCoef(funcStr,coefStr)) 42 Abort "The coefficients and function type do not match. Please correct the selections in the popup menus." 43 endif 44 45 Make/O/D/N=10 root:Packages:NIST:SAS:results 46 Make/O/T/N=10 root:Packages:NIST:SAS:results_desc 47 48 RunMonte(imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,funcStr,coefStr,results) 49 50 End 51 52 Function RunMonte(imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,funcStr,coefStr) 53 Variable imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh 54 String funcStr,coefStr 55 WAVE results 56 57 FUNCREF SANSModelAAO_MCproto func=$funcStr 58 59 Monte_SANS(imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,func,$coefStr,results) 60 End 32 //  make sure that the ratio of scattering coherent/incoherent is properly adjusted for the sample composition 33 // or the volume fraction of solvent. 34 // 35 //  add to the results the fraction of coherently scattered neutrons that are singly scattered, different than 36 // the overall fraction of singly scattered, and maybe more important to know. 37 // 38 //  change the fraction reaching the detector to exclude those that don't interact. These transmitted neutrons 39 // aren't counted. Is the # that interact a better number? 40 // 41 //  do we want to NOT offset the data by a multiplicative factor as it is "frozen" , so that the 42 // effects on the absolute scale can be seen? 43 // 44 //  why is "pure" incoherent scattering giving me a q^1 slope, even with the detector all the way back? 45 //  can I speed up by assuming everything interacts? This would compromise the ability to calculate multiple scattering 46 //  ask John how to verify what is going on 47 //  a number of models are now found to be illbehaved when q=1e10. Then the random deviate calculation blows up. 48 // a warning has been added  but the models are better fixed with the limiting value. 49 // 50 // 51 52 61 53 62 54 ////////// … … 80 72 // 81 73 82 Function Monte_SANS(imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,func,coef,results) 83 Variable imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh 84 FUNCREF SANSModelAAO_MCproto func 85 WAVE coef,results 86 74 Function Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,MC_linear_data,results) 75 WAVE inputWave,ran_dev,nt,j1,j2,nn,MC_linear_data,results 76 77 Variable imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,sig_sas 87 78 Variable NUM_BINS,N_INDEX 88 79 Variable RHO,SIGSAS,SIGABS_0 … … 94 85 // at 0, making things much easier for me... 95 86 //DIMENSION NT(0:5000),J1(0:5000),J2(0:5000),NN(0:100) 96 Make/O/N=5000 root:Packages:NIST:SAS:nt,root:Packages:NIST:SAS:j1,root:Packages:NIST:SAS:j2 97 Make/O/N=100 root:Packages:NIST:SAS:nn 98 WAVE nt = root:Packages:NIST:SAS:nt 99 WAVE j1 = root:Packages:NIST:SAS:j1 100 WAVE j2 = root:Packages:NIST:SAS:j2 101 WAVE nn = root:Packages:NIST:SAS:nn 102 87 103 88 Variable G0,E_NT,E_NN,TRANS_th,Trans_exp,rat 104 89 Variable GG,GG_ED,dS_dW,ds_dw_double,ds_dw_single 105 90 Variable DONE,FIND_THETA,err //used as logicals 106 91 107 Variable Vx,Vy,Vz,Theta_z,qq ,SIG_SAS92 Variable Vx,Vy,Vz,Theta_z,qq 108 93 Variable Sig_scat,Sig_abs,Ratio,Sig_total 109 94 Variable isOn=0,testQ,testPhi,xPixel,yPixel 110 95 111 // detector information 112 // Variable pixSize,sdd,xCtr,yCtr 113 // pixSize = 0.5 //cm 114 // sdd = 400 //cm 115 // xCtr = 100 //pixels 116 // yCtr = 64 117 96 imon = inputWave[0] 97 r1 = inputWave[1] 98 r2 = inputWave[2] 99 xCtr = inputWave[3] 100 yCtr = inputWave[4] 101 sdd = inputWave[5] 102 pixSize = inputWave[6] 103 thick = inputWave[7] 104 wavelength = inputWave[8] 105 sig_incoh = inputWave[9] 106 sig_sas = inputWave[10] 107 118 108 // SetRandomSeed 0.1 //to get a reproduceable sequence 119 109 120 110 //scattering power and maximum qvalue to bin 121 111 // zpow = .1 //scattering power, calculated below 122 qmax = 0.1 //maximum Q to bin 1D data. (A1) (not really used)112 qmax = 4*pi/wavelength //maximum Q to bin 1D data. (A1) (not really used, so set to a big value) 123 113 sigabs_0 = 0.0 // ignore absorption cross section/wavelength [1/(cm A)] 124 114 N_INDEX = 50 // maximum number of scattering events per neutron … … 128 118 // and calculate the scattering power from the model function 129 119 // 130 CalculateRandomDeviate(func,coef,wavelength,"root:Packages:NIST:SAS:ran_dev",SIG_SAS)131 WAVE ran_dev=$"root:Packages:NIST:SAS:ran_dev"132 120 Variable left = leftx(ran_dev) 133 121 Variable delta = deltax(ran_dev) 134 135 // arrays for the data  these should already be there, created by SASCALC initializing the space136 Make/O/D/N=(128,128) root:Packages:NIST:SAS:data,root:Packages:NIST:SAS:linear_data137 WAVE MC_data = root:Packages:NIST:SAS:data138 WAVE MC_linear_data = root:Packages:NIST:SAS:linear_data139 MC_data = 0 //initialize140 MC_linear_data = 0141 122 142 123 //c total SAS crosssection … … 156 137 // RATIO = SIG_ABS / SIG_TOTAL 157 138 RATIO = sig_incoh / SIG_TOTAL 158 //!!!! the ratio is not ye r porperly weighted for the volume fractions of each component!!!139 //!!!! the ratio is not yet properly weighted for the volume fractions of each component!!! 159 140 160 141 //c assuming theta = sin(theta)...OK … … 174 155 nt = 0 175 156 nn=0 176 177 //vector to carry direction information178 Make/O/N=3 root:Packages:NIST:SAS:d_vector179 WAVE d_vector = root:Packages:NIST:SAS:d_vector180 157 181 158 //C MONITOR LOOP  looping over the number of incedent neutrons … … 187 164 Vy = 0.0 188 165 Vz = 1.0 189 d_vector[0] = 0190 d_vector[1] = 0191 d_vector[2] = 1192 166 193 167 Theta = 0.0 // Initialize scattering angle. … … 197 171 INDEX = 0 // Set counter for number of scattering events. 198 172 zz = 0.0 // Set entering dimension of sample. 199 //RR = 2.0*R1 // Make sure next loop runs at least once.200 173 201 174 do // Makes sure position is within circle. … … 212 185 ll = PATH_len(ran,Sig_total) 213 186 //Determine new scattering direction vector. 214 err = NewDirection(d_vector,Theta,Phi) //d_vector is updated, theta, phi unchanged by function 215 vx = d_vector[0] //reassign to local variables 216 vy = d_vector[1] 217 vz = d_vector[2] 187 err = NewDirection(vx,vy,vz,Theta,Phi) //vx,vy,vz is updated, theta, phi unchanged by function 188 218 189 //X,Y,ZPOSITION OF SCATTERING EVENT. 219 190 xx += ll*vx … … 223 194 224 195 //Check whether interaction occurred within sample volume. 225 IF (((zz > 0.0) %& (zz < THICK)) %& (rr < r2))196 IF (((zz > 0.0) && (zz < THICK)) && (rr < r2)) 226 197 //NEUTRON INTERACTED. 227 198 INDEX += 1 //Increment counter of scattering events. … … 234 205 FIND_THETA = 0 //false 235 206 DO 236 ran = abs(enoise(1)) //[0,1] 237 207 //ran = abs(enoise(1)) //[0,1] 238 208 //theta = Scat_angle(Ran,R_DAB,wavelength) // CHOOSE DAB ANGLE  this is 2Theta 239 //theta = Scat_angle(Ran,100,wavelength) // CHOOSE DAB ANGLE  this is 2Theta240 209 //Q0 = 2*PI*THETA/WAVElength // John chose theta, calculated Q 241 210 242 211 // pick a qvalue from the deviate function 243 212 // pnt2x truncates the point to an integer before returning the x 244 //Q0 = pnt2x(ran_dev,binarysearchinterp(ran_dev,abs(enoise(1))))245 213 // so get it from the wave scaling instead 246 214 Q0 =left + binarysearchinterp(ran_dev,abs(enoise(1)))*delta 247 215 theta = Q0/2/Pi*wavelength //SAS approximation 248 216 249 // always accept 250 FIND_THETA = 1 251 // ran = abs(enoise(1)) //[0,1] 252 // BOUND = inten_sphere(Q0,R0) / inten_dab(Q0,R_DAB) 253 // IF(BOUND > 1.0) 254 // Print "****BOUND ERROR.." 255 // ENDIF 256 // IF(Ran <= BOUND) 257 // FIND_THETA = 1 //accept attempt 258 // ENDIF 259 // 217 //Print "q0, theta = ",q0,theta 218 219 FIND_THETA = 1 //always accept 220 260 221 while(!find_theta) 261 222 ran = abs(enoise(1)) //[0,1] … … 289 250 // since the scattering is isotropic, I can safely pick a new, random value 290 251 // this would not be true if simulating anisotropic scattering. 291 testPhi = abs(enoise(1 ,2))*2*Pi252 testPhi = abs(enoise(1))*2*Pi 292 253 // is it on the detector? 293 254 FindPixel(testQ,testPhi,wavelength,sdd,pixSize,xCtr,yCtr,xPixel,yPixel) 255 256 // if(xPixel != xCtr && yPixel != yCtr) 257 // Print "testQ,testPhi,xPixel,yPixel",testQ,testPhi,xPixel,yPixel 258 // endif 259 294 260 if(xPixel != 1 && yPixel != 1) 295 261 isOn += 1 … … 328 294 results[3] = isOn/iMon 329 295 330 MC_data = MC_linear_data331 // MC_data = log(MC_data)332 // for display ONLY, since most of the neutrons just pass through with theta=0333 334 MC_linear_data[xCtr][yCtr]=0335 MC_data[xCtr][yCtr]=0336 337 296 338 297 // OUTPUT of the 1D data, not necessary now since I want the 2D … … 413 372 WAVE Gq = root:Packages:NIST:SAS:gQ 414 373 WAVE xw = root:Packages:NIST:SAS:xw 415 // Make/O/N=(nPts_ran)/D iWave1,iWave2 416 SetScale/I x (0+1e10),qu*(11e10),"", Gq,xw //don't start at zero or run up all the way to qu to avoid numerical errors 374 SetScale/I x (0+1e6),qu*(11e10),"", Gq,xw //don't start at zero or run up all the way to qu to avoid numerical errors 417 375 xw=x //for the AAO 418 376 func(coef,Gq,xw) //call as AAO … … 428 386 429 387 Duplicate/O Gq_INT $outWave 430 431 //Display Gq_INT432 // SetScale/I x 0,qu*(11e10),"", iWave2 //qu = 4Pi/lam (4Pi/6 = 2.09)433 //// iWave2 = DAB_modelX(coef_dab,x)434 // iWave2 = Peak_Gauss_modelX(coef_peak_gauss,x)435 // duplicate/O iWave2 Gqu436 //// Gqu = x*Gqu437 // Gqu = Gqu*sin(2*asin(x/qu))/sqrt(1(x/qu))438 // Integrate/METH=1 Gqu/D=Gqu_INT439 // //Appendtograph Gqu_INT440 // //ModifyGraph lsize(Gq_INT)=3441 // Gq_INT /= gqu_int[nPts_ran2]442 388 443 389 return(0) … … 539 485 tmp = "ASStandardFunction;Ann_1D_Graph;Avg_1D_Graph;BStandardFunction;CStandardFunction;Draw_Plot1D;MyMat2XYZ;NewDirection;SANSModelAAO_MCproto;" 540 486 list = RemoveFromList(tmp, list ,";") 487 list = RemoveFromList("Monte_SANS", list) 541 488 542 489 tmp = FunctionList("f*",";","NPARAMS:2") //point calculations … … 560 507 list = SortList(list) 561 508 return(list) 562 End 563 564 565 // CALCULATES SCATTERING FOR SPHERES, WITH I(0) = 1.0 566 //Function inten_sphere(qq,R0) 567 // Variable qq,r0 568 // 569 // Variable xx,i_sphere 570 // xx = qq*R0 571 // IF (xx < 1.0e6) 572 // I_SPHERE = 1.0 573 // ELSE 574 // I_SPHERE = 9.0*( (SIN(xx)xx*COS(xx)) / xx^3 )^2 575 // ENDIF 576 // RETURN (i_sphere) 577 //END 578 579 580 //CALCULATES SCATTERING FOR DAB MODEL, WITH I(0) = 1.0 581 //FUNCTION inten_dab(qq,R_DAB) 582 // Variable qq,r_dab 583 // 584 // Variable i_dab 585 // 586 // I_DAB = 1.0 / (1.0+(qq*R_DAB)^2)^2 587 // RETURN (i_dab) 588 //END 509 End 589 510 590 511 … … 605 526 //calculates new direction (xyz) from an old direction 606 527 //theta and phi don't change 607 Function NewDirection( d_vector,theta,phi)608 Wave d_vector528 Function NewDirection(vx,vy,vz,theta,phi) 529 Variable &vx,&vy,&vz 609 530 Variable theta,phi 610 531 611 532 Variable err=0,vx0,vy0,vz0 612 Variable Vx,Vy,Vz,nx,ny,mag_xy,tx,ty,tz 613 Vx = d_vector[0] 614 Vy = d_vector[1] 615 Vz = d_vector[2] 533 Variable nx,ny,mag_xy,tx,ty,tz 616 534 617 535 //store old direction vector … … 641 559 Vz = cos(phi)*sin(theta)*Tz + cos(theta)*Vz0 642 560 643 //reassign d_vector before exiting644 d_vector[0] = vx645 d_vector[1] = vy646 d_vector[2] = vz647 648 561 Return(err) 649 562 End … … 654 567 Variable retval 655 568 656 retval = 1*l og(1aval)/sig_tot569 retval = 1*ln(1aval)/sig_tot 657 570 658 571 return(retval) … … 661 574 Window MC_SASCALC() : Panel 662 575 PauseUpdate; Silent 1 // building window... 663 NewPanel /W=(787,44,1088,563) /K=1/N=MC_SASCALC as "SANS Simulator"576 NewPanel /W=(787,44,1088,563) /N=MC_SASCALC as "SANS Simulator" 664 577 CheckBox MC_check0,pos={11,11},size={98,14},title="Use MC Simulation" 665 578 CheckBox MC_check0,variable= root:Packages:NIST:SAS:gDoMonteCarlo … … 674 587 PopupMenu MC_popup0,pos={11,63},size={162,20},proc=MC_ModelPopMenuProc,title="Model Function" 675 588 PopupMenu MC_popup0,mode=1,value= #"MC_FunctionPopupList()" 589 Button MC_button0,pos={17,455},size={50,20},proc=MC_DoItButtonProc,title="Do It" 590 Button MC_button1,pos={15,484},size={80,20},proc=MC_Display2DButtonProc,title="Show 2D" 676 591 677 592 SetDataFolder root:Packages:NIST:SAS: … … 698 613 return 0 699 614 End 615 616 Function MC_DoItButtonProc(ba) : ButtonControl 617 STRUCT WMButtonAction &ba 618 619 switch( ba.eventCode ) 620 case 2: // mouse up 621 // click code here 622 ReCalculateInten(1) 623 break 624 endswitch 625 626 return 0 627 End 628 629 630 Function MC_Display2DButtonProc(ba) : ButtonControl 631 STRUCT WMButtonAction &ba 632 633 switch( ba.eventCode ) 634 case 2: // mouse up 635 // click code here 636 Execute "ChangeDisplay(\"SAS\")" 637 break 638 endswitch 639 640 return 0 641 End 642 643 /////UNUSED, testing routines that have note been updated to work with SASCALC 644 // 645 //Macro Simulate2D_MonteCarlo(imon,r1,r2,xCtr,yCtr,sdd,thick,wavelength,sig_incoh,funcStr) 646 // Variable imon=100000,r1=0.6,r2=0.8,xCtr=100,yCtr=64,sdd=400,thick=0.1,wavelength=6,sig_incoh=0.1 647 // String funcStr 648 // Prompt funcStr, "Pick the model function", popup, MC_FunctionPopupList() 649 // 650 // String coefStr = MC_getFunctionCoef(funcStr) 651 // Variable pixSize = 0.5 // can't have 11 parameters in macro! 652 // 653 // if(!MC_CheckFunctionAndCoef(funcStr,coefStr)) 654 // Abort "The coefficients and function type do not match. Please correct the selections in the popup menus." 655 // endif 656 // 657 // Make/O/D/N=10 root:Packages:NIST:SAS:results 658 // Make/O/T/N=10 root:Packages:NIST:SAS:results_desc 659 // 660 // RunMonte(imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,funcStr,coefStr,results) 661 // 662 //End 663 // 664 //Function RunMonte(imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,funcStr,coefStr) 665 // Variable imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh 666 // String funcStr,coefStr 667 // WAVE results 668 // 669 // FUNCREF SANSModelAAO_MCproto func=$funcStr 670 // 671 // Monte_SANS(imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,sig_sas,ran_dev,linear_data,results) 672 //End 673 // 674 ////// END UNUSED BLOCK
Note: See TracChangeset
for help on using the changeset viewer.