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

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

event mode processing - splitting the binned data into 4 panels, speedup of binning, display of split panels.

fitting routines for fitting 5 peaks for nonlinear corrections.

all needs to documented, then expanded

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