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

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

updates to the beam center corrections that take into account the lateral variation of zero points of the tubes. This is an important correction to the reference beam center, especially when the graphite monochrmoator is used.

various bug fixes included as part of reduction , largely to catch missing or incorrect information in the file header.

New "patch" functions have been added for number of guides, beam stops, source aperture, etc. which were missing in the file and caused the resolution calculation to be totally incorrect.

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