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

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

corrected the behavior of the annular averaging and added a button to the main panel to allow quick access.

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