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

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

updating the IgorVersion? pragma to v7.0 for all files to be consistent.

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