Page History

Composition Setup
cloak.toggle.type = text cloak.toggle.open = + show cloak.toggle.close = - hide

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.

Toggle Cloak

id	toc

outline Outline of contents:

...

Table of Contents

...

maxLevel

...

3

The Basics

Python

http://docs.python.org/tutorial/^{Image Removed}

Pyana

Analysis Workbook. Python-based Analysis

...

Panel

title	Open a terminal at pslogin or psana, and type:

Code Block
newrel ana-current xpptutorial cd xpptutorial ls -lla less .sit_release sit_setup

...

Panel

title	Try:

Code Block

	none
	none

pyxtcreader /reg/d/psdm/xpp/xppi0310/xtc/e81-r0098-s0* | less

Try the same with different verbosity levels:

Code Block

	none
	none


pyxtcreader -v /reg/d/psdm/xpp/xppi0310/xtc/e81-r0098-s0* | less
pyxtcreader -vv /reg/d/psdm/xpp/xppi0310/xtc/e81-r0098-s0* | less

xtcscanner

Code Block

	none
	none

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

...

Code Block

	none
	none

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 - Old - GUI interface that builds pyana modules for you.

...

Panel

title	Try something else:

Code Block

	none
	none


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

('kinit' will ask you for your password and give you a Kerberos ticket valid for 24 h that you need to access our afs software repository)

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 for later:
Edit XtcExplorer/src/pyana_ipimb.py to make a loglog plot of channel1 vs channel0.

...

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:

Toggle Cloak

id	kappe

code

...

id	kappe

...

But before getting to the pyana modules, I'll briefly touch on a few items general to python that may be useful: saving files, matplotlib for plotting, and IPython for interactive work.

NumPy, SciPy and MatPlotLib

These are packages that you may want to look into. Pretty much all our examples here are using them:

http://numpy.scipy.org/ - Numerical package for python (arrays etc)
http://www.scipy.org/ - Scientific tools
http://matplotlib.sourceforge.net/ - plotting package

Saving data arrays

Here are a few examples of how you can save data arrays in python.

Panel

title	saving numpy arrays to numpy file

Code Block
import numpy as np myarray = np.arange(0,100) np.save( "output_file_name.npy", myarray) np.savetxt( "output_file_name.txt", myarray)

Panel

	saving to MatLab file
	saving to MatLab file

Code Block
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 } )

Panel

	saving to HDF5 file
	saving to HDF5 file

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:

Run pyana (start with 200 events):

Panel

Code Block

none


[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

Code Block

	none
	none

pyana -n 200 /reg/d/psdm/XPP/xppi0310/xtc/e81-r0098-s0*

Toggle Cloak

id	delayscan

Highlighting of some code snippets from xppt_delayscan.py:

...

id	delayscan

...

Image peak finding

Here are a collection of useful algorithms for image analysis: http://docs.scipy.org/doc/scipy/reference/ndimage.html^{Image Removed}

This particular example is done with a CSPad image, but only a single section is available. For more typical CSPad module, see next section.

...


   def event( self, evt, env ) :

        elements = evt.getCsPadQuads(self.source, env)
        image = elements[0].data().reshape(185, 388)

...

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/

Plotting with MatPlotLib

One of the most commonly used tools for plotting in python: matplotlib. Other alternatives exist too.

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.
Code Block
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:
Code Block
from matplotlib.pyplot import * ion() plot(array) draw()

Interactive analysis with IPython

The LCLS offline analysis group does have plans for a real interactive pyana, but currently this is not available.
2011-11-04 iPsana Interactive Analysis Framework.pdf

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/

Panel

title	Loading your arrays into (I)Python and plotting interactively:

This example reads in a file produced by the "point detector delay scan" example below.

Code Block


[ofte@psana0106 xpptutorial]$ 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]: ipm3 = load('point_scan_delay.npy')

In [2]: ipm3.shape
Out[2]: (200, 3)

In [3]: ion()

In [4]: delay = ipm3[:,0]

In [5]: ipmraw = ipm3[:,1]

In [6]: ipmnorm = ipm3[:,2]

In [7]: plot(delay,ipmnorm,'ro')
Out[7]: [<matplotlib.lines.Line2D object at 0x59c4c10>]

In [9]: draw()

In [10]:

Sometimes you need to issue the draw() command twice, for some reason. After drawing you can keep working on the arrays and plot more...

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:

Toggle Cloak

id	kappe

...

Cloak

id	kappe

Code Block


# 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

...


        # 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

...


        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])

...

"""

...

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

...

CSPad images and tile arangements

CSPad data structure

CSPad data in xtc is a list of elements. In pyana get the list from the evt (event) object (notice the need for the env (environment) object too!):

...


elements = evt.getCsPadQuads(self.source, env)

elements here is a python list of ElementV1 or ElementV2 (or later versions) objects, each representing one quadrant. The list is not ordered, so to know which quadrant you have, you have to check with element.quad(). To store a local array of the whole CSPad detector, you can do the following.

...

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

...

"""

...

What we have so far gives you a 4d numpy array of all pixels. And if you want to store it in e.g. a numpy array, you can reshape it down to 2 dimensions (this is the format of the official pedestal files made by the translator):

...


pixels = pixel_array.reshape(1480,388)
np.save("pixel_pedestal_file.npy", pixels )

pass

Cloak

Panel

For the following two examples, check out the latest version of the pyana_examples package:

Code Block

	none
	none


addpkg pyana_examples HEAD
scons

(Note, if you don't already have a Kerberos ticket, you need to issue a 'kinit' command before 'addpkg'. You will be prompted for your unix password.)

Datatypes, and how to find data from your detector/device in the xtc file.

Pyana and psana has follows this naming scheme for labeling the datatypes from various devices. You can find the
names by investigating the xtc file with the above-mentioned tools (pyxtcreader, xtcscanner, xtcexplorer).
To see some examples of how to fetch the various data types in pyana (or psana), look at Devices and Datatypes.

Point detector delay scan

The python code for this pyana module resides in pyana_examples/src/xppt_delayscan.py. In this example, we do a point detector delay scan, where we get the time as scan points via a control PV, and where time rebinning based on phase cavity measurement is used to improve the time resolution. One IPIMB device (a.k.a. IPM3) is used for normalization (i0, I Zero) (parameter name ipimb_norm) and another IPIMB device (a.k.a. PIM3) channel 1 is used as the signal (parameter name ipimb_sig).

Panel

Open an editor and save the following in a file named pyana.cfg:

Code Block

	none
	none


[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.1
outputfile = point_scan_delay.npy

If you look at the code (pyana_examples/src/xppt_delayscan.py) you'll notice there are no detector names in there. The names of the detectors in the XTC file are passed as parameters from the configuration file above. The ipimb_norm parameter represents the IPIMB diode used for normalization, and the configuration files set its value to "XppSb3Ipm-1|Ipimb-0" (a.k.a. IPM3). Similarly, the IPIMB diode used for signal is represented by the pipmb_sig and its value set to "XppSb3Pim-1|Ipimb-0" (a.k.a. PIM3). By changing these parameter values, the pyana_examples/src/xppt_delayscan.py module can easily be used for other experiments or instruments.

Run pyana (start with 200 events):

Code Block

	none
	none


pyana -n 200 /reg/d/psdm/XPP/xppi0310/xtc/e81-r0098-s0*

Toggle Cloak

id	delayscan

...

Cloak

id	delayscan

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:

Code Block

	none
	none


    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:

Code Block

	none
	none


        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:

Code Block

	none
	none


    def event( self, evt, env ) :

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

Code Block

	none
	none


        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

Code Block

	none
	none


        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:

Code Block

	none
	none


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

Compute delay time and fill histograms

Code Block

	none
	none


        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 )

Cloak

Image peak finding

Here are a collection of useful algorithms for image analysis: http://docs.scipy.org/doc/scipy/reference/ndimage.html

The python code for this pyana module example resides in pyana_examples/src/xppt_image_analysis.py.
This particular example is done with a CSPad image, but only a single section is available. For more typical CSPad module, see next section.

Panel

Edit pyana.cfg to include configuration for xppt_image_analysis, and comment out the delay_scan module:

Code Block

	none
	none


[pyana]

modules = pyana_examples.xppt_image_analysis
#modules = pyana_examples.xppt_delayscan

[pyana_examples.xppt_image_analysis]
source = XppGon-0|Cspad-0
region = 127.3, 188.4, 95.1, 126.9

[pyana_examples.xppt_delayscan]
controlpv = fs2:ramp_angsft_target
ipimb_norm = XppSb3Ipm-1|Ipimb-0
ipimb_sig = XppSb3Pim-1|Ipimb-0
threshold = 0.1
outputfile = point_scan_delay.npy

Then run the xppt_image_analysis pyana module on xppi0112 run 55:

Code Block

	none
	none


pyana -n 10 /reg/d/psdm/XPP/xppi0112/xtc/e162-r0055-s00-c00.xtc

*Edit pyana.cfg again and comment out the region parameter (add a semicolon ";" to the beginning of the line).
Run again a single event and try to select a region by mouse clicks instead:*

Code Block

	none
	none


pyana -n 1 /reg/d/psdm/XPP/xppi0112/xtc/e162-r0055-s00-c00.xtc

Hit the "Zoom to rectangle" button in the matplotlib toolbar.
Zoom in on a rectangle around the bright spot in the "Region of interest" plot to the right.
You should now see the region marked out in the left window.
Hit the "Zoom" button once more, to go back to normal mode.
Click on the red rectangle in the left plot to print the region parameters and new Center of mass to screen.

Here are some code snippet highlights from the xppt_image_analysis.py module:

For each event, fetch the CsPad information, and get the image array:

Code Block

	none
	none


   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.

Code Block

	none
	none


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

Code Block

	none
	none


        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:

Code Block

	none
	none


       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:

Code Block

	none
	none


        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

The python code for this pyana module example resides in XtcExplorer/src/pyana_image.py.

Panel

Try some plotting of CSPad data using xtcexplorer. Launch the explorer and load xpp48712 run 66 (a dark run):

Code Block

	none
	none


xtcexplorer /reg/d/psdm/XPP/xpp48712/xtc/e153-r0066-s00-c00.xtc

*Look through a couple of events, then "Quit Pyana" and edit the configuration file. Add an output file name, and switch to "NoDisplay" and run 100 events to collect an average of dark images.
With darks collected, load another file from the same experiment: run 141. Edit the pyana configuration file to use the file you just generated to subtract darks. Run the explorer in "SlideShow" mode again.
Change the color scale of the plot by left and right clicking on the colorbar.

CSPad data structure

CSPad data in xtc is a list of elements. In pyana get the list from the evt (event) object (notice the need for the env (environment) object too!):

Code Block

	none
	none


elements = evt.getCsPadQuads(self.source, env)

elements here is a python list of ElementV1 or ElementV2 (or later versions) objects, each representing one quadrant. The list is not ordered, so to know which quadrant you have, you have to check with element.quad(). To store a local array of the whole CSPad detector, you can do the following.

In beginjob, find out from the configuration object what part of the CSPad was in use (sometimes sections are missing):
Code Block
none
none
def beginjob ( self, evt, env ) : config = env.getConfig(xtc.TypeId.Type.Id_CspadConfig, self.source) if not config:

CSPad tile arrangement

To get a rough picture of the full detector, here's an example of how XtcExplorer/src/cspad.py does it:

For each Quadrant (cspad_layout.txt):

Code Block

none


    def get_quad_image( self, data3d, qn) :
        """get_quad_image
        Get an image for this quad (qn)

        @param data3d           3d data array (row vs. col vs. section)
        @param qn               quad number
        """    
        pairs = []
        for i in range (8) :
        
            # 1) insert gap between asics in the 2x1
            asics = np.hsplit( data3d[i], 2)print '*** cspad config object is missing ***'
            return

 gap = np.zeros( (185,3), dtype=data3d.dtype )
  quads = range(4)

        #
 memorize this list of sections for later
     # gap should beself.sections 3 pixels wide
  = map(config.sections, quads)

In each event, get the current CSPad data:

Code Block

	none
	none


def event(self, evt, env):
    elements = evt.getCsPadQuads(self.source,env)

    pairpixel_array = np.hstackzeros( (asics[0], gap, asics[1]) (4,8,185,388), dtype="uint16")

    for element in elements:
     # all sections aredata originally 185 (rows) x 388 (columns) = element.data() # the 3-dimensional data array (list of 2d sections)
        quad =  element.quad() # Re-orientcurrent eachquadrant sectionnumber in the quad(integer value)

        # if any sections if i==0 or i==1 :are missing, insert zeros
        if len( data ) < 8 :
  pair  = pair[:,::-1].T   # reverse columns, switch columnszsec to= rows. np.zeros( (185,388), dtype=data.dtype)
            iffor i==4 or i==5 in range (8) :
                if i pairnot =in pair[::-1,:].Tself.sections[quad] :
    # reverse rows, switch rows to columns
          data = pairsnp.appendinsert( pairdata, i, zsec, axis=0 )

        pixel_array[quad] = data

What we have so far gives you a 4d numpy array of all pixels. And if you want to store it in e.g. a numpy array, you can reshape it down to 2 dimensions (this is the format of the official pedestal files made by the translator):

Code Block

	none
	none


pixels = pixel_array.reshape(1480,388)
np.save("pixel_pedestal_file.npy", pixels )

CSPad tile arrangement

To get a rough picture of the full detector, here's an example of how XtcExplorer/src/cspad.py does it:

For each Quadrant (cspad_layout.txt):

Code Block

	none
	none


    def get_quad_image( self, data3d, qn) :   if self.small_angle_tilt :
                pair = scipy.ndimage.interpolation.rotate(pair,self.tilt_array[qn][i])

        # make the array for this quadrant
        quadrant = np.zeros( (self.npix_quad, self.npix_quad), dtype=data3d.dtype )

"""get_quad_image
        #Get insertan theimage 2x1for sectionsthis according toquad (qn)

        for@param secdata3d in range (8):
        3d data array (row nrows,vs. ncolscol = pairs[sec].shape

vs. section)
        @param qn    # colp,rowp are where the top-left corner of a section should bequad placednumber
        """
    rowp = self.npix_quad - self.sec_offset[0] - (self.section_centers[0][qn][sec] + nrows/2) pairs = []
        for i in  colp = self.npix_quad - self.sec_offset[1] - (self.section_centers[1][qn][sec] + ncols/2)
range (8) :

            # 1) insert gap between asics in the 2x1
            quadrant[rowp:rowp+nrows, colp:colp+ncols] = pairs[sec][0:nrows,0:ncols]
asics = np.hsplit( data3d[i], 2)
            gap if= np.zeros(rowp+nrows > self.npix_quad) or (colp+ncols > self.npix_quad) :
 (185,3), dtype=data3d.dtype )
            #
            print "ERROR"
     # gap should be 3 pixels wide
           print rowp, ":", rowp+nrows, ", ", colp, ":",colp+ncolspair = np.hstack( (asics[0], gap, asics[1]) )

        return quadrant

Then combine all four quadrant images into the full detector image:

...

# all sections are originally

...

185 (rows) x 388 (columns)
            # Re-orient

...

each

...

section in

...

the quad

...

if i=

...

=0 or i==1 :

...

  pair = pair[:,::-1].T   # reverse columns, switch columns

...

to

...

rows. if i==4 or

...

i==5 :

...

pair =

...

pair[::-1,:].T   # reverse rows, switch rows to columns

...

pairs.append( pair )

...

if self.small_angle_tilt :
                pair

...

= scipy.ndimage.interpolation.rotate(pair,self.tilt_array[qn][i])

        # make the array for this quadrant
        quadrant = np.zeros( (850, 850), dtype=data3d.dtype )

        # insert the 2x1 sections according to
        for sec in range (8):
            nrows, ncols = pairs[sec].shape

            # colp,rowp are where the top-left corner of a section should be placed
            rowp = 850 - self.sec_offset[0] - (self.section_centers[0][qn][sec] + nrows/2)
            colp = 850 - self.sec_offset[1] - (self.section_centers[1][qn][sec] + ncols/2)

            quadrant[rowp:rowp+nrows, colp:colp+ncols] = pairs[sec][0:nrows,0:ncols]

        return quadrant

Then combine all four quadrant images into the full detector image:

Code Block

	none
	none


        self.image = np.zeros((2*850+100, 2*850+100 ), dtype="float64")
        for quad in xrange (4):

Fine tuning

Notice that the code snippets above make use of some predefined quantities which it reads in from "calibration files". The files contains calibrated numerical values for individual sections' and quads' rotations and shifts. All of these files are located in the experiment's 'calib' folder, but is not generated automatically. The XtcExplorer currently has a local version which is not correct but which is close enough to display a reasonable image. For how the files have been read in, you can take a look at XtcExplorer/src/cspad.py's read_alignment function.

CSPad alignment

Saving data arrays

saving numpy arrays to numpy file

Code Block
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

Code Block
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/^{Image Removed}

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/^{Image Removed}

[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

quad_image

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

Panel

title	Loading your arrays into (I)Python and plotting interactively:

Code Block

Plotting with MatPlotLib

Matplotlib:

http://matplotlib.sourceforge.net/users/pyplot_tutorial.html^{Image Removed}

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.
Code Block
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:
Code Block
from matplotlib.pyplot import * ion() plot(array) draw()

Fine tuning

Notice that the code snippets above make use of some predefined quantities which it reads in from "calibration files". The files contains calibrated numerical values for individual sections' and quads' rotations and shifts. All of these files are located in the experiment's 'calib' folder, but is not generated automatically. The XtcExplorer currently has a local version which is not correct but which is close enough to display a reasonable image. For how the files have been read in, you can take a look at XtcExplorer/src/cspad.py's read_alignment function.

For how to find the correct constants for each experiment, look at the CSPAD Alignment page.

...

Non-interactive batch analysis

...

Page tree

Versions Compared

Old Version 38

New Version Current

Key

The Basics

Python

Pyana

xtcscanner

xtcexplorer

Extracting the data with pyana, some examples

Outline of a pyana module

NumPy, SciPy and MatPlotLib

Saving data arrays

Point detector delay scan

Image peak finding

Plotting with MatPlotLib

Interactive analysis with IPython

Extracting the data with pyana, some examples

Outline of a pyana module

CSPad images and tile arangements

CSPad data structure

Datatypes, and how to find data from your detector/device in the xtc file.

Point detector delay scan

Image peak finding

CSPad images and tile arangements

CSPad data structure

CSPad tile arrangement

CSPad tile arrangement

Fine tuning

Saving data arrays

saving numpy arrays to numpy file

saving to MatLab file

saving to HDF5 file

Interactive analysis with IPython

Plotting with MatPlotLib

Fine tuning

Non-interactive batch analysis