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

Last change on this file since 1239 was 1239, checked in by srkline, 3 years ago

adding downstream window transmission correction -- added function to detector corrections, and added item to preference panel. Transmission value is currently set as a global since the value is not part of the VSANS header. Global value defaults to T=1= no correction.

Other changes are cleanup of TODO items that were already done and tested. inline commnents have been updated.

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