Page History

import sys
import logging
import time

#-----------------------------
# Imports for other modules --
#-----------------------------
import numpy as np
import matplotlib.pyplot as plt
import h5py

#---------------------
#  Class definition --
#---------------------
class pyana_esort_mpl (object) :
    """Class whose instance will be used as a user analysis module. """

    #----------------
    #  Constructor --
    #----------------
    def __init__ ( self, nenergy, e1, e2, source = "AmoITof-0|Acqiris-0"):
histEnrg = hmgr.h1f("histEnrg",";Energy",self.nebins,self.emin,self.emax)
        self.histBin = hmgr.h1i("histBin",";Bin",self.nebins,0.0, self.nebins)
        self.histEnrg.SetMarkerStyle(4)
        self.histEnrg.SetMarkerColor(4)
        self.histBin.SetMarkerStyle(4)
        self.histBin.SetMarkerColor(4)

        name = "profWfAll"
        title = "Average Waveform;time [s];voltage [V]"
        self.prof_wfAll = hmgr.prof(name,title,self.numSamplesITof,
                            """Class constructor. The parameters to the constructor are passed0.0,
        from pyana configuration file. If parameters do not have default
        values  here then the must be defined in pyana.cfg. All parameters self.numSamplesITof*sampleIntervalITof,
        are passed as strings, convert to correct type before use.
        """
        self.source = source
 "")
        self.nebinsprof_wf = int(nenergy)
[]
        for bin self.eminin =range float(e10,self.nebins):
            self.emaxname = float(e2)"profWf%d"%bin
        self.bwidth = (self.emax-self.emin)/self.nebins

  title = "Average Waveform in Energy self.nev = 0

    #-------------------Bin %d;time [s];voltage [V]" % bin
    #  Public methods --
    #-------------------
self.prof_wf.append( hmgr.prof(name,title,self.numSamplesITof,
    def beginjob( self, evt, env ) :
        """This method is called once at the beginning of the job.
        It will also be called in any Xtc Configure transition0.0, if
        configurations have changed since last it was run.

        @param evt    event data object
        @param env    environment objectself.numSamplesITof*sampleIntervalITof,
        """

        # Preferred way to log information is via logging package
        logging.info( "pyana_esort_mpl.beginjob() called " )

      ""))


  cfg = env.getAcqConfig(self.source)
    # prof.FillN requires a #numChannels = cfg.nbrChannels()weight array...
        self.numSamplesITofwg = cfgnp.horizones().nbrSamples(self.numSamplesITof,dtype=float)
        sampleIntervalITofself.ts = cfg.horiz().sampInterval()None

    def beginrun( self, evt, #env for) matplotlib,:
 we store the array till the end instead of a histogram...
        """ energy / bin histogram """
     """This optional method is called if present at the beginning
        of the new run.

     self.axisEnrg = np.arange(self.emin,self.emax,self.bwidth, dtype=float)
  @param evt    event data object
 self.arrayEnrg = np.zeros(self.nebins, dtype=float)
    @param env   """ time wave"""environment object
        self.ts = None

"""
        logging.info( """ voltage wave"""
pyana_esort_hmgr.beginrun() called ")

        print  self.wf = None"beginrun is done, starting timer..."
        self.wf2start = Nonetime.clock()

    def begincalibcycle( self,  """ voltage wave in bins of energyevt, env ) :
        """
This optional method is called if present at self.bwf = []the beginning
        self.bwf2of =the []
new calibration cycle.

      for bin in@param rangeevt (0, self.nebins):
  event data object
        @param env   self.bwf.append( None )
 environment object
        """
      self.bwf2.append( None )  logging.info( "pyana_esort_hmgr.begincalibcycle() called ")




    def beginrunevent( self, evt, env ) :
        """This optional method is called iffor presentevery at the beginning
        of the new runL1Accept transition.

        @param evt    event data object
        @param env    environment object
        """
        logging.info( "pyana_esort_mplhmgr.beginrunevent() called ")

        print "beginrun is done, starting timer..."self.nev+=1

        self.startshotEnergy = timeevt.clockgetFeeGasDet()


     def begincalibcycle( self, evt, env ) :
 energy = 0.0
        bin  """This optional method is called if present at the beginning
   = 0
        if shotEnergy is None :
     of the new calibration cycle.

   print "No Gas Detector"
  @param evt    event dataelse object:
        @param env   energy environment object= (shotEnergy[2]+shotEnergy[3])/2
        """
        logging.info( "pyana_esort_mpl.begincalibcycle() called ")
self.histEnrg.Fill(energy);

    def event( self, evt, envacqData ) :
  = evt.getAcqValue(self.source, 0, env)
      """This method isif called for every L1Accept transition.

acqData :
         @param evt  wf  event data object
= acqData.waveform()

          @param env if self.ts is environmentNone object:
        """
        logging.info( "pyana_esort_mpl.event() called ")

self.ts = acqData.timestamps()

            self.nev+=1

prof_wfAll.FillN(self.numSamplesITof,self.ts,wf,self.wg,1)

         shotEnergy  = evt.getFeeGasDet( bin = int( ((energy-self.emin)/(self.emax-self.emin))*self.nebins+0.5)
        energy = 0.0
  if bin < 1  : bin = 0
            if shotEnergybin is> Noneself.nebins :
 bin = self.nebins-1
         print "No Gas Detector" self.histBin.Fill(bin)
        else :
   self.prof_wf[bin].FillN(self.numSamplesITof,self.ts,wf,self.wg,1)


    def endcalibcycle( self, env  energy = (shotEnergy[2]+shotEnergy[3])/2
) :
        acqData = evt.getAcqValue(self.source, 0, env)
        if acqData :"""This optional method is called if present at the end of the
            awf = acqData.waveform()calibration cycle.

        @param env   if self.ts is None :environment object
        """

        selflogging.ts = acqData.timestamps()info( "pyana_esort_hmgr.endcalibcycle() called" )

    def endrun( self, env ) :
   # average waveform
   """This optional method is called if present at the ifend self.wfof is None:the run.

        @param env       self.wf = awfenvironment object
        """
        selflogging.wf2 = (awf*awf)
  info( "pyana_esort_hmgr.endrun() called" )

    def endjob( self, env   else) :
        """This method is called at the end of self.wf += awfthe job. It should do
        final cleanup, e.g. close all open files.

   self.wf2 += (awf*awf)

   @param env    environment  object
  # binned waveform
    """
        bin = int( ((energy-self.emin)/(self.emax-self.emin))*self.nebins+0.5)logging.info( "pyana_esort_hmgr.endjob() called" )

        duration = time.clock() - self.start
   if bin < 1 : binprint ="duration 0
in seconds : ", duration
        ifprint bin"Number >of self.nebinsevents processed: bin", = self.nebins-1nev

          print  self.arrayEnrg[bin] += 1
self.prof_wfAll.GetNbinsX()
            if self.bwf[bin] is None for bin in range (20,40):
               print self.prof_wfAll.GetBinContent(bin+1), " +- ", self.prof_wfAll.GetBinError(bin+1)

//-----------------------
// This Class's Header --
//-----------------------
#include "speed/psana_esort.h"

//-----------------
// C/C++ Headers --
//-----------------
#include <time.h>

//-------------------------------
// Collaborating Class Headers --
//-------------------------------
#include "MsgLogger/MsgLogger.h"

//--------------------------------------------------
// to work with detector data include corresponding
// header from psddl_psana package
//--------------------------------------------------
#include "psddl_psana/acqiris.ddl.h"
#include "psddl_psana/bld.ddl.h"

//-----------------------------------------------------------------------
// Local Macros, Typedefs, Structures, Unions and Forward Declarations --
//-----------------------------------------------------------------------

// This declares this class as psana module
using namespace speed;
PSANA_MODULE_FACTORY(psana_esort)

//              ----------------------------------------
//              -- Public Function Member Definitions --
// @param env    environment object
        """

        logging.info( "pyana_esort_mpl.endcalibcycle() called" )

    def endrun( self, env ) :
        """This optional method is called if present at the end of the run.

        @param env    environment object
        """
        logging.info( "pyana_esort_mpl.endrun() called" )

    def endjob( self, env ) :
        """This method is called at the end of the job. It should do
        final cleanup, e.g. close all  ----------------------------------------

namespace speed {

//----------------
// Constructors --
//----------------
psana_esort::psana_esort (const std::string& name)
  : Module(name)
  , m_acqSrc()
  , m_feeSrc()
  , m_maxEvents()
  , m_filter()
  , m_nenergy()
  , m_e1()
  , m_e2()
{
  // get the values from configuration or use defaults
  m_acqSrc = configStr("acqSource", "DetInfo(:Acqiris)");
  m_feeSrc = configStr("feeSource", "BldInfo(FEEGasDetEnergy)");
  m_maxEvents = config("events", 32U);
  m_filter = config("filter", false);
  m_nenergy = config("nenergy",20);
  m_e1 = config("e1",0.0);
  m_e2 = config("e2",4.0);
}

//--------------
// Destructor --
//--------------
psana_esort::~psana_esort ()
{
}

/// Method which is called once at the beginning of the job
void
psana_esort::beginJob(Env& env)
{
  m_count = 0;

  RootHistoManager::RootHMgr& rhmgr = env.rhmgr();

  shared_ptr<Psana::Acqiris::ConfigV1> acqConfig = env.configStore().get(m_acqSrc);
  if (acqConfig.get()) {
    int channel = 0;
    m_nSamples       = acqConfig->horiz().nbrSamples();
    m_sampleInterval = acqConfig->horiz().sampInterval();
    m_slopeopen files.

        @param env    environment object
        """
        logging.info( "pyana_esort_mpl.endjob() called" )

        # have collected the sum so far, but we
        # want to display the average

        self.wf = self.wf / self.nev
        self.wf2 = self.wf2 / self.nev
        self.wf_rms = np.sqrt( self.wf2 - self.wf*self.wf ) / np.sqrt(self.nev)

        self.bwf_rms = []
        for bin in range (0, self.nebins):
            if self.bwf[bin] is None:
                self.bwf_rms.append(None)
                continue

            self.bwf[bin] = self.bwf[bin] / self.arrayEnrg[bin]
            self.bwf2[bin] = self.bwf2[bin] / self.arrayEnrg[bin]
            = acqConfig->vert(channel).slope();
self.bwf_rms.append(np.sqrt( self.bwf2[bin] \
      m_offset         = acqConfig->vert(channel).offset();

    m_timevector = new double[m_nSamples];
    m_wfvector = new double[m_nSamples];
    m_wgvector = new double[m_nSamples];
    for (int i(0);i<m_nSamples;i++){- self.bwf[bin]*self.bwf[bin] ) \
      m_timevector[i]= i * m_sampleInterval;
      m_wfvector[i]=0.0;
       m_wgvector[i]=1.0;
    }

    m_profWfAll = rhmgr.profile("profWfAll","Average Waveform;time [s];voltage [V]", / np.sqrt(self.arrayEnrg[bin]))

        duration = time.clock() - self.start
        print "duration in seconds : ", duration
      AxisDef(0.0, m_nSamples*m_sampleInterval,m_nSamples)
      print "nevents processed: ", self.nev
        duration = time.clock() - self.start
        print "duration in seconds : ",  );duration

    char name[32];
   print char title[232];
    for (int i=0; i<m_nenergy;i++){
"nevents processed: ", self.nev

        # save sprintf(name,"profWf%dhistograms",i);
 to an hdf5 file
  sprintf(title,"Average Waveform in Energy bin %d; time [s]; voltage [V]",i);
      m_profWf[i] = rhmgr.profile(name, title,ofile = h5py.File("outputfile.hdf5", 'w') # open for writing (overwrites existing file)
        group = ofile.create_group("HistogramArrays")
        group.create_dataset("wf_time",data=self.ts)
        group.create_dataset("wf_avg",data=self.wf)
        AxisDef(0.0,m_nSamples*m_sampleInterval,m_nSamplesgroup.create_dataset("wf_avg_rms",data=self.wf_rms)
          group.create_dataset("bin_enrg",data=self.arrayEnrg)
        group.create_dataset("bin_axis",data=self.axisEnrg)
        for bin in range     );
    }
  }
  // create a histogram to display the energy bins
  m_histEnrg = rhmgr.h1f("histEnrg","Energy bins;Energy; Event Count",AxisDef(m_e1,m_e2,m_nenergy));
  m_histBin = rhmgr.h1i("histBin","Energy bins;Bin; Event count",AxisDef(0.0,m_nenergy,m_nenergy));
  m_histEnrg->SetMarkerStyle(4); m_histEnrg->SetMarkerColor(4); // blue circles
  m_histBin->SetMarkerStyle(4);  m_histBin->SetMarkerColor(4); // blue circles

  //m_histEnrg = new TH1F("histEnrg",";Energy",m_nenergy,m_e1,m_e2);
  //m_histBin = new TH1I("histBin",";Bin",m_nenergy,0.0,m_nenergy);

}


/// Method which is called at the beginning of the run
void
psana_esort::beginRun(Env& env)
{
  printf("beginRun is done, starting timer...");
  m_start_time = clock();
}

/// Method which is called at the beginning of the calibration cycle
void
psana_esort::beginCalibCycle(Env& env)
{
}

/// Method which is called with event data, this is the only required
/// method, all other methods are optional
void
psana_esort::event(Event& evt, Env& env)
{

//   // this is how to skip event (all downstream modules will not be called)
//   if (m_filter && m_count % 10 == 0) skip();

//   // this is how to gracefully stop analysis job
//   if (m_count >= m_maxEvents) stop();

  //MsgLog(name(),info,"Event:" << m_count);

  /*
   * Processing FeeGasDet Value
   */
  shared_ptr<Psana::Bld::BldDataFEEGasDetEnergy> fee = evt.get(m_feeSrc);
  double shotEnergy[4] = {0.,0.,0.,0.};
  if (fee.get()) {
    shotEnergy[0] = fee->f_11_ENRC();
    shotEnergy[1] = fee->f_12_ENRC();
    shotEnergy[2] = fee->f_21_ENRC();
    shotEnergy[3] = fee->f_22_ENRC();
  }
  double energy = (shotEnergy[2]+shotEnergy[3])/2;
  if ( energy < m_e1 ) {
    m_histEnrg->Fill(m_e1);
  } else if ( energy > m_e2 ) {
    m_histEnrg->Fill(m_e2);
  } else {
    m_histEnrg->Fill(energy);
  }

  /*
   * Processing ITof waveform
   */
  shared_ptr<Psana::Acqiris::DataDescV1> acqData = evt.get(m_acqSrc);
  if (acqData.get()) {

    m_count++;

    int channel = 0;
    const Psana::Acqiris::DataDescV1Elem& data = acqData->data(channel);

    // copy myana main.cc getAcqValue()
    const int16_t* wfdata = data.waveforms();
    wfdata += data.indexFirstPoint();

    shared_ptr<Psana::Acqiris::ConfigV1> acqConfig = env.configStore().get(m_acqSrc);
    if (acqConfig.get()) {

//       const Psana::Acqiris::HorizV1& hcfg = acqConfig->horiz();
//       const Psana::Acqiris::VertV1& vcfg = acqConfig->vert(channel);
//       float slope = vcfg.slope();
//       float offset = vcfg.offset();
//       int nbrSamples = hcfg.nbrSamples();
//       for (int j=0;j<nbrSamples;j++){
//      int16_t swap = (wfdata[j]&0xff<<8) | (wfdata[j]&0xff00>>8);
//      m_wfvector[j] = swap*slope-offset;

      for (int j=0;j<m_nSamples;j++){
        int16_t swap = (wfdata[j]&0xff<<8) | (wfdata[j]&0xff00>>8);
        m_wfvector[j] = swap* m_slope - m_offset;
      }


      int numSamples = data.nbrSegments() * data.nbrSamplesInSeg();
      if (numSamples != m_nSamples) std::cout << "nSamples mismatch" << std::endl;

      int bin = int(( (energy-m_e1)/(m_e2-m_e1) )*m_nenergy+0.5);
      if (bin<1) bin=0;
      if (bin>m_nenergy-1) bin=int(m_nenergy-1);


      m_histBin->Fill(bin);
      m_profWf[bin]->FillN(numSamples,m_timevector,m_wfvector,m_wgvector,1);
      m_profWfAll->FillN(numSamples,m_timevector,m_wfvector,m_wgvector,1);
    }
  }
}

/// Method which is called at the end of the calibration cycle
void
psana_esort::endCalibCycle(Env& env)
{
  delete[] m_timevector;
}

/// Method which is called at the end of the run
void
psana_esort::endRun(Env& env)
{
}

/// Method which is called once at the end of the job
void
psana_esort::endJob(Env& env)
{

  clock_t finish_time = clock();
  double duration = (finish_time - m_start_time)/CLOCKS_PER_SEC;
  printf("Time elapsed: %.2f s \n",(duration));
  printf("Number of events processed: %d \n", m_count);

  std::cout << m_profWfAll->GetNbinsX() << std::endl;
  for (int bin = 20; bin<40; bin++) {
    std::cout << m_profWfAll->GetBinContent(bin+1) << " +- " << m_profWfAll->GetBinError(bin+1) << std::endl;
  }
}

} // namespace speed(0, self.nebins):
            if self.bwf[bin] is not None:
                group.create_dataset("wf_bin%d"%bin,data=self.bwf[bin])
                group.create_dataset("wf_bin%d_rms"%bin,data=self.bwf_rms[bin])
        ofile.close()


        duration = time.clock() - self.start
        print "duration till after save hdf5, in seconds : ", duration
        # plot
        fig = plt.figure(num=100, figsize=(7.5,5))
        ax = fig.add_subplot(111)
        plt.plot(self.axisEnrg, self.arrayEnrg, 'bo')
        plt.errorbar(self.axisEnrg, self.arrayEnrg, xerr=self.bwidth, yerr=np.sqrt(self.arrayEnrg),fmt=None)
        plt.title("Energy Bins")
        plt.xlabel("Energy")
        plt.ylabel("Event count")
        plt.draw()
        print "energy histogram plotted"
        fig.savefig("pyana_mpl_histEnrg.png")
        print "energy histogram saved"


        fig = plt.figure(num=200, figsize=(7.5,5))
        ax = fig.add_subplot(111)
        plt.clf()
        plt.errorbar(self.ts, self.wf,yerr=self.wf_rms, mew=0.0)
        plt.title("Average Waveform")
        plt.xlabel("time   [s]")
        plt.ylabel("voltage   [V]")
        plt.draw()
        print "wf avg histogram plotted"
        fig.savefig("pyana_mpl_profWfAll.png")
        print "wf avg histogram saved"

        for bin in range (0, self.nebins):
            if self.bwf[bin] is not None:
                # plot
                fig = plt.figure(num=300+bin, figsize=(7.5,5))
                ax = fig.add_subplot(111)
                plt.clf()
                plt.errorbar(self.ts, self.bwf[bin], yerr=self.bwf_rms[bin],mew=0.0)
                plt.title("Waveform bin %d (%.2f<E<%.2f)"\
                          % (bin, self.emin+bin*self.bwidth, self.emin+(bin+1)*self.bwidth))
                plt.xlabel("time    [s]")
                plt.ylabel("voltage    [V]")
                plt.draw()
                print "wf bin %d histogram plotted"%bin
                fig.savefig("pyana_mpl_profWf%d.png"%bin)
                print "wf bin %d histogram saved"%bin


        duration = time.clock() - self.start
        print "duration3 in seconds : ", duration

        # test
        print self.wf.size
        for bin in range (20,40):
            print self.wf[bin], " +-- ", self.wf_rms[bin]