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

Last change on this file since 1141 was 1141, checked in by srkline, 3 years ago

updates to the beam center corrections that take into account the lateral variation of zero points of the tubes. This is an important correction to the reference beam center, especially when the graphite monochrmoator is used.

various bug fixes included as part of reduction , largely to catch missing or incorrect information in the file header.

New "patch" functions have been added for number of guides, beam stops, source aperture, etc. which were missing in the file and caused the resolution calculation to be totally incorrect.

File size: 15.2 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//      ii = trunc(ntube/4)             //random tube number
268//     
269//      Make/O/D/N=(nYpix) tmpTube,tmpMaskTube
270//      tmpTube = qy[ii][p]
271//      tmpMaskTube = mask[ii][p]
272//     
273//      // along the tube, keep the value, or set to NaN if masked
274//      tmpTube = tmpMaskTube == 0 ? tmpTube : NaN
275//      WaveStats/Q tmpTube
276//      Print V_max
277//      Print V_min     
278//      delQy = abs(V_max - V_min)
279
280// DONE: use only dQy for the portion of the detector that was not masked!
281
282// get delQy from the portion of the detector that is not masked out
283        delQy = V_setDeltaQy_Slit(folderStr,detStr)
284       
285
286/// TODO -- this is not necessary, but just for getting the I(Q) display to look "pretty"
287// clean out the near-zero Q point in the T/B  and Back detectors
288//      qBin_qxqy = (abs(qBin_qxqy[p][q]) < 1e-5) ? NaN : qBin_qxqy[p][q]                       
289
290
291// clear out zero data values before exiting...
292        // find the last non-zero point, working backwards
293        val = numpnts(iBin_qxqy)
294        do
295                val -= 1
296        while((iBin_qxqy[val] == 0) && val > 0)
297       
298        DeletePoints val, nq-val, iBin_qxqy,qBin_qxqy,eBin_qxqy
299
300
301// utility function to remove NaN values from the waves
302        V_RemoveNaNsQIS(qBin_qxqy, iBin_qxqy, eBin_qxqy)
303
304// work forwards? this doesn't work...
305//      val = -1
306//      do
307//              val += 1
308//      while(nBin_qxqy[val] == 0 && val < numpnts(iBin_qxqy)-1)       
309//      DeletePoints 0, val, iBin_qxqy,qBin_qxqy,eBin_qxqy
310
311
312        // TODO:
313        // -- This is where I calculate the resolution in SANS (see CircSectAve)
314        // -- use the isVCALC flag to exclude VCALC from the resolution calculation if necessary
315        // -- from the top of the function, folderStr = work folder, detStr = "FLRTB" or other type of averaging
316        //
317
318        nq = numpnts(qBin_qxqy)
319        Make/O/D/N=(nq)  $(folderPath+":"+"sigmaQ_qxqy"+"_"+detStr)
320        Make/O/D/N=(nq)  $(folderPath+":"+"qBar_qxqy"+"_"+detStr)
321        Make/O/D/N=(nq)  $(folderPath+":"+"fSubS_qxqy"+"_"+detStr)
322        Wave sigmaq = $(folderPath+":"+"sigmaQ_qxqy_"+detStr)
323        Wave qbar = $(folderPath+":"+"qBar_qxqy_"+detStr)
324        Wave fsubs = $(folderPath+":"+"fSubS_qxqy_"+detStr)
325
326// TODO:
327// -- calculate the slit resolution here. don't really have any idea how to represent this in VSANS
328//    since it's not an infinite slit, and not anything that I expect could be represented as a Gaussian
329// -- there is also wavelength smearing present
330
331
332
333// ASSUMPTION: As a first approximation, ignore the wavelength smearing component       
334        // TODO : do I use 1/2 of dQy due to the symmetry of my smearing calculation?   (as infinite slit)
335        sigmaq = -delQy/2
336        qbar = -delQy/2
337        fsubs = -delQy/2
338
339
340        SetDataFolder root:
341       
342        return(0)
343End
344
345
346// TODO -- update to new folder structure
347// unused -- update if necessary
348Proc CopyIQWaves()
349
350        SetDataFolder root:Packages:NIST:VSANS:VCALC   
351
352        duplicate/O iBin_qxqy_B iBin_qxqy_B_pin
353        duplicate/O qBin_qxqy_B qBin_qxqy_B_pin
354       
355        duplicate/O iBin_qxqy_MR iBin_qxqy_MR_pin
356        duplicate/O qBin_qxqy_MR qBin_qxqy_MR_pin
357        duplicate/O iBin_qxqy_MT iBin_qxqy_MT_pin
358        duplicate/O qBin_qxqy_MT qBin_qxqy_MT_pin
359        duplicate/O iBin_qxqy_ML iBin_qxqy_ML_pin
360        duplicate/O qBin_qxqy_ML qBin_qxqy_ML_pin
361        duplicate/O iBin_qxqy_MB iBin_qxqy_MB_pin
362        duplicate/O qBin_qxqy_MB qBin_qxqy_MB_pin
363
364        duplicate/O iBin_qxqy_FR iBin_qxqy_FR_pin
365        duplicate/O qBin_qxqy_FR qBin_qxqy_FR_pin
366        duplicate/O iBin_qxqy_FT iBin_qxqy_FT_pin
367        duplicate/O qBin_qxqy_FT qBin_qxqy_FT_pin
368        duplicate/O iBin_qxqy_FL iBin_qxqy_FL_pin
369        duplicate/O qBin_qxqy_FL qBin_qxqy_FL_pin
370        duplicate/O iBin_qxqy_FB iBin_qxqy_FB_pin
371        duplicate/O qBin_qxqy_FB qBin_qxqy_FB_pin
372
373        SetDataFolder root:
374End
375
376// TODO -- update to new folder structure
377// unused -- update if necessary
378Window slit_vs_pin_graph() : Graph
379        PauseUpdate; Silent 1           // building window...
380        String fldrSav0= GetDataFolder(1)
381        SetDataFolder root:Packages:NIST:VSANS:VCALC:
382        Display /W=(1296,44,1976,696) iBin_qxqy_B vs qBin_qxqy_B
383        AppendToGraph iBin_qxqy_B_pin vs qBin_qxqy_B_pin
384        AppendToGraph iBin_qxqy_MR vs qBin_qxqy_MR
385        AppendToGraph iBin_qxqy_MR_pin vs qBin_qxqy_MR_pin
386        AppendToGraph iBin_qxqy_MT vs qBin_qxqy_MT
387        AppendToGraph iBin_qxqy_MT_pin vs qBin_qxqy_MT_pin
388        AppendToGraph iBin_qxqy_FR vs qBin_qxqy_FR
389        AppendToGraph iBin_qxqy_FR_pin vs qBin_qxqy_FR_pin
390        AppendToGraph iBin_qxqy_FT vs qBin_qxqy_FT
391        AppendToGraph iBin_qxqy_FT_pin vs qBin_qxqy_FT_pin
392        SetDataFolder fldrSav0
393        ModifyGraph mode=4
394        ModifyGraph marker=19
395        ModifyGraph lSize=2
396        ModifyGraph rgb(iBin_qxqy_B)=(0,0,0),rgb(iBin_qxqy_B_pin)=(65535,16385,16385),rgb(iBin_qxqy_MR)=(2,39321,1)
397        ModifyGraph rgb(iBin_qxqy_MR_pin)=(0,0,65535),rgb(iBin_qxqy_MT)=(39321,1,31457)
398        ModifyGraph rgb(iBin_qxqy_MT_pin)=(48059,48059,48059),rgb(iBin_qxqy_FR)=(65535,32768,32768)
399        ModifyGraph rgb(iBin_qxqy_FR_pin)=(0,65535,0),rgb(iBin_qxqy_FT)=(16385,65535,65535)
400        ModifyGraph rgb(iBin_qxqy_FT_pin)=(65535,32768,58981)
401        ModifyGraph msize=2
402        ModifyGraph grid=1
403        ModifyGraph log=1
404        ModifyGraph mirror=1
405        SetAxis bottom 1e-05,0.05287132
406        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"
407        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"
408        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"
409EndMacro
410
411
412// finds the delta Qy for the slit height from the data that is not masked
413//
414Function V_setDeltaQy_Slit(folderStr,detStr)
415        String folderStr,detStr
416       
417        Variable delQy,min_qy,max_qy
418
419        String folderPath = "root:Packages:NIST:VSANS:"+folderStr
420        String instPath = ":entry:instrument:detector_"
421
422        WAVE qy = $(folderPath+instPath+detStr+":qy_"+detStr)
423        WAVE mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
424
425        Duplicate/O qy tmp_qy
426        // for the minimum
427        tmp_qy = (mask == 0) ? qY : 1e6
428        min_qy = WaveMin(tmp_qy)
429
430        // for the maximum
431        tmp_qy = (mask == 0) ? qY : -1e6
432        max_qy = WaveMax(tmp_qy)
433
434        KillWaves/Z tmp_qy     
435       
436        delQy = max_qy-min_qy
437       
438        return(delQy)
439end
Note: See TracBrowser for help on using the repository browser.