Ignore:
Timestamp:
May 20, 2016 4:28:56 PM (7 years ago)
Author:
srkline
Message:

changes to a few analysis models to make these Igor 7-ready

adding mask editing utilities

many changes to event mode for easier processing of split lists

updated event mode help file

+ lots more!

File:
1 edited

Legend:

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

    r995 r999  
    370370                if(cmpstr(orientation,"vertical")==0) 
    371371                //      this is vertical tube data dimensioned as (Ntubes,Npix) 
    372                         pixSize = 8                     //V_getDet_y_pixel_size(fname,detStr) 
     372                        pixSize = 8.4           //V_getDet_y_pixel_size(fname,detStr) 
    373373                         
    374374                elseif(cmpstr(orientation,"horizontal")==0) 
     
    396396 
    397397// 
     398// TODO: 
     399// -- MUST VERIFY the definition of SDD and how (if) setback is written to the data files 
     400// -- currently I'm assuming that the SDD is the "nominal" value which is correct for the  
     401//    L/R panels, but is not correct for the T/B panels (must add in the setback) 
     402// 
     403// 
     404// 
    398405// data_realDistX, Y must be previously generated from running NonLinearCorrection() 
    399406// 
     
    411418// get all of the geometry information   
    412419        orientation = V_getDet_tubeOrientation(fname,detStr) 
    413         sdd = V_getDet_distance(fname,detStr) 
     420//      sdd = V_getDet_distance(fname,detStr)           //[cm] 
     421//      sdd += V_getDet_TBSetback(fname,detStr)/10              // written in [mm], convert to [cm], returns 0 for L/R/B panels 
     422 
     423        sdd = V_getDet_ActualDistance(fname,detStr)             //sdd derived, including setback [cm] 
    414424        sdd/=100                // sdd reported in cm, pass in m 
    415425        // this is the ctr in pixels 
     
    552562//    to each pixel 
    553563// 
    554 // not actually used for anything, but here for completeness if anyone asks 
     564// not actually used for any calculations, but here for completeness if anyone asks, or for 2D data export 
    555565// 
    556566// this properly accounts for qz 
     
    579589End 
    580590 
     591 
     592// 
     593// TODO -- VERIFY calculations 
     594// -- This is the actual solid angle per pixel, not a ratio vs. some "unit SA"  
     595//    Do I just correct for the different area vs. the "nominal" central area? 
     596// -- decide how to implement - either directly change the data values (as was done in the past) 
     597//    or use this as a weighting for when the data is binned to I(q). In the second method, 2D data 
     598//    would need this to be applied before exporting 
     599// -- do I keep a wave note indicating that this correction has been applied to the data 
     600//    so that it can be "un-applied"? 
     601// -- do I calculate theta from geometry directly, or get it from Q (Assuming it's present?) 
     602//    (probably just from geometry, since I need SDD and dx and dy values...) 
     603// 
     604// 
     605Function SolidAngleCorrection(w,w_err,fname,detStr,destPath) 
     606        Wave w,w_err 
     607        String fname,detStr,destPath 
     608 
     609        Variable sdd,xCtr,yCtr,lambda 
     610 
     611// get all of the geometry information   
     612//      orientation = V_getDet_tubeOrientation(fname,detStr) 
     613        sdd = V_getDet_ActualDistance(fname,detStr) 
     614        sdd/=100                // sdd in cm, pass in m 
     615 
     616        // this is ctr in mm 
     617        xCtr = V_getDet_beam_center_x_mm(fname,detStr) 
     618        yCtr = V_getDet_beam_center_y_mm(fname,detStr) 
     619        lambda = V_getWavelength(fname) 
     620         
     621        SetDataFolder $(destPath + ":entry:instrument:detector_"+detStr) 
     622         
     623        Wave data_realDistX = data_realDistX 
     624        Wave data_realDistY = data_realDistY 
     625 
     626        Duplicate/O w solid_angle,tmp_theta,tmp_dist            //in the current df 
     627 
     628//// calculate the scattering angle 
     629//      dx = (distX - xctr)             //delta x in mm 
     630//      dy = (distY - yctr)             //delta y in mm 
     631        tmp_dist = sqrt((data_realDistX - xctr)^2 + (data_realDistY - yctr)^2) 
     632         
     633        tmp_dist /= 10  // convert mm to cm 
     634        sdd *=100               //convert to cm 
     635 
     636        tmp_theta = atan(tmp_dist/sdd)          //this is two_theta, the scattering angle 
     637 
     638        Variable ii,jj,numx,numy,dx,dy 
     639        numx = DimSize(tmp_theta,0) 
     640        numy = DimSize(tmp_theta,1) 
     641         
     642        for(ii=0        ;ii<numx;ii+=1) 
     643                for(jj=0;jj<numy;jj+=1) 
     644                         
     645                        if(ii==0)               //do a forward difference if ii==0 
     646                                dx = (data_realDistX[ii+1][jj] - data_realDistX[ii][jj])        //delta x for the pixel 
     647                        else 
     648                                dx = (data_realDistX[ii][jj] - data_realDistX[ii-1][jj])        //delta x for the pixel 
     649                        endif 
     650                         
     651                         
     652                        if(jj==0) 
     653                                dy = (data_realDistY[ii][jj+1] - data_realDistY[ii][jj])        //delta y for the pixel 
     654                        else 
     655                                dy = (data_realDistY[ii][jj] - data_realDistY[ii][jj-1])        //delta y for the pixel 
     656                        endif 
     657         
     658                        dx /= 10 
     659                        dy /= 10                // convert mm to cm (since sdd is in cm) 
     660                        solid_angle[ii][jj] = dx*dy             //this is in cm^2 
     661                endfor 
     662        endfor 
     663         
     664        // to cover up any issues w/negative dx or dy 
     665        solid_angle = abs(solid_angle) 
     666         
     667        // solid_angle correction 
     668        // == dx*dy*cos^3/sdd^2 
     669        solid_angle *= (cos(tmp_theta))^3 
     670        solid_angle /= sdd^2 
     671         
     672        // Here it is! Apply the correction to the intensity (I divide -- to get the counts per solid angle!!) 
     673        w /= solid_angle 
     674         
     675         
     676        // TODO: 
     677        // correctly apply the correction to the error wave (assume a perfect value?) 
     678        // w_err /= solid_angle         //is this correct?? 
     679 
     680// TODO -- clean up after I'm satisfied computations are correct                 
     681//      KillWaves/Z tmp_theta,tmp_dist 
     682         
     683        return(0) 
     684end 
    581685 
    582686 
     
    668772                        // so divide here to get the correct answer (5/22/08 SRK) 
    669773                        if(doEfficiency) 
    670                                 data[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd) 
    671                                 data_err[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd) 
     774//                              data[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd) 
     775//                              data_err[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd) 
    672776//                              solidAngle[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd)           //testing only 
    673777                        endif 
     
    718822 
    719823 
    720 //distances passed in are in mm 
    721 // dtdist is SDD 
    722 // xd and yd are distances from the beam center to the current pixel 
    723 // 
    724 // TODO: 
    725 //   --         DoAlert 0,"This has not yet been updated for VSANS" 
    726 // 
    727 Function DetEffCorr(lambda,dtdist,xd,yd) 
    728         Variable lambda,dtdist,xd,yd 
    729  
    730         DoAlert 0,"This has not yet been updated for VSANS" 
    731          
    732         Variable theta,cosT,ff,stAl,stHe 
    733          
    734         theta = atan( (sqrt(xd^2 + yd^2))/dtdist ) 
    735         cosT = cos(theta) 
    736          
    737         stAl = 0.00967*lambda*0.8               //dimensionless, constants from JGB memo 
    738         stHe = 0.146*lambda*2.5 
    739          
    740         ff = exp(-stAl/cosT)*(1-exp(-stHe/cosT)) / ( exp(-stAl)*(1-exp(-stHe)) ) 
    741                  
    742         return(ff) 
    743 End 
    744824 
    745825// DIVIDE the intensity by this correction to get the right answer 
     
    10311111         
    10321112        Variable err 
    1033         err = DIVCorrection(type)               //returns err = 1 if data doesn't exist in specified folders 
     1113        err = V_DIVCorrection(type)             //returns err = 1 if data doesn't exist in specified folders 
    10341114         
    10351115        if(err) 
    1036                 Abort "error in DIVCorrection()" 
     1116                Abort "error in V_DIVCorrection()" 
    10371117        endif 
    10381118         
     
    10561136// TODO: 
    10571137//   --         DoAlert 0,"This has not yet been updated for VSANS" 
    1058 //  -- how is the error propagation handled? 
     1138//   -- how is the error propagation handled? 
    10591139// 
    10601140//function will divide the contents of "workType" folder with the contents of  
     
    10621142// all data is linear scale for the calculation 
    10631143// 
    1064 Function DIVCorrection(data,data_err,detStr,workType) 
     1144Function V_DIVCorrection(data,data_err,detStr,workType) 
    10651145        Wave data,data_err 
    10661146        String detStr,workType 
     
    10711151 
    10721152        if(WaveExists(data) == 0) 
    1073                 Print "The data wave does not exist in DIVCorrection()" 
     1153                Print "The data wave does not exist in V_DIVCorrection()" 
    10741154                Return(1)               //error condition 
    10751155        Endif 
     
    10801160        WAVE/Z div_data = $("root:Packages:NIST:VSANS:DIV:entry:instrument:detector_"+detStr+":data") 
    10811161        if(WaveExists(div_data) == 0) 
    1082                 Print "The DIV wave does not exist in DIVCorrection()" 
     1162                Print "The DIV wave does not exist in V_DIVCorrection()" 
    10831163                Return(1)               //error condition 
    10841164        Endif 
Note: See TracChangeset for help on using the changeset viewer.