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/Analysis/Models/Models_2D
Files:
6 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 
Note: See TracChangeset for help on using the changeset viewer.