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

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

corrected compiling issue with version 7.87

Fixed logic in defining sector mask
Added mask for Annular average

File size: 27.8 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 degree circle
9//
10// qWidth is +/- around the q-center
11//
12//
13// TODO
14// x- add error bars to the plot of phi
15// x- data writer to export annular data (3-column)
16// x- loader to re-read annular data (will the normal loader work?) -- yes
17// -- integrate this with the protocol
18// -- integrate this with a more general "average panel"
19// -- draw the q-center and width on the image (as a contour? - if no, use a draw layer...)
20// -- Can I "flag" the pixels that contribute to the annualr average, and overlay them like a
21//    Mask(translucent?), like a "thresholding" operation, but act on the Q-value, not the Z- value.
22//
23// TODO!!!
24// x- none of these procedures are aware of the BACK detector
25//
26
27//
28//
29Proc Annular_Binning(folderStr,detGroup,qCtr_Ann,qWidth)
30        String folderStr="SAM",detGroup="F"
31        Variable qCtr_Ann=0.1,qWidth=0.01       // +/- in A^-1
32       
33        V_QBinAllPanels_Annular(folderStr,detGroup,qCtr_Ann,qWidth)
34
35        V_Phi_Graph_Proc(folderStr,detGroup)
36
37End
38
39
40// TODO -- binType == 4 (slit mode) should never end up here, as it makes no sense
41//
42// -- really, the only binning that makes any sense is "one", treating each panel individually,
43// so I may scrap the parameter, or ignore it. so don't count on it in the future.
44//
45Function V_QBinAllPanels_Annular(folderStr,detGroup,qCtr_Ann,qWidth)
46        String folderStr,detGroup
47        Variable qCtr_Ann,qWidth        // +/- in A^-1
48
49        Variable ii,delQ
50        String detStr
51       
52        strswitch(detGroup)
53                case "F":
54                        detStr = "FLRTB"
55                        break
56                case "M":
57                        detStr = "MLRTB"
58                        break
59                case "B":
60                        detStr = "B"
61                        break
62                default:
63                        DoAlert 0,"No detGroup match in V_QBinAllPanels_Annular"
64                        return(0)
65        endswitch
66
67       
68// right now, use all of the detectors. There is a lot of waste in this and it could be
69// done a lot faster, but...
70
71// TODO         
72                // detStr = "FLRTB" or "MLRTB", depending which panel the q-ring is centered on/
73                // for now, no crossing of the rings onto different carriages
74               
75        V_fDoAnnularBin_QxQy2D(folderStr,detStr,qCtr_Ann,qWidth)
76
77
78        return(0)
79End
80
81Function V_Phi_Graph_Proc(folderStr,detGroup)
82        String folderStr,detGroup
83       
84        PauseUpdate; Silent 1           // building window...
85        String fldrSav0= GetDataFolder(1)
86        SetDataFolder $("root:Packages:NIST:VSANS:"+folderStr)
87
88        DoWindow/F V_Phi_Graph
89        if(V_flag == 0)
90                Display /W=(35,45,572,419)/N=V_Phi_Graph /K=1
91        else
92                RemoveFromGraph/Z iPhiBin_qxqy_FLRTB
93                RemoveFromGraph/Z iPhiBin_qxqy_MLRTB
94                RemoveFromGraph/Z iPhiBin_qxqy_B
95        endif
96
97        if(cmpstr(detGroup,"F") == 0)
98                wave iPhiBin_qxqy_FLRTB=iPhiBin_qxqy_FLRTB
99                wave phiBin_qxqy_FLRTB=phiBin_qxqy_FLRTB
100                wave ePhiBin_qxqy_FLRTB=ePhiBin_qxqy_FLRTB
101               
102                AppendToGraph iPhiBin_qxqy_FLRTB vs phiBin_qxqy_FLRTB
103//              Display /W=(35,45,572,419)/N=V_Phi_Graph /K=1 iPhiBin_qxqy_FLRTB vs phiBin_qxqy_FLRTB
104                ModifyGraph mode=4
105                ModifyGraph marker=19
106                ErrorBars iPhiBin_qxqy_FLRTB Y,wave=(ePhiBin_qxqy_FLRTB,ePhiBin_qxqy_FLRTB)
107                Label left "Counts"
108                Label bottom "Phi (degrees)"
109                Legend
110        endif
111        if(cmpstr(detGroup,"M") == 0)
112                wave iPhiBin_qxqy_MLRTB=iPhiBin_qxqy_MLRTB
113                wave phiBin_qxqy_MLRTB=phiBin_qxqy_MLRTB
114                wave ePhiBin_qxqy_MLRTB=ePhiBin_qxqy_MLRTB
115               
116                AppendToGraph iPhiBin_qxqy_MLRTB vs phiBin_qxqy_MLRTB
117                ModifyGraph mode=4
118                ModifyGraph marker=19
119                ErrorBars iPhiBin_qxqy_MLRTB Y,wave=(ePhiBin_qxqy_MLRTB,ePhiBin_qxqy_MLRTB)
120                Label left "Counts"
121                Label bottom "Phi (degrees)"
122                Legend
123        endif
124        if(cmpstr(detGroup,"B") == 0)
125                wave iPhiBin_qxqy_B=iPhiBin_qxqy_B
126                wave phiBin_qxqy_B=phiBin_qxqy_B
127                wave ePhiBin_qxqy_B=ePhiBin_qxqy_B
128               
129                AppendToGraph iPhiBin_qxqy_B vs phiBin_qxqy_B
130                ModifyGraph mode=4
131                ModifyGraph marker=19
132                ErrorBars iPhiBin_qxqy_B Y,wave=(ePhiBin_qxqy_B,ePhiBin_qxqy_B)
133                Label left "Counts"
134                Label bottom "Phi (degrees)"
135                Legend
136        endif
137       
138        SetDataFolder fldrSav0
139EndMacro
140
141
142//////////
143//
144//
145//
146// TODO
147// -- "iErr" is not always defined correctly since it doesn't really apply here for data that is not 2D simulation
148//
149// ** Currently, type is being passed in as "" and ignored (looping through all of the detector panels
150// to potentially add to the annular bins)
151//
152// folderStr = WORK folder, type = the binning type (may include multiple detectors)
153//
154//      Variable qCtr_Ann = 0.1
155//      Variable qWidth = 0.02          // +/- in A^-1
156//
157Function V_fDoAnnularBin_QxQy2D(folderStr,type,qCtr_Ann,qWidth)
158        String folderStr,type
159        Variable qCtr_Ann,qWidth
160       
161        Variable nSets = 0
162        Variable xDim,yDim
163        Variable ii,jj
164        Variable qVal,nq,var,avesq,aveisq
165        Variable binIndex,val,isVCALC=0,maskMissing
166
167        String folderPath = "root:Packages:NIST:VSANS:"+folderStr
168        String instPath = ":entry:instrument:detector_"
169        String detStr
170               
171        if(cmpstr(folderStr,"VCALC") == 0)
172                isVCALC = 1
173        endif
174       
175// now switch on the type to determine which waves to declare and create
176// since there may be more than one panel to step through. There may be two, there may be four
177//
178
179// TODO:
180// -- Solid_Angle -- waves will be present for WORK data other than RAW, but not for RAW
181//
182// assume that the mask files are missing unless we can find them. If VCALC data,
183//  then the Mask is missing by definition
184        maskMissing = 1
185
186        strswitch(type) // string switch
187
188               
189                case "FLRTB":
190                        if(isVCALC)
191                                WAVE inten = $(folderPath+instPath+"FL"+":det_"+"FL")
192                                WAVE/Z iErr = $("iErr_"+"FL")                   // 2D errors -- may not exist, especially for simulation               
193                                WAVE inten2 = $(folderPath+instPath+"FR"+":det_"+"FR")
194                                WAVE/Z iErr2 = $("iErr_"+"FR")                  // 2D errors -- may not exist, especially for simulation       
195                                WAVE inten3 = $(folderPath+instPath+"FT"+":det_"+"FT")
196                                WAVE/Z iErr3 = $("iErr_"+"FT")                  // 2D errors -- may not exist, especially for simulation               
197                                WAVE inten4 = $(folderPath+instPath+"FB"+":det_"+"FB")
198                                WAVE/Z iErr4 = $("iErr_"+"FB")                  // 2D errors -- may not exist, especially for simulation       
199                        else
200                                Wave inten = V_getDetectorDataW(folderStr,"FL")
201                                Wave iErr = V_getDetectorDataErrW(folderStr,"FL")
202                                Wave inten2 = V_getDetectorDataW(folderStr,"FR")
203                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"FR")
204                                Wave inten3 = V_getDetectorDataW(folderStr,"FT")
205                                Wave iErr3 = V_getDetectorDataErrW(folderStr,"FT")
206                                Wave inten4 = V_getDetectorDataW(folderStr,"FB")
207                                Wave iErr4 = V_getDetectorDataErrW(folderStr,"FB")
208                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FL"+":data")
209                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FR"+":data")
210                                Wave/Z mask3 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FT"+":data")
211                                Wave/Z mask4 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FB"+":data")
212                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1 && WaveExists(mask3) == 1 && WaveExists(mask4) == 1)
213                                        maskMissing = 0
214                                endif
215                        endif   
216//                      NVAR delQ = $(folderPath+instPath+"FL"+":gDelQ_FL")
217                       
218                        Wave qTotal = $(folderPath+instPath+"FL"+":qTot_"+"FL")                 // 2D q-values 
219                        Wave qTotal2 = $(folderPath+instPath+"FR"+":qTot_"+"FR")                        // 2D q-values 
220                        Wave qTotal3 = $(folderPath+instPath+"FT"+":qTot_"+"FT")                        // 2D q-values 
221                        Wave qTotal4 = $(folderPath+instPath+"FB"+":qTot_"+"FB")                        // 2D q-values 
222
223                        Wave qx = $(folderPath+instPath+"FL"+":qx_"+"FL")                       // 2D qx-values
224                        Wave qx2 = $(folderPath+instPath+"FR"+":qx_"+"FR")                      // 2D qx-values
225                        Wave qx3 = $(folderPath+instPath+"FT"+":qx_"+"FT")                      // 2D qx-values
226                        Wave qx4 = $(folderPath+instPath+"FB"+":qx_"+"FB")                      // 2D qx-values
227
228                        Wave qy = $(folderPath+instPath+"FL"+":qy_"+"FL")                       // 2D qy-values
229                        Wave qy2 = $(folderPath+instPath+"FR"+":qy_"+"FR")                      // 2D qy-values
230                        Wave qy3 = $(folderPath+instPath+"FT"+":qy_"+"FT")                      // 2D qy-values
231                        Wave qy4 = $(folderPath+instPath+"FB"+":qy_"+"FB")                      // 2D qy-values
232                                                               
233                        nSets = 4
234                        break           
235                       
236       
237               
238                case "MLRTB":
239                        if(isVCALC)
240                                WAVE inten = $(folderPath+instPath+"ML"+":det_"+"ML")
241                                WAVE/Z iErr = $("iErr_"+"ML")                   // 2D errors -- may not exist, especially for simulation               
242                                WAVE inten2 = $(folderPath+instPath+"MR"+":det_"+"MR")
243                                WAVE/Z iErr2 = $("iErr_"+"MR")                  // 2D errors -- may not exist, especially for simulation       
244                                WAVE inten3 = $(folderPath+instPath+"MT"+":det_"+"MT")
245                                WAVE/Z iErr3 = $("iErr_"+"MT")                  // 2D errors -- may not exist, especially for simulation               
246                                WAVE inten4 = $(folderPath+instPath+"MB"+":det_"+"MB")
247                                WAVE/Z iErr4 = $("iErr_"+"MB")                  // 2D errors -- may not exist, especially for simulation       
248                        else
249                                Wave inten = V_getDetectorDataW(folderStr,"ML")
250                                Wave iErr = V_getDetectorDataErrW(folderStr,"ML")
251                                Wave inten2 = V_getDetectorDataW(folderStr,"MR")
252                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"MR")
253                                Wave inten3 = V_getDetectorDataW(folderStr,"MT")
254                                Wave iErr3 = V_getDetectorDataErrW(folderStr,"MT")
255                                Wave inten4 = V_getDetectorDataW(folderStr,"MB")
256                                Wave iErr4 = V_getDetectorDataErrW(folderStr,"MB")
257                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"ML"+":data")
258                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MR"+":data")
259                                Wave/Z mask3 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MT"+":data")
260                                Wave/Z mask4 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MB"+":data")
261                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1 && WaveExists(mask3) == 1 && WaveExists(mask4) == 1)
262                                        maskMissing = 0
263                                endif
264                        endif   
265//                      NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML")
266                       
267                        Wave qTotal = $(folderPath+instPath+"ML"+":qTot_"+"ML")                 // 2D q-values 
268                        Wave qTotal2 = $(folderPath+instPath+"MR"+":qTot_"+"MR")                        // 2D q-values 
269                        Wave qTotal3 = $(folderPath+instPath+"MT"+":qTot_"+"MT")                        // 2D q-values 
270                        Wave qTotal4 = $(folderPath+instPath+"MB"+":qTot_"+"MB")                        // 2D q-values 
271
272                        Wave qx = $(folderPath+instPath+"ML"+":qx_"+"ML")                       // 2D qx-values
273                        Wave qx2 = $(folderPath+instPath+"MR"+":qx_"+"MR")                      // 2D qx-values
274                        Wave qx3 = $(folderPath+instPath+"MT"+":qx_"+"MT")                      // 2D qx-values
275                        Wave qx4 = $(folderPath+instPath+"MB"+":qx_"+"MB")                      // 2D qx-values
276
277                        Wave qy = $(folderPath+instPath+"ML"+":qy_"+"ML")                       // 2D qy-values
278                        Wave qy2 = $(folderPath+instPath+"MR"+":qy_"+"MR")                      // 2D qy-values
279                        Wave qy3 = $(folderPath+instPath+"MT"+":qy_"+"MT")                      // 2D qy-values
280                        Wave qy4 = $(folderPath+instPath+"MB"+":qy_"+"MB")                      // 2D qy-values
281                                       
282                        nSets = 4
283                        break                                                                   
284                                       
285                case "B":
286                        if(isVCALC)
287                                WAVE inten = $(folderPath+instPath+"B"+":det_"+"B")
288                                WAVE/Z iErr = $("iErr_"+"B")                    // 2D errors -- may not exist, especially for simulation               
289                        else
290                                Wave inten = V_getDetectorDataW(folderStr,"B")
291                                Wave iErr = V_getDetectorDataErrW(folderStr,"B")
292
293                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"B"+":data")
294                               
295                                if(WaveExists(mask) == 1)
296                                        maskMissing = 0
297                                endif
298                        endif   
299//                      NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML")
300                       
301                        Wave qTotal = $(folderPath+instPath+"B"+":qTot_"+"B")                   // 2D q-values 
302       
303                        Wave qx = $(folderPath+instPath+"B"+":qx_"+"B")                 // 2D qx-values
304
305                        Wave qy = $(folderPath+instPath+"B"+":qy_"+"B")                 // 2D qy-values
306                                       
307                        nSets = 1
308                        break                                                           
309
310                                       
311                default:
312                        nSets = 0                                                       // optional default expression executed
313                        Print "ERROR   ---- type is not recognized "
314        endswitch
315
316//      Print "delQ = ",delQ," for ",type
317
318        if(nSets == 0)
319                SetDataFolder root:
320                return(0)
321        endif
322
323
324//TODO: properly define the 2D errors here - I'll have this if I do the simulation
325// -- need to propagate the 2D errors up to this point
326//
327        if(WaveExists(iErr)==0  && WaveExists(inten) != 0)
328                Duplicate/O inten,iErr
329                Wave iErr=iErr
330//              iErr = 1+sqrt(inten+0.75)                       // can't use this -- it applies to counts, not intensity (already a count rate...)
331                iErr = sqrt(inten+0.75)                 // TODO -- here I'm just using some fictional value
332        endif
333        if(WaveExists(iErr2)==0 && WaveExists(inten2) != 0)
334                Duplicate/O inten2,iErr2
335                Wave iErr2=iErr2
336//              iErr2 = 1+sqrt(inten2+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
337                iErr2 = sqrt(inten2+0.75)                       // TODO -- here I'm just using some fictional value
338        endif
339        if(WaveExists(iErr3)==0  && WaveExists(inten3) != 0)
340                Duplicate/O inten3,iErr3
341                Wave iErr3=iErr3
342//              iErr3 = 1+sqrt(inten3+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
343                iErr3 = sqrt(inten3+0.75)                       // TODO -- here I'm just using some fictional value
344        endif
345        if(WaveExists(iErr4)==0  && WaveExists(inten4) != 0)
346                Duplicate/O inten4,iErr4
347                Wave iErr4=iErr4
348//              iErr4 = 1+sqrt(inten4+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
349                iErr4 = sqrt(inten4+0.75)                       // TODO -- here I'm just using some fictional value
350        endif
351
352        // TODO -- nq will need to be larger, once the back detector is installed
353        //
354        // note that the back panel of 320x320 (1mm res) results in 447 data points!
355        // - so I upped nq to 600
356
357        nq = 600
358
359//******TODO****** -- where to put the averaged data -- right now, folderStr is forced to ""   
360//      SetDataFolder $("root:"+folderStr)              //should already be here, but make sure...     
361        Make/O/D/N=(nq)  $(folderPath+":"+"iPhiBin_qxqy"+"_"+type)
362//      Make/O/D/N=(nq)  $(folderPath+":"+"qBin_qxqy"+"_"+type)
363        Make/O/D/N=(nq)  $(folderPath+":"+"phiBin_qxqy"+"_"+type)
364        Make/O/D/N=(nq)  $(folderPath+":"+"nPhiBin_qxqy"+"_"+type)
365        Make/O/D/N=(nq)  $(folderPath+":"+"iPhiBin2_qxqy"+"_"+type)
366        Make/O/D/N=(nq)  $(folderPath+":"+"ePhiBin_qxqy"+"_"+type)
367        Make/O/D/N=(nq)  $(folderPath+":"+"ePhiBin2D_qxqy"+"_"+type)
368       
369        Wave iPhiBin_qxqy = $(folderPath+":"+"iPhiBin_qxqy_"+type)
370//      Wave qBin_qxqy = $(folderPath+":"+"qBin_qxqy"+"_"+type)
371        Wave phiBin_qxqy = $(folderPath+":"+"phiBin_qxqy"+"_"+type)
372        Wave nPhiBin_qxqy = $(folderPath+":"+"nPhiBin_qxqy"+"_"+type)
373        Wave iPhiBin2_qxqy = $(folderPath+":"+"iPhiBin2_qxqy"+"_"+type)
374        Wave ePhiBin_qxqy = $(folderPath+":"+"ePhiBin_qxqy"+"_"+type)
375        Wave ePhiBin2D_qxqy = $(folderPath+":"+"ePhiBin2D_qxqy"+"_"+type)
376       
377       
378// TODO:
379// -- get the q-binning values from the VSANS equivalent
380//      SVAR keyListStr = root:myGlobals:Protocols:gAvgInfoStr 
381//      Variable qc = NumberByKey("QCENTER",keyListStr,"=",";")         // q-center
382//      Variable nw = NumberByKey("QDELTA",keyListStr,"=",";")          // (in SANS - number of pixels wide)
383
384
385        Variable nphi,dphi,isIn,phiij,iphi
386
387// TODO: define nphi
388//      dr = 1                  //minimum annulus width, keep this fixed at one
389        NVAR numPhiSteps = root:Packages:NIST:VSANS:Globals:gNPhiSteps
390        nphi = numPhiSteps              //number of anular sectors is set by users
391        dphi = 360/nphi
392       
393
394        iPhiBin_qxqy = 0
395        iPhiBin2_qxqy = 0
396        ePhiBin_qxqy = 0
397        ePhiBin2D_qxqy = 0
398        nPhiBin_qxqy = 0        //number of intensities added to each bin
399
400
401// 4 panels      is currently the only situation
402//
403// this needs to be a double loop now...
404// TODO:
405// -- the iErr (=2D) wave and accumulation of error is NOT CALCULATED CORRECTLY YET
406// -- the solid angle per pixel is not completely implemented.
407//    it will be present for WORK data other than RAW, but not for RAW
408
409// if any of the masks don't exist, display the error, and proceed with the averaging, using all data
410        if(maskMissing == 1)
411                Print "Mask file not found for at least one detector - so all data is used"
412        endif
413
414        Variable mask_val
415// use set 1 (no number) only
416        if(nSets >= 1)
417                xDim=DimSize(inten,0)
418                yDim=DimSize(inten,1)
419       
420                for(ii=0;ii<xDim;ii+=1)
421                        for(jj=0;jj<yDim;jj+=1)
422                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
423                                qVal = qTotal[ii][jj]
424                               
425                                isIn = V_CloseEnough(qVal,qCtr_Ann,qWidth)
426                               
427                                if(isIn)                // it's in the annulus somewhere, do something
428                                        // now I need the qx and qy to find phi
429                                        if (qy[ii][jj] >= 0)
430                                                //phiij is in degrees
431                                                phiij = atan2(qy[ii][jj],qx[ii][jj])*180/Pi             //0 to 180 deg
432                                        else
433                                                phiij = 360 + atan2(qy[ii][jj],qx[ii][jj])*180/Pi               //180 to 360 deg
434                                        Endif
435                                        if (phiij > (360-0.5*dphi))
436                                                phiij -= 360
437                                        Endif
438                                        iphi = trunc(phiij/dphi + 1.501)                        // TODO: why the value of 1.501????
439                                                       
440                                        val = inten[ii][jj]
441                                       
442                                        if(isVCALC || maskMissing)              // mask_val == 0 == keep, mask_val == 1 = YES, mask out the point
443                                                mask_val = 0
444                                        else
445                                                mask_val = mask[ii][jj]
446                                        endif
447                                        if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
448                                                iPhiBin_qxqy[iphi-1] += val
449                                                iPhiBin2_qxqy[iphi-1] += val*val
450                                                ePhiBin2D_qxqy[iphi-1] += iErr[ii][jj]*iErr[ii][jj]
451                                                nPhiBin_qxqy[iphi-1] += 1
452                                        endif
453                               
454                                endif // isIn
455                               
456                        endfor
457                endfor
458               
459        endif
460
461// add in set 2 (set 1 already done)
462        if(nSets >= 2)
463                xDim=DimSize(inten2,0)
464                yDim=DimSize(inten2,1)
465       
466                for(ii=0;ii<xDim;ii+=1)
467                        for(jj=0;jj<yDim;jj+=1)
468                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
469                                qVal = qTotal2[ii][jj]
470                               
471                                isIn = V_CloseEnough(qVal,qCtr_Ann,qWidth)
472                               
473                                if(isIn)                // it's in the annulus somewhere, do something
474                                        // now I need the qx and qy to find phi
475                                        if (qy2[ii][jj] >= 0)
476                                                //phiij is in degrees
477                                                phiij = atan2(qy2[ii][jj],qx2[ii][jj])*180/Pi           //0 to 180 deg
478                                        else
479                                                phiij = 360 + atan2(qy2[ii][jj],qx2[ii][jj])*180/Pi             //180 to 360 deg
480                                        Endif
481                                        if (phiij > (360-0.5*dphi))
482                                                phiij -= 360
483                                        Endif
484                                        iphi = trunc(phiij/dphi + 1.501)                        // TODO: why the value of 1.501????
485                                                       
486                                        val = inten2[ii][jj]
487                                       
488                                        if(isVCALC || maskMissing)              // mask_val == 0 == keep, mask_val == 1 = YES, mask out the point
489                                                mask_val = 0
490                                        else
491                                                mask_val = mask2[ii][jj]
492                                        endif
493                                        if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
494                                                iPhiBin_qxqy[iphi-1] += val
495                                                iPhiBin2_qxqy[iphi-1] += val*val
496                                                ePhiBin2D_qxqy[iphi-1] += iErr2[ii][jj]*iErr2[ii][jj]
497                                                nPhiBin_qxqy[iphi-1] += 1
498                                        endif
499                               
500                                endif // isIn
501
502                        endfor          //jj
503                endfor          //ii
504               
505        endif           // set 2
506
507
508// add in set 3 and 4 (set 1 and 2 already done)
509        if(nSets == 4)
510                xDim=DimSize(inten3,0)
511                yDim=DimSize(inten3,1)
512       
513                for(ii=0;ii<xDim;ii+=1)
514                        for(jj=0;jj<yDim;jj+=1)
515                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
516                                qVal = qTotal3[ii][jj]
517                               
518                                isIn = V_CloseEnough(qVal,qCtr_Ann,qWidth)
519                               
520                                if(isIn)                // it's in the annulus somewhere, do something
521                                        // now I need the qx and qy to find phi
522                                        if (qy3[ii][jj] >= 0)
523                                                //phiij is in degrees
524                                                phiij = atan2(qy3[ii][jj],qx3[ii][jj])*180/Pi           //0 to 180 deg
525                                        else
526                                                phiij = 360 + atan2(qy3[ii][jj],qx3[ii][jj])*180/Pi             //180 to 360 deg
527                                        Endif
528                                        if (phiij > (360-0.5*dphi))
529                                                phiij -= 360
530                                        Endif
531                                        iphi = trunc(phiij/dphi + 1.501)                        // TODO: why the value of 1.501????
532                                                       
533                                        val = inten3[ii][jj]
534                                       
535                                        if(isVCALC || maskMissing)              // mask_val == 0 == keep, mask_val == 1 = YES, mask out the point
536                                                mask_val = 0
537                                        else
538                                                mask_val = mask3[ii][jj]
539                                        endif
540                                        if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
541                                                iPhiBin_qxqy[iphi-1] += val
542                                                iPhiBin2_qxqy[iphi-1] += val*val
543                                                ePhiBin2D_qxqy[iphi-1] += iErr3[ii][jj]*iErr3[ii][jj]
544                                                nPhiBin_qxqy[iphi-1] += 1
545                                        endif
546                               
547                                endif // isIn
548                               
549                       
550                        endfor
551                endfor          // end of ij loop over set 3
552               
553               
554                xDim=DimSize(inten4,0)
555                yDim=DimSize(inten4,1)
556       
557                for(ii=0;ii<xDim;ii+=1)
558                        for(jj=0;jj<yDim;jj+=1)
559                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
560                                qVal = qTotal4[ii][jj]
561                               
562                                isIn = V_CloseEnough(qVal,qCtr_Ann,qWidth)
563                               
564                                if(isIn)                // it's in the annulus somewhere, do something
565                                        // now I need the qx and qy to find phi
566                                        if (qy4[ii][jj] >= 0)
567                                                //phiij is in degrees
568                                                phiij = atan2(qy4[ii][jj],qx4[ii][jj])*180/Pi           //0 to 180 deg
569                                        else
570                                                phiij = 360 + atan2(qy4[ii][jj],qx4[ii][jj])*180/Pi             //180 to 360 deg
571                                        Endif
572                                        if (phiij > (360-0.5*dphi))
573                                                phiij -= 360
574                                        Endif
575                                        iphi = trunc(phiij/dphi + 1.501)                        // TODO: why the value of 1.501????
576                                                       
577                                        val = inten4[ii][jj]
578                                       
579                                        if(isVCALC || maskMissing)              // mask_val == 0 == keep, mask_val == 1 = YES, mask out the point
580                                                mask_val = 0
581                                        else
582                                                mask_val = mask4[ii][jj]
583                                        endif
584                                        if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
585                                                iPhiBin_qxqy[iphi-1] += val
586                                                iPhiBin2_qxqy[iphi-1] += val*val
587                                                ePhiBin2D_qxqy[iphi-1] += iErr4[ii][jj]*iErr4[ii][jj]
588                                                nPhiBin_qxqy[iphi-1] += 1
589                                        endif
590                               
591                                endif // isIn                           
592                               
593                               
594                        endfor
595                endfor          // end ij loop over set 4
596               
597        endif   // adding sets 3 and 4
598
599
600// after looping through all of the data on the panels, calculate errors on I(q),
601// just like in CircSectAve.ipf
602// TODO:
603// -- 2D Errors were NOT properly acculumated through reduction, so this loop of calculations is NOT MEANINGFUL (yet)
604// x- the error on the 1D intensity, is correctly calculated as the standard error of the mean.
605        for(ii=0;ii<nphi;ii+=1)
606       
607                phiBin_qxqy[ii] = dphi*ii
608               
609                if(nPhiBin_qxqy[ii] == 0)
610                        //no pixels in annuli, data unknown
611                        iPhiBin_qxqy[ii] = 0
612                        ePhiBin_qxqy[ii] = 1
613                        ePhiBin2D_qxqy[ii] = NaN
614                else
615                        if(nPhiBin_qxqy[ii] <= 1)
616                                //need more than one pixel to determine error
617                                iPhiBin_qxqy[ii] /= nPhiBin_qxqy[ii]
618                                ePhiBin_qxqy[ii] = 1
619                                ePhiBin2D_qxqy[ii] /= (nPhiBin_qxqy[ii])^2
620                        else
621                                //assume that the intensity in each pixel in annuli is normally distributed about mean...
622                                //  -- this is correctly calculating the error as the standard error of the mean, as
623                                //    was always done for SANS as well.
624                                iPhiBin_qxqy[ii] /= nPhiBin_qxqy[ii]
625                                avesq = iPhiBin_qxqy[ii]^2
626                                aveisq = iPhiBin2_qxqy[ii]/nPhiBin_qxqy[ii]
627                                var = aveisq-avesq
628                                if(var<=0)
629                                        ePhiBin_qxqy[ii] = 1e-6
630                                else
631                                        ePhiBin_qxqy[ii] = sqrt(var/(nPhiBin_qxqy[ii] - 1))
632                                endif
633                                // and calculate as it is propagated pixel-by-pixel
634                                ePhiBin2D_qxqy[ii] /= (nPhiBin_qxqy[ii])^2
635                        endif
636                endif
637        endfor
638       
639        ePhiBin2D_qxqy = sqrt(ePhiBin2D_qxqy)           // as equation (3) of John's memo
640       
641        // find the last non-zero point, working backwards
642        val=nq
643        do
644                val -= 1
645        while((nPhiBin_qxqy[val] == 0) && val > 0)
646       
647//      print val, nPhiBin_qxqy[val]
648        DeletePoints val, nq-val, iPhiBin_qxqy,phiBin_qxqy,nPhiBin_qxqy,iPhiBin2_qxqy,ePhiBin_qxqy,ePhiBin2D_qxqy
649
650        if(val == 0)
651                // all the points were deleted
652                return(0)
653        endif
654       
655       
656        // since the beam center is not always on the detector, many of the low Q bins will have zero pixels
657        // find the first non-zero point, working forwards
658        val = -1
659        do
660                val += 1
661        while(nPhiBin_qxqy[val] == 0)   
662        DeletePoints 0, val, iPhiBin_qxqy,phiBin_qxqy,nPhiBin_qxqy,iPhiBin2_qxqy,ePhiBin_qxqy,ePhiBin2D_qxqy
663
664        // ?? there still may be a point in the q-range that gets zero pixel contribution - so search this out and get rid of it
665        val = numpnts(nPhiBin_qxqy)-1
666        do
667                if(nPhiBin_qxqy[val] == 0)
668                        DeletePoints val, 1, iPhiBin_qxqy,phiBin_qxqy,nPhiBin_qxqy,iPhiBin2_qxqy,ePhiBin_qxqy,ePhiBin2D_qxqy
669                endif
670                val -= 1
671        while(val>0)
672       
673        // TODO:
674        // -- is this where I do the resolution calculation? This is where I calculate the resolution in SANS (see CircSectAve)
675        // -- or do I do it as a separate call?
676        // -- use the isVCALC flag to exclude VCALC from the resolution calculation if necessary
677        //
678       
679        SetDataFolder root:
680       
681        return(0)
682End
683
684
685Proc V_Write1DAnnular(pathStr,folderStr,detGroup,saveName)
686        String pathStr="root:Packages:NIST:VSANS:"
687        String folderStr="SAM"
688        String detGroup = "F"
689        String saveName = "Annular_Data.dat"
690       
691       
692        V_fWrite1DAnnular(pathStr,folderStr,detGroup,saveName)
693       
694end
695
696// TODO:
697// -- this is a temporary solution before a real writer is created
698// -- resolution is not generated here (and it shouldn't be) since resolution is not known yet.
699// -- but a real writer will need to be aware of resolution, and there may be different forms
700//
701// this will bypass save dialogs
702// -- AND WILL OVERWITE DATA WITH THE SAME NAME
703//
704//                      V_Write1DData_ITX("root:Packages:NIST:VSANS:",type,saveName,binType)
705//
706Function V_fWrite1DAnnular(pathStr,folderStr,detGroup,saveName)
707        String pathStr,folderStr,detGroup,saveName
708       
709        String formatStr="",fullpath=""
710        Variable refnum,dialog=1
711
712        SetDataFolder $(pathStr+folderStr)
713
714        if(cmpstr(detGroup,"F") == 0)
715                Wave/Z pw = phiBin_qxqy_FLRTB
716                Wave/Z iw = iPhiBin_qxqy_FLRTB
717                Wave/Z sw = ePhiBin_qxqy_FLRTB
718        endif
719        if(cmpstr(detGroup,"M") == 0)
720                Wave/Z pw = phiBin_qxqy_MLRTB
721                Wave/Z iw = iPhiBin_qxqy_MLRTB
722                Wave/Z sw = ePhiBin_qxqy_MLRTB
723        endif
724        if(cmpstr(detGroup,"B") == 0)
725                Wave/Z pw = phiBin_qxqy_B
726                Wave/Z iw = iPhiBin_qxqy_B
727                Wave/Z sw = ePhiBin_qxqy_B
728        endif
729       
730        String dataSetFolderParent,basestr
731       
732        SVAR gProtoStr = root:Packages:NIST:VSANS:Globals:Protocols:gProtoStr
733        Wave/T proto=$("root:Packages:NIST:VSANS:Globals:Protocols:"+gProtoStr)
734
735        SVAR samFiles = root:Packages:NIST:VSANS:Globals:Protocols:gSAM
736       
737        //make sure the waves exist
738       
739        if(WaveExists(pw) == 0)
740                SetDataFolder root:
741                Abort "phi is missing"
742        endif
743        if(WaveExists(iw) == 0)
744                SetDataFolder root:
745                Abort "i is missing"
746        endif
747        if(WaveExists(sw) == 0)
748                SetDataFolder root:
749                Abort "s is missing"
750        endif
751        if(WaveExists(proto) == 0)
752                SetDataFolder root:
753                Abort "protocol information is missing."
754        endif
755//      if(WaveExists(resw) == 0)
756//              SetDataFolder root:
757//              Abort "Resolution information is missing."
758//      endif
759       
760//      Duplicate/O qw qbar,sigQ,fs
761//      if(dimsize(resW,1) > 4)
762//              //it's USANS put -dQv back in the last 3 columns
763//              NVAR/Z dQv = USANS_dQv
764//              if(NVAR_Exists(dQv) == 0)
765//                      SetDataFolder root:
766//                      Abort "It's USANS data, and I don't know what the slit height is."
767//              endif
768//              sigQ = -dQv
769//              qbar = -dQv
770//              fs = -dQv
771//      else
772//              //it's SANS
773//              sigQ = resw[p][0]
774//              qbar = resw[p][1]
775//              fs = resw[p][2]
776//      endif
777//     
778
779
780// no beg/end trimming is used for annular data
781
782//// if the "default" trimming is used, the proto[] values will be null
783//// fill them in with the default values
784//      String protoStr7,protoStr8
785//      if(strlen(proto[7]) == 0)
786//              protoStr7 = "(Default) "+ ksBinTrimBegDefault
787//      else
788//              protoStr7 = proto[7]
789//      endif
790//      if(strlen(proto[8]) == 0)
791//              protoStr8 = "(Default) "+ ksBinTrimEndDefault
792//      else
793//              protoStr8 = proto[8]
794//      endif   
795
796
797
798        PathInfo catPathName
799        fullPath = S_Path + saveName
800
801        Open refnum as fullpath
802
803        fprintf refnum,"Annular data written from folder %s on %s\r\n",folderStr,(date()+" "+time())
804
805
806        //insert protocol information here
807        //-1 list of sample files
808        //0 - bkg
809        //1 - emp
810        //2 - div
811        //3 - mask
812        //4 - abs params c2-c5
813        //5 - average params
814        //6 - DRK (unused in VSANS)
815        //7 - beginning trim points
816        //8 - end trim points
817        fprintf refnum, "SAM: %s\r\n",samFiles
818        fprintf refnum, "BGD: %s\r\n",proto[0]
819        fprintf refnum, "EMP: %s\r\n",Proto[1]
820        fprintf refnum, "DIV: %s\r\n",Proto[2]
821        fprintf refnum, "MASK: %s\r\n",Proto[3]
822        fprintf refnum, "ABS Parameters (3-6): %s\r\n",Proto[4]
823        fprintf refnum, "Average Choices: %s\r\n",Proto[5]
824//      fprintf refnum, "Beginning Trim Points: %s\r\n",ProtoStr7
825//      fprintf refnum, "End Trim Points: %s\r\n",ProtoStr8
826
827
828// TODO -- make this work for 6-columns (or??)
829//      formatStr = "%15.4g %15.4g %15.4g %15.4g %15.4g %15.4g\r\n"     
830//      fprintf refnum, "The 6 columns are | Q (1/A) | I(Q) (1/cm) | std. dev. I(Q) (1/cm) | sigmaQ | meanQ | ShadowFactor|\r\n"       
831//      wfprintf refnum,formatStr,qw,iw,sw,sigQ,qbar,fs
832
833        //currently, only three columns
834        formatStr = "%15.4g %15.4g %15.4g\r\n" 
835        fprintf refnum, "The 3 columns are | Phi (degrees) | I(phi) (1/cm) | std. dev. I(phi) (1/cm)\r\n"       
836
837        wfprintf refnum,formatStr,pw,iw,sw
838        Close refnum
839       
840//      KillWaves/Z sigQ,qbar,fs
841        Print "Data written to: ",fullpath
842       
843        SetDataFolder root:
844        return(0)
845End
846
847
848//
849// x- I want to mask out everything that is "out" of the annulus
850//
851// 0 = keep the point
852// 1 = yes, mask the point
853Function V_MarkAnnularOverlayPixels(qTotal,overlay,qCtr_ann,qWidth)
854        Wave qTotal,overlay
855        Variable qCtr_ann,qWidth
856               
857       
858        Variable xDim=DimSize(qTotal, 0)
859        Variable yDim=DimSize(qTotal, 1)
860
861        Variable ii,jj,exclude,qVal
862       
863        // initialize the mask to == 1 == exclude everything
864        overlay = 1
865
866// now give every opportunity to keep pixel in
867        for(ii=0;ii<xDim;ii+=1)
868                for(jj=0;jj<yDim;jj+=1)
869                        //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
870                        qval = qTotal[ii][jj]
871                        exclude = 1
872               
873                        // annulus as defined
874                        if(V_CloseEnough(qval,qCtr_ann,qWidth))
875                                exclude = 0
876                        endif
877                       
878                        // set the mask value
879                        overlay[ii][jj] = exclude
880                endfor
881        endfor
882
883
884        return(0)
885End
Note: See TracBrowser for help on using the repository browser.