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

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

cleaning up a lot of the TODO items from the code.

File size: 30.4 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//
14//
15// x- need a way to generate the known, physical dimensions of the slots
16// Make/O/D/N=5 peak_spacing_mm_ctr
17// peak_spacing_mm_ctr = {-350,-190,0,190,350} (to be filled in with the correct measurements,
18//   possibly different for each panel)
19//
20// x- a 128 point wave "tube_pixel" (=p) is made in V_ArrayToTubes(), and is needed for the WM
21//   procedures to identify the peak positions.
22//
23// x- fit with either gauss or lor function to get non-integer pixel values for the peak locations
24//
25// x- do I fit each individually to "tweak" the located values, or fit all 5 at once with a
26//    custom fit function and guess some good starting values for peak height, location, etc.
27//
28//
29// x- find a way to display all of the results - in a way that can quickly identify any fits
30//    that may be incorrect
31//
32//
33
34
35
36// the main routines are:
37
38//(1)
39//
40// to get from an array to individual tubes
41// V_ArrayToTubes(detector)
42//
43// (not needed) to get from individual tubes to an array
44//      V_Tubes_to_Array()                     
45
46//(2)
47// then to locate all of the peak positions
48//      V_MakeTableForPeaks(numTube,numPeak)           
49//      V_Identify_AllPeaks()
50//              AutoFindPeaksCustom()           // if Identify_AllPeaks  doesn't work -try this, setting the "noise" to 1 and smoothing to 2
51
52// (3) Refine the fitted peak positions
53//
54
55//(4)
56// fit to find all of the quadratic coefficients
57//      MakeTableForFitCoefs(numTube,numCoef)
58//      PlotFit_AllPeaks()
59
60
61//(5)
62// then pick a display method
63//
64//      Make_mm_tubes()
65//      Append_Adjusted_mm()
66//
67//      MakeMatrix_PixCenter()
68//      FillMatrix_Pix_Center(ind)
69//
70//
71// -or- (note that the pack_image wave that is generated here is for display ONLY)
72// --since it is interpolated
73//
74// Interpolate_mm_tubes()
75
76
77// The function most used externally is:
78// V_TubePix_to_mm(c0,c1,c2,pix)
79//
80// which will return the real space distance (in mm?) for a given pixel
81// and the tube's coefficients. The distance is relative to the zero position on the
82// detector (which was defined when the coefficients were determined)
83
84
85
86//
87// (0) -- what I start with:
88// -- a table of the mm spacing of the slots (20 of them)
89// -- masked data from each of the (8) tubes
90// -- the table of slots may need to be corrected for parallax, depending on the geometry of the test
91// ** In the table of slots, pick a slot near the center, and SET that to ZERO. Then all of the other
92//   distances are relative to that zero point. This is a necessary reference point.
93//
94
95//
96// x- need a routine to set up the actual measurements of the slot positions
97//
98//
99//
100// x- the slot positioning is different for the L/R and T/B detectors
101//
102Proc V_SetupSlotDimensions()
103        Make/O/D/N=5 peak_spacing_mm_ctr_TB,peak_spacing_mm_ctr_LR
104        peak_spacing_mm_ctr_TB = {-159.54,-80.17,0,80.17,159.54}
105        peak_spacing_mm_ctr_LR = {-379.4,-189.7,0,189.7,380.2}
106        DoWindow/F Real_mm_Table
107        if(V_Flag == 0)
108                Edit/N=Real_mm_Table peak_spacing_mm_ctr_TB,peak_spacing_mm_ctr_LR
109        endif
110End
111
112
113
114//
115// (1) -- get the individual tubes into an array
116//
117//
118Proc V_Tubes_to_Array()
119        Make/O/D/N=(8,1127) pack
120        edit pack
121        display;appendimage pack
122        pack[0][] = tube1[q]
123        pack[1][] = tube2[q]
124        pack[2][] = tube3[q]
125        pack[3][] = tube4[q]
126        pack[4][] = tube5[q]
127        pack[5][] = tube6[q]
128        pack[6][] = tube7[q]
129        pack[7][] = tube8[q]
130        ModifyImage pack ctab= {*,*,ColdWarm,0}
131End
132
133// or the other way around
134Proc V_ArrayToTubes(wStr)
135        String wStr
136       
137        String/G root:detUsed = wStr
138       
139        Variable ii,numTubes=48
140        String str="tube"
141       
142        Variable dim0,dim1
143        dim0 = DimSize($wStr,0)
144        dim1 = DimSize($wStr,1)
145
146       
147        Make/O/D/N=128 tube_pixel
148        tube_pixel = p
149       
150       
151        ii=0
152        do
153                Make/O/D/N=128 $(str+num2str(ii))
154               
155                if(dim0 == 128)
156                        $(str+num2str(ii)) = $(wStr)[p][ii]
157                else
158                        $(str+num2str(ii)) = $(wStr)[ii][p]
159                endif
160               
161                ii+=1
162        while(ii < numTubes)
163
164End
165
166
167// (2) -- for each of the tubes, find the x-position (in pixels) of each of the (20) peaks
168// -- load the Analysis Package "MultiPeakFit 2"
169//
170// automatically find the peaks (after including MultiPeakFit 2)
171//              AutomaticallyFindPeaks()
172//
173//-- or if having difficulty
174//              AutoFindPeaksCustom()           // try this, setting the "noise" to 1 and smoothing to 2
175//
176// -- if really having difficulty, you can do the "full" MultiPeak Fit
177//
178// -- If (hopefully) using the easy way, the results are in:
179// root:WA_PeakCentersY,root:WA_PeakCentersX
180//
181// -- so to see the results:
182//¥Edit/K=0  root:WA_PeakCentersY,root:WA_PeakCentersX
183//
184// -- then sort the results - they seem to be in no real order...
185//¥Sort WA_PeakCentersX WA_PeakCentersY,WA_PeakCentersX
186//
187Proc V_MakeTableForPeaks(numTube,numPeak)
188        Variable numTube=48,numPeak=5
189       
190        Make/O/D/N=(numPeak,numTube) PeakTableX,PeakTableY              //*2 to store x-location and peak height (y)
191       
192        DoWindow/F Peak_Pixel_Loc
193        if(V_flag == 0)
194                Edit/N=Peak_Pixel_Loc peakTableX
195        endif
196        DoAlert 0, "Load the Package: Analysis->MultiPeak Fitting->MultiPeak Fitting 2"
197End
198
199Proc V_Identify_AllPeaks()
200
201        Variable ii,numTubes=48
202        String str="tube"
203       
204        ii=0
205        do
206                V_Identify_Peaks(str+num2str(ii),ii)
207                ii+=1
208        while(ii < numTubes)
209
210End
211
212Proc V_Identify_Peaks(tubeStr,ind)
213        String tubeStr
214        Variable ind
215       
216        // must use a wave of pixels rather than "calculated" -- if calculated is used it only
217        // returns integer values for the peak locations
218       
219//      AutomaticallyFindPeaks() //-- this is a function that doesn't take any parameters - so
220// I need to pull this from the WM function to call the worker directly
221        Variable pBegin=0, pEnd= numpnts($(tubeStr))-1
222        Variable/C estimates= EstPeakNoiseAndSmfact($(tubeStr),pBegin, pEnd)
223        Variable noiselevel=real(estimates)
224        Variable smoothingFactor=imag(estimates)
225        Variable maxPeaks = 20
226        Variable minPeakPercent = 10
227       
228        AutoFindPeaksWorker($(tubeStr), $("tube_pixel"), pBegin, pEnd, maxPeaks, minPeakPercent, noiseLevel, smoothingFactor)
229// end WM function call
230
231        Sort WA_PeakCentersX WA_PeakCentersY,WA_PeakCentersX
232       
233        peakTableX[][ind] = WA_PeakCentersX[p]          // the peak position
234        peakTableY[][ind] = WA_PeakCentersY[p]          // the peak height
235       
236End
237
238
239
240
241// ADD
242// a step to refine the peak positioning - currently an integer value
243//  fit with a gauss or lorentzian
244
245// CurveFit/M=2/W=0/TBOX=(0x310) lor, tube47[29,53]/X=tube_pixel[29,53]/D
246
247//CurveFit/M=2/W=0 lor, tube47[29,53]/X=tube_pixel[29,53]/D
248//fit_tube47= W_coef[0]+W_coef[1]/((x-W_coef[2])^2+W_coef[3])
249//W_coef={-20.37,876.94,40.078,0.5201}
250//W_sigma={6.52,47.3,0.0241,0.0308}
251
252Proc V_MakeTableForRefinedFit(numTube,numPeak)
253        Variable numTube=48,numPeak=5
254       
255        Make/O/D/N=(numPeak,numTube) position_refined,position_refined_err              //
256       
257        DoWindow/F Refined_Positions
258        if(V_flag == 0)
259                Edit/N=Refined_Positions position_refined
260        endif
261       
262End
263
264Proc V_Refine_All_PeakPos()
265
266        Variable ii,numTubes=48
267       
268        ii=0
269        do
270                V_Refine_PeakPos(ii)
271                ii+=1
272        while(ii<numTubes)
273
274End
275
276
277//CurveFit/M=2/W=0 lor, tube47[29,53]/X=tube_pixel[29,53]/D
278//fit_tube47= W_coef[0]+W_coef[1]/((x-W_coef[2])^2+W_coef[3])
279
280Proc V_Refine_PeakPos(ind)
281        Variable ind
282       
283//
284// x- hard-wired for 5 peaks
285
286        Variable ii,lo,hi
287       
288       
289        ii=0
290        do
291       
292                if(ii==0)
293                // 1st peak
294                // define fitting range pixels (integer)
295                        lo = 0
296                else
297                        lo = trunc(0.5*(peakTableX[ii-1][ind] + peakTableX[ii][ind]))
298                endif
299               
300                if(ii==4)
301                        hi = numpnts(tube_pixel)-1
302                else
303                        hi = trunc(0.5*(peakTableX[ii][ind] + peakTableX[ii+1][ind]))
304                endif
305               
306                // do I need initial guesses?
307                CurveFit/M=0/W=2 lor, $("tube"+num2str(ind))[lo,hi]/X=tube_pixel[lo,hi]/D
308               
309                position_refined[ii][ind] = W_coef[2]
310                position_refined_err[ii][ind] = W_sigma[2]
311
312                ii += 1
313
314        while(ii < 5)
315       
316End
317
318
319
320
321
322// -- save a copy of the root:WA_PeakCentersY,root:WA_PeakCentersX values
323//    for later in case the fitting failed, then you can go back and re-do
324//
325// -- then plot:
326//
327//      Display peak_spacing_mm_ctr vs WA_PeakCentersX
328//
329//      Then do a "QuickFit" of the peak position to the data using a polynomial of order 3 (= quadratic)
330//
331// result is in W_coef, W_sigma
332//
333// -- an example of the "quickFit" command is below, so it can be programmed rather than the menu every time
334//¥Display peak_spacing_mm_ctr vs WA_PeakCentersX3
335//¥CurveFit/M=2/W=0/TBOX=(0x310) poly 3, peak_spacing_mm_ctr/X=WA_PeakCentersX3/D
336//  fit_peak_spacing_mm_ctr= poly(W_coef,x)
337//  W_coef={-571.42,1.1135,-4.2444e-05}
338//  V_chisq= 8.5841;V_npnts= 20;V_numNaNs= 0;V_numINFs= 0;
339//  V_startRow= 0;V_endRow= 19;
340//  W_sigma={0.595,0.00246,2.15e-06}
341//  Coefficient values ± one standard deviation
342//      K0      =-571.42 ± 0.595
343//      K1      =1.1135 ± 0.00246
344//      K2      =-4.2444e-05 ± 2.15e-06
345//
346//
347//
348// for (8) tubes, keep all of the fit coefficients
349//
350//¥make/O/D/N=(3,8) fit_coef
351//¥edit fit_coef
352//¥make/O/D/N=(3,8) fit_sigma
353//¥edit fit_sigma
354//
355// -- copy and paste in the W_coef and W_sigma values (or by a command)
356//
357
358
359Proc V_MakeTableForFitCoefs(numTube,numCoef)
360        Variable numTube=48,numCoef=3
361       
362        Make/O/D/N=(numTube,numCoef) TubeCoefTable,TubeSigmaTable               //
363       
364        DoWindow/F Quad_Coefficients
365        if(V_flag == 0)
366                Edit/N=Quad_Coefficients TubeCoefTable
367        endif
368
369        String detUsed = root:detUsed
370       
371        if(strsearch(detUsed,"L",0) >= 0 || strsearch(detUsed,"R",0) >= 0)
372                Duplicate/O     peak_spacing_mm_ctr_LR, peak_spacing_mm_ctr
373                DoAlert 0,"Using peak spacing for L/R"
374        else
375                Duplicate/O     peak_spacing_mm_ctr_TB, peak_spacing_mm_ctr
376                DoAlert 0,"Using peak spacing for T/B"
377        endif
378End
379
380Proc V_PlotFit_AllPeakPosition()
381
382        Variable ii,numTubes=48
383       
384        ii=0
385        do
386                V_PlotFit_PeakPosition(ii)
387                ii+=1
388        while(ii<numTubes)
389
390End
391
392Proc V_PlotFit_PeakPosition(ind)
393        Variable ind
394       
395        Duplicate/O WA_PeakCentersX, tmpX
396       
397//      tmpX = peakTableX[p][ind]
398        tmpX = position_refined[p][ind]
399//      Display peak_spacing_mm_ctr vs tmpX
400       
401//      CurveFit/M=2/W=0/TBOX=(0x310) poly 3, peak_spacing_mm_ctr/X=tmpX/D
402        CurveFit/M=0/W=2 poly 3, peak_spacing_mm_ctr/X=tmpX/D
403       
404        TubeCoefTable[ind][] = W_coef[q]
405        TubeSigmaTable[ind][] = W_sigma[q]
406       
407End
408
409
410
411
412//¥Duplicate tube1 tube1_mm
413//¥tube1_mm = V_TubePix_to_mm(fit_coef[0][0],fit_coef[1][0],fit_coef[2][0],p)
414
415
416////////
417// then there are various display options:
418
419// adjust the center (pixel) of the tube:
420// - measCtr is the pixel location of the DEFINED "zero" peak
421// nominal Ctr is the pixel number of this DEFINED zer0 position, nominally nPix/2-1
422//
423// ( be sure to pick better names, and use a loop...)
424//      adj_tube = raw_tube[p+(measCtr-nominalCtr)]
425//
426// then fill and display a new matrix. The center will be reasonably well aligned, and will
427// get worse towards the ends of the tubes
428// (this may be the "preferred" way of displaying the raw data if the centers are far off)
429// -- this also may be what I need to start with to fit the data to locate the beam center.
430
431
432// I can also display the fully corrected tubes, where the y-axis is now in real-space mm
433// rather than arbitrary pixels. The x-axis is still tube nubmer.
434// -- do this with the procedure"
435//   Append_Adjusted_mm()               // name may change...
436//
437// -- the point of the appending is that it allows each tube to be plotted on an image plot
438// with its own y-axis. Every one will be different and will be non-linear. These conditions
439// BOTH prevent using any of the "normal" image plotting or manipulation routines.
440// - the gist is below:
441//
442//      duplicate tube1_mm adjusted_mm_edge
443//      InsertPoints 0,1, adjusted_mm_edge
444//      // be sure to use the correct set of coefficients
445//      adjusted_mm_edge[0] = V_TubePix_to_mm(fit_coef[0][0],fit_coef[1][0],fit_coef[2][0],-1)
446//     
447//      Display;AppendImage adjusted_pack vs {*,adjusted_mm_edge}
448
449
450Proc V_Make_mm_tubes()
451
452        Variable ii,numTubes=8
453       
454        ii=1
455        do
456                Duplicate $("tube"+num2str(ii)) $("tube"+num2str(ii)+"_mm")
457                $("tube"+num2str(ii)+"_mm") = V_TubePix_to_mm(TubeCoefTable[ii-1][0],TubeCoefTable[ii-1][1],TubeCoefTable[ii-1][2],p)
458                ii+=1
459        while(ii<=numTubes)
460       
461End
462
463
464Proc V_Append_Adjusted_mm()
465
466// a blank base image
467        Duplicate/O pack junk
468        junk=0
469        SetScale/I y -600,600,"", junk          // -600,600 is tooo large and not general
470        Display;Appendimage junk
471
472        Variable ii
473       
474        ii=1
475        do
476                make/O/D/N=(1,1127) $("tube"+num2str(ii)+"_mm_mat")=0   
477       
478                $("tube"+num2str(ii)+"_mm_mat")[0][] = $("tube"+num2str(ii))
479                SetScale/I x (ii-1),ii,"", $("tube"+num2str(ii)+"_mm_mat")              //builds up the x-axis
480               
481                duplicate/O $("tube"+num2str(ii)+"_mm") $("edge"+num2str(ii)+"_mm")
482                InsertPoints 0,1, $("edge"+num2str(ii)+"_mm")           //needs to be one point longer
483        // be sure to use the correct set of coefficients
484                $("edge"+num2str(ii)+"_mm")[0] = V_TubePix_to_mm(TubeCoefTable[0][0],TubeCoefTable[0][1],TubeCoefTable[0][2],-1)
485       
486                AppendImage $("tube"+num2str(ii)+"_mm_mat") vs {*,$("edge"+num2str(ii)+"_mm")}
487                ModifyImage $("tube"+num2str(ii)+"_mm_mat") ctab= {*,*,ColdWarm,0}
488       
489                ii+=1
490        while(ii < 9)
491       
492end
493
494
495
496////////////////////////////////
497
498Proc V_MakeMatrix_PixCenter()
499        Duplicate/O pack pack_centered
500       
501        Variable ii,numTubes=8
502       
503        ii=1
504        do
505                V_FillMatrix_Pix_Center(ii)
506                ii+=1
507        while(ii<=numTubes)
508       
509        Display;AppendImage pack_centered
510        ModifyImage pack_centered ctab= {*,*,ColdWarm,0}
511
512end
513
514//
515// this fills a matrix with the tubes center aligned ONLY, with the y-axis in pixels
516//
517// adj_tube = raw_tube[p+(measCtr-nominalCtr)]
518// finds the center automatically from the tubeN_mm wave, where it crosses zero
519//
520// Tube #1 is the "base" ans others are shifted to match that one's "zero"
521//
522// FindRoots/P=W_coef           can also be used to find the roots (but which one?)
523//
524Function V_FillMatrix_Pix_Center(ind)
525        Variable ind
526       
527        Wave adj=root:pack_centered
528        Wave tube1_mm = $("root:tube1_mm")
529        Wave tube = $("root:tube"+num2str(ind))
530        wave tube_mm = $("root:tube"+num2str(ind)+"_mm")
531
532        Variable base,shift,ii,num,pt
533       
534        num=numpnts(tube)
535       
536        FindLevel tube1_mm 0
537        base = round(V_LevelX)
538       
539       
540        FindLevel tube_mm 0
541        shift = round(V_LevelX)
542       
543        for(ii=0;ii<num;ii+=1)
544                pt = ii + (shift-base)
545                if(pt >= 0 && pt < num)
546                        adj[ind-1][ii] = tube[pt]
547                endif
548        endfor
549       
550        return(0)
551End
552
553
554// this fills a matrix with the tubes center aligned ONLY, with the y-axis in mm
555// -- there seems to be little reason to do this --
556// -- either keep pixels and align centers
557// -- OR -- use mm and append each tube with its own y-axis
558//
559Function V_FillAdjusted(ind)
560        Variable ind
561       
562        Wave adj=root:adjusted_pack
563        Wave tube1_mm
564        Wave tube = $("root:tube"+num2str(ind))
565        wave tube_mm = $("root:tube"+num2str(ind)+"_mm")
566
567        Variable base,shift,ii,num,pt
568       
569        num=numpnts(tube1_mm)
570       
571        FindLevel tube1_mm 0
572        base = round(V_LevelX)
573       
574       
575        FindLevel tube_mm 0
576        shift = round(V_LevelX)
577       
578        for(ii=0;ii<num;ii+=1)
579                pt = ii + (shift-base)
580                if(pt >= 0 && pt < num)
581                        adj[ind-1][ii] = tube[pt]
582                endif
583        endfor
584       
585        return(0)
586End
587
588
589// c0,c1,c2,pix
590// c0-c2 are the fit coefficients
591// pix is the test pixel
592//
593// returns the distance in mm (relative to ctr pixel)
594// ctr is the center pixel, as defined when fitting to quadratic
595//
596Function V_TubePix_to_mm(c0,c1,c2,pix)
597        Variable c0,c1,c2,pix
598       
599        Variable dist
600        dist = c0 + c1*pix + c2*pix*pix
601       
602        return(dist)
603End
604
605////
606
607
608
609// set the (linear) range of the y-axis of the matrix to be the
610// range of the 1st tube. This is completely arbitrary
611//
612Proc V_Interpolate_mm_tubes()
613
614        Duplicate/O pack pack_image
615
616        Variable ii,numTubes=8
617        Variable p1,p2
618        p1 = tube1_mm[0]
619        p2 = tube1_mm[numpnts(tube1_mm)-1]
620       
621        SetScale/I y p1,p2,"", pack_image
622       
623        // then make a temporary 1D wave to help with the interpolation
624        Duplicate/O tube1_mm lin_mm,lin_val
625        SetScale/I x p1,p2,"", lin_mm
626        lin_mm = x                      //fill with the linear mm spacing
627        lin_val=0
628       
629        ii=1
630        do
631                lin_val = interp(lin_mm, $("tube"+num2str(ii)+"_mm"), $("tube"+num2str(ii)))
632                pack_image[ii-1][] = lin_val[q]
633               
634                ii+=1
635        while(ii<=numTubes)
636       
637        display;appendimage pack_image
638        ModifyImage pack_image ctab= {*,*,ColdWarm,0}
639       
640End
641
642
643
644
645
646
647
648
649
650
651
652
653
654// this doesn't work - the interpolation step doesn't do what I want
655// and the plot of the triplet with f(z) for color doesn't fill space like I want
656Proc V_AnotherExample()
657
658        Concatenate/O/NP {tube1_mm,tube2_mm,tube3_mm,tube4_mm,tube5_mm,tube6_mm,tube7_mm,tube8_mm},cat_mm
659        Concatenate/O/NP {tube1,tube2,tube3,tube4,tube5,tube6,tube7,tube8},cat_tubes
660        Duplicate/O cat_mm,cat_num
661        Variable num=1127
662        cat_num[0,num-1]=1
663        cat_num[num,2*num-1]=2
664        cat_num[2*num,3*num-1]=3
665        cat_num[3*num,4*num-1]=4
666        cat_num[4*num,5*num-1]=5
667        cat_num[5*num,6*num-1]=6
668        cat_num[6*num,7*num-1]=7
669        cat_num[7*num,8*num-1]=8
670
671        Display cat_mm vs cat_num
672        ModifyGraph mode=3,marker=9
673        ModifyGraph zColor(cat_mm)={cat_tubes,*,*,ColdWarm,0}
674
675        Concatenate/O {cat_num,cat_mm,cat_tubes}, tripletWave
676        ImageInterpolate Kriging tripletWave
677        AppendImage M_InterpolatedImage
678
679//      Make/O/N=20 xWave=enoise(4),yWave=enoise(5),zWave=enoise(6)  // Random points
680//      Display yWave vs xWave
681//      ModifyGraph mode=3,marker=19
682//      ModifyGraph zColor(yWave)={zWave,*,*,Rainbow,0}
683//
684//      Concatenate/O {xWave,yWave,zWave}, tripletWave
685//      ImageInterpolate/S={-5,0.1,5,-5,0.1,5} voronoi tripletWave
686//      AppendImage M_InterpolatedImage
687
688end
689
690// this desn't work either...
691// (same y-axis for the entire image, which is not the case for the tubes)
692//
693// from the WM help file:
694// Plotting a 2D Z Wave With 1D X and Y Center Data
695//
696Function V_MakeEdgesWave(centers, edgesWave)
697        Wave centers    // Input
698        Wave edgesWave  // Receives output
699       
700        Variable N=numpnts(centers)
701        Redimension/N=(N+1) edgesWave
702
703        edgesWave[0]=centers[0]-0.5*(centers[1]-centers[0])
704        edgesWave[N]=centers[N-1]+0.5*(centers[N-1]-centers[N-2])
705        edgesWave[1,N-1]=centers[p]-0.5*(centers[p]-centers[p-1])
706End
707
708//This function demonstrates the use of MakeEdgesWave:
709Function V_DemoPlotXYZAsImage()
710        Make/O mat={{0,1,2},{2,3,4},{3,4,5}}    // Matrix containing Z values
711        Make/O centersX = {1, 2.5, 5}           // X centers wave
712        Make/O centersY = {300, 400, 600}               // Y centers wave
713        Make/O edgesX; V_MakeEdgesWave(centersX, edgesX)        // Create X edges wave
714        Make/O edgesY; V_MakeEdgesWave(centersY, edgesY)        // Create Y edges wave
715        Display; AppendImage mat vs {edgesX,edgesY}
716End
717
718
719
720////////////////////////////
721//
722// Main entry - open the panel and go through
723// each of the steps for each of the detector panels
724//
725Proc V_TubeCoefPanel() : Panel
726        PauseUpdate; Silent 1           // building window...
727        NewPanel /W=(973,45,1156,535)/K=1
728        DoWindow/C V_TubeCoefPanel
729//      ShowTools/A
730
731        SetDrawLayer UserBack
732        SetDrawEnv fsize= 14,fstyle= 1
733        DrawText 5,58,"(1)"
734        SetDrawEnv fsize= 14,fstyle= 1
735        DrawText 5,108,"(2)"
736        SetDrawEnv fsize= 14,fstyle= 1
737        DrawText 5,158,"(3)"
738        SetDrawEnv fsize= 14,fstyle= 1
739        DrawText 5,208,"(4)"
740        SetDrawEnv fsize= 14,fstyle= 1
741        DrawText 5,258,"(5)"
742        SetDrawEnv fsize= 14,fstyle= 1
743        DrawText 5,308,"(6)"
744        SetDrawEnv fsize= 14,fstyle= 1
745        DrawText 5,358,"(7)"
746        SetDrawEnv fsize= 14,fstyle= 1
747        DrawText 5,408,"(8)"
748        SetDrawEnv fsize= 14,fstyle= 1
749        DrawText 5,458,"(9)"
750                       
751        Button button_0,pos={30.00,40.00},size={120.00,20.00},proc=V_Setup_MasksButton,title="Setup"
752        Button button_1,pos={30.00,90.00},size={120.00,20.00},proc=V_ArrayToTubesButton,title="Array to Tubes"
753        Button button_2,pos={30.00,140.00},size={120.00,20.00},proc=V_TableForPeaksButton,title="Table for Peaks"
754        Button button_3,pos={30.00,190.00},size={120.00,20.00},proc=V_IdentifyPeaksButton,title="Identify Peaks"
755        Button button_4,pos={30.00,240.00},size={120.00,20.00},proc=V_RefineTableButton,title="Refine Peak Table"
756        Button button_5,pos={30.00,290.00},size={120.00,20.00},proc=V_RefinePeaksButton,title="Refine Peaks"
757
758        Button button_6,pos={30.00,340.00},size={120.00,20.00},proc=V_QuadFitTableButton,title="Table for Quad"
759        Button button_7,pos={30.00,390.00},size={120.00,20.00},proc=V_QuadFitButton,title="Fit to Quad"
760        Button button_8,pos={30.00,440},size={120.00,20.00},proc=V_PeakPlotButton,title="Plot Peaks"
761       
762EndMacro
763
764
765// a simple display of the refined results
766//
767Function V_PeakPlotButton(ba) : ButtonControl
768        STRUCT WMButtonAction &ba
769
770        switch( ba.eventCode )
771                case 2: // mouse up
772                        // click code here
773                        Execute "V_OpenPeakResultsGraph()"
774                        break
775                case -1: // control being killed
776                        break
777        endswitch
778
779        return 0
780End
781
782// generate the waves and the table
783//
784Function V_TableForPeaksButton(ba) : ButtonControl
785        STRUCT WMButtonAction &ba
786
787        switch( ba.eventCode )
788                case 2: // mouse up
789                        // click code here
790                        Execute "V_MakeTableForPeaks()"
791                        break
792                case -1: // control being killed
793                        break
794        endswitch
795
796        return 0
797End
798
799// use the WM procedures to quickly identify the peak position (and height)
800// to be used in the refining fits
801//
802Function V_IdentifyPeaksButton(ba) : ButtonControl
803        STRUCT WMButtonAction &ba
804
805        switch( ba.eventCode )
806                case 2: // mouse up
807                        // click code here
808                        Execute "V_Identify_AllPeaks()"
809                        break
810                case -1: // control being killed
811                        break
812        endswitch
813
814        return 0
815End
816
817// generate the waves and the table
818//
819Function V_RefineTableButton(ba) : ButtonControl
820        STRUCT WMButtonAction &ba
821
822        switch( ba.eventCode )
823                case 2: // mouse up
824                        // click code here
825                        Execute "V_MakeTableForRefinedFit()"
826                        break
827                case -1: // control being killed
828                        break
829        endswitch
830
831        return 0
832End
833
834// using the initial peak locations from WM, refine the values
835// by fitting each individual peak
836//
837Function V_RefinePeaksButton(ba) : ButtonControl
838        STRUCT WMButtonAction &ba
839
840        switch( ba.eventCode )
841                case 2: // mouse up
842                        // click code here
843                        Execute "V_Refine_All_PeakPos()"
844                        break
845                case -1: // control being killed
846                        break
847        endswitch
848
849        return 0
850End
851
852// finally, with the peak positions, make waves and a table for the
853// quadratic coefficients
854//
855Function V_QuadFitTableButton(ba) : ButtonControl
856        STRUCT WMButtonAction &ba
857
858        switch( ba.eventCode )
859                case 2: // mouse up
860                        // click code here
861                        Execute "V_MakeTableForFitCoefs()"
862                        break
863                case -1: // control being killed
864                        break
865        endswitch
866
867        return 0
868End
869
870// fit all of the peak positions per tube vs. the actual mm locations to
871// obtain the quadratic coefficients
872//
873Function V_QuadFitButton(ba) : ButtonControl
874        STRUCT WMButtonAction &ba
875
876        switch( ba.eventCode )
877                case 2: // mouse up
878                        // click code here
879                        Execute "V_PlotFit_AllPeakPosition()"
880                        break
881                case -1: // control being killed
882                        break
883        endswitch
884
885        return 0
886End
887
888
889// fill the waves and table with the hard-wired slot positions (mm)
890//
891Function V_Setup_MasksButton(ba) : ButtonControl
892        STRUCT WMButtonAction &ba
893
894        switch( ba.eventCode )
895                case 2: // mouse up
896                        // click code here
897                        Execute "V_SetupSlotDimensions()"
898                        break
899                case -1: // control being killed
900                        break
901        endswitch
902
903        return 0
904End
905
906
907// convert the named detector array to 48 individual tube waves
908//
909Function V_ArrayToTubesButton(ba) : ButtonControl
910        STRUCT WMButtonAction &ba
911
912        switch( ba.eventCode )
913                case 2: // mouse up
914                        // click code here
915                        Execute "V_ArrayToTubes()"
916                        break
917                case -1: // control being killed
918                        break
919        endswitch
920
921        return 0
922End
923
924
925///////////////////////////
926//
927// unused - a simple line graph for each tube is much simpler
928//
929Window Gizmo_refinedPositions() : GizmoPlot
930        PauseUpdate; Silent 1           // building window...
931        // Building Gizmo 7 window...
932        NewGizmo/W=(232,448,747,908)
933        ModifyGizmo startRecMacro=700
934        ModifyGizmo scalingOption=63
935        AppendToGizmo Surface=root:position_refined,name=surface0
936        ModifyGizmo ModifyObject=surface0,objectType=surface,property={ srcMode,0}
937        ModifyGizmo ModifyObject=surface0,objectType=surface,property={ surfaceCTab,Rainbow}
938        AppendToGizmo Axes=boxAxes,name=axes0
939        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={0,axisRange,-1,-1,-1,1,-1,-1}
940        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={1,axisRange,-1,-1,-1,-1,1,-1}
941        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={2,axisRange,-1,-1,-1,-1,-1,1}
942        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={3,axisRange,-1,1,-1,-1,1,1}
943        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={4,axisRange,1,1,-1,1,1,1}
944        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={5,axisRange,1,-1,-1,1,-1,1}
945        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={6,axisRange,-1,-1,1,-1,1,1}
946        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={7,axisRange,1,-1,1,1,1,1}
947        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={8,axisRange,1,-1,-1,1,1,-1}
948        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={9,axisRange,-1,1,-1,1,1,-1}
949        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={10,axisRange,-1,1,1,1,1,1}
950        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={11,axisRange,-1,-1,1,1,-1,1}
951        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={-1,axisScalingMode,1}
952        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={-1,axisColor,0,0,0,1}
953        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={0,ticks,2}
954        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={1,ticks,2}
955        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={2,ticks,2}
956        ModifyGizmo modifyObject=axes0,objectType=Axes,property={-1,Clipped,0}
957        ModifyGizmo setDisplayList=0, object=surface0
958        ModifyGizmo setDisplayList=1, object=axes0
959        ModifyGizmo autoscaling=1
960        ModifyGizmo currentGroupObject=""
961        ModifyGizmo showInfo
962        ModifyGizmo infoWindow={651,303,1468,602}
963        ModifyGizmo endRecMacro
964        ModifyGizmo SETQUATERNION={0.573113,-0.115160,-0.275160,0.763255}
965EndMacro
966
967///////////////////////////
968//
969// unused - a simple line graph for each tube is much simpler
970//
971Window Gizmo_DetPanel() : GizmoPlot
972        PauseUpdate; Silent 1           // building window...
973        // Building Gizmo 7 window...
974        NewGizmo/W=(96,290,611,750)
975        ModifyGizmo startRecMacro=700
976        ModifyGizmo scalingOption=63
977        AppendToGizmo Surface=root:slices_L,name=surface0
978        ModifyGizmo ModifyObject=surface0,objectType=surface,property={ srcMode,0}
979        ModifyGizmo ModifyObject=surface0,objectType=surface,property={ surfaceCTab,ColdWarm}
980        AppendToGizmo Axes=boxAxes,name=axes0
981        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={0,axisRange,-1,-1,-1,1,-1,-1}
982        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={1,axisRange,-1,-1,-1,-1,1,-1}
983        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={2,axisRange,-1,-1,-1,-1,-1,1}
984        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={3,axisRange,-1,1,-1,-1,1,1}
985        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={4,axisRange,1,1,-1,1,1,1}
986        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={5,axisRange,1,-1,-1,1,-1,1}
987        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={6,axisRange,-1,-1,1,-1,1,1}
988        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={7,axisRange,1,-1,1,1,1,1}
989        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={8,axisRange,1,-1,-1,1,1,-1}
990        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={9,axisRange,-1,1,-1,1,1,-1}
991        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={10,axisRange,-1,1,1,1,1,1}
992        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={11,axisRange,-1,-1,1,1,-1,1}
993        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={-1,axisScalingMode,1}
994        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={-1,axisColor,0,0,0,1}
995        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={0,ticks,3}
996        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={1,ticks,3}
997        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={2,ticks,3}
998        ModifyGizmo modifyObject=axes0,objectType=Axes,property={-1,Clipped,0}
999        AppendToGizmo Surface=root:position_refined,name=surface1
1000        ModifyGizmo ModifyObject=surface1,objectType=surface,property={ fillMode,4}
1001        ModifyGizmo ModifyObject=surface1,objectType=surface,property={ srcMode,0}
1002        ModifyGizmo ModifyObject=surface1,objectType=surface,property={ surfaceCTab,Rainbow}
1003        ModifyGizmo setDisplayList=0, object=axes0
1004        ModifyGizmo setDisplayList=1, object=surface0
1005        ModifyGizmo autoscaling=1
1006        ModifyGizmo currentGroupObject=""
1007        ModifyGizmo showInfo
1008        ModifyGizmo infoWindow={550,23,1367,322}
1009        ModifyGizmo endRecMacro
1010        ModifyGizmo SETQUATERNION={0.499484,-0.278571,-0.448869,0.686609}
1011EndMacro
1012
1013
1014////////////////////////////////////
1015//
1016// An easy way to see the fit results to check if the peak locations all make sense.
1017//
1018Proc V_OpenPeakResultsGraph()
1019
1020        DoWindow/F V_PeakResultsGraph
1021        if(V_flag == 0)
1022                Make/O/D/N=5 tmpPeak,dummyLevel
1023                Make/O/D/N=128 tmpTube
1024               
1025                tmpPeak = position_refined[p][0]
1026                dummyLevel = WaveMax(tube0)
1027                tmpTube = tube0
1028               
1029                V_PeakResultsGraph()
1030        endif
1031
1032End
1033
1034
1035///////////////
1036//
1037// An easy way to see the fit results to check if the peak locations all make sense.
1038//
1039Window V_PeakResultsGraph() : Graph
1040        PauseUpdate; Silent 1           // building window...
1041        Display /W=(750,45,1161,376)/K=1 tmpTube vs tube_pixel
1042       
1043        ControlBar 50
1044       
1045       
1046        AppendToGraph dummyLevel vs tmpPeak
1047        ModifyGraph mode(dummyLevel)=3
1048        ModifyGraph marker(dummyLevel)=19
1049        ModifyGraph rgb(dummyLevel)=(1,16019,65535)
1050       
1051        SetVariable setvar0,pos={10.00,10.00},size={120.00,14.00},proc=V_TubePeakSetVarProc,title="Tube"
1052        SetVariable setvar0,limits={0,47,1},value= _NUM:0
1053       
1054        Label left "Counts"
1055        Label bottom "Pixel Number"
1056EndMacro
1057
1058//
1059// cycle through the tubes
1060//
1061Function V_TubePeakSetVarProc(sva) : SetVariableControl
1062        STRUCT WMSetVariableAction &sva
1063
1064        switch( sva.eventCode )
1065                case 1: // mouse up
1066                case 2: // Enter key
1067                case 3: // Live update
1068                        Variable dval = sva.dval
1069                        String sval = sva.sval
1070                       
1071                        Wave tmpPeak = tmpPeak
1072                        Wave dummyLevel = dummyLevel
1073                        Wave tmpTube = tmpTube
1074                       
1075                        Wave pos_ref = position_refined
1076                        Wave tube = $("tube"+num2str(dval))
1077                       
1078                        tmpPeak = pos_ref[p][dval]
1079                        dummyLevel = WaveMax(tube)
1080                        tmpTube = tube
1081               
1082                        break
1083                case -1: // control being killed
1084                        break
1085        endswitch
1086
1087        return 0
1088End
1089
1090////////////////////////////////////
1091//
1092// TODO
1093// -- document the "simple" save of the detector panels for import and subsequent fitting.
1094//
1095// takes the data from RAW, by default. This is OK, since even though whatever is in the calibration data
1096// of the file is used when loading into RAW, it is only used for the calculation of q. The actual data
1097// array is unchanged. Alternatively, the data could be pulled from the RawVSANS folder after a
1098// file catalog operation
1099//
1100Macro V_CopyDetectorsToRoot()
1101
1102        Duplicate/O root:Packages:NIST:VSANS:RAW:entry:instrument:detector_B:data data_B
1103
1104        Duplicate/O root:Packages:NIST:VSANS:RAW:entry:instrument:detector_ML:data data_ML
1105        Duplicate/O root:Packages:NIST:VSANS:RAW:entry:instrument:detector_MR:data data_MR
1106        Duplicate/O root:Packages:NIST:VSANS:RAW:entry:instrument:detector_MT:data data_MT
1107        Duplicate/O root:Packages:NIST:VSANS:RAW:entry:instrument:detector_MB:data data_MB
1108
1109        Duplicate/O root:Packages:NIST:VSANS:RAW:entry:instrument:detector_FL:data data_FL
1110        Duplicate/O root:Packages:NIST:VSANS:RAW:entry:instrument:detector_FR:data data_FR
1111        Duplicate/O root:Packages:NIST:VSANS:RAW:entry:instrument:detector_FT:data data_FT
1112        Duplicate/O root:Packages:NIST:VSANS:RAW:entry:instrument:detector_FB:data data_FB
1113       
1114End
1115
1116//
1117//
1118Macro V_SaveDetectorsITX()
1119// binary save makes each wave an individual file. Igor text groups them all into one file.
1120//      Save/C data_B,data_FB,data_FL,data_FR,data_FT,data_MB,data_ML,data_MR,data_MT
1121        Save/T/M="\r\n" data_B,data_FB,data_FL,data_FR,data_FT,data_MB,data_ML,data_MR,data_MT as "data_B++.itx"
1122
1123End
Note: See TracBrowser for help on using the repository browser.