source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_TubeAdjustments.ipf @ 964

Last change on this file since 964 was 964, checked in by srkline, 7 years ago

Adding a procedure file with routines that allow the fitting of calibration data from the detector tubes. Also allows different ways to plot the data from tubes, shifting the pixels in different ways for display.

File size: 12.5 KB
Line 
1#pragma rtGlobals=3             // Use modern global access method and strict wave access.
2
3
4
5
6
7// the main routines are:
8
9//(1)
10//to get from individual tubes to an array
11//      Tubes_to_Array()                       
12
13//(2)
14// then to locate all of the peak positions
15//      MakeTableForPeaks(numTube,numPeak)             
16//      Identify_AllPeaks()
17//              AutoFindPeaksCustom()           // if Identify_AllPeaks  doesn't work -try this, setting the "noise" to 1 and smoothing to 2
18
19//(3)
20// fit to find all of the quadratic coefficients
21//      MakeTableForFitCoefs(numTube,numCoef)
22//      PlotFit_AllPeaks()
23
24
25//(4)
26// then pick a display method
27//
28//      Make_mm_tubes()
29//      Append_Adjusted_mm()
30//
31//      MakeMatrix_PixCenter()
32//      FillMatrix_Pix_Center(ind)
33//
34
35// The function most used externally is:
36// V_TubePix_to_mm(c0,c1,c2,pix)
37//
38// which will return the real space distance (in mm?) for a given pixel
39// and the tube's coefficients. The distance is relative to the zero position on the
40// detector (which was defined when the coefficients were determined)
41
42
43
44
45
46
47
48// (0) -- what I start with:
49// -- a table of the mm spacing of the slots (20 of them)
50// -- masked data from each of the (8) tubes
51// -- the table of slots may need to be corrected for parallax, depending on the geometry of the test
52// ** In the table of slots, pick a slot near the center, and SET that to ZERO. Then all of the other
53//   distances are relative to that zero point. This is a necessary reference point.
54//
55
56
57// (1) -- get the individual tubes into an array
58//
59//
60Proc Tubes_to_Array()
61        Make/O/D/N=(8,1127) pack
62        edit pack
63        display;appendimage pack
64        pack[0][] = tube1[q]
65        pack[1][] = tube2[q]
66        pack[2][] = tube3[q]
67        pack[3][] = tube4[q]
68        pack[4][] = tube5[q]
69        pack[5][] = tube6[q]
70        pack[6][] = tube7[q]
71        pack[7][] = tube8[q]
72        ModifyImage pack ctab= {*,*,ColdWarm,0}
73End
74
75// (2) -- for each of the tubes, find the x-position (in pixels) of each of the (20) peaks
76// -- load the Analysis Package "MultiPeakFit 2"
77//
78// automatically find the peaks (after including MultiPeakFit 2)
79//              AutomaticallyFindPeaks()
80//
81//-- or if having difficulty
82//              AutoFindPeaksCustom()           // try this, setting the "noise" to 1 and smoothing to 2
83//
84// -- if really having difficulty, you can do the "full" MultiPeak Fit
85//
86// -- If (hopefully) using the easy way, the results are in:
87// root:WA_PeakCentersY,root:WA_PeakCentersX
88//
89// -- so to see the results:
90//¥Edit/K=0  root:WA_PeakCentersY,root:WA_PeakCentersX
91//
92// -- then sort the results - they seem to be in no real order...
93//¥Sort WA_PeakCentersX WA_PeakCentersY,WA_PeakCentersX
94//
95Proc MakeTableForPeaks(numTube,numPeak)
96        Variable numTube,numPeak
97       
98        Make/O/D/N=(numPeak,numTube) PeakTableX,peakTableY              //*2 to store x-location and peak height (y)
99        Edit peakTableX
100End
101
102Proc Identify_AllPeaks()
103
104        Variable ii,numTubes=8
105        String str="tube"
106       
107        ii=1
108        do
109                Identify_Peaks(str+num2str(ii),ii-1)
110                ii+=1
111        while(ii<=numTubes)
112
113End
114
115Proc Identify_Peaks(tubeStr,ind)
116        String tubeStr
117        Variable ind
118       
119        // must use a wave of pixels rather than "calculated" -- if calculated is used it only
120        // returns integer values for the peak locations
121       
122//      AutomaticallyFindPeaks() //-- this is a function that doesn't take any parameters - so
123// I need to pull this from the WM function to call the worker directly
124        Variable pBegin=0, pEnd= numpnts($(tubeStr))-1
125        Variable/C estimates= EstPeakNoiseAndSmfact($(tubeStr),pBegin, pEnd)
126        Variable noiselevel=real(estimates)
127        Variable smoothingFactor=imag(estimates)
128        Variable maxPeaks = 20
129        Variable minPeakPercent = 10
130       
131        AutoFindPeaksWorker($(tubeStr), $("tube_pixel"), pBegin, pEnd, maxPeaks, minPeakPercent, noiseLevel, smoothingFactor)
132// end WM function call
133
134        Sort WA_PeakCentersX WA_PeakCentersY,WA_PeakCentersX
135       
136        peakTableX[][ind] = WA_PeakCentersX[p]
137        peakTableY[][ind] = WA_PeakCentersY[p]
138       
139End
140
141// -- save a copy of the root:WA_PeakCentersY,root:WA_PeakCentersX values
142//    for later in case the fitting failed, then you can go back and re-do
143//
144// -- then plot:
145//
146//      Display peak_spacing_mm_ctr vs WA_PeakCentersX
147//
148//      Then do a "QuickFit" of the peak position to the data using a polynomial of order 3 (= quadratic)
149//
150// result is in W_coef, W_sigma
151//
152// -- an example of the "quickFit" command is below, so it can be programmed rather than the menu every time
153//¥Display peak_spacing_mm_ctr vs WA_PeakCentersX3
154//¥CurveFit/M=2/W=0/TBOX=(0x310) poly 3, peak_spacing_mm_ctr/X=WA_PeakCentersX3/D
155//  fit_peak_spacing_mm_ctr= poly(W_coef,x)
156//  W_coef={-571.42,1.1135,-4.2444e-05}
157//  V_chisq= 8.5841;V_npnts= 20;V_numNaNs= 0;V_numINFs= 0;
158//  V_startRow= 0;V_endRow= 19;
159//  W_sigma={0.595,0.00246,2.15e-06}
160//  Coefficient values ± one standard deviation
161//      K0      =-571.42 ± 0.595
162//      K1      =1.1135 ± 0.00246
163//      K2      =-4.2444e-05 ± 2.15e-06
164//
165//
166//
167// for (8) tubes, keep all of the fit coefficients
168//
169//¥make/O/D/N=(3,8) fit_coef
170//¥edit fit_coef
171//¥make/O/D/N=(3,8) fit_sigma
172//¥edit fit_sigma
173//
174// -- copy and paste in the W_coef and W_sigma values (or by a command)
175//
176
177
178Proc MakeTableForFitCoefs(numTube,numCoef)
179        Variable numTube,numCoef
180       
181        Make/O/D/N=(numCoef,numTube) TubeCoefTable,TubeSigmaTable               //
182        Edit TubeCoefTable
183End
184
185Proc PlotFit_AllPeaks()
186
187        Variable ii,numTubes=8
188       
189        ii=1
190        do
191                PlotFit_Peaks(ii-1)
192                ii+=1
193        while(ii<=numTubes)
194
195End
196
197Proc PlotFit_Peaks(ind)
198        Variable ind
199       
200        //hopefully 20 points - need better control of this
201        Duplicate/O WA_PeakCentersX, tmpX
202       
203        tmpX = peakTableX[p][ind]
204        Display peak_spacing_mm_ctr vs tmpX
205       
206        CurveFit/M=2/W=0/TBOX=(0x310) poly 3, peak_spacing_mm_ctr/X=tmpX/D
207       
208        TubeCoefTable[][ind] = W_coef[p]
209        TubeSigmaTable[][ind] = W_sigma[p]
210       
211End
212
213
214
215
216
217//¥Duplicate tube1 tube1_mm
218//¥tube1_mm = V_TubePix_to_mm(fit_coef[0][0],fit_coef[1][0],fit_coef[2][0],p)
219
220
221
222
223////////
224// then there are various display options:
225
226// adjust the center (pixel) of the tube:
227// - measCtr is the pixel location of the DEFINED "zero" peak
228// nominal Ctr is the pixel number of this DEFINED zer0 position, nominally nPix/2-1
229//
230// ( be sure to pick better names, and use a loop...)
231//      adj_tube = raw_tube[p+(measCtr-nominalCtr)]
232//
233// then fill and display a new matrix. The center will be reasonably well aligned, and will
234// get worse towards the ends of the tubes
235// (this may be the "preferred" way of displaying the raw data if the centers are far off)
236// -- this also may be what I need to start with to fit the data to locate the beam center.
237
238
239// I can also display the fully corrected tubes, where the y-axis is now in real-space mm
240// rather than arbitrary pixels. The x-axis is still tube nubmer.
241// -- do this with the procedure"
242//   Append_Adjusted_mm()               // name may change...
243//
244// -- the point of the appending is that it allows each tube to be plotted on an image plot
245// with its own y-axis. Every one will be different and will be non-linear. These conditions
246// BOTH prevent using any of the "normal" image plotting or manipulation routines.
247// - the gist is below:
248//
249//      duplicate tube1_mm adjusted_mm_edge
250//      InsertPoints 0,1, adjusted_mm_edge
251//      // be sure to use the correct set of coefficients
252//      adjusted_mm_edge[0] = V_TubePix_to_mm(fit_coef[0][0],fit_coef[1][0],fit_coef[2][0],-1)
253//     
254//      Display;AppendImage adjusted_pack vs {*,adjusted_mm_edge}
255
256
257Proc Make_mm_tubes()
258
259        Variable ii,numTubes=8
260       
261        ii=1
262        do
263                Duplicate $("tube"+num2str(ii)) $("tube"+num2str(ii)+"_mm")
264                $("tube"+num2str(ii)+"_mm") = V_TubePix_to_mm(TubeCoefTable[0][ii-1],TubeCoefTable[1][ii-1],TubeCoefTable[2][ii-1],p)
265                ii+=1
266        while(ii<=numTubes)
267       
268End
269
270
271Proc Append_Adjusted_mm()
272
273// a blank base image
274        Duplicate/O pack junk
275        junk=0
276        SetScale/I y -600,600,"", junk          // -600,600 is tooo large and not general
277        Display;Appendimage junk
278
279        Variable ii
280       
281        ii=1
282        do
283                make/O/D/N=(1,1127) $("tube"+num2str(ii)+"_mm_mat")=0   
284       
285                $("tube"+num2str(ii)+"_mm_mat")[0][] = $("tube"+num2str(ii))
286                SetScale/I x (ii-1),ii,"", $("tube"+num2str(ii)+"_mm_mat")              //builds up the x-axis
287               
288                duplicate/O $("tube"+num2str(ii)+"_mm") $("edge"+num2str(ii)+"_mm")
289                InsertPoints 0,1, $("edge"+num2str(ii)+"_mm")           //needs to be one point longer
290        // be sure to use the correct set of coefficients
291                $("edge"+num2str(ii)+"_mm")[0] = V_TubePix_to_mm(TubeCoefTable[0][0],TubeCoefTable[1][0],TubeCoefTable[2][0],-1)
292       
293                AppendImage $("tube"+num2str(ii)+"_mm_mat") vs {*,$("edge"+num2str(ii)+"_mm")}
294                ModifyImage $("tube"+num2str(ii)+"_mm_mat") ctab= {*,*,ColdWarm,0}
295       
296                ii+=1
297        while(ii < 9)
298       
299end
300
301
302
303////////////////////////////////
304
305Proc MakeMatrix_PixCenter()
306        Duplicate/O pack pack_centered
307       
308        Variable ii,numTubes=8
309       
310        ii=1
311        do
312                FillMatrix_Pix_Center(ii)
313                ii+=1
314        while(ii<=numTubes)
315       
316        Display;AppendImage pack_centered
317        ModifyImage pack_centered ctab= {*,*,ColdWarm,0}
318
319end
320
321//
322// this fills a matrix with the tubes center aligned ONLY, with the y-axis in pixels
323//
324// adj_tube = raw_tube[p+(measCtr-nominalCtr)]
325// finds the center automatically from the tubeN_mm wave, where it crosses zero
326//
327// Tube #1 is the "base" ans others are shifted to match that one's "zero"
328//
329// FindRoots/P=W_coef           can also be used to find the roots (but which one?)
330//
331Function FillMatrix_Pix_Center(ind)
332        Variable ind
333       
334        Wave adj=root:pack_centered
335        Wave tube1_mm = $("root:tube1_mm")
336        Wave tube = $("root:tube"+num2str(ind))
337        wave tube_mm = $("root:tube"+num2str(ind)+"_mm")
338
339        Variable base,shift,ii,num,pt
340       
341        num=numpnts(tube)
342       
343        FindLevel tube1_mm 0
344        base = round(V_LevelX)
345       
346       
347        FindLevel tube_mm 0
348        shift = round(V_LevelX)
349       
350        for(ii=0;ii<num;ii+=1)
351                pt = ii + (shift-base)
352                if(pt >= 0 && pt < num)
353                        adj[ind-1][ii] = tube[pt]
354                endif
355        endfor
356       
357        return(0)
358End
359
360
361// this fills a matrix with the tubes center aligned ONLY, with the y-axis in mm
362// -- there seems to be little reason to do this --
363// -- either keep pixels and align centers
364// -- OR -- use mm and append each tube with its own y-axis
365//
366Function FillAdjusted(ind)
367        Variable ind
368       
369        Wave adj=root:adjusted_pack
370        Wave tube1_mm
371        Wave tube = $("root:tube"+num2str(ind))
372        wave tube_mm = $("root:tube"+num2str(ind)+"_mm")
373
374        Variable base,shift,ii,num,pt
375       
376        num=numpnts(tube1_mm)
377       
378        FindLevel tube1_mm 0
379        base = round(V_LevelX)
380       
381       
382        FindLevel tube_mm 0
383        shift = round(V_LevelX)
384       
385        for(ii=0;ii<num;ii+=1)
386                pt = ii + (shift-base)
387                if(pt >= 0 && pt < num)
388                        adj[ind-1][ii] = tube[pt]
389                endif
390        endfor
391       
392        return(0)
393End
394
395
396// c0,c1,c2,pix
397// c0-c2 are the fit coefficients
398// pix is the test pixel
399//
400// returns the distance in mm (relative to ctr pixel)
401// ctr is the center pixel, as defined when fitting to quadratic
402//
403Function V_TubePix_to_mm(c0,c1,c2,pix)
404        Variable c0,c1,c2,pix
405       
406        Variable dist
407        dist = c0 + c1*pix + c2*pix*pix
408       
409        return(dist)
410End
411
412////
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434// this doesn't work - the interpolation step doesn't do what I want
435// and the plot of the triplet with f(z) for color doesn't fill space like I want
436Proc AnotherExample()
437
438        Concatenate/O/NP {tube1_mm,tube2_mm,tube3_mm,tube4_mm,tube5_mm,tube6_mm,tube7_mm,tube8_mm},cat_mm
439        Concatenate/O/NP {tube1,tube2,tube3,tube4,tube5,tube6,tube7,tube8},cat_tubes
440        Duplicate/O cat_mm,cat_num
441        Variable num=1127
442        cat_num[0,num-1]=1
443        cat_num[num,2*num-1]=2
444        cat_num[2*num,3*num-1]=3
445        cat_num[3*num,4*num-1]=4
446        cat_num[4*num,5*num-1]=5
447        cat_num[5*num,6*num-1]=6
448        cat_num[6*num,7*num-1]=7
449        cat_num[7*num,8*num-1]=8
450
451        Display cat_mm vs cat_num
452        ModifyGraph mode=3,marker=9
453        ModifyGraph zColor(cat_mm)={cat_tubes,*,*,ColdWarm,0}
454
455        Concatenate/O {cat_num,cat_mm,cat_tubes}, tripletWave
456        ImageInterpolate Kriging tripletWave
457        AppendImage M_InterpolatedImage
458
459//      Make/O/N=20 xWave=enoise(4),yWave=enoise(5),zWave=enoise(6)  // Random points
460//      Display yWave vs xWave
461//      ModifyGraph mode=3,marker=19
462//      ModifyGraph zColor(yWave)={zWave,*,*,Rainbow,0}
463//
464//      Concatenate/O {xWave,yWave,zWave}, tripletWave
465//      ImageInterpolate/S={-5,0.1,5,-5,0.1,5} voronoi tripletWave
466//      AppendImage M_InterpolatedImage
467
468end
469
470// this desn't work either...
471// (same y-axis for the entire image, which is not the case for the tubes)
472//
473// from the WM help file:
474// Plotting a 2D Z Wave With 1D X and Y Center Data
475//
476Function MakeEdgesWave(centers, edgesWave)
477        Wave centers    // Input
478        Wave edgesWave  // Receives output
479       
480        Variable N=numpnts(centers)
481        Redimension/N=(N+1) edgesWave
482
483        edgesWave[0]=centers[0]-0.5*(centers[1]-centers[0])
484        edgesWave[N]=centers[N-1]+0.5*(centers[N-1]-centers[N-2])
485        edgesWave[1,N-1]=centers[p]-0.5*(centers[p]-centers[p-1])
486End
487
488//This function demonstrates the use of MakeEdgesWave:
489Function DemoPlotXYZAsImage()
490        Make/O mat={{0,1,2},{2,3,4},{3,4,5}}    // Matrix containing Z values
491        Make/O centersX = {1, 2.5, 5}           // X centers wave
492        Make/O centersY = {300, 400, 600}               // Y centers wave
493        Make/O edgesX; MakeEdgesWave(centersX, edgesX)  // Create X edges wave
494        Make/O edgesY; MakeEdgesWave(centersY, edgesY)  // Create Y edges wave
495        Display; AppendImage mat vs {edgesX,edgesY}
496End
497
Note: See TracBrowser for help on using the repository browser.