#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));
  }
}

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

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)

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]