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

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

updating the IgorVersion? pragma to v7.0 for all files to be consistent.

File size: 15.2 KB
Line 
1#pragma rtGlobals=3             // Use modern global access method and strict wave access.
2#pragma IgorVersion = 7.00
3
4
5//////////////////
6//
7// Procedures to average data taken in "slit" aperture geometry. The fake data on the detector panels
8//  is as would be collected in PINHOLE geometry - I do not currently have a simulation for slit
9//  apertures (this would need to be MonteCarlo) - the I(q) averaging here gives the I(q) that you
10//  would measure in 1D using slit geometry. It is done by simply summing the columns of data on each detector.
11//
12//////////////////
13//
14//
15// TODO - big question about averaging in this way...
16// can the T/B panels really be used at all for slit mode? - since there's a big "hole" in the scattering data
17// collected -- you're not getting the full column of data covering a wide range of Qy. L/R panels should be fine,
18//  although the different carriages cover very different "chunks" of Qy
19//
20// - best answer so far is to skip the T/B panels, and simply not use them
21//
22// TODO -- be sure that the absolute scaling of this is correct. Right now, something is off.
23//  -- so write up the documentation for this operation, and try to find the error in the process...
24//
25//
26/////////////////
27
28
29// TODO:
30// -- verify the error calculation
31// -- ? add in functionality to handle FLR and MLR cases (2 panels of data)
32//
33// seems backwards to call this "byRows", but this is the way that Igor indexes
34// LR banks are defined as (48,256) (n,m), sumRows gives sum w/ dimension (n x 1)
35//
36// updated to new folder structure Feb 2016
37// folderStr = RAW,SAM, VCALC or other
38// detStr is the panel identifer "ML", etc.
39Function VC_fBinDetector_byRows(folderStr,detStr)
40        String folderStr,detStr
41       
42//      SetDataFolder root:Packages:NIST:VSANS:VCALC   
43       
44        Variable pixSizeX,pixSizeY,delQy
45        Variable isVCALC=0,maskMissing,nsets=0
46       
47        maskMissing = 1         // set to zero if a mask is actually present
48       
49        if(cmpstr(folderStr,"VCALC") == 0)
50                isVCALC = 1
51        endif
52
53        String folderPath = "root:Packages:NIST:VSANS:"+folderStr
54        String instPath = ":entry:instrument:detector_"
55
56        strswitch(detStr)       // string switch       
57// only one panel, simply pick that panel and move on out of the switch
58                case "FL":
59                case "FR":
60                case "ML":
61                case "MR":
62                case "B":
63                        if(isVCALC)
64                                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)               // 2D detector data
65                                WAVE/Z iErr = $("asdf_iErr_"+detStr)                    // TODO: 2D errors -- may not exist, especially for simulation
66                        else
67                                Wave inten = V_getDetectorDataW(folderStr,detStr)
68                                Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
69                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
70                                if(WaveExists(mask) == 1)
71                                        maskMissing = 0
72                                endif
73                        endif
74                       
75                        nsets = 1
76                        break
77
78//              case "FLR":
79//              // detStr has multiple values now, so unfortuntely, I'm hard-wiring things...
80//                      if(isVCALC)
81//                              WAVE inten = $(folderPath+instPath+"FL"+":det_"+"FL")
82//                              WAVE/Z iErr = $("iErr_"+"FL")                   // 2D errors -- may not exist, especially for simulation               
83//                              WAVE inten2 = $(folderPath+instPath+"FR"+":det_"+"FR")
84//                              WAVE/Z iErr2 = $("iErr_"+"FR")                  // 2D errors -- may not exist, especially for simulation       
85//                      else
86//                              Wave inten = V_getDetectorDataW(folderStr,"FL")
87//                              Wave iErr = V_getDetectorDataErrW(folderStr,"FL")
88//                              Wave inten2 = V_getDetectorDataW(folderStr,"FR")
89//                              Wave iErr2 = V_getDetectorDataErrW(folderStr,"FR")
90//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FL"+":data")
91//                              Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FR"+":data")
92//                              if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
93//                                      maskMissing = 0
94//                              endif
95//                      endif   
96//
97//                      nsets = 2
98//                      break
99//                     
100//              case "MLR":
101//                      if(isVCALC)
102//                              WAVE inten = $(folderPath+instPath+"ML"+":det_"+"ML")
103//                              WAVE/Z iErr = $("iErr_"+"ML")                   // 2D errors -- may not exist, especially for simulation               
104//                              WAVE inten2 = $(folderPath+instPath+"MR"+":det_"+"MR")
105//                              WAVE/Z iErr2 = $("iErr_"+"MR")                  // 2D errors -- may not exist, especially for simulation       
106//                      else
107//                              Wave inten = V_getDetectorDataW(folderStr,"ML")
108//                              Wave iErr = V_getDetectorDataErrW(folderStr,"ML")
109//                              Wave inten2 = V_getDetectorDataW(folderStr,"MR")
110//                              Wave iErr2 = V_getDetectorDataErrW(folderStr,"MR")
111//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"ML"+":data")
112//                              Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MR"+":data")
113//                              if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
114//                                      maskMissing = 0
115//                              endif
116//                      endif   
117//             
118//                      nSets = 2
119//                      break                   
120
121                default:
122                        nSets = 0                                                       
123                        Print "ERROR   ---- type is not recognized "
124        endswitch
125
126        if(nSets == 0)
127                SetDataFolder root:
128                return(0)
129        endif
130
131
132        Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
133        Wave qx = $(folderPath+instPath+detStr+":qx_"+detStr)
134        Wave qy = $(folderPath+instPath+detStr+":qy_"+detStr)
135
136
137// delta Qx is set by the pixel X dimension of the detector, which is the limiting resolution
138//      delQx = abs(qx[0][0] - qx[1][0])
139       
140        Variable nq,val
141        nq = DimSize(inten,0)           //nq == the number of columns (x dimension)
142       
143//      SetDataFolder $(folderPath+instPath+detStr)     
144
145        Make/O/D/N=(nq)  $(folderPath+":"+"iBin_qxqy_"+detStr)
146        Make/O/D/N=(nq)  $(folderPath+":"+"qBin_qxqy_"+detStr)
147        Make/O/D/N=(nq)  $(folderPath+":"+"nBin_qxqy_"+detStr)
148        Make/O/D/N=(nq)  $(folderPath+":"+"iBin2_qxqy_"+detStr)
149        Make/O/D/N=(nq)  $(folderPath+":"+"eBin_qxqy_"+detStr)
150        Make/O/D/N=(nq)  $(folderPath+":"+"eBin2D_qxqy_"+detStr)
151       
152        Wave iBin_qxqy = $(folderPath+":"+"iBin_qxqy_"+detStr)
153        Wave qBin_qxqy = $(folderPath+":"+"qBin_qxqy_"+detStr)
154        Wave nBin_qxqy = $(folderPath+":"+"nBin_qxqy_"+detStr)
155        Wave iBin2_qxqy = $(folderPath+":"+"iBin2_qxqy_"+detStr)
156        Wave eBin_qxqy = $(folderPath+":"+"eBin_qxqy_"+detStr)
157        Wave eBin2D_qxqy = $(folderPath+":"+"eBin2D_qxqy_"+detStr)
158
159
160        iBin_qxqy = 0
161        iBin2_qxqy = 0
162        eBin_qxqy = 0
163        eBin2D_qxqy = 0
164        nBin_qxqy = 0   
165
166
167// sum the rows
168
169// MatrixOp would be fast, but I can't figure out how to apply the mask with sumRows???
170//      MatrixOp/O iBin_qxqy = sumRows(inten)   //automatically generates the destination
171//             
172//// how to properly calculate the error?
173//// This MatrixOp gives values way too large -- larger than the intensity (be sure to correct *delQy below...
174//
175//      MatrixOp/O tmp = sqrt(varCols(inten^t))         // variance: no varRows operation, so use the transpose of the matrix
176//      eBin_qxqy = tmp[0][p]
177
178        Variable ii,jj,ntube,nYpix,sum_inten, sum_n, sum_inten2,avesq,aveisq,var,mask_val
179        ntube = DimSize(inten,0)
180        nYpix = DimSize(inten,1)
181       
182        for(ii=0;ii<ntube;ii+=1)                //for each tube...
183                sum_inten = 0                   // initialize the sum
184                sum_n = 0
185                sum_inten2 = 0
186               
187                for(jj=0;jj<nYpix;jj+=1)                        //sum along y...
188                                val = inten[ii][jj]
189                                if(isVCALC || maskMissing)              // mask_val == 0 == keep, mask_val == 1 = YES, mask out the point
190                                        mask_val = 0
191                                else
192                                        mask_val = mask[ii][jj]
193                                endif
194                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
195                                        sum_inten += val
196                                        sum_n += 1
197                                        sum_inten2 += val*val
198                                endif
199                endfor
200                iBin_qxqy[ii] = sum_inten/sum_n         //the average value
201               
202
203                avesq = (sum_inten/sum_n)^2
204                aveisq = sum_inten2/sum_n
205                var = aveisq-avesq
206                if(var<=0)
207                        eBin_qxqy[ii] = 1e-6
208                else
209                        eBin_qxqy[ii] = sqrt(var/(sum_n - 1))
210                endif
211
212        endfor
213
214// x- use only the Qx component in the y-center of the detector, not Qtotal
215// if the detectors are "L", then the values are all negative, so take the absolute value here
216        qBin_qxqy =  abs(qx[p][nYpix/2])               
217
218
219// for the L panels, sort the q-values (and data) after the abs() step, otherwise the data is reversed
220// won't hurt to sort the R panel data
221        Sort qBin_qxqy, qBin_qxqy,iBin_qxqy,eBin_qxqy
222
223// average and get rid of the duplicate q-values from the L and R sides
224        Variable q1,q2,tol
225        tol = 0.001             // 0.1 %
226        q1 = qBin_qxqy[0]
227        ii=0
228        do
229                q2 = qBin_qxqy[ii+1]
230                if(V_CloseEnough(q1,q2,q1*tol))
231                        // check to be sure that both values are actually real numbers before trying to average
232                        if(numtype(iBin_qxqy[ii])==0 && numtype(iBin_qxqy[ii+1])==0)            //==0 => real number
233                                iBin_qxqy[ii] = (iBin_qxqy[ii] + iBin_qxqy[ii+1])/2             //both OK
234                        endif
235                        if(numtype(iBin_qxqy[ii])==0 && numtype(iBin_qxqy[ii+1])!=0)            //==0 => real number
236                                iBin_qxqy[ii] = iBin_qxqy[ii]           //one OK
237                        endif
238                        if(numtype(iBin_qxqy[ii])!=0 && numtype(iBin_qxqy[ii+1])==0)            //==0 => real number
239                                iBin_qxqy[ii] = iBin_qxqy[ii+1]         //other OK
240                        endif
241                        if(numtype(iBin_qxqy[ii])!=0 && numtype(iBin_qxqy[ii+1])!=0)            //==0 => real number
242                                iBin_qxqy[ii] = (iBin_qxqy[ii])         // both NaN, get rid of it later
243                        endif
244               
245                        if(numtype(eBin_qxqy[ii])==0 && numtype(eBin_qxqy[ii+1])==0)            //==0 => real number
246                                eBin_qxqy[ii] = sqrt(eBin_qxqy[ii]^2 + eBin_qxqy[ii+1]^2)               //both OK
247                        endif
248                        if(numtype(eBin_qxqy[ii])==0 && numtype(eBin_qxqy[ii+1])!=0)            //==0 => real number
249                                eBin_qxqy[ii] = eBin_qxqy[ii]           //one OK
250                        endif
251                        if(numtype(eBin_qxqy[ii])!=0 && numtype(eBin_qxqy[ii+1])==0)            //==0 => real number
252                                eBin_qxqy[ii] = eBin_qxqy[ii+1]         //other OK
253                        endif
254                        if(numtype(eBin_qxqy[ii])!=0 && numtype(eBin_qxqy[ii+1])!=0)            //==0 => real number
255                                eBin_qxqy[ii] = (eBin_qxqy[ii])         // both NaN, get rid of it later
256                        endif
257                       
258                        DeletePoints ii+1, 1, qBin_qxqy,iBin_qxqy,eBin_qxqy,iBin2_qxqy,nBin_qxqy,eBin2D_qxqy
259                else
260                        ii+=1
261                        q1 = q2
262                endif
263        while(ii<numpnts(qBin_qxqy)-2)
264
265        // TODO: do I use dQ for the height of the panel?
266        // TODO : do I use 1/2 of dQ due to the symmetry of my smearing calculation?
267
268//      ii = trunc(ntube/4)             //random tube number
269//     
270//      Make/O/D/N=(nYpix) tmpTube,tmpMaskTube
271//      tmpTube = qy[ii][p]
272//      tmpMaskTube = mask[ii][p]
273//     
274//      // along the tube, keep the value, or set to NaN if masked
275//      tmpTube = tmpMaskTube == 0 ? tmpTube : NaN
276//      WaveStats/Q tmpTube
277//      Print V_max
278//      Print V_min     
279//      delQy = abs(V_max - V_min)
280
281// DONE: use only dQy for the portion of the detector that was not masked!
282
283// get delQy from the portion of the detector that is not masked out
284        delQy = V_setDeltaQy_Slit(folderStr,detStr)
285       
286
287/// TODO -- this is not necessary, but just for getting the I(Q) display to look "pretty"
288// clean out the near-zero Q point in the T/B  and Back detectors
289//      qBin_qxqy = (abs(qBin_qxqy[p][q]) < 1e-5) ? NaN : qBin_qxqy[p][q]                       
290
291
292// clear out zero data values before exiting...
293        // find the last non-zero point, working backwards
294        val = numpnts(iBin_qxqy)
295        do
296                val -= 1
297        while((iBin_qxqy[val] == 0) && val > 0)
298       
299        DeletePoints val, nq-val, iBin_qxqy,qBin_qxqy,eBin_qxqy
300
301
302// utility function to remove NaN values from the waves
303        V_RemoveNaNsQIS(qBin_qxqy, iBin_qxqy, eBin_qxqy)
304
305// work forwards? this doesn't work...
306//      val = -1
307//      do
308//              val += 1
309//      while(nBin_qxqy[val] == 0 && val < numpnts(iBin_qxqy)-1)       
310//      DeletePoints 0, val, iBin_qxqy,qBin_qxqy,eBin_qxqy
311
312
313        // TODO:
314        // -- This is where I calculate the resolution in SANS (see CircSectAve)
315        // -- use the isVCALC flag to exclude VCALC from the resolution calculation if necessary
316        // -- from the top of the function, folderStr = work folder, detStr = "FLRTB" or other type of averaging
317        //
318
319        nq = numpnts(qBin_qxqy)
320        Make/O/D/N=(nq)  $(folderPath+":"+"sigmaQ_qxqy"+"_"+detStr)
321        Make/O/D/N=(nq)  $(folderPath+":"+"qBar_qxqy"+"_"+detStr)
322        Make/O/D/N=(nq)  $(folderPath+":"+"fSubS_qxqy"+"_"+detStr)
323        Wave sigmaq = $(folderPath+":"+"sigmaQ_qxqy_"+detStr)
324        Wave qbar = $(folderPath+":"+"qBar_qxqy_"+detStr)
325        Wave fsubs = $(folderPath+":"+"fSubS_qxqy_"+detStr)
326
327// TODO:
328// -- calculate the slit resolution here. don't really have any idea how to represent this in VSANS
329//    since it's not an infinite slit, and not anything that I expect could be represented as a Gaussian
330// -- there is also wavelength smearing present
331
332
333
334// ASSUMPTION: As a first approximation, ignore the wavelength smearing component       
335        // TODO : do I use 1/2 of dQy due to the symmetry of my smearing calculation?   (as infinite slit)
336        sigmaq = -delQy/2
337        qbar = -delQy/2
338        fsubs = -delQy/2
339
340
341        SetDataFolder root:
342       
343        return(0)
344End
345
346
347// TODO -- update to new folder structure
348// unused -- update if necessary
349Proc CopyIQWaves()
350
351        SetDataFolder root:Packages:NIST:VSANS:VCALC   
352
353        duplicate/O iBin_qxqy_B iBin_qxqy_B_pin
354        duplicate/O qBin_qxqy_B qBin_qxqy_B_pin
355       
356        duplicate/O iBin_qxqy_MR iBin_qxqy_MR_pin
357        duplicate/O qBin_qxqy_MR qBin_qxqy_MR_pin
358        duplicate/O iBin_qxqy_MT iBin_qxqy_MT_pin
359        duplicate/O qBin_qxqy_MT qBin_qxqy_MT_pin
360        duplicate/O iBin_qxqy_ML iBin_qxqy_ML_pin
361        duplicate/O qBin_qxqy_ML qBin_qxqy_ML_pin
362        duplicate/O iBin_qxqy_MB iBin_qxqy_MB_pin
363        duplicate/O qBin_qxqy_MB qBin_qxqy_MB_pin
364
365        duplicate/O iBin_qxqy_FR iBin_qxqy_FR_pin
366        duplicate/O qBin_qxqy_FR qBin_qxqy_FR_pin
367        duplicate/O iBin_qxqy_FT iBin_qxqy_FT_pin
368        duplicate/O qBin_qxqy_FT qBin_qxqy_FT_pin
369        duplicate/O iBin_qxqy_FL iBin_qxqy_FL_pin
370        duplicate/O qBin_qxqy_FL qBin_qxqy_FL_pin
371        duplicate/O iBin_qxqy_FB iBin_qxqy_FB_pin
372        duplicate/O qBin_qxqy_FB qBin_qxqy_FB_pin
373
374        SetDataFolder root:
375End
376
377// TODO -- update to new folder structure
378// unused -- update if necessary
379Window slit_vs_pin_graph() : Graph
380        PauseUpdate; Silent 1           // building window...
381        String fldrSav0= GetDataFolder(1)
382        SetDataFolder root:Packages:NIST:VSANS:VCALC:
383        Display /W=(1296,44,1976,696) iBin_qxqy_B vs qBin_qxqy_B
384        AppendToGraph iBin_qxqy_B_pin vs qBin_qxqy_B_pin
385        AppendToGraph iBin_qxqy_MR vs qBin_qxqy_MR
386        AppendToGraph iBin_qxqy_MR_pin vs qBin_qxqy_MR_pin
387        AppendToGraph iBin_qxqy_MT vs qBin_qxqy_MT
388        AppendToGraph iBin_qxqy_MT_pin vs qBin_qxqy_MT_pin
389        AppendToGraph iBin_qxqy_FR vs qBin_qxqy_FR
390        AppendToGraph iBin_qxqy_FR_pin vs qBin_qxqy_FR_pin
391        AppendToGraph iBin_qxqy_FT vs qBin_qxqy_FT
392        AppendToGraph iBin_qxqy_FT_pin vs qBin_qxqy_FT_pin
393        SetDataFolder fldrSav0
394        ModifyGraph mode=4
395        ModifyGraph marker=19
396        ModifyGraph lSize=2
397        ModifyGraph rgb(iBin_qxqy_B)=(0,0,0),rgb(iBin_qxqy_B_pin)=(65535,16385,16385),rgb(iBin_qxqy_MR)=(2,39321,1)
398        ModifyGraph rgb(iBin_qxqy_MR_pin)=(0,0,65535),rgb(iBin_qxqy_MT)=(39321,1,31457)
399        ModifyGraph rgb(iBin_qxqy_MT_pin)=(48059,48059,48059),rgb(iBin_qxqy_FR)=(65535,32768,32768)
400        ModifyGraph rgb(iBin_qxqy_FR_pin)=(0,65535,0),rgb(iBin_qxqy_FT)=(16385,65535,65535)
401        ModifyGraph rgb(iBin_qxqy_FT_pin)=(65535,32768,58981)
402        ModifyGraph msize=2
403        ModifyGraph grid=1
404        ModifyGraph log=1
405        ModifyGraph mirror=1
406        SetAxis bottom 1e-05,0.05287132
407        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"
408        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"
409        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"
410EndMacro
411
412
413// finds the delta Qy for the slit height from the data that is not masked
414//
415Function V_setDeltaQy_Slit(folderStr,detStr)
416        String folderStr,detStr
417       
418        Variable delQy,min_qy,max_qy
419
420        String folderPath = "root:Packages:NIST:VSANS:"+folderStr
421        String instPath = ":entry:instrument:detector_"
422
423        WAVE qy = $(folderPath+instPath+detStr+":qy_"+detStr)
424        WAVE mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
425
426        Duplicate/O qy tmp_qy
427        // for the minimum
428        tmp_qy = (mask == 0) ? qY : 1e6
429        min_qy = WaveMin(tmp_qy)
430
431        // for the maximum
432        tmp_qy = (mask == 0) ? qY : -1e6
433        max_qy = WaveMax(tmp_qy)
434
435        KillWaves/Z tmp_qy     
436       
437        delQy = max_qy-min_qy
438       
439        return(delQy)
440end
Note: See TracBrowser for help on using the repository browser.