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

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

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

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