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

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:

...

    fee_energy_array = evt.getFeeGasDet()
    gdENRC11 = fee_energy_array[0]
    gdENRC12 = fee_energy_array[1]
    gdENRC21 = fee_energy_array[2]
    gdENRC22 = 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

...

     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

def event(self, evt, env):
    evrdata  try:
  = evt.getEvrData("NoDetector-0|Evr-0")
    
    for i encoderin =range evt(evrdata.get(xtc.TypeId.Type.Id_EncoderData, self.enc_source )numFifoEvents()):
        encoder_value = encoder.value()
    except:
        print "No encoder found in this event"
        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

print "Event code: ", evrdata.fifoEvent(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)

def event(self,evt,env):
    try:
        encoder = evt.get(xtc.TypeId.Type.Id_EncoderData, self.enc_source )
        encoder_value = encoder.value()
    except:
        print "No encoder found in this event"
        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!)

Time data

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

def event ( self, evt,# envEncoder )Parameters :
 to convert to picoseconds
   event delay_timea = 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.

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

def event ( self, evt, env ) :
    # raw data
    ipmRaw 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:

def event(self, evt, env):
    # raw data
    ipmRaw = evt.get(xtc.TypeId.Type.Id_IpimbData, source )
    try:get(xtc.TypeId.Type.Id_IpimbData, source )
    try:
        ch = [ipmRaw.channel0(),
              ipmRaw.channel1(),
              ipmRaw.channel2(),
              ipmRaw.channel3() ]
                
        ch_volt = [ipmRaw.channel0Voltschannel0(),
                   ipmRaw.channel1VoltsipmRaw.channel1(),
                   ipmRaw.channel2Voltschannel2(),
                   ipmRaw.channel3Voltschannel3() ]
    except:
         pass

   
 #  feature-extracted data
    ipmFexch_volt = evt[ipmRaw.get(xtc.TypeId.Type.Id_IpmFex, source )
channel0Volts(),
      try:
         # array of 4 numbers
 ipmRaw.channel1Volts(),
          fex_channel = ipmFex.channel 

         # scalar valuesipmRaw.channel2Volts(),
         fex_sum = ipmFex.sum 
       ipmRaw.channel3Volts()]
  fex_xpos = ipmFex.xposexcept:
         fex_ypos = ipmFex.ypospass

     except:# feature-extracted data
    ipmFex     pass

def event(self, evt, env)= evt.get(xtc.TypeId.Type.Id_IpmFex, source )
    try:
    ipm = evt.getSharedIpimbValue("HFX-DG3-IMB-02")
   # #array of or4 equivalently:numbers
      # ipm   fex_channel = evt.get(xtc.TypeId.Type.Id_SharedIpimb, "HFX-DG3-IMB-02")
    try: ipmFex.channel 

        ### Raw# datascalar ###values
        # arraysfex_sum of= 4ipmFex.sum numbers:
        ch fex_xpos = [ipm.ipimbData.channel0(),ipmFex.xpos
         fex_ypos = ipmFex.ypos

   ipm.ipimbData.channel1(),  except:
         pass

def event(self, evt, env):
    ipm.ipimbData.channel2(), = evt.getSharedIpimbValue("HFX-DG3-IMB-02")
    # or equivalently:
    # ipm =  ipm.ipimbData.channel3()]evt.get(xtc.TypeId.Type.Id_SharedIpimb, "HFX-DG3-IMB-02")
        ch_volt = [^cspad_simple.py]

try: 
        ### Feature-extractedRaw data ###
        # arrayarrays of 4 numbers:
        fex_channelsch = [ipm.ipmFexData.channel ipimbData.channel0(),
        
      ipm.ipimbData.channel1(),
  # scalars:
        fex_sum =  ipm.ipmFexData.sum 
ipimbData.channel2(),
            fex_xpos = ipm.ipmFexData.xposipimbData.channel3()]
        fexch_yposvolt = [ipm.ipmFexData.ypos

ipimbData.channel0Volts(),
          except:
         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:

ipm.ipimbData.channel1Volts(),
     def __init__ ( self ):
        # initialize data ipm.ipimbData.channel2Volts(),
        self.address =  "AmoITof-0|Acqiris-0"
        self.data = []ipm.ipimbData.channel3Volts()]

        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:

### Feature-extracted data ###
     def beginjob ( self,# evt,array envof )4 numbers:
        cfgfex_channels = env.getConfig( _pdsdata.xtc.TypeId.Type.Id_AcqConfig, self.address )ipm.ipmFexData.channel 
        self.num = cfg.nbrChannels()

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

    def event ( self, evt, env ) # scalars:
        channelfex_sum = ipm.ipmFexData.sum 0
        acqDatafex_xpos = evt.getAcqValue( self.address, channel, env)ipm.ipmFexData.xpos
        iffex_ypos acqData := ipm.ipmFexData.ypos

            self.counter+=1except:
            wf = acqData.waveform()   # returns a waveform array of numpy.ndarray type.
    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:

    def __init__ ( self ):
        # initialize data
        self.data.append(wf)

At the end of the job, take the average and plot it:

    def endjob( self, env ) :

address =  "AmoITof-0|Acqiris-0"
        self.data = np.array(self.data)[]
  # this is an array of shape (Nevents, nSamples)
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:

    def beginjob ( self, #evt, takeenv the mean of all events for each sampling time
        xs = np.mean(data, axis=0)

) :
        cfg = env.getConfig( _pdsdata.xtc.TypeId.Type.Id_AcqConfig, self.address )
        self.num = cfg.nbrChannels()

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

plt.plot(xs)

    def event (  plt.xlabel('Seconds')self, evt, env ) :
        channel = 0
        plt.ylabel('Volts'acqData = evt.getAcqValue( self.address, channel, env)
        plt.show()

Which gives you a plot like this
Image Removed

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:

   def __init__ ( self ):
        # initialize data
        self.address =  "SxrEndstation-0|Princeton-0"
if acqData :
            self.counter+=1
            wf = acqData.waveform()   # returns a waveform array of numpy.ndarray type.
            self.data = None
.append(wf)

At the end of the job, take the average and plot itIn each event, we add the image array returned from the getPrincetonValue function:

    def event endjob( self, evt, env ) :

       frame data = evtnp.getPrincetonValuearray( self.address, envdata)
  # this is an array if frame :of shape (Nevents, nSamples)

           # accumulatetake the data
mean of all events for each sampling time
    if self.data is None :
xs = np.mean(data, axis=0)

            self.data = np.float_(frame.data())plt.plot(xs)

           else :plt.xlabel('Seconds')
        plt.ylabel('Volts')
       self.data += frame.data plt.show()

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

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:

   def endjob( self, env ) __init__ ( self ):
        plt.imshow( self.data/self.countpass, origin='lower')# initialize data
        plt.colorbar()
        plt.show()

self.address =  "SxrEndstation-0|Princeton-0"
        #self.data save the full image to a png file= None

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

  def event ( self, evt, env ) :

       frame = pltevt.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

PnCCD image

def event(self,evt,env):
    trygetPrincetonValue( self.address, env)
       if frame :
        frame = evt.getPnCcdValue( self.source, env )  # accumulate the data
           if self.data is None :
        image       self.data = np.float_(frame.data())
    except       else :
              pass

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

These all use the generic getFrameValue function:

 self.data += frame.data()

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

   def endjob( self, env ) 

def event(self,evt,env):
    try:
        frame = evt.getFrameValueplt.imshow( self.source )data/self.countpass, origin='lower')
        plt.colorbar()
      image = frameplt.datashow()

    except:    # save the full image to a png file
        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.

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

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

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

These all use the generic getFrameValue function:

def event(self,evt,env):
   # converttry:
 to 16-bit integer
    image.dtype frame = np.uint16

CsPad data

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

def event(self,evt,env):evt.getFrameValue( self.source )
    quads    image = evtframe.getCsPadQuads(self.img_source, envdata()
    if not quads except:
        print '*** cspad information is missing ***'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.

def event(self,evt,env):
    try:
        return
frame = evt.getFrameValue(      
    # dump information about quadrants
self.source )
        printimage "Number of quadrants: %d" % len(quads= frame.data()
    except:
    for  q in quads:pass

    # convert to 16-bit printinteger
 "  Quadrant image.dtype = np.uint16

CsPad data

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

def event(self,evt,env):%d" % q.quad()
    quads = evt.getCsPadQuads(self.img_source, env)
 print  " if not quads virtual_channel:
  %s"   % q.virtual_channel()
  print '*** cspad information is missing print***'
 "    lane: %s" % q.lane()return
        print
 "   # tid:dump %s"information % q.tid()about quadrants
        print "Number of   acq_countquadrants: %s%d" % q.acq_countlen(quads)
    
    print "for q in quads:
        op_code: %sprint "  Quadrant %d" % q.op_codequad()
        print "    seqvirtual_countchannel: %s" % q.seqvirtual_countchannel()
        print "    tickslane: %s" % q.tickslane()
        print "    fiducialstid: %s" % q.fiducialstid()
        print "    frameacq_typecount: %s" % q.frameacq_typecount()
        print "    sbop_tempcode: %s" % map(q.sb_temp, range(4)op_code()
        print "   
 seq_count: %s" % q.seq_count()
     # image data as 3-dimentional array
   print "    ticks: %s" % q.ticks()
         data =print "    fiducials: %s" % q.data()
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 expectedSo 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 arrayrowcolsec. 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:

...

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 XtcExplorer/src/cspad.py or in cspad_simple.pyin 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:

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 

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