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).

Location:
sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS
Files:
4 edited

Legend:

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

    r1227 r1235  
    243243//      String reducedFolderType="COR",carriageStr="F" 
    244244//      Prompt reducedFolderType, "reduced data folder" 
    245         Prompt carriageStr,"panels to group",popup,"All 8;All F All M;Individual;B;" 
     245        Prompt carriageStr,"panels to group",popup,"Individual;B;All 8;All F All M;" 
    246246                 
    247247        Vf_NormalizeDIV_proc(carriageStr) 
  • 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) 
  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_Initialize.ipf

    r1234 r1235  
    275275        Variable/G root:Packages:NIST:VSANS:Globals:isDemoVersion = V_isDemo() 
    276276 
     277 
     278// for testing 
     279// if this is set to 1, the OLD (incorrect) cos^3 solid angle will be applied to the tubes 
     280// and a lot of alerts will pop up... 
     281        Variable/G root:Packages:NIST:VSANS:Globals:gDo_OLD_SolidAngleCor = 0 
    277282         
    278283        //set XML globals 
  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_WorkFolderUtils.ipf

    r1151 r1235  
    651651                endfor 
    652652 
    653                 if(gIgnoreDetB==1)               
     653                if(gIgnoreDetB==0)               
    654654                        //"B" is separate 
    655655                        detStr = "B" 
     
    738738        // the solid angle correction is calculated for ALL detector panels. 
    739739        NVAR gDoSolidAngleCor = root:Packages:NIST:VSANS:Globals:gDoSolidAngleCor 
     740        NVAR/Z gDo_OLD_SolidAngleCor = root:Packages:NIST:VSANS:Globals:gDo_OLD_SolidAngleCor 
     741        // for older experiments, this won't exist, so generate it and default to zero 
     742        // so the old calculation is not done 
     743        if(NVAR_Exists(gDo_OLD_SolidAngleCor)==0) 
     744                Variable/G root:Packages:NIST:VSANS:Globals:gDo_OLD_SolidAngleCor=0 
     745        endif 
    740746        if (gDoSolidAngleCor == 1) 
    741747                Print "Doing Solid Angle correction"// for "+ detStr 
     
    744750                         
    745751                        if(cmpstr(detStr,"B") == 0  && gIgnoreDetB == 1) 
    746                                 // do nothing 
     752                                // do nothing if the back is ignored 
    747753                        else 
    748754                                Wave w = V_getDetectorDataW(fname,detStr) 
    749755                                Wave w_err = V_getDetectorDataErrW(fname,detStr) 
    750756                                // any other dimensions to pass in? 
    751                                 V_SolidAngleCorrection(w,w_err,fname,detStr,destPath) 
     757                                 
     758                                if(gDo_OLD_SolidAngleCor == 0) 
     759                                        V_SolidAngleCorrection(w,w_err,fname,detStr,destPath) 
     760                                else 
     761                                        // for testing ONLY -- the cos^3 correction is incorrect for tubes, and the normal 
     762                                        // function call above   correctly handles either high-res grid or tubes. This COS3 function 
     763                                        // will incorrectly treat tubes as a grid        
     764                                        //                              Print "TESTING -- using incorrect COS^3 solid angle !"           
     765                                        V_SolidAngleCorrection_COS3(w,w_err,fname,detStr,destPath) 
     766                                endif 
     767                                 
    752768                        endif 
    753769                         
    754770                endfor 
     771                if(gDo_OLD_SolidAngleCor == 1) 
     772                        DoAlert 0,"TESTING -- using incorrect COS^3 solid angle !"               
     773                endif 
     774 
    755775        else 
    756776                Print "Solid Angle correction NOT DONE" 
     
    838858                        // TODO 
    839859                        // -- do I want to update and save the integrated detector count? 
    840                         Variable integratedCount = V_getDet_IntegratedCount(fname,detStr) 
    841                         V_writeDet_IntegratedCount(fname,detStr,integratedCount*scale) 
     860                        // I can't trust the integrated count for the back detector (ever) - since the  
     861                        // data is shifted for registry and some (~ 10% of the detector) is lost 
     862                        // also, ML panel has the wrong value (Oct-nov 2019) and I don't know why. so sum the data 
     863                        // 
     864                        Variable integratedCount = sum(w) 
     865                        V_writeDet_IntegratedCount(fname,detStr,integratedCount)                //already the scaled value for counts 
     866//                      Variable integratedCount = V_getDet_IntegratedCount(fname,detStr) 
     867//                      V_writeDet_IntegratedCount(fname,detStr,integratedCount*scale) 
    842868 
    843869                endfor 
     
    968994        V_putControlMonitorCount(newType,saved_mon_dest+saved_mon_tmp) 
    969995 
     996// TODO (DONE) 
     997// the new, unscaled monitor count was written to the control block, but it needs to be  
     998// written to the BeamMonNormSaved_count field instead, since this is where I read it from. 
     999// - so this worked in the past for adding two files, but fails on 3+ 
     1000// x- write to the NormSaved_count field... 
     1001        V_writeBeamMonNormSaved_count(newType,saved_mon_dest+saved_mon_tmp)                     // save the true count 
     1002 
    9701003 
    9711004// now loop over all of the detector panels 
     
    9951028                data_err_tmp *= scale_tmp 
    9961029                linear_data_tmp *= scale_tmp 
    997                                  
     1030 
     1031// TODO SRK-ERROR? 
     1032// the integrated count may not be correct (ML error Oct/Nov 2019) 
     1033// and is not correct for the back detector since some pixels were lost due to shifting the image registration 
     1034//                               
    9981035                // add them together, the dest is a wave so it is automatically changed in the "dest" folder 
    999                 V_putDet_IntegratedCount(tmpType,detStr,det_integrated_ct_dest+det_integrated_ct_tmp) 
     1036                V_putDet_IntegratedCount(tmpType,detStr,sum(data_dest)+sum(data_tmp))           // adds the unscaled data sums 
     1037//              V_putDet_IntegratedCount(tmpType,detStr,det_integrated_ct_dest+det_integrated_ct_tmp)           // wrong for "B", may be wrong for ML 
    10001038                data_dest += data_tmp 
    10011039                data_err_dest = sqrt(data_err_dest^2 + data_err_tmp^2)          // add in quadrature 
Note: See TracChangeset for help on using the changeset viewer.