 Timestamp:
 Nov 3, 2008 12:48:40 PM (14 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/MultScatter_MonteCarlo_2D.ipf
r431 r434 29 29 //  settings for "count for X seconds" or "how long to 1E6 cts on detector" (run short sim, then multiply) 30 30 //  add quartz window scattering to the simulation somehow 31 //  do smeared models make any sense?? 31 //  do smeared models make any sense?? Yes, John agrees that they do, and may be used in a more realistic simulation 32 32 //  make sure that the ratio of scattering coherent/incoherent is properly adjusted for the sample composition 33 33 // or the volume fraction of solvent. … … 42 42 // effects on the absolute scale can be seen? 43 43 // 44 // why is "pure" incoherent scattering giving me a q^1 slope, even with the detector all the way back?44 // X why is "pure" incoherent scattering giving me a q^1 slope, even with the detector all the way back? 45 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 on46 // X ask John how to verify what is going on 47 47 //  a number of models are now found to be illbehaved when q=1e10. Then the random deviate calculation blows up. 48 48 // a warning has been added  but the models are better fixed with the limiting value. 49 49 // 50 50 // 51 52 // threaded call to the main function, adds up the individual runs, and returns what is to be displayed 53 // results is calculated and sent back for display 54 Function Monte_SANS_Threaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) 55 WAVE inputWave,ran_dev,nt,j1,j2,nn,linear_data,results 56 57 //initialize ran1 in the XOP by passing a negative integer 58 // does nothing in the Igor code 59 Duplicate/O results retWave 60 //results[0] = 1*(datetime) 61 62 Variable NNeutron=inputWave[0] 63 Variable i,nthreads= ThreadProcessorCount 64 if(nthreads>2) //only support 2 processors until I can figure out how to properly thread the XOP and to loop it 65 nthreads=2 66 endif 67 68 // nthreads = 1 69 70 variable mt= ThreadGroupCreate(nthreads) 71 72 inputWave[0] = NNeutron/nthreads //split up the number of neutrons 73 74 for(i=0;i<nthreads;i+=1) 75 Duplicate/O nt $("nt"+num2istr(i)) //new instance for each thread 76 Duplicate/O j1 $("j1"+num2istr(i)) 77 Duplicate/O j2 $("j2"+num2istr(i)) 78 Duplicate/O nn $("nn"+num2istr(i)) 79 Duplicate/O linear_data $("linear_data"+num2istr(i)) 80 Duplicate/O retWave $("retWave"+num2istr(i)) 81 Duplicate/O inputWave $("inputWave"+num2istr(i)) 82 Duplicate/O ran_dev $("ran_dev"+num2istr(i)) 83 84 // ?? I need explicit wave references? 85 if(i==0) 86 WAVE inputWave0,ran_dev0,nt0,j10,j20,nn0,linear_data0,retWave0 87 retWave0[0] = 1*datetime //to initialize ran3 88 ThreadStart mt,i,Monte_SANS_W1(inputWave0,ran_dev0,nt0,j10,j20,nn0,linear_data0,retWave0) 89 Print "started thread 0" 90 endif 91 if(i==1) 92 WAVE inputWave1,ran_dev1,nt1,j11,j21,nn1,linear_data1,retWave1 93 //retWave1[0] = 1*datetime //to initialize ran3 94 ThreadStart mt,i,Monte_SANS_W1(inputWave1,ran_dev1,nt1,j11,j21,nn1,linear_data1,retWave1) 95 Print "started thread 1" 96 endif 97 // if(i==2) 98 // WAVE inputWave2,ran_dev2,nt2,j12,j22,nn2,linear_data2,retWave2 99 // retWave2[0] = 1*datetime //to initialize ran3 100 // ThreadStart mt,i,Monte_SANS_W(inputWave2,ran_dev2,nt2,j12,j22,nn2,linear_data2,retWave2) 101 // endif 102 // if(i==3) 103 // WAVE inputWave3,ran_dev3,nt3,j13,j23,nn3,linear_data3,retWave3 104 // retWave3[0] = 1*datetime //to initialize ran3 105 // ThreadStart mt,i,Monte_SANS_W(inputWave3,ran_dev3,nt3,j13,j23,nn3,linear_data3,retWave3) 106 // endif 107 endfor 108 109 // wait until done 110 do 111 variable tgs= ThreadGroupWait(mt,100) 112 while( tgs != 0 ) 113 variable dummy= ThreadGroupRelease(mt) 114 Print "done with all threads" 115 116 // calculate all of the bits for the results 117 if(nthreads == 1) 118 nt = nt0 // add up each instance 119 j1 = j10 120 j2 = j20 121 nn = nn0 122 linear_data = linear_data0 123 retWave = retWave0 124 endif 125 if(nthreads == 2) 126 nt = nt0+nt1 // add up each instance 127 j1 = j10+j11 128 j2 = j20+j21 129 nn = nn0+nn1 130 linear_data = linear_data0+linear_data1 131 retWave = retWave0+retWave1 132 endif 133 // if(nthreads == 3) 134 // nt = nt0+nt1+nt2 // add up each instance 135 // j1 = j10+j11+j12 136 // j2 = j20+j21+j22 137 // nn = nn0+nn1+nn2 138 // linear_data = linear_data0+linear_data1+linear_data2 139 // retWave = retWave0+retWave1+retWave2 140 // endif 141 // if(nthreads == 4) 142 // nt = nt0+nt1+nt2+nt3 // add up each instance 143 // j1 = j10+j11+j12+j13 144 // j2 = j20+j21+j22+j23 145 // nn = nn0+nn1+nn2+nn3 146 // linear_data = linear_data0+linear_data1+linear_data2+linear_data3 147 // retWave = retWave0+retWave1+retWave2+retWave3 148 // endif 149 150 // fill up the results wave 151 Variable xc,yc 152 xc=inputWave[3] 153 yc=inputWave[4] 154 results[0] = inputWave[9]+inputWave[10] //total XS 155 results[1] = inputWave[10] //SAS XS 156 results[2] = retWave[1] //number that interact n2 157 results[3] = retWave[2]  linear_data[xc][yc] //# reaching detector minus Q(0) 158 results[4] = retWave[3]/retWave[1] //avg# times scattered 159 results[5] = retWave[4]/retWave[1] //single coherent fraction 160 results[6] = retWave[5]/retWave[1] //double coherent fraction 161 results[7] = retWave[6]/retWave[1] //multiple scatter fraction 162 results[8] = (retWave[0]retWave[1])/retWave[0] //trasnmitted fraction 163 164 return(0) 165 End 166 167 // worker function for threads, does nothing except switch between XOP and Igor versions 168 ThreadSafe Function Monte_SANS_W1(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) 169 WAVE inputWave,ran_dev,nt,j1,j2,nn,linear_data,results 170 171 #if exists("Monte_SANSX") 172 Monte_SANSX(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) 173 #else 174 Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) 175 #endif 176 177 return (0) 178 End 179 // worker function for threads, does nothing except switch between XOP and Igor versions 180 ThreadSafe Function Monte_SANS_W2(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) 181 WAVE inputWave,ran_dev,nt,j1,j2,nn,linear_data,results 182 183 #if exists("xxxxMonte_SANSX") 184 Monte_SANSX(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) 185 #else 186 Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) 187 #endif 188 189 return (0) 190 End 191 192 // NONthreaded call to the main function returns what is to be displayed 193 // results is calculated and sent back for display 194 Function Monte_SANS_NotThreaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) 195 WAVE inputWave,ran_dev,nt,j1,j2,nn,linear_data,results 196 197 //initialize ran1 in the XOP by passing a negative integer 198 // does nothing in the Igor code, enoise is already initialized 199 Duplicate/O results retWave 200 WAVE retWave 201 retWave[0] = 1*abs(trunc(100000*enoise(1))) 202 203 #if exists("Monte_SANSX") 204 Monte_SANSX(inputWave,ran_dev,nt,j1,j2,nn,linear_data,retWave) 205 #else 206 Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,retWave) 207 #endif 208 209 // fill up the results wave 210 Variable xc,yc 211 xc=inputWave[3] 212 yc=inputWave[4] 213 results[0] = inputWave[9]+inputWave[10] //total XS 214 results[1] = inputWave[10] //SAS XS 215 results[2] = retWave[1] //number that interact n2 216 results[3] = retWave[2]  linear_data[xc][yc] //# reaching detector minus Q(0) 217 results[4] = retWave[3]/retWave[1] //avg# times scattered 218 results[5] = retWave[4]/retWave[1] //single coherent fraction 219 results[6] = retWave[5]/retWave[1] //double coherent fraction 220 results[7] = retWave[6]/retWave[1] //multiple scatter fraction 221 results[8] = (retWave[0]retWave[1])/retWave[0] //trasnmitted fraction 222 223 return(0) 224 End 51 225 52 226 … … 72 246 // 73 247 74 Function Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,MC_linear_data,results)248 ThreadSafe Function Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,MC_linear_data,results) 75 249 WAVE inputWave,ran_dev,nt,j1,j2,nn,MC_linear_data,results 76 250 … … 82 256 Variable N1,N2,N3,DTH,zz,tt,SIG_SINGLE,xx,yy,PHI,UU,SIG 83 257 Variable THETA,Ran,ll,D_OMEGA,RR,Tabs,Ttot,I1_sumI 84 //in John's implementation, he dimensioned the indices of the arrays to begin85 // at 0, making things much easier for me...86 //DIMENSION NT(0:5000),J1(0:5000),J2(0:5000),NN(0:100)87 88 258 Variable G0,E_NT,E_NN,TRANS_th,Trans_exp,rat 89 259 Variable GG,GG_ED,dS_dW,ds_dw_double,ds_dw_single … … 93 263 Variable Sig_scat,Sig_abs,Ratio,Sig_total 94 264 Variable isOn=0,testQ,testPhi,xPixel,yPixel 265 Variable NSingleIncoherent,NSingleCoherent,NScatterEvents,incoherentEvent,coherentEvent 266 Variable NDoubleCoherent,NMultipleScatter 95 267 96 268 imon = inputWave[0] … … 122 294 123 295 //c total SAS crosssection 124 //125 // input a test value for the incoherent scattering from water126 //127 // sig_incoh = 5.6 //[1/cm] as calculated for H2O, now a parameter128 //129 296 // SIG_SAS = zpow/thick 130 297 zpow = sig_sas*thick //since I now calculate the sig_sas from the model … … 133 300 // Print "The TOTAL XSECTION. (CM1) is ",sig_total 134 301 // Print "The TOTAL SAS XSECTION. (CM1) is ",sig_sas 135 results[0] = sig_total 136 results[1] = sig_sas302 // results[0] = sig_total //assign these after everything's done 303 // results[1] = sig_sas 137 304 // RATIO = SIG_ABS / SIG_TOTAL 138 305 RATIO = sig_incoh / SIG_TOTAL 139 //!!!! the ratio is not yet properly weighted for the volume fractions of each component!!!140 306 141 307 //c assuming theta = sin(theta)...OK … … 149 315 N2 = 0 150 316 N3 = 0 317 NSingleIncoherent = 0 318 NSingleCoherent = 0 319 NDoubleCoherent = 0 320 NMultipleScatter = 0 321 NScatterEvents = 0 151 322 152 323 //C INITIALIZE ARRAYS. … … 171 342 INDEX = 0 // Set counter for number of scattering events. 172 343 zz = 0.0 // Set entering dimension of sample. 344 incoherentEvent = 0 345 coherentEvent = 0 173 346 174 347 do // Makes sure position is within circle. … … 203 376 //Split neutron interactions into scattering and absorption events 204 377 IF(ran > ratio ) //C NEUTRON SCATTERED coherently 378 coherentEvent = 1 205 379 FIND_THETA = 0 //false 206 380 DO … … 228 402 // phi and theta are random over the entire sphere of scattering 229 403 // !can't just choose random theta and phi, won't be random over sphere solid angle 404 incoherentEvent = 1 230 405 231 406 ran = abs(enoise(1)) //[0,1] … … 255 430 FindPixel(testQ,testPhi,wavelength,sdd,pixSize,xCtr,yCtr,xPixel,yPixel) 256 431 257 // if(xPixel != xCtr && yPixel != yCtr)258 // Print "testQ,testPhi,xPixel,yPixel",testQ,testPhi,xPixel,yPixel259 // endif260 261 432 if(xPixel != 1 && yPixel != 1) 262 isOn += 1263 433 //if(index==1) // only the single scattering events 264 434 MC_linear_data[xPixel][yPixel] += 1 //this is the total scattering, including multiple scattering 265 435 //endif 436 isOn += 1 // scattered neutron that lands on detector 266 437 endif 267 438 … … 281 452 EndIf 282 453 454 // increment all of the counters now since done==1 here and I'm sure to exit and get another neutron 455 NScatterEvents += index //total number of scattering events 456 if(index == 1 && incoherentEvent == 1) 457 NSingleIncoherent += 1 458 endif 459 if(index == 1 && coherentEvent == 1) 460 NSingleCoherent += 1 461 endif 462 if(index == 2 && coherentEvent == 1 && incoherentEvent == 0) 463 NDoubleCoherent += 1 464 endif 465 if(index > 1) 466 NMultipleScatter += 1 467 endif 283 468 ENDIF 284 469 while (!done) … … 286 471 287 472 // Print "Monte Carlo Done" 288 trans_th = exp(sig_total*thick) 473 results[0] = n1 474 results[1] = n2 475 results[2] = isOn 476 results[3] = NScatterEvents //sum of # of times that neutrons scattered 477 results[4] = NSingleCoherent //# of events that are single, coherent 478 results[5] = NDoubleCoherent 479 results[6] = NMultipleScatter //# of multiple scattering events 480 481 482 // trans_th = exp(sig_total*thick) 289 483 // TRANS_exp = (N1N2) / N1 // Transmission 290 484 // dsigma/domega assuming isotropic scattering, with no absorption. … … 292 486 // Print "total # of neutrons reaching 2D detector",isOn 293 487 // Print "fraction of incident neutrons reaching detector ",isOn/iMon 294 results[2] = isOn 295 results[3] = isOn/iMon 296 297 298 // OUTPUT of the 1D data, not necessary now since I want the 2D 299 // Make/O/N=(num_bins) qvals,int_ms,sig_ms,int_sing,int_doub 300 // ii=1 301 // Print "binning" 302 // do 303 // //CALCULATE MEAN THETA IN BIN. 304 // THETA_z = (ii0.5)*DTH // Mean scattering angle of bin. 305 // //Solid angle of Ith bin. 306 // D_OMEGA = 2*PI*ABS( COS(ii*DTH)  COS((ii1)*DTH) ) 307 // //SOLID ANGLE CORRECTION: YIELDING CROSSSECTION. 308 // dS_dW = NT[ii]/(N1*THICK*Trans_th*D_OMEGA) 309 // SIG = NT[ii]^0.5/(N1*THICK*Trans_th*D_OMEGA) 310 // ds_dw_single = J1[ii]/(N1*THICK*Trans_th*D_OMEGA) 311 // ds_dw_double = J2[ii]/(N1*THICK*Trans_th*D_OMEGA) 312 // //Deviation from isotropic model. 313 // qq = 4*pi*sin(0.5*theta_z)/wavelength 314 // qvals[ii1] = qq 315 // int_ms[ii1] = dS_dW 316 // sig_ms[ii1] = sig 317 // int_sing[ii1] = ds_dw_single 318 // int_doub[ii1] = ds_dw_double 319 // //140 WRITE (7,145) qq,dS_dW,SIG,ds_dw_single,ds_dw_double 320 // ii+=1 321 // while(ii<=num_bins) 322 // Print "done binning" 323 324 325 // Write(7,*) 326 // Write(7,*) '#Times Sc. #Events ' 327 // DO 150 I = 1,N_INDEX 328 //150 WRITE (7,146) I,NN(I) 329 //146 Format (I5,T10,F8.0) 330 331 /// Write(7,171) N1 332 // Write(7,172) N2 333 // Write(7,173) N3 334 // Write(7,174) N2N3 335 336 //171 Format('Total number of neutrons: N1= ',E10.5) 337 ///172 Format('Number of neutrons that interact: N2= ',E10.5) 338 //173 Format('Number of absorption events: N3= ',E10.5) 339 //174 Format('# of neutrons that scatter out:(N2N3)= ',E10.5) 340 488 341 489 // Print "Total number of neutrons = ",N1 342 490 // Print "Total number of neutrons that interact = ",N2 343 491 // Print "Fraction of singly scattered neutrons = ",sum(j1,inf,inf)/N2 344 results[4] = N2 345 results[5] = sum(j1,inf,inf)/N2 346 347 Tabs = (N1N3)/N1 348 Ttot = (N1N2)/N1 349 I1_sumI = NN[0]/(N2N3) 492 // results[2] = N2 //number that scatter 493 // results[3] = isOn  MC_linear_data[xCtr][yCtr] //# scattered reaching detector minus zero angle 494 495 496 // Tabs = (N1N3)/N1 497 // Ttot = (N1N2)/N1 498 // I1_sumI = NN[0]/(N2N3) 350 499 // Print "Tabs = ",Tabs 351 500 // Print "Transmitted neutrons = ",Ttot 352 results[6] = Ttot501 // results[8] = Ttot 353 502 // Print "I1 / all I1 = ", I1_sumI 354 // Print "DONE!"355 503 356 504 End … … 391 539 End 392 540 393 Function FindPixel(testQ,testPhi,lam,sdd,pixSize,xCtr,yCtr,xPixel,yPixel)541 ThreadSafe Function FindPixel(testQ,testPhi,lam,sdd,pixSize,xCtr,yCtr,xPixel,yPixel) 394 542 Variable testQ,testPhi,lam,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel 395 543 … … 484 632 485 633 // SANS Reduction bits 486 tmp = "ASStandardFunction;Ann_1D_Graph;Avg_1D_Graph;BStandardFunction;CStandardFunction;Draw_Plot1D;MyMat2XYZ;NewDirection;SANSModelAAO_MCproto; "634 tmp = "ASStandardFunction;Ann_1D_Graph;Avg_1D_Graph;BStandardFunction;CStandardFunction;Draw_Plot1D;MyMat2XYZ;NewDirection;SANSModelAAO_MCproto;Monte_SANS_Threaded;Monte_SANS_NotThreaded;Monte_SANS_W1;Monte_SANS_W2;" 487 635 list = RemoveFromList(tmp, list ,";") 488 636 list = RemoveFromList("Monte_SANS", list) … … 507 655 508 656 list = SortList(list) 657 658 list = "default;"+list 509 659 return(list) 510 660 End … … 527 677 //calculates new direction (xyz) from an old direction 528 678 //theta and phi don't change 529 Function NewDirection(vx,vy,vz,theta,phi)679 ThreadSafe Function NewDirection(vx,vy,vz,theta,phi) 530 680 Variable &vx,&vy,&vz 531 681 Variable theta,phi … … 563 713 End 564 714 565 Function path_len(aval,sig_tot)715 ThreadSafe Function path_len(aval,sig_tot) 566 716 Variable aval,sig_tot 567 717 … … 642 792 End 643 793 644 /////UNUSED, testing routines that have not ebeen updated to work with SASCALC794 /////UNUSED, testing routines that have not been updated to work with SASCALC 645 795 // 646 796 //Macro Simulate2D_MonteCarlo(imon,r1,r2,xCtr,yCtr,sdd,thick,wavelength,sig_incoh,funcStr)
Note: See TracChangeset
for help on using the changeset viewer.