#include "UserAlg.h"
#include "RunStatisticsAlg.h"

#include "TF1.h"
#include "TROOT.h"
#include "TStyle.h"

#include <cmath>
#include <map>
#include <sstream>
#include <string>
#include <vector>

using std::cin;
using std::cout;
using std::endl;

using std::map;
using std::string;
using std::vector;

UserAlg::UserAlg(Event* event, Geometry* geo, Headers* runHeader):SystemAlg( event,geo,runHeader){
  m_event = event;
  m_geo = geo;
  m_runHeader = runHeader;

  m_calibration = new CalibrationAlg(m_event, m_geo, m_runHeader);
  m_canvasmanager = new CanvasManagerAlg(m_event, m_geo, m_runHeader);
  m_clustering = new ClusteringAlg(m_event, m_geo, m_runHeader);
  m_noise = new NoiseAlg(m_event, m_geo, m_runHeader);
  RunStatisticsAlg* runstatistics = new RunStatisticsAlg(m_event, m_geo, m_runHeader);
  vector<SystemAlg*> algVector;

  algVector.push_back(m_calibration);
  algVector.push_back(m_canvasmanager);
  algVector.push_back(m_clustering);
  algVector.push_back(m_noise);
  algVector.push_back(runstatistics);
  SystemAlg::setAlgs(algVector);
}

// initialize any variables here
bool UserAlg::initialize(){
  SystemAlg::initialize(); //do NOT remove this line

  // hec
  // enable clustering
  string hecNoiseFileName = "";
  m_hecClusteringIs2D = true;
  m_hecSortingAtEachIteration = false;
  m_hecCellThreshold = 2.;
  m_hecNeighborThreshold = 3.;
  m_hecSeedThreshold = 4.;
  m_clustering->hecEnable(hecNoiseFileName, m_hecClusteringIs2D, m_hecSortingAtEachIteration, m_hecCellThreshold, m_hecNeighborThreshold, m_hecSeedThreshold);
  // load calibration file
  m_calibration->hecLoad("/afs/cern.ch/user/h/hectbmon/public/tb/aug02/dig/coeff_hec_aug02_020903.dat");
  // read noise file
  m_noise->hecReadNoiseFile();

  // emec
  // enable clustering
  string emecNoiseFileName = "";
  m_emecClusteringIs2D = true;
  m_emecSortingAtEachIteration = false;
  m_emecCellThreshold = 2.;
  m_emecNeighborThreshold = 3.;
  m_emecSeedThreshold = 4.;
  m_clustering->emecEnable(emecNoiseFileName, m_emecClusteringIs2D, m_emecSortingAtEachIteration, m_emecCellThreshold, m_emecNeighborThreshold, m_emecSeedThreshold);
  // load calibration file
  m_calibration->emecLoad("/afs/cern.ch/user/h/hectbmon/public/tb/aug02/dig/coeff_emec_high_aug02_031203.dat");
  // read noise file
  m_noise->emecReadNoiseFile();

  // book histograms
  float Emin  = -10.;
  float Emax  = 240.;

  std::ostringstream os;
  os << m_runHeader->global()->getRunNo();

  // hec
  string hecHTitle = "run "  + os.str();
  if (m_hecClusteringIs2D) {
    hecHTitle += " 2D";
  } else {
    hecHTitle += " 3D";
  }
  if (m_hecSortingAtEachIteration) hecHTitle += " sorted";
  m_hecEH  = new TH1F("hecEH", hecHTitle.c_str(), 100, Emin, Emax);
  m_hecClustersEH  = new TH1F("hecClustersEH", hecHTitle.c_str(), 100, Emin, Emax);
  m_hecIterationsH = new TH1F("hecIterationsH", hecHTitle.c_str(), 20, -0.5, 19.5);
  m_hecNumberOfClustersH = new TH1F("hecNumberOfClustersH", hecHTitle.c_str(), 20, -0.5, 19.5);
  m_hecNumberOfCellsH = new TH1F("hecNumberOfCellsH", hecHTitle.c_str(), 50, -0.5, 49.5);
  m_hecNumberOfCellsUsingClustersH = new TH1F("hecNumberOfCellsUsingClustersH", hecHTitle.c_str(), 50, -0.5, 49.5);
  m_hecNumberOfCellsPerClusterH = new TH1F("hecNumberOfCellsPerClusterH", hecHTitle.c_str(), 50, -0.5, 49.5);
  for (int layer = 1; layer < 4; layer++) {
    std::ostringstream osl;
    osl << layer;
    string name = "hecNumberOfCellsPerClusterPerLayerHMap";
    name += osl.str();
    TH1F* h = new TH1F(name.c_str(), hecHTitle.c_str(), 25, -0.5, 24.5);
    m_hecNumberOfCellsPerClusterPerLayerHMap[layer] = h;
  }

  // emec
  string emecHTitle = "run "  + os.str();
  if (m_emecClusteringIs2D) {
    emecHTitle += " 2D";
  } else {
    emecHTitle += " 3D";
  }
  if (m_emecSortingAtEachIteration) emecHTitle += " sorted";
  m_emecEH  = new TH1F("emecEH", emecHTitle.c_str(), 100, Emin, Emax);
  m_emecClustersEH  = new TH1F("emecClustersEH", emecHTitle.c_str(), 100, Emin, Emax);
  m_emecIterationsH = new TH1F("emecIterationsH", emecHTitle.c_str(), 20, -0.5, 19.5);
  m_emecNumberOfClustersH = new TH1F("emecNumberOfClustersH", emecHTitle.c_str(), 50, -0.5, 49.5);
  m_emecNumberOfCellsPerClusterH = new TH1F("emecNumberOfCellsPerClusterH", emecHTitle.c_str(), 1000, -0.5, 999.5);
  m_emecNumberOfCellsH = new TH1F("emecNumberOfCellsH", emecHTitle.c_str(), 1000, -0.5, 999.5);
  m_emecNumberOfCellsUsingClustersH = new TH1F("emecNumberOfCellsUsingClustersH", emecHTitle.c_str(), 1000, -0.5, 999.5);
  for (int layer = 0; layer < 4; layer++) {
    std::ostringstream osl;
    osl << layer;
    string name = "emecNumberOfCellsPerClusterPerLayerHMap";
    name += osl.str();
    TH1F* h = new TH1F(name.c_str(), emecHTitle.c_str(), 500, -0.5, 499.5);
    m_emecNumberOfCellsPerClusterPerLayerHMap[layer] = h;
  }

  // combined
  string HTitle = "run "  + os.str();
  m_EH = new TH1F("EH", HTitle.c_str(), 100, Emin, Emax);
  m_ClustersEH = new TH1F("ClustersEH", HTitle.c_str(), 100, Emin, Emax);

  // find out what timing type
  if (m_runHeader->hec()->getTimingType() > 0) { // global cubic time
    cout << "Timing: global cubic time used" << endl;
  } else {
    cout << "Timing: TDC time used" << endl;
  }

  return true;
}

// put any code you want executed per event here
bool UserAlg::execute(){
  SystemAlg::execute(); //do NOT remove this line

  // only keep physics trigger events
  // Triggers (physics,electron,muon,pi,random)
  // Particle (fELECTRON=1, fMUON=2, fPION=3)
  if (m_event->getHecTriggerBit(0) == 0 || m_event->getHecTriggerBit(4) == 1) {  // not physics or a random
    return true;
  }

  // in the case of global cubic fit timing, only keep good events
  if (m_runHeader->hec()->getTimingType() > 0) { // global cubic time
    if (m_event->globalCubicTimeSigma() == 0) { // bad global timing
      return true;
    }
  }

  // choose an event number for debugging
  int evtN = 25;

  // hec
  // draw an event
  if (m_event->hecEventNumber() == evtN) {
    map<int, float> hecSignalOverNoiseMap;
    map<int, float> hecThresholdMap;
    map<int, float> hecClusterColorMap;
    map<int, int> hecClusterIDToColorMap;  // cluster ID keyed color index map
    // get the cluster ID map
    map<int, int> hecClusterIDMap = m_clustering->getHecClusterIDMap();
    // get the cluster map
    map<int, vector<int> > hecClusterMap = m_clustering->getHecClusterMap();
    // dump the clusters and get the cluster ID to color map for plotting
    int color = 0;
    for (map<int, vector<int> >::iterator it = hecClusterMap.begin(); it != hecClusterMap.end(); it++) {
      vector<int> cluster = it->second;
      cout << "hec  cluster " << it->first << ": ";
      for (vector<int>::iterator ic = cluster.begin(); ic != cluster.end(); ic++) {
        cout << *ic << " ";
      }
      cout << endl;
      hecClusterIDToColorMap[it->first] = color++;
    }
    // loop over connected cells
    for (int febno = m_runHeader->hec()->getFirstCh(); febno <= m_runHeader->hec()->getLastCh(); ++febno) {
      if (m_geo->hec_isConnected(febno)) {
        // get signal
        float s = m_event->hecCell(febno)->signal();
        if (m_runHeader->hec()->getEUnit() == 1) s = m_calibration->hecADCTonA(febno, s);  // ADC to nA
        // get noise
        float noise = m_noise->hecNoise(febno);
        // fill map of signal/noise for plotting
        if (noise > 0.) hecSignalOverNoiseMap[febno] = s/noise;
        // fill map of threshold category for plotting
        if (s > m_hecSeedThreshold*noise) {
          hecThresholdMap[febno] = 3.;
        } else if (fabs(s) > m_hecNeighborThreshold*noise) {
          hecThresholdMap[febno] = 2.;
        } else if (fabs(s) > m_hecCellThreshold*noise) {
          hecThresholdMap[febno] = 1.; 
        }
        // fill febno keyed cluster color map
        if (hecClusterIDMap[febno] >= 0) hecClusterColorMap[febno] = hecClusterIDToColorMap[hecClusterIDMap[febno]];
      }
    }
    m_geo->draw2DHec(hecSignalOverNoiseMap, 0, false, 0, 0, "S/N all cells");
    m_geo->draw2DHec(hecThresholdMap, 0, true, 0, 0, "thresholds");
    m_geo->draw2DHec(hecClusterColorMap, 0, true, 0, 0, "clusters");
  }

  float hecSignal = 0.;                  // total signal
  map<int, float> hecSignalPerLayerMap;  // layer keyed signal per layer map
  int hecNumberOfCells = 0;              // total number of cells used
  for (int febno = m_runHeader->hec()->getFirstCh(); febno <= m_runHeader->hec()->getLastCh(); ++febno) {
    if (m_geo->hec_isConnected(febno)) {
      // get signal
      float s = m_event->hecCell(febno)->signal();
      if (m_runHeader->hec()->getEUnit() == 1) s = m_calibration->hecADCTonA(febno, s);  // ADC to nA
      // get noise
      float noise = m_noise->hecNoise(febno);
      // only consider energy if above noise cut
      if (fabs(s) > m_hecCellThreshold*noise) {
        hecNumberOfCells++;
        // multiply by 2 for second wheel
        if (m_geo->hec_iz(febno) > 2) s *= 2.;
        // convert to GeV EM scale
        s *= 0.003266;  // 3.266 MeV/nA from NIM paper draft
        // cumulate in GeV EM scale 
        hecSignal += s;  
        hecSignalPerLayerMap[m_geo->hec_iz(febno)] += s;
      }
    }
  }
  m_hecEH->Fill(hecSignal);

  float hecClustersSignal = 0.;                  // total signal
  map<int, float> hecClustersSignalPerLayerMap;  // layer keyed clusters signal per layer map
  int hecNumberOfCellsUsingClusters = 0;         // total number of cells used
  m_hecIterationsH->Fill(m_clustering->getHecNumberOfIterations());
  map<int, vector<int> > hecClusterMap = m_clustering->getHecClusterMap();
  m_hecNumberOfClustersH->Fill((int)hecClusterMap.size());
  for (map<int, vector<int> >::iterator it = hecClusterMap.begin(); it != hecClusterMap.end(); it++) {
    vector<int> cluster = it->second;
    hecNumberOfCellsUsingClusters += (int)cluster.size();
    m_hecNumberOfCellsPerClusterH->Fill((int)cluster.size());
    map<int, int> hecNumberOfCellsPerClusterPerLayerMap;
    for (vector<int>::iterator ic = cluster.begin(); ic != cluster.end(); ic++) {
      int febno = *ic;
      hecNumberOfCellsPerClusterPerLayerMap[m_geo->hec_iz(febno)] += 1;
      // get signal
      float s = m_event->hecCell(febno)->signal();
      if (m_runHeader->hec()->getEUnit() == 1) s = m_calibration->hecADCTonA(febno, s);  // ADC to nA
      // get noise
      float noise = m_noise->hecNoise(febno);
      // multiply by 2 for second wheel
      if (m_geo->hec_iz(febno) > 2) s *= 2.;
      // convert to GeV EM scale
      s *= 0.003266;  // 3.266 MeV/nA from NIM paper draft
      // cumulate in GeV EM scale 
      hecClustersSignal += s;  
      hecClustersSignalPerLayerMap[m_geo->hec_iz(febno)] += s;
    }
    for (map<int, TH1F*>::iterator it = m_hecNumberOfCellsPerClusterPerLayerHMap.begin(); it != m_hecNumberOfCellsPerClusterPerLayerHMap.end(); it++) {
      // zero number of cells in a layer only make sense for 3d clusters
      if (m_hecClusteringIs2D) {
        if (hecNumberOfCellsPerClusterPerLayerMap[it->first] > 0) it->second->Fill(hecNumberOfCellsPerClusterPerLayerMap[it->first]);
      } else {
        it->second->Fill(hecNumberOfCellsPerClusterPerLayerMap[it->first]);
      } 
    }
  }
  m_hecClustersEH->Fill(hecClustersSignal);
  if (hecNumberOfCellsUsingClusters > 0) {
    m_hecNumberOfCellsUsingClustersH->Fill(hecNumberOfCellsUsingClusters);
    m_hecNumberOfCellsH->Fill(hecNumberOfCells);
  }

  // emec 
  // draw an event
  if (m_event->emecEventNumber() == evtN) {
    map<int, float> emecSignalOverNoiseMap;
    map<int, float> emecThresholdMap;
    map<int, float> emecClusterColorMap;
    map<int, int> emecClusterIDToColorMap;  // cluster ID keyed color index map
    // get the cluster ID map
    map<int, int> emecClusterIDMap = m_clustering->getEmecClusterIDMap();
    // get the cluster map
    map<int, vector<int> > emecClusterMap = m_clustering->getEmecClusterMap();
    // dump the clusters and get the cluster ID to color map for plotting
    int color = 0;
    for (map<int, vector<int> >::iterator it = emecClusterMap.begin(); it != emecClusterMap.end(); it++) {
      vector<int> cluster = it->second;
      cout << "emec cluster " << it->first << ": ";
      for (vector<int>::iterator ic = cluster.begin(); ic != cluster.end(); ic++) {
        cout << *ic << " ";
      }
      cout << endl;
      emecClusterIDToColorMap[it->first] = color++;
    }
    // loop over connected cells
    for (int febno = m_runHeader->emec()->getFirstCh(); febno <= m_runHeader->emec()->getLastCh(); ++febno) {
      if (m_geo->emec_isConnected(febno)) {
        // get signal
        float s = m_event->emecCell(febno)->signal();
        if (m_runHeader->emec()->getEUnit() == 1) s = m_calibration->emecADCTonA(febno, s);  // ADC to nA
        // get noise
        float noise = m_noise->emecNoise(febno);
        // fill map of signal/noise for plotting
        if (noise > 0.) emecSignalOverNoiseMap[febno] = s/noise;
        // fill map of threshold category for plotting
        if (s > m_emecSeedThreshold*noise) {
          emecThresholdMap[febno] = 3.;
        } else if (fabs(s) > m_emecNeighborThreshold*noise) {
          emecThresholdMap[febno] = 2.;
        } else if (fabs(s) > m_emecCellThreshold*noise) {
          emecThresholdMap[febno] = 1.; 
        }
        // fill febno keyed cluster color map
        if (emecClusterIDMap[febno] >= 0) emecClusterColorMap[febno] = emecClusterIDToColorMap[emecClusterIDMap[febno]];
      }
    }
    m_geo->draw2DEmec(emecSignalOverNoiseMap, 0, false, 0, 0, "S/N all cells");
    m_geo->draw2DEmec(emecThresholdMap, 0, true, 0, 0, "thresholds");
    m_geo->draw2DEmec(emecClusterColorMap, 0, true, 0, 0, "Clusters");
  }

  float emecSignal = 0.;                  // total signal
  map<int, float> emecSignalPerLayerMap;  // layer keyed signal per layer map
  int emecNumberOfCells = 0;              // total number of cells used
  for (int febno = m_runHeader->emec()->getFirstCh(); febno <= m_runHeader->emec()->getLastCh(); ++febno) {
    if (m_geo->emec_isConnected(febno)) {
      // get signal
      float s = m_event->emecCell(febno)->signal();
      if (m_runHeader->emec()->getEUnit() == 1) s = m_calibration->emecADCTonA(febno, s);  // ADC to nA
      // get noise
      float noise = m_noise->emecNoise(febno);
      // only consider energy if above noise cut
      if (fabs(s) > m_emecCellThreshold*noise) {
        emecNumberOfCells++;
        if (m_geo->emec_iz(febno) > 0) {  // don't use presampler for energy 
          // convert to GeV EM scale
          s *= 0.000430;  // 0.430 MeV/nA from NIM paper draft
          // cumulate in GeV EM scale 
          emecSignal += s;  
          emecSignalPerLayerMap[m_geo->emec_iz(febno)] += s;
        }
      }
    }
  }
  m_emecEH->Fill(emecSignal);

  float emecClustersSignal = 0.;                  // total signal
  map<int, float> emecClustersSignalPerLayerMap;  // layer keyed clusters signal per layer map
  int emecNumberOfCellsUsingClusters = 0;         // total number of cells used
  m_emecIterationsH->Fill(m_clustering->getEmecNumberOfIterations());
  map<int, vector<int> > emecClusterMap = m_clustering->getEmecClusterMap();
  m_emecNumberOfClustersH->Fill((int)emecClusterMap.size());
  for (map<int, vector<int> >::iterator it = emecClusterMap.begin(); it != emecClusterMap.end(); it++) {
    vector<int> cluster = it->second;
    emecNumberOfCellsUsingClusters += (int)cluster.size();
    m_emecNumberOfCellsPerClusterH->Fill((int)cluster.size());
    map<int, int> emecNumberOfCellsPerClusterPerLayerMap;
    for (vector<int>::iterator ic = cluster.begin(); ic != cluster.end(); ic++) {
      int febno = *ic;
      emecNumberOfCellsPerClusterPerLayerMap[m_geo->emec_iz(febno)] += 1;
      // get signal
      float s = m_event->emecCell(febno)->signal();
      if (m_runHeader->emec()->getEUnit() == 1) s = m_calibration->emecADCTonA(febno, s);  // ADC to nA
      // get noise
      float noise = m_noise->emecNoise(febno);
      if (m_geo->emec_iz(febno) > 0) {  // don't use presampler for energy
        // convert to GeV EM scale
        s *= 0.000430;  // 0.430 MeV/nA from NIM paper draft
        // cumulate in GeV EM scale 
        emecClustersSignal += s;  
        emecClustersSignalPerLayerMap[m_geo->emec_iz(febno)] += s;
      }
    }
    for (map<int, TH1F*>::iterator it = m_emecNumberOfCellsPerClusterPerLayerHMap.begin(); it != m_emecNumberOfCellsPerClusterPerLayerHMap.end(); it++) {
      // zero number of cells in a layer only make sense for 3d clusters
      if (m_emecClusteringIs2D) {
        if (emecNumberOfCellsPerClusterPerLayerMap[it->first] > 0) it->second->Fill(emecNumberOfCellsPerClusterPerLayerMap[it->first]);
      } else {
        it->second->Fill(emecNumberOfCellsPerClusterPerLayerMap[it->first]);
      } 
    }
  }
  m_emecClustersEH->Fill(emecClustersSignal);
  if (emecNumberOfCellsUsingClusters > 0) {
    m_emecNumberOfCellsUsingClustersH->Fill(emecNumberOfCellsUsingClusters);
    m_emecNumberOfCellsH->Fill(emecNumberOfCells);
  }
  // combined
  m_EH->Fill(hecSignal+emecSignal);
  m_ClustersEH->Fill(hecClustersSignal+emecClustersSignal);

  return true;
}

//put any code here that you want executed after all the events have been processed.
bool UserAlg::finalize(){
  SystemAlg::finalize();//do NOT remove this line

  // draw
  //  gROOT->SetStyle("Plain");
  gStyle->SetOptStat(111110);
  gStyle->SetOptFit(1111);
  gStyle->SetStatW(0.25); // statistics box (% of pad width?)

  int xCanvasSize = 800;
  int yCanvasSize = int(xCanvasSize*8.5/11.); // letter proportions
  {
    std::ostringstream os;
    os << "run_" << m_runHeader->global()->getRunNo() << ".ps";
    string psFile = os.str();
    m_canvasmanager->setPostScript(psFile.c_str());
  }

  m_canvasmanager->setCanvas("all cells", 2, 2, xCanvasSize, yCanvasSize);

  m_canvasmanager->getNextPad();  
  m_hecEH->GetXaxis()->SetTitle("hec energy (GeV EM scale)");
  m_hecEH->Draw(); m_canvasmanager->updateCanvas();
  m_canvasmanager->getNextPad();  
  m_emecEH->GetXaxis()->SetTitle("emec energy (GeV EM scale)");
  m_emecEH->Draw(); m_canvasmanager->updateCanvas();
  m_canvasmanager->getNextPad();  
  m_EH->GetXaxis()->SetTitle("hec+emec energy (GeV EM scale)");
  m_EH->Draw(); m_canvasmanager->updateCanvas();
  m_canvasmanager->getNextPad();  
  {
    // decide number of sigmas to fit
    int nSigma = 3;
    // histogram
    float EMean   = m_EH->GetMean();
    float ESigma  = m_EH->GetRMS();
    cout << "Histogram" << endl;
    cout << "  mean  = " << EMean << endl;
    cout << "  sigma = " << ESigma << endl;
    // first fit
    float EFitmin = EMean - nSigma*ESigma;
    if (EFitmin < m_EH->GetXaxis()->GetXmin()) EFitmin = m_EH->GetXaxis()->GetXmin();
    float EFitmax = EMean + nSigma*ESigma;
    if (EFitmax > m_EH->GetXaxis()->GetXmax()) EFitmax = m_EH->GetXaxis()->GetXmax();
    cout << "First fit" << endl;
    cout << "  fitting between " << EFitmin << " and " << EFitmax << endl;
    TF1* fGauss = new TF1("fGauss", "gaus", EFitmin, EFitmax);
    m_EH->Fit("fGauss", "QR0");
    EMean  = m_EH->GetFunction("fGauss")->GetParameter(1);
    ESigma = m_EH->GetFunction("fGauss")->GetParameter(2);
    float EMeanError  = m_EH->GetFunction("fGauss")->GetParError(1);
    float ESigmaError = m_EH->GetFunction("fGauss")->GetParError(2);
    cout << "  mean  = " << EMean << " +/-" << EMeanError << endl;
    cout << "  sigma = " << ESigma << " +/- " << ESigmaError << endl;
    // second fit
    EFitmin = EMean - nSigma*ESigma;
    if (EFitmin < m_EH->GetXaxis()->GetXmin()) EFitmin = m_EH->GetXaxis()->GetXmin();
    EFitmax = EMean + nSigma*ESigma;
    if (EFitmax > m_EH->GetXaxis()->GetXmax()) EFitmax = m_EH->GetXaxis()->GetXmax();
    cout << "Second fit" << endl;
    cout << "  fitting between " << EFitmin << " and " << EFitmax << endl;
    fGauss = new TF1("fGauss", "gaus", EFitmin, EFitmax);
    m_EH->Fit("fGauss", "QR0");
    EMean  = m_EH->GetFunction("fGauss")->GetParameter(1);
    ESigma = m_EH->GetFunction("fGauss")->GetParameter(2);
    EMeanError  = m_EH->GetFunction("fGauss")->GetParError(1);
    ESigmaError = m_EH->GetFunction("fGauss")->GetParError(2);
    cout << "  mean  = " << EMean << " +/-" << EMeanError << endl;
    cout << "  sigma = " << ESigma << " +/- " << ESigmaError << endl;
    // draw fit
    m_EH->Draw();
    fGauss->SetLineWidth(1);
    fGauss->Draw("LSAME"); m_canvasmanager->updateCanvas();
  }

  m_canvasmanager->setCanvas("cells in clusters", 2, 2, xCanvasSize, yCanvasSize);

  m_canvasmanager->getNextPad();  
  m_hecClustersEH->GetXaxis()->SetTitle("hec clusters energy (GeV EM scale)");
  m_hecClustersEH->Draw(); m_canvasmanager->updateCanvas();
  m_canvasmanager->getNextPad();  
  m_emecClustersEH->GetXaxis()->SetTitle("emec clusters energy (GeV EM scale)");
  m_emecClustersEH->Draw(); m_canvasmanager->updateCanvas();
  m_canvasmanager->getNextPad();  
  m_ClustersEH->GetXaxis()->SetTitle("hec+emec clusters energy (GeV EM scale)");
  m_ClustersEH->Draw(); m_canvasmanager->updateCanvas();
  m_canvasmanager->getNextPad();
  {  
    // decide number of sigmas to fit
    int nSigma = 3;
    // histogram
    float EMean   = m_ClustersEH->GetMean();
    float ESigma  = m_ClustersEH->GetRMS();
    cout << "Histogram" << endl;
    cout << "  mean  = " << EMean << endl;
    cout << "  sigma = " << ESigma << endl;
    // first fit
    float EFitmin = EMean - nSigma*ESigma;
    if (EFitmin < m_ClustersEH->GetXaxis()->GetXmin()) EFitmin = m_ClustersEH->GetXaxis()->GetXmin();
    float EFitmax = EMean + nSigma*ESigma;
    if (EFitmax > m_ClustersEH->GetXaxis()->GetXmax()) EFitmax = m_ClustersEH->GetXaxis()->GetXmax();
    cout << "First fit" << endl;
    cout << "  fitting between " << EFitmin << " and " << EFitmax << endl;
    TF1* fGauss = new TF1("fGauss", "gaus", EFitmin, EFitmax);
    m_ClustersEH->Fit("fGauss", "QR0");
    EMean  = m_ClustersEH->GetFunction("fGauss")->GetParameter(1);
    ESigma = m_ClustersEH->GetFunction("fGauss")->GetParameter(2);
    float EMeanError  = m_ClustersEH->GetFunction("fGauss")->GetParError(1);
    float ESigmaError = m_ClustersEH->GetFunction("fGauss")->GetParError(2);
    cout << "  mean  = " << EMean << " +/-" << EMeanError << endl;
    cout << "  sigma = " << ESigma << " +/- " << ESigmaError << endl;
    // second fit
    EFitmin = EMean - nSigma*ESigma;
    if (EFitmin < m_ClustersEH->GetXaxis()->GetXmin()) EFitmin = m_ClustersEH->GetXaxis()->GetXmin();
    EFitmax = EMean + nSigma*ESigma;
    if (EFitmax > m_ClustersEH->GetXaxis()->GetXmax()) EFitmax = m_ClustersEH->GetXaxis()->GetXmax();
    cout << "Second fit" << endl;
    cout << "  fitting between " << EFitmin << " and " << EFitmax << endl;
    fGauss = new TF1("fGauss", "gaus", EFitmin, EFitmax);
    m_ClustersEH->Fit("fGauss", "QR0");
    EMean  = m_ClustersEH->GetFunction("fGauss")->GetParameter(1);
    ESigma = m_ClustersEH->GetFunction("fGauss")->GetParameter(2);
    EMeanError  = m_ClustersEH->GetFunction("fGauss")->GetParError(1);
    ESigmaError = m_ClustersEH->GetFunction("fGauss")->GetParError(2);
    cout << "  mean  = " << EMean << " +/-" << EMeanError << endl;
    cout << "  sigma = " << ESigma << " +/- " << ESigmaError << endl;
    // draw fit
    m_ClustersEH->Draw();
    fGauss->SetLineWidth(1);
    fGauss->Draw("LSAME"); m_canvasmanager->updateCanvas();
  }

  m_canvasmanager->setCanvas("hec clusters", 2, 2, xCanvasSize, yCanvasSize);

  m_canvasmanager->getNextPad(); gPad->SetLogy();
  m_hecIterationsH->GetXaxis()->SetTitle("hec number of clustering iterations");
  m_hecIterationsH->Draw(); m_canvasmanager->updateCanvas();
  m_canvasmanager->getNextPad(); gPad->SetLogy(); 
  m_hecNumberOfClustersH->GetXaxis()->SetTitle("hec number of clusters");
  m_hecNumberOfClustersH->Draw(); m_canvasmanager->updateCanvas();
  m_canvasmanager->getNextPad(); gPad->SetLogy();
  m_hecNumberOfCellsPerClusterH->GetXaxis()->SetTitle("hec number of cells per cluster");
  m_hecNumberOfCellsPerClusterH->Draw(); m_canvasmanager->updateCanvas();
  m_canvasmanager->getNextPad(); gPad->SetLogy(); 
  m_hecNumberOfCellsUsingClustersH->GetXaxis()->SetTitle("hec number of cluster cells");
  m_hecNumberOfCellsUsingClustersH->Draw(); 
  m_hecNumberOfCellsH->SetLineStyle(3);  // dot-dot
  m_hecNumberOfCellsH->Draw("SAME"); m_canvasmanager->updateCanvas();
  for (map<int, TH1F*>::iterator it = m_hecNumberOfCellsPerClusterPerLayerHMap.begin(); it != m_hecNumberOfCellsPerClusterPerLayerHMap.end(); it++) {
    m_canvasmanager->getNextPad(); gPad->SetLogy(); 
    std::ostringstream osl;
    osl << it->first;
    string title = "hec number of cells per cluster, layer ";
    title += osl.str();
    it->second->GetXaxis()->SetTitle(title.c_str());
    it->second->Draw(); m_canvasmanager->updateCanvas();
  }

  m_canvasmanager->setCanvas("emec clusters", 2, 2, xCanvasSize, yCanvasSize);

  m_canvasmanager->getNextPad(); gPad->SetLogy();  
  m_emecIterationsH->GetXaxis()->SetTitle("emec number of clustering iterations");
  m_emecIterationsH->Draw(); m_canvasmanager->updateCanvas();
  m_canvasmanager->getNextPad(); gPad->SetLogy(); 
  m_emecNumberOfClustersH->GetXaxis()->SetTitle("emec number of clusters");
  m_emecNumberOfClustersH->Draw(); m_canvasmanager->updateCanvas();
  m_canvasmanager->getNextPad(); gPad->SetLogy(); 
  m_emecNumberOfCellsPerClusterH->GetXaxis()->SetTitle("emec number of cells per cluster");
  m_emecNumberOfCellsPerClusterH->Draw(); m_canvasmanager->updateCanvas();
  m_canvasmanager->getNextPad(); gPad->SetLogy(); 
  m_emecNumberOfCellsUsingClustersH->GetXaxis()->SetTitle("emec number of cluster cells");
  m_emecNumberOfCellsUsingClustersH->Draw(); 
  m_emecNumberOfCellsH->SetLineStyle(3);  // dot-dot
  m_emecNumberOfCellsH->Draw("SAME"); m_canvasmanager->updateCanvas();
  for (map<int, TH1F*>::iterator it = m_emecNumberOfCellsPerClusterPerLayerHMap.begin(); it != m_emecNumberOfCellsPerClusterPerLayerHMap.end(); it++) {
    m_canvasmanager->getNextPad(); gPad->SetLogy(); 
    std::ostringstream osl;
    osl << it->first;
    string title = "emec number of cells per cluster, layer ";
    title += osl.str();
    it->second->GetXaxis()->SetTitle(title.c_str());
    it->second->Draw(); m_canvasmanager->updateCanvas();
  }

  m_canvasmanager->closePostScript();

  return true;
}