Ignore:
Timestamp:
Apr 4, 2011 11:12:38 AM (11 years ago)
Author:
srkline
Message:

Lots of changes:
-2D resolution smearing
-2D error propagation

1) 2D resolution smearing has been corrected to use sigma (perp) correctly
rather than phi. This shows up in the quadrature loop in all of the 2D models
and in the Smear_2D "generic" function.

2) 1D resolutionn smearing is now corrected to account for the integration
range of +/- 3 sigma (99.73% of distribution). Now the result is divided by
0.9973 to rescale it to the correct value.

3) Smeared models are now AAO to improve speed and to allow easier use with
functions that are inherently AAO. No changes are needed, since the call is
behind the scenes, replacing Smear_Model_N() with Smear_Model_N_AAO().

4) in PlotUtils_2D added functions to re-bin the QxQy? data into a 1D format
BinQxQy_to_1D(). This also re-bins the errors in two ways, adding the per-pixel
errors in quadrature, or the deviation from the mean of the intensity. Manually
editing the intensity declaration allows 2D->1D binning of smeared models.

5) Per-pixel propagation of errors has been added through the entire data
reduction sequence. Data errors are generated on loading using Poisson
statistics (specifically tailored for accuracy at low counts) and then is
propagated through each manipulation of the data using standard error
propagation. The error matrix is linear_data_error. As a by-product, all of
the math operations on data are explicitly done on linear_data, to avoid
any potential mistakes of log/linear scaling. Results of this propagation
largely match J. Barker's /ERROR routines from the VAX, with some differences
at low data count values (as expected) and at higher count values near the
beam stop (somewhat unexpected). This per-pixel error is ouput in the QxQy_ASCII
data files. NO CHANGE has been made to the 1D data, which uses the deviation from
the mean as the error - since this is correct.

6) Added tables for the uncertainty in attenuator transmission (from JGB)

7) added two more REAL values to the VAX header to store information
necessary for error propagation. These are couting error that are part of
the transmission error and of the absolute scaling error. These add Read/Write?
functions in NCNR_DataReadWrite

The transmission error (one standard deviation) is written at byte 396 (4 bytes)

RealsRead?[41]

The Box Sum error (one standard deviation) is written at byte 400 (4 bytes)

RealsRead?[42]

File:
1 edited

Legend:

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

    r788 r794  
    151151        SigmaQ = QBar*sqrt(v_r/rmd^2 +v_lambda) 
    152152 
     153 
     154// more readable method for calculating the variance in Q 
     155// EXCEPT - this is calculated for Qo, NOT qBar 
     156// (otherwise, they are nearly equivalent, except for close to the beam stop) 
     157//      Variable kap,a_val,a_val_l2,m_h 
     158//      g = 981.0                               //gravity acceleration [cm/s^2] 
     159//      m_h     = 252.8                 // m/h [=] s/cm^2 
     160//       
     161//      kap = 2*pi/lambda 
     162//      a_val = L2*(L1+L2)*g/2*(m_h)^2 
     163//      a_val_L2 = a_val/L2*1e-16               //convert 1/cm^2 to 1/A^2 
     164// 
     165//      sigmaQ = 0 
     166//      sigmaQ = 3*(S1/L1)^2 
     167//       
     168//      if(usingLenses != 0) 
     169//              sigmaQ += 2*(S2/lp)^2*(lambdaWidth)^2   //2nd term w/ lenses 
     170//      else     
     171//              sigmaQ += 2*(S2/lp)^2                                           //2nd term w/ no lenses 
     172//      endif 
     173//       
     174//      sigmaQ += (del_r/L2)^2 
     175//      sigmaQ += 2*(r0/L2)^2*(lambdaWidth)^2 
     176//      sigmaQ += 4*(a_val_l2)^2*lambda^4*(lambdaWidth)^2 
     177//       
     178//      sigmaQ *= kap^2/12 
     179//      sigmaQ = sqrt(sigmaQ) 
     180 
     181 
    153182        results = "success" 
    154183        Return results 
     
    179208// phi is the azimuthal angle, CCW from +x axis 
    180209// r_dist is the real-space distance from ctr of detector to QxQy pixel location 
     210// 
     211// MAR 2011 - removed the del_r terms, they don't apply since no bining is done to the 2D data 
    181212// 
    182213Function/S get2DResolution(inQ,phi,lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,r_dist,SigmaQX,SigmaQY,fSubS) 
     
    255286//      endif 
    256287 
    257         Variable kap,a_val,a_val_L2 
     288        Variable kap,a_val,a_val_L2,proj_DDet 
    258289         
    259290        kap = 2*pi/lambda 
     
    270301// 
    271302//      a_val = 0 
     303//      a_val_l2 = 0 
    272304 
    273305        // the detector pixel is square, so correct for phi 
    274         DDet = DDet*cos(phi) + DDet*sin(phi) 
     306        proj_DDet = DDet*cos(phi) + DDet*sin(phi) 
    275307         
    276308        // this is really sigma_Q_parallel 
    277         SigmaQX = kap*kap/12* (3*(S1/L1)^2 + 3*(S2/LP)^2 + (DDet/L2)^2 + (del_r/L2)^2 + (sin(phi))^2*8*(a_val_L2)^2*lambda^4*lambdaWidth^2) 
     309        SigmaQX = kap*kap/12* (3*(S1/L1)^2 + 3*(S2/LP)^2 + (proj_DDet/L2)^2 + (sin(phi))^2*8*(a_val_L2)^2*lambda^4*lambdaWidth^2) 
    278310        SigmaQX += inQ*inQ*v_lambda 
    279311 
    280312        //this is really sigma_Q_perpendicular 
    281         SigmaQY = 3*(S1/L1)^2 + 3*(S2/LP)^2 + (DDet/L2)^2 + (del_r/L2)^2 + (cos(phi))^2*8*(a_val_L2)^2*lambda^4*lambdaWidth^2 
     313        proj_DDet = DDet*sin(phi) + DDet*cos(phi)               //not necessary, since DDet is the same in both X and Y directions 
     314 
     315        SigmaQY = 3*(S1/L1)^2 + 3*(S2/LP)^2 + (proj_DDet/L2)^2 + (cos(phi))^2*8*(a_val_L2)^2*lambda^4*lambdaWidth^2 
    282316        SigmaQY *= kap*kap/12 
    283317         
     
    11941228        Make/O/N=(num) root:myGlobals:Attenuators:ng3att10 
    11951229         
     1230        // and a wave for the errors at each attenuation factor 
     1231        Make/O/N=(num) root:myGlobals:Attenuators:ng3att0_err 
     1232        Make/O/N=(num) root:myGlobals:Attenuators:ng3att1_err 
     1233        Make/O/N=(num) root:myGlobals:Attenuators:ng3att2_err 
     1234        Make/O/N=(num) root:myGlobals:Attenuators:ng3att3_err 
     1235        Make/O/N=(num) root:myGlobals:Attenuators:ng3att4_err 
     1236        Make/O/N=(num) root:myGlobals:Attenuators:ng3att5_err 
     1237        Make/O/N=(num) root:myGlobals:Attenuators:ng3att6_err 
     1238        Make/O/N=(num) root:myGlobals:Attenuators:ng3att7_err 
     1239        Make/O/N=(num) root:myGlobals:Attenuators:ng3att8_err 
     1240        Make/O/N=(num) root:myGlobals:Attenuators:ng3att9_err 
     1241        Make/O/N=(num) root:myGlobals:Attenuators:ng3att10_err 
     1242         
     1243         
    11961244        //each wave has 10 elements, the transmission of att# at the wavelengths  
    11971245        //lambda = 4,5,6,7,8,10,12,14,17,20 (4 A and 20 A are extrapolated values) 
     
    12101258        root:myGlobals:Attenuators:ng3att9 = {6.27682e-05,3.69e-05,1.908e-05,1.196e-05,8.738e-06,6.996e-06,6.2901e-07,2.60221e-07,7.49748e-08,2.08029e-08} 
    12111259        root:myGlobals:Attenuators:ng3att10 = {1.40323e-05,8.51e-06,5.161e-06,4.4e-06,4.273e-06,1.88799e-07,5.87021e-08,2.08169e-08,4.8744e-09,1.08687e-09} 
     1260   
     1261  // percent errors as measured, May 2007 values 
     1262  // zero error for zero attenuators, appropriate average values put in for unknown values 
     1263        root:myGlobals:Attenuators:ng3att0_err = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } 
     1264        root:myGlobals:Attenuators:ng3att1_err = {0.15,0.142,0.154,0.183,0.221,0.328,0.136,0.13,0.163,0.15} 
     1265        root:myGlobals:Attenuators:ng3att2_err = {0.25,0.257,0.285,0.223,0.271,0.405,0.212,0.223,0.227,0.25} 
     1266        root:myGlobals:Attenuators:ng3att3_err = {0.3,0.295,0.329,0.263,0.323,0.495,0.307,0.28,0.277,0.3} 
     1267        root:myGlobals:Attenuators:ng3att4_err = {0.35,0.331,0.374,0.303,0.379,0.598,0.367,0.322,0.33,0.35} 
     1268        root:myGlobals:Attenuators:ng3att5_err = {0.4,0.365,0.418,0.355,0.454,0.745,0.411,0.367,0.485,0.4} 
     1269        root:myGlobals:Attenuators:ng3att6_err = {0.45,0.406,0.473,0.385,0.498,0.838,0.454,0.49,0.5,0.5} 
     1270        root:myGlobals:Attenuators:ng3att7_err = {0.6,0.554,0.692,0.425,0.562,0.991,0.715,0.8,0.8,0.8} 
     1271        root:myGlobals:Attenuators:ng3att8_err = {0.7,0.705,0.927,0.503,0.691,1.27,1,1,1,1} 
     1272        root:myGlobals:Attenuators:ng3att9_err = {1,0.862,1.172,0.799,1.104,1.891,1.5,1.5,1.5,1.5} 
     1273        root:myGlobals:Attenuators:ng3att10_err = {1.5,1.054,1.435,1.354,1.742,2,2,2,2,2} 
     1274   
    12121275   
    12131276  //old tables, pre-June 2007 
     
    12461309        Make/O/N=(num) root:myGlobals:Attenuators:ng7att10 
    12471310         
     1311        // and a wave for the errors at each attenuation factor 
     1312        Make/O/N=(num) root:myGlobals:Attenuators:ng7att0_err 
     1313        Make/O/N=(num) root:myGlobals:Attenuators:ng7att1_err 
     1314        Make/O/N=(num) root:myGlobals:Attenuators:ng7att2_err 
     1315        Make/O/N=(num) root:myGlobals:Attenuators:ng7att3_err 
     1316        Make/O/N=(num) root:myGlobals:Attenuators:ng7att4_err 
     1317        Make/O/N=(num) root:myGlobals:Attenuators:ng7att5_err 
     1318        Make/O/N=(num) root:myGlobals:Attenuators:ng7att6_err 
     1319        Make/O/N=(num) root:myGlobals:Attenuators:ng7att7_err 
     1320        Make/O/N=(num) root:myGlobals:Attenuators:ng7att8_err 
     1321        Make/O/N=(num) root:myGlobals:Attenuators:ng7att9_err 
     1322        Make/O/N=(num) root:myGlobals:Attenuators:ng7att10_err   
     1323         
    12481324        //NG7 wave has 10 elements, the transmission of att# at the wavelengths  
    12491325        //lambda =4, 5,6,7,8,10,12,14,17,20 
     
    12651341        root:myGlobals:Attenuators:ng7att10 = {2.1607e-05,7.521e-06,2.91221e-06,1.45252e-06,7.93451e-07,1.92309e-07,5.99279e-08,2.87928e-08,2.19862e-08,1.7559e-08} 
    12661342 
     1343  // percent errors as measured, May 2007 values 
     1344  // zero error for zero attenuators, appropriate average values put in for unknown values 
     1345        root:myGlobals:Attenuators:ng7att0_err = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } 
     1346        root:myGlobals:Attenuators:ng7att1_err = {0.2,0.169,0.1932,0.253,0.298,0.4871,0.238,0.245,0.332,0.3} 
     1347        root:myGlobals:Attenuators:ng7att2_err = {0.3,0.305,0.3551,0.306,0.37,0.6113,0.368,0.413,0.45,0.4} 
     1348        root:myGlobals:Attenuators:ng7att3_err = {0.4,0.355,0.4158,0.36,0.4461,0.7643,0.532,0.514,0.535,0.5} 
     1349        root:myGlobals:Attenuators:ng7att4_err = {0.45,0.402,0.4767,0.415,0.5292,0.9304,0.635,0.588,0.623,0.6} 
     1350        root:myGlobals:Attenuators:ng7att5_err = {0.5,0.447,0.5376,0.487,0.6391,1.169,0.708,0.665,0.851,0.8} 
     1351        root:myGlobals:Attenuators:ng7att6_err = {0.6,0.501,0.6136,0.528,0.8796,1.708,0.782,0.874,1,1} 
     1352        root:myGlobals:Attenuators:ng7att7_err = {0.8,0.697,0.9149,0.583,1.173,2.427,1.242,2,2,2} 
     1353        root:myGlobals:Attenuators:ng7att8_err = {1,0.898,1.24,0.696,1.577,3.412,3,3,3,3} 
     1354        root:myGlobals:Attenuators:ng7att9_err = {1.5,1.113,1.599,1.154,2.324,4.721,5,5,5,5} 
     1355        root:myGlobals:Attenuators:ng7att10_err = {1.5,1.493,5,5,5,5,5,5,5,5}   
     1356   
     1357 
    12671358// Pre-June 2007 calibration values - do not use these anymore   
    12681359////    root:myGlobals:Attenuators:ng7att0 = {1, 1, 1, 1, 1, 1, 1, 1 ,1} 
     
    12851376// 
    12861377// Mar 2010 - abs() added to attStr to account for ICE reporting -0.0001 as an attenuator position, which truncates to "-0" 
    1287 Function LookupAttenNG3(lambda,attenNo) 
    1288         Variable lambda, attenNo 
     1378Function LookupAttenNG3(lambda,attenNo,atten_err) 
     1379        Variable lambda, attenNo, &atten_err 
    12891380         
    12901381        Variable trans 
    12911382        String attStr="root:myGlobals:Attenuators:ng3att"+num2str(trunc(abs(attenNo))) 
     1383        String attErrWStr="root:myGlobals:Attenuators:ng3att"+num2str(trunc(abs(attenNo)))+"_err" 
    12921384        String lamStr = "root:myGlobals:Attenuators:ng3lambda" 
    12931385         
     
    13001392        Endif 
    13011393         
    1302         if(!(WaveExists($attStr)) || !(WaveExists($lamStr)) ) 
     1394        if(!(WaveExists($attStr)) || !(WaveExists($lamStr)) || !(WaveExists($attErrWStr))) 
    13031395                Execute "MakeNG3AttenTable()" 
    13041396        Endif 
     
    13111403        //the attenuator must always be an integer 
    13121404        Wave att = $attStr 
     1405        Wave attErrW = $attErrWStr 
    13131406        Wave lam = $lamstr 
    13141407        trans = interp(lambda,lam,att) 
    1315          
     1408        atten_err = interp(lambda,lam,attErrW) 
     1409 
     1410// the error in the tables is % error. return the standard deviation instead 
     1411        atten_err = trans*atten_err/100 
     1412                 
    13161413//      Print "trans = ",trans 
     1414//      Print "trans err = ",atten_err 
    13171415         
    13181416        return trans 
     
    13291427// 
    13301428// Mar 2010 - abs() added to attStr to account for ICE reporting -0.0001 as an attenuator position, which truncates to "-0" 
    1331 Function LookupAttenNG7(lambda,attenNo) 
    1332         Variable lambda, attenNo 
     1429Function LookupAttenNG7(lambda,attenNo,atten_err) 
     1430        Variable lambda, attenNo, &atten_err 
    13331431         
    13341432        Variable trans 
    13351433        String attStr="root:myGlobals:Attenuators:ng7att"+num2str(trunc(abs(attenNo))) 
     1434        String attErrWStr="root:myGlobals:Attenuators:ng7att"+num2str(trunc(abs(attenNo)))+"_err" 
    13361435        String lamStr = "root:myGlobals:Attenuators:ng7lambda" 
    13371436         
     
    13441443        Endif 
    13451444         
    1346         if(!(WaveExists($attStr)) || !(WaveExists($lamStr)) ) 
     1445        if(!(WaveExists($attStr)) || !(WaveExists($lamStr)) || !(WaveExists($attErrWStr))) 
    13471446                Execute "MakeNG7AttenTable()" 
    13481447        Endif 
     
    13551454        //the attenuator must always be an integer 
    13561455        Wave att = $attStr 
     1456        Wave attErrW = $attErrWStr 
    13571457        Wave lam = $lamstr 
    13581458        trans = interp(lambda,lam,att) 
    1359          
    1360         //Print "trans = ",trans 
     1459        atten_err = interp(lambda,lam,attErrW) 
     1460 
     1461// the error in the tables is % error. return the standard deviation instead 
     1462        atten_err = trans*atten_err/100 
     1463         
     1464//      Print "trans = ",trans 
     1465//      Print "trans err = ",atten_err 
    13611466         
    13621467        return trans 
     
    13791484// called by Correct.ipf, ProtocolAsPanel.ipf, Transmission.ipf 
    13801485// 
    1381 Function AttenuationFactor(fileStr,lam,attenNo) 
     1486// 
     1487// as of March 2011, returns the error (one standard deviation) in the attenuation factor as the last parameter, by reference 
     1488Function AttenuationFactor(fileStr,lam,attenNo,atten_err) 
    13821489        String fileStr 
    1383         Variable lam,attenNo 
     1490        Variable lam,attenNo, &atten_err 
    13841491         
    13851492        Variable attenFactor=1,loc 
     
    13881495        strswitch(instr) 
    13891496                case "NG3": 
    1390                         attenFactor = LookupAttenNG3(lam,attenNo) 
     1497                        attenFactor = LookupAttenNG3(lam,attenNo,atten_err) 
    13911498                        break 
    13921499                case "NG5": 
    13931500                        //using NG7 lookup Table 
    1394                         attenFactor = LookupAttenNG7(lam,attenNo) 
     1501                        attenFactor = LookupAttenNG7(lam,attenNo,atten_err) 
    13951502                        break 
    13961503                case "NG7": 
    1397                         attenFactor = LookupAttenNG7(lam,attenNo) 
     1504                        attenFactor = LookupAttenNG7(lam,attenNo,atten_err) 
    13981505                        break 
    13991506                default:                                                         
Note: See TracChangeset for help on using the changeset viewer.