Page History

...

This page holds a few example code-snippets for use in pyana analysis. The analysis is written in python and uses MatPlotLib.PyPlot for plotting of data. Compare with myana user examples to see how (some of) the same things can be done using the myana analysis framework.

BeamLine Data: EBeam

To read out energy, charge and position of the beam from the beamline data, use getEBeam(). It returns a class/structure that has the following members/fields:

The most reliable place for up-to-date information about all the event getters in pyana, see: https://confluence.slac.stanford.edu/display/PCDS/Pyana+Reference+Manual#PyanaReferenceManual-Classpyana.event.Event

For all the examples, you may assume that the pyana module contains a class with at least 'beginjob', 'event' and 'endjob' functions that starts something like this:

Code Block

	none
	none
title	outline of a pyana module


import numpy as np
import matplotlib.pyplot as plt
from pypdsdata import xtc

class mypyana(object):
    def __init__(self,source="")

Code Block

none


def event(self,evt,env):

    ebeam = evt.getEBeam()
    try :
        beamChrgself.source = ebeam.fEbeamChargesource
        beamEnrgself.counter = ebeam.fEbeamL3EnergyNone
        beamPosXself.array = ebeam.fEbeamLTUPosX
[]   # really just a list

 beamPosY   = ebeam.fEbeamLTUPosYdef beginjob(self,evt,env):
        beamAngXself.counter = ebeam.fEbeamLTUAngX 0

    def event(self,evt,env):
        beamAngYself.counter += ebeam.fEbeamLTUAngY1

        # beamPkCrsnippet code =goes ebeam.fEbeamPkCurrBC2here
        printthedata "ebeam: ", beamChrg, beamEnrg, beamPosX, beamPosY, beamAngX, beamAngY, beamPkCr= evt.get(xtc.TypeId.Type.Id_SomeType, self.source )
        self.array.append( thedata.somevalue )

    exceptdef endjob(self,evt,env):
        print "No EBeam object found"

BeamLine Data: Fee Gas Detector

To read out the energy from the front end enclosure (FEE) gas detector, use getFeeGasDet(). This returns and array of 4 numbers:

Code Block

none


    fee_energy_array = evt.getFeeGasDet()
    gdENRC11 = fee_energy_array[0]Job done! Processed %d events. " % self.counter

       # place for plotting etc

    gdENRC12 = fee_energy_array[1]
    gdENRC21 = fee_energy_array[2]
    gdENRC22 = fee_energy_array[3]

    energy = (gdENRC21 - gdENRC22) / 2
    # or use the first two that has a different gain:
    energy = (gdENRC11 - gdENRC12) / 2

BeamLine Data: Phase Cavity

 # convert from python list to a numpy array
       self.array = np.array( self.array )

       # plot graph
       plt.plot(self.array)

BeamLine Data: EBeam

To read out energy, charge and position of the beam from the beamline data, use getEBeam(). It returns a class/structure that has the following members/To read out fit time and charge of the phase cavity, use getPhaseCavity() which returns a structure with the following fields:

Code Block

	none
	none
title	getEBeam


def event(self,evt,env):


     pcebeam = evt.getPhaseCavitygetEBeam()
    try try:
        beamChrg pcFitTime1 = pcebeam.fFitTime1fEbeamCharge
        beamEnrg pcFitTime2 = pcebeam.fFitTime2fEbeamL3Energy
         pcCharge1beamPosX = pcebeam.fCharge1fEbeamLTUPosX
         pcCharge2beamPosY = pcebeam.fCharge2fEbeamLTUPosY
        beamAngX print "PhaseCavity: ", pcFitTime1,  pcFitTime2, pcCharge1, pcCharge2
      except := ebeam.fEbeamLTUAngX
        beamAngY = ebeam.fEbeamLTUAngY
        beamPkCr = ebeam.fEbeamPkCurrBC2
        print print"ebeam: "No, PhasebeamChrg, CavitybeamEnrg, object found"

Encoder data (delay scanner)

Code Block
none	none	def event(self,evt,env):beamPosX, beamPosY, beamAngX, beamAngY, beamPkCr tryexcept: encoderprint = evt.get(xtc.TypeId.Type.Id_EncoderData, self.enc_source )"No EBeam object found"

BeamLine Data: FEE Gas Detector

To read out the energy from the front end enclosure (FEE) gas detector, use getFeeGasDet(). This returns and array of 4 numbers:

Code Block

	none
	none
title	getFeeGasDet

        encoder_valuefee_energy_array = encoderevt.valuegetFeeGasDet()
    except:
gdENRC11 = fee_energy_array[0]
    gdENRC12  print "No encoder found in this event"= fee_energy_array[1]
    gdENRC21 = fee_energy_array[2]
    gdENRC22    return

You could combine it with phase cavity time, and compute a time delay from it, for example (I don't know the origin of these parameters!):

...


     # Encoder Parameters to convert to picoseconds
     delay_a = -80.0e-6;
     delay_b = 0.52168;
     delay_c = 299792458;
     delay_0 = 0;

    delay_time = (delay_a * encoder_value + delay_b)*1.e-3 / delay_c) 
    delay_time = 2 * delay_time / 1.0e-12 + delay_0 + pcFitTime1

Time data

The time of the event can be obtained within the event function:

...


def event ( self, evt, env ) :
    event_time = evt.getTime().seconds() + 1.0e-9*evt.getTime().nanoseconds() )

= fee_energy_array[3]

    energy = (gdENRC21 + gdENRC22) / 2.0
    # or use the first two that has a different gain:
    energy = (gdENRC11 + gdENRC12) / 2.0

BeamLine Data: Phase Cavity

To read out fit time and charge of the phase cavity, use getPhaseCavity() which returns a structure with the following fields:

Code Block

	none
	none
title	getPhaseCavity


     pc = evt.getPhaseCavity()
     try:
         pcFitTime1 = pc.fFitTime1
         pcFitTime2 = pc.fFitTime2
         pcCharge1 = pc.fCharge1
         pcCharge2 = pc.fCharge2
         print "PhaseCavity: ", pcFitTime1,  pcFitTime2, pcCharge1, pcCharge2
      except :
         print "No Phase Cavity object found"

Event code

Code Block

	none
	none
title	EnvData


def event(self, evt, env):
    evrdata = evt.getEvrData("NoDetector-0|Evr-0")
    
    for i in range (evrdata.numFifoEvents()):

IPIMB diode data

Currently there are two data structures that holds data from the same type of devices. Depending on DAQ
configuration, they are either DetInfo type or BldInfo type. Here are examples for extracting both types
in the user module event function:

Code Block

none


def event(self, evt, env):
    # raw data
    ipmRaw = evt.get(xtc.TypeId.Type.Id_IpimbData, source )
    try:
        ch = [ipmRaw.channel0(),
              ipmRaw.channel1(),
        print "Event code: ",   ipmRawevrdata.channel2fifoEvent(),
              ipmRaw.channel3() ]
                
        ch_volt = [ipmRaw.channel0Volts(),
                   ipmRaw.channel1Volts(),
                   ipmRaw.channel2Volts(),
    i).EventCode

In the example above, the address of the EvrData object is given as "NoDetector-0|Evr-0". The address may be different in other cases, so make sure you have the correct address. If you don't know what it is, you can use 'pyxtcreader -vv <xtcfile> | less' to browse your xtcfile and look for it. Look for a lines with 'contains=EvrConfig_V' or 'contains=EvrData_V'. The address will be found on the same line in 'src=DetInfo(<address>)'

Encoder data (delay scanner)

Code Block

	none
	none
title	EncoderData


def event(self,evt,env):
    try:
        encoder       ipmRaw.channel3Volts() ]= evt.get(xtc.TypeId.Type.Id_EncoderData, self.enc_source )
        encoder_value = encoder.value()
    except:
        pass

print "No encoder found #in feature-extracted datathis event"
    ipmFex = evt.get(xtc.TypeId.Type.Id_IpmFex, source ) return

You could combine it with phase cavity time, and compute a time delay from it, for example (I don't know the origin of these parameters!):

Code Block

	none
	none

    try:
# Encoder Parameters to convert to picoseconds
   # arraydelay_a of 4 numbers= -80.0e-6;
    delay_b = 0.52168;
    fexdelay_channelc = ipmFex.channel 
299792458;
    delay_0     # scalar values= 0;

    delay_time = (delay_a *  fexencoder_sumvalue = ipmFex.sum 
  + delay_b)*1.e-3 / delay_c) 
    delay_time = 2 fex* delay_xpostime =/ ipmFex.xpos
         fex_ypos = ipmFex.ypos

     except:
         pass

1.0e-12 + delay_0 + pcFitTime1

Time data

The time of the event can be obtained within the event function:

Code Block

	none
	none
title	getTime


def event ( self, evt, env ) :
    event_time = evt.getTime().seconds() + 1.0e-9*evt.getTime().nanoseconds() )

IPIMB diode data

This is data from sets of 4 diodes around the beam line (Intensity Position, Intensity Monitoring Boards)
that measures the beam intensity in four spots, from which we can also deduce the position of the beam.

Currently there are two data structures that holds data from the same type of devices. Depending on DAQ
configuration, they are either DetInfo type or BldInfo type. Here are examples for extracting both types
in the user module event function:

Code Block

	none
	none
title	DetInfo


def event(self, evt, env):
    # raw data
    ipmRaw = evt.get(xtc.TypeId.Type.Id_IpimbData, source )
    try:

Code Block

none


def event(self, evt, env):
    ipm = evt.getSharedIpimbValue("HFX-DG3-IMB-02")
    # or equivalently:
    # ipm = evt.get(xtc.TypeId.Type.Id_SharedIpimb, "HFX-DG3-IMB-02")
    try: 
        ### Raw data ###
        # arrays of 4 numbers:
        ch = [ipm.ipimbData.channel0(),
              ipm.ipimbData.channel1(),
            ch = ipm.ipimbData.channel2[ipmRaw.channel0(),
              ipmipmRaw.ipimbData.channel3channel1()],
            ch_volt = [ipm.ipimbData.channel0VoltsipmRaw.channel2(),
                   ipm.ipimbData.channel1Volts(),ipmRaw.channel3() ]
                
   ipm.ipimbData.channel2Volts     ch_volt = [ipmRaw.channel0Volts(),
                   ipmipmRaw.ipimbData.channel3Voltschannel1Volts()]
,
        ### Feature-extracted data ###
        # array of 4 numbers:
ipmRaw.channel2Volts(),
            fex_channels = ipm.ipmFexData.channel 
    ipmRaw.channel3Volts()]
    except:
        # scalars:
pass

    # feature-extracted data
    fex_sumipmFex = ipm.ipmFexData.sum evt.get(xtc.TypeId.Type.Id_IpmFex, source )
    try:
        fex_xpos = ipm.ipmFexData.xpos # array of 4 numbers
         fex_yposchannel = ipm.ipmFexData.yposipmFex.channel 

     except:
    # scalar values
   pass

Acqiris waveform data

This method can be used for any detector/device that has Acqiris waveform data. Edit the self.address field to get the detector of your choice.

Initialize class members:

Code Block
def fex__init__ ( self ): sum = ipmFex.sum # initialize data fex_xpos = ipmFex.xpos fex_ypos = self.address = "AmoITof-0\|Acqiris-0"ipmFex.ypos except: self.data = [] pass

Code Block

	none
	none
title	BldInfo


def event(self, evt, env):
    ipm self.counter = 0

If you're curious to see any of the Acqiris configuration, e.g. how many channels do we have, you can inspect the AcqConfig object:

Code Block

= evt.getSharedIpimbValue("HFX-DG3-IMB-02")
    def# beginjob ( self, evt, env ) or equivalently:
    # ipm   cfg = envevt.getConfigget( _pdsdata.xtc.TypeId.Type.Id_AcqConfig, self.address )SharedIpimb, "HFX-DG3-IMB-02")
    try: 
        self.num### = cfg.nbrChannels()

The read the event waveform data (an array) and append it to the self.data list:

Code Block

Raw data ###
    def event ( self, evt, env )  # arrays of 4 numbers:
        channelch = 0[ipm.ipimbData.channel0(),
        acqData  = evt.getAcqValue( self.address, channel, env) ipm.ipimbData.channel1(),
        if acqData :
    ipm.ipimbData.channel2(),
        self.counter+=1
      ipm.ipimbData.channel3()]
      wf  ch_volt = acqData.waveform()   # returns a waveform array of numpy.ndarray type.
[ipm.ipimbData.channel0Volts(),
                   selfipm.dataipimbData.appendchannel1Volts(wf)

...

Code Block

    def endjob( self, env ) :

        data = npipm.ipimbData.arraychannel2Volts(self.data)),
    # this is an array of shape (Nevents, nSamples)

        # take the mean of all events for each sampling timeipm.ipimbData.channel3Volts()]

        ### Feature-extracted data ###
        xs# array = np.mean(data, axis=0)

of 4 numbers:
        plt.plot(xs)
fex_channels = ipm.ipmFexData.channel 
        
        plt.xlabel('Seconds')
# scalars:
        fex_sum = plt.ylabel('Volts')
ipm.ipmFexData.sum 
        fex_xpos = plt.show()

Which gives you a plot like this
Image Removed

Display images from princeton camera

When plotting with MatPlotLib, we don't need to set the limits of the histogram manually, thus we don't need to read the Princeton configuration for this. If we want to sum the image from several events, we must first define and initialize some variables:

Code Block
def __init__ ( self ):ipm.ipmFexData.xpos fex_ypos = ipm.ipmFexData.ypos # initialize dataexcept: self.address = "SxrEndstation-0\|Princeton-0" self.data = None

In each event, we add the image array returned from the getPrincetonValue function:

 pass

Acqiris waveform data

This method can be used for any detector/device that has Acqiris waveform data. Edit the self.address field to get the detector of your choice.

Initialize class members:

Code Block

	none
	none


    def __init__ ( self ):

Code Block


  def event ( self, evt, env ) :

       frame = evt.getPrincetonValue( self.address, env)# initialize data
       if frame :
           # accumulate the data
    self.address =  "AmoITof-0|Acqiris-0"
       if self.data is None :
 = []
              self.datacounter = np.float_(frame.data())
  0

If you're curious to see any of the Acqiris configuration, e.g. how many channels do we have, you can inspect the AcqConfig object:

Code Block

	none
	none


    def beginjob ( self, evt, env else) :
        cfg = env.getConfig( _pdsdata.xtc.TypeId.Type.Id_AcqConfig, self.address )
  self.data += frame.data      self.num = cfg.nbrChannels()

At the end of the job, display/save the array:

The read the event waveform data (an array) and append it to the self.data list:

Code Block

	none
	none

Code Block


   def endjobevent ( self, evt, env ) :
        plt.imshow( self.data/self.countpass, origin='lower')
 channel = 0
       plt.colorbar( acqData = evt.getAcqValue( self.address, channel, env)
        plt.show()

if acqData :
          # save the full image to a png file
 self.counter+=1
            wf = pltacqData.imsave(fname="pyana_princ_image.png",arr=self.data, origin='lower')

Note that imsave saves the image only, pixel by pixel. If you want a view of the figure itself, lower resolution, you can save it from the interactive window you get from plt.show().
Image Removed

CsPad data

waveform()   # returns a waveform array of numpy.ndarray type.
            self.data.append(wf)

At the end of the job, take the average and plot itHere's an example of getting CsPad data from an event:

Code Block

	none
	none

    def eventendjob( self,evt, env ) :

    quads    data = evtnp.getCsPadQuadsarray(self.img_source, env)
data)  # this if not quads :is an array of shape (Nevents, nSamples)

        print '*** cspad information is missing ***'
    # take the mean of all events for each sampling time
    return
    xs =   np.mean(data, axis=0)

    # dump information about quadrants plt.plot(xs)

    print "Number of quadrants: %d" % len(quads plt.xlabel('Seconds')
    
    for q in quads:plt.ylabel('Volts')
        print "  Quadrant %d" % q.quad()
        print "    virtual_channel: %s" % q.virtual_channel()plt.show()

Which gives you a plot like this
Image Added

Princeton camera image

When plotting with MatPlotLib, we don't need to set the limits of the histogram manually, thus we don't need to read the Princeton configuration for this. If we want to sum the image from several events, we must first define and initialize some variables:

Code Block

	none
	none


   def __init__ ( self ):
        print# "initialize data
   lane: %s" % q.lane()
   self.address  =   print "    tid: %s" % q.tid()"SxrEndstation-0|Princeton-0"
        printself.data = None

In each event, we add the image array returned from the getPrincetonValue function:

Code Block

	none
	none
title	getPrincetonValue


  def event ( self, evt, env ) :

       frame = evt.getPrincetonValue( self.address, env"    acq_count: %s" % q.acq_count()
        print "    op_code: %s" % q.op_code()
       if print "    seq_count: %s" % q.seq_count()frame :
        print "  # accumulate ticks: %s" % q.ticks()
the data
          print "if self.data is None fiducials:
 %s" % q.fiducials()
        print "    frame_type: %s" % q.frame_type()
        print "    sb_temp: %s" % map(q.sb_temp, range(4self.data = np.float_(frame.data())
           else :
        # image data as 3-dimentional array
  self.data      data = q+= frame.data()

data2 will give you the third section stored, but be aware that sections sometimes are missing,
and in this case you'll need to check with the configuration information that you can obtain in beginjobAt the end of the job, display/save the array:

Code Block

none



   def beginjobendjob( self,evt, env ) :
    config =   envplt.getConfigimshow(xtc.TypeId.Type.Id_CspadConfig, self.img_source self.data/self.countpass, origin='lower')
    if not config:
  plt.colorbar()
      print '*** cspad config object is missing ***'
 plt.show()

        # save the full image return
to a png file
     
    quads = range(4)

    print 
    print "Cspad configuration"plt.imsave(fname="pyana_princ_image.png",arr=self.data, origin='lower')

Note that imsave saves the image only, pixel by pixel. If you want a view of the figure itself, lower resolution, you can save it from the interactive window you get from plt.show().
Image Added

PnCCD image

Code Block

	none
	none
title	getPnCcdValue


def event(self,evt,env):
    try:
    print "  N quadrantsframe = evt.getPnCcdValue( self.source, env )
        image = frame.data()
    except:
        pass

Other image (FCCD*,Opal,PIM (TM6740), ... )

These all use the generic getFrameValue function:

Code Block

	none
	none
title	getFrameValue


def event(self,evt,env):
    try:
        frame = evt.getFrameValue( self.source )
        image = frame.data()
    except:
        pass

FCCD (Fast CCD) image

The Fast CCD is read out as two 8-bit images, therefore you need this extra line to convert it to the right format.

Code Block

	none
	none
title	getFrameValue


def event(self,evt,env):
    try:
        frame = evt.getFrameValue( self.source )
        image = frame.data()
    except:
        pass

    # convert to 16-bit integer
    image.dtype = np.uint16

CsPad data

Here's an example of getting CsPad data from an event:

Code Block

	none
	none
title	getCsPadQuads


def event(self,evt,env):
    quads = evt.getCsPadQuads(self.img_source, env)
    if not quads :
        print '*** cspad information is missing ***'
        return
        
    # dump information about quadrants
    print "Number of quadrants: %d" % len(quads)
    
    for q in quads:
        print "  Quadrant %d" % q.quad()
        print "    virtual_channel: %s" % q.virtual_channel()
        print "    lane: %s" % q.lane()
        print "    tid: %s" % q.tid()
        print "    acq_count: %s" % q.acq_count()
        print "    op_code: %s" % q.op_code()
        print "    seq_count: %s" % q.seq_count()
        print "    ticks: %s" % q.ticks()
        print "    fiducials: %s" % q.fiducials()
        print "    frame_type: %s" % q.frame_type()
        print "    sb_temp: %s" % map(q.sb_temp, range(4))
            
        # image data as 3-dimentional array
        data = q.data()

So far so good. 'quads' is a list of CsPad Element objects, and not necessarily ordered in the expected way. So you'll need to use q.quad() to obtain the quad number.
q.data() gives you a 3D numpy array [row][col][sec]. Here sections will be ordered as expected, but be aware in case of missing sections, that you may need to check the
configuration object. You can get that from the env object, typically something you do in beginjob:

Code Block

	none
	none


def beginjob(self,evt,env):
    config = env.getConfig(xtc.TypeId.Type.Id_CspadConfig, self.img_source)
    if not config:
        print '*** cspad config object is missing ***'
        return        
    print "Cspad configuration"
    print "  N quadrants   : %d" % config.numQuads()
    print "  Quad mask     : %#x" % config.quadMask()
    print "  payloadSize   : %d" % config.payloadSize()
    print "  badAsicMask0  : %#x" % config.badAsicMask0()
    print "  badAsicMask1  : %#x" % config.badAsicMask1()
    print "  asicMask      : %#x" % config.asicMask()
    print "  numAsicsRead  : %d" % config.numAsicsRead()

   # get the indices of sections in use:
   qn = range(0,config.numQuads())               
   self.sections = map(config.sections, qn )

If you want to draw the whole CsPad image, there's currently no pyana function that does this. Pyana only supplies the pixels in a numpy array, and the
exact location of each pixel depends on the conditions at the time of data collection. A simplified way of making the image can be seen in cspad_simple.py(newer version (cspad.py) available if you check out the XtcExplorer package).

CSPad pixel coordinates.

The CSPad detector image can be drawn by positioning the sections from the data array into a large image array. This is done in cspad_simple.py above. The positions are extracted from optical meterology measurements and additional calibrations. Alternatively one can find the coordinate of each individual pixel from a pixel map, based on the same optical metrology measurements. This is described in details here

Epics Process Variables and ControlConfig

EPICS data is different from DAQ event data. It stores the conditions and settings of the instruments, but values typically change more slowly than your
average shot-by-shot data, and EPICS data is typically updated only when it changes, or every second, or similar. It is not stored in the 'evt' (event) object,
but in the 'env' (environment) object. You typically would read it only at the beginning of each job or if your doing a scan, you'd read it in every calibration cycle:

Code Block

	none
	none
title	env.epicsStore()


def begincalibcycle(self,evt,env):

    ## The returned value should be of the type epics.EpicsPvTime.
    pv = env.epicsStore().value( pv_name )
    if not pv:
        logging.warning('EPICS PV %s does not exist', pv_name)
    else:
        value = pv.value 
        status = pv.status 
        alarm_severity = pv.severity

Code Block

	none
	none
title	ControlConfig


def begincalibcycle(self,evt,env):
    ctrl_config = env.getConfig(xtc.TypeId.Type.Id_ControlConfig)
    
    nControls = ctrl_config.npvControls()
    for ic in range (0, nControls ):

        cpv = ctrl_config.pvControl(ic)
        name = cpv.name()
        value = cpv.value() : %d" % config.numQuads()
    print "  Quad mask     : %#x" % config.quadMask()
    print "  payloadSize   : %d" % config.payloadSize()
    print "  badAsicMask0  : %#x" % config.badAsicMask0()
    print "  badAsicMask1  : %#x" % config.badAsicMask1()
    print "  asicMask      : %#x" % config.asicMask()
    print "  numAsicsRead  : %d" % config.numAsicsRead()
    try:
        # older versions may not have all methods
        print "  roiMask       : [%s]" % ', '.join([hex(config.roiMask(q)) for q in quads])
        print "  numAsicsStored: %s" % str(map(config.numAsicsStored, quads))
    except:
        pass
    print "  sections      : %s" % str(map(config.sections, quads))
    print

Page tree

Versions Compared

Old Version 15

New Version Current

Key

BeamLine Data: EBeam

BeamLine Data: Fee Gas Detector

BeamLine Data: Phase Cavity

BeamLine Data: EBeam

Encoder data (delay scanner)

BeamLine Data: FEE Gas Detector

Time data

BeamLine Data: Phase Cavity

Event code

IPIMB diode data

Encoder data (delay scanner)

Time data

IPIMB diode data

Acqiris waveform data

Display images from princeton camera

Acqiris waveform data

CsPad data

Princeton camera image

PnCCD image

Other image (FCCD*,Opal,PIM (TM6740), ... )

FCCD (Fast CCD) image

CsPad data

CSPad pixel coordinates.

Epics Process Variables and ControlConfig