Python framework: pyana

The most pain-free way to access LCLS XTC data files from python is through LCLS's python framework, pyana. It is a non-interactive framework, but to some extent you can work interactively with the data it produces:

• The XTC Explorer gives you an "interactive" way to configure your analysis.

• Embedded IPython can give you interactive access to the data at regular intervals throughout your analysis.

The LCLS HDF5 data files can be worked with truely interactively from python, for details see How to access HDF5 data from Python. Be aware, though, that the framework (psana or pyana) does the job of syncronizing event data for you, and the lack of syncronization of arrays in the HDF5 files is the biggest drawback of working on datafiles outside of our framework.

Interactive data analysis with IPython

• 'IPython' (http://ipython.org/) is an enhanced python shell for interactive use. Many of the examples here would work equally well with a 'regular' python shell.
• plotting is done with 'matplotlib' (http://matplotlib.sourceforge.net/)
• If you're looking for an IDE to work with, consider 'Spyder' (http://code.google.com/p/spyderlib/).

Interactive access to XTC data

You cannot currently open an IPython or python shell and read in XTC data. Such an interactive framework for LCLS data is currently being explored, but does not exist yet.

What you can do is load data arrays into an interactive ipython session. Since you cannot load the XTC file directly into ipython, you'll need to run pyana or similar to create the arrays first. Or, you also have the option loading arrays from the LCLS HDF5 files, but be aware that when you work with HDF5 files, arrays from different sources (detectors) may not be synchronized. You will need to time-order and/or synchronize them yourself if you want to correlate data from different sources! See e.g. How to access HDF5 data from Python for how to do this.

Benefits of working directly with the XTC files is that they are available immediately, and you can start analyzing before the run is done collecting. The benefits of using the LCLS framework(s) is that each event is easily extracted and you don't have to worry about time-ordering or synchronizing data from different devices.

If you'd like to analyze XTC files with iPython, the options that exist are:

• ipython in combination with pyana. The XTC Explorer is an example of how this can be done (XtcExplorer#IPython).
  Note that currently there is no way to run pyana from IPython, but you can run a pyana job and launch ipython at the end to play with the plots/arrays.

• Write your own application, or be patient and wait for our interactive framework solution.

Because of this limitation, this page isn't really what it pretends to be ("How to access XTC data from Python"). But this page is a placeholder, and attempts to explore some of the functionality that we can use with XTC files later. The interactive analysis in (I)Python is the same in the end.

Numpy and HDF5 files

You can store numpy arrays from a pyana job (reads XTC) and store them in simple numpy files or HDF5 files. Here are some examples:

• saving and loading

  import numpy as np
  
  np.save("filename.npy", array)
  array = np.load("filename.npy")
  
  np.savetxt("filename.dat", array)
  array = loadtxt("filename.dat")
  

  This example shows saving and loading of a binary numpy file (.npy) and an ascii file (.dat).
  This only works with single arrays (max 2 dimensions).
  If you need to save multiple events/shots in the same file you will need to do some tricks (e.g. flatten the array and stack 1d arrays into 2d arrays where axis2 represent event number). Or you could save as an HDF5 file.

• saving to HDF5

  import h5py
  
  def beginjob(self,evt,env):
      self.ofile = h5py.File("outputfile.hdf5", 'w') # open for writing (overwrites existing file)
      self.shot_counter = 0
  
  def event(self,evt,env)
      # example: store several arrays from one shot in a group labeled with shot (event) number
      self.shot_counter += 1
      group = self.ofile.create_group("Shot%d" % self.shot_counter)
  
      image1_source = "CxiSc1-0|TM6740-1"
      image2_source = "CxiSc1-0|TM6740-2"
  
      frame = evt.getFrameValue(image1_source)
      image1 = frame.data()
      frame = evt.getFrameValue(image2_source)
      image2 = frame.data()
  
      dataset1 = group.create_dataset("%s"%image1_source,data=image1)
      dataset2 = group.create_dataset("%s"%image2_source,data=image2)
  
  def endjob(self,env)
      self.ofile.close()
  

  This example is shown in a pyana setting. The HDF5 file is declared and opened in beginjob, datasets created for each event, and the file is closed in the endjob method.
  Or you can group your datasets any other way you find useful, of course.

IPython used "like" MATLAB

Of course MATLAB is much more than this, but here's what we've started with. Here are some examples with IPython based on matlab functions provided by XPP. Thanks to H. Lemke for matlab examples and advice. A python module pymatlab.py defines a number of functions to use in this analysis example.

Starting an interactive session

[ofte@psana0XXX myrelease]$ ipython -pylab
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 pymatlab import *

Generally, it is recommended to load library modules with 'import pymatlab' and access all its methods and classes with pyamatlab.function. In an interactive session it may be easier to have access to the contents of pymatlab in your immediate workspace by doing 'from pymatlab import *'.

In [2]: who
H5getobjnames   ScanInput       ScanOutput      filtvec findmovingmotor
getSTDMEANfrac_from_startpoint  get_filter get_limits       get_limits_automatic
get_limits_channelhist  get_limits_correlation  get_limits_corrfrac
h5py    np      plt     rdXPPdata       runexpNO2fina       scan       scaninput

In [3]: whos
Variable                         Type          Data/Info
--------------------------------------------------------
H5getobjnames                    function      <function H5getobjnames at 0x2b57de8>
ScanInput                        type          <class 'pymatlab.ScanInput'>
ScanOutput                       type          <class 'pymatlab.ScanOutput'>
filtvec                          function      <function filtvec at 0x2b57f50>
findmovingmotor                  function      <function findmovingmotor at 0x2b57d70>
getSTDMEANfrac_from_startpoint   function      <function getSTDMEANfrac_<...>_startpoint at 0x2b581b8>
get_filter                       function      <function get_filter at 0x2b57ed8>
get_limits                       function      <function get_limits at 0x2b58050>
get_limits_automatic             function      <function get_limits_automatic at 0x2b58230>
get_limits_channelhist           function      <function get_limits_channelhist at 0x2b582a8>
get_limits_correlation           function      <function get_limits_correlation at 0x2b580c8>
get_limits_corrfrac              function      <function get_limits_corrfrac at 0x2b58140>
h5py                             module        <module 'h5py' from '/reg<...>ython/h5py/__init__.pyc'>
np                               module        <module 'numpy' from '/re<...>thon/numpy/__init__.pyc'>
plt                              module        <module 'matplotlib.pyplo<...>n/matplotlib/pyplot.pyc'>
rdXPPdata                        function      <function rdXPPdata at 0x2b57c80>
runexpNO2fina                    function      <function runexpNO2fina at 0x2b57e60>
scan                             ScanOutput    <pymatlab.ScanOutput object at 0x2b60536bee90>
scaninput                        ScanInput     <pymatlab.ScanInput object at 0x2b60536b4e90>

Like in MATLAB, who gives you a short list of workspace contents, whos gives you a longer list of workspace contents.

Plot filtered IPIMB data with limits from graphical input:

Here's a log from a session that produces a loglog plot (blue dots) of two IPIMB channels, selects limits from graphical inpu (mouse click),
draws the selected events with red dots.

In [3]: scaninput = ScanInput()

In [4]: scaninput.fina = "/reg/d/psdm/XPP/xpp23410/hdf5/xpp23410-r0107.h5"

In [5]: scan = rdXPPdata(scaninput)
Reading XPP data from  /reg/d/psdm/XPP/xpp23410/hdf5/xpp23410-r0107.h5
Found pv control object  fs2:ramp_angsft_target
Found scan vector  [ 2800120.  2800240.  2800360.  2800480.  2800600.  2800720.  2800840.
  2800960.  2801080.  2801200.  2801320.  2801440.  2801560.  2801680.
  2801800.  2801920.  2802040.  2802160.  2802280.  2802400.  2802520.
  2802640.  2802760.  2802880.  2803000.  2803120.  2803240.  2803360.
  2803480.  2803600.  2803720.  2803840.  2803960.  2804080.  2804200.
  2804320.  2804440.  2804560.  2804680.  2804800.  2804920.  2805040.
  2805160.  2805280.]
Fetching data to correlate with motor
['IPM1', 'IPM2']
(44, 120, 4)

In [6]: channels = np.concatenate(scan.scandata,axis=0)

In [7]: channels.shape
Out[7]: (5280, 4)

In [8]: get_limits(channels,1,"correlation")
4 channels a 5280 events
indexes that pass filter:  (array([   1,    5,    8, ..., 5266, 5272, 5273]),)
Out[8]:
array([[ 0.00086654,  0.01604564],
       [ 0.67172102,  0.71968567],
       [ 0.00194716,  0.01447819],
       [ 0.80365403,  0.73463468]])

In [9]: plt.draw()

Table of comparison (MATLAB vs MatPlotLib)

See also http://www.scipy.org/NumPy_for_Matlab_Users

%
%
%
a1 = subplot(121);
loglog(channels(:,1),channels(:,2),'o')
xlabel('CH0')
ylabel('CH1')
a2 = subplot(122);
loglog(channels(:,3),channels(:,4),'o')
xlabel('CH2')
ylabel('CH3')

import matplotlib.pyplot as plt
import numpy as np

a1 = plt.subplot(221)
plt.loglog(channels[:,0],channels[:,1], 'o' )
plt.xlabel('CH0')
plt.ylabel('CH1')
a2 = plt.subplot(222)
plt.loglog(channels[:,2],channels[:,3], 'o' )
plt.xlabel('CH2')
plt.ylabel('CH3')

axes(a1)
hold on
lims(1:2,:) = ginput(2);

axes(a2)
hold on
lims(3:4,:) = ginput(2);

lims = np.zeros((4,2),dtype="float")

plt.axes(a1)
plt.hold(True)
lims[0:2,:] = plt.ginput(2)

plt.axes(a2)
plt.hold(True)
lims[2:4,:] = plt.ginput(2)

fbool1 = (channels(:,1)>min(lims(1:2,1)))&(channels(:,1)<max(lims(1:2,1)))
fbool2 = (channels(:,2)>min(lims(1:2,2)))&(channels(:,2)<max(lims(1:2,2)));
fbool = fbool1&fbool2
loglog(channels(fbool,1),channels(fbool,2),'or')

fbool3 = (channels(:,3)>min(lims(3:4,3)))&(channels(:,3)<max(lims(3:4,3)))
fbool4 = (channels(:,4)>min(lims(3:4,4)))&(channels(:,4)<max(lims(3:4,4)));
fbool = fbool3&fbool4
loglog(channels(fbool,3),channels(fbool,4),'or') 

fbools0 = (channels[:,0]>lims[:,0].min())&(channels[:,0]<lims[:,0].max())
fbools1 = (channels[:,1]>lims[:,1].min())&(channels[:,1]<lims[:,1].max())
fbools = fbools0 & fbools1

fbools2 = (channels[:,2]>lims[:,2].min())&(channels[:,2]<lims[:,2].max())
fbools3 = (channels[:,3]>lims[:,3].min())&(channels[:,3]<lims[:,3].max())
fbools = fbools2&fbools3

Data visualization with NumPy (arrays) and MatPlotLib (plots).

This is not meant to be documentation or a tutorial for matplotlib or numpy. Just a place to document stuff that I have a hard time finding explained elsewhere.

• Inspecting objects

  for attr_name in dir(obj):
      attribute = getattr(obj, attr_name)
      print attr_name, ": ", attribute
  

  This is really just simple python. But since 'matplotlib' documentation can be frustratingly non-verbose about what functions and attributes are available for its various classes/objects, I found this is a useful way to inspect what an object knows about itself.