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

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

Added two model functions for white beam smearing.

Many other small changes for processing of the back detector, shuffling of VSANS menu items, and consistent naming of V_ procedures.

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