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

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

corrected compiling issue with version 7.87

Fixed logic in defining sector mask
Added mask for Annular average

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