Ignore:
Timestamp:
Apr 4, 2011 11:12:38 AM (12 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/Correct.ipf

    r665 r794  
    103103        //copy SAM information to COR, wiping out the old contents of the COR folder first 
    104104        //do this even if no correction is dispatched (if incorrect mode) 
     105        // the Copy converts "SAM" to linear scale, so "COR" is then linear too 
    105106        err = CopyWorkContents("SAM","COR")      
    106107        if(err==1) 
     
    234235// data exists, checked by dispatch routine 
    235236// 
     237// this is the most common use 
     238// March 2011 added error propagation 
     239//                                      added explicit reference to use linear_data, instead of trusting that data 
     240//                                      was freshly loaded. added final copy of cor result to cor:data and cor:linear_data 
     241// 
    236242Function CorrectMode_1() 
    237243         
    238244        //create the necessary wave references 
    239         WAVE sam_data=$"root:Packages:NIST:SAM:data" 
     245        WAVE sam_data=$"root:Packages:NIST:SAM:linear_data" 
    240246        WAVE sam_reals=$"root:Packages:NIST:SAM:realsread" 
    241247        WAVE sam_ints=$"root:Packages:NIST:SAM:integersread" 
    242248        WAVE/T sam_text=$"root:Packages:NIST:SAM:textread" 
    243         WAVE bgd_data=$"root:Packages:NIST:BGD:data" 
     249        WAVE bgd_data=$"root:Packages:NIST:BGD:linear_data" 
    244250        WAVE bgd_reals=$"root:Packages:NIST:BGD:realsread" 
    245251        WAVE bgd_ints=$"root:Packages:NIST:BGD:integersread" 
    246252        WAVE/T bgd_text=$"root:Packages:NIST:BGD:textread" 
    247         WAVE emp_data=$"root:Packages:NIST:EMP:data" 
     253        WAVE emp_data=$"root:Packages:NIST:EMP:linear_data" 
    248254        WAVE emp_reals=$"root:Packages:NIST:EMP:realsread" 
    249255        WAVE emp_ints=$"root:Packages:NIST:EMP:integersread" 
    250256        WAVE/T emp_text=$"root:Packages:NIST:EMP:textread" 
    251         WAVE cor_data=$"root:Packages:NIST:COR:data" 
     257        WAVE cor_data=$"root:Packages:NIST:COR:linear_data" 
    252258        WAVE/T cor_text=$"root:Packages:NIST:COR:textread" 
     259         
     260        // needed to propagate error 
     261        WAVE cor_data_display=$"root:Packages:NIST:COR:data"            //just for the final copy 
     262        WAVE sam_err =$"root:Packages:NIST:SAM:linear_data_error" 
     263        WAVE bgd_err =$"root:Packages:NIST:BGD:linear_data_error" 
     264        WAVE emp_err =$"root:Packages:NIST:EMP:linear_data_error" 
     265        WAVE cor_err =$"root:Packages:NIST:COR:linear_data_error" 
     266         
     267        Variable sam_trans_err,emp_trans_err 
     268        sam_trans_err = sam_reals[41] 
     269        emp_trans_err = emp_reals[41] 
     270         
    253271         
    254272        //get sam and bgd attenuation factors 
     
    257275        Variable tmonsam,fsam,fbgd,xshift,yshift,rsam,csam,rbgd,cbgd,tmonbgd 
    258276        Variable wcen=0.001,tsam,temp,remp,cemp,tmonemp,femp 
     277        Variable sam_atten_err,emp_atten_err,bgd_atten_err 
    259278        fileStr = sam_text[3] 
    260279        lambda = sam_reals[26] 
    261280        attenNo = sam_reals[3] 
    262         sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo) 
     281        sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,sam_atten_err) 
    263282        fileStr = bgd_text[3] 
    264283        lambda = bgd_reals[26] 
    265284        attenNo = bgd_reals[3] 
    266         bgd_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo) 
     285        bgd_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,bgd_atten_err) 
    267286        fileStr = emp_text[3] 
    268287        lambda = emp_reals[26] 
    269288        attenNo = emp_reals[3] 
    270         emp_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo) 
     289        emp_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,emp_atten_err) 
    271290         
    272291        //get relative monitor counts (should all be 10^8, since normalized in add step) 
     
    320339        cor1 -= (tsam/temp)*(femp*emp_temp/emp_attenFactor - fbgd*bgd_temp/bgd_attenFactor) 
    321340        cor1 *= noadd_bgd*noadd_emp             //zero out the array mismatch values 
     341 
     342// do the error propagation piecewise    
     343        Duplicate/O sam_err, tmp_a, tmp_b, tmp_c, tmp_d,c_val,d_val 
     344        tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2              //sig a ^2 
     345         
     346        tmp_b = (bgd_err/bgd_attenFactor)^2*(tsam/temp - 1)^2 + (bgd_atten_err*bgd_data/bgd_attenFactor^2)^2*(1-tsam/temp)^2            //sig b ^2 
     347 
     348        tmp_c = (sam_trans_err/temp)^2*(emp_data/emp_attenFactor-bgd_data/bgd_attenFactor)^2 
     349        tmp_c += (tsam/temp^2)^2*emp_trans_err^2*(emp_data/emp_attenFactor-bgd_data/bgd_attenFactor)^2 
     350         
     351        tmp_d = (tsam/(temp*emp_attenFactor))^2*(emp_err)^2 + (tsam*emp_data/(temp*emp_attenFactor^2))^2*(emp_atten_err)^2 
     352 
     353        cor_err = sqrt(tmp_a + tmp_b + tmp_c + tmp_d) 
    322354         
    323355        //we're done, get out w/no error 
    324         //set the COR data to the result 
     356        //set the COR data and linear_data to the result 
    325357        cor_data = cor1 
     358        cor_data_display = cor1 
     359         
    326360        //update COR header 
    327361        cor_text[1] = date() + " " + time()             //date + time stamp 
    328362 
    329363        KillWaves/Z cor1,bgd_temp,noadd_bgd,emp_temp,noadd_emp 
     364        Killwaves/Z tmp_a,tmp_b,tmp_c,tmp_d,c_val,d_val 
    330365        SetDataFolder root: 
    331366        Return(0) 
     
    338373 
    339374        //create the necessary wave references 
    340         WAVE sam_data=$"root:Packages:NIST:SAM:data" 
     375        WAVE sam_data=$"root:Packages:NIST:SAM:linear_data" 
    341376        WAVE sam_reals=$"root:Packages:NIST:SAM:realsread" 
    342377        WAVE sam_ints=$"root:Packages:NIST:SAM:integersread" 
    343378        WAVE/T sam_text=$"root:Packages:NIST:SAM:textread" 
    344         WAVE bgd_data=$"root:Packages:NIST:BGD:data" 
     379        WAVE bgd_data=$"root:Packages:NIST:BGD:linear_data" 
    345380        WAVE bgd_reals=$"root:Packages:NIST:BGD:realsread" 
    346381        WAVE bgd_ints=$"root:Packages:NIST:BGD:integersread" 
    347382        WAVE/T bgd_text=$"root:Packages:NIST:BGD:textread" 
    348         WAVE cor_data=$"root:Packages:NIST:COR:data" 
     383        WAVE cor_data=$"root:Packages:NIST:COR:linear_data" 
    349384        WAVE/T cor_text=$"root:Packages:NIST:COR:textread" 
     385 
     386        // needed to propagate error 
     387        WAVE cor_data_display=$"root:Packages:NIST:COR:data"            //just for the final copy 
     388        WAVE sam_err =$"root:Packages:NIST:SAM:linear_data_error" 
     389        WAVE bgd_err =$"root:Packages:NIST:BGD:linear_data_error" 
     390        WAVE cor_err =$"root:Packages:NIST:COR:linear_data_error" 
     391         
     392        Variable sam_trans_err 
     393        sam_trans_err = sam_reals[41] 
     394 
    350395         
    351396        //get sam and bgd attenuation factors 
     
    354399        Variable tmonsam,fsam,fbgd,xshift,yshift,rsam,csam,rbgd,cbgd,tmonbgd 
    355400        Variable wcen=0.001 
     401        Variable sam_atten_err,bgd_atten_err 
    356402        fileStr = sam_text[3] 
    357403        lambda = sam_reals[26] 
    358404        attenNo = sam_reals[3] 
    359         sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo) 
     405        sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,sam_atten_err) 
    360406        fileStr = bgd_text[3] 
    361407        lambda = bgd_reals[26] 
    362408        attenNo = bgd_reals[3] 
    363         bgd_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo) 
     409        bgd_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,bgd_atten_err) 
    364410         
    365411        //Print "atten = ",sam_attenFactor,bgd_attenFactor 
     
    398444        cor1 *= noadd_bgd               //zeros out regions where arrays do not overlap, one otherwise 
    399445 
     446// do the error propagation piecewise    
     447        Duplicate/O sam_err, tmp_a, tmp_b 
     448        tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2              //sig a ^2 
     449         
     450        tmp_b = (bgd_err/bgd_attenFactor)^2 + (bgd_atten_err*bgd_data/bgd_attenFactor^2)^2              //sig b ^2 
     451 
     452        cor_err = sqrt(tmp_a + tmp_b) 
     453 
     454 
    400455        //we're done, get out w/no error 
    401456        //set the COR_data to the result 
    402457        cor_data = cor1 
     458        cor_data_display = cor1 
     459 
    403460        //update COR header 
    404461        cor_text[1] = date() + " " + time()             //date + time stamp 
    405462 
    406463        KillWaves/Z cor1,bgd_temp,noadd_bgd 
     464        Killwaves/Z tmp_a,tmp_b 
     465 
    407466        SetDataFolder root: 
    408467        Return(0) 
     
    414473Function CorrectMode_3() 
    415474        //create the necessary wave references 
    416         WAVE sam_data=$"root:Packages:NIST:SAM:data" 
     475        WAVE sam_data=$"root:Packages:NIST:SAM:linear_data" 
    417476        WAVE sam_reals=$"root:Packages:NIST:SAM:realsread" 
    418477        WAVE sam_ints=$"root:Packages:NIST:SAM:integersread" 
    419478        WAVE/T sam_text=$"root:Packages:NIST:SAM:textread" 
    420         WAVE emp_data=$"root:Packages:NIST:EMP:data" 
     479        WAVE emp_data=$"root:Packages:NIST:EMP:linear_data" 
    421480        WAVE emp_reals=$"root:Packages:NIST:EMP:realsread" 
    422481        WAVE emp_ints=$"root:Packages:NIST:EMP:integersread" 
    423482        WAVE/T emp_text=$"root:Packages:NIST:EMP:textread" 
    424         WAVE cor_data=$"root:Packages:NIST:COR:data" 
     483        WAVE cor_data=$"root:Packages:NIST:COR:linear_data" 
    425484        WAVE/T cor_text=$"root:Packages:NIST:COR:textread" 
     485         
     486        // needed to propagate error 
     487        WAVE cor_data_display=$"root:Packages:NIST:COR:data"            //just for the final copy 
     488        WAVE sam_err =$"root:Packages:NIST:SAM:linear_data_error" 
     489        WAVE emp_err =$"root:Packages:NIST:EMP:linear_data_error" 
     490        WAVE cor_err =$"root:Packages:NIST:COR:linear_data_error" 
     491         
     492        Variable sam_trans_err,emp_trans_err 
     493        sam_trans_err = sam_reals[41] 
     494        emp_trans_err = emp_reals[41]    
    426495         
    427496        //get sam and bgd attenuation factors 
     
    430499        Variable tmonsam,fsam,femp,xshift,yshift,rsam,csam,remp,cemp,tmonemp 
    431500        Variable wcen=0.001,tsam,temp 
     501        Variable sam_atten_err,emp_atten_err 
    432502        fileStr = sam_text[3] 
    433503        lambda = sam_reals[26] 
    434504        attenNo = sam_reals[3] 
    435         sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo) 
     505        sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,sam_atten_err) 
    436506        fileStr = emp_text[3] 
    437507        lambda = emp_reals[26] 
    438508        attenNo = emp_reals[3] 
    439         emp_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo) 
     509        emp_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,emp_atten_err) 
    440510         
    441511        //get relative monitor counts (should all be 10^8, since normalized in add step) 
     
    478548        cor1 *= noadd_emp               //zeros out regions where arrays do not overlap, one otherwise 
    479549 
     550// do the error propagation piecewise    
     551        Duplicate/O sam_err, tmp_a, tmp_c ,c_val 
     552        tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2              //sig a ^2 
     553         
     554        tmp_c = (sam_trans_err*emp_data/(temp*emp_attenFactor))^2 + (emp_err*tsam/(temp*emp_attenFactor))^2 
     555        tmp_c += (tsam*emp_data*emp_trans_err/(temp*temp*emp_attenFactor))^2 + (tsam*emp_data*emp_atten_err/(temp*emp_attenFactor^2))^2//total of 6 terms 
     556 
     557        cor_err = sqrt(tmp_a + tmp_c) 
     558         
    480559        //we're done, get out w/no error 
    481560        //set the COR data to the result 
    482561        cor_data = cor1 
     562        cor_data_display = cor1 
     563 
    483564        //update COR header 
    484565        cor_text[1] = date() + " " + time()             //date + time stamp 
    485566 
    486567        KillWaves/Z cor1,emp_temp,noadd_emp 
     568        Killwaves/Z tmp_a,tmp_c,c_val 
     569 
    487570        SetDataFolder root: 
    488571        Return(0) 
     
    496579Function CorrectMode_4() 
    497580        //create the necessary wave references 
    498         WAVE sam_data=$"root:Packages:NIST:SAM:data" 
     581        WAVE sam_data=$"root:Packages:NIST:SAM:linear_data" 
    499582        WAVE sam_reals=$"root:Packages:NIST:SAM:realsread" 
    500583        WAVE sam_ints=$"root:Packages:NIST:SAM:integersread" 
    501584        WAVE/T sam_text=$"root:Packages:NIST:SAM:textread" 
    502585 
    503         WAVE cor_data=$"root:Packages:NIST:COR:data" 
     586        WAVE cor_data=$"root:Packages:NIST:COR:linear_data" 
    504587        WAVE/T cor_text=$"root:Packages:NIST:COR:textread" 
    505          
     588 
     589        // needed to propagate error 
     590        WAVE cor_data_display=$"root:Packages:NIST:COR:data"            //just for the final copy 
     591        WAVE sam_err =$"root:Packages:NIST:SAM:linear_data_error" 
     592        WAVE cor_err =$"root:Packages:NIST:COR:linear_data_error" 
     593                 
    506594        //get sam and bgd attenuation factors 
    507595        String fileStr="" 
    508         Variable lambda,attenNo,sam_AttenFactor 
     596        Variable lambda,attenNo,sam_AttenFactor,sam_atten_err 
    509597        fileStr = sam_text[3] 
    510598        lambda = sam_reals[26] 
    511599        attenNo = sam_reals[3] 
    512         sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo) 
     600        sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,sam_atten_err) 
    513601 
    514602        NVAR pixelsX = root:myGlobals:gNPixelsX 
     
    518606        cor1 = sam_data/sam_AttenFactor         //simply rescale the data 
    519607 
     608// do the error propagation piecewise    
     609        Duplicate/O sam_err, tmp_a 
     610        tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2              //sig a ^2 
     611 
     612        cor_err = sqrt(tmp_a) 
     613         
     614 
    520615        //we're done, get out w/no error 
    521616        //set the COR data to the result 
    522617        cor_data = cor1 
     618        cor_data_display = cor1 
     619 
    523620        //update COR header 
    524621        cor_text[1] = date() + " " + time()             //date + time stamp 
    525622 
    526623        KillWaves/Z cor1 
     624        Killwaves/Z tmp_a 
     625 
    527626        SetDataFolder root: 
    528627        Return(0) 
     
    531630Function CorrectMode_11() 
    532631        //create the necessary wave references 
    533         WAVE sam_data=$"root:Packages:NIST:SAM:data" 
     632        WAVE sam_data=$"root:Packages:NIST:SAM:linear_data" 
    534633        WAVE sam_reals=$"root:Packages:NIST:SAM:realsread" 
    535634        WAVE sam_ints=$"root:Packages:NIST:SAM:integersread" 
    536635        WAVE/T sam_text=$"root:Packages:NIST:SAM:textread" 
    537         WAVE bgd_data=$"root:Packages:NIST:BGD:data" 
     636        WAVE bgd_data=$"root:Packages:NIST:BGD:linear_data" 
    538637        WAVE bgd_reals=$"root:Packages:NIST:BGD:realsread" 
    539638        WAVE bgd_ints=$"root:Packages:NIST:BGD:integersread" 
    540639        WAVE/T bgd_text=$"root:Packages:NIST:BGD:textread" 
    541         WAVE emp_data=$"root:Packages:NIST:EMP:data" 
     640        WAVE emp_data=$"root:Packages:NIST:EMP:linear_data" 
    542641        WAVE emp_reals=$"root:Packages:NIST:EMP:realsread" 
    543642        WAVE emp_ints=$"root:Packages:NIST:EMP:integersread" 
    544643        WAVE/T emp_text=$"root:Packages:NIST:EMP:textread" 
    545         WAVE drk_data=$"root:Packages:NIST:DRK:data" 
     644        WAVE drk_data=$"root:Packages:NIST:DRK:linear_data" 
    546645        WAVE drk_reals=$"root:Packages:NIST:DRK:realsread" 
    547646        WAVE drk_ints=$"root:Packages:NIST:DRK:integersread" 
    548647        WAVE/T drk_text=$"root:Packages:NIST:DRK:textread" 
    549         WAVE cor_data=$"root:Packages:NIST:COR:data" 
     648        WAVE cor_data=$"root:Packages:NIST:COR:linear_data" 
    550649        WAVE/T cor_text=$"root:Packages:NIST:COR:textread" 
     650 
     651        // needed to propagate error 
     652        WAVE cor_data_display=$"root:Packages:NIST:COR:data"            //just for the final copy 
     653        WAVE sam_err =$"root:Packages:NIST:SAM:linear_data_error" 
     654        WAVE bgd_err =$"root:Packages:NIST:BGD:linear_data_error" 
     655        WAVE emp_err =$"root:Packages:NIST:EMP:linear_data_error" 
     656        WAVE drk_err =$"root:Packages:NIST:DRK:linear_data_error" 
     657        WAVE cor_err =$"root:Packages:NIST:COR:linear_data_error" 
     658         
     659        Variable sam_trans_err,emp_trans_err 
     660        sam_trans_err = sam_reals[41] 
     661        emp_trans_err = emp_reals[41] 
    551662         
    552663        //get sam and bgd attenuation factors 
     
    555666        Variable tmonsam,fsam,fbgd,xshift,yshift,rsam,csam,rbgd,cbgd,tmonbgd 
    556667        Variable wcen=0.001,tsam,temp,remp,cemp,tmonemp,femp,time_sam,time_drk,savmon_sam 
     668        Variable sam_atten_err,bgd_atten_err,emp_atten_err 
    557669        fileStr = sam_text[3] 
    558670        lambda = sam_reals[26] 
    559671        attenNo = sam_reals[3] 
    560         sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo) 
     672        sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,sam_atten_err) 
    561673        fileStr = bgd_text[3] 
    562674        lambda = bgd_reals[26] 
    563675        attenNo = bgd_reals[3] 
    564         bgd_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo) 
     676        bgd_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,bgd_atten_err) 
    565677        fileStr = emp_text[3] 
    566678        lambda = emp_reals[26] 
    567679        attenNo = emp_reals[3] 
    568         emp_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo) 
     680        emp_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,emp_atten_err) 
    569681         
    570682        //get relative monitor counts (should all be 10^8, since normalized in add step) 
     
    587699        NVAR pixelsY = root:myGlobals:gNPixelsY 
    588700        //rescale drk to sam cnt time and then multiply by the same monitor scaling as SAM 
    589         Make/D/O/N=(pixelsX,pixelsY) drk_temp 
     701        Make/D/O/N=(pixelsX,pixelsY) drk_temp, drk_tmp_err 
    590702        drk_temp = drk_data*(time_sam/time_drk)*(tmonsam/savmon_sam) 
     703        drk_tmp_err *= drk_err*(time_sam/time_drk)*(tmonsam/savmon_sam)                 //temporarily rescale the error of DRK 
    591704         
    592705        if(temp==0) 
     
    628741        cor1 *= noadd_bgd*noadd_emp             //zero out the array mismatch values 
    629742         
     743// do the error propagation piecewise    
     744        Duplicate/O sam_err, tmp_a, tmp_b, tmp_c, tmp_d,c_val,d_val 
     745        tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2              //sig a ^2 
     746         
     747        tmp_b = (bgd_err/bgd_attenFactor)^2*(tsam/temp - 1)^2 + (bgd_atten_err*bgd_data/bgd_attenFactor^2)^2*(1-tsam/temp)^2            //sig b ^2 
     748 
     749        tmp_c = (sam_trans_err/temp)^2*(emp_data/emp_attenFactor-bgd_data/bgd_attenFactor)^2 
     750        tmp_c += (tsam/temp^2)^2*emp_trans_err^2*(emp_data/emp_attenFactor-bgd_data/bgd_attenFactor)^2 
     751         
     752        tmp_d = (tsam/(temp*emp_attenFactor))^2*(emp_err)^2 + (tsam*emp_data/(temp*emp_attenFactor^2))^2*(emp_atten_err)^2 
     753 
     754        cor_err = sqrt(tmp_a + tmp_b + tmp_c + tmp_d + drk_tmp_err^2) 
     755         
    630756        //we're done, get out w/no error 
    631757        //set the COR data to the result 
    632758        cor_data = cor1 
     759        cor_data_display = cor1 
     760 
    633761        //update COR header 
    634762        cor_text[1] = date() + " " + time()             //date + time stamp 
    635763 
    636764        KillWaves/Z cor1,bgd_temp,noadd_bgd,emp_temp,noadd_emp,drk_temp 
     765        Killwaves/Z tmp_a,tmp_b,tmp_c,tmp_d,c_val,d_val,drk_tmp_err 
     766 
    637767        SetDataFolder root: 
    638768        Return(0) 
     
    643773Function CorrectMode_12() 
    644774        //create the necessary wave references 
    645         WAVE sam_data=$"root:Packages:NIST:SAM:data" 
     775        WAVE sam_data=$"root:Packages:NIST:SAM:linear_data" 
    646776        WAVE sam_reals=$"root:Packages:NIST:SAM:realsread" 
    647777        WAVE sam_ints=$"root:Packages:NIST:SAM:integersread" 
    648778        WAVE/T sam_text=$"root:Packages:NIST:SAM:textread" 
    649         WAVE bgd_data=$"root:Packages:NIST:BGD:data" 
     779        WAVE bgd_data=$"root:Packages:NIST:BGD:linear_data" 
    650780        WAVE bgd_reals=$"root:Packages:NIST:BGD:realsread" 
    651781        WAVE bgd_ints=$"root:Packages:NIST:BGD:integersread" 
    652782        WAVE/T bgd_text=$"root:Packages:NIST:BGD:textread" 
    653         WAVE drk_data=$"root:Packages:NIST:DRK:data" 
     783        WAVE drk_data=$"root:Packages:NIST:DRK:linear_data" 
    654784        WAVE drk_reals=$"root:Packages:NIST:DRK:realsread" 
    655785        WAVE drk_ints=$"root:Packages:NIST:DRK:integersread" 
    656786        WAVE/T drk_text=$"root:Packages:NIST:DRK:textread" 
    657         WAVE cor_data=$"root:Packages:NIST:COR:data" 
     787        WAVE cor_data=$"root:Packages:NIST:COR:linear_data" 
    658788        WAVE/T cor_text=$"root:Packages:NIST:COR:textread" 
     789 
     790        // needed to propagate error 
     791        WAVE cor_data_display=$"root:Packages:NIST:COR:data"            //just for the final copy 
     792        WAVE sam_err =$"root:Packages:NIST:SAM:linear_data_error" 
     793        WAVE bgd_err =$"root:Packages:NIST:BGD:linear_data_error" 
     794        WAVE drk_err =$"root:Packages:NIST:DRK:linear_data_error" 
     795        WAVE cor_err =$"root:Packages:NIST:COR:linear_data_error" 
     796         
     797        Variable sam_trans_err 
     798        sam_trans_err = sam_reals[41] 
     799         
    659800         
    660801        //get sam and bgd attenuation factors 
     
    663804        Variable tmonsam,fsam,fbgd,xshift,yshift,rsam,csam,rbgd,cbgd,tmonbgd 
    664805        Variable wcen=0.001,time_drk,time_sam,savmon_sam,tsam 
     806        Variable sam_atten_err,bgd_atten_err 
    665807        fileStr = sam_text[3] 
    666808        lambda = sam_reals[26] 
    667809        attenNo = sam_reals[3] 
    668         sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo) 
     810        sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,sam_atten_err) 
    669811        fileStr = bgd_text[3] 
    670812        lambda = bgd_reals[26] 
    671813        attenNo = bgd_reals[3] 
    672         bgd_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo) 
     814        bgd_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,bgd_atten_err) 
    673815         
    674816        //get relative monitor counts (should all be 10^8, since normalized in add step) 
     
    687829        NVAR pixelsY = root:myGlobals:gNPixelsY 
    688830        //rescale drk to sam cnt time and then multiply by the same monitor scaling as SAM 
    689         Make/D/O/N=(pixelsX,pixelsY) drk_temp 
     831        Make/D/O/N=(pixelsX,pixelsY) drk_temp,drk_tmp_err 
    690832        drk_temp = drk_data*(time_sam/time_drk)*(tmonsam/savmon_sam) 
    691          
     833        drk_tmp_err *= drk_err*(time_sam/time_drk)*(tmonsam/savmon_sam)                 //temporarily rescale the error of DRK 
     834 
    692835        // set up beamcenter shift, relative to SAM 
    693836        xshift = cbgd-csam 
     
    712855        cor1 *= noadd_bgd               //zeros out regions where arrays do not overlap, one otherwise 
    713856 
     857// do the error propagation piecewise    
     858        Duplicate/O sam_err, tmp_a, tmp_b 
     859        tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2              //sig a ^2 
     860         
     861        tmp_b = (bgd_err/bgd_attenFactor)^2 + (bgd_atten_err*bgd_data/bgd_attenFactor^2)^2              //sig b ^2 
     862 
     863        cor_err = sqrt(tmp_a + tmp_b + drk_tmp_err^2) 
     864 
    714865        //we're done, get out w/no error 
    715866        //set the COR_data to the result 
    716867        cor_data = cor1 
     868        cor_data_display = cor1 
     869 
    717870        //update COR header 
    718871        cor_text[1] = date() + " " + time()             //date + time stamp 
    719872 
    720 //      KillWaves/Z cor1,bgd_temp,noadd_bgd,drk_temp 
     873        KillWaves/Z cor1,bgd_temp,noadd_bgd,drk_temp 
     874        Killwaves/Z tmp_a,tmp_b,drk_tmp_err 
     875 
    721876        SetDataFolder root: 
    722877        Return(0) 
     
    729884Function CorrectMode_13() 
    730885        //create the necessary wave references 
    731         WAVE sam_data=$"root:Packages:NIST:SAM:data" 
     886        WAVE sam_data=$"root:Packages:NIST:SAM:linear_data" 
    732887        WAVE sam_reals=$"root:Packages:NIST:SAM:realsread" 
    733888        WAVE sam_ints=$"root:Packages:NIST:SAM:integersread" 
    734889        WAVE/T sam_text=$"root:Packages:NIST:SAM:textread" 
    735         WAVE emp_data=$"root:Packages:NIST:EMP:data" 
     890        WAVE emp_data=$"root:Packages:NIST:EMP:linear_data" 
    736891        WAVE emp_reals=$"root:Packages:NIST:EMP:realsread" 
    737892        WAVE emp_ints=$"root:Packages:NIST:EMP:integersread" 
    738893        WAVE/T emp_text=$"root:Packages:NIST:EMP:textread" 
    739         WAVE drk_data=$"root:DRK:data" 
    740         WAVE drk_reals=$"root:DRK:realsread" 
    741         WAVE drk_ints=$"root:DRK:integersread" 
    742         WAVE/T drk_text=$"root:DRK:textread" 
    743         WAVE cor_data=$"root:Packages:NIST:COR:data" 
     894        WAVE drk_data=$"root:Packages:NIST:DRK:linear_data" 
     895        WAVE drk_reals=$"root:Packages:NIST:DRK:realsread" 
     896        WAVE drk_ints=$"root:Packages:NIST:DRK:integersread" 
     897        WAVE/T drk_text=$"root:Packages:NIST:DRK:textread" 
     898        WAVE cor_data=$"root:Packages:NIST:COR:linear_data" 
    744899        WAVE/T cor_text=$"root:Packages:NIST:COR:textread" 
     900 
     901        // needed to propagate error 
     902        WAVE cor_data_display=$"root:Packages:NIST:COR:data"            //just for the final copy 
     903        WAVE sam_err =$"root:Packages:NIST:SAM:linear_data_error" 
     904        WAVE emp_err =$"root:Packages:NIST:EMP:linear_data_error" 
     905        WAVE drk_err =$"root:Packages:NIST:DRK:linear_data_error" 
     906        WAVE cor_err =$"root:Packages:NIST:COR:linear_data_error" 
     907         
     908        Variable sam_trans_err,emp_trans_err 
     909        sam_trans_err = sam_reals[41] 
     910        emp_trans_err = emp_reals[41] 
    745911         
    746912        //get sam and bgd attenuation factors (DRK irrelevant) 
     
    749915        Variable tmonsam,fsam,femp,xshift,yshift,rsam,csam,remp,cemp,tmonemp 
    750916        Variable wcen=0.001,tsam,temp,savmon_sam,time_sam,time_drk 
     917        Variable sam_atten_err,emp_atten_err 
    751918        fileStr = sam_text[3] 
    752919        lambda = sam_reals[26] 
    753920        attenNo = sam_reals[3] 
    754         sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo) 
     921        sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,sam_atten_err) 
    755922        fileStr = emp_text[3] 
    756923        lambda = emp_reals[26] 
    757924        attenNo = emp_reals[3] 
    758         emp_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo) 
     925        emp_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,emp_atten_err) 
    759926         
    760927        //get relative monitor counts (should all be 10^8, since normalized in add step) 
     
    774941        NVAR pixelsY = root:myGlobals:gNPixelsY 
    775942        //rescale drk to sam cnt time and then multiply by the same monitor scaling as SAM 
    776         Make/D/O/N=(pixelsX,pixelsY) drk_temp 
     943        Make/D/O/N=(pixelsX,pixelsY) drk_temp,drk_tmp_err 
    777944        drk_temp = drk_data*(time_sam/time_drk)*(tmonsam/savmon_sam) 
     945        drk_tmp_err *= drk_err*(time_sam/time_drk)*(tmonsam/savmon_sam)                 //temporarily rescale the error of DRK 
     946 
    778947         
    779948        if(temp==0) 
     
    805974        cor1 *= noadd_emp               //zeros out regions where arrays do not overlap, one otherwise 
    806975 
     976// do the error propagation piecewise    
     977        Duplicate/O sam_err, tmp_a, tmp_c, c_val 
     978        tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2              //sig a ^2 
     979         
     980        tmp_c = (sam_trans_err*emp_data/(temp*emp_attenFactor))^2 + (emp_err*tsam/(temp*emp_attenFactor))^2 
     981        tmp_c += (tsam*emp_data*emp_trans_err/(temp*temp*emp_attenFactor))^2 + (tsam*emp_data*emp_atten_err/(temp*emp_attenFactor^2))^2//total of 6 terms 
     982         
     983        cor_err = sqrt(tmp_a + tmp_c + drk_tmp_err^2) 
     984         
    807985        //we're done, get out w/no error 
    808986        //set the COR data to the result 
    809987        cor_data = cor1 
     988        cor_data_display = cor1 
     989 
    810990        //update COR header 
    811991        cor_text[1] = date() + " " + time()             //date + time stamp 
    812992 
    813993        KillWaves/Z cor1,emp_temp,noadd_emp,drk_temp 
     994        Killwaves/Z tmp_a,tmp_c,c_val,drk_tmp_err 
     995 
    814996        SetDataFolder root: 
    815997        Return(0) 
     
    8201002Function CorrectMode_14() 
    8211003        //create the necessary wave references 
    822         WAVE sam_data=$"root:Packages:NIST:SAM:data" 
     1004        WAVE sam_data=$"root:Packages:NIST:SAM:linear_data" 
    8231005        WAVE sam_reals=$"root:Packages:NIST:SAM:realsread" 
    8241006        WAVE sam_ints=$"root:Packages:NIST:SAM:integersread" 
    8251007        WAVE/T sam_text=$"root:Packages:NIST:SAM:textread" 
    826  
    827         WAVE drk_data=$"root:DRK:data" 
    828         WAVE drk_reals=$"root:DRK:realsread" 
    829         WAVE drk_ints=$"root:DRK:integersread" 
    830         WAVE/T drk_text=$"root:DRK:textread" 
    831         WAVE cor_data=$"root:Packages:NIST:COR:data" 
     1008        WAVE drk_data=$"root:Packages:NIST:DRK:linear_data" 
     1009        WAVE drk_reals=$"root:Packages:NIST:DRK:realsread" 
     1010        WAVE drk_ints=$"root:Packages:NIST:DRK:integersread" 
     1011        WAVE/T drk_text=$"root:Packages:NIST:DRK:textread" 
     1012        WAVE cor_data=$"root:Packages:NIST:COR:linear_data" 
    8321013        WAVE/T cor_text=$"root:Packages:NIST:COR:textread" 
     1014 
     1015        // needed to propagate error 
     1016        WAVE cor_data_display=$"root:Packages:NIST:COR:data"            //just for the final copy 
     1017        WAVE sam_err =$"root:Packages:NIST:SAM:linear_data_error" 
     1018        WAVE drk_err =$"root:Packages:NIST:DRK:linear_data_error" 
     1019        WAVE cor_err =$"root:Packages:NIST:COR:linear_data_error" 
     1020         
     1021        Variable sam_trans_err 
     1022        sam_trans_err = sam_reals[41] 
     1023         
    8331024         
    8341025        //get sam and bgd attenuation factors 
     
    8371028        Variable tmonsam,fsam,fbgd,xshift,yshift,rsam,csam,rbgd,cbgd,tmonbgd 
    8381029        Variable wcen=0.001,time_drk,time_sam,savmon_sam,tsam 
     1030        Variable sam_atten_err 
    8391031        fileStr = sam_text[3] 
    8401032        lambda = sam_reals[26] 
    8411033        attenNo = sam_reals[3] 
    842         sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo) 
     1034        sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,sam_atten_err) 
    8431035         
    8441036        //get relative monitor counts (should all be 10^8, since normalized in add step) 
     
    8551047        NVAR pixelsY = root:myGlobals:gNPixelsY 
    8561048        //rescale drk to sam cnt time and then multiply by the same monitor scaling as SAM 
    857         Make/D/O/N=(pixelsX,pixelsY) drk_temp 
     1049        Make/D/O/N=(pixelsX,pixelsY) drk_temp,drk_tmp_err 
    8581050        drk_temp = drk_data*(time_sam/time_drk)*(tmonsam/savmon_sam) 
    859          
     1051        drk_tmp_err *= drk_err*(time_sam/time_drk)*(tmonsam/savmon_sam)                 //temporarily rescale the error of DRK 
     1052 
    8601053        Make/D/O/N=(pixelsX,pixelsY) cor1       //temp arrays 
    8611054        //always ignore the DRK center shift 
     
    8681061        cor1 = fsam*sam_data/sam_AttenFactor  - drk_temp/sam_attenFactor 
    8691062 
     1063// do the error propagation piecewise    
     1064        Duplicate/O sam_err, tmp_a 
     1065        tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2              //sig a ^2 
     1066 
     1067        cor_err = sqrt(tmp_a + drk_tmp_err^2) 
     1068         
    8701069        //we're done, get out w/no error 
    8711070        //set the COR_data to the result 
    8721071        cor_data = cor1 
     1072        cor_data_display = cor1 
     1073 
    8731074        //update COR header 
    8741075        cor_text[1] = date() + " " + time()             //date + time stamp 
    8751076 
    876 //      KillWaves/Z cor1,bgd_temp,noadd_bgd,drk_temp 
     1077        KillWaves/Z cor1,bgd_temp,noadd_bgd,drk_temp 
     1078        Killwaves/Z tmp_a,tmp_b,tmp_c,tmp_d,c_val,d_val,drk_tmp_err 
     1079 
    8771080        SetDataFolder root: 
    8781081        Return(0) 
     
    11381341        Variable SamAttenFactor,lambda,attenNo,err=0 
    11391342        String samfileStr="" 
     1343        Variable sam_atten_err 
    11401344        samfileStr = sam_text[3] 
    11411345        lambda = sam_reals[26] 
    11421346        attenNo = sam_reals[3] 
    1143         SamAttenFactor = AttenuationFactor(samFileStr,lambda,AttenNo) 
     1347        SamAttenFactor = AttenuationFactor(samFileStr,lambda,AttenNo,sam_atten_err) 
    11441348        //if sample trans is zero, do only SAM-BGD subtraction (notify the user) 
    11451349        Variable sam_trans = sam_reals[4] 
Note: See TracChangeset for help on using the changeset viewer.