source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/Marquee.ipf @ 418

Last change on this file since 418 was 418, checked in by ajj, 14 years ago

Moving data folders into root:Packages:NIST

This could be hairy.

File size: 15.4 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma version=5.0
3#pragma IgorVersion=6.0
4
5////////////
6// vers 1.21 3 may 06 total transmission incorporated (BSG)
7//
8//**************************
9// Vers 1.2 091901
10//
11// marquee functions are used to:
12//
13// locate the beamcenter
14// set the (x,y) box coordinates for the sum for transmission calculation
15// read out coordinates
16// do a "box sum" of the same box over a range of files
17// do a 2D Gaussian fit over a selected range
18// do a save of the current image (with colorBar) - as a graphics image
19//
20//***************************
21
22//sums the data counts in the box specified by (x1,y1) to (x2,y2)
23//assuming that x1<x2, and y1<y2
24//the x,y values must also be in axis coordinates[0,127] NOT (1,128) detector coords.
25//
26// accepts arbitrary detector coordinates. calling function is responsible for
27// keeping selection in bounds
28Function SumCountsInBox(x1,x2,y1,y2,type)
29        Variable x1,x2,y1,y2
30        String type
31       
32        Variable counts = 0,ii,jj
33       
34        String dest =  "root:Packages:NIST:"+type
35       
36        //check for logscale data, but don't change the data
37        NVAR gIsLogScale = $(dest + ":gIsLogScale")
38        if (gIsLogScale)
39                wave w=$(dest + ":linear_data")
40        else
41                wave w=$(dest + ":data")
42        endif
43       
44        ii=x1
45        jj=y1
46        do
47                do
48                        counts += w[ii][jj]
49                        jj+=1
50                while(jj<=y2)
51                jj=y1
52                ii+=1
53        while(ii<=x2)
54       
55        Return (counts)
56End
57
58
59//from a marquee selection:
60//calculates the sum of counts in the box, and records this count value in
61//globals for both patch and Trans routines, then records the box coordinates
62//and the count value for the box in the header of the (currently displayed)
63//empty beam file for use later when calculating transmissions of samples
64//these values are written to unused (analysis) fields in the data header
65//4 integers and one real value are written
66//
67//re-written to work with Transmission panel as well as PatchPanel
68//which now bothe work from the same folder (:Patch)
69//
70Function SetXYBoxCoords() :  GraphMarquee
71
72        GetMarquee left,bottom
73        if(V_flag == 0)
74                Abort "There is no Marquee"
75        Endif
76        SVAR dispType=root:myGlobals:gDataDisplayType
77        if(cmpstr(dispType,"SAM")!=0)
78                DoAlert 0, "You can only use SetXYBox on SAM data files"
79                return(1)
80        endif
81        //printf "marquee left in bottom axis terms: %g\r",round(V_left)
82        //printf "marquee right in bottom axis terms: %g\r",round(V_right)
83        //printf "marquee top in left axis terms: %g\r",round(V_top)
84        //printf "marquee bottom in left axis terms: %g\r",round(V_bottom)
85        Variable x1,x2,y1,y2
86        x1 = round(V_left)
87        x2 = round(V_right)
88        y1 = round(V_bottom)
89        y2 = round(V_top)
90       
91        KeepSelectionInBounds(x1,x2,y1,y2)
92       
93        //check to make sure that Patch and Trans data folders exist for writing of global variables
94        If( ! (DataFolderExists("root:myGlobals:Patch"))  )
95                Execute "InitializePatchPanel()"
96        Endif
97        //check to make sure that Patch and Trans data folders exist for writing of global variables
98        If( ! (DataFolderExists("root:myGlobals:TransHeaderInfo"))  )
99                Execute "InitializeTransPanel()"
100        Endif
101       
102        //write string as keyword-packed string, to use IGOR parsing functions
103        String msgStr = "X1="+num2str(x1)+";"
104        msgStr += "X2="+num2str(x2)+";"
105        msgStr += "Y1="+num2str(y1)+";"
106        msgStr += "Y2="+num2str(y2)+";"
107        String/G root:myGlobals:Patch:gPS3 = msgStr
108        String/G root:myGlobals:Patch:gEmpBox = msgStr
109        //changing this global wil update the display variable on the TransPanel
110        String/G root:myGlobals:TransHeaderInfo:gBox = msgStr
111       
112        //sum the counts in the patch - working on the SAM data, to be sure that it's normalized
113        //to the same monitor counts and corrected for detector deadtime
114        String type = "SAM"
115        Variable counts
116        counts = SumCountsInBox(x1,x2,y1,y2,type)
117        Print " marquee counts =",counts
118        //Set the global gTransCts
119        Variable/G root:myGlobals:Patch:gTransCts = counts
120       
121        //now change the extra variables in the empty beam file
122        //get the filename from the SAM folder (there will only be one file)
123        SVAR partialName = root:Packages:NIST:SAM:FileList
124        //construct valid filename, then prepend path
125        String tempName = FindValidFilename(partialName)
126        Print "in marquee",partialName
127        //Print tempName
128        if(cmpstr(tempName,"")==0)
129                //file not found, get out
130                Abort "file not found, marquee"
131        Endif
132        //name is ok, prepend path to tempName for read routine
133        PathInfo catPathName
134        if (V_flag == 0)
135                //path does not exist - no folder selected
136                Abort "no path selected"
137        else
138                String filename = S_path + tempName
139        endif
140       
141        if(cmpstr(filename,"no file selected")==0)
142                Abort "no file selected"
143        Endif
144       
145        WriteXYBoxToHeader(filename,x1,x2,y1,y2)
146       
147        Print counts, " counts in XY box"
148        WriteBoxCountsToHeader(filename,counts)
149       
150End
151
152//finds the beam center (the centroid) of the selected region
153//and simply prints out the results to the history window
154//values are printed out in detector coordinates, not IGOR coords.
155//
156Function FindBeamCenter() :  GraphMarquee
157
158        //get the current displayed data (so the correct folder is used)
159        SVAR cur_folder=root:myGlobals:gDataDisplayType
160        String dest = "root:Packages:NIST:" + cur_folder
161       
162        Variable xzsum,yzsum,zsum,xctr,yctr
163        Variable left,right,bottom,top,ii,jj,counts
164       
165        // data wave is hard-wired in as the displayed data
166        NVAR dataIsLog=$(dest + ":gIsLogScale")         //check for log-scaling in current data folder
167        if (dataIsLog)
168                wave data=$(dest + ":linear_data")
169        else
170                wave data=$(dest + ":data")
171        endif
172       
173        GetMarquee left,bottom
174        if(V_flag == 0)
175                Print "There is no Marquee"
176        else
177                left = round(V_left)            //get integer values for selection limits
178                right = round(V_right)
179                top = round(V_top)
180                bottom = round(V_bottom)
181               
182                KeepSelectionInBounds(left,right,bottom,top)
183               
184                // selection valid now, calculate beamcenter
185                xzsum = 0
186                yzsum = 0
187                zsum = 0
188                // count over rectangular selection, doing each row, L-R, bottom to top
189                ii = bottom -1
190                do
191                        ii +=1
192                        jj = left-1
193                        do
194                                jj += 1
195                                counts = data[jj][ii]
196                                xzsum += jj*counts
197                                yzsum += ii*counts
198                                zsum += counts
199                        while(jj<right)
200                while(ii<top)
201               
202                xctr = xzsum/zsum
203                yctr = yzsum/zsum
204               
205                // add 1 to each to get to detector coordinates (1,128)
206                // rather than the data array which is [0,127]
207                xctr+=1
208                yctr+=1
209               
210                Print "X-center (in detector coordinates) = ",xctr
211                Print "Y-center (in detector coordinates) = ",yctr
212        endif
213       
214        //back to root folder (redundant)
215        SetDataFolder root:
216       
217End
218
219//still need to error check - out-of-bounds...waves exist.
220// allows a 2D Gaussian fit to a selected region of data in a SANS_Data window
221//puts up a new graph with the fitted contour
222Function Do_2D_Gaussian_Fit() :  GraphMarquee
223        String topWin=WinName(0,1)              //top *graph* window
224        //exit nicely if not in the Data display window
225        if(cmpstr(topWin,"SANS_Data") != 0)
226                DoAlert 0,"2D Gaussian fitting is only available from the Data Display Window"
227                return(1)
228        Endif
229       
230        GetMarquee/K left,bottom
231        Variable x1,x2,y1,y2,qxlo,qxhi,qylo,qyhi
232        if(V_flag == 0)
233                Print "There is no Marquee"
234        else
235                String junk="",df=""
236               
237                //**hard-wired info about the x-y q-scales
238                qxlo = DimOffset(root:myGlobals:q_x_axis,0)
239                qxhi = DimDelta(root:myGlobals:q_x_axis,0) + qxlo
240//              Print "qxlo,qxhi = ",qxlo,qxhi
241                Wave w=$"root:myGlobals:q_y_axis"
242                qylo=w[0]
243                qyhi=w[1]
244//              print "qylo,qyhi = ",qylo,qyhi
245               
246                junk=ImageInfo("SANS_Data","data",0)
247                df=StringByKey("ZWAVEDF", junk,":",";")
248//              print df
249                Duplicate/O $(df+"data") data,data_err
250                data_err=sqrt(data)             //for weighting
251               
252// comment out the SetScale lines if you want the result in terms of pixels as a way of
253// measuring the beam center. Note that you need to ADD ONE to fitted x0 and y0 to get detector
254// coordinates rather than the zero-indexed array. 2D fitting does have the benefit of
255// reporting error bars on the xy (if you believe that 2D gaussian is correct)         
256                SetScale/I x qxlo,qxhi,"",data
257                SetScale/I y qylo,qyhi,"",data
258               
259                Display /W=(10,50,361,351) /K=1
260                AppendImage data
261                ModifyImage data ctab= {*,*,Grays,0}
262                ModifyGraph width={Plan,1,bottom,left}
263                ModifyGraph mirror=2
264                ModifyGraph lowTrip=1e-04
265                ModifyImage data cindex=$"root:myGlobals:NIHColors"
266                SVAR angst = root:myGlobals:gAngstStr
267                Label bottom "Qx ("+angst+"\\S-1\\M)"
268                Label left "Qy ("+angst+"\\S-1\\M)"
269
270                //keep selection in-bounds
271                x1=V_left
272                x2=V_right
273                y1=V_bottom
274                y2=V_top
275                KeepSelectionInBounds(x1,x2,y1,y2)
276
277                //cross correlation coefficent (K6) must be between 0 and 1, need constraints
278                Make/O/T/N=2 temp_constr
279                temp_constr = {"K6>0","K6<1"}
280               
281                CurveFit/N Gauss2D data[x1,x2][y1,y2] /I=1 /W=data_err /D /R /A=0 /C=temp_constr
282               
283                Killwaves/Z temp_constr
284        endif
285End
286
287// to save the image, simply invoke the IGOR menu item for saving graphics
288//
289Function SaveSANSGraphic() : GraphMarquee
290       
291        NVAR isDemoVersion=root:myGlobals:isDemoVersion
292        if(isDemoVersion==1)
293                //      comment out in DEMO_MODIFIED version, and show the alert
294                DoAlert 0,"This operation is not available in the Demo version of IGOR"
295        else
296                DoIGORMenu "File","Save Graphics"
297        endif
298End
299
300//does a sum over each of the files in the list over the specified range
301// x,y are assumed to already be in-bounds of the data array
302// output is dumped to the command window
303//
304Function DoBoxSum(fileStr,x1,x2,y1,y2)
305        String fileStr
306        Variable x1,x2,y1,y2
307       
308        //parse the list of file numbers
309        String fileList="",item="",pathStr="",fullPath=""
310        Variable ii,num,err,cts
311       
312        PathInfo catPathName
313        If(V_Flag==0)
314                Abort "no path selected"
315        Endif
316        pathStr = S_Path
317       
318        fileList=ParseRunNumberList(fileStr)
319        num=ItemsInList(fileList,",")
320       
321        //loop over the list
322        //add each file to SAM (to normalize to monitor counts)
323        //sum over the box
324        //print the results
325        Make/O/N=(num) FileID,BoxCounts
326        Print "Results are stored in root:FileID and root:BoxCounts waves"
327        for(ii=0;ii<num;ii+=1)
328                item=StringFromList(ii,fileList,",")
329                FileID[ii] = GetRunNumFromFile(item)            //do this here, since the list is now valid
330                fullPath = pathStr+item
331                ReadHeaderAndData(fullPath)
332//              String/G root:myGlobals:gDataDisplayType="RAW"
333//              fRawWindowHook()
334                err = Raw_to_work("SAM")
335                String/G root:myGlobals:gDataDisplayType="SAM" 
336                fRawWindowHook()
337                cts=SumCountsInBox(x1,x2,y1,y2,"SAM")
338                BoxCounts[ii]=cts
339                Print item+" counts = ",cts
340        endfor
341       
342        DoBoxGraph(FileID,BoxCounts)
343       
344        return(0)
345End
346
347Function DoBoxGraph(FileID,BoxCounts)
348        Wave FileID,BoxCounts
349       
350        Sort FileID BoxCounts,FileID            //sort the waves, in case the run numbers were entered out of numerical order
351       
352        Display /W=(5,44,383,306) BoxCounts vs FileID
353        ModifyGraph mode=4
354        ModifyGraph marker=8
355        ModifyGraph grid=2
356        ModifyGraph mirror=2
357        Label left "Counts (per 10^8 monitor counts)"
358        Label bottom "Run Number"
359       
360        return(0)
361End
362
363//
364// promts the user for a range of file numbers to perform the sum over
365// list must be comma delimited numbers (or dashes) just as in the BuildProtocol panel
366// the (x,y) range is already selected from the marquee
367//
368Function BoxSum() :  GraphMarquee
369        GetMarquee left,bottom
370        if(V_flag == 0)
371                Abort "There is no Marquee"
372        Endif
373        SVAR dispType=root:myGlobals:gDataDisplayType
374        if(cmpstr(dispType,"RealTime")==0)
375                Print "Can't do a BoxSum for a RealTime file"
376                return(1)
377        endif
378        Variable x1,x2,y1,y2
379        x1 = V_left
380        x2 = V_right
381        y1 = V_bottom
382        y2 = V_top
383        KeepSelectionInBounds(x1,x2,y1,y2)
384       
385        String fileStr="",msgStr="Enter a comma-delimited list of run numbers, use dashes for ranges"
386        Prompt fileStr,msgStr
387        DoPrompt "Pick the file range",fileStr
388        Print "fileStr = ",fileStr
389        printf "(x1,x2) (y1,y2) = (%d,%d) (%d,%d)\r",x1,x2,y1,y2
390       
391        DoBoxSum(fileStr,x1,x2,y1,y2)
392       
393        return(0)
394End     
395
396//function that keeps the marquee selection in the range [0,127] inclusive
397// (igor coordinate system)
398// uses pass-by reference!
399//
400// x1 = left
401// x2 = right
402// y1 = bottom
403// y2 = top
404//
405// accepts any detector size
406Function KeepSelectionInBounds(x1,x2,y1,y2)
407        Variable &x1,&x2,&y1,&y2
408       
409        NVAR pixelsX = root:myGlobals:gNPixelsX
410        NVAR pixelsY = root:myGlobals:gNPixelsY
411       
412        //keep selection in-bounds
413        x1 = (round(x1) >= 0) ? round(x1) : 0
414        x2 = (round(x2) <= (pixelsX-1)) ? round(x2) : (pixelsX-1)
415        y1 = (round(y1) >= 0) ? round(y1) : 0
416        y2 = (round(y2) <= (pixelsY-1)) ? round(y2) : (pixelsY-1)
417        return(0)
418End
419
420//testing function, not used
421Function testKeepInBounds(x1,x2,y1,y2)
422        Variable x1,x2,y1,y2
423       
424        KeepSelectionInBounds(x1,x2,y1,y2)
425        Print x1,x2,y1,y2
426        return(0)
427End
428
429// generates a histogram of the data as defined by the marquee. The longer dimension of the marquee
430// becomes the x-axis of the histogram (this may need to be changed for some odd case). Pixel range specified
431// by the marquee is inclusive, and is automatically kept in-bounds
432//
433// The counts over the (short) dimension are averaged, and plotted vs. the pixel position.
434// Pixel position is reported as Detector coordinates (1,128). Counts are whatever the current display
435// happens to be.
436//
437Function SANS_Histogram() :  GraphMarquee
438        GetMarquee left,bottom
439        if(V_flag == 0)
440                Abort "There is no Marquee"
441        endif
442        //
443        Variable count,x1,x2,y1,y2,xwidth,ywidth,vsX=1,xx,yy
444        x1 = V_left
445        x2 = V_right
446        y1 = V_bottom
447        y2 = V_top
448        KeepSelectionInBounds(x1,x2,y1,y2)
449        Print "x1,x2,y1,y2 (det) =",x1+1,x2+1,y1+1,y2+1
450        //determine whether to do x vs y or y vs x
451        xwidth=x2-x1
452        ywidth=y2-y1
453        if(xwidth < ywidth)
454                vsX=0           //sum and graph vs Y
455        endif
456        SVAR cur_folder=root:myGlobals:gDataDisplayType
457        WAVE data=$("root:Packages:NIST:"+cur_folder+":data")           //don't care if it's log or linear scale
458        Make/O/N=(max(xwidth,ywidth)+1) Position,AvgCounts
459        AvgCounts=0
460        //set position wave
461        if(vsX)
462                position=p+x1
463        else
464                position=p+y1
465        endif
466        //convert the position to Detector coordinates
467        position += 1
468       
469        //Compute the histogram (manually)
470        if(vsX)
471                for(xx=x1;xx<=x2;xx+=1)         //outer loop is the "x-axis"
472                        for(yy=y1;yy<=y2;yy+=1)
473                                AvgCounts[xx-x1] += data[xx][yy]
474                        endfor
475                endfor
476                AvgCounts /= (ywidth+1)
477        else
478                for(yy=y1;yy<=y2;yy+=1)
479                        for(xx=x1;xx<=x2;xx+=1)
480                                AvgCounts[yy-y1] += data[xx][yy]
481                        endfor
482                endfor
483                AvgCounts /= (xwidth+1)
484        endif
485        GetMarquee/K            //to keep from drawing the marquee on the new histo graph
486        //draw the graph, or just bring to the front with the new data
487        DoWindow/F SANS_Histo
488        if(V_Flag != 1)
489                Draw_Histo()
490        endif
491       
492        return(0)
493End
494
495//draws the histogram of the 2d data as specified by AvgCounts and Position
496//both wave are assumed to exist in the data folder. The SANS_Histogram() marquee
497//operation is responsible for creating them.
498//
499Function Draw_Histo()
500        Display /W=(197,329,567,461)/K=1 AvgCounts vs Position
501        DoWindow/C SANS_Histo
502        DoWindow/T SANS_Histo,"Histogram"
503        ModifyGraph mode=5,grid=1,mirror=2
504        ModifyGraph rgb=(21845,21845,21845)
505        ModifyGraph standoff=0
506        ModifyGraph hbFill=2
507        ModifyGraph useNegPat=1
508        ModifyGraph usePlusRGB=1
509        ModifyGraph useNegRGB=1
510        ModifyGraph hBarNegFill=2
511        ModifyGraph negRGB=(0,0,65535)
512        SetAxis/A/N=2 left
513        Label left "Counts"
514        Label bottom "Pixel (detector coordinates)"
515End
516
517//function will print marquee coordinates in axis terms, not detector terms
518//since IGOR is [0][127] and detector is (1,128)
519Function PrintMarqueeCoords() :  GraphMarquee
520        GetMarquee left,bottom
521        if(V_flag == 0)
522                Print "There is no Marquee"
523        else
524                Variable count,x1,x2,y1,y2
525                x1 = V_left
526                x2 = V_right
527                y1 = V_bottom
528                y2 = V_top
529                printf "marquee left in bottom axis terms: %g\r",round(V_left)
530                printf "marquee right in bottom axis terms: %g\r",round(V_right)
531                printf "marquee bottom in left axis terms: %g\r",round(V_bottom)
532                printf "marquee top in left axis terms: %g\r",round(V_top)
533                printf "**note that you must add 1 to each axis coordinate to get detector coordinates\r"
534               
535                KeepSelectionInBounds(x1,x2,y1,y2)
536                SVAR cur_folder=root:myGlobals:gDataDisplayType
537                count = SumCountsInBox(x1,x2,y1,y2,cur_folder)
538                Print "counts = "+ num2str(count)
539        endif
540End
Note: See TracBrowser for help on using the repository browser.