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]

Location:
sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS
Files:
9 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] 

  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/Marquee.ipf

    r788 r794  
    2626// accepts arbitrary detector coordinates. calling function is responsible for  
    2727// keeping selection in bounds 
    28 Function SumCountsInBox(x1,x2,y1,y2,type) 
    29         Variable x1,x2,y1,y2 
     28Function SumCountsInBox(x1,x2,y1,y2,ct_err,type) 
     29        Variable x1,x2,y1,y2,&ct_err 
    3030        String type 
    3131         
    32         Variable counts = 0,ii,jj 
     32        Variable counts = 0,ii,jj,err2_sum 
    3333         
    3434        String dest =  "root:Packages:NIST:"+type 
    … …  
    4141                wave w=$(dest + ":data") 
    4242        endif 
    43          
     43        wave data_err = $(dest + ":linear_data_error") 
     44         
     45        err2_sum = 0            // running total of the squared error 
    4446        ii=x1 
    4547        jj=y1 
    … …  
    4749                do 
    4850                        counts += w[ii][jj] 
     51                        err2_sum += data_err[ii][jj]*data_err[ii][jj] 
    4952                        jj+=1 
    5053                while(jj<=y2) 
    … …  
    5255                ii+=1 
    5356        while(ii<=x2) 
     57         
     58        err2_sum = sqrt(err2_sum) 
     59        ct_err = err2_sum 
     60         
     61//      Print "error = ",err2_sum 
     62//      Print "error/counts = ",err2_sum/counts 
     63         
    5464         
    5565        Return (counts) 
    … …  
    113123        //to the same monitor counts and corrected for detector deadtime 
    114124        String type = "SAM" 
    115         Variable counts 
    116         counts = SumCountsInBox(x1,x2,y1,y2,type) 
    117         Print " marquee counts =",counts 
     125        Variable counts,ct_err 
     126        counts = SumCountsInBox(x1,x2,y1,y2,ct_err,type) 
     127//      Print "marquee counts =",counts 
     128//      Print "relative error = ",ct_err/counts 
     129         
    118130        //Set the global gTransCts 
    119131        Variable/G root:myGlobals:Patch:gTransCts = counts 
    … …  
    148160        WriteBoxCountsToHeader(filename,counts) 
    149161         
     162        WriteBoxCountsErrorToHeader(filename,ct_err) 
     163         
     164        return(0) 
    150165End 
    151166 
    … …  
    315330        //parse the list of file numbers 
    316331        String fileList="",item="",pathStr="",fullPath="" 
    317         Variable ii,num,err,cts 
     332        Variable ii,num,err,cts,ct_err 
    318333         
    319334        PathInfo catPathName 
    … …  
    330345        //sum over the box 
    331346        //print the results 
    332         Make/O/N=(num) FileID,BoxCounts 
     347        Make/O/N=(num) FileID,BoxCounts,BoxCount_err 
    333348        Print "Results are stored in root:FileID and root:BoxCounts waves" 
    334349        for(ii=0;ii<num;ii+=1) 
    … …  
    344359                String/G root:myGlobals:gDataDisplayType=type 
    345360                fRawWindowHook() 
    346                 cts=SumCountsInBox(x1,x2,y1,y2,type) 
     361                cts=SumCountsInBox(x1,x2,y1,y2,ct_err,type) 
    347362                BoxCounts[ii]=cts 
     363                BoxCount_err[ii]=ct_err 
    348364                Print item+" counts = ",cts 
    349365        endfor 
    350366         
    351         DoBoxGraph(FileID,BoxCounts) 
     367        DoBoxGraph(FileID,BoxCounts,BoxCount_err) 
    352368         
    353369        return(0) 
    354370End 
    355371 
    356 Function DoBoxGraph(FileID,BoxCounts) 
    357         Wave FileID,BoxCounts 
     372Function DoBoxGraph(FileID,BoxCounts,BoxCount_err) 
     373        Wave FileID,BoxCounts,BoxCount_err 
    358374         
    359375        Sort FileID BoxCounts,FileID            //sort the waves, in case the run numbers were entered out of numerical order 
    … …  
    364380        ModifyGraph grid=2 
    365381        ModifyGraph mirror=2 
     382        ErrorBars/T=0 BoxCounts Y,wave=(BoxCount_err,BoxCount_err) 
    366383        Label left "Counts (per 10^8 monitor counts)" 
    367384        Label bottom "Run Number" 
    … …  
    561578                Print "There is no Marquee" 
    562579        else 
    563                 Variable count,x1,x2,y1,y2 
     580                Variable count,x1,x2,y1,y2,ct_err 
    564581                x1 = V_left 
    565582                x2 = V_right 
    … …  
    574591                KeepSelectionInBounds(x1,x2,y1,y2) 
    575592                SVAR cur_folder=root:myGlobals:gDataDisplayType 
    576                 count = SumCountsInBox(x1,x2,y1,y2,cur_folder) 
    577                 Print "counts = "+ num2str(count) 
     593                count = SumCountsInBox(x1,x2,y1,y2,ct_err,cur_folder) 
     594                Print "counts = ",count 
     595                Print "err/counts = ",ct_err/count 
    578596        endif 
    579597End 

  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/NCNR_DataReadWrite.ipf

    r776 r794  
    369369         
    370370        Duplicate/O data linear_data                    // at this point, the data is still the raw data, and is linear_data 
     371         
     372        // proper error for counting statistics, good for low count values too 
     373        // rather than just sqrt(n) 
     374        // see N. Gehrels, Astrophys. J., 303 (1986) 336-346, equation (7) 
     375        // for S = 1 in eq (7), this corresponds to one sigma error bars 
     376        Duplicate/O linear_data linear_data_error 
     377        linear_data_error = 1 + sqrt(linear_data + 0.75)                                 
    371378        // 
    372379         
    … …  
    759766         
    760767        //read in the data 
    761          GBLoadStr="GBLoadWave/O/N=tempGBwave/T={2,2}/J=2/W=1/Q" 
     768        GBLoadStr="GBLoadWave/O/N=tempGBwave/T={2,2}/J=2/W=1/Q" 
    762769 
    763770        curPath = "root:Packages:NIST:"+cur_folder 
    … …  
    904911        Redimension/N=(pixelsX,pixelsY) data            //,linear_data 
    905912         
     913        Duplicate/O data linear_data_error 
     914        linear_data_error = 1 + sqrt(data + 0.75) 
     915         
     916        //just in case there are odd inputs to this, like negative intensities 
     917        WaveStats/Q linear_data_error 
     918        linear_data_error = numtype(linear_data_error[p]) == 0 ? linear_data_error[p] : V_avg 
     919        linear_data_error = linear_data_error[p] != 0 ? linear_data_error[p] : V_avg 
     920         
    906921        //linear_data = data 
    907922         
    … …  
    10821097End 
    10831098 
     1099//sample transmission error (one sigma) is a real value at byte 396 
     1100Function WriteTransmissionErrorToHeader(fname,transErr) 
     1101        String fname 
     1102        Variable transErr 
     1103         
     1104        WriteVAXReal(fname,transErr,396)                //transmission start byte is 396 
     1105        return(0) 
     1106End 
     1107 
     1108 
    10841109//whole transmission is a real value at byte 392 
    10851110Function WriteWholeTransToHeader(fname,trans) 
    … …  
    10971122         
    10981123        WriteVAXReal(fname,counts,494)          // start byte is 494 
     1124        return(0) 
     1125End 
     1126 
     1127//box sum counts error is is a real value at byte 400 
     1128Function WriteBoxCountsErrorToHeader(fname,rel_err) 
     1129        String fname 
     1130        Variable rel_err 
     1131         
     1132        WriteVAXReal(fname,rel_err,400)         // start byte is 400 
    10991133        return(0) 
    11001134End 
    … …  
    14931527end 
    14941528 
     1529//transmission error (one sigma) is at byte 396 
     1530Function getSampleTransError(fname) 
     1531        String fname 
     1532         
     1533        return(getRealValueFromHeader(fname,396)) 
     1534end 
     1535 
    14951536//box counts are stored at byte 494 
    14961537Function getBoxCounts(fname) 
    … …  
    14991540        return(getRealValueFromHeader(fname,494)) 
    15001541end 
     1542 
     1543//box counts error are stored at byte 400 
     1544Function getBoxCountsError(fname) 
     1545        String fname 
     1546         
     1547        return(getRealValueFromHeader(fname,400)) 
     1548end 
     1549 
    15011550 
    15021551//whole detector trasmission is at byte 392 

  • 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:                                                         

  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/ProDiv.ipf

    r665 r794  
    3838         
    3939        WAVE data=$("root:Packages:NIST:"+type+":data") 
     40        WAVE data_lin=$("root:Packages:NIST:"+type+":linear_data") 
     41        WAVE data_err=$("root:Packages:NIST:"+type+":linear_data_error") 
     42         
    4043        Variable totCts=sum(data,Inf,-Inf)              //sum all of the data 
    4144        NVAR pixelX = root:myGlobals:gNPixelsX 
    … …  
    4649        data *= pixelX*pixelY 
    4750         
     51        data_lin /= totCts 
     52        data_lin *= pixelX*pixelY 
     53         
     54        data_err /= totCts 
     55        data_err *= pixelX*pixelY 
     56                 
    4857        return(0) 
    4958End 
    … …  
    132141        WAVE ctrData=$("root:Packages:NIST:"+ctrtype+":data") 
    133142        WAVE offData=$("root:Packages:NIST:"+offtype+":data") 
     143         
     144        WAVE ctrData_lin=$("root:Packages:NIST:"+ctrtype+":linear_data") 
     145        WAVE offData_lin=$("root:Packages:NIST:"+offtype+":linear_data") 
     146         
     147        WAVE ctrData_err=$("root:Packages:NIST:"+ctrtype+":linear_data_error") 
     148        WAVE offData_err=$("root:Packages:NIST:"+offtype+":linear_data_error") 
     149         
    134150        Variable ii,jj 
    135151         
    … …  
    137153                for(jj=y1;jj<=y2;jj+=1) 
    138154                        ctrData[ii][jj] = offData[ii][jj] 
     155                        ctrData_lin[ii][jj] = offData_lin[ii][jj] 
     156                        ctrData_err[ii][jj] = offData_err[ii][jj] 
    139157                endfor 
    140158        endfor 

  • 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         

  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/Transmission.ipf

    r670 r794  
    8080        Make/O/D/N=0 $"root:myGlobals:TransHeaderInfo:S_Lambda" 
    8181        Make/O/D/N=0 $"root:myGlobals:TransHeaderInfo:S_Transmission" 
     82        Make/O/D/N=0 $"root:myGlobals:TransHeaderInfo:S_Trans_Error" 
     83 
    8284        Make/O/D/N=0 $"root:myGlobals:TransHeaderInfo:S_Whole" 
    8385        If(V_Flag==0) 
    … …  
    177179        Wave   S_Lambda = $"root:myGlobals:TransHeaderInfo:S_Lambda" 
    178180        Wave   S_Transmission = $"root:myGlobals:TransHeaderInfo:S_Transmission" 
    179          
    180         Edit/K=1/W=(50,50,520,270) S_TRANS_Filenames, S_Filenames, S_Labels, S_SDD, S_Lambda, S_Transmission as "ScatteringFiles" 
     181        Wave   S_GTrans_err =  $"root:myGlobals:TransHeaderInfo:S_Trans_Error" 
     182 
     183         
     184        Edit/K=1/W=(50,50,520,270) S_TRANS_Filenames, S_Filenames, S_Labels, S_SDD, S_Lambda, S_Transmission, S_GTrans_err as "ScatteringFiles" 
    181185        String name="ScatterFileTable" 
    182186        DoWindow/C $name 
    … …  
    207211        Wave   S_GLambda = $"root:myGlobals:TransHeaderInfo:S_Lambda" 
    208212        Wave   S_GTransmission = $"root:myGlobals:TransHeaderInfo:S_Transmission" 
     213        Wave     S_GTrans_err =  $"root:myGlobals:TransHeaderInfo:S_Trans_Error" 
    209214        Wave   S_Whole= $"root:myGlobals:TransHeaderInfo:S_Whole" 
    210215         
    211216        Sort T_GSuffix, T_GEMP_Filenames, T_GSuffix, T_GSuffices, T_GFilenames, T_GLabels, T_GSDD, T_GLambda, T_GTransmission, T_Whole 
    212         Sort S_GSuffix, S_GTRANS_Filenames, S_GSuffix, S_GSuffices, S_GFilenames, S_GLabels, S_GSDD, S_GLambda, S_GTransmission, S_Whole 
     217        Sort S_GSuffix, S_GTRANS_Filenames, S_GSuffix, S_GSuffices, S_GFilenames, S_GLabels, S_GSDD, S_GLambda, S_GTransmission, S_Whole, S_GTrans_err 
    213218         
    214219End 
    … …  
    236241        Wave   S_GLambda = $"root:myGlobals:TransHeaderInfo:S_Lambda" 
    237242        Wave   S_GTransmission = $"root:myGlobals:TransHeaderInfo:S_Transmission" 
     243        Wave   S_GTrans_err =  $"root:myGlobals:TransHeaderInfo:S_Trans_Error" 
    238244        Wave   S_Whole= $"root:myGlobals:TransHeaderInfo:S_Whole" 
    239245         
    240246        Sort T_GLabels, T_GEMP_Filenames, T_GSuffix, T_GSuffices, T_GFilenames, T_GLabels, T_GSDD, T_GLambda, T_GTransmission, T_Whole 
    241         Sort S_GLabels, S_GTRANS_Filenames, S_GSuffix, S_GSuffices, S_GFilenames, S_GLabels, S_GSDD, S_GLambda, S_GTransmission, S_Whole 
     247        Sort S_GLabels, S_GTRANS_Filenames, S_GSuffix, S_GSuffices, S_GFilenames, S_GLabels, S_GSDD, S_GLambda, S_GTransmission, S_Whole, S_GTrans_err 
    242248         
    243249End 
    … …  
    279285                Wave GLambda = $"root:myGlobals:TransHeaderInfo:S_Lambda" 
    280286                Wave GTransmission = $"root:myGlobals:TransHeaderInfo:S_Transmission" 
     287                Wave S_GTrans_err =  $"root:myGlobals:TransHeaderInfo:S_Trans_Error" 
    281288                Wave GWhole = $"root:myGlobals:TransHeaderInfo:S_Whole" 
     289                 
     290                        //Transmission error 
     291                        InsertPoints lastPoint,1,S_GTrans_err 
     292                        S_GTrans_err[lastPoint] = getSampleTransError(fname) 
     293                         
    282294        Endif 
    283295 
    … …  
    306318        InsertPoints lastPoint,1,GTransmission 
    307319        GTransmission[lastPoint]=getSampleTrans(fname) 
     320         
     321 
    308322         
    309323        //Whole detector Transmission 
    … …  
    477491        Make/O/D/N=0 $"root:myGlobals:TransHeaderInfo:S_Lambda" 
    478492        Make/O/D/N=0 $"root:myGlobals:TransHeaderInfo:S_Transmission" 
     493        Make/O/D/N=0 $"root:myGlobals:TransHeaderInfo:S_Trans_Error" 
    479494        Make/O/D/N=0 $"root:myGlobals:TransHeaderInfo:S_Whole" 
    480495End 
    … …  
    605620        Wave/T S_GSuffix = $"root:myGlobals:TransHeaderInfo:S_Suffix" 
    606621        Wave S_GTransmission =  $"root:myGlobals:TransHeaderInfo:S_Transmission" 
     622        Wave S_GTrans_err =  $"root:myGlobals:TransHeaderInfo:S_Trans_Error" 
    607623        Wave/T T_EMP_Filenames = $"root:myGlobals:TransHeaderInfo:T_EMP_Filenames" 
    608624        Wave/T T_GFilenames = $"root:myGlobals:TransHeaderInfo:T_Filenames" 
    … …  
    643659                        WriteWholeTransToHeader(filename,1)             //WholeTrans start byte is 392 
    644660                         
     661                        WriteTransmissionErrorToHeader(filename,0)                      //reset transmission error to zero 
     662                        WriteBoxCountsErrorToHeader(filename,0) 
     663                         
    645664                        //then update the table that is displayed 
    646665                        T_EMP_Filenames[ii] = "" 
    … …  
    663682                        //write a trans==1 to the file header of the raw data (open/close done in function) 
    664683                        WriteTransmissionToHeader(filename,1)           //sample trans 
     684                        WriteTransmissionErrorToHeader(filename,0)                      //reset transmission error to zero 
     685                        WriteBoxCountsErrorToHeader(filename,0) 
    665686                         
    666687                        //then update the table that is displayed 
    667688                        S_TRANS_Filenames[ii] = "" 
    668689                        S_GTransmission[ii] = 1 
     690                        S_GTrans_err[ii] = 0 
    669691                         
    670692                        ii+=1 
    … …  
    710732        Wave/T S_GFilenames = $"root:myGlobals:TransHeaderInfo:S_Filenames" 
    711733        Wave S_GTransmission =  $"root:myGlobals:TransHeaderInfo:S_Transmission" 
    712          
     734        Wave S_GTrans_err =  $"root:myGlobals:TransHeaderInfo:S_Trans_Error" 
     735 
    713736        Variable num_s_files, num_t_files, ii, jj 
    714         Variable refnum, transCts, emptyCts,attenRatio,lambda,trans 
    715         Variable x1,x2,y1,y2,err,attenEmp,attenSam 
     737        Variable refnum, transCts, emptyCts,attenRatio,lambda,trans,sam_atten_err,emp_atten_err 
     738        Variable x1,x2,y1,y2,err,attenEmp,attenSam,empty_ct_err,sam_ct_err,samAttenFactor,empAttenFactor,trans_err 
    716739        String suffix = "",pathName,textStr,abortStr,emptyFile,transFile,samFileStr 
    717740        String GBLoadStr="GBLoadWave/O/N=tempGBwave/T={2,2}/J=2/W=1/Q" 
    … …  
    742765                                                //read the real count value 
    743766                                                emptyCts = getBoxCounts(emptyFile) 
     767                                                // get the error in count value 
     768                                                empty_ct_err = getBoxCountsError(emptyFile) 
     769                                                 
    744770                                                // read the attenuator number of the empty beam file 
    745771                                                attenEmp = getAttenNumber(emptyFile) 
    … …  
    757783                                                err = Raw_to_work("SAM") 
    758784                                                //sum region in SAM 
    759                                                 transCts =  SumCountsInBox(x1,x2,y1,y2,"SAM")    
     785                                                transCts =  SumCountsInBox(x1,x2,y1,y2,sam_ct_err,"SAM")         
    760786                                                // get the attenuator, lambda, and sample string (to get the instrument) 
    761787                                                WAVE/T samText = $"root:Packages:NIST:SAM:textRead" 
    … …  
    764790                                                lambda = samReals[26] 
    765791                                                attenSam = samReals[3] 
     792                                                 
    766793                                                //calculate the ratio of attenuation factors - assumes that same instrument used for each, AND same lambda 
    767                                                 AttenRatio = AttenuationFactor(samFileStr,lambda,attenEmp)/AttenuationFactor(samFileStr,lambda,attenSam) 
     794                                                samAttenFactor = AttenuationFactor(samFileStr,lambda,attenSam,sam_atten_err) 
     795                                                empAttenFactor = AttenuationFactor(samFileStr,lambda,attenEmp,emp_atten_err) 
     796                                                AttenRatio = empAttenFactor/samAttenFactor 
    768797                                                //calculate trans based on empty beam value and rescale by attenuation ratio 
    769798                                                trans= transCts/emptyCts * AttenRatio 
    770799                                                 
     800                                                // squared, relative error 
     801                                                if(AttenRatio == 1) 
     802                                                        trans_err = (sam_ct_err/transCts)^2 + (empty_ct_err/emptyCts)^2         //same atten, att_err drops out 
     803                                                else 
     804                                                        trans_err = (sam_ct_err/transCts)^2 + (empty_ct_err/emptyCts)^2 + (sam_atten_err/samAttenFactor)^2 + (emp_atten_err/empAttenFactor)^2 
     805                                                endif 
     806                                                trans_err = sqrt(trans_err) 
     807                                                trans_err *= trans              // now, one std deviation 
     808                                                 
    771809                                                //write out counts and transmission to history window, showing the attenuator ratio, if it is not unity 
    772810                                                If(attenRatio==1) 
    773                                                         Printf "%s\t\tTrans Counts = %g\tTrans = %g\r",S_GFilenames[ii], transCts,trans 
     811                                                        Printf "%s\t\tTrans Counts = %g\tTrans = %g +/- %g\r",S_GFilenames[ii], transCts,trans,trans_err 
    774812                                                else 
    775                                                         Printf "%s\t\tTrans Counts = %g\tTrans = %g\tAttenuatorRatio = %g\r",S_GFilenames[ii], transCts,trans,attenRatio 
     813                                                        Printf "%s\t\tTrans Counts = %g\tTrans = %g +/- %g\tAttenuatorRatio = %g\r",S_GFilenames[ii], transCts,trans,trans_err,attenRatio 
    776814                                                endif 
    777815                                                //write the trans to the file header of the raw data (open/close done in function) 
    778816                                                WriteTransmissionToHeader(filename,trans) 
    779817                                                 
     818                                                WriteTransmissionErrorToHeader(filename,trans_err) 
     819/////                                            
    780820                                                //then update the global that is displayed 
    781821                                                S_GTransmission[ii] = trans 
     822                                                S_GTrans_err[ii] = trans_err 
     823 
    782824                                                 
    783825                                        else  // There is no empty assigned 
    … …  
    821863         
    822864        Variable num_t_files, ii, jj 
    823         Variable refnum, transCts, emptyCts,attenRatio,lambda,trans 
    824         Variable x1,x2,y1,y2,err,attenEmp,attenSam 
     865        Variable refnum, transCts, emptyCts,attenRatio,lambda,trans,samAttenFactor,empAttenFactor,trans_err 
     866        Variable x1,x2,y1,y2,err,attenEmp,attenSam,empty_ct_err,sam_ct_err,emp_atten_err,sam_atten_err 
    825867        String suffix = "",pathName,textStr,abortStr,emptyFile,transFile,samFileStr 
    826868        String GBLoadStr="GBLoadWave/O/N=tempGBwave/T={2,2}/J=2/W=1/Q" 
    … …  
    851893                                                //read the real count value 
    852894                                                emptyCts = getBoxCounts(emptyFile) 
     895                                                // get the count error 
     896                                                empty_ct_err = getBoxCountsError(emptyFile) 
     897                                                 
    853898                                                // read the attenuator number of the empty beam file 
    854899                                                attenEmp = getAttenNumber(emptyFile) 
     900                                                 
    855901                                                // 
    856902                                                if( ((x1-x2)==0) || ((y1-y2)==0) )              //zero width marquee in either direction 
    … …  
    866912                                                err = Raw_to_work("SAM") 
    867913                                                //sum region in SAM 
    868                                                 transCts =  SumCountsInBox(x1,x2,y1,y2,"SAM")    
     914                                                transCts =  SumCountsInBox(x1,x2,y1,y2,sam_ct_err,"SAM")         
    869915                                                // get the attenuator, lambda, and sample string (to get the instrument) 
    870916                                                WAVE/T samText = $"root:Packages:NIST:SAM:textRead" 
    … …  
    873919                                                lambda = samReals[26] 
    874920                                                attenSam = samReals[3] 
     921                                                 
    875922                                                //calculate the ratio of attenuation factors - assumes that same instrument used for each, AND same lambda 
    876                                                 AttenRatio = AttenuationFactor(samFileStr,lambda,attenEmp)/AttenuationFactor(samFileStr,lambda,attenSam) 
     923                                                samAttenFactor = AttenuationFactor(samFileStr,lambda,attenSam,sam_atten_err) 
     924                                                empAttenFactor = AttenuationFactor(samFileStr,lambda,attenEmp,emp_atten_err) 
     925                                                AttenRatio = empAttenFactor/samAttenFactor 
    877926                                                //calculate trans based on empty beam value and rescale by attenuation ratio 
    878927                                                trans= transCts/emptyCts * AttenRatio 
     928                                                 
     929                                                // squared, relative error 
     930                                                if(AttenRatio == 1) 
     931                                                        trans_err = (sam_ct_err/transCts)^2 + (empty_ct_err/emptyCts)^2         //same atten, att_err drops out 
     932                                                else 
     933                                                        trans_err = (sam_ct_err/transCts)^2 + (empty_ct_err/emptyCts)^2 + (sam_atten_err/samAttenFactor)^2 + (emp_atten_err/empAttenFactor)^2 
     934                                                endif 
     935                                                trans_err = sqrt(trans_err) 
     936                                                trans_err *= trans              // now, one std deviation                                                
    879937                                                 
    880938                                                //write out counts and transmission to history window, showing the attenuator ratio, if it is not unity 
    881939                                                If(attenRatio==1) 
    882940                                                        //Printf "%s\t\tTrans Counts = %g\t Actual Trans = %g\r",T_GFilenames[ii], transCts,trans 
    883                                                         Printf "%s\t\tBox Counts = %g\t Trans = %g\r",T_GFilenames[ii], transCts,trans 
     941                                                        Printf "%s\t\tBox Counts = %g\t Trans = %g +/- %g\r",T_GFilenames[ii], transCts,trans,trans_err 
    884942                                                else 
    885943                                                        //Printf "%s\t\tTrans Counts = %g\t Trans = %g\tAttenuatorRatio = %g\r",T_GFilenames[ii], transCts,trans,attenRatio 
    886                                                         Printf "%s\t\tBox Counts = %g\t Trans = %g\t AttenuatorRatio = %g\r",T_GFilenames[ii], transCts,trans,attenRatio 
     944                                                        Printf "%s\t\tBox Counts = %g\t Trans = %g +/- %g\t AttenuatorRatio = %g\r",T_GFilenames[ii], transCts,trans,trans_err,attenRatio 
    887945                                                endif 
    888946                                                //write the trans to the file header of the raw data (open/close done in function) 
    889947                                                WriteTransmissionToHeader(filename,trans)               //transmission start byte is 158 
    890948                                                 
     949                                                WriteTransmissionErrorToHeader(filename,trans_err) 
     950                                                                                                 
    891951                                                //then update the global that is displayed 
    892952                                                T_GTransmission[ii] = trans 
    … …  
    936996        Variable num_t_files, ii, jj 
    937997        Variable refnum, transCts, emptyCts,attenRatio,lambda,trans 
    938         Variable x1,x2,y1,y2,err,attenEmp,attenSam 
     998        Variable x1,x2,y1,y2,err,attenEmp,attenSam,empty_ct_err,sam_ct_err,emp_atten_err,sam_atten_err 
    939999        String suffix = "",pathName,textStr,abortStr,emptyFile,transFile,samFileStr 
    9401000        String GBLoadStr="GBLoadWave/O/N=tempGBwave/T={2,2}/J=2/W=1/Q" 
    … …  
    9651025                                                //read the real count value 
    9661026                                                emptyCts = getBoxCounts(emptyFile) 
     1027                                                // get the box count error 
     1028                                                empty_ct_err = getBoxCountsError(emptyFile) 
     1029                                                 
    9671030                                                // read the attenuator number of the empty beam file 
    9681031                                                attenEmp = getAttenNumber(emptyFile) 
     1032                                                 
    9691033                                                // 
    9701034                                                if( ((x1-x2)==0) || ((y1-y2)==0) )              //zero width marquee in either direction 
    … …  
    9801044                                                err = Raw_to_work("SAM") 
    9811045                                                //sum region in SAM 
    982                                                 transCts =  SumCountsInBox(0,pixelsX-1,0,pixelsY-1,"SAM")        
     1046                                                transCts =  SumCountsInBox(0,pixelsX-1,0,pixelsY-1,sam_ct_err,"SAM")     
    9831047                                                // get the attenuator, lambda, and sample string (to get the instrument) 
    9841048                                                WAVE/T samText = $"root:Packages:NIST:SAM:textRead" 
    … …  
    9871051                                                lambda = samReals[26] 
    9881052                                                attenSam = samReals[3] 
     1053                                                 
    9891054                                                //calculate the ratio of attenuation factors - assumes that same instrument used for each, AND same lambda 
    990                                                 AttenRatio = AttenuationFactor(samFileStr,lambda,attenEmp)/AttenuationFactor(samFileStr,lambda,attenSam) 
     1055                                                AttenRatio = AttenuationFactor(samFileStr,lambda,attenEmp,emp_atten_err)/AttenuationFactor(samFileStr,lambda,attenSam,sam_atten_err) 
    9911056                                                //calculate trans based on empty beam value and rescale by attenuation ratio 
    9921057                                                trans= transCts/emptyCts * AttenRatio 

  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/WorkFileUtils.ipf

    r788 r794  
    2020//*************************** 
    2121 
     22// 
     23// 
     24// MARCH 2011 - changed the references that manipulate the data to explcitly work on the linear_data wave 
     25//                                      at the end of routines that maniplulate linear_data, it is copied back to "data" which is displayed 
     26// 
     27// 
    2228 
    2329//testing procedure, not called anymore 
    … …  
    6369         
    6470        // if the desired workfile doesn't exist, let the user know, and just make a new one 
    65         destPath = "root:Packages:NIST:"+newType + ":data" 
    66         if(WaveExists($destpath) == 0) 
     71        if(WaveExists($("root:Packages:NIST:"+newType + ":data")) == 0) 
    6772                Print "There is no old work file to add to - a new one will be created" 
    6873                //call Raw_to_work(), then return from this function 
    … …  
    7681        //now make references to data in newType folder 
    7782        DestPath="root:Packages:NIST:"+newType   
    78         WAVE data=$(destPath +":data")                  // these wave references point to the EXISTING work data 
     83        WAVE data=$(destPath +":linear_data")                   // these wave references point to the EXISTING work data 
     84        WAVE data_copy=$(destPath +":data")                     // these wave references point to the EXISTING work data 
     85        WAVE dest_data_err=$(destPath +":linear_data_error")                    // these wave references point to the EXISTING work data 
    7986        WAVE/T textread=$(destPath + ":textread") 
    8087        WAVE integersread=$(destPath + ":integersread") 
    … …  
    110117        //then unscale the data array 
    111118        data *= uscale 
     119        dest_data_err *= uscale 
    112120         
    113121        //DetCorr() has not been applied to the data in RAW , do it now in a local reference to the raw data 
    114         WAVE raw_data = $"root:Packages:NIST:RAW:data" 
     122        WAVE raw_data = $"root:Packages:NIST:RAW:linear_data" 
     123        WAVE raw_data_err = $"root:Packages:NIST:RAW:linear_data_error" 
    115124        WAVE raw_reals =  $"root:Packages:NIST:RAW:realsread" 
    116125        WAVE/T raw_text = $"root:Packages:NIST:RAW:textread" 
    … …  
    128137        endif    
    129138         
    130         DetCorr(raw_data,raw_reals,doEfficiency,doTrans)        //applies correction to raw_data, and overwrites it 
     139        DetCorr(raw_data,raw_data_err,raw_reals,doEfficiency,doTrans)   //applies correction to raw_data, and overwrites it 
    131140         
    132141        //if RAW data is ILL type detector, correct raw_data for same counts being written to 4 pixels 
    133142        if(cmpstr(raw_text[9], "ILL   ") == 0 )         //text field in header is 6 characters "ILL---" 
    134143                raw_data /= 4 
     144                raw_data_err /= 4 
    135145        endif 
    136146         
    … …  
    155165                dscale = 1/(1-deadTime*cntrate) 
    156166                // multiply data[ii][] by the dead time 
    157                 data[][ii] *= dscale 
     167                raw_data[][ii] *= dscale 
     168                raw_data_err[][ii] *= dscale 
    158169        endfor 
     170#else 
     171        // dead time correction on all other RAW data, including NCNR 
     172        raw_data *= dscale 
     173        raw_data_err *= dscale 
    159174#endif 
    160175 
    … …  
    162177        total_mon += raw_reals[0] 
    163178#if (exists("ILL_D22")==6) 
    164         total_det += sum(data,-inf,inf)                 //add the newly scaled detector array 
     179        total_det += sum(raw_data,-inf,inf)                     //add the newly scaled detector array 
    165180#else 
    166181        total_det += dscale*raw_reals[2] 
    … …  
    184199         
    185200        If((xshift == 0) && (yshift == 0))              //no shift, just add them 
    186                 data += dscale*raw_data         //do the deadtime correction on RAW here 
     201                data += raw_data                //deadtime correction has already been done to the raw data 
     202                dest_data_err = sqrt(dest_data_err^2 + raw_data_err^2)                  // error of the sum 
    187203        else 
    188204                //shift the beamcenter, then add 
    … …  
    202218                                else 
    203219                                        //add the raw_data + shifted sum (and do the deadtime correction on both) 
    204                                         data[ii][jj] += dscale*(raw_data[ii][jj]+sh_sum)                //do the deadtime correction on RAW here 
     220                                        data[ii][jj] += (raw_data[ii][jj]+sh_sum)               //do the deadtime correction on RAW here 
     221                                        dest_data_err[ii][jj] = sqrt(dest_data_err[ii][jj]^2 + raw_data_err[ii][jj]^2)                  // error of the sum 
    205222                                Endif 
    206223                                jj+=1 
    … …  
    213230        scale = defmon/total_mon 
    214231        data *= scale 
     232        dest_data_err *= scale 
     233         
     234        // keep "data" and linear_data in sync in the destination folder 
     235        data_copy = data 
    215236         
    216237        //all is done, except for the bookkeeping of updating the header info in the work folder 
    … …  
    277298        DestPath = "root:Packages:NIST:"+newType 
    278299        Duplicate/O $"root:Packages:NIST:RAW:data",$(destPath + ":data") 
     300        Duplicate/O $"root:Packages:NIST:RAW:linear_data_error",$(destPath + ":linear_data_error") 
     301        Duplicate/O $"root:Packages:NIST:RAW:linear_data",$(destPath + ":linear_data") 
    279302//      Duplicate/O $"root:Packages:NIST:RAW:vlegend",$(destPath + ":vlegend") 
    280303        Duplicate/O $"root:Packages:NIST:RAW:textread",$(destPath + ":textread") 
    … …  
    282305        Duplicate/O $"root:Packages:NIST:RAW:realsread",$(destPath + ":realsread") 
    283306        Variable/G $(destPath + ":gIsLogscale")=0                       //overwite flag in newType folder, data converted (above) to linear scale 
    284          
    285         WAVE data=$(destPath + ":data")                         // these wave references point to the data in work 
    286         WAVE/T textread=$(destPath + ":textread")                       //that are to be directly operated on 
     307 
     308        WAVE data=$(destPath + ":linear_data")                                                  // these wave references point to the data in work 
     309        WAVE data_copy=$(destPath + ":data")                                                    // these wave references point to the data in work 
     310        WAVE data_err=$(destPath + ":linear_data_error") 
     311        WAVE/T textread=$(destPath + ":textread")                               //that are to be directly operated on 
    287312        WAVE integersread=$(destPath + ":integersread") 
    288313        WAVE realsread=$(destPath + ":realsread") 
    … …  
    298323        endif 
    299324         
    300         DetCorr(data,realsread,doEfficiency,doTrans)            //the parameters are waves, and will be changed by the function 
    301          
     325        DetCorr(data,data_err,realsread,doEfficiency,doTrans)           //the parameters are waves, and will be changed by the function 
     326 
    302327        //if ILL type detector, correct for same counts being written to 4 pixels 
    303328        if(cmpstr(textread[9], "ILL   ") == 0 )         //text field in header is 6 characters "ILL---" 
    304329                data /= 4 
     330                data_err /= 4           //rescale error 
    305331        endif 
    306332         
    … …  
    326352                // multiply data[ii][] by the dead time 
    327353                data[][ii] *= dscale 
     354                data_err[][ii] *= dscale 
    328355        endfor 
    329356#endif 
    330          
     357 
     358        // NO xcenter,ycenter shifting is done - this is the first (and only) file in the work folder 
     359         
     360        //only ONE data file- no addition of multiple runs in this function, so data is 
     361        //just simply corrected for deadtime. 
     362#if (exists("ILL_D22")==0)              //ILL correction done tube-by-tube above 
     363        data *= dscale          //deadtime correction for everyone else, including NCNR 
     364        data_err *= dscale 
     365#endif 
    331366         
    332367        //update totals to put in the work header (at the end of the function) 
    … …  
    341376        total_numruns +=1 
    342377         
    343         // NO xcenter,ycenter shifting is done - this is the first (and only) file in the work folder 
    344          
    345         //only ONE data file- no addition of multiple runs in this function, so data is 
    346         //just simply corrected for deadtime. 
    347 #ifndef ILL_D22         //correction done tube-by-tube above 
    348         data *= dscale          //deadtime correction 
    349 #endif 
    350  
     378         
    351379        //scale the data to the default montor counts 
    352380        scale = defmon/total_mon 
    353381        data *= scale 
     382        data_err *= scale               //assumes total monitor count is so large there is essentially no error 
     383         
     384        // to keep "data" and linear_data in sync 
     385        data_copy = data 
    354386         
    355387        //all is done, except for the bookkeeping, updating the header information in the work folder 
    … …  
    383415        Raw_to_work(type) 
    384416        //data is now in "type" folder 
    385         WAVE data=$("root:Packages:NIST:"+type+":data") 
     417        WAVE data=$("root:Packages:NIST:"+type+":linear_data") 
     418        WAVE data_copy=$("root:Packages:NIST:"+type+":data") 
     419        WAVE data_err=$("root:Packages:NIST:"+type+":linear_data_error") 
    386420        WAVE new_reals=$("root:Packages:NIST:"+type+":realsread") 
    387421         
    … …  
    393427         
    394428        data /= scale           //unscale the data 
     429        data_err /= scale 
     430         
     431        // to keep "data" and linear_data in sync 
     432        data_copy = data 
    395433         
    396434        return(0) 
    … …  
    409447        Add_Raw_to_work(type) 
    410448        //data is now in "type" folder 
    411         WAVE data=$("root:Packages:NIST:"+type+":data") 
     449        WAVE data=$("root:Packages:NIST:"+type+":linear_data") 
     450        WAVE data_copy=$("root:Packages:NIST:"+type+":data") 
     451        WAVE data_err=$("root:Packages:NIST:"+type+":linear_data_error") 
    412452        WAVE new_reals=$("root:Packages:NIST:"+type+":realsread") 
    413453         
    … …  
    419459         
    420460        data /= scale           //unscale the data 
     461        data_err /= scale 
     462         
     463        // to keep "data" and linear_data in sync 
     464        data_copy = data 
    421465         
    422466        return(0) 
    … …  
    427471//works on the actual data array, assumes that is is already on LINEAR scale 
    428472// 
    429 Function DetCorr(data,realsread,doEfficiency,doTrans) 
    430         Wave data,realsread 
     473Function DetCorr(data,data_err,realsread,doEfficiency,doTrans) 
     474        Wave data,data_err,realsread 
    431475        Variable doEfficiency,doTrans 
    432476         
    … …  
    434478        Variable ii,jj,dtdist,dtdis2 
    435479        Variable xi,xd,yd,rad,ratio,domega,xy 
    436         Variable lambda,trans 
     480        Variable lambda,trans,trans_err,lat_err,tmp_err,lat_corr 
    437481         
    438482//      Print "...doing jacobian and non-linear corrections" 
    … …  
    457501        lambda = realsRead[26] 
    458502        trans = RealsRead[4] 
     503        trans_err = RealsRead[41]               //new, March 2011 
    459504         
    460505        xx0 = dc_fx(x0,sx,sx3,xcenter) 
    … …  
    492537                        data[ii][jj] *= xy*ratio 
    493538                        solidAngle[ii][jj] = xy*ratio           //testing only   
    494  
     539                        data_err[ii][jj] *= xy*ratio                    //error propagation assumes that SA and Jacobian are exact, so simply scale error 
     540                         
    495541                         
    496542                        // correction factor for detector efficiency JBG memo det_eff_cor2.doc 3/20/07 
    … …  
    501547#if (exists("ILL_D22")==6) 
    502548                                data[ii][jj]  /= DetEffCorrILL(lambda,dtdist,xd)                //tube-by-tube corrections  
    503                   solidAngle[ii][jj] = DetEffCorrILL(lambda,dtdist,xd) 
     549                                data_err[ii][jj] /= DetEffCorrILL(lambda,dtdist,xd)                     //assumes correction is exact 
     550//                solidAngle[ii][jj] = DetEffCorrILL(lambda,dtdist,xd) 
    504551#else 
    505552                                data[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd) 
     553                                data_err[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd) 
    506554//                              solidAngle[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd)           //testing only 
    507555#endif 
    … …  
    522570                                        trans = 1 
    523571                                endif 
    524                                          
    525                                 data[ii][jj] /= LargeAngleTransmissionCorr(trans,dtdist,xd,yd)          //moved from 1D avg SRK 11/2007 
     572                                 
     573                                // pass in the transmission error, and the error in the correction is returned as the last parameter 
     574                                lat_corr = LargeAngleTransmissionCorr(trans,dtdist,xd,yd,trans_err,lat_err)             //moved from 1D avg SRK 11/2007 
     575                                data[ii][jj] /= lat_corr                        //divide by the correction factor 
     576                                // 
     577                                // 
     578                                // 
     579                                // relative errors add in quadrature 
     580                                tmp_err = (data_err[ii][jj]/lat_corr)^2 + (lat_err/lat_corr)^2*data[ii][jj]*data[ii][jj]/lat_corr^2 
     581                                tmp_err = sqrt(tmp_err) 
     582                                 
     583                                data_err[ii][jj] = tmp_err 
     584                                 
     585                                solidAngle[ii][jj] = lat_err 
     586 
     587                                 
    526588                                //solidAngle[ii][jj] = LargeAngleTransmissionCorr(trans,dtdist,xd,yd)           //testing only 
    527589                        endif 
    … …  
    599661 
    600662// DIVIDE the intensity by this correction to get the right answer 
    601 Function LargeAngleTransmissionCorr(trans,dtdist,xd,yd) 
    602         Variable trans,dtdist,xd,yd 
     663Function LargeAngleTransmissionCorr(trans,dtdist,xd,yd,trans_err,err) 
     664        Variable trans,dtdist,xd,yd,trans_err,&err 
    603665 
    604666        //angle dependent transmission correction  
    … …  
    618680        //optical thickness 
    619681        uval = -ln(trans)               //use natural logarithm 
    620  
    621 //      theta = 2*asin(lambda*qval/(4*pi)) 
    622  
    623682        cos_th = cos(theta) 
    624683        arg = (1-cos_th)/cos_th 
    625         if((uval<0.01) || (cos_th>0.99))                //OR 
     684         
     685        // a Taylor series around uval*arg=0 only needs about 4 terms for very good accuracy 
     686        //                      correction= 1 - 0.5*uval*arg + (uval*arg)^2/6 - (uval*arg)^3/24 + (uval*arg)^4/120 
     687        // OR 
     688        if((uval<0.01) || (cos_th>0.99))         
    626689                //small arg, approx correction 
    627690                correction= 1-0.5*uval*arg 
    … …  
    631694        endif 
    632695 
     696        Variable tmp 
     697         
     698        if(trans == 1) 
     699                err = 0         //no correction, no error 
     700        else 
     701                //sigT, calculated from the Taylor expansion 
     702                tmp = (1/trans)*(arg/2-arg^2/3*uval+arg^3/8*uval^2-arg^4/30*uval^3) 
     703                tmp *= tmp 
     704                tmp *= trans_err^2 
     705                tmp = sqrt(tmp)         //sigT 
     706                 
     707                err = tmp 
     708        endif 
     709         
     710//      Printf "trans error = %g\r",trans_err 
     711//      Printf "correction = %g +/- %g\r", correction, err 
     712         
    633713        //end of transmission/pathlength correction 
    634714 
    … …  
    806886        // if the desired workfile doesn't exist, let the user know, and abort 
    807887        String destPath="" 
    808         destPath = "root:Packages:NIST:"+Type + ":data" 
    809         if(WaveExists($destpath) == 0) 
     888 
     889        if(WaveExists($("root:Packages:NIST:"+Type + ":data")) == 0) 
    810890                Print "There is no work file in "+type+"--Aborting" 
    811891                Return(1)               //error condition 
    … …  
    813893        //check for DIV 
    814894        // if the DIV workfile doesn't exist, let the user know,and abort 
    815         destPath = "root:Packages:NIST:DIV:data" 
    816         if(WaveExists($destpath) == 0) 
     895 
     896        if(WaveExists($"root:Packages:NIST:DIV:data") == 0) 
    817897                Print "There is no work file in DIV --Aborting" 
    818898                Return(1)               //error condition 
    … …  
    836916        destPath = "root:Packages:NIST:" + type 
    837917        Duplicate/O $(destPath + ":data"),$"root:Packages:NIST:CAL:data" 
     918        Duplicate/O $(destPath + ":linear_data"),$"root:Packages:NIST:CAL:linear_data" 
     919        Duplicate/O $(destPath + ":linear_data_error"),$"root:Packages:NIST:CAL:linear_data_error" 
    838920//      Duplicate/O $(destPath + ":vlegend"),$"root:Packages:NIST:CAL:vlegend" 
    839921        Duplicate/O $(destPath + ":textread"),$"root:Packages:NIST:CAL:textread" 
    … …  
    846928        destPath = "root:Packages:NIST:CAL" 
    847929        //make appropriate wave references 
    848         Wave data=$(destPath + ":data")                                 // these wave references point to the data in CAL 
     930        Wave data=$(destPath + ":linear_data")                                  // these wave references point to the data in CAL 
     931//      Wave data_err=$(destPath + ":linear_data_err")                                  // these wave references point to the data in CAL 
     932        Wave data_copy=$(destPath + ":data")                                    // these wave references point to the data in CAL 
    849933        Wave/t textread=$(destPath + ":textread")                       //that are to be directly operated on 
    850934        Wave integersread=$(destPath + ":integersread") 
    … …  
    857941        //do the division, changing data in CAL 
    858942        data /= div_data 
     943         
     944//      data_err /= div_data 
     945         
     946        // keep "data" in sync with linear_data 
     947        data_copy = data 
    859948         
    860949        //update CAL header 
    … …  
    904993//w_ is the "work" file 
    905994//both are work files and should already be normalized to 10^8 monitor counts 
    906 Function Absolute_Scale(type,w_trans,w_thick,s_trans,s_thick,s_izero,s_cross) 
     995Function Absolute_Scale(type,w_trans,w_thick,s_trans,s_thick,s_izero,s_cross,kappa_err) 
    907996        String type 
    908         Variable w_trans,w_thick,s_trans,s_thick,s_izero,s_cross 
     997        Variable w_trans,w_thick,s_trans,s_thick,s_izero,s_cross,kappa_err 
    909998         
    910999        //convert the "type" data to absolute scale using the given standard information 
    … …  
    9141003        String destPath 
    9151004        //check for "type" 
    916         destPath = "root:Packages:NIST:"+Type + ":data" 
    917         if(WaveExists($destpath) == 0) 
     1005        if(WaveExists($("root:Packages:NIST:"+Type + ":data")) == 0) 
    9181006                Print "There is no work file in "+type+"--Aborting" 
    9191007                Return(1)               //error condition 
    … …  
    9341022        //copy from current dir (type) to ABS, defined by destPath 
    9351023        Duplicate/O $(oldType + ":data"),$"root:Packages:NIST:ABS:data" 
     1024        Duplicate/O $(oldType + ":linear_data"),$"root:Packages:NIST:ABS:linear_data" 
     1025        Duplicate/O $(oldType + ":linear_data_error"),$"root:Packages:NIST:ABS:linear_data_error" 
    9361026//      Duplicate/O $(oldType + ":vlegend"),$"root:Packages:NIST:ABS:vlegend" 
    9371027        Duplicate/O $(oldType + ":textread"),$"root:Packages:NIST:ABS:textread" 
    … …  
    9451035        //now switch to ABS folder 
    9461036        //make appropriate wave references 
    947         WAVE data=$"root:Packages:NIST:ABS:data"                                        // these wave references point to the "type" data in ABS 
     1037        WAVE data=$"root:Packages:NIST:ABS:linear_data"                                 // these wave references point to the "type" data in ABS 
     1038        WAVE data_err=$"root:Packages:NIST:ABS:linear_data_error"                                       // these wave references point to the "type" data in ABS 
     1039        WAVE data_copy=$"root:Packages:NIST:ABS:data"                                   // just for display 
    9481040        WAVE/T textread=$"root:Packages:NIST:ABS:textread"                      //that are to be directly operated on 
    9491041        WAVE integersread=$"root:Packages:NIST:ABS:integersread" 
    … …  
    9621054         
    9631055        //calculate scale factor 
     1056        Variable scale,trans_err 
    9641057        s1 = defmon/realsread[0]                //[0] is monitor count (s1 should be 1) 
    9651058        s2 = s_thick/w_thick 
    … …  
    9671060        s4 = s_cross/s_izero 
    9681061         
     1062        // kappa comes in as s_izero, so be sure to use 1/kappa_err 
     1063         
    9691064        data *= s1*s2*s3*s4 
     1065         
     1066        scale = s1*s2*s3*s4 
     1067        trans_err = realsRead[41] 
     1068         
     1069//      print scale 
     1070//      print data[0][0] 
     1071         
     1072        data_err = sqrt(scale^2*data_err^2 + scale^2*data^2*(kappa_err^2/s_izero^2 +trans_err^2/w_trans^2)) 
     1073 
     1074//      print data_err[0][0] 
     1075         
     1076// keep "data" in sync with linear_data  
     1077        data_copy = data 
    9701078         
    9711079        //********* 15APR02 
    … …  
    9961104        // if the desired workfile doesn't exist, let the user know, and abort 
    9971105        String destPath="" 
    998         destPath = "root:Packages:NIST:"+oldType + ":data" 
    999         if(WaveExists($destpath) == 0) 
     1106        if(WaveExists($("root:Packages:NIST:"+oldType + ":data")) == 0) 
    10001107                Print "There is no work file in "+oldtype+"--Aborting" 
    10011108                Return(1)               //error condition 
    … …  
    10101117        destPath = "root:Packages:NIST:" + oldtype 
    10111118        Duplicate/O $(destPath + ":data"),$("root:Packages:NIST:"+newtype+":data") 
     1119        Duplicate/O $(destPath + ":linear_data"),$("root:Packages:NIST:"+newtype+":linear_data") 
     1120        Duplicate/O $(destPath + ":linear_data_error"),$("root:Packages:NIST:"+newtype+":linear_data_error") 
    10121121        Duplicate/O $(destPath + ":textread"),$("root:Packages:NIST:"+newtype+":textread") 
    10131122        Duplicate/O $(destPath + ":integersread"),$("root:Packages:NIST:"+newtype+":integersread") 
    … …  
    10151124        // 
    10161125        // be sure to get rid of the linear_data if it exists in the destination folder 
    1017         KillWaves/Z $("root:Packages:NIST:"+newtype+":linear_data") 
     1126//      KillWaves/Z $("root:Packages:NIST:"+newtype+":linear_data") 
     1127 
    10181128        //need to save a copy of filelist string too (from the current type folder) 
    10191129        SVAR oldFileList = $(destPath + ":fileList") 
    … …  
    11081218         
    11091219        WAVE/Z data1=$("root:Packages:NIST:"+workMathStr+"File_1:data") 
     1220        WAVE/Z err1=$("root:Packages:NIST:"+workMathStr+"File_1:linear_data_error") 
     1221         
     1222        // set # 2 
    11101223        If(cmpstr(str2,"UNIT MATRIX")==0) 
    11111224                Make/O/N=(pixelsX,pixelsY) root:Packages:NIST:WorkMath:data             //don't put in File_2 folder 
    11121225                Wave/Z data2 =  root:Packages:NIST:WorkMath:data                        //it's not real data! 
    11131226                data2=1 
     1227                Duplicate/O data2 err2 
     1228                err2 = 0 
    11141229        else 
    11151230                //Load set #2 
    11161231                Load_NamedASC_File(pathStr+str2,workMathStr+"File_2") 
    11171232                WAVE/Z data2=$("root:Packages:NIST:"+workMathStr+"File_2:data") 
     1233                WAVE/Z err2=$("root:Packages:NIST:"+workMathStr+"File_2:linear_data_error") 
    11181234        Endif 
    11191235 
    … …  
    11281244        CopyWorkContents(workMathStr+"File_1",workMathStr+dest) 
    11291245        WAVE/Z destData=$("root:Packages:NIST:"+workMathStr+dest+":data") 
     1246        WAVE/Z destErr=$("root:Packages:NIST:"+workMathStr+dest+":linear_data_error") 
    11301247         
    11311248        //dispatch 
    … …  
    11331250                case "*":               //multiplication 
    11341251                        destData = const1*data1 * const2*data2 
     1252                        destErr = const1^2*const2^2*(err1^2*data2^2 + err2^2*data1^2) 
     1253                        destErr = sqrt(destErr) 
    11351254                        break    
    11361255                case "_":               //subtraction 
    11371256                        destData = const1*data1 - const2*data2 
     1257                        destErr = const1^2*err1^2 + const2^2*err2^2 
     1258                        destErr = sqrt(destErr) 
    11381259                        break 
    11391260                case "/":               //division 
    11401261                        destData = (const1*data1) / (const2*data2) 
     1262                        destErr = const1^2/const2^2*(err1^2/data2^2 + err2^2*data1^2/data2^4) 
     1263                        destErr = sqrt(destErr) 
    11411264                        break 
    11421265                case "+":               //addition 
    11431266                        destData = const1*data1 + const2*data2 
     1267                        destErr = const1^2*err1^2 + const2^2*err2^2 
     1268                        destErr = sqrt(destErr) 
    11441269                        break                    
    11451270        endswitch 

  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/WriteQIS.ipf

    r762 r794  
    419419        //labelWave[19] = "ASCII data created " +date()+" "+time() 
    420420        PathInfo catPathName 
    421         String sfPath = S_Path+StringFromList(0,samfiles,";") 
     421        String sfPath = S_Path+RemoveAllSpaces(StringFromList(0,samfiles,";"))          //make sure there are no leading spaces in the file name 
    422422        print sfPath 
    423423        labelWave[19] = "RAW SAM FILE "+StringFromList(0, samfiles  , ";")+ " TIMESTAMP: "+getFileCreationDate(sfPath) 
    … …  
    675675         
    676676        Wave data=$(destStr+typeStr) 
     677        Wave data_err=$(destStr+":linear_data_error") 
    677678        WAVE intw=$(destStr + ":integersRead") 
    678679        WAVE rw=$(destStr + ":realsRead") 
    … …  
    691692        If(!(WaveExists(data))) 
    692693                Abort "data DNExist QxQy_Export()" 
     694        Endif 
     695        If(!(WaveExists(data_err))) 
     696                Abort "linear data error DNExist QxQy_Export()" 
    693697        Endif 
    694698        If(!(WaveExists(intw))) 
    … …  
    813817//********************* 
    814818 
    815         // generate my own error wave for I(qx,qy) 
    816         Duplicate/O z_val sw 
    817         sw = sqrt(z_val)                //assumes Poisson statistics for each cell (counter) 
    818         //      sw = 0.05*sw            // uniform 5% error? tends to favor the low intensity too strongly 
    819         // get rid of the "bad" errorsby replacing the NaN, Inf, and zero with V_avg 
    820         // THIS IS EXTREMEMLY IMPORTANT - if this is not done, there are some "bad" values in the  
    821         // error wave (things that are not numbers) - and this wrecks the smeared model fitting. 
    822         // It appears to have no effect on the unsmeared model. 
    823         WaveStats/Q sw 
    824         sw = numtype(sw[p]) == 0 ? sw[p] : V_avg 
    825         sw = sw[p] != 0 ? sw[p] : V_avg 
    826          
     819//      // generate my own error wave for I(qx,qy) 
     820//      Duplicate/O z_val sw 
     821//      sw = sqrt(z_val)                //assumes Poisson statistics for each cell (counter) 
     822//      //      sw = 0.05*sw            // uniform 5% error? tends to favor the low intensity too strongly 
     823//      // get rid of the "bad" errorsby replacing the NaN, Inf, and zero with V_avg 
     824//      // THIS IS EXTREMEMLY IMPORTANT - if this is not done, there are some "bad" values in the  
     825//      // error wave (things that are not numbers) - and this wrecks the smeared model fitting. 
     826//      // It appears to have no effect on the unsmeared model. 
     827//      WaveStats/Q sw 
     828//      sw = numtype(sw[p]) == 0 ? sw[p] : V_avg 
     829//      sw = sw[p] != 0 ? sw[p] : V_avg 
     830         
     831        // now use the properly propagated 2D error 
     832        Duplicate/O data_err sw 
     833        Redimension/N=(pixelsX*pixelsY) sw 
     834 
     835 
    827836 
    828837        //not demo-compatible, but approx 8x faster!!    
Note: See TracChangeset for help on using the changeset viewer.