EcalPerEvtLaserAnalyzer.cc

Go to the documentation of this file.

/* 
 *  \class EcalPerEvtLaserAnalyzer
 *
 *  primary author: Julie Malcles - CEA/Saclay
 *  author: Gautier Hamel De Monchenault - CEA/Saclay
 */
 
#include <TFile.h>
#include <TTree.h>
 
#include "EcalPerEvtLaserAnalyzer.h"
 
#include <sstream>
#include <iostream>
#include <iomanip>
#include <ctime>
 
#include <FWCore/MessageLogger/interface/MessageLogger.h>
#include <FWCore/Utilities/interface/Exception.h>
 
#include <CalibCalorimetry/EcalLaserAnalyzer/interface/MEEEGeom.h>
#include <CalibCalorimetry/EcalLaserAnalyzer/interface/MEEBGeom.h>
 
#include <FWCore/Framework/interface/Event.h>
#include <FWCore/Framework/interface/MakerMacros.h>
#include <FWCore/ParameterSet/interface/ParameterSet.h>
#include <FWCore/Framework/interface/EventSetup.h>
#include <FWCore/Framework/interface/ESHandle.h>
 
#include <Geometry/EcalMapping/interface/EcalElectronicsMapping.h>
#include <Geometry/EcalMapping/interface/EcalMappingRcd.h>
 
#include <DataFormats/EcalDigi/interface/EcalDigiCollections.h>
#include <DataFormats/EcalDetId/interface/EcalElectronicsId.h>
#include <DataFormats/EcalDetId/interface/EcalDetIdCollections.h>
#include <DataFormats/EcalRawData/interface/EcalRawDataCollections.h>
 
#include <vector>
 
#include <CalibCalorimetry/EcalLaserAnalyzer/interface/TPNFit.h>
#include <CalibCalorimetry/EcalLaserAnalyzer/interface/PulseFitWithFunction.h>
 
using namespace std;
 
//========================================================================
EcalPerEvtLaserAnalyzer::EcalPerEvtLaserAnalyzer(const edm::ParameterSet& iConfig)
    //========================================================================
    : iEvent(0),
 
      // framework parameters with default values
      _nsamples(iConfig.getUntrackedParameter<unsigned int>("nSamples", 10)),
      _presample(iConfig.getUntrackedParameter<unsigned int>("nPresamples", 3)),
      _firstsample(iConfig.getUntrackedParameter<unsigned int>("firstSample", 1)),
      _lastsample(iConfig.getUntrackedParameter<unsigned int>("lastSample", 2)),
      _samplemin(iConfig.getUntrackedParameter<unsigned int>("sampleMin", 3)),
      _samplemax(iConfig.getUntrackedParameter<unsigned int>("sampleMax", 9)),
      _nsamplesPN(iConfig.getUntrackedParameter<unsigned int>("nSamplesPN", 50)),
      _presamplePN(iConfig.getUntrackedParameter<unsigned int>("nPresamplesPN", 6)),
      _firstsamplePN(iConfig.getUntrackedParameter<unsigned int>("firstSamplePN", 7)),
      _lastsamplePN(iConfig.getUntrackedParameter<unsigned int>("lastSamplePN", 8)),
      _timingcutlow(iConfig.getUntrackedParameter<unsigned int>("timingCutLow", 3)),
      _timingcuthigh(iConfig.getUntrackedParameter<unsigned int>("timingCutHigh", 7)),
      _niter(iConfig.getUntrackedParameter<unsigned int>("nIter", 3)),
      _fedid(iConfig.getUntrackedParameter<unsigned int>("fedID", 0)),
      _tower(iConfig.getUntrackedParameter<unsigned int>("tower", 1)),
      _channel(iConfig.getUntrackedParameter<unsigned int>("channel", 1)),
      _ecalPart(iConfig.getUntrackedParameter<std::string>("ecalPart", "EB")),
      nCrys(NCRYSEB),
      nTT(NTTEB),
      nSides(NSIDES),
      IsFileCreated(0),
      nCh(0),
      runType(-1),
      runNum(0),
      dccID(-1),
      lightside(2),
      doesRefFileExist(0),
      ttMat(-1),
      peakMat(-1),
      peak(-1),
      evtMat(-1),
      colMat(-1)
//========================================================================
 
{
  //now do what ever initialization is needed
 
  resdir_ = iConfig.getUntrackedParameter<std::string>("resDir");
  refalphabeta_ = iConfig.getUntrackedParameter<std::string>("refAlphaBeta");
 
  digiCollection_ = iConfig.getParameter<std::string>("digiCollection");
  digiPNCollection_ = iConfig.getParameter<std::string>("digiPNCollection");
  digiProducer_ = iConfig.getParameter<std::string>("digiProducer");
 
  eventHeaderCollection_ = iConfig.getParameter<std::string>("eventHeaderCollection");
  eventHeaderProducer_ = iConfig.getParameter<std::string>("eventHeaderProducer");
 
  // Define geometrical constants
 
  if (_ecalPart == "EB") {
    nCrys = NCRYSEB;
    nTT = NTTEB;
  } else {
    nCrys = NCRYSEE;
    nTT = NTTEE;
  }
}
 
//========================================================================
EcalPerEvtLaserAnalyzer::~EcalPerEvtLaserAnalyzer() {
  //========================================================================
 
  // do anything here that needs to be done at desctruction time
  // (e.g. close files, deallocate resources etc.)
}
 
//========================================================================
void EcalPerEvtLaserAnalyzer::beginJob() {
  //========================================================================
 
  // Define temporary files' names
 
  stringstream namefile1;
  namefile1 << resdir_ << "/ADCSamples.root";
 
  ADCfile = namefile1.str();
 
  // Create temporary file and trees to save adc samples
 
  ADCFile = new TFile(ADCfile.c_str(), "RECREATE");
 
  stringstream name;
  name << "ADCTree";
 
  ADCtrees = new TTree(name.str().c_str(), name.str().c_str());
 
  ADCtrees->Branch("iphi", &phi, "phi/I");
  ADCtrees->Branch("ieta", &eta, "eta/I");
  ADCtrees->Branch("dccID", &dccID, "dccID/I");
  ADCtrees->Branch("event", &event, "event/I");
  ADCtrees->Branch("color", &color, "color/I");
  ADCtrees->Branch("adc", &adc, "adc[10]/D");
  ADCtrees->Branch("ttrigMatacq", &ttrig, "ttrig/D");
  ADCtrees->Branch("peakMatacq", &peak, "peak/D");
  ADCtrees->Branch("pn0", &pn0, "pn0/D");
  ADCtrees->Branch("pn1", &pn1, "pn1/D");
 
  ADCtrees->SetBranchAddress("ieta", &eta);
  ADCtrees->SetBranchAddress("iphi", &phi);
  ADCtrees->SetBranchAddress("dccID", &dccID);
  ADCtrees->SetBranchAddress("event", &event);
  ADCtrees->SetBranchAddress("color", &color);
  ADCtrees->SetBranchAddress("adc", adc);
  ADCtrees->SetBranchAddress("ttrigMatacq", &ttrig);
  ADCtrees->SetBranchAddress("peakMatacq", &peak);
  ADCtrees->SetBranchAddress("pn0", &pn0);
  ADCtrees->SetBranchAddress("pn1", &pn1);
 
  IsFileCreated = 0;
  nCh = nCrys;
}
 
//========================================================================
void EcalPerEvtLaserAnalyzer::analyze(const edm::Event& e, const edm::EventSetup& c) {
  //========================================================================
 
  ++iEvent;
 
  // retrieving DCC header
  edm::Handle<EcalRawDataCollection> pDCCHeader;
  const EcalRawDataCollection* DCCHeader = nullptr;
  try {
    e.getByLabel(eventHeaderProducer_, eventHeaderCollection_, pDCCHeader);
    DCCHeader = pDCCHeader.product();
  } catch (std::exception& ex) {
    std::cerr << "Error! can't get the product  retrieving DCC header" << eventHeaderCollection_.c_str() << std::endl;
  }
 
  // retrieving crystal data from Event
  edm::Handle<EBDigiCollection> pEBDigi;
  const EBDigiCollection* EBDigi = nullptr;
  edm::Handle<EEDigiCollection> pEEDigi;
  const EEDigiCollection* EEDigi = nullptr;
 
  if (_ecalPart == "EB") {
    try {
      e.getByLabel(digiProducer_, digiCollection_, pEBDigi);
      EBDigi = pEBDigi.product();
    } catch (std::exception& ex) {
      std::cerr << "Error! can't get the product retrieving EB crystal data " << digiCollection_.c_str() << std::endl;
    }
  } else {
    try {
      e.getByLabel(digiProducer_, digiCollection_, pEEDigi);
      EEDigi = pEEDigi.product();
    } catch (std::exception& ex) {
      std::cerr << "Error! can't get the product retrieving EE crystal data " << digiCollection_.c_str() << std::endl;
    }
  }
 
  // retrieving crystal PN diodes from Event
  edm::Handle<EcalPnDiodeDigiCollection> pPNDigi;
  const EcalPnDiodeDigiCollection* PNDigi = nullptr;
  try {
    e.getByLabel(digiProducer_, digiPNCollection_, pPNDigi);
    PNDigi = pPNDigi.product();
  } catch (std::exception& ex) {
    std::cerr << "Error! can't get the product " << digiPNCollection_.c_str() << std::endl;
  }
 
  // retrieving electronics mapping
  edm::ESHandle<EcalElectronicsMapping> ecalmapping;
  const EcalElectronicsMapping* TheMapping = nullptr;
  try {
    c.get<EcalMappingRcd>().get(ecalmapping);
    TheMapping = ecalmapping.product();
  } catch (std::exception& ex) {
    std::cerr << "Error! can't get the product EcalMappingRcd" << std::endl;
  }
 
  // ====================================
  // Decode Basic DCCHeader Information
  // ====================================
 
  for (EcalRawDataCollection::const_iterator headerItr = DCCHeader->begin(); headerItr != DCCHeader->end();
       ++headerItr) {
    // Get run type and run number
 
    int fed = headerItr->fedId();
 
    if (fed != _fedid && _fedid != -999)
      continue;
 
    runType = headerItr->getRunType();
    runNum = headerItr->getRunNumber();
    event = headerItr->getLV1();
    dccID = headerItr->getDccInTCCCommand();
    fedID = headerItr->fedId();
 
    // take event only if the fed corresponds to the DCC in TCC
    if (600 + dccID != fedID)
      continue;
 
    // Cut on runType
    if (runType != EcalDCCHeaderBlock::LASER_STD && runType != EcalDCCHeaderBlock::LASER_GAP)
      return;
 
    // Define output results files' names
 
    if (IsFileCreated == 0) {
      stringstream namefile2;
      namefile2 << resdir_ << "/APDAmpl_Run" << runNum << "_" << _fedid << "_" << _tower << "_" << _channel << ".root";
      resfile = namefile2.str();
 
      // Get Matacq ttrig
 
      stringstream namefile;
      namefile << resdir_ << "/Matacq-Run" << runNum << ".root";
 
      doesRefFileExist = 0;
 
      FILE* test;
      test = fopen(namefile.str().c_str(), "r");
      if (test)
        doesRefFileExist = 1;
 
      if (doesRefFileExist == 1) {
        matacqFile = new TFile((namefile.str().c_str()));
        matacqTree = (TTree*)matacqFile->Get("MatacqShape");
 
        matacqTree->SetBranchAddress("event", &evtMat);
        matacqTree->SetBranchAddress("color", &colMat);
        matacqTree->SetBranchAddress("peak", &peakMat);
        matacqTree->SetBranchAddress("ttrig", &ttMat);
      }
 
      IsFileCreated = 1;
    }
 
    // Retrieve laser color and event number
 
    EcalDCCHeaderBlock::EcalDCCEventSettings settings = headerItr->getEventSettings();
    int color = settings.wavelength;
 
    vector<int>::iterator iter = find(colors.begin(), colors.end(), color);
    if (iter == colors.end()) {
      colors.push_back(color);
      cout << " new color found " << color << " " << colors.size() << endl;
    }
  }
 
  // cut on fedID
 
  if (fedID != _fedid && _fedid != -999)
    return;
 
  // ======================
  // Decode PN Information
  // ======================
 
  TPNFit* pnfit = new TPNFit();
  pnfit->init(_nsamplesPN, _firstsamplePN, _lastsamplePN);
 
  double chi2pn = 0;
  double ypnrange[50];
  double dsum = 0.;
  double dsum1 = 0.;
  double bl = 0.;
  double bl1 = 0.;
  double val_max = 0.;
  unsigned int samplemax = 0;
  unsigned int k;
 
  std::vector<double> allPNAmpl;
 
  for (EcalPnDiodeDigiCollection::const_iterator pnItr = PNDigi->begin(); pnItr != PNDigi->end();
       ++pnItr) {  // Loop on PNs
 
    for (int samId = 0; samId < (*pnItr).size(); samId++) {  // Loop on PN samples
      pn[samId] = (*pnItr).sample(samId).adc();
    }
 
    for (dsum = 0., k = 0; k < _presamplePN; k++) {
      dsum += pn[k];
    }
    bl = dsum / ((double)_presamplePN);
 
    for (val_max = 0., k = 0; k < _nsamplesPN; k++) {
      ypnrange[k] = pn[k] - bl;
 
      if (ypnrange[k] > val_max) {
        val_max = ypnrange[k];
        samplemax = k;
      }
    }
 
    chi2pn = pnfit->doFit(samplemax, &ypnrange[0]);
 
    if (chi2pn == 101 || chi2pn == 102 || chi2pn == 103)
      pnAmpl = 0.;
    else
      pnAmpl = pnfit->getAmpl();
 
    allPNAmpl.push_back(pnAmpl);
  }
 
  // ===========
  // Get Matacq
  // ===========
 
  ttrig = -1.;
  peak = -1.;
 
  if (doesRefFileExist == 1) {
    // FIXME
    if (color == 0)
      matacqTree->GetEntry(event - 1);
    else if (color == 3)
      matacqTree->GetEntry(matacqTree->GetEntries("color==0") + event - 1);
    ttrig = ttMat;
    peak = peakMat;
  }
 
  // ===========================
  // Decode EBDigis Information
  // ===========================
 
  double yrange[10];
  int adcGain = 0;
  int side = 0;
 
  if (EBDigi) {
    for (EBDigiCollection::const_iterator digiItr = EBDigi->begin(); digiItr != EBDigi->end();
         ++digiItr) {  // Loop on crystals
 
      EBDetId id_crystal(digiItr->id());
      EBDataFrame df(*digiItr);
 
      int etaG = id_crystal.ieta();  // global
      int phiG = id_crystal.iphi();  // global
 
      int etaL;  // local
      int phiL;  // local
      std::pair<int, int> LocalCoord = MEEBGeom::localCoord(etaG, phiG);
 
      etaL = LocalCoord.first;
      phiL = LocalCoord.second;
 
      eta = etaG;
      phi = phiG;
 
      side = MEEBGeom::side(etaG, phiG);
 
      EcalElectronicsId elecid_crystal = TheMapping->getElectronicsId(id_crystal);
 
      int towerID = elecid_crystal.towerId();
      // int channelID=elecid_crystal.channelId()-1;  // FIXME so far for endcap only
      int strip = elecid_crystal.stripId();
      int xtal = elecid_crystal.xtalId();
      int channelID = 5 * (strip - 1) + xtal - 1;  // FIXME
 
      int module = MEEBGeom::lmmod(etaG, phiG);
 
      std::pair<int, int> pnpair = MEEBGeom::pn(module);
      unsigned int MyPn0 = pnpair.first;
      unsigned int MyPn1 = pnpair.second;
 
      unsigned int channel = MEEBGeom::electronic_channel(etaL, phiL);
      assert(channel < nCrys);
 
      double adcmax = 0.0;
 
      if (towerID != int(_tower) || channelID != int(_channel) || dccID != int(_fedid - 600))
        continue;
      else
        channelNumber = channel;
 
      for (unsigned int i = 0; i < (*digiItr).size(); ++i) {  // Loop on adc samples
 
        EcalMGPASample samp_crystal(df.sample(i));
        adc[i] = samp_crystal.adc();
        adcG[i] = samp_crystal.gainId();
 
        if (i == 0)
          adcGain = adcG[i];
        if (i > 0)
          adcGain = TMath::Max(adcG[i], adcGain);
 
        if (adc[i] > adcmax) {
          adcmax = adc[i];
        }
      }
 
      for (dsum = 0., dsum1 = 0., k = 0; k < _presample; k++) {
        dsum += adc[k];
        if (k < _presample - 1)
          dsum1 += adc[k];
      }
 
      bl = dsum / ((double)_presample);
      bl1 = dsum1 / ((double)_presample - 1);
 
      for (val_max = 0., k = 0; k < _nsamples; k++) {
        yrange[k] = adc[k] - bl;
        if (yrange[k] > val_max) {
          val_max = yrange[k];
          samplemax = k;
        }
      }
 
      if (samplemax == 4 || samplemax == 5) {
        val_max = val_max + bl - bl1;
        for (k = 0; k < _nsamples; k++) {
          yrange[k] = yrange[k] + bl - bl1;
        }
      }
 
      for (unsigned int k = 0; k < _nsamples; k++) {
        adc[k] = yrange[k];
      }
 
      pn0 = allPNAmpl[MyPn0];
      pn1 = allPNAmpl[MyPn1];
 
      if (samplemax >= _timingcutlow && samplemax <= _timingcuthigh && lightside == side)
        ADCtrees->Fill();
    }
 
  } else {
    for (EEDigiCollection::const_iterator digiItr = EEDigi->begin(); digiItr != EEDigi->end();
         ++digiItr) {  // Loop on crystals
 
      EEDetId id_crystal(digiItr->id());
      EEDataFrame df(*digiItr);
 
      phi = id_crystal.ix() - 1;
      eta = id_crystal.iy() - 1;
      side = 0;  // FIXME
 
      // Recover the TT id and the electronic crystal numbering from EcalElectronicsMapping
 
      EcalElectronicsId elecid_crystal = TheMapping->getElectronicsId(id_crystal);
 
      int towerID = elecid_crystal.towerId();
      int channelID = elecid_crystal.channelId() - 1;
 
      int module = MEEEGeom::lmmod(phi, eta);
 
      std::pair<int, int> pnpair = MEEEGeom::pn(module, _fedid);
      unsigned int MyPn0 = pnpair.first;
      unsigned int MyPn1 = pnpair.second;
 
      unsigned int channel = MEEEGeom::crystal(phi, eta);
      assert(channel < nCrys);
 
      double adcmax = 0.0;
 
      if (towerID != int(_tower) || channelID != int(_channel) || dccID != int(_fedid - 600))
        continue;
      else
        channelNumber = channel;
 
      for (unsigned int i = 0; i < (*digiItr).size(); ++i) {  // Loop on adc samples
 
        EcalMGPASample samp_crystal(df.sample(i));
        adc[i] = samp_crystal.adc();
        adcG[i] = samp_crystal.gainId();
 
        if (i == 0)
          adcGain = adcG[i];
        if (i > 0)
          adcGain = TMath::Max(adcG[i], adcGain);
 
        if (adc[i] > adcmax) {
          adcmax = adc[i];
        }
      }
 
      for (dsum = 0., dsum1 = 0., k = 0; k < _presample; k++) {
        dsum += adc[k];
        if (k < _presample - 1)
          dsum1 += adc[k];
      }
 
      bl = dsum / ((double)_presample);
      bl1 = dsum1 / ((double)_presample - 1);
 
      for (val_max = 0., k = 0; k < _nsamples; k++) {
        yrange[k] = adc[k] - bl;
        if (yrange[k] > val_max) {
          val_max = yrange[k];
          samplemax = k;
        }
      }
 
      if (samplemax == 4 || samplemax == 5) {
        val_max = val_max + bl - bl1;
        for (k = 0; k < _nsamples; k++) {
          yrange[k] = yrange[k] + bl - bl1;
        }
      }
 
      for (unsigned int k = 0; k < _nsamples; k++) {
        adc[k] = yrange[k];
      }
 
      pn0 = allPNAmpl[MyPn0];
      pn1 = allPNAmpl[MyPn1];
 
      if (samplemax >= _timingcutlow && samplemax <= _timingcuthigh && lightside == side)
        ADCtrees->Fill();
    }
  }
 
}  // analyze
 
//========================================================================
void EcalPerEvtLaserAnalyzer::endJob() {
  //========================================================================
 
  assert(colors.size() <= nColor);
  unsigned int nCol = colors.size();
 
  ADCtrees->Write();
 
  cout << "\n\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << endl;
  cout << "\t+=+       Analyzing laser data: getting per event               +=+" << endl;
  cout << "\t+=+            APD Amplitudes and ADC samples                   +=+" << endl;
  cout << "\t+=+    for fed:" << _fedid << ", tower:" << _tower << ", and channel:" << _channel << endl;
 
  // Define temporary tree to save APD amplitudes
 
  APDFile = new TFile(resfile.c_str(), "RECREATE");
 
  int ieta, iphi, channelID, towerID, flag;
  double alpha, beta;
 
  colors.push_back(color);
 
  for (unsigned int i = 0; i < nCol; i++) {
    stringstream name1;
    name1 << "headerCol" << colors[i];
 
    header[i] = new TTree(name1.str().c_str(), name1.str().c_str());
 
    header[i]->Branch("alpha", &alpha, "alpha/D");
    header[i]->Branch("beta", &beta, "beta/D");
    header[i]->Branch("iphi", &iphi, "iphi/I");
    header[i]->Branch("ieta", &ieta, "ieta/I");
    header[i]->Branch("dccID", &dccID, "dccID/I");
    header[i]->Branch("towerID", &towerID, "towerID/I");
    header[i]->Branch("channelID", &channelID, "channelID/I");
 
    header[i]->SetBranchAddress("alpha", &alpha);
    header[i]->SetBranchAddress("beta", &beta);
    header[i]->SetBranchAddress("iphi", &iphi);
    header[i]->SetBranchAddress("ieta", &ieta);
    header[i]->SetBranchAddress("dccID", &dccID);
    header[i]->SetBranchAddress("towerID", &towerID);
    header[i]->SetBranchAddress("channelID", &channelID);
  }
 
  stringstream name2;
  name2 << "APDTree";
  APDtrees = new TTree(name2.str().c_str(), name2.str().c_str());
 
  //List of branches
 
  APDtrees->Branch("event", &event, "event/I");
  APDtrees->Branch("color", &color, "color/I");
  APDtrees->Branch("adc", &adc, "adc[10]/D");
  APDtrees->Branch("peakMatacq", &peak, "peak/D");
  APDtrees->Branch("ttrigMatacq", &ttrig, "ttrig/D");
  APDtrees->Branch("apdAmpl", &apdAmpl, "apdAmpl/D");
  APDtrees->Branch("apdTime", &apdTime, "apdTime/D");
  APDtrees->Branch("flag", &flag, "flag/I");
  APDtrees->Branch("pn0", &pn0, "pn0/D");
  APDtrees->Branch("pn1", &pn1, "pn1/D");
 
  APDtrees->SetBranchAddress("event", &event);
  APDtrees->SetBranchAddress("color", &color);
  APDtrees->SetBranchAddress("adc", adc);
  APDtrees->SetBranchAddress("peakMatacq", &peak);
  APDtrees->SetBranchAddress("ttrigMatacq", &ttrig);
  APDtrees->SetBranchAddress("apdAmpl", &apdAmpl);
  APDtrees->SetBranchAddress("apdTime", &apdTime);
  APDtrees->SetBranchAddress("flag", &flag);
  APDtrees->SetBranchAddress("pn0", &pn0);
  APDtrees->SetBranchAddress("pn1", &pn1);
 
  // Retrieve alpha and beta for APD amplitudes calculation
 
  TFile* alphaFile = new TFile(refalphabeta_.c_str());
  TTree* alphaTree[2];
 
  Double_t alphaRun, betaRun;
  int ietaRun, iphiRun, channelIDRun, towerIDRun, dccIDRun, flagRun;
 
  for (unsigned int i = 0; i < nCol; i++) {
    stringstream name3;
    name3 << "ABCol" << i;
    alphaTree[i] = (TTree*)alphaFile->Get(name3.str().c_str());
    alphaTree[i]->SetBranchStatus("*", false);
    alphaTree[i]->SetBranchStatus("alpha", true);
    alphaTree[i]->SetBranchStatus("beta", true);
    alphaTree[i]->SetBranchStatus("iphi", true);
    alphaTree[i]->SetBranchStatus("ieta", true);
    alphaTree[i]->SetBranchStatus("dccID", true);
    alphaTree[i]->SetBranchStatus("towerID", true);
    alphaTree[i]->SetBranchStatus("channelID", true);
    alphaTree[i]->SetBranchStatus("flag", true);
 
    alphaTree[i]->SetBranchAddress("alpha", &alphaRun);
    alphaTree[i]->SetBranchAddress("beta", &betaRun);
    alphaTree[i]->SetBranchAddress("iphi", &iphiRun);
    alphaTree[i]->SetBranchAddress("ieta", &ietaRun);
    alphaTree[i]->SetBranchAddress("dccID", &dccIDRun);
    alphaTree[i]->SetBranchAddress("towerID", &towerIDRun);
    alphaTree[i]->SetBranchAddress("channelID", &channelIDRun);
    alphaTree[i]->SetBranchAddress("flag", &flagRun);
  }
 
  PulseFitWithFunction* pslsfit = new PulseFitWithFunction();
 
  double chi2;
 
  for (unsigned int icol = 0; icol < nCol; icol++) {
    IsThereDataADC[icol] = 1;
    stringstream cut;
    cut << "color==" << colors.at(icol);
    if (ADCtrees->GetEntries(cut.str().c_str()) < 10)
      IsThereDataADC[icol] = 0;
    IsHeaderFilled[icol] = 0;
  }
 
  // Define submodule and channel number inside the submodule (as Patrice)
 
  Long64_t nbytes = 0, nb = 0;
  for (Long64_t jentry = 0; jentry < ADCtrees->GetEntriesFast(); jentry++) {  // Loop on events
    nb = ADCtrees->GetEntry(jentry);
    nbytes += nb;
 
    int iCry = channelNumber;
 
    // get back color
 
    unsigned int iCol = 0;
    for (unsigned int i = 0; i < nCol; i++) {
      if (color == colors[i]) {
        iCol = i;
        i = colors.size();
      }
    }
 
    alphaTree[iCol]->GetEntry(iCry);
 
    flag = flagRun;
    iphi = iphiRun;
    ieta = ietaRun;
    towerID = towerIDRun;
    channelID = channelIDRun;
    alpha = alphaRun;
    beta = betaRun;
 
    if (IsHeaderFilled[iCol] == 0) {
      header[iCol]->Fill();
      IsHeaderFilled[iCol] = 1;
    }
    // Amplitude calculation
 
    apdAmpl = 0;
    apdTime = 0;
 
    pslsfit->init(_nsamples, _firstsample, _lastsample, _niter, alphaRun, betaRun);
    chi2 = pslsfit->doFit(&adc[0]);
 
    if (chi2 < 0. || chi2 == 102 || chi2 == 101) {
      apdAmpl = 0;
      apdTime = 0;
 
    } else {
      apdAmpl = pslsfit->getAmpl();
      apdTime = pslsfit->getTime();
    }
 
    APDtrees->Fill();
  }
 
  alphaFile->Close();
 
  ADCFile->Close();
 
  APDFile->Write();
  APDFile->Close();
 
  // Remove unwanted files
 
  stringstream del;
  del << "rm " << ADCfile;
  system(del.str().c_str());
 
  cout << "\t+=+    .................................................. done  +=+" << endl;
  cout << "\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << endl;
}
 
DEFINE_FWK_MODULE(EcalPerEvtLaserAnalyzer);