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

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

more utilities and bug fixes to handle:
(1) generation of DIV files
(2) generation and loading of slit-smeared VSANS data sets

Loading of VSANS data sets and generating the resolution smearing matrix required adding a conditional compile statement to the PlotUtilsMacro? file, to distinguish VSANS data - crudely done now by checking to see that VSANS is defined (which will NOT be the case for a "normal" Analysis experiment. This will need to be revisited in the future.

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// ASSUMPTION: As a first approximation, ignore the wavelength smearing component       
271        // TODO : do I use 1/2 of dQy due to the symmetry of my smearing calculation?   
272        sigmaq = -delQy/2
273        qbar = -delQy/2
274        fsubs = -delQy/2
275
276
277        SetDataFolder root:
278       
279        return(0)
280End
281
282
283// TODO -- update to new folder structure
284// unused -- update if necessary
285Proc CopyIQWaves()
286
287        SetDataFolder root:Packages:NIST:VSANS:VCALC   
288
289        duplicate/O iBin_qxqy_B iBin_qxqy_B_pin
290        duplicate/O qBin_qxqy_B qBin_qxqy_B_pin
291       
292        duplicate/O iBin_qxqy_MR iBin_qxqy_MR_pin
293        duplicate/O qBin_qxqy_MR qBin_qxqy_MR_pin
294        duplicate/O iBin_qxqy_MT iBin_qxqy_MT_pin
295        duplicate/O qBin_qxqy_MT qBin_qxqy_MT_pin
296        duplicate/O iBin_qxqy_ML iBin_qxqy_ML_pin
297        duplicate/O qBin_qxqy_ML qBin_qxqy_ML_pin
298        duplicate/O iBin_qxqy_MB iBin_qxqy_MB_pin
299        duplicate/O qBin_qxqy_MB qBin_qxqy_MB_pin
300
301        duplicate/O iBin_qxqy_FR iBin_qxqy_FR_pin
302        duplicate/O qBin_qxqy_FR qBin_qxqy_FR_pin
303        duplicate/O iBin_qxqy_FT iBin_qxqy_FT_pin
304        duplicate/O qBin_qxqy_FT qBin_qxqy_FT_pin
305        duplicate/O iBin_qxqy_FL iBin_qxqy_FL_pin
306        duplicate/O qBin_qxqy_FL qBin_qxqy_FL_pin
307        duplicate/O iBin_qxqy_FB iBin_qxqy_FB_pin
308        duplicate/O qBin_qxqy_FB qBin_qxqy_FB_pin
309
310        SetDataFolder root:
311End
312
313// TODO -- update to new folder structure
314// unused -- update if necessary
315Window slit_vs_pin_graph() : Graph
316        PauseUpdate; Silent 1           // building window...
317        String fldrSav0= GetDataFolder(1)
318        SetDataFolder root:Packages:NIST:VSANS:VCALC:
319        Display /W=(1296,44,1976,696) iBin_qxqy_B vs qBin_qxqy_B
320        AppendToGraph iBin_qxqy_B_pin vs qBin_qxqy_B_pin
321        AppendToGraph iBin_qxqy_MR vs qBin_qxqy_MR
322        AppendToGraph iBin_qxqy_MR_pin vs qBin_qxqy_MR_pin
323        AppendToGraph iBin_qxqy_MT vs qBin_qxqy_MT
324        AppendToGraph iBin_qxqy_MT_pin vs qBin_qxqy_MT_pin
325        AppendToGraph iBin_qxqy_FR vs qBin_qxqy_FR
326        AppendToGraph iBin_qxqy_FR_pin vs qBin_qxqy_FR_pin
327        AppendToGraph iBin_qxqy_FT vs qBin_qxqy_FT
328        AppendToGraph iBin_qxqy_FT_pin vs qBin_qxqy_FT_pin
329        SetDataFolder fldrSav0
330        ModifyGraph mode=4
331        ModifyGraph marker=19
332        ModifyGraph lSize=2
333        ModifyGraph rgb(iBin_qxqy_B)=(0,0,0),rgb(iBin_qxqy_B_pin)=(65535,16385,16385),rgb(iBin_qxqy_MR)=(2,39321,1)
334        ModifyGraph rgb(iBin_qxqy_MR_pin)=(0,0,65535),rgb(iBin_qxqy_MT)=(39321,1,31457)
335        ModifyGraph rgb(iBin_qxqy_MT_pin)=(48059,48059,48059),rgb(iBin_qxqy_FR)=(65535,32768,32768)
336        ModifyGraph rgb(iBin_qxqy_FR_pin)=(0,65535,0),rgb(iBin_qxqy_FT)=(16385,65535,65535)
337        ModifyGraph rgb(iBin_qxqy_FT_pin)=(65535,32768,58981)
338        ModifyGraph msize=2
339        ModifyGraph grid=1
340        ModifyGraph log=1
341        ModifyGraph mirror=1
342        SetAxis bottom 1e-05,0.05287132
343        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"
344        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"
345        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"
346EndMacro
Note: See TracBrowser for help on using the repository browser.