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