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

Last change on this file since 1121 was 1121, checked in by srkline, 4 years ago

adding procedures for calculating sector averages of the data

still need to add overlay of sector angles.

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