| 1 | #pragma TextEncoding = "MacRoman" |
|---|
| 2 | #pragma rtGlobals=3 // Use modern global access method and strict wave access. |
|---|
| 3 | |
|---|
| 4 | // |
|---|
| 5 | // routines to do a sector average |
|---|
| 6 | // |
|---|
| 7 | // possible ideas: |
|---|
| 8 | // - start with calculation of the phi matrix for each panel |
|---|
| 9 | // - with the definition of the angle +/-, I can decide which points to keep during the average. |
|---|
| 10 | // - I can also make this an actual mask (sum with the protocol mask) and use it this way |
|---|
| 11 | // so that I can still run it through the circular average routine. |
|---|
| 12 | // - I can also then overlay the sector mask onto the data (once I figure out how to overlay masks on |
|---|
| 13 | // the regular data display |
|---|
| 14 | |
|---|
| 15 | |
|---|
| 16 | |
|---|
| 17 | Function/WAVE V_MakePhiMatrix(qTotal,folderStr,detStr,folderPath) |
|---|
| 18 | Wave qTotal |
|---|
| 19 | String folderStr,detStr,folderPath |
|---|
| 20 | |
|---|
| 21 | |
|---|
| 22 | Variable xctr = V_getDet_beam_center_x_pix(folderStr,detStr) |
|---|
| 23 | Variable yctr = V_getDet_beam_center_y_pix(folderStr,detStr) |
|---|
| 24 | |
|---|
| 25 | Duplicate/O qTotal,$(folderPath+":phi") |
|---|
| 26 | Wave phi = $(folderPath+":phi") |
|---|
| 27 | Variable pixSizeX,pixSizeY |
|---|
| 28 | pixSizeX = V_getDet_x_pixel_size(folderStr,detStr) |
|---|
| 29 | pixSizeY = V_getDet_y_pixel_size(folderStr,detStr) |
|---|
| 30 | phi = V_FindPhi( pixSizeX*((p+1)-xctr) , pixSizeY*((q+1)-yctr)) //(dx,dy) |
|---|
| 31 | |
|---|
| 32 | return phi |
|---|
| 33 | End |
|---|
| 34 | |
|---|
| 35 | |
|---|
| 36 | // |
|---|
| 37 | // x- I want to mask out everything that is "out" of the sector |
|---|
| 38 | // |
|---|
| 39 | // 0 = keep the point |
|---|
| 40 | // 1 = yes, mask the point |
|---|
| 41 | // |
|---|
| 42 | // |
|---|
| 43 | // phiCtr is in the range (-90,90) degrees |
|---|
| 44 | // delta is in the range (0,90) for a total width of 2*delta = 180 degrees |
|---|
| 45 | // |
|---|
| 46 | Function V_MarkSectorOverlayPixels(phi,overlay,phiCtr,delta,side) |
|---|
| 47 | Wave phi,overlay |
|---|
| 48 | Variable phiCtr,delta |
|---|
| 49 | String side |
|---|
| 50 | |
|---|
| 51 | Variable phiVal |
|---|
| 52 | |
|---|
| 53 | // convert the imput from degrees to radians , since phi is in radians |
|---|
| 54 | phiCtr *= pi/180 |
|---|
| 55 | delta *= pi/180 |
|---|
| 56 | |
|---|
| 57 | Variable xDim=DimSize(phi, 0) |
|---|
| 58 | Variable yDim=DimSize(phi, 1) |
|---|
| 59 | |
|---|
| 60 | Variable ii,jj,exclude,mirror_phiCtr,crossZero,keepPix |
|---|
| 61 | |
|---|
| 62 | // initialize the mask to == 1 == exclude everything |
|---|
| 63 | overlay = 1 |
|---|
| 64 | |
|---|
| 65 | // now give every opportunity to keep pixel in |
|---|
| 66 | // comparisons use a modified phiCtr to match the definition of the phi field (0= +x-axis) |
|---|
| 67 | // |
|---|
| 68 | for(ii=0;ii<xDim;ii+=1) |
|---|
| 69 | for(jj=0;jj<yDim;jj+=1) |
|---|
| 70 | //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2) |
|---|
| 71 | phiVal = phi[ii][jj] |
|---|
| 72 | keepPix = 0 //start with not keeping |
|---|
| 73 | |
|---|
| 74 | // if within the right or left, flag to keep the pixel |
|---|
| 75 | if(cmpstr(side,"right")==0) |
|---|
| 76 | //right, when 0->pi/2 |
|---|
| 77 | if(V_CloseEnough(phiVal,phiCtr,delta)) |
|---|
| 78 | keepPix = 1 |
|---|
| 79 | endif |
|---|
| 80 | // condition here to get the 3pi/2 -> 2pi region |
|---|
| 81 | if(V_CloseEnough(phiVal,phiCtr+2*pi,delta)) |
|---|
| 82 | keepPix = 1 |
|---|
| 83 | endif |
|---|
| 84 | endif |
|---|
| 85 | |
|---|
| 86 | if(cmpstr(side,"left")==0) |
|---|
| 87 | if(V_CloseEnough(phiVal,phiCtr+pi,delta)) |
|---|
| 88 | keepPix = 1 |
|---|
| 89 | endif |
|---|
| 90 | endif |
|---|
| 91 | |
|---|
| 92 | // both sides, duplicates the conditions above |
|---|
| 93 | if(cmpstr(side,"both")==0) |
|---|
| 94 | //right, when 0->pi/2 |
|---|
| 95 | if(V_CloseEnough(phiVal,phiCtr,delta)) |
|---|
| 96 | keepPix = 1 |
|---|
| 97 | endif |
|---|
| 98 | // right, when 3pi/2 -> 2pi |
|---|
| 99 | if(V_CloseEnough(phiVal,phiCtr+2*pi,delta)) |
|---|
| 100 | keepPix = 1 |
|---|
| 101 | endif |
|---|
| 102 | |
|---|
| 103 | //left |
|---|
| 104 | if(V_CloseEnough(phiVal,phiCtr+pi,delta)) |
|---|
| 105 | keepPix = 1 |
|---|
| 106 | endif |
|---|
| 107 | |
|---|
| 108 | endif |
|---|
| 109 | |
|---|
| 110 | // set the mask value (entire overlay initialized to 1 to start) |
|---|
| 111 | if(keepPix > 0) |
|---|
| 112 | overlay[ii][jj] = 0 |
|---|
| 113 | endif |
|---|
| 114 | |
|---|
| 115 | endfor |
|---|
| 116 | endfor |
|---|
| 117 | |
|---|
| 118 | |
|---|
| 119 | return(0) |
|---|
| 120 | End |
|---|
| 121 | |
|---|
| 122 | |
|---|
| 123 | |
|---|
| 124 | // |
|---|
| 125 | // TODO -- binType == 4 (slit mode) should never end up here |
|---|
| 126 | // -- new logic in calling routines to dispatch to proper routine |
|---|
| 127 | // -- AND need to write the routine for binning_SlitMode |
|---|
| 128 | // |
|---|
| 129 | // side = one of "left;right;both;" |
|---|
| 130 | // phi_rad = center of sector in radians |
|---|
| 131 | // dphi_rad = half-width of sector, also in radians |
|---|
| 132 | Function V_QBinAllPanels_Sector(folderStr,binType,collimationStr,side,phi_rad,dphi_rad) |
|---|
| 133 | String folderStr |
|---|
| 134 | Variable binType |
|---|
| 135 | String collimationStr,side |
|---|
| 136 | Variable phi_rad,dphi_rad |
|---|
| 137 | |
|---|
| 138 | // do the back, middle, and front separately |
|---|
| 139 | |
|---|
| 140 | // figure out the binning type (where is it set?) |
|---|
| 141 | Variable ii,delQ |
|---|
| 142 | String detStr |
|---|
| 143 | |
|---|
| 144 | // binType = V_GetBinningPopMode() |
|---|
| 145 | |
|---|
| 146 | // set delta Q for binning (used later inside VC_fDoBinning_QxQy2D) |
|---|
| 147 | for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1) |
|---|
| 148 | detStr = StringFromList(ii, ksDetectorListAll, ";") |
|---|
| 149 | |
|---|
| 150 | delQ = V_SetDeltaQ(folderStr,detStr) // this sets (overwrites) the global value for each panel type |
|---|
| 151 | endfor |
|---|
| 152 | |
|---|
| 153 | |
|---|
| 154 | switch(binType) |
|---|
| 155 | case 1: |
|---|
| 156 | V_fDoSectorBin_QxQy2D(folderStr,"FL",collimationStr,side,phi_rad,dphi_rad) |
|---|
| 157 | V_fDoSectorBin_QxQy2D(folderStr,"FR",collimationStr,side,phi_rad,dphi_rad) |
|---|
| 158 | V_fDoSectorBin_QxQy2D(folderStr,"FT",collimationStr,side,phi_rad,dphi_rad) |
|---|
| 159 | V_fDoSectorBin_QxQy2D(folderStr,"FB",collimationStr,side,phi_rad,dphi_rad) |
|---|
| 160 | V_fDoSectorBin_QxQy2D(folderStr,"ML",collimationStr,side,phi_rad,dphi_rad) |
|---|
| 161 | V_fDoSectorBin_QxQy2D(folderStr,"MR",collimationStr,side,phi_rad,dphi_rad) |
|---|
| 162 | V_fDoSectorBin_QxQy2D(folderStr,"MT",collimationStr,side,phi_rad,dphi_rad) |
|---|
| 163 | V_fDoSectorBin_QxQy2D(folderStr,"MB",collimationStr,side,phi_rad,dphi_rad) |
|---|
| 164 | V_fDoSectorBin_QxQy2D(folderStr, "B",collimationStr,side,phi_rad,dphi_rad) |
|---|
| 165 | |
|---|
| 166 | break |
|---|
| 167 | case 2: |
|---|
| 168 | V_fDoSectorBin_QxQy2D(folderStr,"FLR",collimationStr,side,phi_rad,dphi_rad) |
|---|
| 169 | V_fDoSectorBin_QxQy2D(folderStr,"FTB",collimationStr,side,phi_rad,dphi_rad) |
|---|
| 170 | V_fDoSectorBin_QxQy2D(folderStr,"MLR",collimationStr,side,phi_rad,dphi_rad) |
|---|
| 171 | V_fDoSectorBin_QxQy2D(folderStr,"MTB",collimationStr,side,phi_rad,dphi_rad) |
|---|
| 172 | V_fDoSectorBin_QxQy2D(folderStr, "B",collimationStr,side,phi_rad,dphi_rad) |
|---|
| 173 | |
|---|
| 174 | break |
|---|
| 175 | case 3: |
|---|
| 176 | V_fDoSectorBin_QxQy2D(folderStr,"MLRTB",collimationStr,side,phi_rad,dphi_rad) |
|---|
| 177 | V_fDoSectorBin_QxQy2D(folderStr,"FLRTB",collimationStr,side,phi_rad,dphi_rad) |
|---|
| 178 | V_fDoSectorBin_QxQy2D(folderStr, "B",collimationStr,side,phi_rad,dphi_rad) |
|---|
| 179 | |
|---|
| 180 | break |
|---|
| 181 | case 4: /// this is for a tall, narrow slit mode |
|---|
| 182 | VC_fBinDetector_byRows(folderStr,"FL") |
|---|
| 183 | VC_fBinDetector_byRows(folderStr,"FR") |
|---|
| 184 | VC_fBinDetector_byRows(folderStr,"ML") |
|---|
| 185 | VC_fBinDetector_byRows(folderStr,"MR") |
|---|
| 186 | VC_fBinDetector_byRows(folderStr,"B") |
|---|
| 187 | |
|---|
| 188 | break |
|---|
| 189 | case 5: |
|---|
| 190 | V_fDoSectorBin_QxQy2D(folderStr,"FTB",collimationStr,side,phi_rad,dphi_rad) |
|---|
| 191 | V_fDoSectorBin_QxQy2D(folderStr,"FLR",collimationStr,side,phi_rad,dphi_rad) |
|---|
| 192 | V_fDoSectorBin_QxQy2D(folderStr,"MLRTB",collimationStr,side,phi_rad,dphi_rad) |
|---|
| 193 | V_fDoSectorBin_QxQy2D(folderStr, "B",collimationStr,side,phi_rad,dphi_rad) |
|---|
| 194 | |
|---|
| 195 | break |
|---|
| 196 | case 6: |
|---|
| 197 | V_fDoSectorBin_QxQy2D(folderStr,"FLRTB",collimationStr,side,phi_rad,dphi_rad) |
|---|
| 198 | V_fDoSectorBin_QxQy2D(folderStr,"MLR",collimationStr,side,phi_rad,dphi_rad) |
|---|
| 199 | V_fDoSectorBin_QxQy2D(folderStr, "B",collimationStr,side,phi_rad,dphi_rad) |
|---|
| 200 | |
|---|
| 201 | break |
|---|
| 202 | case 7: |
|---|
| 203 | V_fDoSectorBin_QxQy2D(folderStr,"FTB",collimationStr,side,phi_rad,dphi_rad) |
|---|
| 204 | V_fDoSectorBin_QxQy2D(folderStr,"FLR",collimationStr,side,phi_rad,dphi_rad) |
|---|
| 205 | V_fDoSectorBin_QxQy2D(folderStr,"MLR",collimationStr,side,phi_rad,dphi_rad) |
|---|
| 206 | V_fDoSectorBin_QxQy2D(folderStr, "B",collimationStr,side,phi_rad,dphi_rad) |
|---|
| 207 | |
|---|
| 208 | break |
|---|
| 209 | |
|---|
| 210 | default: |
|---|
| 211 | Abort "Binning mode not found in V_QBinAllPanels_Circular"// when no case matches |
|---|
| 212 | endswitch |
|---|
| 213 | |
|---|
| 214 | |
|---|
| 215 | return(0) |
|---|
| 216 | End |
|---|
| 217 | |
|---|
| 218 | |
|---|
| 219 | ////////// |
|---|
| 220 | // |
|---|
| 221 | // Function that bins a 2D detctor panel into I(q) based on the q-value of the pixel |
|---|
| 222 | // - each pixel QxQyQz has been calculated beforehand |
|---|
| 223 | // - if multiple panels are selected to be combined, it is done here during the binning |
|---|
| 224 | // - the setting of deltaQ step is still a little suspect (TODO) |
|---|
| 225 | // |
|---|
| 226 | // |
|---|
| 227 | // see the equivalent function in PlotUtils2D_v40.ipf |
|---|
| 228 | // |
|---|
| 229 | //Function fDoBinning_QxQy2D(inten,qx,qy,qz) |
|---|
| 230 | // |
|---|
| 231 | // this has been modified to accept different detector panels and to take arrays |
|---|
| 232 | // -- type = FL or FR or...other panel identifiers |
|---|
| 233 | // |
|---|
| 234 | // TODO "iErr" is not always defined correctly since it doesn't really apply here for data that is not 2D simulation |
|---|
| 235 | // |
|---|
| 236 | // |
|---|
| 237 | // updated Feb2016 to take new folder structure |
|---|
| 238 | // TODO |
|---|
| 239 | // -- VERIFY |
|---|
| 240 | // -- figure out what the best location is to put the averaged data? currently @ top level of WORK folder |
|---|
| 241 | // but this is a lousy choice. |
|---|
| 242 | // x- binning is now Mask-aware. If mask is not present, all data is used. If data is from VCALC, all data is used |
|---|
| 243 | // x- Where do I put the solid angle correction? In here as a weight for each point, or later on as |
|---|
| 244 | // a blanket correction (matrix multiply) for an entire panel? (Solid Angle correction is done in the |
|---|
| 245 | // step where data is added to a WORK file (see Raw_to_Work()) |
|---|
| 246 | // |
|---|
| 247 | // |
|---|
| 248 | // TODO: |
|---|
| 249 | // -- some of the input parameters for the resolution calcuation are either assumed (apOff) or are currently |
|---|
| 250 | // hard-wired. these need to be corrected before even the pinhole resolution is correct |
|---|
| 251 | // x- resolution calculation is in the correct place. The calculation is done per-panel (specified by TYPE), |
|---|
| 252 | // and then the unwanted points can be discarded (all 6 columns) as the data is trimmed and concatenated |
|---|
| 253 | // is separate functions that are resolution-aware. |
|---|
| 254 | // |
|---|
| 255 | // |
|---|
| 256 | // folderStr = WORK folder, type = the binning type (may include multiple detectors) |
|---|
| 257 | // |
|---|
| 258 | // side = one of "left;right;both;" |
|---|
| 259 | // phi_rad = center of sector in radians |
|---|
| 260 | // dphi_rad = half-width of sector, also in radians |
|---|
| 261 | // |
|---|
| 262 | Function V_fDoSectorBin_QxQy2D(folderStr,type,collimationStr,side,phi_rad,dphi_rad) |
|---|
| 263 | String folderStr,type,collimationStr,side |
|---|
| 264 | Variable phi_rad,dphi_rad |
|---|
| 265 | |
|---|
| 266 | Variable nSets = 0 |
|---|
| 267 | Variable xDim,yDim |
|---|
| 268 | Variable ii,jj |
|---|
| 269 | Variable qVal,nq,var,avesq,aveisq |
|---|
| 270 | Variable binIndex,val,isVCALC=0,maskMissing |
|---|
| 271 | |
|---|
| 272 | String folderPath = "root:Packages:NIST:VSANS:"+folderStr |
|---|
| 273 | String instPath = ":entry:instrument:detector_" |
|---|
| 274 | String detStr |
|---|
| 275 | |
|---|
| 276 | if(cmpstr(folderStr,"VCALC") == 0) |
|---|
| 277 | isVCALC = 1 |
|---|
| 278 | endif |
|---|
| 279 | |
|---|
| 280 | detStr = type |
|---|
| 281 | |
|---|
| 282 | // now switch on the type to determine which waves to declare and create |
|---|
| 283 | // since there may be more than one panel to step through. There may be two, there may be four |
|---|
| 284 | // |
|---|
| 285 | |
|---|
| 286 | // TODO: |
|---|
| 287 | // -- Solid_Angle -- waves will be present for WORK data other than RAW, but not for RAW |
|---|
| 288 | // |
|---|
| 289 | // assume that the mask files are missing unless we can find them. If VCALC data, |
|---|
| 290 | // then the Mask is missing by definition |
|---|
| 291 | maskMissing = 1 |
|---|
| 292 | |
|---|
| 293 | strswitch(type) // string switch |
|---|
| 294 | |
|---|
| 295 | // only one panel, simply pick that panel and move on out of the switch |
|---|
| 296 | case "FL": |
|---|
| 297 | case "FR": |
|---|
| 298 | case "FT": |
|---|
| 299 | case "FB": |
|---|
| 300 | case "ML": |
|---|
| 301 | case "MR": |
|---|
| 302 | case "MT": |
|---|
| 303 | case "MB": |
|---|
| 304 | case "B": |
|---|
| 305 | if(isVCALC) |
|---|
| 306 | WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr) |
|---|
| 307 | WAVE/Z iErr = $("iErr_"+detStr) // 2D errors -- may not exist, especially for simulation |
|---|
| 308 | else |
|---|
| 309 | Wave inten = V_getDetectorDataW(folderStr,detStr) |
|---|
| 310 | Wave iErr = V_getDetectorDataErrW(folderStr,detStr) |
|---|
| 311 | Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data") |
|---|
| 312 | if(WaveExists(mask) == 1) |
|---|
| 313 | maskMissing = 0 |
|---|
| 314 | endif |
|---|
| 315 | endif |
|---|
| 316 | NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr) |
|---|
| 317 | Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr) // 2D q-values |
|---|
| 318 | Wave phi = V_MakePhiMatrix(qTotal,folderStr,detStr,folderPath+instPath+detStr) |
|---|
| 319 | nSets = 1 |
|---|
| 320 | break |
|---|
| 321 | |
|---|
| 322 | case "FLR": |
|---|
| 323 | // detStr has multiple values now, so unfortuntely, I'm hard-wiring things... |
|---|
| 324 | // TODO |
|---|
| 325 | // -- see if I can un-hard-wire some of this below when more than one panel is combined |
|---|
| 326 | if(isVCALC) |
|---|
| 327 | WAVE inten = $(folderPath+instPath+"FL"+":det_"+"FL") |
|---|
| 328 | WAVE/Z iErr = $("iErr_"+"FL") // 2D errors -- may not exist, especially for simulation |
|---|
| 329 | WAVE inten2 = $(folderPath+instPath+"FR"+":det_"+"FR") |
|---|
| 330 | WAVE/Z iErr2 = $("iErr_"+"FR") // 2D errors -- may not exist, especially for simulation |
|---|
| 331 | else |
|---|
| 332 | Wave inten = V_getDetectorDataW(folderStr,"FL") |
|---|
| 333 | Wave iErr = V_getDetectorDataErrW(folderStr,"FL") |
|---|
| 334 | Wave inten2 = V_getDetectorDataW(folderStr,"FR") |
|---|
| 335 | Wave iErr2 = V_getDetectorDataErrW(folderStr,"FR") |
|---|
| 336 | Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FL"+":data") |
|---|
| 337 | Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FR"+":data") |
|---|
| 338 | if(WaveExists(mask) == 1 && WaveExists(mask2) == 1) |
|---|
| 339 | maskMissing = 0 |
|---|
| 340 | endif |
|---|
| 341 | endif |
|---|
| 342 | NVAR delQ = $(folderPath+instPath+"FL"+":gDelQ_FL") |
|---|
| 343 | |
|---|
| 344 | Wave qTotal = $(folderPath+instPath+"FL"+":qTot_"+"FL") // 2D q-values |
|---|
| 345 | Wave qTotal2 = $(folderPath+instPath+"FR"+":qTot_"+"FR") // 2D q-values |
|---|
| 346 | Wave phi = V_MakePhiMatrix(qTotal,folderStr,"FL",folderPath+instPath+"FL") |
|---|
| 347 | Wave phi2 = V_MakePhiMatrix(qTotal2,folderStr,"FR",folderPath+instPath+"FR") |
|---|
| 348 | |
|---|
| 349 | nSets = 2 |
|---|
| 350 | break |
|---|
| 351 | |
|---|
| 352 | case "FTB": |
|---|
| 353 | if(isVCALC) |
|---|
| 354 | WAVE inten = $(folderPath+instPath+"FT"+":det_"+"FT") |
|---|
| 355 | WAVE/Z iErr = $("iErr_"+"FT") // 2D errors -- may not exist, especially for simulation |
|---|
| 356 | WAVE inten2 = $(folderPath+instPath+"FB"+":det_"+"FB") |
|---|
| 357 | WAVE/Z iErr2 = $("iErr_"+"FB") // 2D errors -- may not exist, especially for simulation |
|---|
| 358 | else |
|---|
| 359 | Wave inten = V_getDetectorDataW(folderStr,"FT") |
|---|
| 360 | Wave iErr = V_getDetectorDataErrW(folderStr,"FT") |
|---|
| 361 | Wave inten2 = V_getDetectorDataW(folderStr,"FB") |
|---|
| 362 | Wave iErr2 = V_getDetectorDataErrW(folderStr,"FB") |
|---|
| 363 | Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FT"+":data") |
|---|
| 364 | Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FB"+":data") |
|---|
| 365 | if(WaveExists(mask) == 1 && WaveExists(mask2) == 1) |
|---|
| 366 | maskMissing = 0 |
|---|
| 367 | endif |
|---|
| 368 | endif |
|---|
| 369 | NVAR delQ = $(folderPath+instPath+"FT"+":gDelQ_FT") |
|---|
| 370 | |
|---|
| 371 | Wave qTotal = $(folderPath+instPath+"FT"+":qTot_"+"FT") // 2D q-values |
|---|
| 372 | Wave qTotal2 = $(folderPath+instPath+"FB"+":qTot_"+"FB") // 2D q-values |
|---|
| 373 | |
|---|
| 374 | Wave phi = V_MakePhiMatrix(qTotal,folderStr,"FT",folderPath+instPath+"FT") |
|---|
| 375 | Wave phi2 = V_MakePhiMatrix(qTotal2,folderStr,"FB",folderPath+instPath+"FB") |
|---|
| 376 | |
|---|
| 377 | nSets = 2 |
|---|
| 378 | break |
|---|
| 379 | |
|---|
| 380 | case "FLRTB": |
|---|
| 381 | if(isVCALC) |
|---|
| 382 | WAVE inten = $(folderPath+instPath+"FL"+":det_"+"FL") |
|---|
| 383 | WAVE/Z iErr = $("iErr_"+"FL") // 2D errors -- may not exist, especially for simulation |
|---|
| 384 | WAVE inten2 = $(folderPath+instPath+"FR"+":det_"+"FR") |
|---|
| 385 | WAVE/Z iErr2 = $("iErr_"+"FR") // 2D errors -- may not exist, especially for simulation |
|---|
| 386 | WAVE inten3 = $(folderPath+instPath+"FT"+":det_"+"FT") |
|---|
| 387 | WAVE/Z iErr3 = $("iErr_"+"FT") // 2D errors -- may not exist, especially for simulation |
|---|
| 388 | WAVE inten4 = $(folderPath+instPath+"FB"+":det_"+"FB") |
|---|
| 389 | WAVE/Z iErr4 = $("iErr_"+"FB") // 2D errors -- may not exist, especially for simulation |
|---|
| 390 | else |
|---|
| 391 | Wave inten = V_getDetectorDataW(folderStr,"FL") |
|---|
| 392 | Wave iErr = V_getDetectorDataErrW(folderStr,"FL") |
|---|
| 393 | Wave inten2 = V_getDetectorDataW(folderStr,"FR") |
|---|
| 394 | Wave iErr2 = V_getDetectorDataErrW(folderStr,"FR") |
|---|
| 395 | Wave inten3 = V_getDetectorDataW(folderStr,"FT") |
|---|
| 396 | Wave iErr3 = V_getDetectorDataErrW(folderStr,"FT") |
|---|
| 397 | Wave inten4 = V_getDetectorDataW(folderStr,"FB") |
|---|
| 398 | Wave iErr4 = V_getDetectorDataErrW(folderStr,"FB") |
|---|
| 399 | Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FL"+":data") |
|---|
| 400 | Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FR"+":data") |
|---|
| 401 | Wave/Z mask3 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FT"+":data") |
|---|
| 402 | Wave/Z mask4 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FB"+":data") |
|---|
| 403 | if(WaveExists(mask) == 1 && WaveExists(mask2) == 1 && WaveExists(mask3) == 1 && WaveExists(mask4) == 1) |
|---|
| 404 | maskMissing = 0 |
|---|
| 405 | endif |
|---|
| 406 | endif |
|---|
| 407 | NVAR delQ = $(folderPath+instPath+"FL"+":gDelQ_FL") |
|---|
| 408 | |
|---|
| 409 | Wave qTotal = $(folderPath+instPath+"FL"+":qTot_"+"FL") // 2D q-values |
|---|
| 410 | Wave qTotal2 = $(folderPath+instPath+"FR"+":qTot_"+"FR") // 2D q-values |
|---|
| 411 | Wave qTotal3 = $(folderPath+instPath+"FT"+":qTot_"+"FT") // 2D q-values |
|---|
| 412 | Wave qTotal4 = $(folderPath+instPath+"FB"+":qTot_"+"FB") // 2D q-values |
|---|
| 413 | |
|---|
| 414 | Wave phi = V_MakePhiMatrix(qTotal,folderStr,"FL",folderPath+instPath+"FL") |
|---|
| 415 | Wave phi2 = V_MakePhiMatrix(qTotal2,folderStr,"FR",folderPath+instPath+"FR") |
|---|
| 416 | Wave phi3 = V_MakePhiMatrix(qTotal3,folderStr,"FT",folderPath+instPath+"FT") |
|---|
| 417 | Wave phi4 = V_MakePhiMatrix(qTotal4,folderStr,"FB",folderPath+instPath+"FB") |
|---|
| 418 | |
|---|
| 419 | nSets = 4 |
|---|
| 420 | break |
|---|
| 421 | |
|---|
| 422 | case "MLR": |
|---|
| 423 | if(isVCALC) |
|---|
| 424 | WAVE inten = $(folderPath+instPath+"ML"+":det_"+"ML") |
|---|
| 425 | WAVE/Z iErr = $("iErr_"+"ML") // 2D errors -- may not exist, especially for simulation |
|---|
| 426 | WAVE inten2 = $(folderPath+instPath+"MR"+":det_"+"MR") |
|---|
| 427 | WAVE/Z iErr2 = $("iErr_"+"MR") // 2D errors -- may not exist, especially for simulation |
|---|
| 428 | else |
|---|
| 429 | Wave inten = V_getDetectorDataW(folderStr,"ML") |
|---|
| 430 | Wave iErr = V_getDetectorDataErrW(folderStr,"ML") |
|---|
| 431 | Wave inten2 = V_getDetectorDataW(folderStr,"MR") |
|---|
| 432 | Wave iErr2 = V_getDetectorDataErrW(folderStr,"MR") |
|---|
| 433 | Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"ML"+":data") |
|---|
| 434 | Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MR"+":data") |
|---|
| 435 | if(WaveExists(mask) == 1 && WaveExists(mask2) == 1) |
|---|
| 436 | maskMissing = 0 |
|---|
| 437 | endif |
|---|
| 438 | endif |
|---|
| 439 | NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML") |
|---|
| 440 | |
|---|
| 441 | Wave qTotal = $(folderPath+instPath+"ML"+":qTot_"+"ML") // 2D q-values |
|---|
| 442 | Wave qTotal2 = $(folderPath+instPath+"MR"+":qTot_"+"MR") // 2D q-values |
|---|
| 443 | |
|---|
| 444 | Wave phi = V_MakePhiMatrix(qTotal,folderStr,"ML",folderPath+instPath+"ML") |
|---|
| 445 | Wave phi2 = V_MakePhiMatrix(qTotal2,folderStr,"MR",folderPath+instPath+"MR") |
|---|
| 446 | |
|---|
| 447 | nSets = 2 |
|---|
| 448 | break |
|---|
| 449 | |
|---|
| 450 | case "MTB": |
|---|
| 451 | if(isVCALC) |
|---|
| 452 | WAVE inten = $(folderPath+instPath+"MT"+":det_"+"MT") |
|---|
| 453 | WAVE/Z iErr = $("iErr_"+"MT") // 2D errors -- may not exist, especially for simulation |
|---|
| 454 | WAVE inten2 = $(folderPath+instPath+"MB"+":det_"+"MB") |
|---|
| 455 | WAVE/Z iErr2 = $("iErr_"+"MB") // 2D errors -- may not exist, especially for simulation |
|---|
| 456 | else |
|---|
| 457 | Wave inten = V_getDetectorDataW(folderStr,"MT") |
|---|
| 458 | Wave iErr = V_getDetectorDataErrW(folderStr,"MT") |
|---|
| 459 | Wave inten2 = V_getDetectorDataW(folderStr,"MB") |
|---|
| 460 | Wave iErr2 = V_getDetectorDataErrW(folderStr,"MB") |
|---|
| 461 | Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MT"+":data") |
|---|
| 462 | Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MB"+":data") |
|---|
| 463 | if(WaveExists(mask) == 1 && WaveExists(mask2) == 1) |
|---|
| 464 | maskMissing = 0 |
|---|
| 465 | endif |
|---|
| 466 | endif |
|---|
| 467 | NVAR delQ = $(folderPath+instPath+"MT"+":gDelQ_MT") |
|---|
| 468 | |
|---|
| 469 | Wave qTotal = $(folderPath+instPath+"MT"+":qTot_"+"MT") // 2D q-values |
|---|
| 470 | Wave qTotal2 = $(folderPath+instPath+"MB"+":qTot_"+"MB") // 2D q-values |
|---|
| 471 | |
|---|
| 472 | Wave phi = V_MakePhiMatrix(qTotal,folderStr,"MT",folderPath+instPath+"MT") |
|---|
| 473 | Wave phi2 = V_MakePhiMatrix(qTotal2,folderStr,"MB",folderPath+instPath+"MB") |
|---|
| 474 | |
|---|
| 475 | nSets = 2 |
|---|
| 476 | break |
|---|
| 477 | |
|---|
| 478 | case "MLRTB": |
|---|
| 479 | if(isVCALC) |
|---|
| 480 | WAVE inten = $(folderPath+instPath+"ML"+":det_"+"ML") |
|---|
| 481 | WAVE/Z iErr = $("iErr_"+"ML") // 2D errors -- may not exist, especially for simulation |
|---|
| 482 | WAVE inten2 = $(folderPath+instPath+"MR"+":det_"+"MR") |
|---|
| 483 | WAVE/Z iErr2 = $("iErr_"+"MR") // 2D errors -- may not exist, especially for simulation |
|---|
| 484 | WAVE inten3 = $(folderPath+instPath+"MT"+":det_"+"MT") |
|---|
| 485 | WAVE/Z iErr3 = $("iErr_"+"MT") // 2D errors -- may not exist, especially for simulation |
|---|
| 486 | WAVE inten4 = $(folderPath+instPath+"MB"+":det_"+"MB") |
|---|
| 487 | WAVE/Z iErr4 = $("iErr_"+"MB") // 2D errors -- may not exist, especially for simulation |
|---|
| 488 | else |
|---|
| 489 | Wave inten = V_getDetectorDataW(folderStr,"ML") |
|---|
| 490 | Wave iErr = V_getDetectorDataErrW(folderStr,"ML") |
|---|
| 491 | Wave inten2 = V_getDetectorDataW(folderStr,"MR") |
|---|
| 492 | Wave iErr2 = V_getDetectorDataErrW(folderStr,"MR") |
|---|
| 493 | Wave inten3 = V_getDetectorDataW(folderStr,"MT") |
|---|
| 494 | Wave iErr3 = V_getDetectorDataErrW(folderStr,"MT") |
|---|
| 495 | Wave inten4 = V_getDetectorDataW(folderStr,"MB") |
|---|
| 496 | Wave iErr4 = V_getDetectorDataErrW(folderStr,"MB") |
|---|
| 497 | Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"ML"+":data") |
|---|
| 498 | Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MR"+":data") |
|---|
| 499 | Wave/Z mask3 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MT"+":data") |
|---|
| 500 | Wave/Z mask4 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MB"+":data") |
|---|
| 501 | if(WaveExists(mask) == 1 && WaveExists(mask2) == 1 && WaveExists(mask3) == 1 && WaveExists(mask4) == 1) |
|---|
| 502 | maskMissing = 0 |
|---|
| 503 | endif |
|---|
| 504 | endif |
|---|
| 505 | NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML") |
|---|
| 506 | |
|---|
| 507 | Wave qTotal = $(folderPath+instPath+"ML"+":qTot_"+"ML") // 2D q-values |
|---|
| 508 | Wave qTotal2 = $(folderPath+instPath+"MR"+":qTot_"+"MR") // 2D q-values |
|---|
| 509 | Wave qTotal3 = $(folderPath+instPath+"MT"+":qTot_"+"MT") // 2D q-values |
|---|
| 510 | Wave qTotal4 = $(folderPath+instPath+"MB"+":qTot_"+"MB") // 2D q-values |
|---|
| 511 | |
|---|
| 512 | Wave phi = V_MakePhiMatrix(qTotal,folderStr,"ML",folderPath+instPath+"ML") |
|---|
| 513 | Wave phi2 = V_MakePhiMatrix(qTotal2,folderStr,"MR",folderPath+instPath+"MR") |
|---|
| 514 | Wave phi3 = V_MakePhiMatrix(qTotal3,folderStr,"MT",folderPath+instPath+"MT") |
|---|
| 515 | Wave phi4 = V_MakePhiMatrix(qTotal4,folderStr,"MB",folderPath+instPath+"MB") |
|---|
| 516 | |
|---|
| 517 | nSets = 4 |
|---|
| 518 | break |
|---|
| 519 | |
|---|
| 520 | default: |
|---|
| 521 | nSets = 0 |
|---|
| 522 | Print "ERROR ---- type is not recognized " |
|---|
| 523 | endswitch |
|---|
| 524 | |
|---|
| 525 | // Print "delQ = ",delQ," for ",type |
|---|
| 526 | |
|---|
| 527 | if(nSets == 0) |
|---|
| 528 | SetDataFolder root: |
|---|
| 529 | return(0) |
|---|
| 530 | endif |
|---|
| 531 | |
|---|
| 532 | |
|---|
| 533 | // RAW data is currently read in and the 2D error wave is correctly generated |
|---|
| 534 | // 2D error is propagated through all reduction steps, but I have not |
|---|
| 535 | // verified that it is an exact duplication of the 1D error |
|---|
| 536 | // |
|---|
| 537 | // |
|---|
| 538 | // |
|---|
| 539 | // IF ther is no 2D error wave present for some reason, make a fake one |
|---|
| 540 | if(WaveExists(iErr)==0 && WaveExists(inten) != 0) |
|---|
| 541 | Duplicate/O inten,iErr |
|---|
| 542 | Wave iErr=iErr |
|---|
| 543 | // iErr = 1+sqrt(inten+0.75) // can't use this -- it applies to counts, not intensity (already a count rate...) |
|---|
| 544 | iErr = sqrt(inten+0.75) // TODO -- here I'm just using some fictional value |
|---|
| 545 | endif |
|---|
| 546 | if(WaveExists(iErr2)==0 && WaveExists(inten2) != 0) |
|---|
| 547 | Duplicate/O inten2,iErr2 |
|---|
| 548 | Wave iErr2=iErr2 |
|---|
| 549 | // iErr2 = 1+sqrt(inten2+0.75) // can't use this -- it applies to counts, not intensity (already a count rate...) |
|---|
| 550 | iErr2 = sqrt(inten2+0.75) // TODO -- here I'm just using some fictional value |
|---|
| 551 | endif |
|---|
| 552 | if(WaveExists(iErr3)==0 && WaveExists(inten3) != 0) |
|---|
| 553 | Duplicate/O inten3,iErr3 |
|---|
| 554 | Wave iErr3=iErr3 |
|---|
| 555 | // iErr3 = 1+sqrt(inten3+0.75) // can't use this -- it applies to counts, not intensity (already a count rate...) |
|---|
| 556 | iErr3 = sqrt(inten3+0.75) // TODO -- here I'm just using some fictional value |
|---|
| 557 | endif |
|---|
| 558 | if(WaveExists(iErr4)==0 && WaveExists(inten4) != 0) |
|---|
| 559 | Duplicate/O inten4,iErr4 |
|---|
| 560 | Wave iErr4=iErr4 |
|---|
| 561 | // iErr4 = 1+sqrt(inten4+0.75) // can't use this -- it applies to counts, not intensity (already a count rate...) |
|---|
| 562 | iErr4 = sqrt(inten4+0.75) // TODO -- here I'm just using some fictional value |
|---|
| 563 | endif |
|---|
| 564 | |
|---|
| 565 | // TODO -- nq will need to be larger, once the back detector is installed |
|---|
| 566 | // |
|---|
| 567 | if(cmpstr(type,"B") == 0) |
|---|
| 568 | nq = 8000 |
|---|
| 569 | else |
|---|
| 570 | nq=600 |
|---|
| 571 | endif |
|---|
| 572 | |
|---|
| 573 | //******TODO****** -- where to put the averaged data -- right now, folderStr is forced to "" |
|---|
| 574 | // SetDataFolder $("root:"+folderStr) //should already be here, but make sure... |
|---|
| 575 | Make/O/D/N=(nq) $(folderPath+":"+"iBin_qxqy"+"_"+type) |
|---|
| 576 | Make/O/D/N=(nq) $(folderPath+":"+"qBin_qxqy"+"_"+type) |
|---|
| 577 | Make/O/D/N=(nq) $(folderPath+":"+"nBin_qxqy"+"_"+type) |
|---|
| 578 | Make/O/D/N=(nq) $(folderPath+":"+"iBin2_qxqy"+"_"+type) |
|---|
| 579 | Make/O/D/N=(nq) $(folderPath+":"+"eBin_qxqy"+"_"+type) |
|---|
| 580 | Make/O/D/N=(nq) $(folderPath+":"+"eBin2D_qxqy"+"_"+type) |
|---|
| 581 | |
|---|
| 582 | Wave iBin_qxqy = $(folderPath+":"+"iBin_qxqy_"+type) |
|---|
| 583 | Wave qBin_qxqy = $(folderPath+":"+"qBin_qxqy"+"_"+type) |
|---|
| 584 | Wave nBin_qxqy = $(folderPath+":"+"nBin_qxqy"+"_"+type) |
|---|
| 585 | Wave iBin2_qxqy = $(folderPath+":"+"iBin2_qxqy"+"_"+type) |
|---|
| 586 | Wave eBin_qxqy = $(folderPath+":"+"eBin_qxqy"+"_"+type) |
|---|
| 587 | Wave eBin2D_qxqy = $(folderPath+":"+"eBin2D_qxqy"+"_"+type) |
|---|
| 588 | |
|---|
| 589 | |
|---|
| 590 | // delQ = abs(sqrt(qx[2]^2+qy[2]^2+qz[2]^2) - sqrt(qx[1]^2+qy[1]^2+qz[1]^2)) //use bins of 1 pixel width |
|---|
| 591 | // TODO: not sure if I want to set dQ in x or y direction... |
|---|
| 592 | // the short dimension is the 8mm tubes, use this direction as dQ? |
|---|
| 593 | // but don't use the corner of the detector, since dQ will be very different on T/B or L/R due to the location of [0,0] |
|---|
| 594 | // WRT the beam center. use qx or qy directly. Still not happy with this way... |
|---|
| 595 | |
|---|
| 596 | |
|---|
| 597 | qBin_qxqy[] = p*delQ |
|---|
| 598 | SetScale/P x,0,delQ,"",qBin_qxqy //allows easy binning |
|---|
| 599 | |
|---|
| 600 | iBin_qxqy = 0 |
|---|
| 601 | iBin2_qxqy = 0 |
|---|
| 602 | eBin_qxqy = 0 |
|---|
| 603 | eBin2D_qxqy = 0 |
|---|
| 604 | nBin_qxqy = 0 //number of intensities added to each bin |
|---|
| 605 | |
|---|
| 606 | // now there are situations of: |
|---|
| 607 | // 1 panel |
|---|
| 608 | // 2 panels |
|---|
| 609 | // 4 panels |
|---|
| 610 | // |
|---|
| 611 | // this needs to be a double loop now... |
|---|
| 612 | // TODO: |
|---|
| 613 | // -- the iErr (=2D) wave and accumulation of error is NOT CALCULATED CORRECTLY YET |
|---|
| 614 | // -- verify the 2D error propagation by reducing it to 1D error |
|---|
| 615 | // |
|---|
| 616 | // |
|---|
| 617 | // The 1D error does not use iErr, and IS CALCULATED CORRECTLY |
|---|
| 618 | // |
|---|
| 619 | // x- the solid angle per pixel will be present for WORK data other than RAW, but not for RAW |
|---|
| 620 | |
|---|
| 621 | // |
|---|
| 622 | // if any of the masks don't exist, display the error, and proceed with the averaging, using all data |
|---|
| 623 | if(maskMissing == 1) |
|---|
| 624 | Print "Mask file not found for at least one detector - so all data is used" |
|---|
| 625 | endif |
|---|
| 626 | |
|---|
| 627 | NVAR gIgnoreDetB = root:Packages:NIST:VSANS:Globals:gIgnoreDetB |
|---|
| 628 | if(gIgnoreDetB && cmpstr(type,"B") == 0) |
|---|
| 629 | maskMissing = 1 |
|---|
| 630 | Print "Mask skipped for B due to possible mismatch (Panel B ignored in preferences)" |
|---|
| 631 | endif |
|---|
| 632 | |
|---|
| 633 | Variable mask_val,phiVal,isIn |
|---|
| 634 | // use set 1 (no number) only |
|---|
| 635 | if(nSets >= 1) |
|---|
| 636 | xDim=DimSize(inten,0) |
|---|
| 637 | yDim=DimSize(inten,1) |
|---|
| 638 | |
|---|
| 639 | for(ii=0;ii<xDim;ii+=1) |
|---|
| 640 | for(jj=0;jj<yDim;jj+=1) |
|---|
| 641 | //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2) |
|---|
| 642 | qVal = qTotal[ii][jj] |
|---|
| 643 | binIndex = trunc(x2pnt(qBin_qxqy, qVal)) |
|---|
| 644 | val = inten[ii][jj] |
|---|
| 645 | |
|---|
| 646 | if(isVCALC || maskMissing) // mask_val == 0 == keep, mask_val == 1 = YES, mask out the point |
|---|
| 647 | mask_val = 0 |
|---|
| 648 | else |
|---|
| 649 | mask_val = mask[ii][jj] |
|---|
| 650 | endif |
|---|
| 651 | |
|---|
| 652 | phiVal = phi[ii][jj] |
|---|
| 653 | isIn = 0 // start with exclude, now see if it's a keeper |
|---|
| 654 | |
|---|
| 655 | // if within the right or left, flag to keep the pixel |
|---|
| 656 | if(cmpstr(side,"right")==0) |
|---|
| 657 | //right, when 0->pi/2 |
|---|
| 658 | if(V_CloseEnough(phiVal,phi_rad,dphi_rad)) |
|---|
| 659 | isIn = 1 |
|---|
| 660 | endif |
|---|
| 661 | // condition here to get the 3pi/2 -> 2pi region |
|---|
| 662 | if(V_CloseEnough(phiVal,phi_rad+2*pi,dphi_rad)) |
|---|
| 663 | isIn = 1 |
|---|
| 664 | endif |
|---|
| 665 | endif |
|---|
| 666 | |
|---|
| 667 | if(cmpstr(side,"left")==0) |
|---|
| 668 | if(V_CloseEnough(phiVal,phi_rad+pi,dphi_rad)) |
|---|
| 669 | isIn = 1 |
|---|
| 670 | endif |
|---|
| 671 | endif |
|---|
| 672 | |
|---|
| 673 | // both sides, duplicates the conditions above |
|---|
| 674 | if(cmpstr(side,"both")==0) |
|---|
| 675 | //right, when 0->pi/2 |
|---|
| 676 | if(V_CloseEnough(phiVal,phi_rad,dphi_rad)) |
|---|
| 677 | isIn = 1 |
|---|
| 678 | endif |
|---|
| 679 | // right, when 3pi/2 -> 2pi |
|---|
| 680 | if(V_CloseEnough(phiVal,phi_rad+2*pi,dphi_rad)) |
|---|
| 681 | isIn = 1 |
|---|
| 682 | endif |
|---|
| 683 | |
|---|
| 684 | //left |
|---|
| 685 | if(V_CloseEnough(phiVal,phi_rad+pi,dphi_rad)) |
|---|
| 686 | isIn = 1 |
|---|
| 687 | endif |
|---|
| 688 | |
|---|
| 689 | endif //end the check of phiVal within sector and side |
|---|
| 690 | |
|---|
| 691 | |
|---|
| 692 | if (numType(val)==0 && mask_val == 0 && isIn > 0) //count only the good points, in the sector, and ignore Nan or Inf |
|---|
| 693 | iBin_qxqy[binIndex] += val |
|---|
| 694 | iBin2_qxqy[binIndex] += val*val |
|---|
| 695 | eBin2D_qxqy[binIndex] += iErr[ii][jj]*iErr[ii][jj] |
|---|
| 696 | nBin_qxqy[binIndex] += 1 |
|---|
| 697 | endif |
|---|
| 698 | |
|---|
| 699 | |
|---|
| 700 | endfor |
|---|
| 701 | endfor |
|---|
| 702 | |
|---|
| 703 | endif |
|---|
| 704 | |
|---|
| 705 | // add in set 2 (set 1 already done) |
|---|
| 706 | if(nSets >= 2) |
|---|
| 707 | xDim=DimSize(inten2,0) |
|---|
| 708 | yDim=DimSize(inten2,1) |
|---|
| 709 | |
|---|
| 710 | for(ii=0;ii<xDim;ii+=1) |
|---|
| 711 | for(jj=0;jj<yDim;jj+=1) |
|---|
| 712 | //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2) |
|---|
| 713 | qVal = qTotal2[ii][jj] |
|---|
| 714 | binIndex = trunc(x2pnt(qBin_qxqy, qVal)) |
|---|
| 715 | val = inten2[ii][jj] |
|---|
| 716 | |
|---|
| 717 | if(isVCALC || maskMissing) |
|---|
| 718 | mask_val = 0 |
|---|
| 719 | else |
|---|
| 720 | mask_val = mask2[ii][jj] |
|---|
| 721 | endif |
|---|
| 722 | |
|---|
| 723 | phiVal = phi2[ii][jj] |
|---|
| 724 | isIn = 0 // start with exclude, now see if it's a keeper |
|---|
| 725 | |
|---|
| 726 | // if within the right or left, flag to keep the pixel |
|---|
| 727 | if(cmpstr(side,"right")==0) |
|---|
| 728 | //right, when 0->pi/2 |
|---|
| 729 | if(V_CloseEnough(phiVal,phi_rad,dphi_rad)) |
|---|
| 730 | isIn = 1 |
|---|
| 731 | endif |
|---|
| 732 | // condition here to get the 3pi/2 -> 2pi region |
|---|
| 733 | if(V_CloseEnough(phiVal,phi_rad+2*pi,dphi_rad)) |
|---|
| 734 | isIn = 1 |
|---|
| 735 | endif |
|---|
| 736 | endif |
|---|
| 737 | |
|---|
| 738 | if(cmpstr(side,"left")==0) |
|---|
| 739 | if(V_CloseEnough(phiVal,phi_rad+pi,dphi_rad)) |
|---|
| 740 | isIn = 1 |
|---|
| 741 | endif |
|---|
| 742 | endif |
|---|
| 743 | |
|---|
| 744 | // both sides, duplicates the conditions above |
|---|
| 745 | if(cmpstr(side,"both")==0) |
|---|
| 746 | //right, when 0->pi/2 |
|---|
| 747 | if(V_CloseEnough(phiVal,phi_rad,dphi_rad)) |
|---|
| 748 | isIn = 1 |
|---|
| 749 | endif |
|---|
| 750 | // right, when 3pi/2 -> 2pi |
|---|
| 751 | if(V_CloseEnough(phiVal,phi_rad+2*pi,dphi_rad)) |
|---|
| 752 | isIn = 1 |
|---|
| 753 | endif |
|---|
| 754 | |
|---|
| 755 | //left |
|---|
| 756 | if(V_CloseEnough(phiVal,phi_rad+pi,dphi_rad)) |
|---|
| 757 | isIn = 1 |
|---|
| 758 | endif |
|---|
| 759 | |
|---|
| 760 | endif //end the check of phiVal within sector and side |
|---|
| 761 | |
|---|
| 762 | |
|---|
| 763 | if (numType(val)==0 && mask_val == 0 && isIn > 0) //count only the good points, ignore Nan or Inf |
|---|
| 764 | iBin_qxqy[binIndex] += val |
|---|
| 765 | iBin2_qxqy[binIndex] += val*val |
|---|
| 766 | eBin2D_qxqy[binIndex] += iErr2[ii][jj]*iErr2[ii][jj] |
|---|
| 767 | nBin_qxqy[binIndex] += 1 |
|---|
| 768 | endif |
|---|
| 769 | endfor |
|---|
| 770 | endfor |
|---|
| 771 | |
|---|
| 772 | endif |
|---|
| 773 | |
|---|
| 774 | // add in set 3 and 4 (set 1 and 2 already done) |
|---|
| 775 | if(nSets == 4) |
|---|
| 776 | xDim=DimSize(inten3,0) |
|---|
| 777 | yDim=DimSize(inten3,1) |
|---|
| 778 | |
|---|
| 779 | for(ii=0;ii<xDim;ii+=1) |
|---|
| 780 | for(jj=0;jj<yDim;jj+=1) |
|---|
| 781 | //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2) |
|---|
| 782 | qVal = qTotal3[ii][jj] |
|---|
| 783 | binIndex = trunc(x2pnt(qBin_qxqy, qVal)) |
|---|
| 784 | val = inten3[ii][jj] |
|---|
| 785 | |
|---|
| 786 | if(isVCALC || maskMissing) |
|---|
| 787 | mask_val = 0 |
|---|
| 788 | else |
|---|
| 789 | mask_val = mask3[ii][jj] |
|---|
| 790 | endif |
|---|
| 791 | |
|---|
| 792 | phiVal = phi3[ii][jj] |
|---|
| 793 | isIn = 0 // start with exclude, now see if it's a keeper |
|---|
| 794 | |
|---|
| 795 | // if within the right or left, flag to keep the pixel |
|---|
| 796 | if(cmpstr(side,"right")==0) |
|---|
| 797 | //right, when 0->pi/2 |
|---|
| 798 | if(V_CloseEnough(phiVal,phi_rad,dphi_rad)) |
|---|
| 799 | isIn = 1 |
|---|
| 800 | endif |
|---|
| 801 | // condition here to get the 3pi/2 -> 2pi region |
|---|
| 802 | if(V_CloseEnough(phiVal,phi_rad+2*pi,dphi_rad)) |
|---|
| 803 | isIn = 1 |
|---|
| 804 | endif |
|---|
| 805 | endif |
|---|
| 806 | |
|---|
| 807 | if(cmpstr(side,"left")==0) |
|---|
| 808 | if(V_CloseEnough(phiVal,phi_rad+pi,dphi_rad)) |
|---|
| 809 | isIn = 1 |
|---|
| 810 | endif |
|---|
| 811 | endif |
|---|
| 812 | |
|---|
| 813 | // both sides, duplicates the conditions above |
|---|
| 814 | if(cmpstr(side,"both")==0) |
|---|
| 815 | //right, when 0->pi/2 |
|---|
| 816 | if(V_CloseEnough(phiVal,phi_rad,dphi_rad)) |
|---|
| 817 | isIn = 1 |
|---|
| 818 | endif |
|---|
| 819 | // right, when 3pi/2 -> 2pi |
|---|
| 820 | if(V_CloseEnough(phiVal,phi_rad+2*pi,dphi_rad)) |
|---|
| 821 | isIn = 1 |
|---|
| 822 | endif |
|---|
| 823 | |
|---|
| 824 | //left |
|---|
| 825 | if(V_CloseEnough(phiVal,phi_rad+pi,dphi_rad)) |
|---|
| 826 | isIn = 1 |
|---|
| 827 | endif |
|---|
| 828 | |
|---|
| 829 | endif //end the check of phiVal within sector and side |
|---|
| 830 | |
|---|
| 831 | if (numType(val)==0 && mask_val == 0 && isIn > 0) //count only the good points, ignore Nan or Inf |
|---|
| 832 | iBin_qxqy[binIndex] += val |
|---|
| 833 | iBin2_qxqy[binIndex] += val*val |
|---|
| 834 | eBin2D_qxqy[binIndex] += iErr3[ii][jj]*iErr3[ii][jj] |
|---|
| 835 | nBin_qxqy[binIndex] += 1 |
|---|
| 836 | endif |
|---|
| 837 | endfor |
|---|
| 838 | endfor |
|---|
| 839 | |
|---|
| 840 | |
|---|
| 841 | xDim=DimSize(inten4,0) |
|---|
| 842 | yDim=DimSize(inten4,1) |
|---|
| 843 | |
|---|
| 844 | for(ii=0;ii<xDim;ii+=1) |
|---|
| 845 | for(jj=0;jj<yDim;jj+=1) |
|---|
| 846 | //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2) |
|---|
| 847 | qVal = qTotal4[ii][jj] |
|---|
| 848 | binIndex = trunc(x2pnt(qBin_qxqy, qVal)) |
|---|
| 849 | val = inten4[ii][jj] |
|---|
| 850 | |
|---|
| 851 | if(isVCALC || maskMissing) |
|---|
| 852 | mask_val = 0 |
|---|
| 853 | else |
|---|
| 854 | mask_val = mask4[ii][jj] |
|---|
| 855 | endif |
|---|
| 856 | |
|---|
| 857 | phiVal = phi4[ii][jj] |
|---|
| 858 | isIn = 0 // start with exclude, now see if it's a keeper |
|---|
| 859 | |
|---|
| 860 | // if within the right or left, flag to keep the pixel |
|---|
| 861 | if(cmpstr(side,"right")==0) |
|---|
| 862 | //right, when 0->pi/2 |
|---|
| 863 | if(V_CloseEnough(phiVal,phi_rad,dphi_rad)) |
|---|
| 864 | isIn = 1 |
|---|
| 865 | endif |
|---|
| 866 | // condition here to get the 3pi/2 -> 2pi region |
|---|
| 867 | if(V_CloseEnough(phiVal,phi_rad+2*pi,dphi_rad)) |
|---|
| 868 | isIn = 1 |
|---|
| 869 | endif |
|---|
| 870 | endif |
|---|
| 871 | |
|---|
| 872 | if(cmpstr(side,"left")==0) |
|---|
| 873 | if(V_CloseEnough(phiVal,phi_rad+pi,dphi_rad)) |
|---|
| 874 | isIn = 1 |
|---|
| 875 | endif |
|---|
| 876 | endif |
|---|
| 877 | |
|---|
| 878 | // both sides, duplicates the conditions above |
|---|
| 879 | if(cmpstr(side,"both")==0) |
|---|
| 880 | //right, when 0->pi/2 |
|---|
| 881 | if(V_CloseEnough(phiVal,phi_rad,dphi_rad)) |
|---|
| 882 | isIn = 1 |
|---|
| 883 | endif |
|---|
| 884 | // right, when 3pi/2 -> 2pi |
|---|
| 885 | if(V_CloseEnough(phiVal,phi_rad+2*pi,dphi_rad)) |
|---|
| 886 | isIn = 1 |
|---|
| 887 | endif |
|---|
| 888 | |
|---|
| 889 | //left |
|---|
| 890 | if(V_CloseEnough(phiVal,phi_rad+pi,dphi_rad)) |
|---|
| 891 | isIn = 1 |
|---|
| 892 | endif |
|---|
| 893 | |
|---|
| 894 | endif //end the check of phiVal within sector and side |
|---|
| 895 | |
|---|
| 896 | if (numType(val)==0 && mask_val == 0 && isIn > 0) //count only the good points, ignore Nan or Inf |
|---|
| 897 | iBin_qxqy[binIndex] += val |
|---|
| 898 | iBin2_qxqy[binIndex] += val*val |
|---|
| 899 | eBin2D_qxqy[binIndex] += iErr4[ii][jj]*iErr4[ii][jj] |
|---|
| 900 | nBin_qxqy[binIndex] += 1 |
|---|
| 901 | endif |
|---|
| 902 | endfor |
|---|
| 903 | endfor |
|---|
| 904 | |
|---|
| 905 | endif |
|---|
| 906 | |
|---|
| 907 | |
|---|
| 908 | // after looping through all of the data on the panels, calculate errors on I(q), |
|---|
| 909 | // just like in CircSectAve.ipf |
|---|
| 910 | // TODO: |
|---|
| 911 | // -- 2D Errors were (maybe) properly acculumated through reduction, so this loop of calculations is NOT VERIFIED (yet) |
|---|
| 912 | // x- the error on the 1D intensity, is correctly calculated as the standard error of the mean. |
|---|
| 913 | for(ii=0;ii<nq;ii+=1) |
|---|
| 914 | if(nBin_qxqy[ii] == 0) |
|---|
| 915 | //no pixels in annuli, data unknown |
|---|
| 916 | iBin_qxqy[ii] = 0 |
|---|
| 917 | eBin_qxqy[ii] = 1 |
|---|
| 918 | eBin2D_qxqy[ii] = NaN |
|---|
| 919 | else |
|---|
| 920 | if(nBin_qxqy[ii] <= 1) |
|---|
| 921 | //need more than one pixel to determine error |
|---|
| 922 | iBin_qxqy[ii] /= nBin_qxqy[ii] |
|---|
| 923 | eBin_qxqy[ii] = 1 |
|---|
| 924 | eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2 |
|---|
| 925 | else |
|---|
| 926 | //assume that the intensity in each pixel in annuli is normally distributed about mean... |
|---|
| 927 | // -- this is correctly calculating the error as the standard error of the mean, as |
|---|
| 928 | // was always done for SANS as well. |
|---|
| 929 | iBin_qxqy[ii] /= nBin_qxqy[ii] |
|---|
| 930 | avesq = iBin_qxqy[ii]^2 |
|---|
| 931 | aveisq = iBin2_qxqy[ii]/nBin_qxqy[ii] |
|---|
| 932 | var = aveisq-avesq |
|---|
| 933 | if(var<=0) |
|---|
| 934 | eBin_qxqy[ii] = 1e-6 |
|---|
| 935 | else |
|---|
| 936 | eBin_qxqy[ii] = sqrt(var/(nBin_qxqy[ii] - 1)) |
|---|
| 937 | endif |
|---|
| 938 | // and calculate as it is propagated pixel-by-pixel |
|---|
| 939 | eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2 |
|---|
| 940 | endif |
|---|
| 941 | endif |
|---|
| 942 | endfor |
|---|
| 943 | |
|---|
| 944 | eBin2D_qxqy = sqrt(eBin2D_qxqy) // as equation (3) of John's memo |
|---|
| 945 | |
|---|
| 946 | // find the last non-zero point, working backwards |
|---|
| 947 | val=nq |
|---|
| 948 | do |
|---|
| 949 | val -= 1 |
|---|
| 950 | while((nBin_qxqy[val] == 0) && val > 0) |
|---|
| 951 | |
|---|
| 952 | // print val, nBin_qxqy[val] |
|---|
| 953 | DeletePoints val, nq-val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy |
|---|
| 954 | |
|---|
| 955 | if(val == 0) |
|---|
| 956 | // all the points were deleted, make dummy waves for resolution |
|---|
| 957 | Make/O/D/N=0 $(folderPath+":"+"sigmaQ_qxqy"+"_"+type) |
|---|
| 958 | Make/O/D/N=0 $(folderPath+":"+"qBar_qxqy"+"_"+type) |
|---|
| 959 | Make/O/D/N=0 $(folderPath+":"+"fSubS_qxqy"+"_"+type) |
|---|
| 960 | return(0) |
|---|
| 961 | endif |
|---|
| 962 | |
|---|
| 963 | |
|---|
| 964 | // since the beam center is not always on the detector, many of the low Q bins will have zero pixels |
|---|
| 965 | // find the first non-zero point, working forwards |
|---|
| 966 | val = -1 |
|---|
| 967 | do |
|---|
| 968 | val += 1 |
|---|
| 969 | while(nBin_qxqy[val] == 0) |
|---|
| 970 | DeletePoints 0, val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy |
|---|
| 971 | |
|---|
| 972 | // ?? there still may be a point in the q-range that gets zero pixel contribution - so search this out and get rid of it |
|---|
| 973 | val = numpnts(nBin_qxqy)-1 |
|---|
| 974 | do |
|---|
| 975 | if(nBin_qxqy[val] == 0) |
|---|
| 976 | DeletePoints val, 1, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy |
|---|
| 977 | endif |
|---|
| 978 | val -= 1 |
|---|
| 979 | while(val>0) |
|---|
| 980 | |
|---|
| 981 | // utility function to remove NaN values from the waves |
|---|
| 982 | V_RemoveNaNsQIS(qBin_qxqy, iBin_qxqy, eBin_qxqy) |
|---|
| 983 | |
|---|
| 984 | |
|---|
| 985 | // TODO: |
|---|
| 986 | // -- This is where I calculate the resolution in SANS (see CircSectAve) |
|---|
| 987 | // -- use the isVCALC flag to exclude VCALC from the resolution calculation if necessary |
|---|
| 988 | // -- from the top of the function, folderStr = work folder, type = "FLRTB" or other type of averaging |
|---|
| 989 | // |
|---|
| 990 | nq = numpnts(qBin_qxqy) |
|---|
| 991 | Make/O/D/N=(nq) $(folderPath+":"+"sigmaQ_qxqy"+"_"+type) |
|---|
| 992 | Make/O/D/N=(nq) $(folderPath+":"+"qBar_qxqy"+"_"+type) |
|---|
| 993 | Make/O/D/N=(nq) $(folderPath+":"+"fSubS_qxqy"+"_"+type) |
|---|
| 994 | Wave sigmaq = $(folderPath+":"+"sigmaQ_qxqy_"+type) |
|---|
| 995 | Wave qbar = $(folderPath+":"+"qBar_qxqy_"+type) |
|---|
| 996 | Wave fsubs = $(folderPath+":"+"fSubS_qxqy_"+type) |
|---|
| 997 | Variable ret1,ret2,ret3 |
|---|
| 998 | |
|---|
| 999 | // all of the different collimation conditions are handled within the V_getResolution function |
|---|
| 1000 | // which is responsible for switching based on the different collimation types (white beam, slit, Xtal, etc) |
|---|
| 1001 | // to calculate the correct resolution, or fill the waves with the correct "flags" |
|---|
| 1002 | // |
|---|
| 1003 | |
|---|
| 1004 | // For white beam data, the wavelength distribution can't be represented as a gaussian, but all of the other |
|---|
| 1005 | // geometric corrections still apply. Passing zero for the lambdaWidth will return the geometry contribution, |
|---|
| 1006 | // as long as the wavelength can be handled separately. It appears to be correct to do as a double integral, |
|---|
| 1007 | // with the inner(lambda) calculated first, then the outer(geometry). |
|---|
| 1008 | // |
|---|
| 1009 | |
|---|
| 1010 | // possible values are: |
|---|
| 1011 | // |
|---|
| 1012 | // pinhole |
|---|
| 1013 | // pinhole_whiteBeam |
|---|
| 1014 | // convergingPinholes |
|---|
| 1015 | // |
|---|
| 1016 | // *slit data should be reduced using the slit routine, not here, proceed but warn |
|---|
| 1017 | // narrowSlit |
|---|
| 1018 | // narrowSlit_whiteBeam |
|---|
| 1019 | // |
|---|
| 1020 | ii=0 |
|---|
| 1021 | do |
|---|
| 1022 | V_getResolution(qBin_qxqy[ii],folderStr,type,collimationStr,ret1,ret2,ret3) |
|---|
| 1023 | sigmaq[ii] = ret1 |
|---|
| 1024 | qbar[ii] = ret2 |
|---|
| 1025 | fsubs[ii] = ret3 |
|---|
| 1026 | ii+=1 |
|---|
| 1027 | while(ii<nq) |
|---|
| 1028 | |
|---|
| 1029 | |
|---|
| 1030 | SetDataFolder root: |
|---|
| 1031 | |
|---|
| 1032 | return(0) |
|---|
| 1033 | End |
|---|
| 1034 | |
|---|
