source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_IQ_Annular.ipf @ 1056

Last change on this file since 1056 was 1056, checked in by srkline, 5 years ago

added Annular averaging routine. still need writer, and link with protocols.

replaced tic/toc with v_tic/v_toc to avoid missing functions

File size: 28.4 KB
Line 
1#pragma TextEncoding = "MacRoman"
2#pragma rtGlobals=3             // Use modern global access method and strict wave access.
3
4
5// Procedures to do an annular binning of the data
6//
7// As for SANS, needs a Q-center and Q-delta to define the annular ring,
8// and the number of bins to divide the 360 circle
9//
10//
11//
12// TODO
13// x- add error bars to the plot of phi
14// -- data writer to export annular data (3-column)
15// -- loader to re-read annular data (will the normal loader work?
16// -- integrate this with the protocol
17// -- integrate this with a more general "average panel"
18// -- draw the q-center and width on the image (as a contour?)
19
20Proc Annular_Binning(folderStr,detGroup,qCtr_Ann,qWidth)
21        String folderStr="SAM",detGroup="F"
22        Variable qCtr_Ann=0.1,qWidth=0.01       // +/- in A^-1
23       
24        V_QBinAllPanels_Annular(folderStr,detGroup,qCtr_Ann,qWidth)
25       
26        Phi_Graph(folderStr)
27       
28End
29
30
31// TODO -- binType == 4 (slit mode) should never end up here, as it makes no sense
32//
33// -- really, the onle binning that makes any sense is "one", treating each panel individually,
34// so I may scrap the parameter, or ignore it. so don't count on it in the future.
35//
36Function V_QBinAllPanels_Annular(folderStr,detGroup,qCtr_Ann,qWidth)
37        String folderStr,detGroup
38        Variable qCtr_Ann,qWidth        // +/- in A^-1
39
40        Variable ii,delQ
41        String detStr
42       
43        if(cmpstr(detGroup,"F") == 0)
44                detStr = "FLRTB"
45        else
46                detStr = "MLRTB"
47        endif
48       
49// right now, use all of the detectors. There is a lot of waste in this and it could be
50// done a lot faster, but...
51
52
53// TODO         
54                // detStr = "FLRTB" or "MLRTB", depending which panel the q-ring is centered on/
55                // for now, no crossing of the rings onto different panels
56               
57        V_fDoAnnularBin_QxQy2D(folderStr,detStr,qCtr_Ann,qWidth)
58
59
60       
61
62        return(0)
63End
64
65Window Phi_graph(folderStr) : Graph
66        String folderStr
67       
68        PauseUpdate; Silent 1           // building window...
69        String fldrSav0= GetDataFolder(1)
70        SetDataFolder $("root:Packages:NIST:VSANS:"+folderStr)
71        Display /W=(35,45,572,419) iBin_qxqy_FLRTB vs phiBin_qxqy_FLRTB
72        ModifyGraph mode=4
73        ModifyGraph marker=19
74        ErrorBars iBin_qxqy_FLRTB Y,wave=(eBin_qxqy_FLRTB,eBin_qxqy_FLRTB)
75
76       
77        SetDataFolder fldrSav0
78EndMacro
79
80
81//////////
82//
83//
84//
85// TODO
86// -- "iErr" is not always defined correctly since it doesn't really apply here for data that is not 2D simulation
87//
88// ** Currently, type is being passed in as "" and ignored (looping through all of the detector panels
89// to potentially add to the annular bins)
90//
91// folderStr = WORK folder, type = the binning type (may include multiple detectors)
92//
93//      Variable qCtr_Ann = 0.1
94//      Variable qWidth = 0.02          // +/- in A^-1
95//
96Function V_fDoAnnularBin_QxQy2D(folderStr,type,qCtr_Ann,qWidth)
97        String folderStr,type
98        Variable qCtr_Ann,qWidth
99       
100        Variable nSets = 0
101        Variable xDim,yDim
102        Variable ii,jj
103        Variable qVal,nq,var,avesq,aveisq
104        Variable binIndex,val,isVCALC=0,maskMissing
105
106        String folderPath = "root:Packages:NIST:VSANS:"+folderStr
107        String instPath = ":entry:instrument:detector_"
108        String detStr
109               
110        if(cmpstr(folderStr,"VCALC") == 0)
111                isVCALC = 1
112        endif
113       
114// now switch on the type to determine which waves to declare and create
115// since there may be more than one panel to step through. There may be two, there may be four
116//
117
118// TODO:
119// -- Solid_Angle -- waves will be present for WORK data other than RAW, but not for RAW
120//
121// assume that the mask files are missing unless we can find them. If VCALC data,
122//  then the Mask is missing by definition
123        maskMissing = 1
124
125        strswitch(type) // string switch
126//              case "FL":              // execute if case matches expression
127//              case "FR":
128//                      detStr = type
129//                      if(isVCALC)
130//                              WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
131//                              WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation
132//                      else
133//                              Wave inten = V_getDetectorDataW(folderStr,detStr)
134//                              Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
135//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
136//                              if(WaveExists(mask) == 1)
137//                                      maskMissing = 0
138//                              endif
139//                             
140//                      endif
141//                      NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
142//                      Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
143//                      nSets = 1
144//                      break   
145//                                                             
146//              case "FT":             
147//              case "FB":
148//                      detStr = type
149//                      if(isVCALC)
150//                              WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
151//                              WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
152//                      else
153//                              Wave inten = V_getDetectorDataW(folderStr,detStr)
154//                              Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
155//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
156//                              if(WaveExists(mask) == 1)
157//                                      maskMissing = 0
158//                              endif
159//                      endif
160//                      NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
161//                      Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
162//                      nSets = 1
163//                      break
164//                     
165//              case "ML":             
166//              case "MR":
167//                      detStr = type
168//                      if(isVCALC)
169//                              WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
170//                              WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
171//                      else
172//                              Wave inten = V_getDetectorDataW(folderStr,detStr)
173//                              Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
174//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
175//                              if(WaveExists(mask) == 1)
176//                                      maskMissing = 0
177//                              endif
178//                      endif   
179//                      //TODO:
180//                      // -- decide on the proper deltaQ for binning. either nominal value for LR, or one
181//                      //    determined specifically for that panel (currently using one tube width as deltaQ)
182//                      // -- this is repeated multiple times in this switch
183//                      NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
184////                    NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML")
185//                      Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
186//                      nSets = 1
187//                      break   
188//                                     
189//              case "MT":             
190//              case "MB":
191//                      detStr = type
192//                      if(isVCALC)
193//                              WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
194//                              WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
195//                      else
196//                              Wave inten = V_getDetectorDataW(folderStr,detStr)
197//                              Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
198//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
199//                              if(WaveExists(mask) == 1)
200//                                      maskMissing = 0
201//                              endif
202//                      endif   
203//                      NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
204//                      Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
205//                      nSets = 1
206//                      break   
207//                                     
208//              case "B":       
209//                      detStr = type
210//                      if(isVCALC)
211//                              WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
212//                              WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
213//                      else
214//                              Wave inten = V_getDetectorDataW(folderStr,detStr)
215//                              Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
216//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
217//                              if(WaveExists(mask) == 1)
218//                                      maskMissing = 0
219//                              endif
220//                      endif   
221//                      NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_B")
222//                      Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values 
223//                      nSets = 1
224//                      break   
225//                     
226//              case "FLR":
227//              // detStr has multiple values now, so unfortuntely, I'm hard-wiring things...
228//              // TODO
229//              // -- see if I can un-hard-wire some of this below when more than one panel is combined
230//                      if(isVCALC)
231//                              WAVE inten = $(folderPath+instPath+"FL"+":det_"+"FL")
232//                              WAVE/Z iErr = $("iErr_"+"FL")                   // 2D errors -- may not exist, especially for simulation               
233//                              WAVE inten2 = $(folderPath+instPath+"FR"+":det_"+"FR")
234//                              WAVE/Z iErr2 = $("iErr_"+"FR")                  // 2D errors -- may not exist, especially for simulation       
235//                      else
236//                              Wave inten = V_getDetectorDataW(folderStr,"FL")
237//                              Wave iErr = V_getDetectorDataErrW(folderStr,"FL")
238//                              Wave inten2 = V_getDetectorDataW(folderStr,"FR")
239//                              Wave iErr2 = V_getDetectorDataErrW(folderStr,"FR")
240//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FL"+":data")
241//                              Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FR"+":data")
242//                              if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
243//                                      maskMissing = 0
244//                              endif
245//                      endif   
246//                      NVAR delQ = $(folderPath+instPath+"FL"+":gDelQ_FL")
247//                     
248//                      Wave qTotal = $(folderPath+instPath+"FL"+":qTot_"+"FL")                 // 2D q-values 
249//                      Wave qTotal2 = $(folderPath+instPath+"FR"+":qTot_"+"FR")                        // 2D q-values 
250//             
251//                      nSets = 2
252//                      break                   
253//             
254//              case "FTB":
255//                      if(isVCALC)
256//                              WAVE inten = $(folderPath+instPath+"FT"+":det_"+"FT")
257//                              WAVE/Z iErr = $("iErr_"+"FT")                   // 2D errors -- may not exist, especially for simulation               
258//                              WAVE inten2 = $(folderPath+instPath+"FB"+":det_"+"FB")
259//                              WAVE/Z iErr2 = $("iErr_"+"FB")                  // 2D errors -- may not exist, especially for simulation       
260//                      else
261//                              Wave inten = V_getDetectorDataW(folderStr,"FT")
262//                              Wave iErr = V_getDetectorDataErrW(folderStr,"FT")
263//                              Wave inten2 = V_getDetectorDataW(folderStr,"FB")
264//                              Wave iErr2 = V_getDetectorDataErrW(folderStr,"FB")
265//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FT"+":data")
266//                              Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FB"+":data")
267//                              if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
268//                                      maskMissing = 0
269//                              endif
270//                      endif   
271//                      NVAR delQ = $(folderPath+instPath+"FT"+":gDelQ_FT")
272//                     
273//                      Wave qTotal = $(folderPath+instPath+"FT"+":qTot_"+"FT")                 // 2D q-values 
274//                      Wave qTotal2 = $(folderPath+instPath+"FB"+":qTot_"+"FB")                        // 2D q-values 
275//     
276//                      nSets = 2
277//                      break           
278               
279                case "FLRTB":
280                        if(isVCALC)
281                                WAVE inten = $(folderPath+instPath+"FL"+":det_"+"FL")
282                                WAVE/Z iErr = $("iErr_"+"FL")                   // 2D errors -- may not exist, especially for simulation               
283                                WAVE inten2 = $(folderPath+instPath+"FR"+":det_"+"FR")
284                                WAVE/Z iErr2 = $("iErr_"+"FR")                  // 2D errors -- may not exist, especially for simulation       
285                                WAVE inten3 = $(folderPath+instPath+"FT"+":det_"+"FT")
286                                WAVE/Z iErr3 = $("iErr_"+"FT")                  // 2D errors -- may not exist, especially for simulation               
287                                WAVE inten4 = $(folderPath+instPath+"FB"+":det_"+"FB")
288                                WAVE/Z iErr4 = $("iErr_"+"FB")                  // 2D errors -- may not exist, especially for simulation       
289                        else
290                                Wave inten = V_getDetectorDataW(folderStr,"FL")
291                                Wave iErr = V_getDetectorDataErrW(folderStr,"FL")
292                                Wave inten2 = V_getDetectorDataW(folderStr,"FR")
293                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"FR")
294                                Wave inten3 = V_getDetectorDataW(folderStr,"FT")
295                                Wave iErr3 = V_getDetectorDataErrW(folderStr,"FT")
296                                Wave inten4 = V_getDetectorDataW(folderStr,"FB")
297                                Wave iErr4 = V_getDetectorDataErrW(folderStr,"FB")
298                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FL"+":data")
299                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FR"+":data")
300                                Wave/Z mask3 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FT"+":data")
301                                Wave/Z mask4 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FB"+":data")
302                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1 && WaveExists(mask3) == 1 && WaveExists(mask4) == 1)
303                                        maskMissing = 0
304                                endif
305                        endif   
306//                      NVAR delQ = $(folderPath+instPath+"FL"+":gDelQ_FL")
307                       
308                        Wave qTotal = $(folderPath+instPath+"FL"+":qTot_"+"FL")                 // 2D q-values 
309                        Wave qTotal2 = $(folderPath+instPath+"FR"+":qTot_"+"FR")                        // 2D q-values 
310                        Wave qTotal3 = $(folderPath+instPath+"FT"+":qTot_"+"FT")                        // 2D q-values 
311                        Wave qTotal4 = $(folderPath+instPath+"FB"+":qTot_"+"FB")                        // 2D q-values 
312
313                        Wave qx = $(folderPath+instPath+"FL"+":qx_"+"FL")                       // 2D qx-values
314                        Wave qx2 = $(folderPath+instPath+"FR"+":qx_"+"FR")                      // 2D qx-values
315                        Wave qx3 = $(folderPath+instPath+"FT"+":qx_"+"FT")                      // 2D qx-values
316                        Wave qx4 = $(folderPath+instPath+"FB"+":qx_"+"FB")                      // 2D qx-values
317
318                        Wave qy = $(folderPath+instPath+"FL"+":qy_"+"FL")                       // 2D qy-values
319                        Wave qy2 = $(folderPath+instPath+"FR"+":qy_"+"FR")                      // 2D qy-values
320                        Wave qy3 = $(folderPath+instPath+"FT"+":qy_"+"FT")                      // 2D qy-values
321                        Wave qy4 = $(folderPath+instPath+"FB"+":qy_"+"FB")                      // 2D qy-values
322                                                               
323                        nSets = 4
324                        break           
325                       
326//              case "MLR":
327//                      if(isVCALC)
328//                              WAVE inten = $(folderPath+instPath+"ML"+":det_"+"ML")
329//                              WAVE/Z iErr = $("iErr_"+"ML")                   // 2D errors -- may not exist, especially for simulation               
330//                              WAVE inten2 = $(folderPath+instPath+"MR"+":det_"+"MR")
331//                              WAVE/Z iErr2 = $("iErr_"+"MR")                  // 2D errors -- may not exist, especially for simulation       
332//                      else
333//                              Wave inten = V_getDetectorDataW(folderStr,"ML")
334//                              Wave iErr = V_getDetectorDataErrW(folderStr,"ML")
335//                              Wave inten2 = V_getDetectorDataW(folderStr,"MR")
336//                              Wave iErr2 = V_getDetectorDataErrW(folderStr,"MR")
337//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"ML"+":data")
338//                              Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MR"+":data")
339//                              if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
340//                                      maskMissing = 0
341//                              endif
342//                      endif   
343//                      NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML")
344//                     
345//                      Wave qTotal = $(folderPath+instPath+"ML"+":qTot_"+"ML")                 // 2D q-values 
346//                      Wave qTotal2 = $(folderPath+instPath+"MR"+":qTot_"+"MR")                        // 2D q-values 
347//             
348//                      nSets = 2
349//                      break                   
350//             
351//              case "MTB":
352//                      if(isVCALC)
353//                              WAVE inten = $(folderPath+instPath+"MT"+":det_"+"MT")
354//                              WAVE/Z iErr = $("iErr_"+"MT")                   // 2D errors -- may not exist, especially for simulation               
355//                              WAVE inten2 = $(folderPath+instPath+"MB"+":det_"+"MB")
356//                              WAVE/Z iErr2 = $("iErr_"+"MB")                  // 2D errors -- may not exist, especially for simulation       
357//                      else
358//                              Wave inten = V_getDetectorDataW(folderStr,"MT")
359//                              Wave iErr = V_getDetectorDataErrW(folderStr,"MT")
360//                              Wave inten2 = V_getDetectorDataW(folderStr,"MB")
361//                              Wave iErr2 = V_getDetectorDataErrW(folderStr,"MB")
362//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MT"+":data")
363//                              Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MB"+":data")
364//                              if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
365//                                      maskMissing = 0
366//                              endif
367//                      endif   
368//                      NVAR delQ = $(folderPath+instPath+"MT"+":gDelQ_MT")
369//                     
370//                      Wave qTotal = $(folderPath+instPath+"MT"+":qTot_"+"MT")                 // 2D q-values 
371//                      Wave qTotal2 = $(folderPath+instPath+"MB"+":qTot_"+"MB")                        // 2D q-values 
372//             
373//                      nSets = 2
374//                      break                           
375               
376                case "MLRTB":
377                        if(isVCALC)
378                                WAVE inten = $(folderPath+instPath+"ML"+":det_"+"ML")
379                                WAVE/Z iErr = $("iErr_"+"ML")                   // 2D errors -- may not exist, especially for simulation               
380                                WAVE inten2 = $(folderPath+instPath+"MR"+":det_"+"MR")
381                                WAVE/Z iErr2 = $("iErr_"+"MR")                  // 2D errors -- may not exist, especially for simulation       
382                                WAVE inten3 = $(folderPath+instPath+"MT"+":det_"+"MT")
383                                WAVE/Z iErr3 = $("iErr_"+"MT")                  // 2D errors -- may not exist, especially for simulation               
384                                WAVE inten4 = $(folderPath+instPath+"MB"+":det_"+"MB")
385                                WAVE/Z iErr4 = $("iErr_"+"MB")                  // 2D errors -- may not exist, especially for simulation       
386                        else
387                                Wave inten = V_getDetectorDataW(folderStr,"ML")
388                                Wave iErr = V_getDetectorDataErrW(folderStr,"ML")
389                                Wave inten2 = V_getDetectorDataW(folderStr,"MR")
390                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"MR")
391                                Wave inten3 = V_getDetectorDataW(folderStr,"MT")
392                                Wave iErr3 = V_getDetectorDataErrW(folderStr,"MT")
393                                Wave inten4 = V_getDetectorDataW(folderStr,"MB")
394                                Wave iErr4 = V_getDetectorDataErrW(folderStr,"MB")
395                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"ML"+":data")
396                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MR"+":data")
397                                Wave/Z mask3 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MT"+":data")
398                                Wave/Z mask4 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MB"+":data")
399                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1 && WaveExists(mask3) == 1 && WaveExists(mask4) == 1)
400                                        maskMissing = 0
401                                endif
402                        endif   
403//                      NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML")
404                       
405                        Wave qTotal = $(folderPath+instPath+"ML"+":qTot_"+"ML")                 // 2D q-values 
406                        Wave qTotal2 = $(folderPath+instPath+"MR"+":qTot_"+"MR")                        // 2D q-values 
407                        Wave qTotal3 = $(folderPath+instPath+"MT"+":qTot_"+"MT")                        // 2D q-values 
408                        Wave qTotal4 = $(folderPath+instPath+"MB"+":qTot_"+"MB")                        // 2D q-values 
409
410                        Wave qx = $(folderPath+instPath+"ML"+":qx_"+"ML")                       // 2D qx-values
411                        Wave qx2 = $(folderPath+instPath+"MR"+":qx_"+"MR")                      // 2D qx-values
412                        Wave qx3 = $(folderPath+instPath+"MT"+":qx_"+"MT")                      // 2D qx-values
413                        Wave qx4 = $(folderPath+instPath+"MB"+":qx_"+"MB")                      // 2D qx-values
414
415                        Wave qy = $(folderPath+instPath+"ML"+":qy_"+"ML")                       // 2D qy-values
416                        Wave qy2 = $(folderPath+instPath+"MR"+":qy_"+"MR")                      // 2D qy-values
417                        Wave qy3 = $(folderPath+instPath+"MT"+":qy_"+"MT")                      // 2D qy-values
418                        Wave qy4 = $(folderPath+instPath+"MB"+":qy_"+"MB")                      // 2D qy-values
419                                       
420                        nSets = 4
421                        break                                                                   
422                                       
423                default:
424                        nSets = 0                                                       // optional default expression executed
425                        Print "ERROR   ---- type is not recognized "
426        endswitch
427
428//      Print "delQ = ",delQ," for ",type
429
430        if(nSets == 0)
431                SetDataFolder root:
432                return(0)
433        endif
434
435
436//TODO: properly define the 2D errors here - I'll have this if I do the simulation
437// -- need to propagate the 2D errors up to this point
438//
439        if(WaveExists(iErr)==0  && WaveExists(inten) != 0)
440                Duplicate/O inten,iErr
441                Wave iErr=iErr
442//              iErr = 1+sqrt(inten+0.75)                       // can't use this -- it applies to counts, not intensity (already a count rate...)
443                iErr = sqrt(inten+0.75)                 // TODO -- here I'm just using some fictional value
444        endif
445        if(WaveExists(iErr2)==0 && WaveExists(inten2) != 0)
446                Duplicate/O inten2,iErr2
447                Wave iErr2=iErr2
448//              iErr2 = 1+sqrt(inten2+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
449                iErr2 = sqrt(inten2+0.75)                       // TODO -- here I'm just using some fictional value
450        endif
451        if(WaveExists(iErr3)==0  && WaveExists(inten3) != 0)
452                Duplicate/O inten3,iErr3
453                Wave iErr3=iErr3
454//              iErr3 = 1+sqrt(inten3+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
455                iErr3 = sqrt(inten3+0.75)                       // TODO -- here I'm just using some fictional value
456        endif
457        if(WaveExists(iErr4)==0  && WaveExists(inten4) != 0)
458                Duplicate/O inten4,iErr4
459                Wave iErr4=iErr4
460//              iErr4 = 1+sqrt(inten4+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
461                iErr4 = sqrt(inten4+0.75)                       // TODO -- here I'm just using some fictional value
462        endif
463
464        // TODO -- nq will need to be larger, once the back detector is installed
465        //
466        // note that the back panel of 320x320 (1mm res) results in 447 data points!
467        // - so I upped nq to 600
468
469        nq = 600
470
471//******TODO****** -- where to put the averaged data -- right now, folderStr is forced to ""   
472//      SetDataFolder $("root:"+folderStr)              //should already be here, but make sure...     
473        Make/O/D/N=(nq)  $(folderPath+":"+"iBin_qxqy"+"_"+type)
474//      Make/O/D/N=(nq)  $(folderPath+":"+"qBin_qxqy"+"_"+type)
475        Make/O/D/N=(nq)  $(folderPath+":"+"phiBin_qxqy"+"_"+type)
476        Make/O/D/N=(nq)  $(folderPath+":"+"nBin_qxqy"+"_"+type)
477        Make/O/D/N=(nq)  $(folderPath+":"+"iBin2_qxqy"+"_"+type)
478        Make/O/D/N=(nq)  $(folderPath+":"+"eBin_qxqy"+"_"+type)
479        Make/O/D/N=(nq)  $(folderPath+":"+"eBin2D_qxqy"+"_"+type)
480       
481        Wave iBin_qxqy = $(folderPath+":"+"iBin_qxqy_"+type)
482//      Wave qBin_qxqy = $(folderPath+":"+"qBin_qxqy"+"_"+type)
483        Wave phiBin_qxqy = $(folderPath+":"+"phiBin_qxqy"+"_"+type)
484        Wave nBin_qxqy = $(folderPath+":"+"nBin_qxqy"+"_"+type)
485        Wave iBin2_qxqy = $(folderPath+":"+"iBin2_qxqy"+"_"+type)
486        Wave eBin_qxqy = $(folderPath+":"+"eBin_qxqy"+"_"+type)
487        Wave eBin2D_qxqy = $(folderPath+":"+"eBin2D_qxqy"+"_"+type)
488       
489       
490// TODO:
491// -- get the q-binning values from the VSANS equivalent
492//      SVAR keyListStr = root:myGlobals:Protocols:gAvgInfoStr 
493//      Variable qc = NumberByKey("QCENTER",keyListStr,"=",";")         // q-center
494//      Variable nw = NumberByKey("QDELTA",keyListStr,"=",";")          // (in SANS - number of pixels wide)
495
496
497        Variable nphi,dphi,isIn,phiij,iphi
498
499// TODO: define nphi
500//      dr = 1                  //minimum annulus width, keep this fixed at one
501        NVAR numPhiSteps = root:Packages:NIST:VSANS:Globals:gNPhiSteps
502        nphi = numPhiSteps              //number of anular sectors is set by users
503        dphi = 360/nphi
504       
505
506        iBin_qxqy = 0
507        iBin2_qxqy = 0
508        eBin_qxqy = 0
509        eBin2D_qxqy = 0
510        nBin_qxqy = 0   //number of intensities added to each bin
511
512
513// 4 panels      is currently the only situation
514//
515// this needs to be a double loop now...
516// TODO:
517// -- the iErr (=2D) wave and accumulation of error is NOT CALCULATED CORRECTLY YET
518// -- the solid angle per pixel is not completely implemented.
519//    it will be present for WORK data other than RAW, but not for RAW
520
521// if any of the masks don't exist, display the error, and proceed with the averaging, using all data
522        if(maskMissing == 1)
523                Print "Mask file not found for at least one detector - so all data is used"
524        endif
525
526        Variable mask_val
527// use set 1 (no number) only
528        if(nSets >= 1)
529                xDim=DimSize(inten,0)
530                yDim=DimSize(inten,1)
531       
532                for(ii=0;ii<xDim;ii+=1)
533                        for(jj=0;jj<yDim;jj+=1)
534                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
535                                qVal = qTotal[ii][jj]
536                               
537                                isIn = V_CloseEnough(qVal,qCtr_Ann,qWidth)
538                               
539                                if(isIn)                // it's in the annulus somewhere, do something
540                                        // now I need the qx and qy to find phi
541                                        if (qy[ii][jj] >= 0)
542                                                //phiij is in degrees
543                                                phiij = atan2(qy[ii][jj],qx[ii][jj])*180/Pi             //0 to 180 deg
544                                        else
545                                                phiij = 360 + atan2(qy[ii][jj],qx[ii][jj])*180/Pi               //180 to 360 deg
546                                        Endif
547                                        if (phiij > (360-0.5*dphi))
548                                                phiij -= 360
549                                        Endif
550                                        iphi = trunc(phiij/dphi + 1.501)                        // TODO: why the value of 1.501????
551                                                       
552                                        val = inten[ii][jj]
553                                       
554                                        if(isVCALC || maskMissing)              // mask_val == 0 == keep, mask_val == 1 = YES, mask out the point
555                                                mask_val = 0
556                                        else
557                                                mask_val = mask[ii][jj]
558                                        endif
559                                        if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
560                                                iBin_qxqy[iphi-1] += val
561                                                iBin2_qxqy[iphi-1] += val*val
562                                                eBin2D_qxqy[iphi-1] += iErr[ii][jj]*iErr[ii][jj]
563                                                nBin_qxqy[iphi-1] += 1
564                                        endif
565                               
566                                endif // isIn
567                               
568                        endfor
569                endfor
570               
571        endif
572
573// add in set 2 (set 1 already done)
574        if(nSets >= 2)
575                xDim=DimSize(inten2,0)
576                yDim=DimSize(inten2,1)
577       
578                for(ii=0;ii<xDim;ii+=1)
579                        for(jj=0;jj<yDim;jj+=1)
580                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
581                                qVal = qTotal2[ii][jj]
582                               
583                                isIn = V_CloseEnough(qVal,qCtr_Ann,qWidth)
584                               
585                                if(isIn)                // it's in the annulus somewhere, do something
586                                        // now I need the qx and qy to find phi
587                                        if (qy2[ii][jj] >= 0)
588                                                //phiij is in degrees
589                                                phiij = atan2(qy2[ii][jj],qx2[ii][jj])*180/Pi           //0 to 180 deg
590                                        else
591                                                phiij = 360 + atan2(qy2[ii][jj],qx2[ii][jj])*180/Pi             //180 to 360 deg
592                                        Endif
593                                        if (phiij > (360-0.5*dphi))
594                                                phiij -= 360
595                                        Endif
596                                        iphi = trunc(phiij/dphi + 1.501)                        // TODO: why the value of 1.501????
597                                                       
598                                        val = inten2[ii][jj]
599                                       
600                                        if(isVCALC || maskMissing)              // mask_val == 0 == keep, mask_val == 1 = YES, mask out the point
601                                                mask_val = 0
602                                        else
603                                                mask_val = mask2[ii][jj]
604                                        endif
605                                        if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
606                                                iBin_qxqy[iphi-1] += val
607                                                iBin2_qxqy[iphi-1] += val*val
608                                                eBin2D_qxqy[iphi-1] += iErr2[ii][jj]*iErr2[ii][jj]
609                                                nBin_qxqy[iphi-1] += 1
610                                        endif
611                               
612                                endif // isIn
613
614                        endfor          //jj
615                endfor          //ii
616               
617        endif           // set 2
618
619
620// add in set 3 and 4 (set 1 and 2 already done)
621        if(nSets == 4)
622                xDim=DimSize(inten3,0)
623                yDim=DimSize(inten3,1)
624       
625                for(ii=0;ii<xDim;ii+=1)
626                        for(jj=0;jj<yDim;jj+=1)
627                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
628                                qVal = qTotal3[ii][jj]
629                               
630                                isIn = V_CloseEnough(qVal,qCtr_Ann,qWidth)
631                               
632                                if(isIn)                // it's in the annulus somewhere, do something
633                                        // now I need the qx and qy to find phi
634                                        if (qy3[ii][jj] >= 0)
635                                                //phiij is in degrees
636                                                phiij = atan2(qy3[ii][jj],qx3[ii][jj])*180/Pi           //0 to 180 deg
637                                        else
638                                                phiij = 360 + atan2(qy3[ii][jj],qx3[ii][jj])*180/Pi             //180 to 360 deg
639                                        Endif
640                                        if (phiij > (360-0.5*dphi))
641                                                phiij -= 360
642                                        Endif
643                                        iphi = trunc(phiij/dphi + 1.501)                        // TODO: why the value of 1.501????
644                                                       
645                                        val = inten3[ii][jj]
646                                       
647                                        if(isVCALC || maskMissing)              // mask_val == 0 == keep, mask_val == 1 = YES, mask out the point
648                                                mask_val = 0
649                                        else
650                                                mask_val = mask3[ii][jj]
651                                        endif
652                                        if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
653                                                iBin_qxqy[iphi-1] += val
654                                                iBin2_qxqy[iphi-1] += val*val
655                                                eBin2D_qxqy[iphi-1] += iErr3[ii][jj]*iErr3[ii][jj]
656                                                nBin_qxqy[iphi-1] += 1
657                                        endif
658                               
659                                endif // isIn
660                               
661                       
662                        endfor
663                endfor          // end of ij loop over set 3
664               
665               
666                xDim=DimSize(inten4,0)
667                yDim=DimSize(inten4,1)
668       
669                for(ii=0;ii<xDim;ii+=1)
670                        for(jj=0;jj<yDim;jj+=1)
671                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
672                                qVal = qTotal4[ii][jj]
673                               
674                                isIn = V_CloseEnough(qVal,qCtr_Ann,qWidth)
675                               
676                                if(isIn)                // it's in the annulus somewhere, do something
677                                        // now I need the qx and qy to find phi
678                                        if (qy4[ii][jj] >= 0)
679                                                //phiij is in degrees
680                                                phiij = atan2(qy4[ii][jj],qx4[ii][jj])*180/Pi           //0 to 180 deg
681                                        else
682                                                phiij = 360 + atan2(qy4[ii][jj],qx4[ii][jj])*180/Pi             //180 to 360 deg
683                                        Endif
684                                        if (phiij > (360-0.5*dphi))
685                                                phiij -= 360
686                                        Endif
687                                        iphi = trunc(phiij/dphi + 1.501)                        // TODO: why the value of 1.501????
688                                                       
689                                        val = inten4[ii][jj]
690                                       
691                                        if(isVCALC || maskMissing)              // mask_val == 0 == keep, mask_val == 1 = YES, mask out the point
692                                                mask_val = 0
693                                        else
694                                                mask_val = mask4[ii][jj]
695                                        endif
696                                        if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
697                                                iBin_qxqy[iphi-1] += val
698                                                iBin2_qxqy[iphi-1] += val*val
699                                                eBin2D_qxqy[iphi-1] += iErr4[ii][jj]*iErr4[ii][jj]
700                                                nBin_qxqy[iphi-1] += 1
701                                        endif
702                               
703                                endif // isIn                           
704                               
705                               
706                        endfor
707                endfor          // end ij loop over set 4
708               
709        endif   // adding sets 3 and 4
710
711
712// after looping through all of the data on the panels, calculate errors on I(q),
713// just like in CircSectAve.ipf
714// TODO:
715// -- 2D Errors were NOT properly acculumated through reduction, so this loop of calculations is NOT MEANINGFUL (yet)
716// x- the error on the 1D intensity, is correctly calculated as the standard error of the mean.
717        for(ii=0;ii<nphi;ii+=1)
718       
719                phiBin_qxqy[ii] = dphi*ii
720               
721                if(nBin_qxqy[ii] == 0)
722                        //no pixels in annuli, data unknown
723                        iBin_qxqy[ii] = 0
724                        eBin_qxqy[ii] = 1
725                        eBin2D_qxqy[ii] = NaN
726                else
727                        if(nBin_qxqy[ii] <= 1)
728                                //need more than one pixel to determine error
729                                iBin_qxqy[ii] /= nBin_qxqy[ii]
730                                eBin_qxqy[ii] = 1
731                                eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
732                        else
733                                //assume that the intensity in each pixel in annuli is normally distributed about mean...
734                                //  -- this is correctly calculating the error as the standard error of the mean, as
735                                //    was always done for SANS as well.
736                                iBin_qxqy[ii] /= nBin_qxqy[ii]
737                                avesq = iBin_qxqy[ii]^2
738                                aveisq = iBin2_qxqy[ii]/nBin_qxqy[ii]
739                                var = aveisq-avesq
740                                if(var<=0)
741                                        eBin_qxqy[ii] = 1e-6
742                                else
743                                        eBin_qxqy[ii] = sqrt(var/(nBin_qxqy[ii] - 1))
744                                endif
745                                // and calculate as it is propagated pixel-by-pixel
746                                eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
747                        endif
748                endif
749        endfor
750       
751        eBin2D_qxqy = sqrt(eBin2D_qxqy)         // as equation (3) of John's memo
752       
753        // find the last non-zero point, working backwards
754        val=nq
755        do
756                val -= 1
757        while((nBin_qxqy[val] == 0) && val > 0)
758       
759//      print val, nBin_qxqy[val]
760        DeletePoints val, nq-val, iBin_qxqy,phiBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
761
762        if(val == 0)
763                // all the points were deleted
764                return(0)
765        endif
766       
767       
768        // since the beam center is not always on the detector, many of the low Q bins will have zero pixels
769        // find the first non-zero point, working forwards
770        val = -1
771        do
772                val += 1
773        while(nBin_qxqy[val] == 0)     
774        DeletePoints 0, val, iBin_qxqy,phiBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
775
776        // ?? there still may be a point in the q-range that gets zero pixel contribution - so search this out and get rid of it
777        val = numpnts(nBin_qxqy)-1
778        do
779                if(nBin_qxqy[val] == 0)
780                        DeletePoints val, 1, iBin_qxqy,phiBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
781                endif
782                val -= 1
783        while(val>0)
784       
785        // TODO:
786        // -- is this where I do the resolution calculation? This is where I calculate the resolution in SANS (see CircSectAve)
787        // -- or do I do it as a separate call?
788        // -- use the isVCALC flag to exclude VCALC from the resolution calculation if necessary
789        //
790       
791       
792       
793        SetDataFolder root:
794       
795        return(0)
796End
797
798
Note: See TracBrowser for help on using the repository browser.