Ignore:
Timestamp:
Jan 14, 2020 2:15:48 PM (3 years ago)
Author:
srkline
Message:

Changed the calculation of solid angle to properly reflect the tube geometry. This mode is now the (correct) default calculation. A global switch can be applied for testing to use the old method (with many warning alerts).

File:
1 edited

Legend:

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

    r1147 r1235  
    951951//    (YES just from geometry, since I need SDD and dx and dy values...) 
    952952// 
     953//              "B" is passed, so I need to check for "B" and for panel orientation 
     954// -if it is detector "B" (not tubes), then the normal solid angle correction applies 
     955// -if it is a tube panel, then I need to know the orientation, to know which angles 
     956//    and pixel dimensions to use 
     957// 
     958// *** UPDATED 1/2020 SRK 
     959// -using new calculation since the lateral direction of the tubes does not affect the solid angle 
     960// projection (see He (2015) and John's memo) 
     961// 
    953962// 
    954963Function V_SolidAngleCorrection(w,w_err,fname,detStr,destPath) 
     
    957966 
    958967        Variable sdd,xCtr,yCtr,lambda 
    959  
     968        String orientation 
     969         
    960970// get all of the geometry information   
    961 //      orientation = V_getDet_tubeOrientation(fname,detStr) 
     971        orientation = V_getDet_tubeOrientation(fname,detStr) 
    962972        sdd = V_getDet_ActualDistance(fname,detStr) 
    963  
    964973 
    965974        // this is ctr in mm 
     
    973982        Wave data_realDistY = data_realDistY 
    974983 
    975         Duplicate/O w solid_angle,tmp_theta,tmp_dist            //in the current df 
     984        Duplicate/O w solid_angle,tmp_theta,tmp_dist,tmp_theta_i        //in the current df 
    976985 
    977986//// calculate the scattering angle 
     
    983992        // sdd is in [cm] 
    984993 
    985         tmp_theta = atan(tmp_dist/sdd)          //this is two_theta, the scattering angle 
     994        tmp_theta = atan(tmp_dist/sdd)          //this is two_theta, the (total) scattering angle 
    986995 
    987996        Variable ii,jj,numx,numy,dx,dy 
    988997        numx = DimSize(tmp_theta,0) 
    989998        numy = DimSize(tmp_theta,1) 
    990          
     999 
     1000        if(cmpstr(detStr,"B")==0) 
     1001                //detector B is a grid, straightforward cos^3 solid angle 
     1002                for(ii=0        ;ii<numx;ii+=1) 
     1003                        for(jj=0;jj<numy;jj+=1) 
     1004                                 
     1005                                if(ii==0)               //do a forward difference if ii==0 
     1006                                        dx = (data_realDistX[ii+1][jj] - data_realDistX[ii][jj])        //delta x for the pixel 
     1007                                else 
     1008                                        dx = (data_realDistX[ii][jj] - data_realDistX[ii-1][jj])        //delta x for the pixel 
     1009                                endif 
     1010                                 
     1011                                 
     1012                                if(jj==0) 
     1013                                        dy = (data_realDistY[ii][jj+1] - data_realDistY[ii][jj])        //delta y for the pixel 
     1014                                else 
     1015                                        dy = (data_realDistY[ii][jj] - data_realDistY[ii][jj-1])        //delta y for the pixel 
     1016                                endif 
     1017                 
     1018                                dx /= 10 
     1019                                dy /= 10                // convert mm to cm (since sdd is in cm) 
     1020                                solid_angle[ii][jj] = dx*dy             //this is in cm^2 
     1021                        endfor 
     1022                endfor 
     1023                 
     1024                // to cover up any issues w/negative dx or dy 
     1025                solid_angle = abs(solid_angle) 
     1026                 
     1027                // solid_angle correction 
     1028                // == dx*dy*cos^3/sdd^2 
     1029                solid_angle *= (cos(tmp_theta))^3 
     1030                solid_angle /= sdd^2 
     1031                 
     1032                // Here it is! Apply the correction to the intensity (I divide -- to get the counts per solid angle!!) 
     1033                w /= solid_angle 
     1034 
     1035                // correctly apply the correction to the error wave (assume a perfect value?) 
     1036                w_err /= solid_angle            // 
     1037 
     1038        else 
     1039                //               
     1040                //different calculation for the tubes, different calculation based on XY orientation 
     1041                // 
     1042                if(cmpstr(orientation,"vertical")==0) 
     1043                        // L/R panels, tube axis is y-direction 
     1044                        // this is now a different tmp_dist 
     1045                        // convert everything to cm first! 
     1046                        // sdd is in [cm], everything else is in [mm] 
     1047                        tmp_dist = (data_realDistY/10 - yctr/10)/sqrt((data_realDistX/10 - xctr/10)^2 + sdd^2)           
     1048                        tmp_theta_i = atan(tmp_dist)            //this is theta_y 
     1049                         
     1050                else 
     1051                        // horizontal orientation (T/B panels) 
     1052                        // this is now a different tmp_dist 
     1053                        // convert everything to cm first! 
     1054                        // sdd is in [cm], everything else is in [mm] 
     1055                        tmp_dist = (data_realDistX/10 - xctr/10)/sqrt((data_realDistY/10 - yctr/10)^2 + sdd^2)           
     1056                        tmp_theta_i = atan(tmp_dist)            //this is theta_x 
     1057                 
     1058                endif 
     1059                 
     1060                for(ii=0        ;ii<numx;ii+=1) 
     1061                        for(jj=0;jj<numy;jj+=1) 
     1062                                 
     1063                                if(ii==0)               //do a forward difference if ii==0 
     1064                                        dx = (data_realDistX[ii+1][jj] - data_realDistX[ii][jj])        //delta x for the pixel 
     1065                                else 
     1066                                        dx = (data_realDistX[ii][jj] - data_realDistX[ii-1][jj])        //delta x for the pixel 
     1067                                endif 
     1068                                 
     1069                                 
     1070                                if(jj==0) 
     1071                                        dy = (data_realDistY[ii][jj+1] - data_realDistY[ii][jj])        //delta y for the pixel 
     1072                                else 
     1073                                        dy = (data_realDistY[ii][jj] - data_realDistY[ii][jj-1])        //delta y for the pixel 
     1074                                endif 
     1075                 
     1076                                dx /= 10 
     1077                                dy /= 10                // convert mm to cm (since sdd is in cm) 
     1078                                solid_angle[ii][jj] = dx*dy             //this is in cm^2 
     1079                        endfor 
     1080                endfor 
     1081                 
     1082                // to cover up any issues w/negative dx or dy 
     1083                solid_angle = abs(solid_angle) 
     1084                 
     1085                // solid_angle correction 
     1086                // == dx*dy*cos(th)^2*cos(th_i)/sdd^2           using either the theta_x or theta_y value 
     1087                solid_angle *= (cos(tmp_theta))^2*cos(tmp_theta_i) 
     1088                solid_angle /= sdd^2 
     1089                 
     1090                // Here it is! Apply the correction to the intensity (I divide -- to get the counts per solid angle!!) 
     1091                w /= solid_angle 
     1092                 
     1093                // 
     1094                // correctly apply the correction to the error wave (assume a perfect value?) 
     1095                w_err /= solid_angle            // 
     1096         
     1097        endif 
     1098         
     1099 
     1100// DONE x- clean up after I'm satisfied computations are correct                 
     1101        KillWaves/Z tmp_theta,tmp_dist,tmp_theta_i 
     1102         
     1103        return(0) 
     1104end 
     1105 
     1106// this is the incorrect solid angle correction that does not take into  
     1107// account the tube geometry. It is correct for the high-res detector (and the 30m Ordela) 
     1108// 
     1109// -- only for testing to prove that the cos(th)^2 *cos(th_i) is correct 
     1110// 
     1111Function V_SolidAngleCorrection_COS3(w,w_err,fname,detStr,destPath) 
     1112        Wave w,w_err 
     1113        String fname,detStr,destPath 
     1114 
     1115        Variable sdd,xCtr,yCtr,lambda 
     1116        String orientation 
     1117         
     1118// get all of the geometry information   
     1119        orientation = V_getDet_tubeOrientation(fname,detStr) 
     1120        sdd = V_getDet_ActualDistance(fname,detStr) 
     1121 
     1122        // this is ctr in mm 
     1123        xCtr = V_getDet_beam_center_x_mm(fname,detStr) 
     1124        yCtr = V_getDet_beam_center_y_mm(fname,detStr) 
     1125        lambda = V_getWavelength(fname) 
     1126         
     1127        SetDataFolder $(destPath + ":entry:instrument:detector_"+detStr) 
     1128         
     1129        Wave data_realDistX = data_realDistX 
     1130        Wave data_realDistY = data_realDistY 
     1131 
     1132        Duplicate/O w solid_angle,tmp_theta,tmp_dist    //in the current df 
     1133 
     1134//// calculate the scattering angle 
     1135//      dx = (distX - xctr)             //delta x in mm 
     1136//      dy = (distY - yctr)             //delta y in mm 
     1137        tmp_dist = sqrt((data_realDistX - xctr)^2 + (data_realDistY - yctr)^2) 
     1138         
     1139        tmp_dist /= 10  // convert mm to cm 
     1140        // sdd is in [cm] 
     1141 
     1142        tmp_theta = atan(tmp_dist/sdd)          //this is two_theta, the (total) scattering angle 
     1143 
     1144        Variable ii,jj,numx,numy,dx,dy 
     1145        numx = DimSize(tmp_theta,0) 
     1146        numy = DimSize(tmp_theta,1) 
     1147 
    9911148        for(ii=0        ;ii<numx;ii+=1) 
    9921149                for(jj=0;jj<numy;jj+=1) 
     
    10211178        // Here it is! Apply the correction to the intensity (I divide -- to get the counts per solid angle!!) 
    10221179        w /= solid_angle 
    1023          
    1024         // 
     1180 
    10251181        // correctly apply the correction to the error wave (assume a perfect value?) 
    10261182        w_err /= solid_angle            // 
     1183         
    10271184 
    10281185// DONE x- clean up after I'm satisfied computations are correct                 
    1029         KillWaves/Z tmp_theta,tmp_dist 
     1186        KillWaves/Z tmp_theta,tmp_dist,tmp_theta_i 
    10301187         
    10311188        return(0) 
    10321189end 
    1033  
    10341190 
    10351191 
     
    14621618 
    14631619// this is necessary for some old data with the 150x150 back (dummy) panel 
    1464         NVAR gIgnoreDetB = root:Packages:NIST:VSANS:Globals:gIgnoreDetB 
    1465         if(gIgnoreDetB == 1) 
     1620// the proper back detector has an x-dimension of 680 pixels. Don't do the shift 
     1621// if the dimensions are incorrect. 
     1622        if(DimSize(w,0) < 680) 
    14661623                adjW=w 
    14671624                return(0) 
Note: See TracChangeset for help on using the changeset viewer.