source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_Marquee_Operations.ipf @ 1243

Last change on this file since 1243 was 1243, checked in by srkline, 2 years ago

added some test procedures to make it easier to calculate and patch the corrected beam center to a series of files. This new method can be found under the VSANS menu.

Changed the mask drawing objects to be semi-transparent to allow easier masking of bad spots - since you can see what you're doing now.

File size: 15.1 KB
Line 
1#pragma TextEncoding = "MacRoman"
2#pragma rtGlobals=3             // Use modern global access method and strict wave access.
3#pragma IgorVersion = 7.00
4
5
6// TODO:
7// x- get BoxSum Operational, for rudimentary Transmission calculations and ABS scale
8// x- need to be able to identify the detector panel from the Marquee selection
9// x- then do the count in (unscaled) coordinates of the actual data array
10//
11//
12// -- still many operations in (SANS)Marquee.ipf that will be useful to implement:
13//
14// -- writing the box coordinates and the counts (and error) to the data file
15// x- determining the beam center (centroid) of the selection
16//  -- writing beam center (centroid) to the file?
17//  -- a box sum over a range of files (with a plot)
18// -- box sum over annular regions
19// -- box sum over arcs
20// --  saving Graphics image
21// -- histogram of the counts
22//
23//
24//
25
26
27
28
29//
30//function will print marquee coordinates in axis terms, not detector terms
31//
32// will also calculate the sum in the box, and the error, and print it out
33//
34Function V_PrintMarqueeCoords() :  GraphMarquee
35        GetMarquee left,bottom
36        if(V_flag == 0)
37                Print "There is no Marquee"
38        else
39                Variable count,x1,x2,y1,y2,ct_err
40                x1 = V_left
41                x2 = V_right
42                y1 = V_bottom
43                y2 = V_top
44                printf "marquee left in bottom axis terms: %g\r",round(V_left)
45                printf "marquee right in bottom axis terms: %g\r",round(V_right)
46                printf "marquee bottom in left axis terms: %g\r",round(V_bottom)
47                printf "marquee top in left axis terms: %g\r",round(V_top)
48//              printf "**note that you must add 1 to each axis coordinate to get detector coordinates\r"
49               
50                // NOTE:
51                // this function MODIFIES x and y values on return, converting them to panel coordinates
52                // detector panel is identified from the (left,top) coordinate (x1,y2)
53                String detStr = V_FindDetStrFromLoc(x1,x2,y1,y2)               
54//              Printf "Detector = %s\r",detStr
55
56//
57                SVAR gCurDispType = root:Packages:NIST:VSANS:Globals:gCurDispType       
58                // this function will modify the x and y values (passed by reference) as needed to keep on the panel
59                V_KeepSelectionInBounds(x1,x2,y1,y2,detStr,gCurDispType)
60                Printf "%d;%d;%d;%d;\r",x1,x2,y1,y2
61
62
63                count = V_SumCountsInBox(x1,x2,y1,y2,ct_err,gCurDispType,detStr)
64
65               
66                Print "total counts = ",count
67                Print "err/counts = ",ct_err/count
68                print "average counts per pixel = ",count/(x2-x1)/(y2-y1)
69        endif
70End
71
72
73// NOTE:
74// this function MODIFIES x and y values on return, converting them to panel coordinates
75// detector panel is identified from the (left,top) coordinate (x1,y2)
76Function/S V_FindDetStrFromLoc(x1,x2,y1,y2)
77        Variable &x1,&x2,&y1,&y2
78       
79        // which images are here?
80        String detStr="",imStr,carriageStr
81        String currentImageRef,activeSubwindow,imageList
82        Variable ii,nIm,testX,testY,tab
83       
84        // which tab is selected? -this is the main graph panel (subwindow may not be the active one!)
85        ControlInfo/W=VSANS_Data tab0
86        tab = V_Value
87        if(tab == 0)
88                activeSubwindow = "VSANS_Data#det_panelsF"
89        elseif (tab == 1)
90                activeSubwindow = "VSANS_Data#det_panelsM"
91        else
92                activeSubwindow = "VSANS_Data#det_panelsB"
93        endif
94       
95        imageList = ImageNameList(activeSubwindow,";") 
96
97        nIm = ItemsInList(imageList,";")
98        if(nIm==0)
99                return("")              //problem, get out
100        endif
101
102        // images were added in the order TBLR, so look back through in the order RLBT, checking each to see if
103        // the xy value is found on that (scaled) array
104                               
105        // loop backwards through the list of panels (may only be one if on the back)
106        for(ii=nIm-1;ii>=0;ii-=1)
107                Wave w = ImageNameToWaveRef(activeSubwindow,StringFromList(ii, imageList,";"))
108               
109                // which, if any image is the mouse xy location on?
110                // use a multidemensional equivalent to x2pnt: (ScaledDimPos - DimOffset(waveName, dim))/DimDelta(waveName,dim)
111
112               
113                testX = ScaleToIndex(w,x1,0)
114                testY = ScaleToIndex(w,y2,1)
115               
116                if( (testX >= 0 && testX < DimSize(w,0)) && (testY >= 0 && testY < DimSize(w,1)) )
117                        // we're in-bounds on this wave
118                       
119                        // deduce the detector panel
120                        currentImageRef = StringFromList(ii, imageList,";")     //the image instance ##
121                        // string is "data", or "data#2" etc. - so this returns "", "1", "2", or "3"
122                        imStr = StringFromList(1, currentImageRef,"#")         
123                        carriageStr = activeSubWindow[strlen(activeSubWindow)-1]
124                       
125                        if(cmpstr(carriageStr,"B")==0)
126                                detStr = carriageStr
127                        else
128                                if(strlen(imStr)==0)
129                                        imStr = "9"                     // a dummy value so I can replace it later
130                                endif
131                                detStr = carriageStr+imStr              // "F2" or something similar
132                                detStr = ReplaceString("9", detStr, "T")        // ASSUMPTION :::: instances 0123 correspond to TBLR
133                                detStr = ReplaceString("1", detStr, "B")        // ASSUMPTION :::: this is the order that the panels
134                                detStr = ReplaceString("2", detStr, "L")        // ASSUMPTION :::: are ALWAYS added to the graph
135                                detStr = ReplaceString("3", detStr, "R")        // ASSUMPTION ::::
136                        endif
137                       
138                        Printf "Detector panel %s=(%d,%d)\r",detStr,testX,testY
139                       
140                        x1 = ScaleToIndex(w,x1,0)               // get all four marquee values to pass back to the calling function
141                        x2 = ScaleToIndex(w,x2,0)               // converted into detector coordinates
142                        y1 = ScaleToIndex(w,y1,1)
143                        y2 = ScaleToIndex(w,y2,1)
144                       
145                       
146                        ii = -1         //look no further, set ii to bad value to exit the for loop
147
148                endif
149               
150        endfor
151       
152        return(detStr)
153       
154End
155
156//testing function only, not called anywhere
157Function V_testKeepInBounds(x1,x2,y1,y2,detStr,folderStr)
158        Variable x1,x2,y1,y2
159        String detStr,folderStr
160       
161        V_KeepSelectionInBounds(x1,x2,y1,y2,detStr,folderStr)
162        Print x1,x2,y1,y2
163        return(0)
164End
165
166// for the given detector
167Function V_KeepSelectionInBounds(x1,x2,y1,y2,detStr,folderStr)
168        Variable &x1,&x2,&y1,&y2
169        String detStr,folderStr
170       
171        Variable pixelsX = V_getDet_pixel_num_x(folderStr,detStr)
172        Variable pixelsY = V_getDet_pixel_num_y(folderStr,detStr)
173
174        //keep selection in-bounds
175        x1 = (round(x1) >= 0) ? round(x1) : 0
176        x2 = (round(x2) <= (pixelsX-1)) ? round(x2) : (pixelsX-1)
177        y1 = (round(y1) >= 0) ? round(y1) : 0
178        y2 = (round(y2) <= (pixelsY-1)) ? round(y2) : (pixelsY-1)
179       
180        return(0)
181End
182
183
184
185//sums the data counts in the box specified by (x1,y1) to (x2,y2)
186//assuming that x1<x2, and y1<y2
187//the x,y values must also be in array coordinates[0] NOT scaled detector coords.
188//
189// accepts arbitrary detector coordinates. calling function is responsible for
190// keeping selection in bounds
191//
192Function V_SumCountsInBox(x1,x2,y1,y2,ct_err,type,detStr)
193        Variable x1,x2,y1,y2,&ct_err
194        String type,detStr
195       
196        Variable counts = 0,ii,jj,err2_sum
197       
198// get the waves of the data and the data_err
199        Wave w = V_getDetectorDataW(type,detStr)
200        Wave data_err = V_getDetectorDataErrW(type,detStr)
201
202       
203        err2_sum = 0            // running total of the squared error
204        ii=x1
205        jj=y1
206        do
207                do
208                        counts += w[ii][jj]
209                        err2_sum += data_err[ii][jj]*data_err[ii][jj]
210                        jj+=1
211                while(jj<=y2)
212                jj=y1
213                ii+=1
214        while(ii<=x2)
215       
216        err2_sum = sqrt(err2_sum)
217        ct_err = err2_sum
218       
219//      Print "error = ",ct_err
220//      Print "error/counts = ",ct_err/counts
221       
222        Return (counts)
223End
224
225
226//sums the data counts in the box specified by (x1,y1) to (x2,y2)
227//assuming that x1<x2, and y1<y2
228//the x,y values must also be in array coordinates[0] NOT scaled detector coords.
229//
230// accepts arbitrary detector coordinates. calling function is responsible for
231// keeping selection in bounds
232//
233// basically the same as V_SumCountsInBox, except the PBR value has been removed so that the
234// function can be used from the command line
235//
236Function V_SumCountsInBox_Cmd(x1,x2,y1,y2,type,detStr)
237        Variable x1,x2,y1,y2
238        String type,detStr
239       
240        Variable counts = 0,ii,jj,err2_sum,ct_err
241       
242// get the waves of the data and the data_err
243        Wave w = V_getDetectorDataW(type,detStr)
244        Wave data_err = V_getDetectorDataErrW(type,detStr)
245
246                       
247        err2_sum = 0            // running total of the squared error
248        ii=x1
249        jj=y1
250        do
251                do
252                        counts += w[ii][jj]
253                        err2_sum += data_err[ii][jj]*data_err[ii][jj]
254                        jj+=1
255                while(jj<=y2)
256                jj=y1
257                ii+=1
258        while(ii<=x2)
259       
260        err2_sum = sqrt(err2_sum)
261        ct_err = err2_sum
262
263        Print "sum of counts = ",counts
264        Print "error = ",ct_err
265        Print "error/counts = ",ct_err/counts
266       
267        Return (counts)
268End
269
270Proc pV_SumCountsInBox_Cmd(x1,x2,y1,y2,type,detStr)
271        Variable x1=280,x2=430,y1=350,y2=1020
272        String type="RAW",detStr="B"
273       
274        V_SumCountsInBox_Cmd(x1,x2,y1,y2,type,detStr)
275End
276
277
278
279Function V_Find_BeamCentroid() :  GraphMarquee
280
281//      //get the current displayed data (so the correct folder is used)
282//      SVAR cur_folder=root:myGlobals:gDataDisplayType
283//      String dest = "root:Packages:NIST:" + cur_folder
284       
285        Variable xzsum,yzsum,zsum,xctr,yctr
286        Variable left,right,bottom,top,ii,jj,counts
287        Variable x_mm_sum,y_mm_sum,x_mm,y_mm
288        Variable xRef,yRef
289
290       
291        GetMarquee left,bottom
292        if(V_flag == 0)
293                Print "There is no Marquee"
294        else
295                left = round(V_left)            //get integer values for selection limits
296                right = round(V_right)
297                top = round(V_top)
298                bottom = round(V_bottom)
299
300                // detector panel is identified from the (left,top) coordinate (x1,y2)
301                String detStr = V_FindDetStrFromLoc(left,right,bottom,top)             
302        //      Printf "Detector = %s\r",detStr
303       
304        //
305                SVAR gCurDispType = root:Packages:NIST:VSANS:Globals:gCurDispType       
306                // this function will modify the x and y values (passed by reference) as needed to keep on the panel
307                V_KeepSelectionInBounds(left,right,bottom,top,detStr,gCurDispType)
308                Print left,right,bottom,top
309                       
310               
311                // selection valid now, calculate beamcenter
312                // get the waves of the data and the data_err
313                Wave data = V_getDetectorDataW(gCurDispType,detStr)
314                Wave data_err = V_getDetectorDataErrW(gCurDispType,detStr)
315               
316                // get the real-space information
317                String destPath = "root:Packages:NIST:VSANS:"+gCurDispType
318                Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
319                Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
320       
321                xzsum = 0
322                yzsum = 0
323                zsum = 0
324                x_mm_sum = 0
325                y_mm_sum = 0
326               
327                // count over rectangular selection, doing each row, L-R, bottom to top
328                ii = bottom -1
329                do
330                        ii +=1
331                        jj = left-1
332                        do
333                                jj += 1
334                                counts = data[jj][ii]
335                                xzsum += jj*counts
336                                yzsum += ii*counts
337                                zsum += counts
338                               
339                                x_mm_sum += data_realDistX[jj][ii]*counts
340                                y_mm_sum += data_realDistY[jj][ii]*counts
341                        while(jj<right)
342                while(ii<top)
343               
344                xctr = xzsum/zsum
345                yctr = yzsum/zsum
346               
347                x_mm = x_mm_sum/zsum
348                y_mm = y_mm_sum/zsum
349                // add 1 to each to get to detector coordinates (1,128)
350                // rather than the data array which is [0,127]
351//              xctr+=1
352//              yctr+=1
353
354// correct for the zero position (y-position) on the L/R panels not being exactly equal
355// the lateral scan data from Dec 2018 is used to correct this. The span of zero points
356// is relatively small (+- 0.5 pixel) but is significant for data using graphite monochromator
357//     
358        // check that the correction waves exist, if not, generate them V_TubeZeroPointTables()
359                Wave/Z tube_num = $("root:Packages:NIST:VSANS:Globals:tube_" + detStr)
360                Wave/Z yCtr_tube = $("root:Packages:NIST:VSANS:Globals:yCtr_" + detStr)
361                if(!WaveExists(tube_num))
362                        Execute "V_TubeZeroPointTables()"
363                        Wave/Z tube_num = $("root:Packages:NIST:VSANS:Globals:tube_" + detStr)
364                        Wave/Z yCtr_tube = $("root:Packages:NIST:VSANS:Globals:yCtr_" + detStr)
365                endif
366               
367                Variable yCorrection = interp(xCtr,tube_num,yCtr_tube)
368                Variable yPixSize = V_getDet_y_pixel_size(gCurDispType,detStr)
369                yPixSize /= 10          // convert mm to cm
370                // offsets were determined in Dec 2018 using:
371                // FR tube # 7 = 61.70 pix
372                // MR tube # 10 = 61.94 pix
373               
374                Print "X-center (in array coordinates 0->n-1 ) = ",xctr
375                Print "Y-center (in array coordinates 0->n-1 ) = ",yctr
376               
377                Print "X-center (cm) = ",x_mm/10
378                Print "Y-center (cm) = ",y_mm/10
379
380                if(cmpstr(detStr,"FR") == 0)
381                        Print "Reference Y-Center is corrected for tube #7 zero position"               
382
383                        yCorrection = 61.70 - yCorrection
384                        Print "yCorrection (pix) = ",yCorrection
385                        Print "yCorrection (cm) = ",yCorrection*yPixSize
386                        Print "FRONT Reference X-center (cm) = ",x_mm/10
387                        Print "FRONT Reference Y-center (cm) = ",y_mm/10 + yCorrection*yPixSize
388                        xRef = x_mm/10
389                        yRef = y_mm/10 + yCorrection*yPixSize
390                endif
391
392                if(cmpstr(detStr,"MR") == 0)
393                        Print "Reference Y-Center is corrected for tube #10 zero position"             
394
395                        yCorrection = 61.94 - yCorrection
396                        Print "yCorrection (pix) = ",yCorrection
397                        Print "yCorrection (cm) = ",yCorrection*yPixSize
398                        Print "MIDDLE Reference X-center (cm) = ",x_mm/10
399                        Print "MIDDLE Reference Y-center (cm) = ",y_mm/10 + yCorrection*yPixSize
400                        xRef = x_mm/10
401                        yRef = y_mm/10 + yCorrection*yPixSize
402                endif
403               
404// if measured on the LEFT panel, convert to the RIGHT coordinates for the reference value     
405// these corrections are exactly the opposite (subtract, not add) of what is done in V_fDeriveBeamCenters(xFR,yFR,xMR,yMR)
406                if(cmpstr(detStr,"FL") == 0)
407                        Print "FRONT Reference X-center (cm) = ",x_mm/10 - kBCtrOffset_FL_x     // NEW Dec 2018 values
408                        Print "FRONT Reference Y-center (cm) = ",y_mm/10 - kBCtrOffset_FL_y
409                        xRef = x_mm/10 - kBCtrOffset_FL_x
410                        yRef = y_mm/10 - kBCtrOffset_FL_y
411                endif
412               
413                if(cmpstr(detStr,"ML") == 0)
414                        Print "MIDDLE Reference X-center (cm) = ",x_mm/10 - kBCtrOffset_ML_x
415                        Print "MIDDLE Reference Y-center (cm) = ",y_mm/10 - kBCtrOffset_ML_y
416                        xRef = x_mm/10 - kBCtrOffset_ML_x
417                        yRef = y_mm/10 - kBCtrOffset_ML_y
418                endif
419        endif
420
421// TODO
422// ?? store the xy reference values somewhere so that the conversion to proper
423// beam center values can be done automatically, rather than copying numbers into a procedure
424//
425// - either I need 6 globals for the three panels, or I need to store the values in the
426// reduction block of the file (comment?) - but I don't have the fileName here - could I find it
427// somewhere? gFileList in the current data folder?
428//
429        String ctrStr=""
430        if(cmpstr(detStr,"B") == 0)
431                xRef = xCtr
432                yRef = yCtr                     //these are in pixels
433        endif
434        sprintf ctrStr,"XREF=%g;YREF=%g;",xRef,yRef
435        SVAR gFileList = $("root:Packages:NIST:VSANS:"+gCurDispType+":gFileList")
436       
437        V_writeReductionComments(gFileList,ctrStr)
438       
439       
440       
441        //back to root folder (redundant)
442        SetDataFolder root:
443       
444End
445
446
447
448//
449//function writes new box coordinates to the data file
450//
451// also writes the panel where the coordinates were set (non-nice field in /reduction)
452//
453Function V_UpdateBoxCoords() :  GraphMarquee
454        GetMarquee left,bottom
455        if(V_flag == 0)
456                Print "There is no Marquee"
457        else
458                Variable count,x1,x2,y1,y2,ct_err
459                x1 = V_left
460                x2 = V_right
461                y1 = V_bottom
462                y2 = V_top
463
464//              Print x1,x2,y1,y2
465
466                // NOTE:
467                // this function MODIFIES x and y values on return, converting them to panel coordinates
468                // detector panel is identified from the (left,top) coordinate (x1,y2)
469                String detStr = V_FindDetStrFromLoc(x1,x2,y1,y2)               
470//
471                SVAR gCurDispType = root:Packages:NIST:VSANS:Globals:gCurDispType
472                string boxStr
473                // this function will modify the x and y values (passed by reference) as needed to keep on the panel
474                V_KeepSelectionInBounds(x1,x2,y1,y2,detStr,gCurDispType)
475                sprintf boxStr,"%d;%d;%d;%d;",x1,x2,y1,y2
476
477//              Print x1,x2,y1,y2
478               
479                SVAR gCurrentFile = root:Packages:NIST:VSANS:Globals:gLastLoadedFile            //for the status of the display
480
481                V_writeBoxCoordinates(gCurrentFile,V_List2NumWave(boxStr,";","inW"))
482       
483                V_writeReduction_BoxPanel(gCurrentFile,detStr)
484
485//              count = V_SumCountsInBox(x1,x2,y1,y2,ct_err,gCurDispType,detStr)               
486//              Print "counts = ",count
487//              Print "err/counts = ",ct_err/count
488
489                // kill the file from RawVSANS so that the updated box coordinates will be re-read in
490                //
491                V_KillNamedDataFolder(gCurrentFile)
492        endif
493End
Note: See TracBrowser for help on using the repository browser.