Tutorial - python, pyana and matplotlib

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.

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*

addpkg XtcExplorer
scons
xtcexplorer /reg/d/psdm/xpp/xppi0310/xtc/e81-r0098-s0*

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

# useful imports
import numpy as np
import matplotlib.pyplot as plt
from pypdsdata.xtc import TypeId

class mymodule (object) :
    """Class whose instance will be used as a user analysis module. """

    def __init__ ( self,
                   source = ""
                   threshold = "" ):
        """Class constructor.
        The parameters to the constructor are passed from pyana configuration file.
        If parameters do not have default values  here then the must be defined in
        pyana.cfg. All parameters are passed as strings, convert to correct type before use.

        @param source         name of device, format 'Det-ID|Dev-ID'
        @param threshold      threshold value (remember to convert from string)
        """
        self.source = source
        self.threshold = float(threshold)

    def beginjob( self, evt, env ) :
        """This method is called once at the beginning of the job. It should
        do a one-time initialization possible extracting values from event
        data (which is a Configure object) or environment.

        @param evt    event data object
        @param env    environment object
        """
        pass

    def beginrun( self, evt, env ) :
        """This optional method is called if present at the beginning of the new run.

        @param evt    event data object
        @param env    environment object
        """
        pass

    def begincalibcycle( self, evt, env ) :
        """This optional method is called if present at the beginning 
        of the new calibration cycle.

        @param evt    event data object
        @param env    environment object
        """
        pass

    def event( self, evt, env ) :
        """This method is called for every L1Accept transition.

        @param evt    event data object
        @param env    environment object
        """
        pass

    def endcalibcycle( self, env ) :
        """This optional method is called if present at the end of the 
        calibration cycle.
        
        @param env    environment object
        """
        pass
        
    def endrun( self, env ) :
        """This optional method is called if present at the end of the run.
        
        @param env    environment object
        """
        pass

    def endjob( self, env ) :
        """This method is called at the end of the job. It should do 
        final cleanup, e.g. close all open files.
        
        @param env    environment object
        """
        pass

Point detector delay scan

• Fetching the ControlPV information:
  ControlPV is available from the env object, and since it only changes at the beginning
  of each calibration cycle, the begincalibcycle function is the appropriate place to get it:
  

      def begincalibcycle( self, evt, env ) :
  

  
  The ControlConfig object may contain several pvControl and pvMonitor objects. In this case
  there's only one, but make sure the name matches anyway:

          ctrl_config = env.getConfig(TypeId.Type.Id_ControlConfig)
  
          for ic in range (0, ctrl_config.npvControls() ):
              cpv = ctrl_config.pvControl(ic)
              if cpv.name()=="fs2:ramp_angsft_target":
  
                  # store the value in a class variable (visible in every class method)
                  self.current_pv_value = cpv.value() )
  

• Fetching the IPIMB and PhaseCavity information:
  All the other information that we need, is available through the evt object, and
  event member function is the place to get it:
  

      def event( self, evt, env ) :
  

  
  Use "XppSb3Ipm-1|Ipimb-0" (a.k.a. IPM3) sum of all channels for normalization and filtering
  

          ipmN_raw = evt.get(TypeId.Type.Id_IpimbData, "XppSb3Ipm-1|Ipimb-0")
          ipmN_fex = evt.get(TypeId.Type.Id_IpmFex, "XppSb3Ipm-1|Ipimb-0")
  
          ipmN_norm = ipmN_fex.sum
  

  
  Use "XppSb3Pim-1|Ipimb-0" (a.k.a. PIM3) channel 1 as signal
  

          ipmS_raw = evt.get(TypeId.Type.Id_IpimbData, "XppSb3Pim-1|Ipimb-0" )
          ipmS_fex = evt.get(TypeId.Type.Id_IpmFex, "XppSb3Pim-1|Ipimb-0" )
  
          ipm_sig = ipmS_fex.channel[1]
  

  
  Get the phase cavity:
  

          pc = evt.getPhaseCavity()
          phasecav1 = pc.fFitTime1
          phasecav2 = pc.fFitTime2
          charge1 = pc.fCharge1
          charge2 = pc.fCharge2
  

  
  Compute delay time and fill histograms
  

          delaytime = self.current_pv_value + phasecav1*1e3    
  
          # The "histograms" are nothing but python lists. Append to them, and turn them into arrays at the end. 
          self.h_ipm_rsig.append( ipm_sig )
          self.h_ipm_nsig.append( ipm_sig/ipm_norm )
          self.h_delaytime.append( delaytime )
  

Image peak finding

http://docs.scipy.org/doc/scipy/reference/ndimage.html

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...

http://ipython.org/

[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:

• http://matplotlib.sourceforge.net/users/pyplot_tutorial.html

• 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.