Analysis code examples (comparisons of frameworks)
Here's one comparison of an example analysis in three different frameworks:
- myana (online C++ framework), using ROOT histogramming package
- psana (offline C++ framework), using ROOT histogramming package
- pyana (offline python framework)
- using ROOT histogramming package
- alternative without ROOT, using NumPy to store arrays and MatPlotLib to plot those arrays/
They all make histograms/plots of AMO ITof wavefunctions, in 20 energy bins (energy measured by FEE Gas Detector).
Timing and plots from all frameworks (pdf)
myana_esort.cc
// header files for this class
#include "myana.hh"
#include "main.hh"
// these are needed to make ROOT histograms
#include <TROOT.h>
#include <TH1F.h>
#include <TProfile.h>
// for tracking execution time
#include <time.h>
// global variables:
static int numChannelsITof;
static int numSamplesITof;
static double sampleIntervalITof;
static const int maxITofBin = 50;
static TProfile* profWf[maxITofBin];
static TProfile* profWfAll;
static TH1F* histEnrg;
static TH1I* histBin;
static unsigned int shotCountITof = 0;
float e1,e2;
int nenergy;
static int bin;
clock_t start_time;
clock_t finish_time;
using namespace Pds;
// This method is called once per job
void beginjob() {
// get ITof parameters from myana.inp (configuration file)
const char inpname[20]="myana.inp";
FILE *inpfile;
inpfile=fopen(inpname,"r");
fscanf(inpfile,"%d\n",&nenergy); // nbins of histogram
fscanf(inpfile,"%f\n",&e1); // xmin of histogram
fscanf(inpfile,"%f\n",&e2); // xmax of histogram
printf("%d %e %e\n",nenergy,e1,e2);
fclose(inpfile);
printf("calling beginjob()\n"); fflush(stdout);
int i;
int fail = 0;
/*
* Get ITof Config Data at the beginning of the job
*/
fail = getAcqConfig( AmoITof, numChannelsITof, numSamplesITof, sampleIntervalITof);
// create an "averaging" histogram (called "profile")
profWfAll = new TProfile("profWfAll","Average Waveform;Seconds;Volts"
, numSamplesITof
, 0.0,sampleIntervalITof,"");
char name[32];
char title[232];
for (i=0;i<nenergy;i=i+1) {
sprintf(name,"profWf%d",i);
sprintf(title,"Average Waveform in Energy bin %d",i);
profWf[i] = new TProfile(name,title,numSamplesITof,
0.0,sampleIntervalITof,"");
profWf[i]->SetYTitle("Volts"); //optional
profWf[i]->SetXTitle("Seconds"); //optional
}
// create a histogram to display the energy bins
histEnrg = new TH1F("histEnrg","Energy bins;Energy; Event Count",nenergy,e1,e2);
histBin = new TH1I("histBin","Energy bins;Bin; Event Count",nenergy,0.0,nenergy);
histEnrg->SetMarkerStyle(4); histEnrg->SetMarkerColor(4); // blue circles
histBin->SetMarkerStyle(4); histBin->SetMarkerColor(4); // blue circles
printf("beginjob is done, starting timer...");
start_time = clock();
}
// This method is called once per Run
void beginrun()
{
}
// This method is called once per CalibCycle
void begincalib()
{
}
// This is called once every shot. You can ask for
// the individual detector shot data here.
void event()
{
int fail;
/*
* Processing FeeGasDet Value for this event
*/
double shotEnergy[4];
fail = getFeeGasDet( shotEnergy );
double energy = (shotEnergy[2]+shotEnergy[3])/2;
if ( energy < e1 ) histEnrg->Fill(e1);
else if ( energy > e2 ) histEnrg->Fill(e2);
else histEnrg->Fill(energy);
/*
* Processing ITof waveform for this event
*/
double* timeITof;
double* voltageITof;
int channel = 0;
fail = getAcqValue( AmoITof, channel, timeITof, voltageITof);
if ( fail != 0 )
printf( "event(): getAcqValue() failed, code = %d\n", fail );
else
{
bin = int(( (energy-e1)/(e2-e1) )*nenergy+0.5);
if (bin<1) bin=0;
if (bin>nenergy-1) bin=int(nenergy-1);
for (int i=0;i<numSamplesITof;i=i+1)
{
double t = timeITof[i];
double v = voltageITof[i];
profWf[bin]->Fill(t,v);
profWfAll->Fill(t,v);
}
histBin->Fill(bin);
shotCountITof++;
}
}
void endcalib()
{
}
void endrun()
{
printf("User analysis endrun() routine called.\n");
}
void endjob()
{
printf("User analysis endjob() routine called.\n");
finish_time = clock();
double duration = (finish_time - start_time)/CLOCKS_PER_SEC;
printf("Time elapsed: %.2f s \n",(duration));
printf("Number of events processed: %d \n", shotCountITof);
printf("%d \n",profWfAll->GetNbinsX());
for (int bin=20;bin<40;bin++){
printf("%.10f +- %.10f \n",profWfAll->GetBinContent(bin+1),profWfAll->GetBinError(bin+1));
}
}
"psana_esort.cpp"
//-----------------------
// 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 --
// ----------------------------------------
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_slope = acqConfig->vert(channel).slope();
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++){
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]",
AxisDef(0.0, m_nSamples*m_sampleInterval,m_nSamples)
);
char name[32];
char title[232];
for (int i=0; i<m_nenergy;i++){
sprintf(name,"profWf%d",i);
sprintf(title,"Average Waveform in Energy bin %d; time [s]; voltage [V]",i);
m_profWf[i] = rhmgr.profile(name, title,
AxisDef(0.0,m_nSamples*m_sampleInterval,m_nSamples)
);
}
}
// 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
}
/// 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)
{
//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()) {
for (int j=0;j<m_nSamples;j++){
m_wfvector[j] = wfdata[j] * 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
"pyana_esort_hmgr.py"
import sys
import logging
import time
#-----------------------------
# Imports for other modules --
#-----------------------------
import numpy as np
import matplotlib.pyplot as plt
#---------------------
# Class definition --
#---------------------
class pyana_esort_hmgr (object) :
"""Class whose instance will be used as a user analysis module. """
#----------------
# Constructor --
#----------------
def __init__ ( self, nenergy, e1, e2, source = "AmoITof-0|Acqiris-0"):
"""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.
"""
self.source = source
self.nebins = int(nenergy)
self.emin = float(e1)
self.emax = float(e2)
self.nev = 0
#-------------------
# Public methods --
#-------------------
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 transition, if
configurations have changed since last it was run.
@param evt event data object
@param env environment object
"""
# Preferred way to log information is via logging package
logging.info( "pyana_esort_hmgr.beginjob() called " )
cfg = env.getAcqConfig(self.source)
#numChannels = cfg.nbrChannels()
self.numSamplesITof = cfg.horiz().nbrSamples()
sampleIntervalITof = cfg.horiz().sampInterval()
print self.numSamplesITof
print sampleIntervalITof
hmgr = env.hmgr()
self.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,
0.0,
self.numSamplesITof*sampleIntervalITof,
"")
self.prof_wf = []
for bin in range (0,self.nebins):
name = "profWf%d"%bin
title = "Average Waveform in Energy Bin %d;time [s];voltage [V]" % bin
self.prof_wf.append( hmgr.prof(name,title,self.numSamplesITof,
0.0,
self.numSamplesITof*sampleIntervalITof,
""))
# prof.FillN requires a weight array...
self.wg = np.ones(self.numSamplesITof,dtype=float)
self.ts = None
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
"""
logging.info( "pyana_esort_hmgr.beginrun() called ")
print "beginrun is done, starting timer..."
self.start = time.clock()
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
"""
logging.info( "pyana_esort_hmgr.begincalibcycle() called ")
def event( self, evt, env ) :
"""This method is called for every L1Accept transition.
@param evt event data object
@param env environment object
"""
logging.info( "pyana_esort_hmgr.event() called ")
self.nev+=1
shotEnergy = evt.getFeeGasDet()
energy = 0.0
bin = 0
if shotEnergy is None :
print "No Gas Detector"
else :
energy = (shotEnergy[2]+shotEnergy[3])/2
self.histEnrg.Fill(energy);
acqData = evt.getAcqValue(self.source, 0, env)
if acqData :
wf = acqData.waveform()
if self.ts is None :
self.ts = acqData.timestamps()
self.prof_wfAll.FillN(self.numSamplesITof,self.ts,wf,self.wg,1)
bin = int( ((energy-self.emin)/(self.emax-self.emin))*self.nebins+0.5)
if bin < 1 : bin = 0
if bin > self.nebins : bin = self.nebins-1
self.histBin.Fill(bin)
self.prof_wf[bin].FillN(self.numSamplesITof,self.ts,wf,self.wg,1)
def endcalibcycle( self, env ) :
"""This optional method is called if present at the end of the
calibration cycle.
@param env environment object
"""
logging.info( "pyana_esort_hmgr.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_hmgr.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 open files.
@param env environment object
"""
logging.info( "pyana_esort_hmgr.endjob() called" )
duration = time.clock() - self.start
print "duration in seconds : ", duration
print "Number of events processed: ", self.nev
print self.prof_wfAll.GetNbinsX()
for bin in range (20,40):
print self.prof_wfAll.GetBinContent(bin+1), " +- ", self.prof_wfAll.GetBinError(bin+1)
"pyana_esort_mpl.py"
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"):
"""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.
"""
self.source = source
self.nebins = int(nenergy)
self.emin = float(e1)
self.emax = float(e2)
self.bwidth = (self.emax-self.emin)/self.nebins
self.nev = 0
#-------------------
# Public methods --
#-------------------
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 transition, if
configurations have changed since last it was run.
@param evt event data object
@param env environment object
"""
# Preferred way to log information is via logging package
logging.info( "pyana_esort_mpl.beginjob() called " )
cfg = env.getAcqConfig(self.source)
#numChannels = cfg.nbrChannels()
self.numSamplesITof = cfg.horiz().nbrSamples()
sampleIntervalITof = cfg.horiz().sampInterval()
# for matplotlib, we store the array till the end instead of a histogram...
""" energy / bin histogram """
self.axisEnrg = np.arange(self.emin,self.emax,self.bwidth, dtype=float)
self.arrayEnrg = np.zeros(self.nebins, dtype=float)
""" time wave"""
self.ts = None
""" voltage wave"""
self.wf = None
self.wf2 = None
""" voltage wave in bins of energy """
self.bwf = []
self.bwf2 = []
for bin in range (0, self.nebins):
self.bwf.append( None )
self.bwf2.append( None )
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
"""
logging.info( "pyana_esort_mpl.beginrun() called ")
print "beginrun is done, starting timer..."
self.start = time.clock()
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
"""
logging.info( "pyana_esort_mpl.begincalibcycle() called ")
def event( self, evt, env ) :
"""This method is called for every L1Accept transition.
@param evt event data object
@param env environment object
"""
logging.info( "pyana_esort_mpl.event() called ")
self.nev+=1
shotEnergy = evt.getFeeGasDet()
energy = 0.0
bin = 0
if shotEnergy is None :
print "No Gas Detector"
else :
energy = (shotEnergy[2]+shotEnergy[3])/2
acqData = evt.getAcqValue(self.source, 0, env)
if acqData :
awf = acqData.waveform()
if self.ts is None :
self.ts = acqData.timestamps()
# average waveform
if self.wf is None:
self.wf = awf
self.wf2 = (awf*awf)
else :
self.wf += awf
self.wf2 += (awf*awf)
# binned waveform
bin = int( ((energy-self.emin)/(self.emax-self.emin))*self.nebins+0.5)
if bin < 1 : bin = 0
if bin > self.nebins : bin = self.nebins-1
self.arrayEnrg[bin] += 1
if self.bwf[bin] is None :
self.bwf[bin]=awf
self.bwf2[bin]=(awf*awf)
else :
self.bwf[bin]+=awf
self.bwf2[bin]+=(awf*awf)
def endcalibcycle( self, env ) :
"""This optional method is called if present at the end of the
calibration cycle.
@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 open 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]
self.bwf_rms.append(np.sqrt( self.bwf2[bin] \
- self.bwf[bin]*self.bwf[bin] ) \
/ np.sqrt(self.arrayEnrg[bin]))
duration = time.clock() - self.start
print "duration in seconds : ", duration
print "nevents processed: ", self.nev
duration = time.clock() - self.start
print "duration in seconds : ", duration
print "nevents processed: ", self.nev
# save "histograms" to an hdf5 file
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)
group.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 (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]
Overview
Content Tools