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

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

a lot of little changes:

changed the name of the Raw Data display procedure file (removed test)

lots of bug fixes, moving items from the macros menu to proper locations, getting the file status to display properly, some error checking, and cleaning up a few TODO items.

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