source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/VC_DetectorBinning_Slit.ipf @ 1100

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

Significant restructuring of V_ExecuteProtocol to make the logic of the flow more tractable and adaptable in the future.

Added major change to VSANS event mode reduction panel to allow F and M binned events to be saved to a copy of the correspondng data file. A new data reduction panel for multiple slice reduction is now presented for Event data. This allows reduction of a single slice (all 8 F,M panels, back ignored), or all of the slices.

File size: 12.4 KB
Line 
1#pragma rtGlobals=3             // Use modern global access method and strict wave access.
2
3
4//////////////////
5//
6// Procedures to average data taken in "slit" aperture geometry. The fake data on the detector panels
7//  is as would be collected in PINHOLE geometry - I do not currently have a simulation for slit
8//  apertures (this would need to be MonteCarlo) - the I(q) averaging here gives the I(q) that you
9//  would measure in 1D using slit geometry. It is done by simply summing the columns of data on each detector.
10//
11//////////////////
12//
13//
14// TODO - big question about averaging in this way...
15// can the T/B panels really be used at all for slit mode? - since there's a big "hole" in the scattering data
16// collected -- you're not getting the full column of data covering a wide range of Qy. L/R panels should be fine,
17//  although the different carriages cover very different "chunks" of Qy
18//
19// - best answer so far is to skip the T/B panels, and simply not use them
20//
21// TODO -- be sure that the absolute scaling of this is correct. Right now, something is off.
22//  -- so write up the documentation for this operation, and try to find the error in the process...
23//
24//
25/////////////////
26
27
28// TODO:
29// -- verify the error calculation
30// -- ? add in functionality to handle FLR and MLR cases (2 panels of data)
31//
32// seems backwards to call this "byRows", but this is the way that Igor indexes
33// LR banks are defined as (48,256) (n,m), sumRows gives sum w/ dimension (n x 1)
34//
35// updated to new folder structure Feb 2016
36// folderStr = RAW,SAM, VCALC or other
37// detStr is the panel identifer "ML", etc.
38Function VC_fBinDetector_byRows(folderStr,detStr)
39        String folderStr,detStr
40       
41//      SetDataFolder root:Packages:NIST:VSANS:VCALC   
42       
43        Variable pixSizeX,pixSizeY,delQy
44        Variable isVCALC=0,maskMissing,nsets=0
45       
46        maskMissing = 1         // set to zero if a mask is actually present
47       
48        if(cmpstr(folderStr,"VCALC") == 0)
49                isVCALC = 1
50        endif
51
52        String folderPath = "root:Packages:NIST:VSANS:"+folderStr
53        String instPath = ":entry:instrument:detector_"
54
55        strswitch(detStr)       // string switch       
56// only one panel, simply pick that panel and move on out of the switch
57                case "FL":
58                case "FR":
59                case "ML":
60                case "MR":
61                case "B":
62                        if(isVCALC)
63                                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)               // 2D detector data
64                                WAVE/Z iErr = $("asdf_iErr_"+detStr)                    // TODO: 2D errors -- may not exist, especially for simulation
65                        else
66                                Wave inten = V_getDetectorDataW(folderStr,detStr)
67                                Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
68                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
69                                if(WaveExists(mask) == 1)
70                                        maskMissing = 0
71                                endif
72                        endif
73                       
74                        nsets = 1
75                        break
76
77//              case "FLR":
78//              // detStr has multiple values now, so unfortuntely, I'm hard-wiring things...
79//                      if(isVCALC)
80//                              WAVE inten = $(folderPath+instPath+"FL"+":det_"+"FL")
81//                              WAVE/Z iErr = $("iErr_"+"FL")                   // 2D errors -- may not exist, especially for simulation               
82//                              WAVE inten2 = $(folderPath+instPath+"FR"+":det_"+"FR")
83//                              WAVE/Z iErr2 = $("iErr_"+"FR")                  // 2D errors -- may not exist, especially for simulation       
84//                      else
85//                              Wave inten = V_getDetectorDataW(folderStr,"FL")
86//                              Wave iErr = V_getDetectorDataErrW(folderStr,"FL")
87//                              Wave inten2 = V_getDetectorDataW(folderStr,"FR")
88//                              Wave iErr2 = V_getDetectorDataErrW(folderStr,"FR")
89//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FL"+":data")
90//                              Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FR"+":data")
91//                              if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
92//                                      maskMissing = 0
93//                              endif
94//                      endif   
95//
96//                      nsets = 2
97//                      break
98//                     
99//              case "MLR":
100//                      if(isVCALC)
101//                              WAVE inten = $(folderPath+instPath+"ML"+":det_"+"ML")
102//                              WAVE/Z iErr = $("iErr_"+"ML")                   // 2D errors -- may not exist, especially for simulation               
103//                              WAVE inten2 = $(folderPath+instPath+"MR"+":det_"+"MR")
104//                              WAVE/Z iErr2 = $("iErr_"+"MR")                  // 2D errors -- may not exist, especially for simulation       
105//                      else
106//                              Wave inten = V_getDetectorDataW(folderStr,"ML")
107//                              Wave iErr = V_getDetectorDataErrW(folderStr,"ML")
108//                              Wave inten2 = V_getDetectorDataW(folderStr,"MR")
109//                              Wave iErr2 = V_getDetectorDataErrW(folderStr,"MR")
110//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"ML"+":data")
111//                              Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MR"+":data")
112//                              if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
113//                                      maskMissing = 0
114//                              endif
115//                      endif   
116//             
117//                      nSets = 2
118//                      break                   
119
120                default:
121                        nSets = 0                                                       
122                        Print "ERROR   ---- type is not recognized "
123        endswitch
124
125        if(nSets == 0)
126                SetDataFolder root:
127                return(0)
128        endif
129
130
131        Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
132        Wave qx = $(folderPath+instPath+detStr+":qx_"+detStr)
133        Wave qy = $(folderPath+instPath+detStr+":qy_"+detStr)
134
135
136// delta Qx is set by the pixel X dimension of the detector, which is the limiting resolution
137//      delQx = abs(qx[0][0] - qx[1][0])
138       
139        Variable nq,val
140        nq = DimSize(inten,0)           //nq == the number of columns (x dimension)
141       
142//      SetDataFolder $(folderPath+instPath+detStr)     
143
144        Make/O/D/N=(nq)  $(folderPath+":"+"iBin_qxqy_"+detStr)
145        Make/O/D/N=(nq)  $(folderPath+":"+"qBin_qxqy_"+detStr)
146        Make/O/D/N=(nq)  $(folderPath+":"+"nBin_qxqy_"+detStr)
147        Make/O/D/N=(nq)  $(folderPath+":"+"iBin2_qxqy_"+detStr)
148        Make/O/D/N=(nq)  $(folderPath+":"+"eBin_qxqy_"+detStr)
149        Make/O/D/N=(nq)  $(folderPath+":"+"eBin2D_qxqy_"+detStr)
150       
151        Wave iBin_qxqy = $(folderPath+":"+"iBin_qxqy_"+detStr)
152        Wave qBin_qxqy = $(folderPath+":"+"qBin_qxqy_"+detStr)
153        Wave nBin_qxqy = $(folderPath+":"+"nBin_qxqy_"+detStr)
154        Wave iBin2_qxqy = $(folderPath+":"+"iBin2_qxqy_"+detStr)
155        Wave eBin_qxqy = $(folderPath+":"+"eBin_qxqy_"+detStr)
156        Wave eBin2D_qxqy = $(folderPath+":"+"eBin2D_qxqy_"+detStr)
157
158
159        iBin_qxqy = 0
160        iBin2_qxqy = 0
161        eBin_qxqy = 0
162        eBin2D_qxqy = 0
163        nBin_qxqy = 0   
164
165
166// sum the rows
167
168// MatrixOp would be fast, but I can't figure out how to apply the mask with sumRows???
169//      MatrixOp/O iBin_qxqy = sumRows(inten)   //automatically generates the destination
170//             
171//// how to properly calculate the error?
172//// This MatrixOp gives values way too large -- larger than the intensity (be sure to correct *delQy below...
173//
174//      MatrixOp/O tmp = sqrt(varCols(inten^t))         // variance: no varRows operation, so use the transpose of the matrix
175//      eBin_qxqy = tmp[0][p]
176
177        Variable ii,jj,ntube,nYpix,sum_inten, sum_n, sum_inten2,avesq,aveisq,var,mask_val
178        ntube = DimSize(inten,0)
179        nYpix = DimSize(inten,1)
180       
181        for(ii=0;ii<ntube;ii+=1)
182                sum_inten = 0                   // initialize the sum
183                sum_n = 0
184                sum_inten2 = 0
185               
186                for(jj=0;jj<nYpix;jj+=1)
187                                val = inten[ii][jj]
188                                if(isVCALC || maskMissing)              // mask_val == 0 == keep, mask_val == 1 = YES, mask out the point
189                                        mask_val = 0
190                                else
191                                        mask_val = mask[ii][jj]
192                                endif
193                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
194                                        sum_inten += val
195                                        sum_n += 1
196                                        sum_inten2 += val*val
197                                endif
198                endfor
199                iBin_qxqy[ii] = sum_inten/sum_n         //the average value
200               
201//              if(numtype(iBin_qxqy[ii]) == 2)
202//                      print "asdfasdf"
203//              endif
204               
205                avesq = (sum_inten/sum_n)^2
206                aveisq = sum_inten2/sum_n
207                var = aveisq-avesq
208                if(var<=0)
209                        eBin_qxqy[ii] = 1e-6
210                else
211                        eBin_qxqy[ii] = sqrt(var/(sum_n - 1))
212                endif
213
214        endfor
215
216//TODO:  use only the Qx component in the y-center of the detector, not Qtotal
217// if the detectors are "L", then the values are all negative, so take the absolute value here
218        qBin_qxqy =  abs(qx[p][nYpix/2])               
219
220// for the L panels, sort the q-values (and data) after the abs() step, otherwise the data is reversed
221// won't hurt to sort the R panel data
222        Sort qBin_qxqy, qBin_qxqy,iBin_qxqy,eBin_qxqy
223
224
225        // TODO: do I use dQ for the height of the panel?
226        // TODO : do I use 1/2 of dQ due to the symmetry of my smearing calculation?   
227        delQy = abs(qy[0][nYpix-1] - qy[0][0])
228
229//      iBin_qxqy *= delQy
230//      eBin_qxqy *= delQy
231       
232
233/// TODO -- this is not necessary, but just for getting the I(Q) display to look "pretty"
234// clean out the near-zero Q point in the T/B  and Back detectors
235//      qBin_qxqy = (abs(qBin_qxqy[p][q]) < 1e-5) ? NaN : qBin_qxqy[p][q]                       
236
237
238// clear out zero data values before exiting...
239        // find the last non-zero point, working backwards
240        val = numpnts(iBin_qxqy)
241        do
242                val -= 1
243        while((iBin_qxqy[val] == 0) && val > 0)
244       
245        DeletePoints val, nq-val, iBin_qxqy,qBin_qxqy,eBin_qxqy
246
247
248// utility function to remove NaN values from the waves
249        V_RemoveNaNsQIS(qBin_qxqy, iBin_qxqy, eBin_qxqy)
250
251// work forwards? this doesn't work...
252//      val = -1
253//      do
254//              val += 1
255//      while(nBin_qxqy[val] == 0 && val < numpnts(iBin_qxqy)-1)       
256//      DeletePoints 0, val, iBin_qxqy,qBin_qxqy,eBin_qxqy
257
258
259        // TODO:
260        // -- This is where I calculate the resolution in SANS (see CircSectAve)
261        // -- use the isVCALC flag to exclude VCALC from the resolution calculation if necessary
262        // -- from the top of the function, folderStr = work folder, detStr = "FLRTB" or other type of averaging
263        //
264
265        nq = numpnts(qBin_qxqy)
266        Make/O/D/N=(nq)  $(folderPath+":"+"sigmaQ_qxqy"+"_"+detStr)
267        Make/O/D/N=(nq)  $(folderPath+":"+"qBar_qxqy"+"_"+detStr)
268        Make/O/D/N=(nq)  $(folderPath+":"+"fSubS_qxqy"+"_"+detStr)
269        Wave sigmaq = $(folderPath+":"+"sigmaQ_qxqy_"+detStr)
270        Wave qbar = $(folderPath+":"+"qBar_qxqy_"+detStr)
271        Wave fsubs = $(folderPath+":"+"fSubS_qxqy_"+detStr)
272
273// TODO:
274// -- calculate the slit resolution here. don't really have any idea how to represent this in VSANS
275//    since it's not an infinite slit, and not anything that I expect could be represented as a Gaussian
276// -- there is also wavelength smearing present
277
278
279
280// ASSUMPTION: As a first approximation, ignore the wavelength smearing component       
281        // TODO : do I use 1/2 of dQy due to the symmetry of my smearing calculation?   (as infinite slit)
282        sigmaq = -delQy/2
283        qbar = -delQy/2
284        fsubs = -delQy/2
285
286
287        SetDataFolder root:
288       
289        return(0)
290End
291
292
293// TODO -- update to new folder structure
294// unused -- update if necessary
295Proc CopyIQWaves()
296
297        SetDataFolder root:Packages:NIST:VSANS:VCALC   
298
299        duplicate/O iBin_qxqy_B iBin_qxqy_B_pin
300        duplicate/O qBin_qxqy_B qBin_qxqy_B_pin
301       
302        duplicate/O iBin_qxqy_MR iBin_qxqy_MR_pin
303        duplicate/O qBin_qxqy_MR qBin_qxqy_MR_pin
304        duplicate/O iBin_qxqy_MT iBin_qxqy_MT_pin
305        duplicate/O qBin_qxqy_MT qBin_qxqy_MT_pin
306        duplicate/O iBin_qxqy_ML iBin_qxqy_ML_pin
307        duplicate/O qBin_qxqy_ML qBin_qxqy_ML_pin
308        duplicate/O iBin_qxqy_MB iBin_qxqy_MB_pin
309        duplicate/O qBin_qxqy_MB qBin_qxqy_MB_pin
310
311        duplicate/O iBin_qxqy_FR iBin_qxqy_FR_pin
312        duplicate/O qBin_qxqy_FR qBin_qxqy_FR_pin
313        duplicate/O iBin_qxqy_FT iBin_qxqy_FT_pin
314        duplicate/O qBin_qxqy_FT qBin_qxqy_FT_pin
315        duplicate/O iBin_qxqy_FL iBin_qxqy_FL_pin
316        duplicate/O qBin_qxqy_FL qBin_qxqy_FL_pin
317        duplicate/O iBin_qxqy_FB iBin_qxqy_FB_pin
318        duplicate/O qBin_qxqy_FB qBin_qxqy_FB_pin
319
320        SetDataFolder root:
321End
322
323// TODO -- update to new folder structure
324// unused -- update if necessary
325Window slit_vs_pin_graph() : Graph
326        PauseUpdate; Silent 1           // building window...
327        String fldrSav0= GetDataFolder(1)
328        SetDataFolder root:Packages:NIST:VSANS:VCALC:
329        Display /W=(1296,44,1976,696) iBin_qxqy_B vs qBin_qxqy_B
330        AppendToGraph iBin_qxqy_B_pin vs qBin_qxqy_B_pin
331        AppendToGraph iBin_qxqy_MR vs qBin_qxqy_MR
332        AppendToGraph iBin_qxqy_MR_pin vs qBin_qxqy_MR_pin
333        AppendToGraph iBin_qxqy_MT vs qBin_qxqy_MT
334        AppendToGraph iBin_qxqy_MT_pin vs qBin_qxqy_MT_pin
335        AppendToGraph iBin_qxqy_FR vs qBin_qxqy_FR
336        AppendToGraph iBin_qxqy_FR_pin vs qBin_qxqy_FR_pin
337        AppendToGraph iBin_qxqy_FT vs qBin_qxqy_FT
338        AppendToGraph iBin_qxqy_FT_pin vs qBin_qxqy_FT_pin
339        SetDataFolder fldrSav0
340        ModifyGraph mode=4
341        ModifyGraph marker=19
342        ModifyGraph lSize=2
343        ModifyGraph rgb(iBin_qxqy_B)=(0,0,0),rgb(iBin_qxqy_B_pin)=(65535,16385,16385),rgb(iBin_qxqy_MR)=(2,39321,1)
344        ModifyGraph rgb(iBin_qxqy_MR_pin)=(0,0,65535),rgb(iBin_qxqy_MT)=(39321,1,31457)
345        ModifyGraph rgb(iBin_qxqy_MT_pin)=(48059,48059,48059),rgb(iBin_qxqy_FR)=(65535,32768,32768)
346        ModifyGraph rgb(iBin_qxqy_FR_pin)=(0,65535,0),rgb(iBin_qxqy_FT)=(16385,65535,65535)
347        ModifyGraph rgb(iBin_qxqy_FT_pin)=(65535,32768,58981)
348        ModifyGraph msize=2
349        ModifyGraph grid=1
350        ModifyGraph log=1
351        ModifyGraph mirror=1
352        SetAxis bottom 1e-05,0.05287132
353        Legend/C/N=text0/J/X=64.73/Y=7.25 "\\Z12\\s(iBin_qxqy_B) iBin_qxqy_B\r\\s(iBin_qxqy_B_pin) iBin_qxqy_B_pin\r\\s(iBin_qxqy_MR) iBin_qxqy_MR"
354        AppendText "\\s(iBin_qxqy_MR_pin) iBin_qxqy_MR_pin\r\\s(iBin_qxqy_MT) iBin_qxqy_MT\r\\s(iBin_qxqy_MT_pin) iBin_qxqy_MT_pin\r\\s(iBin_qxqy_FR) iBin_qxqy_FR"
355        AppendText "\\s(iBin_qxqy_FR_pin) iBin_qxqy_FR_pin\r\\s(iBin_qxqy_FT) iBin_qxqy_FT\r\\s(iBin_qxqy_FT_pin) iBin_qxqy_FT_pin"
356EndMacro
Note: See TracBrowser for help on using the repository browser.