 Timestamp:
 Jan 14, 2020 2:15:48 PM (3 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_DetectorCorrections.ipf
r1147 r1235 951 951 // (YES just from geometry, since I need SDD and dx and dy values...) 952 952 // 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 // 953 962 // 954 963 Function V_SolidAngleCorrection(w,w_err,fname,detStr,destPath) … … 957 966 958 967 Variable sdd,xCtr,yCtr,lambda 959 968 String orientation 969 960 970 // get all of the geometry information 961 //orientation = V_getDet_tubeOrientation(fname,detStr)971 orientation = V_getDet_tubeOrientation(fname,detStr) 962 972 sdd = V_getDet_ActualDistance(fname,detStr) 963 964 973 965 974 // this is ctr in mm … … 973 982 Wave data_realDistY = data_realDistY 974 983 975 Duplicate/O w solid_angle,tmp_theta,tmp_dist 984 Duplicate/O w solid_angle,tmp_theta,tmp_dist,tmp_theta_i //in the current df 976 985 977 986 //// calculate the scattering angle … … 983 992 // sdd is in [cm] 984 993 985 tmp_theta = atan(tmp_dist/sdd) //this is two_theta, the scattering angle994 tmp_theta = atan(tmp_dist/sdd) //this is two_theta, the (total) scattering angle 986 995 987 996 Variable ii,jj,numx,numy,dx,dy 988 997 numx = DimSize(tmp_theta,0) 989 998 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[ii1][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][jj1]) //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 ydirection 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[ii1][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][jj1]) //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) 1104 end 1105 1106 // this is the incorrect solid angle correction that does not take into 1107 // account the tube geometry. It is correct for the highres detector (and the 30m Ordela) 1108 // 1109 //  only for testing to prove that the cos(th)^2 *cos(th_i) is correct 1110 // 1111 Function 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 991 1148 for(ii=0 ;ii<numx;ii+=1) 992 1149 for(jj=0;jj<numy;jj+=1) … … 1021 1178 // Here it is! Apply the correction to the intensity (I divide  to get the counts per solid angle!!) 1022 1179 w /= solid_angle 1023 1024 // 1180 1025 1181 // correctly apply the correction to the error wave (assume a perfect value?) 1026 1182 w_err /= solid_angle // 1183 1027 1184 1028 1185 // 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 1030 1187 1031 1188 return(0) 1032 1189 end 1033 1034 1190 1035 1191 … … 1462 1618 1463 1619 // 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 xdimension of 680 pixels. Don't do the shift 1621 // if the dimensions are incorrect. 1622 if(DimSize(w,0) < 680) 1466 1623 adjW=w 1467 1624 return(0)
Note: See TracChangeset
for help on using the changeset viewer.