Ignore:
Timestamp:
May 25, 2017 1:29:50 PM (6 years ago)
Author:
srkline
Message:

Added the angle dependent transmission correction to the data correction in the raw_to_work step, in 2D

added a testing file that can generate fake event data, read, write, and decode it. Read is based on GBLoadWave. Hoepfully I'll not need to write an XOP. manipulation of the 64 bit words are done with simple bit shifts and logic.

also added are a number of error checking routines to improve behavior when wave, folders, etc. are missing.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_DetectorCorrections.ipf

    r1027 r1041  
    676676        // TODO: 
    677677        // correctly apply the correction to the error wave (assume a perfect value?) 
    678         // w_err /= solid_angle         //is this correct?? 
     678        w_err /= solid_angle            //is this correct?? 
    679679 
    680680// TODO -- clean up after I'm satisfied computations are correct                 
     
    793793                                 
    794794                                // pass in the transmission error, and the error in the correction is returned as the last parameter 
    795                                 lat_corr = V_LargeAngleTransmissionCorr(trans,dtdist,xd,yd,trans_err,lat_err)           //moved from 1D avg SRK 11/2007 
     795 
     796//                              lat_corr = V_LargeAngleTransmissionCorr(trans,dtdist,xd,yd,trans_err,lat_err)           //moved from 1D avg SRK 11/2007 
     797 
    796798                                data[ii][jj] /= lat_corr                        //divide by the correction factor 
    797799                                // 
     
    822824 
    823825 
    824  
     826// 
    825827// DIVIDE the intensity by this correction to get the right answer 
    826828// TODO: 
     
    828830// 
    829831// 
    830 Function V_LargeAngleTransmissionCorr(trans,dtdist,xd,yd,trans_err,err) 
    831         Variable trans,dtdist,xd,yd,trans_err,&err 
    832  
    833         DoAlert 0,"This has not yet been updated for VSANS" 
    834          
    835         //angle dependent transmission correction  
    836         Variable uval,arg,cos_th,correction,theta 
    837          
    838         ////this section is the trans_correct() VAX routine 
    839 //      if(trans<0.1) 
    840 //              Print "***transmission is less than 0.1*** and is a significant correction" 
    841 //      endif 
    842 //      if(trans==0) 
    843 //              Print "***transmission is ZERO*** and has been reset to 1.0 for the averaging calculation" 
    844 //              trans = 1 
    845 //      endif 
    846          
    847         theta = atan( (sqrt(xd^2 + yd^2))/dtdist )              //theta at the input pixel 
     832 
     833// Apply the large angle transmssion correction as the data is converted to WORK 
     834// so that whether the data is saved as 2D or 1D, the correction has properly been done. 
     835// 
     836// This is, however, a SAMPLE dependent calculation, not purely instrument geometry. 
     837// 
     838Function V_LargeAngleTransmissionCorr(w,w_err,fname,detStr,destPath) 
     839        Wave w,w_err 
     840        String fname,detStr,destPath 
     841 
     842        Variable sdd,xCtr,yCtr,trans,trans_err,uval 
     843 
     844// get all of the geometry information   
     845//      orientation = V_getDet_tubeOrientation(fname,detStr) 
     846        sdd = V_getDet_ActualDistance(fname,detStr) 
     847        sdd/=100                // sdd in cm, pass in m 
     848 
     849        // this is ctr in mm 
     850        xCtr = V_getDet_beam_center_x_mm(fname,detStr) 
     851        yCtr = V_getDet_beam_center_y_mm(fname,detStr) 
     852        trans = V_getSampleTransmission(fname) 
     853        trans_err = V_getSampleTransError(fname) 
     854         
     855        SetDataFolder $(destPath + ":entry:instrument:detector_"+detStr) 
     856         
     857        Wave data_realDistX = data_realDistX 
     858        Wave data_realDistY = data_realDistY 
     859 
     860        Duplicate/O w lat_corr,tmp_theta,tmp_dist,lat_err,tmp_err               //in the current df 
     861 
     862//// calculate the scattering angle 
     863//      dx = (distX - xctr)             //delta x in mm 
     864//      dy = (distY - yctr)             //delta y in mm 
     865        tmp_dist = sqrt((data_realDistX - xctr)^2 + (data_realDistY - yctr)^2) 
     866         
     867        tmp_dist /= 10  // convert mm to cm 
     868        sdd *=100               //convert to cm 
     869 
     870        tmp_theta = atan(tmp_dist/sdd)          //this is two_theta, the scattering angle 
     871 
     872        Variable ii,jj,numx,numy,dx,dy,cos_th,arg,tmp 
     873        numx = DimSize(tmp_theta,0) 
     874        numy = DimSize(tmp_theta,1) 
     875         
    848876         
    849877        //optical thickness 
    850878        uval = -ln(trans)               //use natural logarithm 
    851         cos_th = cos(theta) 
    852         arg = (1-cos_th)/cos_th 
    853          
    854         // a Taylor series around uval*arg=0 only needs about 4 terms for very good accuracy 
    855         //                      correction= 1 - 0.5*uval*arg + (uval*arg)^2/6 - (uval*arg)^3/24 + (uval*arg)^4/120 
    856         // OR 
    857         if((uval<0.01) || (cos_th>0.99))         
    858                 //small arg, approx correction 
    859                 correction= 1-0.5*uval*arg 
    860         else 
    861                 //large arg, exact correction 
    862                 correction = (1-exp(-uval*arg))/(uval*arg) 
    863         endif 
    864  
    865         Variable tmp 
    866          
    867         if(trans == 1) 
    868                 err = 0         //no correction, no error 
    869         else 
    870                 //sigT, calculated from the Taylor expansion 
    871                 tmp = (1/trans)*(arg/2-arg^2/3*uval+arg^3/8*uval^2-arg^4/30*uval^3) 
    872                 tmp *= tmp 
    873                 tmp *= trans_err^2 
    874                 tmp = sqrt(tmp)         //sigT 
    875                  
    876                 err = tmp 
    877         endif 
    878          
    879 //      Printf "trans error = %g\r",trans_err 
    880 //      Printf "correction = %g +/- %g\r", correction, err 
    881          
    882         //end of transmission/pathlength correction 
    883  
    884         return(correction) 
     879         
     880        for(ii=0        ;ii<numx;ii+=1) 
     881                for(jj=0;jj<numy;jj+=1) 
     882                         
     883                        cos_th = cos(tmp_theta[ii][jj]) 
     884                        arg = (1-cos_th)/cos_th 
     885                         
     886                        // a Taylor series around uval*arg=0 only needs about 4 terms for very good accuracy 
     887                        //                      correction= 1 - 0.5*uval*arg + (uval*arg)^2/6 - (uval*arg)^3/24 + (uval*arg)^4/120 
     888                        // OR 
     889                        if((uval<0.01) || (cos_th>0.99))         
     890                                //small arg, approx correction 
     891                                lat_corr[ii][jj] = 1-0.5*uval*arg 
     892                        else 
     893                                //large arg, exact correction 
     894                                lat_corr[ii][jj] = (1-exp(-uval*arg))/(uval*arg) 
     895                        endif 
     896                          
     897                        // TODO 
     898                        // -- properly calculate and apply the 2D error propagation 
     899                        if(trans == 1) 
     900                                lat_err[ii][jj] = 0             //no correction, no error 
     901                        else 
     902                                //sigT, calculated from the Taylor expansion 
     903                                tmp = (1/trans)*(arg/2-arg^2/3*uval+arg^3/8*uval^2-arg^4/30*uval^3) 
     904                                tmp *= tmp 
     905                                tmp *= trans_err^2 
     906                                tmp = sqrt(tmp)         //sigT 
     907                                 
     908                                lat_err[ii][jj] = tmp 
     909                        endif 
     910                          
     911  
     912                endfor 
     913        endfor 
     914         
     915 
     916         
     917        // Here it is! Apply the correction to the intensity (divide -- to get the proper correction) 
     918        w /= lat_corr 
     919 
     920        // relative errors add in quadrature to the current 2D error 
     921        tmp_err = (w_err/lat_corr)^2 + (lat_err/lat_corr)^2*w*w/lat_corr^2 
     922        tmp_err = sqrt(tmp_err) 
     923         
     924        w_err = tmp_err  
     925         
     926        // TODO: 
     927        // correctly apply the correction to the error wave (assume a perfect value?) 
     928        // w_err /= tmp         //is this correct?? 
     929 
     930        // TODO -- clean up after I'm satisfied computations are correct                 
     931        KillWaves/Z tmp_theta,tmp_dist,tmp_err,lat_err 
     932         
     933        return(0) 
    885934end 
    886935 
Note: See TracChangeset for help on using the changeset viewer.