Changeset 794 for sans/Dev/trunk


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
Files:
18 edited

Legend:

Unmodified
Added
Removed

  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/Models_2D/CoreShellCyl2D_v40.ipf

    r771 r794  
    276276        Make/O/D/N=15 CSCyl2D_tmp 
    277277        CSCyl2D_tmp[0,13] = cw 
    278         CSCyl2D_tmp[14] = 25 
     278        CSCyl2D_tmp[14] = 5             // hard-wire the number of integration points, small number since smearing is used 
    279279        CSCyl2D_tmp[7] = 0              //send a zero background to the calculation, add it in later 
    280280 
    … …  
    383383         
    384384        Variable spl,spp,apl,app,bpl,bpp,phi_pt,qpl_pt 
     385        Variable qperp_pt,phi_prime,q_prime 
    385386 
    386387        //loop over q-values 
    … …  
    401402                apl = -numStdDev*spl + qval             //parallel = q integration limits 
    402403                bpl = numStdDev*spl + qval 
    403                 app = -numStdDev*spp + phi              //perpendicular = phi integration limits 
    404                 bpp = numStdDev*spp + phi 
     404                app = -numStdDev*spp + 0                //q_perp = 0 
     405                bpp = numStdDev*spp + 0 
    405406                 
    406407                //make sure the limits are reasonable. 
    … …  
    413414                sumOut = 0 
    414415                for(jj=0;jj<nord;jj+=1)         // call phi the "outer' 
    415                         phi_pt = (xi[jj]*(bpp-app)+app+bpp)/2 
     416                        qperp_pt = (xi[jj]*(bpp-app)+app+bpp)/2         //this is now q_perp 
    416417 
    417418                        sumIn=0 
    … …  
    420421                                qpl_pt = (xi[kk]*(bpl-apl)+apl+bpl)/2 
    421422                                 
    422                                 FindQxQy(qpl_pt,phi_pt,qx_pt,qy_pt)             //find the corresponding QxQy to the Q,phi 
     423                                // find QxQy given Qpl and Qperp on the grid 
     424                                // 
     425                                q_prime = sqrt(qpl_pt^2+qperp_pt^2) 
     426                                phi_prime = phi + qperp_pt/qpl_pt 
     427                                FindQxQy(q_prime,phi_prime,qx_pt,qy_pt) 
     428                                 
    423429                                yPtw[kk] = qy_pt                                        //phi is the same in this loop, but qy is not 
    424430                                xPtW[kk] = qx_pt                                        //qx is different here too, as we're varying Qpl 
    425431                                 
    426                                 res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (phi_pt-phi)^2/spp/spp ) ) 
     432                                res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (qperp_pt)^2/spp/spp ) ) 
    427433                                res_tot[kk] /= normFactor 
    428434//                              res_tot[kk] *= fs 

  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/Models_2D/Cylinder_2D_v40.ipf

    r788 r794  
    221221/// now using the MultiThread keyword. as of Igor 6.20, the manual threading 
    222222// as above gives a wave read error (index out of range). Same code works fine in Igor 6.12 
     223//ThreadSafe Function Cylinder2D(cw,zw,xw,yw) : FitFunc 
    223224Function Cylinder2D(cw,zw,xw,yw) : FitFunc 
    224225        Wave cw,zw,xw,yw 
    … …  
    252253//      Smear_2DModel_PP(Cylinder2D_noThread,s,10) 
    253254 
     255// this is generic, but I need to declare the Cylinder2D threadsafe 
     256// and this calculation is significantly slower than the manually threaded calculation 
     257// if the function is fast to calculate. Once the function has polydispersity on 2 or more parameters 
     258// then this AAO calculation and the manual threading are both painfully slow, and more similar in execution time 
     259// 
     260// see the function for timing details 
     261// 
     262//      Smear_2DModel_PP_AAO(Cylinder2D,s,10) 
    254263 
    255264//// the last param is nord 
    … …  
    385394 
    386395// 
    387 // - worker function for threads of Sphere2D 
     396// - worker function for threads of Cylinder2D 
    388397// 
    389398ThreadSafe Function SmearedCylinder2D_T(coef,qxw,qyw,qzw,sxw,syw,fsw,zw,wt,xi,pt1,pt2,nord) 
    … …  
    409418         
    410419        Variable spl,spp,apl,app,bpl,bpp,phi_pt,qpl_pt 
     420        Variable qperp_pt,phi_prime,q_prime 
    411421 
    412422        //loop over q-values 
    … …  
    427437                apl = -numStdDev*spl + qval             //parallel = q integration limits 
    428438                bpl = numStdDev*spl + qval 
    429                 app = -numStdDev*spp + phi              //perpendicular = phi integration limits 
    430                 bpp = numStdDev*spp + phi 
     439                app = -numStdDev*spp + 0                //q_perp = 0 
     440                bpp = numStdDev*spp + 0 
    431441                 
    432442                //make sure the limits are reasonable. 
    … …  
    439449                sumOut = 0 
    440450                for(jj=0;jj<nord;jj+=1)         // call phi the "outer' 
    441                         phi_pt = (xi[jj]*(bpp-app)+app+bpp)/2 
     451                        qperp_pt = (xi[jj]*(bpp-app)+app+bpp)/2         //this is now q_perp 
    442452 
    443453                        sumIn=0 
    … …  
    446456                                qpl_pt = (xi[kk]*(bpl-apl)+apl+bpl)/2 
    447457                                 
    448                                 FindQxQy(qpl_pt,phi_pt,qx_pt,qy_pt)             //find the corresponding QxQy to the Q,phi 
     458                                // find QxQy given Qpl and Qperp on the grid 
     459                                // 
     460                                q_prime = sqrt(qpl_pt^2+qperp_pt^2) 
     461                                phi_prime = phi + qperp_pt/qpl_pt 
     462                                FindQxQy(q_prime,phi_prime,qx_pt,qy_pt) 
     463                                 
    449464                                yPtw[kk] = qy_pt                                        //phi is the same in this loop, but qy is not 
    450465                                xPtW[kk] = qx_pt                                        //qx is different here too, as we're varying Qpl 
    451466                                 
    452                                 res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (phi_pt-phi)^2/spp/spp ) ) 
     467                                res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (qperp_pt)^2/spp/spp ) ) 
    453468                                res_tot[kk] /= normFactor 
    454469//                              res_tot[kk] *= fs 

  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/Models_2D/Ellipsoid2D_v40.ipf

    r771 r794  
    270270        Make/O/D/N=13 Ellip2D_tmp 
    271271        Ellip2D_tmp[0,11] = cw 
    272         Ellip2D_tmp[12] = 25 
     272        Ellip2D_tmp[12] = 5                     //use a small number of integration points since smearing is used 
    273273        Ellip2D_tmp[5] = 0              //pass in a zero background and add it in later 
    274274         
    … …  
    377377         
    378378        Variable spl,spp,apl,app,bpl,bpp,phi_pt,qpl_pt 
     379        Variable qperp_pt,phi_prime,q_prime 
    379380 
    380381        //loop over q-values 
    … …  
    395396                apl = -numStdDev*spl + qval             //parallel = q integration limits 
    396397                bpl = numStdDev*spl + qval 
    397                 app = -numStdDev*spp + phi              //perpendicular = phi integration limits 
    398                 bpp = numStdDev*spp + phi 
     398                app = -numStdDev*spp + 0                //q_perp = 0 
     399                bpp = numStdDev*spp + 0 
    399400                 
    400401                //make sure the limits are reasonable. 
    … …  
    407408                sumOut = 0 
    408409                for(jj=0;jj<nord;jj+=1)         // call phi the "outer' 
    409                         phi_pt = (xi[jj]*(bpp-app)+app+bpp)/2 
     410                        qperp_pt = (xi[jj]*(bpp-app)+app+bpp)/2         //this is now q_perp 
    410411 
    411412                        sumIn=0 
    … …  
    414415                                qpl_pt = (xi[kk]*(bpl-apl)+apl+bpl)/2 
    415416                                 
    416                                 FindQxQy(qpl_pt,phi_pt,qx_pt,qy_pt)             //find the corresponding QxQy to the Q,phi 
     417                                // find QxQy given Qpl and Qperp on the grid 
     418                                // 
     419                                q_prime = sqrt(qpl_pt^2+qperp_pt^2) 
     420                                phi_prime = phi + qperp_pt/qpl_pt 
     421                                FindQxQy(q_prime,phi_prime,qx_pt,qy_pt) 
     422                                 
    417423                                yPtw[kk] = qy_pt                                        //phi is the same in this loop, but qy is not 
    418424                                xPtW[kk] = qx_pt                                        //qx is different here too, as we're varying Qpl 
    419425                                 
    420                                 res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (phi_pt-phi)^2/spp/spp ) ) 
     426                                res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (qperp_pt)^2/spp/spp ) ) 
    421427                                res_tot[kk] /= normFactor 
    422428//                              res_tot[kk] *= fs 

  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/Models_2D/EllipticalCylinder2D_v40.ipf

    r771 r794  
    273273        Make/O/D/N=15 EllCyl2D_tmp 
    274274        EllCyl2D_tmp[0,13] = cw 
    275         EllCyl2D_tmp[14] = 25 
     275        EllCyl2D_tmp[14] = 5            // small number of integration points since smearing is used 
    276276        EllCyl2D_tmp[6] = 0             //pass in a zero background and add it in later 
    277277         
    … …  
    380380         
    381381        Variable spl,spp,apl,app,bpl,bpp,phi_pt,qpl_pt 
     382        Variable qperp_pt,phi_prime,q_prime 
    382383 
    383384        //loop over q-values 
    … …  
    398399                apl = -numStdDev*spl + qval             //parallel = q integration limits 
    399400                bpl = numStdDev*spl + qval 
    400                 app = -numStdDev*spp + phi              //perpendicular = phi integration limits 
    401                 bpp = numStdDev*spp + phi 
     401                app = -numStdDev*spp + 0                //q_perp = 0 
     402                bpp = numStdDev*spp + 0 
    402403                 
    403404                //make sure the limits are reasonable. 
    … …  
    410411                sumOut = 0 
    411412                for(jj=0;jj<nord;jj+=1)         // call phi the "outer' 
    412                         phi_pt = (xi[jj]*(bpp-app)+app+bpp)/2 
     413                        qperp_pt = (xi[jj]*(bpp-app)+app+bpp)/2         //this is now q_perp 
    413414 
    414415                        sumIn=0 
    … …  
    417418                                qpl_pt = (xi[kk]*(bpl-apl)+apl+bpl)/2 
    418419                                 
    419                                 FindQxQy(qpl_pt,phi_pt,qx_pt,qy_pt)             //find the corresponding QxQy to the Q,phi 
     420                                // find QxQy given Qpl and Qperp on the grid 
     421                                // 
     422                                q_prime = sqrt(qpl_pt^2+qperp_pt^2) 
     423                                phi_prime = phi + qperp_pt/qpl_pt 
     424                                FindQxQy(q_prime,phi_prime,qx_pt,qy_pt) 
     425                                 
    420426                                yPtw[kk] = qy_pt                                        //phi is the same in this loop, but qy is not 
    421427                                xPtW[kk] = qx_pt                                        //qx is different here too, as we're varying Qpl 
    422428                                 
    423                                 res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (phi_pt-phi)^2/spp/spp ) ) 
     429                                res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (qperp_pt)^2/spp/spp ) ) 
    424430                                res_tot[kk] /= normFactor 
    425431//                              res_tot[kk] *= fs 

  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/Models_2D/PeakGauss2D_v40.ipf

    r775 r794  
    328328         
    329329        Variable spl,spp,apl,app,bpl,bpp,phi_pt,qpl_pt 
     330        Variable qperp_pt,phi_prime,q_prime 
    330331 
    331332        //loop over q-values 
    … …  
    346347                apl = -numStdDev*spl + qval             //parallel = q integration limits 
    347348                bpl = numStdDev*spl + qval 
    348                 app = -numStdDev*spp + phi              //perpendicular = phi integration limits 
    349                 bpp = numStdDev*spp + phi 
     349                app = -numStdDev*spp + 0                //q_perp = 0 
     350                bpp = numStdDev*spp + 0 
    350351                 
    351352                //make sure the limits are reasonable. 
    … …  
    358359                sumOut = 0 
    359360                for(jj=0;jj<nord;jj+=1)         // call phi the "outer' 
    360                         phi_pt = (xi[jj]*(bpp-app)+app+bpp)/2 
     361                        qperp_pt = (xi[jj]*(bpp-app)+app+bpp)/2         //this is now q_perp 
    361362 
    362363                        sumIn=0 
    … …  
    365366                                qpl_pt = (xi[kk]*(bpl-apl)+apl+bpl)/2 
    366367                                 
    367                                 FindQxQy(qpl_pt,phi_pt,qx_pt,qy_pt)             //find the corresponding QxQy to the Q,phi 
     368                                // find QxQy given Qpl and Qperp on the grid 
     369                                // 
     370                                q_prime = sqrt(qpl_pt^2+qperp_pt^2) 
     371                                phi_prime = phi + qperp_pt/qpl_pt 
     372                                FindQxQy(q_prime,phi_prime,qx_pt,qy_pt) 
     373                                 
    368374                                yPtw[kk] = qy_pt                                        //phi is the same in this loop, but qy is not 
    369375                                xPtW[kk] = qx_pt                                        //qx is different here too, as we're varying Qpl 
    370376                                 
    371                                 res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (phi_pt-phi)^2/spp/spp ) ) 
     377                                res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (qperp_pt)^2/spp/spp ) ) 
    372378                                res_tot[kk] /= normFactor 
    373379//                              res_tot[kk] *= fs 

  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/Models_2D/Sphere_2D_v40.ipf

    r771 r794  
    347347         
    348348        Variable spl,spp,apl,app,bpl,bpp,phi_pt,qpl_pt 
     349        Variable qperp_pt,phi_prime,q_prime 
    349350 
    350351        //loop over q-values 
    … …  
    358359                spp = syw[ii] 
    359360                fs = fsw[ii] 
     361 
    360362                 
    361363                normFactor = 2*pi*spl*spp 
    … …  
    365367                apl = -numStdDev*spl + qval             //parallel = q integration limits 
    366368                bpl = numStdDev*spl + qval 
    367                 app = -numStdDev*spp + phi              //perpendicular = phi integration limits 
    368                 bpp = numStdDev*spp + phi 
     369///             app = -numStdDev*spp + phi              //perpendicular = phi integration limits (WRONG) 
     370///             bpp = numStdDev*spp + phi 
     371                app = -numStdDev*spp + 0                //q_perp = 0 
     372                bpp = numStdDev*spp + 0 
    369373                 
    370374                //make sure the limits are reasonable. 
    … …  
    377381                sumOut = 0 
    378382                for(jj=0;jj<nord;jj+=1)         // call phi the "outer' 
    379                         phi_pt = (xi[jj]*(bpp-app)+app+bpp)/2 
    380  
     383///                     phi_pt = (xi[jj]*(bpp-app)+app+bpp)/2 
     384                        qperp_pt = (xi[jj]*(bpp-app)+app+bpp)/2         //this is now q_perp 
     385                         
    381386                        sumIn=0 
    382387                        for(kk=0;kk<nord;kk+=1)         //at phi, integrate over Qpl 
    … …  
    384389                                qpl_pt = (xi[kk]*(bpl-apl)+apl+bpl)/2 
    385390                                 
    386                                 FindQxQy(qpl_pt,phi_pt,qx_pt,qy_pt)             //find the corresponding QxQy to the Q,phi 
     391///                             FindQxQy(qpl_pt,phi_pt,qx_pt,qy_pt)             //find the corresponding QxQy to the Q,phi 
     392 
     393                                // find QxQy given Qpl and Qperp on the grid 
     394                                // 
     395                                q_prime = sqrt(qpl_pt^2+qperp_pt^2) 
     396                                phi_prime = phi + qperp_pt/q_prime 
     397                                FindQxQy(q_prime,phi_prime,qx_pt,qy_pt) 
     398                                 
    387399                                yPtw[kk] = qy_pt                                        //phi is the same in this loop, but qy is not 
    388400                                xPtW[kk] = qx_pt                                        //qx is different here too, as we're varying Qpl 
    389                                  
    390                                 res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (phi_pt-phi)^2/spp/spp ) ) 
     401 
     402                                res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (qperp_pt)^2/spp/spp ) ) 
     403///                             res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (phi_pt-phi)^2/spp/spp ) ) 
    391404                                res_tot[kk] /= normFactor 
    392405//                              res_tot[kk] *= fs 

  • sans/Dev/trunk/NCNR_User_Procedures/Common/GaussUtils_v40.ipf

    r708 r794  
    742742// "backwards" wrapped to reduce redundant code 
    743743// there are only 4 choices of N (5,10,20,76) for smearing 
     744// 
     745// 
     746// 4 MAR 2011 
     747// Note: In John's paper, he integrated the Gaussian to +/- 3 sigma and then renormalized 
     748//       to an integral of 1. This "truncated" gaussian was a somewhat better approximation 
     749//       to the triangular resolution function. Here, I integrate to +/- 3 sigma and 
     750//       do not renormalize the integral to 1. Hence the smeared calculation is 0.27% low. 
     751//       This is easily seen by smearing a constant value. 
     752// 
     753// Using 5 quadrature points is not recommended, as it doesn't normalize properly using .9973 
     754//  -- instead, it normalizes to 1.0084,  
     755// 
    744756Function Smear_Model_N(fcn,w,x,resW,wi,zi,nord) 
    745757        FUNCREF SANSModelAAO_proto fcn 
    … …  
    755767// local variables 
    756768        Variable ii,va,vb 
    757         Variable answer,i_shad,i_qbar,i_sigq 
     769        Variable answer,i_shad,i_qbar,i_sigq,normalize=1 
    758770 
    759771        // current x point is the q-value for evaluation 
    … …  
    788800                //integration from va to vb 
    789801         
     802                // for +/- 3 sigma ONLY 
     803                if(nord == 5) 
     804                        normalize = 1.0057              //empirical correction, N=5 shouldn't be any different 
     805                else 
     806                        normalize = 0.9973 
     807                endif 
     808                 
    790809                va = -3*i_sigq + i_qbar 
    791810                if (va<0) 
    792811                        va=0            //to avoid numerical error when  va<0 (-ve q-value) 
    793 //                      Print "truncated Gaussian at nominal q = ",x 
     812                        Print "truncated Gaussian at nominal q = ",x 
    794813                endif 
    795814                vb = 3*i_sigq + i_qbar 
     815                 
    796816                 
    797817                // Using 20 Gauss points                     
    … …  
    815835                // all scaling, background addition... etc. is done in the model calculation 
    816836         
     837                        // renormalize to 1 
     838                        answer /= normalize 
    817839        else 
    818840                //smear with the USANS routine 
    … …  
    884906                endif 
    885907         
    886                 answer = Smear_Model_N(fcn,w,x,resW,weightW,abscissW,nord) 
     908//              answer = Smear_Model_N(fcn,w,x,resW,weightW,abscissW,nord) 
     909                Smear_Model_N_AAO(fcn,w,x,resW,weightW,abscissW,nord,answer) 
     910                 
    887911                Return (0) 
    888912        endif 
    … …  
    932956                endif 
    933957         
    934                 answer = Smear_Model_N(fcn,w,x,resW,weightW,abscissW,nord) 
     958//              answer = Smear_Model_N(fcn,w,x,resW,weightW,abscissW,nord) 
     959                Smear_Model_N_AAO(fcn,w,x,resW,weightW,abscissW,nord,answer) 
     960 
    935961                Return (0) 
    936962        endif 
    … …  
    9861012                endif 
    9871013         
    988                 answer = Smear_Model_N(fcn,w,x,resW,weightW,abscissW,nord) 
     1014//              answer = Smear_Model_N(fcn,w,x,resW,weightW,abscissW,nord) 
     1015                Smear_Model_N_AAO(fcn,w,x,resW,weightW,abscissW,nord,answer) 
     1016 
    9891017                Return (0) 
    9901018        endif 
    … …  
    10321060                endif 
    10331061         
    1034                 answer = Smear_Model_N(fcn,w,x,resW,weightW,abscissW,nord) 
     1062//              answer = Smear_Model_N(fcn,w,x,resW,weightW,abscissW,nord) 
     1063                Smear_Model_N_AAO(fcn,w,x,resW,weightW,abscissW,nord,answer) 
     1064 
    10351065                Return (0) 
    10361066        endif 
    … …  
    15651595        return(0) 
    15661596end 
    1567          
     1597 
     1598 
     1599// 7 MAR 2011 SRK 
     1600// 
     1601// calculate the resolution smearing AAO 
     1602// 
     1603// - many of the form factor calculations are threaded, so they benefit 
     1604// from being passed large numbers of q-values at once, rather than suffering the  
     1605// overhead penalty of setting up threads. 
     1606// 
     1607// In general, single integral functions benefit from this, multiple integrals not so much. 
     1608// As an example, a fit using SmearedCylinderForm took 4.3s passing nord (=20) q-values 
     1609// at a time, but only 1.1s by passing all (Nq*nord) q-values at once. For Cyl_polyRad, 
     1610// the difference was not so large, 16.2s vs. 11.9s. This is due to CylPolyRad being a  
     1611// double integral and slow enough of a calculation that passing even 20 points at once 
     1612// provides some speedup. 
     1613// 
     1614// 
     1615// 
     1616// "backwards" wrapped to reduce redundant code 
     1617// there are only 4 choices of N (5,10,20,76) for smearing 
     1618// 
     1619// 
     1620// 4 MAR 2011 SRK 
     1621// Note: In John's paper, he integrated the Gaussian to +/- 3 sigma and then renormalized 
     1622//       to an integral of 1. This "truncated" gaussian was a somewhat better approximation 
     1623//       to the triangular resolution function. Here, I integrate to +/- 3 sigma and 
     1624//       do not renormalize the integral to 1. Hence the smeared calculation is 0.27% low. 
     1625//       This is easily seen by smearing a constant value. 
     1626// 
     1627// Using 5 quadrature points is not recommended, as it doesn't normalize properly using .9973 
     1628//  -- instead, it normalizes to 1.0084. 
     1629// 
     1630Function Smear_Model_N_AAO(fcn,w,x,resW,wi,zi,nord,sm_ans) 
     1631        FUNCREF SANSModelAAO_proto fcn 
     1632        Wave w                  //coefficients of function fcn(w,x) 
     1633        WAVE x                  //x-value (q) for the calculation THIS IS PASSED IN AS A WAVE 
     1634        Wave resW               // Nx4 or NxN matrix of resolution 
     1635        Wave wi         //weight wave 
     1636        Wave zi         //abscissa wave 
     1637        Variable nord           //order of integration 
     1638        Wave sm_ans             // wave returned with the smeared model 
     1639 
     1640        NVAR dQv = root:Packages:NIST:USANS_dQv 
     1641 
     1642// local variables 
     1643        Variable ii,jj 
     1644        Variable normalize=1 
     1645        Variable nTot,num,block_sum 
     1646 
     1647 
     1648        // current x point is the q-value for evaluation 
     1649        // 
     1650        // * for the input x, the resolution function waves are interpolated to get the correct values for 
     1651        //  sigq, qbar and shad - since the model x-spacing may not be the same as 
     1652        // the experimental QSIG data. This is always the case when curve fitting, since fit_wave is  
     1653        // Igor-defined as 200 points and has its own (linear) q-(x)-scaling which will be quite different 
     1654        // from experimental data. 
     1655        // **note** if the (x) passed in is the experimental q-values, these values are 
     1656        // returned from the interpolation (as expected) 
     1657 
     1658        Make/O/D/N=(DimSize(resW, 0)) sigQ,qbar,shad,qvals,va,vb 
     1659        sigq = resW[p][0]               //std dev of resolution fn 
     1660        qbar = resW[p][1]               //mean q-value 
     1661        shad = resW[p][2]               //beamstop shadow factor 
     1662        qvals = resW[p][3]      //q-values where R(q) is known 
     1663 
     1664        //SKIP the interpolation, points passed in ARE (MUST) be the experimental q-values 
     1665         
     1666         
     1667        // if USANS data, handle separately 
     1668        // -- but this would only ever be used if the calculation was forced to use trapezoid integration 
     1669        if ( ! isSANSResolution(sigq[0]) ) 
     1670                        //smear with the USANS routine 
     1671                // Make global string and local variables 
     1672                // now data folder aware, necessary for GlobalFit = FULL path to wave    
     1673                String/G gTrap_coefStr = GetWavesDataFolder(w, 2 )       
     1674                Variable maxiter=20, tol=1e-4,uva,uvb 
     1675                 
     1676                num=numpnts(x) 
     1677                // set up limits for the integration 
     1678                uva=0 
     1679                uvb=abs(dQv) 
     1680                 
     1681                //loop over the q-values 
     1682                for(jj=0;jj<num;jj+=1) 
     1683                        Variable/G gEvalQval = x[jj] 
     1684 
     1685                        // call qtrap to do actual work 
     1686                        sm_ans[jj] = qtrap_USANS(fcn,uva,uvb,tol,maxiter) 
     1687                        sm_ans[jj] /= (uvb - uva) 
     1688                endfor 
     1689                 
     1690                return(0)        
     1691        endif 
     1692 
     1693 
     1694// now the idea is to calculate a long vector of all of the zi's (Nq * nord) 
     1695// and pass these AAO to the AAO function, to make the most use of the threading 
     1696// passing repeated short lengths of q to the function can actually be slower 
     1697// due to the overhead. 
     1698         
     1699        num = numpnts(x) 
     1700        nTot = nord*num 
     1701         
     1702        Make/O/D/N=(nTot) Resoln,yyy,xGauss,wts 
     1703 
     1704        //loop over q 
     1705        for(jj=0;jj<num;jj+=1) 
     1706         
     1707                //for each q, set up the integration range 
     1708                // end points of integration limits are technically 0-inf, but wisely choose interesting region of q where R() is nonzero 
     1709                // +/- 3 sigq catches 99.73% of distrubution 
     1710                // change limits (and spacing of zi) at each evaluation based on R() 
     1711                //integration from va to vb 
     1712 
     1713                va[jj] = -3*sigq[jj] + qbar[jj] 
     1714                if (va[jj]<0) 
     1715                        va[jj]=0                //to avoid numerical error when  va<0 (-ve q-value) 
     1716//                      Print "truncated Gaussian at nominal q = ",x 
     1717                endif 
     1718                vb[jj] = 3*sigq[jj] + qbar[jj] 
     1719         
     1720                // loop over the Gauss points 
     1721                for(ii=0;ii<nord;ii+=1) 
     1722                        // calculate Gauss points on integration interval (q-value for evaluation) 
     1723                        xGauss[nord*jj+ii] = ( zi[ii]*(vb[jj]-va[jj]) + vb[jj] + va[jj] )/2.0 
     1724                        // calculate resolution function at input q-value (use the interpolated values and zi) 
     1725                        Resoln[nord*jj+ii] = shad[jj]/sqrt(2*pi*sigq[jj]*sigq[jj]) 
     1726                        Resoln[nord*jj+ii] *= exp((-1*(xGauss[nord*jj+ii] - qbar[jj])^2)/(2*sigq[jj]*sigq[jj])) 
     1727//                      Resoln[nord*jj+ii] *= exp((-1*(xGauss[nord*jj+ii] - qvals[jj])^2)/(2*sigq[jj]*sigq[jj]))                //WRONG, but just for testing 
     1728                        // carry a copy of the weights 
     1729                        wts[nord*jj+ii] = wi[ii] 
     1730                endfor          // end of loop over quadrature points 
     1731         
     1732        endfor          //loop over q 
     1733         
     1734        //calculate AAO 
     1735        yyy = 0 
     1736        fcn(w,yyy,xGauss)               //yyy is the return value as a wave 
     1737 
     1738        //multiply by weights 
     1739        yyy *= wts*Resoln               //multiply function by resolution and weights 
     1740         
     1741        //sum up blockwise to get the final answer 
     1742        for(jj=0;jj<num;jj+=1) 
     1743                block_sum = 0 
     1744                for(ii=0;ii<nord;ii+=1) 
     1745                        block_sum += yyy[nord*jj+ii] 
     1746                endfor 
     1747                sm_ans[jj] = (vb[jj]-va[jj])/2.0*block_sum 
     1748        endfor 
     1749         
     1750         
     1751        // then normalize for +/- 3 sigma ONLY 
     1752        if(nord == 5) 
     1753                normalize = 1.0057              //empirical correction, N=5 shouldn't be any different 
     1754        else 
     1755                normalize = 0.9973 
     1756        endif 
     1757         
     1758        sm_ans /= normalize 
     1759         
     1760        return(0) 
     1761         
     1762End 
     1763 

  • sans/Dev/trunk/NCNR_User_Procedures/Common/Packages/PlotManager/PlotUtils2D_v40.ipf

    r771 r794  
    1717//                                              2D model calculations (model_lin) 2D data files as well. 
    1818 
     19// Call the procedure BinQxQy_to_1D() to re-bin the QxQy data into I(q). good for testing to check the 2D error propagation 
     20// 
     21// 
     22 
    1923 
    2024// 
    … …  
    2529Proc LoadQxQy() 
    2630 
     31        SetDataFolder root: 
     32         
    2733        LoadWave/G/D/W/A 
    2834        String fileName = S_fileName 
    … …  
    10181024// be made too, if needed. 
    10191025// 
    1020 // X- need error on I(q) 
    1021 // -- need to set "proper" number of data points (delta and qMax?) 
     1026// X- need error on I(q) - now calculated both ways, from 2D error propagation, and from deviation from the mean 
     1027// X- need to set "proper" number of data points (delta and qMax?) 
    10221028// X- need to remove points at high Q end 
    10231029// 
    … …  
    10271033// the results are in iBin_qxqy, qBin_qxqy, and eBin_qxqy, in the folder passed 
    10281034//  
     1035 
     1036Proc BinQxQy_to_1D(folderStr) 
     1037        String folderStr 
     1038        Prompt folderStr,"Pick the data folder containing 2D data",popup,getAList(4) 
     1039 
     1040        fDoBinning_QxQy2D(folderStr) 
     1041End 
     1042 
     1043 
     1044 
    10291045//Function fDoBinning_QxQy2D(inten,qx,qy,qz) 
     1046//      Wave inten,qx,qy,qz 
     1047 
    10301048Function fDoBinning_QxQy2D(folderStr) 
    10311049        String folderStr 
    10321050 
    1033 //      Wave inten,qx,qy,qz 
    10341051 
    10351052        SetDataFolder $("root:"+folderStr) 
    10361053         
    1037         WAVE inten = $(folderStr + "_i") 
     1054        WAVE inten = $("smeared_sf2D") 
     1055 
     1056//      WAVE inten = $(folderStr + "_i") 
     1057        WAVE iErr = $(folderStr + "_iErr") 
    10381058        WAVE qx = $(folderStr + "_qx") 
    10391059        WAVE qy = $(folderStr + "_qy") 
    … …  
    10481068         
    10491069        yDim = XDim 
    1050         Make/O/D/N=(nq) iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy 
     1070        Make/O/D/N=(nq) iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy 
    10511071        delQ = abs(sqrt(qx[2]^2+qy[2]^2+qz[2]^2) - sqrt(qx[1]^2+qy[1]^2+qz[1]^2))               //use bins of 1 pixel width  
    10521072        qBin_qxqy[] =  p*       delQ     
    … …  
    10561076        iBin2_qxqy = 0 
    10571077        eBin_qxqy = 0 
     1078        eBin2D_qxqy = 0 
    10581079        nBin_qxqy = 0   //number of intensities added to each bin 
    10591080         
    … …  
    10651086                        iBin_qxqy[binIndex] += val 
    10661087                        iBin2_qxqy[binIndex] += val*val 
     1088                        eBin2D_qxqy[binIndex] += iErr[ii]*iErr[ii] 
    10671089                        nBin_qxqy[binIndex] += 1 
    10681090                endif 
    … …  
    10751097                        iBin_qxqy[ii] = 0 
    10761098                        eBin_qxqy[ii] = 1 
     1099                        eBin2D_qxqy[ii] = NaN 
    10771100                else 
    10781101                        if(nBin_qxqy[ii] <= 1) 
    … …  
    10801103                                iBin_qxqy[ii] /= nBin_qxqy[ii] 
    10811104                                eBin_qxqy[ii] = 1 
     1105                                eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2 
    10821106                        else 
    10831107                                //assume that the intensity in each pixel in annuli is normally 
    … …  
    10921116                                        eBin_qxqy[ii] = sqrt(var/(nBin_qxqy[ii] - 1)) 
    10931117                                endif 
     1118                                // and calculate as it is propagated pixel-by-pixel 
     1119                                eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2 
    10941120                        endif 
    10951121                endif 
    10961122        endfor 
     1123         
     1124        eBin2D_qxqy = sqrt(eBin2D_qxqy)         // as equation (3) of John's memo 
    10971125         
    10981126        // find the last non-zero point, working backwards 
    … …  
    11031131         
    11041132//      print val, nBin_qxqy[val] 
    1105         DeletePoints val, nq-val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy 
     1133        DeletePoints val, nq-val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy 
    11061134         
    11071135        SetDataFolder root: 

  • sans/Dev/trunk/NCNR_User_Procedures/Common/Smear_2D.ipf

    r708 r794  
    9595         
    9696        Variable spl,spp,apl,app,bpl,bpp,phi_pt,qpl_pt 
     97        Variable qperp_pt,phi_prime,q_prime 
     98 
    9799        Variable t1=StopMSTimer(-2) 
    98100 
    … …  
    118120                apl = -numStdDev*spl + qval             //parallel = q integration limits 
    119121                bpl = numStdDev*spl + qval 
    120                 app = -numStdDev*spp + phi              //perpendicular = phi integration limits 
    121                 bpp = numStdDev*spp + phi 
     122///             app = -numStdDev*spp + phi              //perpendicular = phi integration limits (WRONG) 
     123///             bpp = numStdDev*spp + phi 
     124                app = -numStdDev*spp + 0                //q_perp = 0 
     125                bpp = numStdDev*spp + 0 
    122126                 
    123127                //make sure the limits are reasonable. 
    … …  
    130134                sumOut = 0 
    131135                for(jj=0;jj<nord;jj+=1)         // call phi the "outer' 
    132                         phi_pt = (xi[jj]*(bpp-app)+app+bpp)/2 
     136///                     phi_pt = (xi[jj]*(bpp-app)+app+bpp)/2 
     137                        qperp_pt = (xi[jj]*(bpp-app)+app+bpp)/2         //this is now q_perp 
    133138 
    134139                        sumIn=0 
    … …  
    137142                                qpl_pt = (xi[kk]*(bpl-apl)+apl+bpl)/2 
    138143                                 
    139                                 FindQxQy(qpl_pt,phi_pt,qx_pt,qy_pt)             //find the corresponding QxQy to the Q,phi 
     144///                             FindQxQy(qpl_pt,phi_pt,qx_pt,qy_pt)             //find the corresponding QxQy to the Q,phi 
     145 
     146                                // find QxQy given Qpl and Qperp on the grid 
     147                                // 
     148                                q_prime = sqrt(qpl_pt^2+qperp_pt^2) 
     149                                phi_prime = phi + qperp_pt/qpl_pt 
     150                                FindQxQy(q_prime,phi_prime,qx_pt,qy_pt) 
     151                                 
    140152                                yPtw[kk] = qy_pt                                        //phi is the same in this loop, but qy is not 
    141153                                xPtW[kk] = qx_pt                                        //qx is different here too, as we're varying Qpl 
    142                                  
    143                                 res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (phi_pt-phi)^2/spp/spp ) ) 
     154 
     155                                res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (qperp_pt)^2/spp/spp ) ) 
     156///                             res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (phi_pt-phi)^2/spp/spp ) ) 
    144157                                res_tot[kk] /= normFactor 
    145158//                              res_tot[kk] *= fs 
    … …  
    164177        return(0) 
    165178end 
     179 
     180 
     181// this is generic, but I need to declare the Cylinder2D function threadsafe 
     182// and this calculation is significantly slower than the manually threaded calculation 
     183// if the function is fast to calculate. Once the function has polydispersity on 2 or more parameters 
     184// then this AAO calculation and the manual threading are both painfully slow, and more similar in execution time 
     185// 
     186// For 128x128 data, and 10x10 smearing, using 25 points for each polydispersity (using DANSE code) 
     187//      Successively making more things polydisperse gives the following timing (in seconds): 
     188// 
     189//                              monodisp                sigR            sigTheta                sigPhi 
     190// manual THR           1.1             3.9                     74                      1844 
     191//      this AAO                        8.6             11.4                    104             1930 
     192// 
     193// and using 5 points for each polydispersity: (I see no visual difference) 
     194// manual THR           1.1             1.6                     3.9             16.1 
     195// 
     196// so clearly -- use 5 points for the polydispersities, unless there's a good reason not to - and  
     197// certainly for a survey, it's the way to go. 
     198// 
     199Function Smear_2DModel_PP_AAO(fcn,s,nord) 
     200        FUNCREF SANS_2D_ModelAAO_proto fcn 
     201        Struct ResSmear_2D_AAOStruct &s 
     202        Variable nord 
     203         
     204        String weightStr,zStr 
     205 
     206        Variable ii,jj,kk,num 
     207        Variable qx,qy,qz,qval,fs 
     208        Variable qy_pt,qx_pt,res_x,res_y,answer,sumIn,sumOut 
     209         
     210        Variable a,b,c,normFactor,phi,theta,maxSig,numStdDev=3 
     211         
     212        switch(nord)     
     213                case 5:          
     214                        weightStr="gauss5wt" 
     215                        zStr="gauss5z" 
     216                        if (WaveExists($weightStr) == 0) 
     217                                Make/O/D/N=(nord) $weightStr,$zStr 
     218                                Make5GaussPoints($weightStr,$zStr)       
     219                        endif 
     220                        break                            
     221                case 10:                 
     222                        weightStr="gauss10wt" 
     223                        zStr="gauss10z" 
     224                        if (WaveExists($weightStr) == 0) 
     225                                Make/O/D/N=(nord) $weightStr,$zStr 
     226                                Make10GaussPoints($weightStr,$zStr)      
     227                        endif 
     228                        break                            
     229                case 20:                 
     230                        weightStr="gauss20wt" 
     231                        zStr="gauss20z" 
     232                        if (WaveExists($weightStr) == 0) 
     233                                Make/O/D/N=(nord) $weightStr,$zStr 
     234                                Make20GaussPoints($weightStr,$zStr)      
     235                        endif 
     236                        break 
     237                default:                                                         
     238                        Abort "Smear_2DModel_PP called with invalid nord value"                                  
     239        endswitch 
     240         
     241        Wave/Z wt = $weightStr 
     242        Wave/Z xi = $zStr 
     243         
     244/// keep these waves local 
     245//      Make/O/D/N=1 yPtW 
     246        Make/O/D/N=(nord*nord) fcnRet,xptW,res_tot,yptW,wts 
     247        Make/O/D/N=(nord) phi_pt,qpl_pt,qperp_pt 
     248                 
     249        // now just loop over the points as specified 
     250        num=numpnts(s.xw[0]) 
     251         
     252        answer=0 
     253         
     254        Variable spl,spp,apl,app,bpl,bpp 
     255        Variable phi_prime,q_prime 
     256 
     257        Variable t1=StopMSTimer(-2) 
     258 
     259        //loop over q-values 
     260        for(ii=0;ii<num;ii+=1) 
     261         
     262//              if(mod(ii, 1000 ) == 0) 
     263//                      Print "ii= ",ii 
     264//              endif 
     265                 
     266                qx = s.xw[0][ii] 
     267                qy = s.xw[1][ii] 
     268                qz = s.qz[ii] 
     269                qval = sqrt(qx^2+qy^2+qz^2) 
     270                spl = s.sQpl[ii] 
     271                spp = s.sQpp[ii] 
     272                fs = s.fs[ii] 
     273                 
     274                normFactor = 2*pi*spl*spp 
     275                 
     276                phi = -1*FindTheta(qx,qy)               //Findtheta is an exact duplicate of FindPhi() * -1 
     277                 
     278                apl = -numStdDev*spl + qval             //parallel = q integration limits 
     279                bpl = numStdDev*spl + qval 
     280///             app = -numStdDev*spp + phi              //perpendicular = phi integration limits (WRONG) 
     281///             bpp = numStdDev*spp + phi 
     282                app = -numStdDev*spp + 0                //q_perp = 0 
     283                bpp = numStdDev*spp + 0 
     284                 
     285                //make sure the limits are reasonable. 
     286                if(apl < 0) 
     287                        apl = 0 
     288                endif 
     289                // do I need to specially handle limits when phi ~ 0? 
     290         
     291                 
     292//              sumOut = 0 
     293                for(jj=0;jj<nord;jj+=1)         // call phi the "outer' 
     294//                      phi_pt[jj] = (xi[jj]*(bpp-app)+app+bpp)/2 
     295                        qperp_pt[jj] = (xi[jj]*(bpp-app)+app+bpp)/2             //this is now q_perp 
     296                         
     297//                      sumIn=0 
     298                        for(kk=0;kk<nord;kk+=1)         //at phi, integrate over Qpl 
     299 
     300                                qpl_pt[kk] = (xi[kk]*(bpl-apl)+apl+bpl)/2 
     301                                 
     302///                             FindQxQy(qpl_pt[kk],phi_pt[jj],qx_pt,qy_pt)             //find the corresponding QxQy to the Q,phi 
     303                                 
     304                                // find QxQy given Qpl and Qperp on the grid 
     305                                // 
     306                                q_prime = sqrt(qpl_pt[kk]^2+qperp_pt[jj]^2) 
     307                                phi_prime = phi + qperp_pt[jj]/qpl_pt[kk] 
     308                                FindQxQy(q_prime,phi_prime,qx_pt,qy_pt) 
     309                                                                 
     310                                yPtw[nord*jj+kk] = qy_pt                                        //phi is the same in this loop, but qy is not 
     311                                xPtW[nord*jj+kk] = qx_pt                                        //qx is different here too, as we're varying Qpl 
     312                                 
     313                                res_tot[nord*jj+kk] = exp(-0.5*( (qpl_pt[kk]-qval)^2/spl/spl + (qperp_pt[jj])^2/spp/spp ) ) 
     314//                              res_tot[nord*jj+kk] = exp(-0.5*( (qpl_pt[kk]-qval)^2/spl/spl + (phi_pt[jj]-phi)^2/spp/spp ) ) 
     315                                res_tot[nord*jj+kk] /= normFactor 
     316//                              res_tot[kk] *= fs 
     317 
     318                                //weighting 
     319                                wts[nord*jj+kk] = wt[jj]*wt[kk] 
     320                        endfor 
     321                         
     322                endfor 
     323                 
     324                fcn(s.coefW,fcnRet,xptw,yptw)                   //calculate nord*nord pts at a time 
     325                 
     326                fcnRet *= wts*res_tot 
     327                // 
     328                answer = (bpl-apl)/2.0*sum(fcnRet)              // get the sum, normalize to parallel direction 
     329                answer *= (bpp-app)/2.0                                         // and normalize to perpendicular direction 
     330                 
     331                s.zw[ii] = answer 
     332        endfor 
     333         
     334        Variable elap = (StopMSTimer(-2) - t1)/1e6 
     335        Print "elapsed time = ",elap 
     336         
     337        return(0) 
     338end 
     339 
    166340 
    167341 

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