Ignore:
Timestamp:
Feb 14, 2020 2:44:16 PM (3 years ago)
Author:
srkline
Message:

adding downstream window transmission correction -- added function to detector corrections, and added item to preference panel. Transmission value is currently set as a global since the value is not part of the VSANS header. Global value defaults to T=1= no correction.

Other changes are cleanup of TODO items that were already done and tested. inline commnents have been updated.

File:
1 edited

Legend:

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

    r1236 r1239  
    22002200///////////// 
    22012201 
     2202// 
     2203// 
     2204// testing function to calculate the correction for the attenuation 
     2205// of the scattered beam by windows downstream of the sample 
     2206// (the back window of the sample block, the Si window) 
     2207// 
     2208// For implementation, this function could be made identical 
     2209// to the large angle transmission correction, since the math is 
     2210// identical - only the Tw value is different (and should ideally be 
     2211// quite close to 1). With Tw near 1, this would be a few percent correction 
     2212// at the largest scattering angles. 
     2213// 
     2214// 
     2215Function V_WindowTransmission(tw) 
     2216        Variable tw 
     2217         
     2218        Make/O/D/N=100 theta,method1,method2,arg 
     2219         
     2220        theta = p/2 
     2221        theta = theta/360*2*pi          //convert to radians 
     2222         
     2223//      method1 = exp( -ln(tw)/cos(theta) )/tw 
     2224         
     2225        Variable tau 
     2226        tau = -ln(tw) 
     2227        arg = (1-cos(theta))/cos(theta) 
     2228         
     2229        if(tau < 0.01) 
     2230                method2 = 1 - 0.5*tau*arg        
     2231        else 
     2232                method2 = ( 1 - exp(-tau*arg) )/(tau*arg) 
     2233        endif 
     2234         
     2235        return(0) 
     2236end 
     2237 
     2238 
     2239// 
     2240// Large angle transmission correction for the downstream window 
     2241// 
     2242// DIVIDE the intensity by this correction to get the right answer 
     2243// 
     2244// -- this is a duplication of the math for the large angle 
     2245// sample tranmission correction. Same situation, but now the  
     2246// scattered neutrons are attenuated by whatever windows are  
     2247// downstream of the scattering event. 
     2248// (the back window of the sample block, the Si window) 
     2249// 
     2250// For implementation, this function is made identical 
     2251// to the large angle transmission correction, since the math is 
     2252// identical - only the Tw value is different (and should ideally be 
     2253// quite close to 1). With Tw near 1, this would be a few percent correction 
     2254// at the largest scattering angles. 
     2255// 
     2256// 
     2257Function V_DownstreamWindowTransmission(w,w_err,fname,detStr,destPath) 
     2258        Wave w,w_err 
     2259        String fname,detStr,destPath 
     2260 
     2261        Variable sdd,xCtr,yCtr,uval 
     2262 
     2263// get all of the geometry information   
     2264//      orientation = V_getDet_tubeOrientation(fname,detStr) 
     2265        sdd = V_getDet_ActualDistance(fname,detStr) 
     2266 
     2267        // this is ctr in mm 
     2268        xCtr = V_getDet_beam_center_x_mm(fname,detStr) 
     2269        yCtr = V_getDet_beam_center_y_mm(fname,detStr) 
     2270 
     2271// get the value of the overall transmission of the downstream components 
     2272// + error if available. 
     2273//      trans = V_getSampleTransmission(fname) 
     2274//      trans_err = V_getSampleTransError(fname) 
     2275// TODO -- HARD WIRED values, need to set a global or find a place in the header (instrument block?) (reduction?) 
     2276// currently globals are forced to one in WorkFolderUtils.ipf as the correction is done 
     2277        NVAR trans = root:Packages:NIST:VSANS:Globals:gDownstreamWinTrans 
     2278        NVAR trans_err = root:Packages:NIST:VSANS:Globals:gDownstreamWinTransErr         
     2279 
     2280        SetDataFolder $(destPath + ":entry:instrument:detector_"+detStr) 
     2281         
     2282        Wave data_realDistX = data_realDistX 
     2283        Wave data_realDistY = data_realDistY 
     2284 
     2285        Duplicate/O w dwt_corr,tmp_theta,tmp_dist,dwt_err,tmp_err               //in the current df 
     2286 
     2287//// calculate the scattering angle 
     2288//      dx = (distX - xctr)             //delta x in mm 
     2289//      dy = (distY - yctr)             //delta y in mm 
     2290        tmp_dist = sqrt((data_realDistX - xctr)^2 + (data_realDistY - yctr)^2) 
     2291         
     2292        tmp_dist /= 10  // convert mm to cm 
     2293        // sdd is in [cm] 
     2294 
     2295        tmp_theta = atan(tmp_dist/sdd)          //this is two_theta, the scattering angle 
     2296 
     2297        Variable ii,jj,numx,numy,dx,dy,cos_th,arg,tmp 
     2298        numx = DimSize(tmp_theta,0) 
     2299        numy = DimSize(tmp_theta,1) 
     2300         
     2301         
     2302        //optical thickness 
     2303        uval = -ln(trans)               //use natural logarithm 
     2304         
     2305        for(ii=0        ;ii<numx;ii+=1) 
     2306                for(jj=0;jj<numy;jj+=1) 
     2307                         
     2308                        cos_th = cos(tmp_theta[ii][jj]) 
     2309                        arg = (1-cos_th)/cos_th 
     2310                         
     2311                        // a Taylor series around uval*arg=0 only needs about 4 terms for very good accuracy 
     2312                        //                      correction= 1 - 0.5*uval*arg + (uval*arg)^2/6 - (uval*arg)^3/24 + (uval*arg)^4/120 
     2313                        // OR 
     2314                        if((uval<0.01) || (cos_th>0.99))         
     2315                                //small arg, approx correction 
     2316                                dwt_corr[ii][jj] = 1-0.5*uval*arg 
     2317                        else 
     2318                                //large arg, exact correction 
     2319                                dwt_corr[ii][jj] = (1-exp(-uval*arg))/(uval*arg) 
     2320                        endif 
     2321                          
     2322                        // (DONE) 
     2323                        // x- properly calculate and apply the 2D error propagation 
     2324                        if(trans == 1) 
     2325                                dwt_err[ii][jj] = 0             //no correction, no error 
     2326                        else 
     2327                                //sigT, calculated from the Taylor expansion 
     2328                                tmp = (1/trans)*(arg/2-arg^2/3*uval+arg^3/8*uval^2-arg^4/30*uval^3) 
     2329                                tmp *= tmp 
     2330                                tmp *= trans_err^2 
     2331                                tmp = sqrt(tmp)         //sigT 
     2332                                 
     2333                                dwt_err[ii][jj] = tmp 
     2334                        endif 
     2335                          
     2336                endfor 
     2337        endfor 
     2338         
     2339        // Here it is! Apply the correction to the intensity (divide -- to get the proper correction) 
     2340        w /= dwt_corr 
     2341 
     2342        // relative errors add in quadrature to the current 2D error 
     2343        tmp_err = (w_err/dwt_corr)^2 + (dwt_err/dwt_corr)^2*w*w/dwt_corr^2 
     2344        tmp_err = sqrt(tmp_err) 
     2345         
     2346        w_err = tmp_err  
     2347         
     2348        // DONE x- clean up after I'm satisfied computations are correct                 
     2349        KillWaves/Z tmp_theta,tmp_dist,tmp_err,dwt_err 
     2350         
     2351        return(0) 
     2352end 
     2353 
Note: See TracChangeset for help on using the changeset viewer.