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

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

renamed white beam smearing models so they would be grouped together on the file list.

Re-worked the logic and flow of the averaging/plotting/saving steps of the reduction protocol so that it would flow cleanly and leave room for changes for the multitude of different collimation conditions. The averaging routines are now aware of the collimation conditions so that the appropriate resolution can be calculated. The collimation string is also written out to the averaged data file as element[9] of the protocol. The hope is that one could key on this collimation string to decide how to proceed with the analysis.

File size: 12.3 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 documentaiton 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                avesq = (sum_inten/sum_n)^2
202                aveisq = sum_inten2/sum_n
203                var = aveisq-avesq
204                if(var<=0)
205                        eBin_qxqy[ii] = 1e-6
206                else
207                        eBin_qxqy[ii] = sqrt(var/(sum_n - 1))
208                endif
209
210        endfor
211
212//TODO:  use only the Qx component in the y-center of the detector, not Qtotal
213// if the detectors are "L", then the values are all negative, so take the absolute value here
214        qBin_qxqy =  abs(qx[p][nYpix/2])               
215
216// for the L panels, sort the q-values (and data) after the abs() step, otherwise the data is reversed
217// won't hurt to sort the R panel data
218        Sort qBin_qxqy, qBin_qxqy,iBin_qxqy,eBin_qxqy
219
220
221        // TODO: do I use dQ for the height of the panel?
222        // TODO : do I use 1/2 of dQ due to the symmetry of my smearing calculation?   
223        delQy = abs(qy[0][nYpix-1] - qy[0][0])
224
225//      iBin_qxqy *= delQy
226//      eBin_qxqy *= delQy
227       
228
229/// TODO -- this is not necessary, but just for getting the I(Q) display to look "pretty"
230// clean out the near-zero Q point in the T/B  and Back detectors
231//      qBin_qxqy = (abs(qBin_qxqy[p][q]) < 1e-5) ? NaN : qBin_qxqy[p][q]                       
232
233
234// clear out zero data values before exiting...
235        // find the last non-zero point, working backwards
236        val = numpnts(iBin_qxqy)
237        do
238                val -= 1
239        while((iBin_qxqy[val] == 0) && val > 0)
240       
241        DeletePoints val, nq-val, iBin_qxqy,qBin_qxqy,eBin_qxqy
242
243// work forwards? this doesn't work...
244//      val = -1
245//      do
246//              val += 1
247//      while(nBin_qxqy[val] == 0 && val < numpnts(iBin_qxqy)-1)       
248//      DeletePoints 0, val, iBin_qxqy,qBin_qxqy,eBin_qxqy
249
250
251        // TODO:
252        // -- This is where I calculate the resolution in SANS (see CircSectAve)
253        // -- use the isVCALC flag to exclude VCALC from the resolution calculation if necessary
254        // -- from the top of the function, folderStr = work folder, detStr = "FLRTB" or other type of averaging
255        //
256
257        nq = numpnts(qBin_qxqy)
258        Make/O/D/N=(nq)  $(folderPath+":"+"sigmaQ_qxqy"+"_"+detStr)
259        Make/O/D/N=(nq)  $(folderPath+":"+"qBar_qxqy"+"_"+detStr)
260        Make/O/D/N=(nq)  $(folderPath+":"+"fSubS_qxqy"+"_"+detStr)
261        Wave sigmaq = $(folderPath+":"+"sigmaQ_qxqy_"+detStr)
262        Wave qbar = $(folderPath+":"+"qBar_qxqy_"+detStr)
263        Wave fsubs = $(folderPath+":"+"fSubS_qxqy_"+detStr)
264
265// TODO:
266// -- calculate the slit resolution here. don't really have any idea how to represent this in VSANS
267//    since it's not an infinite slit, and not anything that I expect could be represented as a Gaussian
268// -- there is also wavelength smearing present
269
270
271
272// ASSUMPTION: As a first approximation, ignore the wavelength smearing component       
273        // TODO : do I use 1/2 of dQy due to the symmetry of my smearing calculation?   (as infinite slit)
274        sigmaq = -delQy/2
275        qbar = -delQy/2
276        fsubs = -delQy/2
277
278
279        SetDataFolder root:
280       
281        return(0)
282End
283
284
285// TODO -- update to new folder structure
286// unused -- update if necessary
287Proc CopyIQWaves()
288
289        SetDataFolder root:Packages:NIST:VSANS:VCALC   
290
291        duplicate/O iBin_qxqy_B iBin_qxqy_B_pin
292        duplicate/O qBin_qxqy_B qBin_qxqy_B_pin
293       
294        duplicate/O iBin_qxqy_MR iBin_qxqy_MR_pin
295        duplicate/O qBin_qxqy_MR qBin_qxqy_MR_pin
296        duplicate/O iBin_qxqy_MT iBin_qxqy_MT_pin
297        duplicate/O qBin_qxqy_MT qBin_qxqy_MT_pin
298        duplicate/O iBin_qxqy_ML iBin_qxqy_ML_pin
299        duplicate/O qBin_qxqy_ML qBin_qxqy_ML_pin
300        duplicate/O iBin_qxqy_MB iBin_qxqy_MB_pin
301        duplicate/O qBin_qxqy_MB qBin_qxqy_MB_pin
302
303        duplicate/O iBin_qxqy_FR iBin_qxqy_FR_pin
304        duplicate/O qBin_qxqy_FR qBin_qxqy_FR_pin
305        duplicate/O iBin_qxqy_FT iBin_qxqy_FT_pin
306        duplicate/O qBin_qxqy_FT qBin_qxqy_FT_pin
307        duplicate/O iBin_qxqy_FL iBin_qxqy_FL_pin
308        duplicate/O qBin_qxqy_FL qBin_qxqy_FL_pin
309        duplicate/O iBin_qxqy_FB iBin_qxqy_FB_pin
310        duplicate/O qBin_qxqy_FB qBin_qxqy_FB_pin
311
312        SetDataFolder root:
313End
314
315// TODO -- update to new folder structure
316// unused -- update if necessary
317Window slit_vs_pin_graph() : Graph
318        PauseUpdate; Silent 1           // building window...
319        String fldrSav0= GetDataFolder(1)
320        SetDataFolder root:Packages:NIST:VSANS:VCALC:
321        Display /W=(1296,44,1976,696) iBin_qxqy_B vs qBin_qxqy_B
322        AppendToGraph iBin_qxqy_B_pin vs qBin_qxqy_B_pin
323        AppendToGraph iBin_qxqy_MR vs qBin_qxqy_MR
324        AppendToGraph iBin_qxqy_MR_pin vs qBin_qxqy_MR_pin
325        AppendToGraph iBin_qxqy_MT vs qBin_qxqy_MT
326        AppendToGraph iBin_qxqy_MT_pin vs qBin_qxqy_MT_pin
327        AppendToGraph iBin_qxqy_FR vs qBin_qxqy_FR
328        AppendToGraph iBin_qxqy_FR_pin vs qBin_qxqy_FR_pin
329        AppendToGraph iBin_qxqy_FT vs qBin_qxqy_FT
330        AppendToGraph iBin_qxqy_FT_pin vs qBin_qxqy_FT_pin
331        SetDataFolder fldrSav0
332        ModifyGraph mode=4
333        ModifyGraph marker=19
334        ModifyGraph lSize=2
335        ModifyGraph rgb(iBin_qxqy_B)=(0,0,0),rgb(iBin_qxqy_B_pin)=(65535,16385,16385),rgb(iBin_qxqy_MR)=(2,39321,1)
336        ModifyGraph rgb(iBin_qxqy_MR_pin)=(0,0,65535),rgb(iBin_qxqy_MT)=(39321,1,31457)
337        ModifyGraph rgb(iBin_qxqy_MT_pin)=(48059,48059,48059),rgb(iBin_qxqy_FR)=(65535,32768,32768)
338        ModifyGraph rgb(iBin_qxqy_FR_pin)=(0,65535,0),rgb(iBin_qxqy_FT)=(16385,65535,65535)
339        ModifyGraph rgb(iBin_qxqy_FT_pin)=(65535,32768,58981)
340        ModifyGraph msize=2
341        ModifyGraph grid=1
342        ModifyGraph log=1
343        ModifyGraph mirror=1
344        SetAxis bottom 1e-05,0.05287132
345        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"
346        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"
347        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"
348EndMacro
Note: See TracBrowser for help on using the repository browser.