You are viewing an old version of this page. View the current version.

Compare with Current View Page History

« Previous Version 8 Next »

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"
    #include <TROOT.h>
    #include <TH1F.h>
    #include <TProfile.h>
    #include <time.h>
    
    #include "myana.hh"
    #include "main.hh"
    
    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;
    
    const char* ParameterFileName = "ParameterFileName";
    
    FILE *LaserParamOut;
    FILE *ShotsOut;
    char ShotsName[200];
    
    using namespace Pds;
    
    // This function is called once at the beginning of the analysis job,
    // You can ask for detector "configuration" information here.
    
    void beginjob() {
    
      // get ITof parameters from myana.inp
      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;
    
    
      /*
       * Processing ITof Config Data
       */
      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
    }
    
    // This function is called once for each run.  You should check to see
    // if detector configuration information has changed.
    
    void beginrun()
    {
      printf("beginrun is done, starting timer...");
      start_time = clock();
    }
    
    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
       */
      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
       */
      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]
    
    
    • No labels