Changeset 535


Ignore:
Timestamp:
Jul 7, 2009 6:35:07 PM (13 years ago)
Author:
srkline
Message:

changes to move 2D MonteCarlo? to a hidden feature and include 1D simulation. Skeleton of functionality only. Math still
needs to be done as well as interface.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/MultScatter_MonteCarlo_2D.ipf

    r494 r535  
    756756         
    757757        // 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;" 
    759759        list = RemoveFromList(tmp, list  ,";") 
    760760 
     
    10871087 
    10881088 
    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() 
     1090Function        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 well-behaved.",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 non-threaded 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)*1e-6 
     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 "non-reentrant". 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)*1e-6 
     1202        Printf  "MC sim time = %g seconds\r",t0 
     1203         
     1204        trans = results[8]                      //(n1-n2)/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((p-xCtr)^2+(q-yCtr)^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        // re-average 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) 
     1253end 
     1254 
     1255 
     1256 
     1257// as a stand-alone panel, extra control bar  (right) and subwindow implementations don't work right  
     1258// for various reasons... 
     1259Window 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         
     1301EndMacro 
     1302 
     1303Function 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 
     1320End 
     1321 
     1322 
     1323Function 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 
     1337End 
     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 
    11081355//End 
    11091356// 
    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 
    11181369//End 
    1119 // 
    1120 ////// END UNUSED BLOCK 
Note: See TracChangeset for help on using the changeset viewer.