- Timestamp:
- Dec 6, 2017 1:36:11 PM (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/VC_DetectorBinning_Slit.ipf
r992 r1076 24 24 25 25 26 // TODO: 27 // -- verify the error calculation 28 // -- add in functionality to handle FLR and MLR cases (2 panels of data) 26 29 // 27 30 // seems backwards to call this "byRows", but this is the way that Igor indexes … … 36 39 // SetDataFolder root:Packages:NIST:VSANS:VCALC 37 40 38 Variable pixSizeX,pixSizeY,delQx, delQy 39 Variable isVCALC=0 41 Variable pixSizeX,pixSizeY,delQy 42 Variable isVCALC=0,maskMissing,nsets=0 43 44 maskMissing = 1 // set to zero if a mask is actually present 40 45 41 46 if(cmpstr(folderStr,"VCALC") == 0) … … 45 50 String folderPath = "root:Packages:NIST:VSANS:"+folderStr 46 51 String instPath = ":entry:instrument:detector_" 47 48 if(isVCALC) 49 WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr) // 2D detector data 50 WAVE/Z iErr = $("asdf_iErr_"+detStr) // TODO: 2D errors -- may not exist, especially for simulation 51 else 52 Wave inten = V_getDetectorDataW(folderStr,detStr) 53 Wave iErr = V_getDetectorDataErrW(folderStr,detStr) 52 53 strswitch(detStr) // string switch 54 // only one panel, simply pick that panel and move on out of the switch 55 case "FL": 56 case "FR": 57 case "ML": 58 case "MR": 59 case "B": 60 if(isVCALC) 61 WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr) // 2D detector data 62 WAVE/Z iErr = $("asdf_iErr_"+detStr) // TODO: 2D errors -- may not exist, especially for simulation 63 else 64 Wave inten = V_getDetectorDataW(folderStr,detStr) 65 Wave iErr = V_getDetectorDataErrW(folderStr,detStr) 66 Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data") 67 if(WaveExists(mask) == 1) 68 maskMissing = 0 69 endif 70 endif 71 72 nsets = 1 73 break 74 75 // case "FLR": 76 // // detStr has multiple values now, so unfortuntely, I'm hard-wiring things... 77 // if(isVCALC) 78 // WAVE inten = $(folderPath+instPath+"FL"+":det_"+"FL") 79 // WAVE/Z iErr = $("iErr_"+"FL") // 2D errors -- may not exist, especially for simulation 80 // WAVE inten2 = $(folderPath+instPath+"FR"+":det_"+"FR") 81 // WAVE/Z iErr2 = $("iErr_"+"FR") // 2D errors -- may not exist, especially for simulation 82 // else 83 // Wave inten = V_getDetectorDataW(folderStr,"FL") 84 // Wave iErr = V_getDetectorDataErrW(folderStr,"FL") 85 // Wave inten2 = V_getDetectorDataW(folderStr,"FR") 86 // Wave iErr2 = V_getDetectorDataErrW(folderStr,"FR") 87 // Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FL"+":data") 88 // Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FR"+":data") 89 // if(WaveExists(mask) == 1 && WaveExists(mask2) == 1) 90 // maskMissing = 0 91 // endif 92 // endif 93 // 94 // nsets = 2 95 // break 96 // 97 // case "MLR": 98 // if(isVCALC) 99 // WAVE inten = $(folderPath+instPath+"ML"+":det_"+"ML") 100 // WAVE/Z iErr = $("iErr_"+"ML") // 2D errors -- may not exist, especially for simulation 101 // WAVE inten2 = $(folderPath+instPath+"MR"+":det_"+"MR") 102 // WAVE/Z iErr2 = $("iErr_"+"MR") // 2D errors -- may not exist, especially for simulation 103 // else 104 // Wave inten = V_getDetectorDataW(folderStr,"ML") 105 // Wave iErr = V_getDetectorDataErrW(folderStr,"ML") 106 // Wave inten2 = V_getDetectorDataW(folderStr,"MR") 107 // Wave iErr2 = V_getDetectorDataErrW(folderStr,"MR") 108 // Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"ML"+":data") 109 // Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MR"+":data") 110 // if(WaveExists(mask) == 1 && WaveExists(mask2) == 1) 111 // maskMissing = 0 112 // endif 113 // endif 114 // 115 // nSets = 2 116 // break 117 118 default: 119 nSets = 0 120 Print "ERROR ---- type is not recognized " 121 endswitch 122 123 if(nSets == 0) 124 SetDataFolder root: 125 return(0) 54 126 endif 127 55 128 56 129 Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr) // 2D q-values … … 58 131 Wave qy = $(folderPath+instPath+detStr+":qy_"+detStr) 59 132 60 // ?? TODO not needed here? 61 // pixSizeX = VCALC_getPixSizeX(detStr) 62 // pixSizeY = VCALC_getPixSizeY(detStr) 63 64 delQx = abs(qx[0][0] - qx[1][0]) 65 delQy = abs(qy[0][1] - qy[0][0]) 66 67 // delta Qx is set by the pixel X dimension of the detector, which is the limiting resolution 68 133 134 // delta Qx is set by the pixel X dimension of the detector, which is the limiting resolution 135 // delQx = abs(qx[0][0] - qx[1][0]) 136 69 137 Variable nq,val 70 138 nq = DimSize(inten,0) //nq == the number of columns (x dimension) … … 86 154 Wave eBin2D_qxqy = $(folderPath+":"+"eBin2D_qxqy_"+detStr) 87 155 156 157 iBin_qxqy = 0 158 iBin2_qxqy = 0 159 eBin_qxqy = 0 160 eBin2D_qxqy = 0 161 nBin_qxqy = 0 162 163 88 164 // sum the rows 89 MatrixOp/O iBin_qxqy = sumRows(inten) //automatically generates the destination 90 165 166 // MatrixOp would be fast, but I can't figure out how to apply the mask with sumRows??? 167 // MatrixOp/O iBin_qxqy = sumRows(inten) //automatically generates the destination 168 // 169 //// how to properly calculate the error? 170 //// This MatrixOp gives values way too large -- larger than the intensity (be sure to correct *delQy below... 171 // 172 // MatrixOp/O tmp = sqrt(varCols(inten^t)) // variance: no varRows operation, so use the transpose of the matrix 173 // eBin_qxqy = tmp[0][p] 174 175 Variable ii,jj,ntube,npix,sum_inten, sum_n, sum_inten2,avesq,aveisq,var,mask_val 176 ntube = DimSize(inten,0) 177 npix = DimSize(inten,1) 178 179 for(ii=0;ii<ntube;ii+=1) 180 sum_inten = 0 // initialize the sum 181 sum_n = 0 182 sum_inten2 = 0 183 184 for(jj=0;jj<npix;jj+=1) 185 val = inten[ii][jj] 186 if(isVCALC || maskMissing) // mask_val == 0 == keep, mask_val == 1 = YES, mask out the point 187 mask_val = 0 188 else 189 mask_val = mask[ii][jj] 190 endif 191 if (numType(val)==0 && mask_val == 0) //count only the good points, ignore Nan or Inf 192 sum_inten += val 193 sum_n += 1 194 sum_inten2 += val*val 195 endif 196 endfor 197 iBin_qxqy[ii] = sum_inten 198 199 avesq = sum_inten 200 aveisq = sum_inten2/sum_n 201 var = aveisq-avesq 202 if(var<=0) 203 eBin_qxqy[ii] = 1e-6 204 else 205 eBin_qxqy[ii] = sqrt(var/(sum_n - 1)) 206 endif 207 208 endfor 209 210 qBin_qxqy = qx[p][npix/2] //TODO: use only the Qx component in the y-center of the detector, not Qtotal 211 91 212 // if the detectors are "L", then the values are all negative... 92 // if the detectors are T/B, then half is negative, and there's a very nearly zero point in the middle...93 // and it may make no sense to use T/B anyways...94 213 qBin_qxqy = abs(qx[p][0]) 95 214 215 // for the L panels, sort the q-values (and data) after the abs() step, otherwise the data is reversed 216 // won't hurt to sort the R data 217 Sort qBin_qxqy, qBin_qxqy,iBin_qxqy,eBin_qxqy 96 218 97 219 //now get the scaling correct 98 220 // q-integration (rectangular), matrixOp simply summed, so I need to multiply by dy (pixelSizeY -> as Qy?) 99 221 222 // delQy = abs(qy[0][1] - qy[0][0]) // this is only one pixel 223 delQy = abs(qy[0][npix-1] - qy[0][0]) // TODO: do I use dQ for the height of the panel? 224 100 225 iBin_qxqy *= delQy 101 226 eBin_qxqy *= delQy 227 228 229 /// TODO -- this is not necessary, but just for getting the I(Q) display to look "pretty" 230 // clean out the near-zero Q point in the T/B and Back detectors 231 qBin_qxqy = (abs(qBin_qxqy[p][q]) < 1e-5) ? NaN : qBin_qxqy[p][q] 232 233 234 // clear out zero data values before exiting... 235 // find the last non-zero point, working backwards 236 val = numpnts(iBin_qxqy) 237 do 238 val -= 1 239 while((iBin_qxqy[val] == 0) && val > 0) 240 241 DeletePoints val, nq-val, iBin_qxqy,qBin_qxqy,eBin_qxqy 242 243 // work forwards? this doesn't work... 244 // val = -1 245 // do 246 // val += 1 247 // while(nBin_qxqy[val] == 0 && val < numpnts(iBin_qxqy)-1) 248 // DeletePoints 0, val, iBin_qxqy,qBin_qxqy,eBin_qxqy 249 250 251 // TODO: 252 // -- calculate the slit resolution here. don't really have any idea how to represent this in VSANS 253 // since it's not an infinite slit, and not anything that I expect could be represented as a Gaussian 254 255 // TODO: 256 // -- This is where I calculate the resolution in SANS (see CircSectAve) 257 // -- use the isVCALC flag to exclude VCALC from the resolution calculation if necessary 258 // -- from the top of the function, folderStr = work folder, detStr = "FLRTB" or other type of averaging 259 // 260 261 nq = numpnts(qBin_qxqy) 262 Make/O/D/N=(nq) $(folderPath+":"+"sigmaQ_qxqy"+"_"+detStr) 263 Make/O/D/N=(nq) $(folderPath+":"+"qBar_qxqy"+"_"+detStr) 264 Make/O/D/N=(nq) $(folderPath+":"+"fSubS_qxqy"+"_"+detStr) 265 Wave sigmaq = $(folderPath+":"+"sigmaQ_qxqy_"+detStr) 266 Wave qbar = $(folderPath+":"+"qBar_qxqy_"+detStr) 267 Wave fsubs = $(folderPath+":"+"fSubS_qxqy_"+detStr) 268 102 269 // TODO 103 // iBin_qxqy *= 4 //why the factor of 4??? -- this is what I needed to do with FFT->USANS. Do I need it here? 104 105 106 /// TODO -- this is not correct, but just for getting the I(Q) display to look "pretty" 107 qBin_qxqy = (abs(qBin_qxqy[p][q]) < 1e-5) ? NaN : qBin_qxqy[p][q] // clean out the near-zero Q point in the T/B and Back detectors 270 // -- these are DUMMY VALUES!!! 271 sigmaq = -delQy 272 qbar = -delQy 273 fsubs = -delQy 108 274 109 275
Note: See TracChangeset
for help on using the changeset viewer.
