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

Last change on this file since 1109 was 1109, checked in by srkline, 4 years ago

Many changes:

Made the VCALC panel aware of all of the binning options
Corrected the behavior of the VCALC preset conditions
Adjusted how the Slit data is binned so that there are not duplicated q-values in the output

Made Absolute scaling aware of the back detector. Now the ABS String in the protocol has a second
set of scaling constants tagged with "_B" for the back detector. There is an added button
on the protocol panel to set the second set of constants. For the back detector, the read noise
is subtracted by reading it from the empty beam file (shifting over to the right by one box width)
All of the associated abs procedures are now aware of this.
More error checking needs to be added.

Back detector image is now shifted upon loading of the data. the default mask takes this into account
and masks out the padded (zero) regions.

in the protocol, DIV and MSK do not use grep any longer. it was just way too slow. Now it depends on

the file name having DIV or MASK respectively.



Raw data files can now be added together, in the usual way from the protocol panel.



File size: 14.5 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)                //for each tube...
182                sum_inten = 0                   // initialize the sum
183                sum_n = 0
184                sum_inten2 = 0
185               
186                for(jj=0;jj<nYpix;jj+=1)                        //sum along y...
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
202                avesq = (sum_inten/sum_n)^2
203                aveisq = sum_inten2/sum_n
204                var = aveisq-avesq
205                if(var<=0)
206                        eBin_qxqy[ii] = 1e-6
207                else
208                        eBin_qxqy[ii] = sqrt(var/(sum_n - 1))
209                endif
210
211        endfor
212
213// x- use only the Qx component in the y-center of the detector, not Qtotal
214// if the detectors are "L", then the values are all negative, so take the absolute value here
215        qBin_qxqy =  abs(qx[p][nYpix/2])               
216
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 panel data
220        Sort qBin_qxqy, qBin_qxqy,iBin_qxqy,eBin_qxqy
221
222// average and get rid of the duplicate q-values from the L and R sides
223        Variable q1,q2,tol
224        tol = 0.001             // 0.1 %
225        q1 = qBin_qxqy[0]
226        ii=0
227        do
228                q2 = qBin_qxqy[ii+1]
229                if(V_CloseEnough(q1,q2,q1*tol))
230                        // check to be sure that both values are actually real numbers before trying to average
231                        if(numtype(iBin_qxqy[ii])==0 && numtype(iBin_qxqy[ii+1])==0)            //==0 => real number
232                                iBin_qxqy[ii] = (iBin_qxqy[ii] + iBin_qxqy[ii+1])/2             //both OK
233                        endif
234                        if(numtype(iBin_qxqy[ii])==0 && numtype(iBin_qxqy[ii+1])!=0)            //==0 => real number
235                                iBin_qxqy[ii] = iBin_qxqy[ii]           //one OK
236                        endif
237                        if(numtype(iBin_qxqy[ii])!=0 && numtype(iBin_qxqy[ii+1])==0)            //==0 => real number
238                                iBin_qxqy[ii] = iBin_qxqy[ii+1]         //other OK
239                        endif
240                        if(numtype(iBin_qxqy[ii])!=0 && numtype(iBin_qxqy[ii+1])!=0)            //==0 => real number
241                                iBin_qxqy[ii] = (iBin_qxqy[ii])         // both NaN, get rid of it later
242                        endif
243               
244                        if(numtype(eBin_qxqy[ii])==0 && numtype(eBin_qxqy[ii+1])==0)            //==0 => real number
245                                eBin_qxqy[ii] = sqrt(eBin_qxqy[ii]^2 + eBin_qxqy[ii+1]^2)               //both OK
246                        endif
247                        if(numtype(eBin_qxqy[ii])==0 && numtype(eBin_qxqy[ii+1])!=0)            //==0 => real number
248                                eBin_qxqy[ii] = eBin_qxqy[ii]           //one OK
249                        endif
250                        if(numtype(eBin_qxqy[ii])!=0 && numtype(eBin_qxqy[ii+1])==0)            //==0 => real number
251                                eBin_qxqy[ii] = eBin_qxqy[ii+1]         //other OK
252                        endif
253                        if(numtype(eBin_qxqy[ii])!=0 && numtype(eBin_qxqy[ii+1])!=0)            //==0 => real number
254                                eBin_qxqy[ii] = (eBin_qxqy[ii])         // both NaN, get rid of it later
255                        endif
256                       
257                        DeletePoints ii+1, 1, qBin_qxqy,iBin_qxqy,eBin_qxqy,iBin2_qxqy,nBin_qxqy,eBin2D_qxqy
258                else
259                        ii+=1
260                        q1 = q2
261                endif
262        while(ii<numpnts(qBin_qxqy)-2)
263
264        // TODO: do I use dQ for the height of the panel?
265        // TODO : do I use 1/2 of dQ due to the symmetry of my smearing calculation?
266
267// TODO: use only dQy for the portion of the detector that was not masked!
268       
269        ii = trunc(ntube/4)             //random tube number
270       
271        Make/O/D/N=(nYpix) tmpTube,tmpMaskTube
272        tmpTube = qy[ii][p]
273        tmpMaskTube = mask[ii][p]
274       
275        // along the tube, keep the value, or set to NaN if masked
276        tmpTube = tmpMaskTube == 0 ? tmpTube : NaN
277        WaveStats/Q tmpTube
278//      Print V_max
279//      Print V_min
280       
281        delQy = abs(V_max - V_min)
282//print delQy
283       
284        // not quite correct - this uses the whole detector height, but there is some masked out
285//      delQy = abs(qy[0][nYpix-1] - qy[0][0])
286
287//      iBin_qxqy *= delQy
288//      eBin_qxqy *= delQy
289       
290
291/// TODO -- this is not necessary, but just for getting the I(Q) display to look "pretty"
292// clean out the near-zero Q point in the T/B  and Back detectors
293//      qBin_qxqy = (abs(qBin_qxqy[p][q]) < 1e-5) ? NaN : qBin_qxqy[p][q]                       
294
295
296// clear out zero data values before exiting...
297        // find the last non-zero point, working backwards
298        val = numpnts(iBin_qxqy)
299        do
300                val -= 1
301        while((iBin_qxqy[val] == 0) && val > 0)
302       
303        DeletePoints val, nq-val, iBin_qxqy,qBin_qxqy,eBin_qxqy
304
305
306// utility function to remove NaN values from the waves
307        V_RemoveNaNsQIS(qBin_qxqy, iBin_qxqy, eBin_qxqy)
308
309// work forwards? this doesn't work...
310//      val = -1
311//      do
312//              val += 1
313//      while(nBin_qxqy[val] == 0 && val < numpnts(iBin_qxqy)-1)       
314//      DeletePoints 0, val, iBin_qxqy,qBin_qxqy,eBin_qxqy
315
316
317        // TODO:
318        // -- This is where I calculate the resolution in SANS (see CircSectAve)
319        // -- use the isVCALC flag to exclude VCALC from the resolution calculation if necessary
320        // -- from the top of the function, folderStr = work folder, detStr = "FLRTB" or other type of averaging
321        //
322
323        nq = numpnts(qBin_qxqy)
324        Make/O/D/N=(nq)  $(folderPath+":"+"sigmaQ_qxqy"+"_"+detStr)
325        Make/O/D/N=(nq)  $(folderPath+":"+"qBar_qxqy"+"_"+detStr)
326        Make/O/D/N=(nq)  $(folderPath+":"+"fSubS_qxqy"+"_"+detStr)
327        Wave sigmaq = $(folderPath+":"+"sigmaQ_qxqy_"+detStr)
328        Wave qbar = $(folderPath+":"+"qBar_qxqy_"+detStr)
329        Wave fsubs = $(folderPath+":"+"fSubS_qxqy_"+detStr)
330
331// TODO:
332// -- calculate the slit resolution here. don't really have any idea how to represent this in VSANS
333//    since it's not an infinite slit, and not anything that I expect could be represented as a Gaussian
334// -- there is also wavelength smearing present
335
336
337
338// ASSUMPTION: As a first approximation, ignore the wavelength smearing component       
339        // TODO : do I use 1/2 of dQy due to the symmetry of my smearing calculation?   (as infinite slit)
340        sigmaq = -delQy/2
341        qbar = -delQy/2
342        fsubs = -delQy/2
343
344
345        SetDataFolder root:
346       
347        return(0)
348End
349
350
351// TODO -- update to new folder structure
352// unused -- update if necessary
353Proc CopyIQWaves()
354
355        SetDataFolder root:Packages:NIST:VSANS:VCALC   
356
357        duplicate/O iBin_qxqy_B iBin_qxqy_B_pin
358        duplicate/O qBin_qxqy_B qBin_qxqy_B_pin
359       
360        duplicate/O iBin_qxqy_MR iBin_qxqy_MR_pin
361        duplicate/O qBin_qxqy_MR qBin_qxqy_MR_pin
362        duplicate/O iBin_qxqy_MT iBin_qxqy_MT_pin
363        duplicate/O qBin_qxqy_MT qBin_qxqy_MT_pin
364        duplicate/O iBin_qxqy_ML iBin_qxqy_ML_pin
365        duplicate/O qBin_qxqy_ML qBin_qxqy_ML_pin
366        duplicate/O iBin_qxqy_MB iBin_qxqy_MB_pin
367        duplicate/O qBin_qxqy_MB qBin_qxqy_MB_pin
368
369        duplicate/O iBin_qxqy_FR iBin_qxqy_FR_pin
370        duplicate/O qBin_qxqy_FR qBin_qxqy_FR_pin
371        duplicate/O iBin_qxqy_FT iBin_qxqy_FT_pin
372        duplicate/O qBin_qxqy_FT qBin_qxqy_FT_pin
373        duplicate/O iBin_qxqy_FL iBin_qxqy_FL_pin
374        duplicate/O qBin_qxqy_FL qBin_qxqy_FL_pin
375        duplicate/O iBin_qxqy_FB iBin_qxqy_FB_pin
376        duplicate/O qBin_qxqy_FB qBin_qxqy_FB_pin
377
378        SetDataFolder root:
379End
380
381// TODO -- update to new folder structure
382// unused -- update if necessary
383Window slit_vs_pin_graph() : Graph
384        PauseUpdate; Silent 1           // building window...
385        String fldrSav0= GetDataFolder(1)
386        SetDataFolder root:Packages:NIST:VSANS:VCALC:
387        Display /W=(1296,44,1976,696) iBin_qxqy_B vs qBin_qxqy_B
388        AppendToGraph iBin_qxqy_B_pin vs qBin_qxqy_B_pin
389        AppendToGraph iBin_qxqy_MR vs qBin_qxqy_MR
390        AppendToGraph iBin_qxqy_MR_pin vs qBin_qxqy_MR_pin
391        AppendToGraph iBin_qxqy_MT vs qBin_qxqy_MT
392        AppendToGraph iBin_qxqy_MT_pin vs qBin_qxqy_MT_pin
393        AppendToGraph iBin_qxqy_FR vs qBin_qxqy_FR
394        AppendToGraph iBin_qxqy_FR_pin vs qBin_qxqy_FR_pin
395        AppendToGraph iBin_qxqy_FT vs qBin_qxqy_FT
396        AppendToGraph iBin_qxqy_FT_pin vs qBin_qxqy_FT_pin
397        SetDataFolder fldrSav0
398        ModifyGraph mode=4
399        ModifyGraph marker=19
400        ModifyGraph lSize=2
401        ModifyGraph rgb(iBin_qxqy_B)=(0,0,0),rgb(iBin_qxqy_B_pin)=(65535,16385,16385),rgb(iBin_qxqy_MR)=(2,39321,1)
402        ModifyGraph rgb(iBin_qxqy_MR_pin)=(0,0,65535),rgb(iBin_qxqy_MT)=(39321,1,31457)
403        ModifyGraph rgb(iBin_qxqy_MT_pin)=(48059,48059,48059),rgb(iBin_qxqy_FR)=(65535,32768,32768)
404        ModifyGraph rgb(iBin_qxqy_FR_pin)=(0,65535,0),rgb(iBin_qxqy_FT)=(16385,65535,65535)
405        ModifyGraph rgb(iBin_qxqy_FT_pin)=(65535,32768,58981)
406        ModifyGraph msize=2
407        ModifyGraph grid=1
408        ModifyGraph log=1
409        ModifyGraph mirror=1
410        SetAxis bottom 1e-05,0.05287132
411        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"
412        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"
413        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"
414EndMacro
Note: See TracBrowser for help on using the repository browser.