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

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

A large number of changes to the size of panels to enable "Laptop Mode" where all of the panels and controls are scaled down so that they fit on screen and are still in correct proportion. For the laptop I'm using for testing, the resolution is 1920x1080. For this, a scaling of 0.7 seems to work. The on/off of the "laptop Mode" is controlled by a checkbox in the preference panel (under the General tab).

There are still more panels to update in the next commit.

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