Changeset 859
 Timestamp:
 Jul 23, 2012 4:54:42 PM (10 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

sans/Dev/trunk/NCNR_User_Procedures/Reduction/Polarization/Pol_PolarizationPanels.ipf
r858 r859 631 631 SetDataFolder root:Packages:NIST:Polarization:Cells 632 632 633 Make/O/D/N=(1, 8) $("Decay_"+popStr)633 Make/O/D/N=(1,9) $("Decay_"+popStr) 634 634 WAVE decay = $("Decay_"+popStr) 635 635 // set the column labels … … 638 638 SetDimLabel 1,2,Blocked,decay 639 639 SetDimLabel 1,3,mu_star,decay 640 SetDimLabel 1,4,Pol_Cell,decay 641 SetDimLabel 1,5,T_Major,decay 642 SetDimLabel 1,6,Include,decay //for a mask wave, nonzero is used in the fit 643 SetDimLabel 1,7,elapsed_hr,decay 644 decay[0][6] = 1 //default to include the point 640 SetDimLabel 1,4,Effective_Pol,decay 641 SetDimLabel 1,5,Atomic_Pol,decay 642 SetDimLabel 1,6,T_Major,decay 643 SetDimLabel 1,7,Include,decay //for a mask wave, nonzero is used in the fit 644 SetDimLabel 1,8,elapsed_hr,decay 645 decay[0][7] = 1 //default to include the point 645 646 646 647 // generate the dummy wave note now, change as needed … … 775 776 // 2 find the mu and Te values for cellStr 776 777 SVAR gCellKW = $("root:Packages:NIST:Polarization:Cells:gCell_"+cellStr) 777 //(moved to a separate function, just pass the string) 778 // Te = NumberByKey("Te", gCellKW, "=", ",", 0) 779 // err_Te = NumberByKey("err_Te", gCellKW, "=", ",", 0) 780 // mu = NumberByKey("mu", gCellKW, "=", ",", 0) 781 // err_mu = NumberByKey("err_mu", gCellKW, "=", ",", 0) 778 782 779 // 783 780 // 3 Calc muPo and error … … 789 786 // 3.5 calc Polarization of cell (no value or error stored in calc wave?) 790 787 PCell = Calc_PCell(muPo,err_muPo,err_PCell) 791 // PCell = Calc_PCell(2,err_muPo,err_PCell) 792 w[selRow][%Pol_Cell] = PCell 788 w[selRow][%Effective_Pol] = PCell 793 789 794 790 // 4 calc Po and error 795 791 Po = Calc_Po(gCellKW,muPo,err_muPo,err_Po) 796 // Po = Calc_Po(gCellKW,2,err_muPo,err_Po)797 792 calc[selRow][%Po] = Po 798 793 calc[selRow][%err_Po] = err_Po 794 w[selRow][%Atomic_Pol] = Po //for display 799 795 800 796 // 5 calc Tmaj and error … … 949 945 tmp = (err_muPo/muPo)^2 + (err_mu/mu)^2 950 946 err_Po = Po * sqrt(tmp) 947 // tmp = 1/mu^2*err_muPo^2 + muPo^2/mu^4*err_mu^2 948 // err_Po = sqrt(tmp) 951 949 952 950 Printf "Po = %g +/ %g (%g%)\r",Po,err_Po,err_Po/Po*100 … … 999 997 1000 998 1001 //Function testCR(num)1002 //Variable num1003 //Variable err_cr1004 //1005 //Variable noNorm=01006 //Variable cr = TotalCR_FromRun(num,err_cr,noNorm)1007 //printf "CR = %g +/ %g (%g%)\r",cr,err_cr,err_cr/cr*1001008 //return(0)1009 //End999 Function testCR(num) 1000 Variable num 1001 Variable err_cr 1002 1003 Variable noNorm=0 1004 Variable cr = TotalCR_FromRun(num,err_cr,noNorm) 1005 printf "CR = %g +/ %g (%g%)\r",cr,err_cr,err_cr/cr*100 1006 return(0) 1007 End 1010 1008 1011 1009 // calculate the total detector CR and its error. … … 1087 1085 1088 1086 String cellStr="" 1089 Variable num 1087 Variable num,plot_P 1088 1089 plot_P = 1 //plot_Pol as Y, or muP if this == 0 1090 1090 1091 1091 switch( ba.eventCode ) … … 1103 1103 // make temp copies for the fit and plot, extra for mask 1104 1104 num = DimSize(calc,0) 1105 Make/O/D/N=(num) tmp_Mask,tmp_hr,tmp_ muP,tmp_err_muP,tmp_muP21105 Make/O/D/N=(num) tmp_Mask,tmp_hr,tmp_Y,tmp_err_Y,tmp_Y2 1106 1106 1107 1107 tmp_Mask = decay[p][%Include] 1108 1108 tmp_hr = decay[p][%elapsed_hr] 1109 tmp_muP = calc[p][%muPo] 1110 tmp_muP2 = tmp_muP 1111 tmp_err_muP = calc[p][%err_muPo] 1112 1113 tmp_muP2 = (tmp_Mask == 1) ? NaN : tmp_muP2 //only excluded points will plot 1114 1109 1110 if(plot_P == 1) 1111 tmp_Y = calc[p][%Po] 1112 tmp_Y2 = tmp_Y 1113 tmp_err_Y = calc[p][%err_Po] 1114 else //plot muP as the yaxis 1115 tmp_Y = calc[p][%muPo] 1116 tmp_Y2 = tmp_Y 1117 tmp_err_Y = calc[p][%err_muPo] 1118 endif 1119 1120 tmp_Y2 = (tmp_Mask == 1) ? NaN : tmp_Y2 //only excluded points will plot 1121 1122 1123 1124 1115 1125 // clear old data, and plot the new 1116 1126 // 1117 CheckDisplayed/W=DecayPanel#G0 tmp_ muP,tmp_muP2,fit_tmp_muP1127 CheckDisplayed/W=DecayPanel#G0 tmp_Y,tmp_Y2,fit_tmp_Y 1118 1128 // if both present, bit 0 + bit 1 = 3 1119 1129 if(V_flag & 2^0) //check bit 0 1120 RemoveFromGraph/W=DecayPanel#G0 tmp_ muP1130 RemoveFromGraph/W=DecayPanel#G0 tmp_Y 1121 1131 endif 1122 1132 if(V_flag & 2^1) 1123 RemoveFromGraph/W=DecayPanel#G0 tmp_ muP21133 RemoveFromGraph/W=DecayPanel#G0 tmp_Y2 1124 1134 endif 1125 1135 if(V_flag & 2^2) 1126 RemoveFromGraph/W=DecayPanel#G0 fit_tmp_ muP1136 RemoveFromGraph/W=DecayPanel#G0 fit_tmp_Y 1127 1137 endif 1128 1138 1129 AppendToGraph/W=DecayPanel#G0 tmp_ muPvs tmp_hr1130 AppendToGraph/W=DecayPanel#G0 tmp_ muP2 vs tmp_hr1139 AppendToGraph/W=DecayPanel#G0 tmp_Y vs tmp_hr 1140 AppendToGraph/W=DecayPanel#G0 tmp_Y2 vs tmp_hr 1131 1141 1132 1142 ModifyGraph/W=DecayPanel#G0 log(left)=1 … … 1134 1144 ModifyGraph/W=DecayPanel#G0 mode=3 1135 1145 ModifyGraph/W=DecayPanel#G0 marker=19 1136 ModifyGraph/W=DecayPanel#G0 rgb(tmp_ muP)=(1,16019,65535),rgb(tmp_muP2)=(65535,0,0)1146 ModifyGraph/W=DecayPanel#G0 rgb(tmp_Y)=(1,16019,65535),rgb(tmp_Y2)=(65535,0,0) 1137 1147 ModifyGraph/W=DecayPanel#G0 msize=3 1138 ErrorBars/W=DecayPanel#G0 tmp_muP,Y wave=(tmp_err_muP,tmp_err_muP) 1139 1140 Label/W=DecayPanel#G0 left "mu*P" 1148 ErrorBars/W=DecayPanel#G0 tmp_Y,Y wave=(tmp_err_Y,tmp_err_Y) 1149 1150 if(plot_P == 1) 1151 Label/W=DecayPanel#G0 left "Atomic Polarization, P" 1152 else 1153 Label/W=DecayPanel#G0 left "mu*P" 1154 endif 1141 1155 Label/W=DecayPanel#G0 bottom "time (h)" 1142 1156 1143 1157 // do the fit 1144 1158 // as long as the constant X0 doesn't stray from zero, exp_XOffset is OK, otherwise I'll need to switch to the exp function 1159 //  use the /K={0} flag to set the constant to zero 1145 1160 1146 1161 SetActiveSubwindow DecayPanel#G0 //to get the automatic fit to show up on the graph 1147 1162 1163 Make/O/D/N=3 fitCoef={0,5,0.05} 1164 CurveFit/H="100"/M=2/W=0/TBOX=(0x310)/K={0} exp_XOffset, kwCWave=fitCoef, tmp_Y /X=tmp_hr /D /I=1 /W=tmp_err_Y /M=tmp_Mask 1165 // gives nice error on gamma, but it's wrong since it doesn't use the errors on muP 1166 // CurveFit/H="100"/M=2/W=0/TBOX=(0x310)/K={0} exp_XOffset, kwCWave=fitCoef, tmp_Y /X=tmp_hr /D /M=tmp_Mask 1167 1168 1148 1169 // Make/O/D/N=3 fitCoef={0,5,0.05} 1149 // CurveFit/H="100"/M=2/W=0/TBOX=(0x310) exp_XOffset, kwCWave=fitCoef, tmp_muP /X=tmp_hr /D /I=1 /W=tmp_err_muP /M=tmp_Mask 1150 1151 1152 Make/O/D/N=3 fitCoef={0,5,0.05} 1153 CurveFit/H="100"/M=2/W=0/TBOX=(0x310) exp, kwCWave=fitCoef, tmp_muP /X=tmp_hr /D /I=1 /W=tmp_err_muP /M=tmp_Mask 1170 // CurveFit/H="100"/M=2/W=0/TBOX=(0x310) exp, kwCWave=fitCoef, tmp_Y /X=tmp_hr /D /I=1 /W=tmp_err_Y /M=tmp_Mask 1154 1171 1155 1172 … … 1171 1188 SVAR gMuPo = root:Packages:NIST:Polarization:Cells:gMuPo 1172 1189 SVAR gPo = root:Packages:NIST:Polarization:Cells:gPo 1173 1174 muPo = fitCoef[1] 1175 err_muPo = W_sigma[1] 1176 1177 Po = Calc_Po(gCellKW,muPo,err_muPo,err_Po) 1190 1191 Variable mu,err_mu 1192 1193 if(plot_P == 1) 1194 Po = fitCoef[1] 1195 err_Po = W_sigma[1] 1196 // calc muPo inline here, since the Calc_muP function uses a different method 1197 // cell constants 1198 mu = NumberByKey("mu", gCellKW, "=", ",", 0) 1199 err_mu = NumberByKey("err_mu", gCellKW, "=", ",", 0) 1200 1201 muPo = Po*mu 1202 err_muPo = muPo*sqrt( (err_Po/Po)^2 + (err_mu/mu)^2 ) 1203 else 1204 muPo = fitCoef[1] 1205 err_muPo = W_sigma[1] 1206 1207 Po = Calc_Po(gCellKW,muPo,err_muPo,err_Po) 1208 endif 1209 1210 1178 1211 1179 1212 // if exp_XOffset used 1180 // gGamma= fitCoef[2]1181 //err_Gamma = W_sigma[2]1213 gamma_val = fitCoef[2] 1214 err_Gamma = W_sigma[2] 1182 1215 1183 1216 // calculating the error using exp is the inverse of coef[2]: 1184 gamma_val = 1/fitCoef[2]1185 err_gamma = W_sigma[2]/(fitCoef[2])^21217 // gamma_val = 1/fitCoef[2] 1218 // err_gamma = W_sigma[2]/(fitCoef[2])^2 1186 1219 1187 1220 … … 1201 1234 1202 1235 // for the panel display 1203 sprintf gMuPo, "%g +/ %g", fitCoef[1],W_sigma[1]1236 sprintf gMuPo, "%g +/ %g",muPo,err_muPo 1204 1237 sprintf gPo, "%g +/ %g",Po,err_Po 1205 1238 sprintf gGamma, "%g +/ %g",gamma_val,err_gamma
Note: See TracChangeset
for help on using the changeset viewer.