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

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

adding the annular average option to the protocol definition and execution.

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