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

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

lots of changes here:
many little fixes to clean up TODO items and marke them DONE

changed the handling of the panel "gap" to split the gap evenly. Q-calculations have been re-verified with this change.

re-named the list of "bin Type" values, and added a few more choices. Streamlined how the averaging and plotting works with this list so that it can be more easily modified as different combinations of binning are envisioned. This resulted in a lot of excess code being cut out and replaced with cleaner logic. This change has also been verified to work as intended.

Attenuation is now always calculated from the table. The table also by (NEW) definition has values for the white beam (one waelength) and graphite (multiple possible wavelengths) where the wavelengths are artificially scaled (*1000) or *1e6) so that the interpolations can be done internally without the need for multiple attenuator tables.

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