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
Open a terminal at pslogin or psana, and type:
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'.
Try:
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).
Try:
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.
Try:
xtcexplorer /reg/d/psdm/xpp/xppi0310/xtc/e81-r0098-s0*
{none}
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.
h2. Extracting the data with pyana, some examples
h3. 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:
- 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
h3. 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:
\\
{code: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:
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 theevtobject, and
eventmember 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
CSPad images and tile arangements
Saving data arrays
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...
Plotting with 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.
import matplotlib.pyplot as plt plt.ion() plt.plot(array) plt.draw()
Non-interactive batch analysis
Overview
Content Tools