Page History

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

...

The Basics

Python

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

Pyana

Analysis Workbook. Python-based Analysis

...

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

• http://numpy.scipy.org/^{Image Removed} - Numerical package for python (arrays etc)
• http://www.scipy.org/^{Image Removed} - Scientific tools
• http://matplotlib.sourceforge.net/^{Image Removed} - plotting package

Other useful links:

• http://www.scipy.org/NumPy_for_Matlab_Users^{Image Removed}
• http://www.scipy.org/^{Image Removed}
• http://www.sagemath.org/^{Image Removed}
• http://code.google.com/p/spyderlib/^{Image Removed}

Saving data arrays

...

import numpy as np

myarray = np.arange(0,100)
np.save( "output_file_name.npy", myarray) 
np.savetxt( "output_file_name.txt", myarray) 

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

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

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^{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()

...

However, the latest IPython has loads of new and interesting features...

http://ipython.org/^{Image Removed}

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

...

addpkg pyana_examples HEAD
scons

...

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 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 (parameter name ipimb_norm) and another IPIMB device (a.k.a. PIM3) channel 1 is used as the signal (parameter name ipimb_sig).

[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

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

...

Image peak finding

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

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.

[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

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

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

...

• 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

...

• 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: print '*** cspad config object is missing ***' return quads = range(4) # memorize this list of sections for later self.sections = 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) pixel_array = np.zeros((4,8,185,388), dtype="uint16") for element in elements: data = element.data() # the 3-dimensional data array (list of 2d sections) quad = element.quad() # current quadrant number (integer value) # if any sections are missing, insert zeros if len( data ) < 8 : zsec = np.zeros( (185,388), dtype=data.dtype) for i in range (8) : if i not in self.sections[quad] : data = np.insert( data, i, zsec, axis=0 ) pixel_array[quad] = data

...

• For each Quadrant (cspad_layout.txt):

  Code Block
  none 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) gap = np.zeros( (185,3), dtype=data3d.dtype ) # # gap should be 3 pixels wide pair = np.hstack( (asics[0], gap, asics[1]) ) # 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

...