source: sans/SANSReduction/branches/kline_29MAR07/Put in User Procedures/SANS_Reduction_v5.00/Marquee.ipf @ 66

Last change on this file since 66 was 40, checked in by srkline, 16 years ago

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