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

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

woking on routines to correctly reduce data colloected with slit apertures. trying to get the data to absolute scale, with proper errro bars and resolution information.

Still not there yet.

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" apertuere 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//
18// - best answer so far is to skip the T/B panels, and simply not use them
19//
20// TODO -- be sure that the absolute scaling of this is correct. I have no guarantee that
21//                      it'll be correct.
22//
23/////////////////
24
25
26// TODO:
27// -- verify the error calculation
28// -- add in functionality to handle FLR and MLR cases (2 panels of data)
29//
30// seems backwards to call this "byRows", but this is the way that Igor indexes
31// LR banks are defined as (48,256) (n,m), sumRows gives sum w/ dimension (n x 1)
32//
33// updated to new folder structure Feb 2016
34// folderStr = RAW,SAM, VCALC or other
35// detStr is the panel identifer "ML", etc.
36Function VC_fBinDetector_byRows(folderStr,detStr)
37        String folderStr,detStr
38       
39//      SetDataFolder root:Packages:NIST:VSANS:VCALC   
40       
41        Variable pixSizeX,pixSizeY,delQy
42        Variable isVCALC=0,maskMissing,nsets=0
43       
44        maskMissing = 1         // set to zero if a mask is actually present
45       
46        if(cmpstr(folderStr,"VCALC") == 0)
47                isVCALC = 1
48        endif
49
50        String folderPath = "root:Packages:NIST:VSANS:"+folderStr
51        String instPath = ":entry:instrument:detector_"
52
53        strswitch(detStr)       // string switch       
54// only one panel, simply pick that panel and move on out of the switch
55                case "FL":
56                case "FR":
57                case "ML":
58                case "MR":
59                case "B":
60                        if(isVCALC)
61                                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)               // 2D detector data
62                                WAVE/Z iErr = $("asdf_iErr_"+detStr)                    // TODO: 2D errors -- may not exist, especially for simulation
63                        else
64                                Wave inten = V_getDetectorDataW(folderStr,detStr)
65                                Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
66                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
67                                if(WaveExists(mask) == 1)
68                                        maskMissing = 0
69                                endif
70                        endif
71                       
72                        nsets = 1
73                        break
74
75//              case "FLR":
76//              // detStr has multiple values now, so unfortuntely, I'm hard-wiring things...
77//                      if(isVCALC)
78//                              WAVE inten = $(folderPath+instPath+"FL"+":det_"+"FL")
79//                              WAVE/Z iErr = $("iErr_"+"FL")                   // 2D errors -- may not exist, especially for simulation               
80//                              WAVE inten2 = $(folderPath+instPath+"FR"+":det_"+"FR")
81//                              WAVE/Z iErr2 = $("iErr_"+"FR")                  // 2D errors -- may not exist, especially for simulation       
82//                      else
83//                              Wave inten = V_getDetectorDataW(folderStr,"FL")
84//                              Wave iErr = V_getDetectorDataErrW(folderStr,"FL")
85//                              Wave inten2 = V_getDetectorDataW(folderStr,"FR")
86//                              Wave iErr2 = V_getDetectorDataErrW(folderStr,"FR")
87//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FL"+":data")
88//                              Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FR"+":data")
89//                              if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
90//                                      maskMissing = 0
91//                              endif
92//                      endif   
93//
94//                      nsets = 2
95//                      break
96//                     
97//              case "MLR":
98//                      if(isVCALC)
99//                              WAVE inten = $(folderPath+instPath+"ML"+":det_"+"ML")
100//                              WAVE/Z iErr = $("iErr_"+"ML")                   // 2D errors -- may not exist, especially for simulation               
101//                              WAVE inten2 = $(folderPath+instPath+"MR"+":det_"+"MR")
102//                              WAVE/Z iErr2 = $("iErr_"+"MR")                  // 2D errors -- may not exist, especially for simulation       
103//                      else
104//                              Wave inten = V_getDetectorDataW(folderStr,"ML")
105//                              Wave iErr = V_getDetectorDataErrW(folderStr,"ML")
106//                              Wave inten2 = V_getDetectorDataW(folderStr,"MR")
107//                              Wave iErr2 = V_getDetectorDataErrW(folderStr,"MR")
108//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"ML"+":data")
109//                              Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MR"+":data")
110//                              if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
111//                                      maskMissing = 0
112//                              endif
113//                      endif   
114//             
115//                      nSets = 2
116//                      break                   
117
118                default:
119                        nSets = 0                                                       
120                        Print "ERROR   ---- type is not recognized "
121        endswitch
122
123        if(nSets == 0)
124                SetDataFolder root:
125                return(0)
126        endif
127
128
129        Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
130        Wave qx = $(folderPath+instPath+detStr+":qx_"+detStr)
131        Wave qy = $(folderPath+instPath+detStr+":qy_"+detStr)
132
133
134// delta Qx is set by the pixel X dimension of the detector, which is the limiting resolution
135//      delQx = abs(qx[0][0] - qx[1][0])
136       
137        Variable nq,val
138        nq = DimSize(inten,0)           //nq == the number of columns (x dimension)
139       
140//      SetDataFolder $(folderPath+instPath+detStr)     
141
142        Make/O/D/N=(nq)  $(folderPath+":"+"iBin_qxqy_"+detStr)
143        Make/O/D/N=(nq)  $(folderPath+":"+"qBin_qxqy_"+detStr)
144        Make/O/D/N=(nq)  $(folderPath+":"+"nBin_qxqy_"+detStr)
145        Make/O/D/N=(nq)  $(folderPath+":"+"iBin2_qxqy_"+detStr)
146        Make/O/D/N=(nq)  $(folderPath+":"+"eBin_qxqy_"+detStr)
147        Make/O/D/N=(nq)  $(folderPath+":"+"eBin2D_qxqy_"+detStr)
148       
149        Wave iBin_qxqy = $(folderPath+":"+"iBin_qxqy_"+detStr)
150        Wave qBin_qxqy = $(folderPath+":"+"qBin_qxqy_"+detStr)
151        Wave nBin_qxqy = $(folderPath+":"+"nBin_qxqy_"+detStr)
152        Wave iBin2_qxqy = $(folderPath+":"+"iBin2_qxqy_"+detStr)
153        Wave eBin_qxqy = $(folderPath+":"+"eBin_qxqy_"+detStr)
154        Wave eBin2D_qxqy = $(folderPath+":"+"eBin2D_qxqy_"+detStr)
155
156
157        iBin_qxqy = 0
158        iBin2_qxqy = 0
159        eBin_qxqy = 0
160        eBin2D_qxqy = 0
161        nBin_qxqy = 0   
162
163
164// sum the rows
165
166// MatrixOp would be fast, but I can't figure out how to apply the mask with sumRows???
167//      MatrixOp/O iBin_qxqy = sumRows(inten)   //automatically generates the destination
168//             
169//// how to properly calculate the error?
170//// This MatrixOp gives values way too large -- larger than the intensity (be sure to correct *delQy below...
171//
172//      MatrixOp/O tmp = sqrt(varCols(inten^t))         // variance: no varRows operation, so use the transpose of the matrix
173//      eBin_qxqy = tmp[0][p]
174
175        Variable ii,jj,ntube,npix,sum_inten, sum_n, sum_inten2,avesq,aveisq,var,mask_val
176        ntube = DimSize(inten,0)
177        npix = DimSize(inten,1)
178       
179        for(ii=0;ii<ntube;ii+=1)
180                sum_inten = 0                   // initialize the sum
181                sum_n = 0
182                sum_inten2 = 0
183               
184                for(jj=0;jj<npix;jj+=1)
185                                val = inten[ii][jj]
186                                if(isVCALC || maskMissing)              // mask_val == 0 == keep, mask_val == 1 = YES, mask out the point
187                                        mask_val = 0
188                                else
189                                        mask_val = mask[ii][jj]
190                                endif
191                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
192                                        sum_inten += val
193                                        sum_n += 1
194                                        sum_inten2 += val*val
195                                endif
196                endfor
197                iBin_qxqy[ii] = sum_inten
198               
199                avesq = sum_inten
200                aveisq = sum_inten2/sum_n
201                var = aveisq-avesq
202                if(var<=0)
203                        eBin_qxqy[ii] = 1e-6
204                else
205                        eBin_qxqy[ii] = sqrt(var/(sum_n - 1))
206                endif
207
208        endfor
209
210        qBin_qxqy =  qx[p][npix/2]              //TODO:  use only the Qx component in the y-center of the detector, not Qtotal
211               
212// if the detectors are "L", then the values are all negative...
213        qBin_qxqy = abs(qx[p][0])
214
215// for the L panels, sort the q-values (and data) after the abs() step, otherwise the data is reversed
216// won't hurt to sort the R data
217        Sort qBin_qxqy, qBin_qxqy,iBin_qxqy,eBin_qxqy
218       
219        //now get the scaling correct
220        // q-integration (rectangular), matrixOp simply summed, so I need to multiply by dy (pixelSizeY -> as Qy?)
221
222//      delQy = abs(qy[0][1] - qy[0][0])        // this is only one pixel
223        delQy = abs(qy[0][npix-1] - qy[0][0])   // TODO: do I use dQ for the height of the panel?
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// -- calculate the slit resolution here. don't really have any idea how to represent this in VSANS
253//    since it's not an infinite slit, and not anything that I expect could be represented as a Gaussian
254
255        // TODO:
256        // -- This is where I calculate the resolution in SANS (see CircSectAve)
257        // -- use the isVCALC flag to exclude VCALC from the resolution calculation if necessary
258        // -- from the top of the function, folderStr = work folder, detStr = "FLRTB" or other type of averaging
259        //
260
261        nq = numpnts(qBin_qxqy)
262        Make/O/D/N=(nq)  $(folderPath+":"+"sigmaQ_qxqy"+"_"+detStr)
263        Make/O/D/N=(nq)  $(folderPath+":"+"qBar_qxqy"+"_"+detStr)
264        Make/O/D/N=(nq)  $(folderPath+":"+"fSubS_qxqy"+"_"+detStr)
265        Wave sigmaq = $(folderPath+":"+"sigmaQ_qxqy_"+detStr)
266        Wave qbar = $(folderPath+":"+"qBar_qxqy_"+detStr)
267        Wave fsubs = $(folderPath+":"+"fSubS_qxqy_"+detStr)
268
269// TODO
270// -- these are DUMMY VALUES!!!
271        sigmaq = -delQy
272        qbar = -delQy
273        fsubs = -delQy
274
275
276        SetDataFolder root:
277       
278        return(0)
279End
280
281
282// TODO -- update to new folder structure
283// unused -- update if necessary
284Proc CopyIQWaves()
285
286        SetDataFolder root:Packages:NIST:VSANS:VCALC   
287
288        duplicate/O iBin_qxqy_B iBin_qxqy_B_pin
289        duplicate/O qBin_qxqy_B qBin_qxqy_B_pin
290       
291        duplicate/O iBin_qxqy_MR iBin_qxqy_MR_pin
292        duplicate/O qBin_qxqy_MR qBin_qxqy_MR_pin
293        duplicate/O iBin_qxqy_MT iBin_qxqy_MT_pin
294        duplicate/O qBin_qxqy_MT qBin_qxqy_MT_pin
295        duplicate/O iBin_qxqy_ML iBin_qxqy_ML_pin
296        duplicate/O qBin_qxqy_ML qBin_qxqy_ML_pin
297        duplicate/O iBin_qxqy_MB iBin_qxqy_MB_pin
298        duplicate/O qBin_qxqy_MB qBin_qxqy_MB_pin
299
300        duplicate/O iBin_qxqy_FR iBin_qxqy_FR_pin
301        duplicate/O qBin_qxqy_FR qBin_qxqy_FR_pin
302        duplicate/O iBin_qxqy_FT iBin_qxqy_FT_pin
303        duplicate/O qBin_qxqy_FT qBin_qxqy_FT_pin
304        duplicate/O iBin_qxqy_FL iBin_qxqy_FL_pin
305        duplicate/O qBin_qxqy_FL qBin_qxqy_FL_pin
306        duplicate/O iBin_qxqy_FB iBin_qxqy_FB_pin
307        duplicate/O qBin_qxqy_FB qBin_qxqy_FB_pin
308
309        SetDataFolder root:
310End
311
312// TODO -- update to new folder structure
313// unused -- update if necessary
314Window slit_vs_pin_graph() : Graph
315        PauseUpdate; Silent 1           // building window...
316        String fldrSav0= GetDataFolder(1)
317        SetDataFolder root:Packages:NIST:VSANS:VCALC:
318        Display /W=(1296,44,1976,696) iBin_qxqy_B vs qBin_qxqy_B
319        AppendToGraph iBin_qxqy_B_pin vs qBin_qxqy_B_pin
320        AppendToGraph iBin_qxqy_MR vs qBin_qxqy_MR
321        AppendToGraph iBin_qxqy_MR_pin vs qBin_qxqy_MR_pin
322        AppendToGraph iBin_qxqy_MT vs qBin_qxqy_MT
323        AppendToGraph iBin_qxqy_MT_pin vs qBin_qxqy_MT_pin
324        AppendToGraph iBin_qxqy_FR vs qBin_qxqy_FR
325        AppendToGraph iBin_qxqy_FR_pin vs qBin_qxqy_FR_pin
326        AppendToGraph iBin_qxqy_FT vs qBin_qxqy_FT
327        AppendToGraph iBin_qxqy_FT_pin vs qBin_qxqy_FT_pin
328        SetDataFolder fldrSav0
329        ModifyGraph mode=4
330        ModifyGraph marker=19
331        ModifyGraph lSize=2
332        ModifyGraph rgb(iBin_qxqy_B)=(0,0,0),rgb(iBin_qxqy_B_pin)=(65535,16385,16385),rgb(iBin_qxqy_MR)=(2,39321,1)
333        ModifyGraph rgb(iBin_qxqy_MR_pin)=(0,0,65535),rgb(iBin_qxqy_MT)=(39321,1,31457)
334        ModifyGraph rgb(iBin_qxqy_MT_pin)=(48059,48059,48059),rgb(iBin_qxqy_FR)=(65535,32768,32768)
335        ModifyGraph rgb(iBin_qxqy_FR_pin)=(0,65535,0),rgb(iBin_qxqy_FT)=(16385,65535,65535)
336        ModifyGraph rgb(iBin_qxqy_FT_pin)=(65535,32768,58981)
337        ModifyGraph msize=2
338        ModifyGraph grid=1
339        ModifyGraph log=1
340        ModifyGraph mirror=1
341        SetAxis bottom 1e-05,0.05287132
342        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"
343        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"
344        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"
345EndMacro
Note: See TracBrowser for help on using the repository browser.