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

Last change on this file since 1242 was 1242, checked in by srkline, 2 years ago

updating the IgorVersion? pragma to v7.0 for all files to be consistent.

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