A quick walk-through of the tools that exist for analysis of xtc files with python.
The main focus is on pyana, and the examples are from and for XPP primarily,
but may be useful examples to other experiments too.
outline of contents
The Basics
Python
http://docs.python.org/tutorial/
Pyana
Analysis Workbook. Python-based Analysis
Setting up a work directory (a.k.a. offline release directory)
Prior to this, you may need to set up your account for offline analysis:
Analysis Workbook. Account Setup
The general version of this is in Analysis Workbook. Quick Tour
newrel ana-current xpptutorial cd xpptutorial ls -l less .sit_release sit_setup
Exploring an xtc file
pyxtcreader
pyxtcreader -h usage: pyxtcreader [options] xtc-files ... options: -h, --help show this help message and exit -v, --verbose -l L1_OFFSET, --l1-offset=L1_OFFSET
Loops through the xtc datagrams and dumps info to screen. I recommend piping it to 'less'.
pyxtcreader /reg/d/psdm/xpp/xppi0310/xtc/e81-r0098-s0* | less
xtcscanner
xtcscanner -h usage: xtcscanner [options] xtc-files ... options: -h, --help show this help message and exit -n NDATAGRAMS, --ndatagrams=NDATAGRAMS -v, --verbose -l L1_OFFSET, --l1-offset=L1_OFFSET -e, --epics
Similar to pyxtcreader in that it loops throug xtc datagrams, but doesn't print to screen. Internally counts the datatypes it finds, and at the end dumps a summary only. Optinally prints out epics information (default no).
xtcscanner /reg/d/psdm/xpp/xppi0310/xtc/e81-r0098-s0*
You should see output similar to this:
Scanning.... Start parsing files: ['/reg/d/psdm/xpp/xppi0310/xtc/e81-r0098-s00-c00.xtc', '/reg/d/psdm/xpp/xppi0310/xtc/e81-r0098-s01-c00.xtc'] 14826 datagrams read in 4.120000 s . . . . . . . ------------------------------------------------------------- XtcScanner information: - 61 calibration cycles. - Events per calib cycle: [240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240, 240] Information from 1 control channels found: fs2:ramp_angsft_target Information from 11 devices found BldInfo:EBeam: EBeamBld_V1 (14641) BldInfo:FEEGasDetEnergy: FEEGasDetEnergy (14563) Any (78) BldInfo:NH2-SB1-IPM-01: SharedIpimb (14641) BldInfo:PhaseCavity: PhaseCavity (14641) DetInfo:EpicsArch-0|NoDevice-0: Epics_V1 (107580) DetInfo:NoDetector-0|Evr-0: EvrConfig_V4 (62) EvrData_V3 (14640) DetInfo:XppSb2Ipm-1|Ipimb-0: IpimbConfig_V1 (1) IpmFexConfig_V1 (1) IpimbData_V1 (14640) IpmFex_V1 (14640) DetInfo:XppSb3Ipm-1|Ipimb-0: IpimbConfig_V1 (1) IpmFexConfig_V1 (1) IpimbData_V1 (14640) IpmFex_V1 (14640) DetInfo:XppSb3Pim-1|Ipimb-0: IpimbConfig_V1 (1) IpmFexConfig_V1 (1) IpimbData_V1 (14640) IpmFex_V1 (14640) DetInfo:XppSb4Pim-1|Ipimb-0: IpimbConfig_V1 (1) IpmFexConfig_V1 (1) IpimbData_V1 (14640) IpmFex_V1 (14640) ProcInfo:: RunControlConfig_V1 (62) XtcScanner is done! -------------------------------------------------------------
xtcexplorer
XTC Explorer - GUI interface that builds pyana modules for you.
xtcexplorer /reg/d/psdm/xpp/xppi0310/xtc/e81-r0098-s0*
- Then hit the "Scan File(s)" button (can you find it?!?)
- What do you see? Compare the GUI that pops up, with the output in the terminal window.
- Checkmark the IPM3 checkbox ('XppSb3Ipm=1|Ipimb-0')
- hit the 'Write configuration to file' button,
- hit the 'Run pyana' button
- hit 'OK' and wait till a plot pops up... Close the window and wait again...
- hit the 'Quit pyana' button
- go to the 'General Settings' tab and switch Display mode to 'SlideShow'
- go back to 'General Settings' again and change the number of events to accumulate to 240
- hit the 'Write configuration to file'
- hit the 'Edit configuration file' button. Edit the line with 'quantities = ': remove 'fex:channels' and add 'fex:ch1' and 'fex:ch0' instead
- hit 'Run pyana' button again (as well as 'OK' when that pops up). Stare at the plot for a while...
addpkg XtcExplorer scons xtcexplorer /reg/d/psdm/xpp/xppi0310/xtc/e81-r0098-s0*
Now you have a local version of the XtcExplorer package in your directory. That allows you to edit the
source code and customize the analysis modules in the XtcExplorer/src
directory.
Exercise:
Edit XtcExplorer/src/pyana_ipimb.py
to make a loglog plot of channel1 vs channel0.
The xtcexplorer has several shortcomings. It tries to be very generic, and thus is sometimes slower than it would have to be. It is also currently only capable of plotting from a single device in each plot, so many correlation plots will need to be added by hand. However, it is a simple tool to just look at the data contents, and provides many examples through its pyana modules.
For some more useful analysis examples, in the following we'll stick to writing customized pyana modules and running pyana from the command line.
Extracting the data with pyana, some examples
Outline of a pyana module
Like the other frameworks, pyana is an executable that loops through the XTC file and calls all
requested user modules at certain transitions. All the analysts need to do is to fill in the
relevant functions in their user analysis module:
code
For the following two examples, check out the latest version of the pyana_examples
package:
addpkg pyana_examples HEAD scons
Point detector delay scan
Open an editor and save the following in a file named pyana.cfg:
[pyana] modules pyana_examples.xppt_delayscan [pyana_examples.xppt_delayscan] controlpv = fs2:ramp_angsft_target ipimb_norm = XppSb3Ipm-1|Ipimb-0 ipimb_sig = XppSb3Pim-1|Ipimb-0 threshold = 0 outputfile = point_scan_delay.npy
Run pyana (start with 200 events):
pyana -n 200 /reg/d/psdm/XPP/xppi0310/xtc/e81-r0098-s0*
Highlighting of some code snippets from xppt_delayscan.py
:
Image peak finding
Here are a collection of useful algorithms for image analysis: http://docs.scipy.org/doc/scipy/reference/ndimage.html
This particular example is done with a CSPad image, but only a single section is available. For more typical CSPad module, see next section.
- For each event, fetch the CsPad information, and get the image array:
def event( self, evt, env ) : elements = evt.getCsPadQuads(self.source, env) image = elements[0].data().reshape(185, 388)
- Select a region of interest. If none is given (optional module parameter), set RoI to be the whole image.
# Region of Interest (RoI) if self.roi is None: self.roi = [ 0, image.shape[1], 0, image.shape[0] ] # [x1,x2,y1,y2] print "ROI [x1, x2, y1, y2] = ", self.roi
- Using only the RoI subset of the image, compute center-of-mass using one of the SciPy.ndimamge algorithms. Add to it the position of the RoI to get the CMS in global pixel coordinates:
roi_array = image[self.roi[2]:self.roi[3],self.roi[0]:self.roi[1]] cms = scipy.ndimage.measurements.center_of_mass(roi_array) print "Center-of-mass of the ROI: (x, y) = (%.2f, %.2f)" %(self.roi[0]+cms[1],self.roi[2]+cms[0])
- Here's an example how you can make an interactive plot and select the Region of Interest with the mouse. Here we plot the image in two axes (subpads on the canvas). The first will always show the full image. In the second axes, you can select a rectangular region in "Zoom" mode (click on the Toolbar's Zoom button). The selected region will be drawn on top of the full image to the left, while the right plot will zoom into the selected region:
fig = plt.figure(1,figsize=(16,5)) axes1 = fig.add_subplot(121) axes2 = fig.add_subplot(122) axim1 = axes1.imshow(image) axes1.set_title("Full image") axim2 = axes2.imshow(roi_array, extent=(self.roi[0],self.roi[1],self.roi[3],self.roi[2])) axes2.set_title("Region of Interest") # rectangular ROI selector rect = UpdatingRect([0, 0], 0, 0, facecolor='None', edgecolor='red', picker=10) rect.set_bounds(*axes2.viewLim.bounds) axes1.add_patch(rect) # Connect for changing the view limits axes2.callbacks.connect('xlim_changed', rect) axes2.callbacks.connect('ylim_changed', rect)
- To compute the center-of-mass of the selected region, revert back to non-zoom mode (hit the 'zoom' button again) and click on the rectangle. The rectangle is connected to the 'onpick' function which updates
self.roi
and computes the center-of-mass:def onpick(event): xrange = axes2.get_xbound() yrange = axes2.get_ybound() self.roi = [ xrange[0], xrange[1], yrange[0], yrange[1]] roi_array = image[self.roi[2]:self.roi[3],self.roi[0]:self.roi[1]] cms = scipy.ndimage.measurements.center_of_mass(roi_array) print "Center-of-mass of the ROI: (x, y) = (%.2f, %.2f)" % (self.roi[0]+cms[1],self.roi[2]+cms[0]) fig.canvas.mpl_connect('pick_event', onpick) plt.draw()
CSPad images and tile arangements
Saving data arrays
saving numpy arrays to numpy file
import numpy as np myarray = np.arange(0,100) np.save( "output_file_name.npy", myarray) np.savetxt( "output_file_name.txt", myarray)
Both of these work with arrays of maximum 2 dimensions. And only one array per file.
saving to MatLab file
import scipy.io N_array = np.arange(0,100) x_array = np.random(100) y_array = np.random(100) scipy.io.savemat( "output_file_name.mat", mdict={'N': N_array, 'x' : x_array, 'y' : y_array } )
saving to HDF5 file
import h5py ofile = h5py.File("output_file_name.h5",'w') group = ofile.create_group("MyGroup") group.create_dataset('delaytime',data=np.array(self.h_delaytime)) group.create_dataset('rawsignal',data=np.array(self.h_ipm_rsig)) group.create_dataset('normsignal',data=np.array(self.h_ipm_nsig)) ofile.close()
For more examples, see How to access HDF5 data from Python and http://code.google.com/p/h5py/
Interactive analysis with IPython
The version available in our offline release system is
IPython 0.9.1 – An enhanced Interactive Python.
so this is the one I've been using in these examples.
Not a whole lot more than a python shell.
However, the latest IPython has loads of new and interesting features...
[ofte@psana0106 xpptutorial]$ ipython Python 2.4.3 (#1, Nov 3 2010, 12:52:40) Type "copyright", "credits" or "license" for more information. IPython 0.9.1 -- An enhanced Interactive Python. ? -> Introduction and overview of IPython's features. %quickref -> Quick reference. help -> Python's own help system. object? -> Details about 'object'. ?object also works, ?? prints more. In [1]: from numpy import * In [2]: from matplotlib.pyplot import * In [3]: ipm3 = load('point_scan_delay.npy') In [4]: ipm3.shape Out[4]: (200, 3) In [5]: ion() In [6]: delay = ipm3[:,0] In [7]: ipmraw = ipm3[:,1] In [8]: ipmnorm = ipm3[:,2] n [9]: plot(delay,ipmnorm,'ro') Out[9]: [<matplotlib.lines.Line2D object at 0x59c4c10>] In [10]: draw() In [11]:
Plotting with MatPlotLib
Matplotlib:
- The plotting can be done directly in the pyana module, but be aware that you need to disable plotting for the
module to run successfully in a batch job.import matplotlib.pyplot as plt plt.plot(array) plt.show()
- Or you can load arrays from a file and interactively plot them in iPython. The same ('recommended') syntax as above can be used, or if you use 'import *' you don't need to prepend the commands with the package name, which is handy when plotting interactively:
from matplotlib.pyplot import * ion() plot(array) draw()
Related useful tools and links
- http://www.scipy.org/NumPy_for_Matlab_Users
- http://www.scipy.org/
- http://www.sagemath.org/
- http://code.google.com/p/spyderlib/
Non-interactive batch analysis
Pyana jobs are designed to do batch analysis, but matplotlib plotting does not go well with this. If you want your job to produce graphics, make sure to use a matplotlib backend that writes the graphics directly to file, e.g. png files.
Multiprocessing
Pyana can make use of multiple core processing. On the command line, add the option '-p N' where N is the number of cores to use.
Extra care needs to be taken when plotting. Also, output files need to be made with the pyana mkfile command. The output will be merged at the end of the job, but may not be in order. So if you need events to be written to a file in chronological order, you're better off using single core processing.