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

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

changes to read in new USANS Raw data structure, based on the file creation date.

Some additional changes to sector averaging and viewing the "avg" masks on the detector panels. still not quite complete.

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