Ignore:
Timestamp:
Jul 23, 2012 4:54:42 PM (11 years ago)
Author:
srkline
Message:

Switched the Decay panel plot to use atomic polariztion as the y-axis, rather than mu*P.

Added atomic polarization to the disaplayed matrix on the panel. It was already part of the calculation behind the scenes.

Switched the decay fit to use exp_XOffset (with an offset of 0). this fits naturally as gamma, not 1/gamma.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/Polarization/Pol_PolarizationPanels.ipf

    r858 r859  
    631631        SetDataFolder root:Packages:NIST:Polarization:Cells 
    632632 
    633         Make/O/D/N=(1,8) $("Decay_"+popStr) 
     633        Make/O/D/N=(1,9) $("Decay_"+popStr) 
    634634        WAVE decay = $("Decay_"+popStr) 
    635635        // set the column labels 
    … …  
    638638        SetDimLabel 1,2,Blocked,decay 
    639639        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, non-zero 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, non-zero is used in the fit 
     644        SetDimLabel 1,8,elapsed_hr,decay 
     645        decay[0][7] = 1                 //default to include the point 
    645646         
    646647        // generate the dummy wave note now, change as needed 
    … …  
    775776                                // 2 find the mu and Te values for cellStr 
    776777                                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 
    782779        //                       
    783780                                // 3 Calc muPo and error 
    … …  
    789786                                // 3.5 calc Polarization of cell (no value or error stored in calc wave?) 
    790787                                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 
    793789         
    794790                                // 4 calc Po and error 
    795791                                Po = Calc_Po(gCellKW,muPo,err_muPo,err_Po) 
    796         //                      Po = Calc_Po(gCellKW,2,err_muPo,err_Po) 
    797792                                calc[selRow][%Po] = Po 
    798793                                calc[selRow][%err_Po] = err_Po 
     794                                w[selRow][%Atomic_Pol] = Po             //for display 
    799795                                 
    800796                                // 5 calc Tmaj and error 
    … …  
    949945        tmp = (err_muPo/muPo)^2 + (err_mu/mu)^2 
    950946        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) 
    951949         
    952950        Printf "Po = %g +/- %g (%g%)\r",Po,err_Po,err_Po/Po*100 
    … …  
    999997 
    1000998 
    1001 //Function testCR(num) 
    1002 //      Variable num 
    1003 //      Variable err_cr 
    1004 //       
    1005 //      Variable noNorm=0 
    1006 //      Variable cr = TotalCR_FromRun(num,err_cr,noNorm) 
    1007 //      printf "CR = %g +/- %g (%g%)\r",cr,err_cr,err_cr/cr*100  
    1008 //      return(0) 
    1009 //End 
     999Function 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) 
     1007End 
    10101008 
    10111009// calculate the total detector CR and its error. 
    … …  
    10871085 
    10881086        String cellStr="" 
    1089         Variable num 
     1087        Variable num,plot_P 
     1088         
     1089        plot_P = 1              //plot_Pol as Y, or muP if this == 0 
    10901090 
    10911091        switch( ba.eventCode ) 
    … …  
    11031103//                      make temp copies for the fit and plot, extra for mask 
    11041104                        num = DimSize(calc,0) 
    1105                         Make/O/D/N=(num)                tmp_Mask,tmp_hr,tmp_muP,tmp_err_muP,tmp_muP2 
     1105                        Make/O/D/N=(num)                tmp_Mask,tmp_hr,tmp_Y,tmp_err_Y,tmp_Y2 
    11061106                         
    11071107                        tmp_Mask = decay[p][%Include] 
    11081108                        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 y-axis 
     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 
    11151125                        // clear old data, and plot the new 
    11161126                        // 
    1117                         CheckDisplayed/W=DecayPanel#G0 tmp_muP,tmp_muP2,fit_tmp_muP 
     1127                        CheckDisplayed/W=DecayPanel#G0 tmp_Y,tmp_Y2,fit_tmp_Y 
    11181128                        // if both present, bit 0 + bit 1 = 3 
    11191129                        if(V_flag & 2^0)                        //check bit 0 
    1120                                 RemoveFromGraph/W=DecayPanel#G0 tmp_muP 
     1130                                RemoveFromGraph/W=DecayPanel#G0 tmp_Y 
    11211131                        endif 
    11221132                        if(V_flag & 2^1) 
    1123                                 RemoveFromGraph/W=DecayPanel#G0 tmp_muP2 
     1133                                RemoveFromGraph/W=DecayPanel#G0 tmp_Y2 
    11241134                        endif 
    11251135                        if(V_flag & 2^2) 
    1126                                 RemoveFromGraph/W=DecayPanel#G0 fit_tmp_muP 
     1136                                RemoveFromGraph/W=DecayPanel#G0 fit_tmp_Y 
    11271137                        endif 
    11281138                         
    1129                         AppendToGraph/W=DecayPanel#G0 tmp_muP vs tmp_hr 
    1130                         AppendToGraph/W=DecayPanel#G0 tmp_muP2 vs tmp_hr 
     1139                        AppendToGraph/W=DecayPanel#G0 tmp_Y vs tmp_hr 
     1140                        AppendToGraph/W=DecayPanel#G0 tmp_Y2 vs tmp_hr 
    11311141 
    11321142                        ModifyGraph/W=DecayPanel#G0 log(left)=1 
    … …  
    11341144                        ModifyGraph/W=DecayPanel#G0 mode=3 
    11351145                        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) 
    11371147                        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                    
    11411155                        Label/W=DecayPanel#G0 bottom "time (h)" 
    11421156                         
    11431157// do the fit 
    11441158//       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 
    11451160 
    11461161                        SetActiveSubwindow DecayPanel#G0                        //to get the automatic fit to show up on the graph 
    11471162 
     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 
    11481169//                      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 
    11541171                         
    11551172 
    … …  
    11711188                        SVAR gMuPo = root:Packages:NIST:Polarization:Cells:gMuPo 
    11721189                        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 
    11781211 
    11791212                        // 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] 
    11821215 
    11831216                        // 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])^2 
     1217//                      gamma_val  = 1/fitCoef[2] 
     1218//                      err_gamma = W_sigma[2]/(fitCoef[2])^2 
    11861219                 
    11871220                         
    … …  
    12011234                         
    12021235                        // for the panel display 
    1203                         sprintf gMuPo, "%g +/- %g",fitCoef[1],W_sigma[1] 
     1236                        sprintf gMuPo, "%g +/- %g",muPo,err_muPo 
    12041237                        sprintf gPo, "%g +/- %g",Po,err_Po 
    12051238                        sprintf gGamma, "%g +/- %g",gamma_val,err_gamma 
Note: See TracChangeset for help on using the changeset viewer.