Changeset 535
 Timestamp:
 Jul 7, 2009 6:35:07 PM (14 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/MultScatter_MonteCarlo_2D.ipf
r494 r535 756 756 757 757 // MOTOFIT/GenFit bits 758 tmp = "GEN_allatoncefitfunc;GEN_fitfunc;GetCheckBoxesState;MOTO_GFFitAllAtOnceTemplate;MOTO_GFFitFuncTemplate;MOTO_NewGF_SetXWaveInList;MOTO_NewGlblFitFunc;MOTO_NewGlblFitFuncAllAtOnce; "758 tmp = "GEN_allatoncefitfunc;GEN_fitfunc;GetCheckBoxesState;MOTO_GFFitAllAtOnceTemplate;MOTO_GFFitFuncTemplate;MOTO_NewGF_SetXWaveInList;MOTO_NewGlblFitFunc;MOTO_NewGlblFitFuncAllAtOnce;GeneticFit_UnSmearedModel;GeneticFit_SmearedModel;" 759 759 list = RemoveFromList(tmp, list ,";") 760 760 … … 1087 1087 1088 1088 1089 /////UNUSED, testing routines that have not been updated to work with SASCALC 1090 // 1091 //Macro Simulate2D_MonteCarlo(imon,r1,r2,xCtr,yCtr,sdd,thick,wavelength,sig_incoh,funcStr) 1092 // Variable imon=100000,r1=0.6,r2=0.8,xCtr=100,yCtr=64,sdd=400,thick=0.1,wavelength=6,sig_incoh=0.1 1093 // String funcStr 1094 // Prompt funcStr, "Pick the model function", popup, MC_FunctionPopupList() 1095 // 1096 // String coefStr = MC_getFunctionCoef(funcStr) 1097 // Variable pixSize = 0.5 // can't have 11 parameters in macro! 1098 // 1099 // if(!MC_CheckFunctionAndCoef(funcStr,coefStr)) 1100 // Abort "The coefficients and function type do not match. Please correct the selections in the popup menus." 1101 // endif 1102 // 1103 // Make/O/D/N=10 root:Packages:NIST:SAS:results 1104 // Make/O/T/N=10 root:Packages:NIST:SAS:results_desc 1105 // 1106 // RunMonte(imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,funcStr,coefStr,results) 1107 // 1089 /// called in SASCALC:ReCalculateInten() 1090 Function Simulate_2D_MC(funcStr,aveint,qval,sigave,sigmaq,qbar,fsubs) 1091 String funcStr 1092 WAVE aveint,qval,sigave,sigmaq,qbar,fsubs 1093 1094 NVAR doMonteCarlo = root:Packages:NIST:SAS:gDoMonteCarlo // == 1 if 2D MonteCarlo set by hidden flag 1095 WAVE rw=root:Packages:NIST:SAS:realsRead 1096 1097 NVAR imon = root:Packages:NIST:SAS:gImon 1098 NVAR thick = root:Packages:NIST:SAS:gThick 1099 NVAR sig_incoh = root:Packages:NIST:SAS:gSig_incoh 1100 NVAR r2 = root:Packages:NIST:SAS:gR2 1101 1102 // do the simulation here, or not 1103 Variable r1,xCtr,yCtr,sdd,pixSize,wavelength 1104 String coefStr,abortStr,str 1105 1106 r1 = rw[24]/2/10 // sample diameter convert diam in [mm] to radius in cm 1107 xCtr = rw[16] 1108 yCtr = rw[17] 1109 sdd = rw[18]*100 //conver header of [m] to [cm] 1110 pixSize = rw[10]/10 // convert pix size in mm to cm 1111 wavelength = rw[26] 1112 coefStr = MC_getFunctionCoef(funcStr) 1113 1114 if(!MC_CheckFunctionAndCoef(funcStr,coefStr)) 1115 doMonteCarlo = 0 //we're getting out now, reset the flag so we don't continually end up here 1116 Abort "The coefficients and function type do not match. Please correct the selections in the popup menus." 1117 endif 1118 1119 Variable sig_sas 1120 // FUNCREF SANSModelAAO_MCproto func=$("fSmeared"+funcStr) //a wrapper for the structure version 1121 FUNCREF SANSModelAAO_MCproto func=$(funcStr) //unsmeared 1122 WAVE results = root:Packages:NIST:SAS:results 1123 WAVE linear_data = root:Packages:NIST:SAS:linear_data 1124 WAVE data = root:Packages:NIST:SAS:data 1125 1126 results = 0 1127 linear_data = 0 1128 1129 CalculateRandomDeviate(func,$coefStr,wavelength,"root:Packages:NIST:SAS:ran_dev",SIG_SAS) 1130 if(sig_sas > 100) 1131 sprintf abortStr,"sig_sas = %g. Please check that the model coefficients have a zero background, or the low q is wellbehaved.",sig_sas 1132 Abort abortStr 1133 endif 1134 1135 WAVE ran_dev=$"root:Packages:NIST:SAS:ran_dev" 1136 1137 Make/O/D/N=5000 root:Packages:NIST:SAS:nt=0,root:Packages:NIST:SAS:j1=0,root:Packages:NIST:SAS:j2=0 1138 Make/O/D/N=100 root:Packages:NIST:SAS:nn=0 1139 Make/O/D/N=11 root:Packages:NIST:SAS:inputWave=0 1140 1141 WAVE nt = root:Packages:NIST:SAS:nt 1142 WAVE j1 = root:Packages:NIST:SAS:j1 1143 WAVE j2 = root:Packages:NIST:SAS:j2 1144 WAVE nn = root:Packages:NIST:SAS:nn 1145 WAVE inputWave = root:Packages:NIST:SAS:inputWave 1146 1147 inputWave[0] = imon 1148 inputWave[1] = r1 1149 inputWave[2] = r2 1150 inputWave[3] = xCtr 1151 inputWave[4] = yCtr 1152 inputWave[5] = sdd 1153 inputWave[6] = pixSize 1154 inputWave[7] = thick 1155 inputWave[8] = wavelength 1156 inputWave[9] = sig_incoh 1157 inputWave[10] = sig_sas 1158 1159 linear_data = 0 //initialize 1160 1161 Variable t0,trans 1162 1163 // get a time estimate, and give the user a chance to exit if they're unsure. 1164 t0 = stopMStimer(2) 1165 inputWave[0] = 1000 1166 NVAR useXOP = root:Packages:NIST:SAS:gUse_MC_XOP //if zero, will use nonthreaded Igor code 1167 1168 if(useXOP) 1169 //use a single thread, otherwise time is dominated by overhead 1170 Monte_SANS_NotThreaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) 1171 else 1172 Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) 1173 endif 1174 1175 t0 = (stopMSTimer(2)  t0)*1e6 1176 t0 *= imon/1000/ThreadProcessorCount //projected time, in seconds (using threads for the calculation) 1177 inputWave[0] = imon //reset 1178 1179 if(t0>10) 1180 sprintf str,"The simulation will take approximately %d seconds.\r Proceed?",t0 1181 DoAlert 1,str 1182 if(V_flag == 2) 1183 doMonteCarlo = 0 1184 reCalculateInten(1) //come back in and do the smeared calculation 1185 return(0) 1186 endif 1187 endif 1188 1189 linear_data = 0 //initialize 1190 // threading crashes!!  there must be some operation in the XOP that is not threadSafe. What, I don't know... 1191 // I think it's the ran() calls, being "nonreentrant". So the XOP now defines two separate functions, that each 1192 // use a different rng. This works. 1.75x speedup. 1193 t0 = stopMStimer(2) 1194 1195 if(useXOP) 1196 Monte_SANS_Threaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) 1197 else 1198 Monte_SANS_NotThreaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) 1199 endif 1200 1201 t0 = (stopMSTimer(2)  t0)*1e6 1202 Printf "MC sim time = %g seconds\r",t0 1203 1204 trans = results[8] //(n1n2)/n1 1205 if(trans == 0) 1206 trans = 1 1207 endif 1208 1209 Print "counts on detector, including transmitted = ",sum(linear_data,inf,inf) 1210 1211 // linear_data[xCtr][yCtr] = 0 //snip out the transmitted spike 1212 // Print "counts on detector not transmitted = ",sum(linear_data,inf,inf) 1213 1214 // or simulate a beamstop 1215 NVAR MC_BS_in = root:Packages:NIST:SAS:gBeamStopIn //if zero, beam stop is "out", as in a transmission measurement 1216 1217 Variable rad=beamstopDiam()/2 //beamstop radius in cm 1218 if(MC_BS_in) 1219 rad /= 0.5 //convert cm to pixels 1220 rad += 0. // (no  it cuts off the low Q artificially) add an extra pixel to each side to account for edge 1221 Duplicate/O linear_data,root:Packages:NIST:SAS:tmp_mask//,root:Packages:NIST:SAS:MC_linear_data 1222 WAVE tmp_mask = root:Packages:NIST:SAS:tmp_mask 1223 tmp_mask = (sqrt((pxCtr)^2+(qyCtr)^2) < rad) ? 0 : 1 //behind beamstop = 0, away = 1 1224 1225 linear_data *= tmp_mask 1226 endif 1227 1228 results[9] = sum(linear_data,inf,inf) 1229 // Print "counts on detector not behind beamstop = ",results[9] 1230 1231 // convert to absolute scale 1232 Variable kappa //,beaminten = beamIntensity() 1233 // kappa = beamInten*pi*r1*r1*thick*(pixSize/sdd)^2*trans*(iMon/beaminten) 1234 kappa = thick*(pixSize/sdd)^2*trans*iMon 1235 1236 //use kappa to get back to counts => linear_data = round(linear_data*kappa) 1237 Note/K linear_data ,"KAPPA="+num2str(kappa)+";" 1238 1239 NVAR rawCts = root:Packages:NIST:SAS:gRawCounts 1240 if(!rawCts) //go ahead and do the linear scaling 1241 linear_data = linear_data / kappa 1242 endif 1243 data = linear_data 1244 1245 // reaverage the 2D data 1246 S_CircularAverageTo1D("SAS") 1247 1248 // put the new result into the simulation folder 1249 Fake1DDataFolder(qval,aveint,sigave,sigmaQ,qbar,fSubs,"Simulation") 1250 1251 1252 return(0) 1253 end 1254 1255 1256 1257 // as a standalone panel, extra control bar (right) and subwindow implementations don't work right 1258 // for various reasons... 1259 Window Sim_1D_Panel() : Panel 1260 PauseUpdate; Silent 1 // building window... 1261 NewPanel /W=(92,556,713,818)/K=1 as "1D SANS Simulator" 1262 SetVariable MC_setvar0,pos={28,73},size={144,15},bodyWidth=80,title="# of neutrons" 1263 SetVariable MC_setvar0,format="%5.4g" 1264 SetVariable MC_setvar0,limits={0,inf,100},value= root:Packages:NIST:SAS:gImon 1265 SetVariable MC_setvar0_1,pos={28,119},size={131,15},bodyWidth=60,title="Thickness (cm)" 1266 SetVariable MC_setvar0_1,limits={0,inf,0.1},value= root:Packages:NIST:SAS:gThick 1267 SetVariable MC_setvar0_2,pos={28,96},size={149,15},bodyWidth=60,title="Incoherent XS (cm)" 1268 SetVariable MC_setvar0_2,limits={0,inf,0.1},value= root:Packages:NIST:SAS:gSig_incoh 1269 // SetVariable MC_setvar0_3,pos={28,142},size={150,15},bodyWidth=60,title="Sample Radius (cm)" 1270 // SetVariable MC_setvar0_3,limits={inf,inf,0.1},value= root:Packages:NIST:SAS:gR2 1271 SetVariable setvar0_3,pos={105,484},size={50,20},disable=1 1272 PopupMenu MC_popup0,pos={13,13},size={165,20},proc=Sim_1D_ModelPopMenuProc,title="Model Function" 1273 PopupMenu MC_popup0,mode=1,value= #"MC_FunctionPopupList()" 1274 // Button MC_button0,pos={17,181},size={130,20},proc=Sim_1D_DoItButtonProc,title="Do 1D Simulation" 1275 // Button MC_button0,fColor=(3,52428,1) 1276 // Button MC_button1,pos={17,208},size={80,20},proc=Sim_1D_Display2DButtonProc,title="Show 2D" 1277 GroupBox group0,pos={15,42},size={267,130},title="Sample Setup" 1278 SetVariable cntVar,pos={190,73},size={80,15},proc=Sim_1D_CountTimeSetVarProc,title="time(s)" 1279 SetVariable cntVar,format="%d" 1280 SetVariable cntVar,limits={1,600,1},value= root:Packages:NIST:SAS:gCntTime 1281 // Button MC_button2,pos={17,234},size={100,20},proc=SaveAsVAXButtonProc,title="Save 2D VAX" 1282 CheckBox check0,pos={216,180},size={68,14},title="Raw counts",variable = root:Packages:NIST:SAS:gRawCounts 1283 CheckBox check0_1,pos={216,199},size={60,14},title="Yes Offset",variable= root:Packages:NIST:SAS:gDoTraceOffset 1284 // CheckBox check0_2,pos={216,199+19},size={60,14},title="Beam Stop in",variable= root:Packages:NIST:SAS:gBeamStopIn 1285 // CheckBox check0_3,pos={216,199+2*19},size={60,14},title="use XOP",variable= root:Packages:NIST:SAS:gUse_MC_XOP 1286 1287 String fldrSav0= GetDataFolder(1) 1288 SetDataFolder root:Packages:NIST:SAS: 1289 Edit/W=(344,23,606,248)/HOST=# results_desc,results 1290 ModifyTable format(Point)=1,width(Point)=0,width(results_desc)=150 1291 SetDataFolder fldrSav0 1292 RenameWindow #,T_results 1293 SetActiveSubwindow ## 1294 1295 // set the flags here  do the simulation, but not 2D 1296 1297 root:Packages:NIST:SAS:doSimulation = 1 // == 1 if 1D simulated data, 0 if other from the checkbox 1298 root:Packages:NIST:SAS:gDoMonteCarlo = 0 // == 1 if 2D MonteCarlo set by hidden flag 1299 1300 1301 EndMacro 1302 1303 Function Sim_1D_CountTimeSetVarProc(sva) : SetVariableControl 1304 STRUCT WMSetVariableAction &sva 1305 1306 switch( sva.eventCode ) 1307 case 1: // mouse up 1308 case 2: // Enter key 1309 case 3: // Live update 1310 Variable dval = sva.dval 1311 1312 // get the neutron flux, multiply, and reset the global for # neutrons 1313 NVAR imon=root:Packages:NIST:SAS:gImon 1314 imon = dval*beamIntensity() 1315 1316 break 1317 endswitch 1318 1319 return 0 1320 End 1321 1322 1323 Function Sim_1D_ModelPopMenuProc(pa) : PopupMenuControl 1324 STRUCT WMPopupAction &pa 1325 1326 switch( pa.eventCode ) 1327 case 2: // mouse up 1328 Variable popNum = pa.popNum 1329 String popStr = pa.popStr 1330 SVAR gStr = root:Packages:NIST:SAS:gFuncStr 1331 gStr = popStr 1332 1333 break 1334 endswitch 1335 1336 return 0 1337 End 1338 1339 1340 1341 //Function Sim_1D_DoItButtonProc(ba) : ButtonControl 1342 // STRUCT WMButtonAction &ba 1343 // 1344 // switch( ba.eventCode ) 1345 // case 2: // mouse up 1346 // // click code here 1347 // NVAR doMC = root:Packages:NIST:SAS:gDoMonteCarlo 1348 // doMC = 1 1349 // ReCalculateInten(1) 1350 // doMC = 0 //so the next time won't be MC 1351 // break 1352 // endswitch 1353 // 1354 // return 0 1108 1355 //End 1109 1356 // 1110 //Function RunMonte(imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,funcStr,coefStr) 1111 // Variable imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh 1112 // String funcStr,coefStr 1113 // WAVE results 1114 // 1115 // FUNCREF SANSModelAAO_MCproto func=$funcStr 1116 // 1117 // Monte_SANS(imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,sig_sas,ran_dev,linear_data,results) 1357 // 1358 //Function Sim_1D_Display2DButtonProc(ba) : ButtonControl 1359 // STRUCT WMButtonAction &ba 1360 // 1361 // switch( ba.eventCode ) 1362 // case 2: // mouse up 1363 // // click code here 1364 // Execute "ChangeDisplay(\"SAS\")" 1365 // break 1366 // endswitch 1367 // 1368 // return 0 1118 1369 //End 1119 //1120 ////// END UNUSED BLOCK
Note: See TracChangeset
for help on using the changeset viewer.