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

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

more attempts to correct the slit binning - to get the absolute scaling correct - still in progress.

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