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

Last change on this file since 1097 was 1097, checked in by srkline, 4 years ago

added procedures to compare files to see if they are from the same configuration, same wavelength, etc. so they can be properly chosen for transmission files, scattering files, and properly identified for the different resolution conditions.

Re-worked the logic of dispatching averaging, plotting, and saving in Execute_Protocol. Hopefully this will alow for easier dispatching for future conditions, including getting the correct resolution calculation.

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