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

Last change on this file since 794 was 794, checked in by srkline, 12 years ago

Lots of changes:
-2D resolution smearing
-2D error propagation

1) 2D resolution smearing has been corrected to use sigma (perp) correctly
rather than phi. This shows up in the quadrature loop in all of the 2D models
and in the Smear_2D "generic" function.

2) 1D resolutionn smearing is now corrected to account for the integration
range of +/- 3 sigma (99.73% of distribution). Now the result is divided by
0.9973 to rescale it to the correct value.

3) Smeared models are now AAO to improve speed and to allow easier use with
functions that are inherently AAO. No changes are needed, since the call is
behind the scenes, replacing Smear_Model_N() with Smear_Model_N_AAO().

4) in PlotUtils_2D added functions to re-bin the QxQy? data into a 1D format
BinQxQy_to_1D(). This also re-bins the errors in two ways, adding the per-pixel
errors in quadrature, or the deviation from the mean of the intensity. Manually
editing the intensity declaration allows 2D->1D binning of smeared models.

5) Per-pixel propagation of errors has been added through the entire data
reduction sequence. Data errors are generated on loading using Poisson
statistics (specifically tailored for accuracy at low counts) and then is
propagated through each manipulation of the data using standard error
propagation. The error matrix is linear_data_error. As a by-product, all of
the math operations on data are explicitly done on linear_data, to avoid
any potential mistakes of log/linear scaling. Results of this propagation
largely match J. Barker's /ERROR routines from the VAX, with some differences
at low data count values (as expected) and at higher count values near the
beam stop (somewhat unexpected). This per-pixel error is ouput in the QxQy_ASCII
data files. NO CHANGE has been made to the 1D data, which uses the deviation from
the mean as the error - since this is correct.

6) Added tables for the uncertainty in attenuator transmission (from JGB)

7) added two more REAL values to the VAX header to store information
necessary for error propagation. These are couting error that are part of
the transmission error and of the absolute scaling error. These add Read/Write?
functions in NCNR_DataReadWrite

The transmission error (one standard deviation) is written at byte 396 (4 bytes)

RealsRead?[41]

The Box Sum error (one standard deviation) is written at byte 400 (4 bytes)

RealsRead?[42]

File size: 22.1 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma version=5.0
3#pragma IgorVersion=6.1
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,ct_err,type)
29        Variable x1,x2,y1,y2,&ct_err
30        String type
31       
32        Variable counts = 0,ii,jj,err2_sum
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        wave data_err = $(dest + ":linear_data_error")
44       
45        err2_sum = 0            // running total of the squared error
46        ii=x1
47        jj=y1
48        do
49                do
50                        counts += w[ii][jj]
51                        err2_sum += data_err[ii][jj]*data_err[ii][jj]
52                        jj+=1
53                while(jj<=y2)
54                jj=y1
55                ii+=1
56        while(ii<=x2)
57       
58        err2_sum = sqrt(err2_sum)
59        ct_err = err2_sum
60       
61//      Print "error = ",err2_sum
62//      Print "error/counts = ",err2_sum/counts
63       
64       
65        Return (counts)
66End
67
68
69//from a marquee selection:
70//calculates the sum of counts in the box, and records this count value in
71//globals for both patch and Trans routines, then records the box coordinates
72//and the count value for the box in the header of the (currently displayed)
73//empty beam file for use later when calculating transmissions of samples
74//these values are written to unused (analysis) fields in the data header
75//4 integers and one real value are written
76//
77//re-written to work with Transmission panel as well as PatchPanel
78//which now bothe work from the same folder (:Patch)
79//
80Function SetXYBoxCoords() :  GraphMarquee
81
82        GetMarquee left,bottom
83        if(V_flag == 0)
84                Abort "There is no Marquee"
85        Endif
86        SVAR dispType=root:myGlobals:gDataDisplayType
87        if(cmpstr(dispType,"SAM")!=0)
88                DoAlert 0, "You can only use SetXYBox on SAM data files"
89                return(1)
90        endif
91        //printf "marquee left in bottom axis terms: %g\r",round(V_left)
92        //printf "marquee right in bottom axis terms: %g\r",round(V_right)
93        //printf "marquee top in left axis terms: %g\r",round(V_top)
94        //printf "marquee bottom in left axis terms: %g\r",round(V_bottom)
95        Variable x1,x2,y1,y2
96        x1 = round(V_left)
97        x2 = round(V_right)
98        y1 = round(V_bottom)
99        y2 = round(V_top)
100       
101        KeepSelectionInBounds(x1,x2,y1,y2)
102       
103        //check to make sure that Patch and Trans data folders exist for writing of global variables
104        If( ! (DataFolderExists("root:myGlobals:Patch"))  )
105                Execute "InitializePatchPanel()"
106        Endif
107        //check to make sure that Patch and Trans data folders exist for writing of global variables
108        If( ! (DataFolderExists("root:myGlobals:TransHeaderInfo"))  )
109                Execute "InitializeTransPanel()"
110        Endif
111       
112        //write string as keyword-packed string, to use IGOR parsing functions
113        String msgStr = "X1="+num2str(x1)+";"
114        msgStr += "X2="+num2str(x2)+";"
115        msgStr += "Y1="+num2str(y1)+";"
116        msgStr += "Y2="+num2str(y2)+";"
117        String/G root:myGlobals:Patch:gPS3 = msgStr
118        String/G root:myGlobals:Patch:gEmpBox = msgStr
119        //changing this global wil update the display variable on the TransPanel
120        String/G root:myGlobals:TransHeaderInfo:gBox = msgStr
121       
122        //sum the counts in the patch - working on the SAM data, to be sure that it's normalized
123        //to the same monitor counts and corrected for detector deadtime
124        String type = "SAM"
125        Variable counts,ct_err
126        counts = SumCountsInBox(x1,x2,y1,y2,ct_err,type)
127//      Print "marquee counts =",counts
128//      Print "relative error = ",ct_err/counts
129       
130        //Set the global gTransCts
131        Variable/G root:myGlobals:Patch:gTransCts = counts
132       
133        //now change the extra variables in the empty beam file
134        //get the filename from the SAM folder (there will only be one file)
135        SVAR partialName = root:Packages:NIST:SAM:FileList
136        //construct valid filename, then prepend path
137        String tempName = FindValidFilename(partialName)
138        Print "in marquee",partialName
139        //Print tempName
140        if(cmpstr(tempName,"")==0)
141                //file not found, get out
142                Abort "file not found, marquee"
143        Endif
144        //name is ok, prepend path to tempName for read routine
145        PathInfo catPathName
146        if (V_flag == 0)
147                //path does not exist - no folder selected
148                Abort "no path selected"
149        else
150                String filename = S_path + tempName
151        endif
152       
153        if(cmpstr(filename,"no file selected")==0)
154                Abort "no file selected"
155        Endif
156       
157        WriteXYBoxToHeader(filename,x1,x2,y1,y2)
158       
159        Print counts, " counts in XY box"
160        WriteBoxCountsToHeader(filename,counts)
161       
162        WriteBoxCountsErrorToHeader(filename,ct_err)
163       
164        return(0)
165End
166
167//finds the beam center (the centroid) of the selected region
168//and simply prints out the results to the history window
169//values are printed out in detector coordinates, not IGOR coords.
170//
171Function FindBeamCenter() :  GraphMarquee
172
173        //get the current displayed data (so the correct folder is used)
174        SVAR cur_folder=root:myGlobals:gDataDisplayType
175        String dest = "root:Packages:NIST:" + cur_folder
176       
177        Variable xzsum,yzsum,zsum,xctr,yctr
178        Variable left,right,bottom,top,ii,jj,counts
179       
180        // data wave is hard-wired in as the displayed data
181        NVAR dataIsLog=$(dest + ":gIsLogScale")         //check for log-scaling in current data folder
182        if (dataIsLog)
183                wave data=$(dest + ":linear_data")
184        else
185                wave data=$(dest + ":data")
186        endif
187       
188        GetMarquee left,bottom
189        if(V_flag == 0)
190                Print "There is no Marquee"
191        else
192                left = round(V_left)            //get integer values for selection limits
193                right = round(V_right)
194                top = round(V_top)
195                bottom = round(V_bottom)
196               
197                KeepSelectionInBounds(left,right,bottom,top)
198               
199                // selection valid now, calculate beamcenter
200                xzsum = 0
201                yzsum = 0
202                zsum = 0
203                // count over rectangular selection, doing each row, L-R, bottom to top
204                ii = bottom -1
205                do
206                        ii +=1
207                        jj = left-1
208                        do
209                                jj += 1
210                                counts = data[jj][ii]
211                                xzsum += jj*counts
212                                yzsum += ii*counts
213                                zsum += counts
214                        while(jj<right)
215                while(ii<top)
216               
217                xctr = xzsum/zsum
218                yctr = yzsum/zsum
219               
220                // add 1 to each to get to detector coordinates (1,128)
221                // rather than the data array which is [0,127]
222                xctr+=1
223                yctr+=1
224               
225                Print "X-center (in detector coordinates) = ",xctr
226                Print "Y-center (in detector coordinates) = ",yctr
227        endif
228       
229        //back to root folder (redundant)
230        SetDataFolder root:
231       
232End
233
234//still need to error check - out-of-bounds...waves exist.
235// allows a 2D Gaussian fit to a selected region of data in a SANS_Data window
236//puts up a new graph with the fitted contour
237Function Do_2D_Gaussian_Fit() :  GraphMarquee
238        String topWin=WinName(0,1)              //top *graph* window
239        //exit nicely if not in the Data display window
240        if(cmpstr(topWin,"SANS_Data") != 0)
241                DoAlert 0,"2D Gaussian fitting is only available from the Data Display Window"
242                return(1)
243        Endif
244       
245        GetMarquee/K left,bottom
246        Variable x1,x2,y1,y2,qxlo,qxhi,qylo,qyhi
247        if(V_flag == 0)
248                Print "There is no Marquee"
249        else
250                String junk="",df=""
251               
252                //**hard-wired info about the x-y q-scales
253                qxlo = DimOffset(root:myGlobals:q_x_axis,0)
254                qxhi = DimDelta(root:myGlobals:q_x_axis,0) + qxlo
255//              Print "qxlo,qxhi = ",qxlo,qxhi
256                Wave w=$"root:myGlobals:q_y_axis"
257                qylo=w[0]
258                qyhi=w[1]
259//              print "qylo,qyhi = ",qylo,qyhi
260               
261                junk=ImageInfo("SANS_Data","data",0)
262                df=StringByKey("ZWAVEDF", junk,":",";")
263//              print df
264                Duplicate/O $(df+"data") data,data_err
265                data_err=sqrt(data)             //for weighting
266               
267// comment out the SetScale lines if you want the result in terms of pixels as a way of
268// measuring the beam center. Note that you need to ADD ONE to fitted x0 and y0 to get detector
269// coordinates rather than the zero-indexed array. 2D fitting does have the benefit of
270// reporting error bars on the xy (if you believe that 2D gaussian is correct)         
271                SetScale/I x qxlo,qxhi,"",data
272                SetScale/I y qylo,qyhi,"",data
273               
274                Display /W=(10,50,361,351) /K=1
275                AppendImage data
276                ModifyImage data ctab= {*,*,Grays,0}
277                ModifyGraph width={Plan,1,bottom,left}
278                ModifyGraph mirror=2
279                ModifyGraph lowTrip=1e-04
280                ModifyImage data cindex=$"root:myGlobals:NIHColors"
281                SVAR/Z angst = root:Packages:NIST:gAngstStr
282                Label bottom "Qx ("+angst+"\\S-1\\M)"
283                Label left "Qy ("+angst+"\\S-1\\M)"
284
285                //keep selection in-bounds
286                x1=V_left
287                x2=V_right
288                y1=V_bottom
289                y2=V_top
290                KeepSelectionInBounds(x1,x2,y1,y2)
291
292                //cross correlation coefficent (K6) must be between 0 and 1, need constraints
293                Make/O/T/N=2 temp_constr
294                temp_constr = {"K6>0","K6<1"}
295               
296                CurveFit/N Gauss2D data[x1,x2][y1,y2] /I=1 /W=data_err /D /R /A=0 /C=temp_constr
297               
298                Killwaves/Z temp_constr
299        endif
300End
301
302// to save the image, simply invoke the IGOR menu item for saving graphics
303//
304Function SaveSANSGraphic() : GraphMarquee
305       
306        NVAR isDemoVersion=root:myGlobals:isDemoVersion
307        if(isDemoVersion==1)
308                //      comment out in DEMO_MODIFIED version, and show the alert
309                DoAlert 0,"This operation is not available in the Demo version of IGOR"
310        else
311                DoAlert 1,"Do you want the controls too?"
312                if(V_flag==1)
313                        GetMarquee/K/Z
314                        SavePICT /E=-5/SNAP=1
315                else
316                        DoIGORMenu "File","Save Graphics"
317                endif
318        endif
319End
320
321//does a sum over each of the files in the list over the specified range
322// x,y are assumed to already be in-bounds of the data array
323// output is dumped to the command window
324//
325Function DoBoxSum(fileStr,x1,x2,y1,y2,type)
326        String fileStr
327        Variable x1,x2,y1,y2
328        String type
329       
330        //parse the list of file numbers
331        String fileList="",item="",pathStr="",fullPath=""
332        Variable ii,num,err,cts,ct_err
333       
334        PathInfo catPathName
335        If(V_Flag==0)
336                Abort "no path selected"
337        Endif
338        pathStr = S_Path
339       
340        fileList=ParseRunNumberList(fileStr)
341        num=ItemsInList(fileList,",")
342       
343        //loop over the list
344        //add each file to SAM (to normalize to monitor counts)
345        //sum over the box
346        //print the results
347        Make/O/N=(num) FileID,BoxCounts,BoxCount_err
348        Print "Results are stored in root:FileID and root:BoxCounts waves"
349        for(ii=0;ii<num;ii+=1)
350                item=StringFromList(ii,fileList,",")
351                FileID[ii] = GetRunNumFromFile(item)            //do this here, since the list is now valid
352                fullPath = pathStr+item
353                ReadHeaderAndData(fullPath)
354//              String/G root:myGlobals:gDataDisplayType="RAW"
355//              fRawWindowHook()
356                if(cmpstr(type,"SAM")==0)
357                        err = Raw_to_work("SAM")
358                endif
359                String/G root:myGlobals:gDataDisplayType=type
360                fRawWindowHook()
361                cts=SumCountsInBox(x1,x2,y1,y2,ct_err,type)
362                BoxCounts[ii]=cts
363                BoxCount_err[ii]=ct_err
364                Print item+" counts = ",cts
365        endfor
366       
367        DoBoxGraph(FileID,BoxCounts,BoxCount_err)
368       
369        return(0)
370End
371
372Function DoBoxGraph(FileID,BoxCounts,BoxCount_err)
373        Wave FileID,BoxCounts,BoxCount_err
374       
375        Sort FileID BoxCounts,FileID            //sort the waves, in case the run numbers were entered out of numerical order
376       
377        Display /W=(5,44,383,306) BoxCounts vs FileID
378        ModifyGraph mode=4
379        ModifyGraph marker=8
380        ModifyGraph grid=2
381        ModifyGraph mirror=2
382        ErrorBars/T=0 BoxCounts Y,wave=(BoxCount_err,BoxCount_err)
383        Label left "Counts (per 10^8 monitor counts)"
384        Label bottom "Run Number"
385       
386        return(0)
387End
388
389//
390// promts the user for a range of file numbers to perform the sum over
391// list must be comma delimited numbers (or dashes) just as in the BuildProtocol panel
392// the (x,y) range is already selected from the marquee
393//
394Function BoxSum() :  GraphMarquee
395        GetMarquee left,bottom
396        if(V_flag == 0)
397                Abort "There is no Marquee"
398        Endif
399        SVAR dispType=root:myGlobals:gDataDisplayType
400        if(cmpstr(dispType,"RealTime")==0)
401                Print "Can't do a BoxSum for a RealTime file"
402                return(1)
403        endif
404        Variable x1,x2,y1,y2
405        x1 = V_left
406        x2 = V_right
407        y1 = V_bottom
408        y2 = V_top
409        KeepSelectionInBounds(x1,x2,y1,y2)
410       
411        String fileStr="",msgStr="Enter a comma-delimited list of run numbers, use dashes for ranges"
412        String type="RAW"
413        Prompt fileStr,msgStr
414        Prompt type,"RAW or Normalized (SAM)",popup,"RAW;SAM;"
415        DoPrompt "Pick the file range",fileStr,type
416        Print "fileStr = ",fileStr
417        printf "(x1,x2) (y1,y2) = (%d,%d) (%d,%d)\r",x1,x2,y1,y2
418       
419        DoBoxSum(fileStr,x1,x2,y1,y2,type)
420       
421        return(0)
422End     
423
424//function that keeps the marquee selection in the range [0,127] inclusive
425// (igor coordinate system)
426// uses pass-by reference!
427//
428// x1 = left
429// x2 = right
430// y1 = bottom
431// y2 = top
432//
433// accepts any detector size
434Function KeepSelectionInBounds(x1,x2,y1,y2)
435        Variable &x1,&x2,&y1,&y2
436       
437        NVAR pixelsX = root:myGlobals:gNPixelsX
438        NVAR pixelsY = root:myGlobals:gNPixelsY
439       
440        //keep selection in-bounds
441        x1 = (round(x1) >= 0) ? round(x1) : 0
442        x2 = (round(x2) <= (pixelsX-1)) ? round(x2) : (pixelsX-1)
443        y1 = (round(y1) >= 0) ? round(y1) : 0
444        y2 = (round(y2) <= (pixelsY-1)) ? round(y2) : (pixelsY-1)
445        return(0)
446End
447
448//testing function, not used
449Function testKeepInBounds(x1,x2,y1,y2)
450        Variable x1,x2,y1,y2
451       
452        KeepSelectionInBounds(x1,x2,y1,y2)
453        Print x1,x2,y1,y2
454        return(0)
455End
456
457Function SANS_Histogram_Pair() :  GraphMarquee
458        GetMarquee left,bottom
459        if(V_flag == 0)
460                Abort "There is no Marquee"
461        endif
462       
463        Cursor/W=SANS_Data/F/I A data 64,64
464        Cursor/M/S=2/H=1/L=0/C=(3,52428,1) A
465                               
466        // if cursor A on graph
467        // Do histogram pair
468        Variable aExists= strlen(CsrInfo(A)) > 0        // A is a name, not a string
469        if(aExists)
470                DoHistogramPair(hcsr(A),vcsr(A))
471        else
472                DoHistogramPair(64,64)
473        endif
474        return(0)
475       
476        //
477End
478// generates a histogram of the data as defined by the marquee. The longer dimension of the marquee
479// becomes the x-axis of the histogram (this may need to be changed for some odd case). Pixel range specified
480// by the marquee is inclusive, and is automatically kept in-bounds
481//
482// The counts over the (short) dimension are averaged, and plotted vs. the pixel position.
483// Pixel position is reported as Detector coordinates (1,128). Counts are whatever the current display
484// happens to be.
485//
486Function SANS_Histogram() :  GraphMarquee
487        GetMarquee left,bottom
488        if(V_flag == 0)
489                Abort "There is no Marquee"
490        endif
491//      // if cursor A on graph
492//      // Do histogram pair
493//      Variable aExists= strlen(CsrInfo(A)) > 0        // A is a name, not a string
494//      if(aExists)
495//              DoHistogramPair(hcsr(A),vcsr(A))
496//              return(0)
497//      endif
498        //
499        Variable count,x1,x2,y1,y2,xwidth,ywidth,vsX=1,xx,yy
500        x1 = V_left
501        x2 = V_right
502        y1 = V_bottom
503        y2 = V_top
504        KeepSelectionInBounds(x1,x2,y1,y2)
505        Print "x1,x2,y1,y2 (det) =",x1+1,x2+1,y1+1,y2+1
506        //determine whether to do x vs y or y vs x
507        xwidth=x2-x1
508        ywidth=y2-y1
509        if(xwidth < ywidth)
510                vsX=0           //sum and graph vs Y
511        endif
512        SVAR cur_folder=root:myGlobals:gDataDisplayType
513        WAVE data=$("root:Packages:NIST:"+cur_folder+":data")           //don't care if it's log or linear scale
514        Make/O/N=(max(xwidth,ywidth)+1) Position,AvgCounts
515        AvgCounts=0
516        //set position wave
517        if(vsX)
518                position=p+x1
519        else
520                position=p+y1
521        endif
522        //convert the position to Detector coordinates
523        position += 1
524       
525        //Compute the histogram (manually)
526        if(vsX)
527                for(xx=x1;xx<=x2;xx+=1)         //outer loop is the "x-axis"
528                        for(yy=y1;yy<=y2;yy+=1)
529                                AvgCounts[xx-x1] += data[xx][yy]
530                        endfor
531                endfor
532                AvgCounts /= (ywidth+1)
533        else
534                for(yy=y1;yy<=y2;yy+=1)
535                        for(xx=x1;xx<=x2;xx+=1)
536                                AvgCounts[yy-y1] += data[xx][yy]
537                        endfor
538                endfor
539                AvgCounts /= (xwidth+1)
540        endif
541        GetMarquee/K            //to keep from drawing the marquee on the new histo graph
542        //draw the graph, or just bring to the front with the new data
543        DoWindow/F SANS_Histo
544        if(V_Flag != 1)
545                Draw_Histo()
546        endif
547       
548        return(0)
549End
550
551//draws the histogram of the 2d data as specified by AvgCounts and Position
552//both wave are assumed to exist in the data folder. The SANS_Histogram() marquee
553//operation is responsible for creating them.
554//
555Function Draw_Histo()
556        Display /W=(197,329,567,461)/K=1 AvgCounts vs Position
557        DoWindow/C SANS_Histo
558        DoWindow/T SANS_Histo,"Histogram"
559        ModifyGraph mode=0,grid=1,mirror=2
560        ModifyGraph rgb=(21845,21845,21845)
561        ModifyGraph standoff=0
562        ModifyGraph hbFill=2
563        ModifyGraph useNegPat=1
564        ModifyGraph usePlusRGB=1
565        ModifyGraph useNegRGB=1
566        ModifyGraph hBarNegFill=2
567        ModifyGraph negRGB=(0,0,65535)
568        SetAxis/A/N=2 left
569        Label left "Counts"
570        Label bottom "Pixel (detector coordinates)"
571End
572
573//function will print marquee coordinates in axis terms, not detector terms
574//since IGOR is [0][127] and detector is (1,128)
575Function PrintMarqueeCoords() :  GraphMarquee
576        GetMarquee left,bottom
577        if(V_flag == 0)
578                Print "There is no Marquee"
579        else
580                Variable count,x1,x2,y1,y2,ct_err
581                x1 = V_left
582                x2 = V_right
583                y1 = V_bottom
584                y2 = V_top
585                printf "marquee left in bottom axis terms: %g\r",round(V_left)
586                printf "marquee right in bottom axis terms: %g\r",round(V_right)
587                printf "marquee bottom in left axis terms: %g\r",round(V_bottom)
588                printf "marquee top in left axis terms: %g\r",round(V_top)
589                printf "**note that you must add 1 to each axis coordinate to get detector coordinates\r"
590               
591                KeepSelectionInBounds(x1,x2,y1,y2)
592                SVAR cur_folder=root:myGlobals:gDataDisplayType
593                count = SumCountsInBox(x1,x2,y1,y2,ct_err,cur_folder)
594                Print "counts = ",count
595                Print "err/counts = ",ct_err/count
596        endif
597End
598
599//
600//
601// The histogram could potentially be done with the Igor ImageLineProfile operation
602// but this seems to be just as efficient to do.
603//
604Function DoHistogramPair(xin,yin)
605        Variable xin,yin
606       
607        Variable count,x1,x2,y1,y2,xwidth,ywidth,pt1,pt2,xx,yy
608        SVAR cur_folder=root:myGlobals:gDataDisplayType
609        WAVE data=$("root:Packages:NIST:"+cur_folder+":data")           //don't care if it's log or linear scale
610       
611
612        pt1 = 1         // extent along the "long" direction of the swath
613        pt2 = 128
614               
615        Make/O/D/N=(pt2-pt1+1) PositionX,AvgCountsX
616        Make/O/D/N=(pt2-pt1+1) PositionY,AvgCountsY
617        AvgCountsX=0
618        AvgCountsY=0
619       
620        //set position wave, in detector coordinates
621        positionX=p+pt1
622        positionY=p+pt1
623       
624        //do the vertical, then the horizontal
625        ControlInfo/W=HistoPair setvar0
626//      Print "width = ",V_Value
627        xwidth = V_Value                //+ -
628        ywidth = V_Value
629        x1 = xin - xwidth
630        x2 = xin + xwidth
631        y1 = pt1
632        y2 = pt2
633       
634        KeepSelectionInBounds(x1,x2,y1,y2)
635//      Print "x1,x2,y1,y2 (det) =",x1+1,x2+1,y1+1,y2+1
636       
637        //Compute the histogram (manually)
638        for(yy=y1;yy<=y2;yy+=1)
639                for(xx=x1;xx<=x2;xx+=1)
640                        AvgCountsY[yy-y1] += data[xx][yy]
641                endfor
642        endfor
643        AvgCountsY /= (xwidth+1)
644
645        // now do the Y
646        y1 = yin - ywidth
647        y2 = yin + ywidth
648        x1 = pt1
649        x2 = pt2
650               
651        KeepSelectionInBounds(x1,x2,y1,y2)
652//      Print "x1,x2,y1,y2 (det) =",x1+1,x2+1,y2+1,y2+1
653        for(xx=x1;xx<=x2;xx+=1)         //outer loop is the "x-axis"
654                for(yy=y1;yy<=y2;yy+=1)
655                        AvgCountsX[xx-x1] += data[xx][yy]
656                endfor
657        endfor
658        AvgCountsX /= (ywidth+1)
659       
660        GetMarquee/K            //to keep from drawing the marquee on the new histo graph
661        //draw the graph, or just bring to the front with the new data
662        Variable step = 30
663        DoWindow/F HistoPair
664        if(V_Flag != 1)
665                Draw_HistoPair()
666        else
667                SetAxis bottom (xin-step),(xin+step)
668                SetAxis bottomY (yin-step),(yin+step)
669        endif
670               
671        // so use a pair of cursors instead (how do I easily get rid of them?) - a "done" button
672        DoWindow SANS_Data
673        if(V_flag)
674                Cursor/W=SANS_Data/K B
675                Cursor/W=SANS_Data/K C
676       
677                Cursor/W=SANS_Data/F/I B data (xin-xwidth), (yin-yWidth)
678                Cursor/W=SANS_Data/M/S=2/H=1/L=1/C=(3,52428,1) B
679        //      Cursor/W=SANS_Data/M/A=0 B
680               
681                Cursor/W=SANS_Data/F/I C data (xin+xwidth), (yin+yWidth)
682                Cursor/W=SANS_Data/M/S=2/H=1/L=1/C=(3,52428,1) C
683        //      Cursor/W=SANS_Data/M/A=0 C
684        endif   
685        return(0)
686end
687
688
689Function Draw_HistoPair()
690        PauseUpdate; Silent 1           // building window...
691        Display /W=(432.75,431.75,903,698.75)/K=2 AvgCountsX vs PositionX as "Histogram Pair"
692        AppendToGraph/L=leftY/B=bottomY AvgCountsY vs PositionY
693        DoWindow/C HistoPair
694       
695        ModifyGraph rgb(AvgCountsX)=(21845,21845,21845)
696        ModifyGraph hbFill(AvgCountsX)=2
697        ModifyGraph useNegPat(AvgCountsX)=1
698        ModifyGraph usePlusRGB(AvgCountsX)=1
699        ModifyGraph useNegRGB(AvgCountsX)=1
700        ModifyGraph hBarNegFill(AvgCountsX)=2
701        ModifyGraph negRGB(AvgCountsX)=(0,0,65535)
702        ModifyGraph grid(left)=1,grid(bottom)=1,grid(leftY)=1
703        ModifyGraph mirror(left)=2,mirror(bottom)=2,mirror(leftY)=2
704        ModifyGraph standoff(left)=0,standoff(bottom)=0,standoff(leftY)=0
705        ModifyGraph lblPos(left)=62,lblPos(bottom)=39
706        ModifyGraph freePos(leftY)=0
707        ModifyGraph freePos(bottomY)={0,leftY}
708        ModifyGraph axisEnab(left)={0,0.4}
709        ModifyGraph axisEnab(leftY)={0.6,1}
710        Label left "Counts"
711        Label bottom "Pixel (detector coordinates)"
712        SetAxis/A/N=0 left
713        TextBox/C/N=text0/X=5.0/Y=5.0 "TOP"
714        TextBox/C/N=text0_1/X=5.0/Y=67.0 "RIGHT"
715        TextBox/C/N=text0_2/X=84.0/Y=67.0 "LEFT"
716        TextBox/C/N=text0_3/X=84.0/Y=5.0 "BOTTOM"
717       
718        ControlBar 40
719//      CheckBox check0,pos={300,11},size={72,14},proc=SH_FreeCursorCheck,title="Free Cursor"
720//      CheckBox check0,value= 0
721        Button button0 title="Update",size={70,20},pos={200,9},proc=SH_RecalcButton
722        SetVariable setvar0,pos={20,11},size={120,16},title="Width (pixels)"
723        SetVariable setvar0,limits={0,64,1},value= _NUM:5,proc=SH_WidthSetVarProc
724       
725        Button button1 title="Done",size={70,20},pos={300,9},proc=SH_DoneButton
726
727EndMacro
728
729
730// not used, just for testing
731Function CursorForHistogram()
732
733        Wave w=root:Packages:NIST:RAW:RealsRead
734       
735        Cursor/W=SANS_Data/F/I A data w[16],w[17]
736        Cursor/M/S=2/H=1/L=0/C=(3,52428,1) A
737       
738End
739
740
741Function SH_FreeCursorCheck(cba) : CheckBoxControl
742        STRUCT WMCheckboxAction &cba
743
744        switch( cba.eventCode )
745                case 2: // mouse up
746                        Variable checked = cba.checked
747                       
748                                //don't move the cursor
749                                //Cursor/W=SANS_Data/F/M/S=2/H=1/L=0/C=(3,52428,1) A data
750                       
751                                Cursor/W=SANS_Data/F/I A data 64,64
752                                Cursor/M/S=2/H=1/L=0/C=(3,52428,1) A
753                               
754                        break
755        endswitch
756
757        return 0
758End
759
760Function SH_RecalcButton(ba) : ButtonControl
761        STRUCT WMButtonAction &ba
762
763        switch( ba.eventCode )
764                case 2: // mouse up
765                        // click code here
766//                      Print "at = ",hcsr(A,"SANS_Data"),vcsr(A,"SANS_Data")
767                        DoHistogramPair(hcsr(A,"SANS_Data"),vcsr(A,"SANS_Data"))
768                        break
769        endswitch
770
771        return 0
772End
773
774Function SH_DoneButton(ba) : ButtonControl
775        STRUCT WMButtonAction &ba
776
777        switch( ba.eventCode )
778                case 2: // mouse up
779                        // click code here
780                        DoWindow/K HistoPair
781                       
782                        DoWindow SANS_Data
783                        if(V_flag)
784                                Cursor/W=SANS_Data/K A
785                                Cursor/W=SANS_Data/K B
786                                Cursor/W=SANS_Data/K C
787                        endif
788                        break
789        endswitch
790
791        return 0
792End
793Function SH_WidthSetVarProc(sva) : SetVariableControl
794        STRUCT WMSetVariableAction &sva
795
796        switch( sva.eventCode )
797                case 1: // mouse up
798                case 2: // Enter key
799                case 3: // Live update
800                        Variable dval = sva.dval
801                        String sval = sva.sval
802                       
803                        DoHistogramPair(hcsr(A,"SANS_Data"),vcsr(A,"SANS_Data"))
804                       
805                        break
806        endswitch
807
808        return 0
809End
Note: See TracBrowser for help on using the repository browser.