source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_Sector_Average.ipf @ 1217

Last change on this file since 1217 was 1217, checked in by srkline, 3 years ago

fixed minor bugs in attenuator table definition (didn't affect calculations), VCALC display, and sector averaging of the back detector on VSANS

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