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/ProtocolAsPanel.ipf

    r764 r794  
    242242        numItems = ItemsInList(list,";") 
    243243        checked = 1 
    244         if(numitems == 4) 
     244        if(numitems == 4 || numitems == 5)              //allow for protocols with no SDEV list item 
    245245                //correct number of parameters, assume ok 
    246246                String/G root:myGlobals:Protocols:gAbsStr = list 
     
    17401740        Endif 
    17411741         
    1742         Variable c2,c3,c4,c5 
     1742        Variable c2,c3,c4,c5,kappa_err 
    17431743        //do absolute scaling if desired 
    17441744        if(cmpstr("none",prot[4])!=0) 
     
    17521752                        c4 = NumberByKey("IZERO", junkAbsStr, "=", ";") 
    17531753                        c5 = NumberByKey("XSECT", junkAbsStr, "=", ";") 
     1754                        kappa_err = NumberByKey("SDEV", junkAbsStr, "=", ";") 
    17541755                else 
    17551756                        //get the parames from the list 
     
    17581759                        c4 = NumberByKey("IZERO", prot[4], "=", ";") 
    17591760                        c5 = NumberByKey("XSECT", prot[4], "=", ";") 
     1761                        kappa_err = NumberByKey("SDEV", prot[4], "=", ";") 
    17601762                Endif 
    17611763                //get the sample trans and thickness from the activeType folder 
     
    17651767                Variable c1 = dest[5]           //sample thickness 
    17661768                 
    1767                 err = Absolute_Scale(activeType,c0,c1,c2,c3,c4,c5) 
     1769                err = Absolute_Scale(activeType,c0,c1,c2,c3,c4,c5,kappa_err) 
    17681770                if(err) 
    17691771                        SetDataFolder root: 
     
    19901992//values are passed back as a global string variable (keyword=value) 
    19911993// 
    1992 Proc AskForAbsoluteParams(c2,c3,c4,c5) 
    1993         Variable c2=0.95,c3=0.1,c4=1,c5=32.0 
     1994Proc AskForAbsoluteParams(c2,c3,c4,c5,err) 
     1995        Variable c2=0.95,c3=0.1,c4=1,c5=32.0,err=0 
    19941996        Prompt c2, "Standard Transmission" 
    19951997        Prompt c3, "Standard Thickness (cm)" 
    19961998        Prompt c4, "I(0) from standard fit (normalized to 1E8 monitor cts)" 
    19971999        Prompt c5, "Standard Cross-Section (cm-1)" 
     2000        Prompt err, "error in I(q=0) (one std dev)" 
    19982001         
    19992002        String/G root:myGlobals:Protocols:gAbsStr="" 
     
    20032006        root:myGlobals:Protocols:gAbsStr +=  ";" + "IZERO="+num2str(c4) 
    20042007        root:myGlobals:Protocols:gAbsStr +=  ";" + "XSECT="+num2str(c5) 
     2008        root:myGlobals:Protocols:gAbsStr +=  ";" + "SDEV="+num2str(err) 
    20052009         
    20062010End 
     
    20592063                Variable lambda = rw[26] 
    20602064                Variable attenNo = rw[3] 
    2061                 attenTrans = AttenuationFactor(acctStr,lambda,attenNo) 
     2065                Variable atten_err 
     2066                attenTrans = AttenuationFactor(acctStr,lambda,attenNo,atten_err) 
    20622067                //Print "attenTrans = ",attenTrans 
    20632068                 
    20642069                //get the XY box, if needed 
    2065                 Variable x1,x2,y1,y2 
     2070                Variable x1,x2,y1,y2,ct_err 
    20662071                String filename=tw[0],tempStr 
    20672072                PathInfo/S catPathName 
     
    20922097                 
    20932098                //need the detector sensitivity file - make a guess, allow to override 
    2094                 String junkStr="" 
     2099                String junkStr="",errStr="" 
    20952100                if(! waveexists($"root:Packages:NIST:DIV:data")) 
    20962101                        junkStr = PromptForPath("Select the detector sensitivity file") 
     
    21242129//              detCnt = sum($"root:Packages:NIST:raw:data", -inf, inf ) 
    21252130//              Print "box is now ",x1,x2,y1,y2 
    2126                 detCnt = SumCountsInBox(x1,x2,y1,y2,"RAW") 
     2131                detCnt = SumCountsInBox(x1,x2,y1,y2,ct_err,"RAW") 
    21272132                if(cmpstr(tw[9],"ILL   ")==0) 
    21282133                        detCnt /= 4             // for cerca detector, header is right, sum(data) is 4x too large this is usually corrected in the Add step 
     
    21312136                //               
    21322137                kappa = detCnt/countTime/attenTrans*1.0e8/(monCnt/countTime)*(pixel/sdd)^2 
     2138                 
     2139                Variable kappa_err 
     2140                kappa_err = (ct_err/detCnt)^2 + (atten_err/attenTrans)^2 
     2141                kappa_err = sqrt(kappa_err) * kappa 
     2142                 
    21332143                junkStr = num2str(kappa) 
     2144                errStr = num2Str(kappa_err) 
    21342145                // set the parameters in the global string 
    2135                 Execute "AskForAbsoluteParams(1,1,"+junkStr+",1)"               //no missing parameters, no dialog 
     2146//              Print "AskForAbsoluteParams(1,1,"+junkStr+",1,"+errStr+")"       
     2147                Execute "AskForAbsoluteParams(1,1,"+junkStr+",1,"+errStr+")"            //no missing parameters, no dialog 
    21362148                 
    21372149                //should wipe out the data in the RAW folder, since it's not really RAW now 
     
    21402152                // - obsucre bug if "ask" in ABS section of protocol clears RAW folder, then Q-axes can't be set from RAW:RealsRead 
    21412153                //ClearDataFolder("RAW") 
    2142                 Print "Kappa was successfully calculated as = ",kappa    
     2154                Printf "Kappa was successfully calculated as = %g +/- %g (%g %)\r",kappa,kappa_err,(kappa_err/kappa)*100 
    21432155        Endif 
    21442156         
Note: See TracChangeset for help on using the changeset viewer.