EcalLaserAnalyzer2.cc

Go to the documentation of this file.

/* 
 *  \class EcalLaserAnalyzer2
 *  primary author: Julie Malcles - CEA/Saclay
 *  author: Gautier Hamel De Monchenault - CEA/Saclay
 */

#include <TAxis.h>
#include <TH1.h>
#include <TProfile.h>
#include <TTree.h>
#include <TChain.h>
#include <TFile.h>
#include <TMath.h>

#include "CalibCalorimetry/EcalLaserAnalyzer/plugins/EcalLaserAnalyzer2.h"

#include <sstream>
#include <iostream>
#include <iomanip>
#include <fstream>

#include <FWCore/MessageLogger/interface/MessageLogger.h>
#include <FWCore/Utilities/interface/Exception.h>

#include <FWCore/Framework/interface/ESHandle.h>
#include <FWCore/Framework/interface/EventSetup.h>

#include <Geometry/EcalMapping/interface/EcalElectronicsMapping.h>
#include <Geometry/EcalMapping/interface/EcalMappingRcd.h>

#include <FWCore/Framework/interface/Event.h>
#include <FWCore/Framework/interface/MakerMacros.h>
#include <FWCore/ParameterSet/interface/ParameterSet.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/TMom.h>
#include <CalibCalorimetry/EcalLaserAnalyzer/interface/TAPD.h>
#include <CalibCalorimetry/EcalLaserAnalyzer/interface/TAPDPulse.h>
#include <CalibCalorimetry/EcalLaserAnalyzer/interface/TPN.h>
#include <CalibCalorimetry/EcalLaserAnalyzer/interface/TPNPulse.h>
#include <CalibCalorimetry/EcalLaserAnalyzer/interface/TPNCor.h>
#include <CalibCalorimetry/EcalLaserAnalyzer/interface/TMem.h>
#include <CalibCalorimetry/EcalLaserAnalyzer/interface/PulseFitWithShape.h>
#include <CalibCalorimetry/EcalLaserAnalyzer/interface/ME.h>
#include <CalibCalorimetry/EcalLaserAnalyzer/interface/MEGeom.h>

using namespace std;

//========================================================================
EcalLaserAnalyzer2::EcalLaserAnalyzer2(const edm::ParameterSet& iConfig)
    //========================================================================
    : iEvent(0),

      // framework parameters with default values
      _nsamples(iConfig.getUntrackedParameter<unsigned int>("nSamples", 10)),
      _presample(iConfig.getUntrackedParameter<unsigned int>("nPresamples", 2)),
      _firstsample(iConfig.getUntrackedParameter<unsigned int>("firstSample", 1)),
      _lastsample(iConfig.getUntrackedParameter<unsigned int>("lastSample", 2)),
      _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", 2)),
      _timingcuthigh(iConfig.getUntrackedParameter<unsigned int>("timingCutHigh", 9)),
      _timingquallow(iConfig.getUntrackedParameter<unsigned int>("timingQualLow", 3)),
      _timingqualhigh(iConfig.getUntrackedParameter<unsigned int>("timingQualHigh", 8)),
      _ratiomincutlow(iConfig.getUntrackedParameter<double>("ratioMinCutLow", 0.4)),
      _ratiomincuthigh(iConfig.getUntrackedParameter<double>("ratioMinCutHigh", 0.95)),
      _ratiomaxcutlow(iConfig.getUntrackedParameter<double>("ratioMaxCutLow", 0.8)),
      _presamplecut(iConfig.getUntrackedParameter<double>("presampleCut", 5.0)),
      _niter(iConfig.getUntrackedParameter<unsigned int>("nIter", 5)),
      _noise(iConfig.getUntrackedParameter<double>("noise", 2.0)),
      _ecalPart(iConfig.getUntrackedParameter<std::string>("ecalPart", "EB")),
      _saveshapes(iConfig.getUntrackedParameter<bool>("saveShapes", true)),
      _docorpn(iConfig.getUntrackedParameter<bool>("doCorPN", false)),
      _fedid(iConfig.getUntrackedParameter<int>("fedID", -999)),
      _saveallevents(iConfig.getUntrackedParameter<bool>("saveAllEvents", false)),
      _qualpercent(iConfig.getUntrackedParameter<double>("qualPercent", 0.2)),
      _debug(iConfig.getUntrackedParameter<int>("debug", 0)),
      nCrys(NCRYSEB),
      nPNPerMod(NPNPERMOD),
      nMod(NMODEB),
      nSides(NSIDES),
      nSamplesShapes(NSAMPSHAPES),
      IsMatacqOK(false),
      runType(-1),
      runNum(0),
      towerID(-1),
      channelID(-1),
      fedID(-1),
      dccID(-1),
      side(2),
      lightside(2),
      iZ(1),
      phi(-1),
      eta(-1),
      event(0),
      color(0),
      pn0(0),
      pn1(0),
      apdAmpl(0),
      apdTime(0),
      pnAmpl(0),
      pnID(-1),
      moduleID(-1),
      flag(0),
      channelIteratorEE(0),
      ShapeCor(0)

//========================================================================

{
  // Initialization from cfg file

  resdir_ = iConfig.getUntrackedParameter<std::string>("resDir");
  elecfile_ = iConfig.getUntrackedParameter<std::string>("elecFile");
  pncorfile_ = iConfig.getUntrackedParameter<std::string>("pnCorFile");

  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");

  // Geometrical constants initialization
  if (_ecalPart == "EB") {
    nCrys = NCRYSEB;
  } else {
    nCrys = NCRYSEE;
  }
  iZ = 1;
  if (_fedid <= 609)
    iZ = -1;
  modules = ME::lmmodFromDcc(_fedid);
  nMod = modules.size();
  nRefChan = NREFCHAN;

  for (unsigned int j = 0; j < nCrys; j++) {
    iEta[j] = -1;
    iPhi[j] = -1;
    iModule[j] = 10;
    iTowerID[j] = -1;
    iChannelID[j] = -1;
    idccID[j] = -1;
    iside[j] = -1;
    wasTimingOK[j] = true;
    wasGainOK[j] = true;
  }

  for (unsigned int j = 0; j < nMod; j++) {
    int ii = modules[j];
    firstChanMod[ii - 1] = 0;
    isFirstChanModFilled[ii - 1] = 0;
  }

  // Quality check flags

  isGainOK = true;
  isTimingOK = true;

  // PN linearity corrector

  pnCorrector = new TPNCor(pncorfile_);

  // Objects dealing with pulses

  APDPulse = new TAPDPulse(_nsamples,
                           _presample,
                           _firstsample,
                           _lastsample,
                           _timingcutlow,
                           _timingcuthigh,
                           _timingquallow,
                           _timingqualhigh,
                           _ratiomincutlow,
                           _ratiomincuthigh,
                           _ratiomaxcutlow);
  PNPulse = new TPNPulse(_nsamplesPN, _presamplePN);

  // Object dealing with MEM numbering

  Mem = new TMem(_fedid);

  // Objects needed for npresample calculation

  Delta01 = new TMom();
  Delta12 = new TMom();
}

//========================================================================
EcalLaserAnalyzer2::~EcalLaserAnalyzer2() {
  //========================================================================

  // do anything here that needs to be done at desctruction time
  // (e.g. close files, deallocate resources etc.)
}

//========================================================================
void EcalLaserAnalyzer2::beginJob() {
  //========================================================================

  // Create temporary files and trees to save adc samples
  //======================================================

  ADCfile = resdir_;
  ADCfile += "/APDSamplesLaser.root";

  APDfile = resdir_;
  APDfile += "/APDPNLaserAllEvents.root";

  ADCFile = new TFile(ADCfile.c_str(), "RECREATE");

  for (unsigned int i = 0; i < nCrys; i++) {
    stringstream name;
    name << "ADCTree" << i + 1;
    ADCtrees[i] = new TTree(name.str().c_str(), name.str().c_str());

    ADCtrees[i]->Branch("iphi", &phi, "phi/I");
    ADCtrees[i]->Branch("ieta", &eta, "eta/I");
    ADCtrees[i]->Branch("towerID", &towerID, "towerID/I");
    ADCtrees[i]->Branch("channelID", &channelID, "channelID/I");
    ADCtrees[i]->Branch("dccID", &dccID, "dccID/I");
    ADCtrees[i]->Branch("side", &side, "side/I");
    ADCtrees[i]->Branch("event", &event, "event/I");
    ADCtrees[i]->Branch("color", &color, "color/I");
    ADCtrees[i]->Branch("adc", &adc, "adc[10]/D");
    ADCtrees[i]->Branch("pn0", &pn0, "pn0/D");
    ADCtrees[i]->Branch("pn1", &pn1, "pn1/D");

    ADCtrees[i]->SetBranchAddress("ieta", &eta);
    ADCtrees[i]->SetBranchAddress("iphi", &phi);
    ADCtrees[i]->SetBranchAddress("towerID", &towerID);
    ADCtrees[i]->SetBranchAddress("channelID", &channelID);
    ADCtrees[i]->SetBranchAddress("dccID", &dccID);
    ADCtrees[i]->SetBranchAddress("side", &side);
    ADCtrees[i]->SetBranchAddress("event", &event);
    ADCtrees[i]->SetBranchAddress("color", &color);
    ADCtrees[i]->SetBranchAddress("adc", adc);
    ADCtrees[i]->SetBranchAddress("pn0", &pn0);
    ADCtrees[i]->SetBranchAddress("pn1", &pn1);
  }

  // Define output results filenames
  //==================================
  stringstream namefile1;
  namefile1 << resdir_ << "/SHAPE_LASER.root";
  shapefile = namefile1.str();

  stringstream namefile2;
  namefile2 << resdir_ << "/APDPN_LASER.root";
  resfile = namefile2.str();

  stringstream namefile3;
  namefile3 << resdir_ << "/MATACQ.root";

  matfile = namefile3.str();

  // Get Pulse Shapes
  //==================

  IsMatacqOK = getShapes();
  if (!IsMatacqOK) {
    cout << " ERROR! No matacq shape available: analysis aborted !" << endl;
    return;
  }

  // Laser events counter

  laserEvents = 0;
}

//========================================================================
void EcalLaserAnalyzer2::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 if (_ecalPart == "EE") {
    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;
    }
  } else {
    cout << " Wrong ecalPart in cfg file " << endl;
    return;
  }

  // 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();
    lightside = headerItr->getRtHalf();

    // Check fed corresponds to the DCC in TCC

    if (600 + dccID != fedID)
      continue;

    // Cut on runType

    if (runType != EcalDCCHeaderBlock::LASER_STD && runType != EcalDCCHeaderBlock::LASER_GAP &&
        runType != EcalDCCHeaderBlock::LASER_POWER_SCAN && runType != EcalDCCHeaderBlock::LASER_DELAY_SCAN)
      return;

    // Retrieve laser color and event number

    EcalDCCHeaderBlock::EcalDCCEventSettings settings = headerItr->getEventSettings();
    color = settings.wavelength;
    if (color < 0)
      return;

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

  // Check Matacq shape exists

  if (!IsMatacqOK)
    return;

  // Cut on fedID

  if (fedID != _fedid && _fedid != -999)
    return;

  // Count laser events

  laserEvents++;

  // ======================
  // Decode PN Information
  // ======================

  TPNFit* pnfit = new TPNFit();
  pnfit->init(_nsamplesPN, _firstsamplePN, _lastsamplePN);

  double chi2pn = 0;
  unsigned int samplemax = 0;
  int pnGain = 0;

  map<int, vector<double> > allPNAmpl;
  map<int, vector<double> > allPNGain;

  // Loop on PNs digis

  for (EcalPnDiodeDigiCollection::const_iterator pnItr = PNDigi->begin(); pnItr != PNDigi->end(); ++pnItr) {
    EcalPnDiodeDetId pnDetId = EcalPnDiodeDetId((*pnItr).id());

    if (_debug == 1)
      cout << "-- debug test -- Inside PNDigi - pnID=" << pnDetId.iPnId() << ", dccID=" << pnDetId.iDCCId() << endl;

    // Skip MEM DCC without relevant data

    bool isMemRelevant = Mem->isMemRelevant(pnDetId.iDCCId());
    if (!isMemRelevant)
      continue;

    // Loop on PN samples

    for (int samId = 0; samId < (*pnItr).size(); samId++) {
      pn[samId] = (*pnItr).sample(samId).adc();
      pnG[samId] = (*pnItr).sample(samId).gainId();
      if (samId == 0)
        pnGain = pnG[samId];
      if (samId > 0)
        pnGain = int(TMath::Max(pnG[samId], pnGain));
    }

    if (pnGain != 1)
      cout << "PN gain different from 1" << endl;

    // Calculate amplitude from pulse

    PNPulse->setPulse(pn);
    pnNoPed = PNPulse->getAdcWithoutPedestal();
    samplemax = PNPulse->getMaxSample();
    chi2pn = pnfit->doFit(samplemax, &pnNoPed[0]);
    if (chi2pn == 101 || chi2pn == 102 || chi2pn == 103)
      pnAmpl = 0.;
    else
      pnAmpl = pnfit->getAmpl();

    // Apply linearity correction

    double corr = 1.0;
    if (_docorpn)
      corr = pnCorrector->getPNCorrectionFactor(pnAmpl, pnGain);
    pnAmpl *= corr;

    // Fill PN ampl vector

    allPNAmpl[pnDetId.iDCCId()].push_back(pnAmpl);

    if (_debug == 1)
      cout << "-- debug -- Inside PNDigi - PNampl=" << pnAmpl << ", PNgain=" << pnGain << endl;
  }

  // ===========================
  // Decode EBDigis Information
  // ===========================

  int adcGain = 0;

  if (EBDigi) {
    // Loop on crystals
    //===================

    for (EBDigiCollection::const_iterator digiItr = EBDigi->begin(); digiItr != EBDigi->end();
         ++digiItr) {  // Loop on crystals

      // Retrieve geometry
      //===================

      EBDetId id_crystal(digiItr->id());
      EBDataFrame df(*digiItr);
      EcalElectronicsId elecid_crystal = TheMapping->getElectronicsId(id_crystal);

      int etaG = id_crystal.ieta();  // global
      int phiG = id_crystal.iphi();  // global

      std::pair<int, int> LocalCoord = MEEBGeom::localCoord(etaG, phiG);

      int etaL = LocalCoord.first;   // local
      int phiL = LocalCoord.second;  // local

      int strip = elecid_crystal.stripId();
      int xtal = elecid_crystal.xtalId();

      int module = MEEBGeom::lmmod(etaG, phiG);
      int tower = elecid_crystal.towerId();

      int apdRefTT = MEEBGeom::apdRefTower(module);

      std::pair<int, int> pnpair = MEEBGeom::pn(module);
      unsigned int MyPn0 = pnpair.first;
      unsigned int MyPn1 = pnpair.second;

      int lmr = MEEBGeom::lmr(etaG, phiG);
      unsigned int channel = MEEBGeom::electronic_channel(etaL, phiL);
      assert(channel < nCrys);

      setGeomEB(etaG, phiG, module, tower, strip, xtal, apdRefTT, channel, lmr);

      if (_debug == 1)
        cout << "-- debug -- Inside EBDigi - towerID:" << towerID << " channelID:" << channelID << " module:" << module
             << " modules:" << modules.size() << endl;

      // APD Pulse
      //===========

      // Loop on adc samples

      for (unsigned int i = 0; i < (*digiItr).size(); ++i) {
        EcalMGPASample samp_crystal(df.sample(i));
        adc[i] = samp_crystal.adc();
        adcG[i] = samp_crystal.gainId();
        adc[i] *= adcG[i];
        if (i == 0)
          adcGain = adcG[i];
        if (i > 0)
          adcGain = TMath::Max(adcG[i], adcGain);
      }

      APDPulse->setPulse(adc);

      // Quality checks
      //================

      if (adcGain != 1)
        nEvtBadGain[channel]++;
      if (!APDPulse->isTimingQualOK())
        nEvtBadTiming[channel]++;
      nEvtTot[channel]++;

      // Associate PN ampl
      //===================

      int mem0 = Mem->Mem(lmr, 0);
      int mem1 = Mem->Mem(lmr, 1);

      if (allPNAmpl[mem0].size() > MyPn0)
        pn0 = allPNAmpl[mem0][MyPn0];
      else
        pn0 = 0;
      if (allPNAmpl[mem1].size() > MyPn1)
        pn1 = allPNAmpl[mem1][MyPn1];
      else
        pn1 = 0;

      // Fill if Pulse is fine
      //=======================

      if (APDPulse->isPulseOK() && lightside == side) {
        ADCtrees[channel]->Fill();

        Delta01->addEntry(APDPulse->getDelta(0, 1));
        Delta12->addEntry(APDPulse->getDelta(1, 2));
      }
    }
  } else if (EEDigi) {
    // Loop on crystals
    //===================

    for (EEDigiCollection::const_iterator digiItr = EEDigi->begin(); digiItr != EEDigi->end(); ++digiItr) {
      // Retrieve geometry
      //===================

      EEDetId id_crystal(digiItr->id());
      EEDataFrame df(*digiItr);
      EcalElectronicsId elecid_crystal = TheMapping->getElectronicsId(id_crystal);

      int etaG = id_crystal.iy();
      int phiG = id_crystal.ix();

      int iX = (phiG - 1) / 5 + 1;
      int iY = (etaG - 1) / 5 + 1;

      int tower = elecid_crystal.towerId();
      int ch = elecid_crystal.channelId() - 1;

      int module = MEEEGeom::lmmod(iX, iY);
      if (module >= 18 && side == 1)
        module += 2;
      int lmr = MEEEGeom::lmr(iX, iY, iZ);
      int dee = MEEEGeom::dee(lmr);
      int apdRefTT = MEEEGeom::apdRefTower(lmr, module);

      std::pair<int, int> pnpair = MEEEGeom::pn(dee, module);
      unsigned int MyPn0 = pnpair.first;
      unsigned int MyPn1 = pnpair.second;

      int hashedIndex = 100000 * eta + phi;
      if (channelMapEE.count(hashedIndex) == 0) {
        channelMapEE[hashedIndex] = channelIteratorEE;
        channelIteratorEE++;
      }
      unsigned int channel = channelMapEE[hashedIndex];
      assert(channel < nCrys);

      setGeomEE(etaG, phiG, iX, iY, iZ, module, tower, ch, apdRefTT, channel, lmr);

      if (_debug == 1)
        cout << "-- debug -- Inside EEDigi - towerID:" << towerID << " channelID:" << channelID << " module:" << module
             << " modules:" << modules.size() << endl;

      // APD Pulse
      //===========

      if ((*digiItr).size() > 10)
        cout << "SAMPLES SIZE > 10!" << (*digiItr).size() << endl;

      // Loop on adc samples

      for (unsigned int i = 0; i < (*digiItr).size(); ++i) {
        EcalMGPASample samp_crystal(df.sample(i));
        adc[i] = samp_crystal.adc();
        adcG[i] = samp_crystal.gainId();
        adc[i] *= adcG[i];

        if (i == 0)
          adcGain = adcG[i];
        if (i > 0)
          adcGain = TMath::Max(adcG[i], adcGain);
      }

      APDPulse->setPulse(adc);

      // Quality checks
      //================

      if (adcGain != 1)
        nEvtBadGain[channel]++;
      if (!APDPulse->isTimingQualOK())
        nEvtBadTiming[channel]++;
      nEvtTot[channel]++;

      // Associate PN ampl
      //===================

      int mem0 = Mem->Mem(lmr, 0);
      int mem1 = Mem->Mem(lmr, 1);

      if (allPNAmpl[mem0].size() > MyPn0)
        pn0 = allPNAmpl[mem0][MyPn0];
      else
        pn0 = 0;
      if (allPNAmpl[mem1].size() > MyPn1)
        pn1 = allPNAmpl[mem1][MyPn1];
      else
        pn1 = 0;

      // Fill if Pulse is fine
      //=======================

      if (APDPulse->isPulseOK() && lightside == side) {
        ADCtrees[channel]->Fill();

        Delta01->addEntry(APDPulse->getDelta(0, 1));
        Delta12->addEntry(APDPulse->getDelta(1, 2));
      }
    }
  }
}
// analyze

//========================================================================
void EcalLaserAnalyzer2::endJob() {
  //========================================================================

  if (!IsMatacqOK) {
    cout << "\n\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << endl;
    cout << "\t+=+     WARNING! NO MATACQ        +=+" << endl;
    cout << "\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << endl;
    return;
  }

  // Adjust channel numbers for EE
  //===============================

  if (_ecalPart == "EE") {
    nCrys = channelMapEE.size();
  }

  // Set presamples number
  //======================

  double delta01 = Delta01->getMean();
  double delta12 = Delta12->getMean();
  if (delta12 > _presamplecut) {
    _presample = 2;
    if (delta01 > _presamplecut)
      _presample = 1;
  }

  APDPulse->setPresamples(_presample);

  // Don't do anything if there is no events
  //=========================================

  if (laserEvents == 0) {
    ADCFile->Close();
    stringstream del;
    del << "rm " << ADCfile;
    system(del.str().c_str());
    cout << " No Laser Events " << endl;
    return;
  }

  // Set quality flags for gains and timing
  //=========================================

  double BadGainEvtPercentage = 0.0;
  double BadTimingEvtPercentage = 0.0;

  int nChanBadGain = 0;
  int nChanBadTiming = 0;

  for (unsigned int i = 0; i < nCrys; i++) {
    if (nEvtTot[i] != 0) {
      BadGainEvtPercentage = double(nEvtBadGain[i]) / double(nEvtTot[i]);
      BadTimingEvtPercentage = double(nEvtBadTiming[i]) / double(nEvtTot[i]);
    }
    if (BadGainEvtPercentage > _qualpercent) {
      wasGainOK[i] = false;
      nChanBadGain++;
    }
    if (BadTimingEvtPercentage > _qualpercent) {
      wasTimingOK[i] = false;
      nChanBadTiming++;
    }
  }

  double BadGainChanPercentage = double(nChanBadGain) / double(nCrys);
  double BadTimingChanPercentage = double(nChanBadTiming) / double(nCrys);

  if (BadGainChanPercentage > _qualpercent)
    isGainOK = false;
  if (BadTimingChanPercentage > _qualpercent)
    isTimingOK = false;

  // Analyze adc samples to get amplitudes
  //=======================================

  cout << "\n\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << endl;
  cout << "\t+=+     Analyzing laser data: getting APD, PN, APD/PN, PN/PN    +=+" << endl;

  if (!isGainOK)
    cout << "\t+=+ ............................ WARNING! APD GAIN WAS NOT 1    +=+" << endl;
  if (!isTimingOK)
    cout << "\t+=+ ............................ WARNING! TIMING WAS BAD        +=+" << endl;

  APDFile = new TFile(APDfile.c_str(), "RECREATE");

  int ieta, iphi;

  int flagfit;
  for (unsigned int i = 0; i < nCrys; i++) {
    stringstream name;
    name << "APDTree" << i + 1;

    APDtrees[i] = new TTree(name.str().c_str(), name.str().c_str());

    //List of branches

    APDtrees[i]->Branch("event", &event, "event/I");
    APDtrees[i]->Branch("color", &color, "color/I");
    APDtrees[i]->Branch("iphi", &iphi, "iphi/I");
    APDtrees[i]->Branch("ieta", &ieta, "ieta/I");
    APDtrees[i]->Branch("side", &side, "side/I");
    APDtrees[i]->Branch("dccID", &dccID, "dccID/I");
    APDtrees[i]->Branch("towerID", &towerID, "towerID/I");
    APDtrees[i]->Branch("channelID", &channelID, "channelID/I");
    APDtrees[i]->Branch("apdAmpl", &apdAmpl, "apdAmpl/D");
    APDtrees[i]->Branch("apdTime", &apdTime, "apdTime/D");
    if (_saveallevents)
      APDtrees[i]->Branch("adc", &adc, "adc[10]/D");
    APDtrees[i]->Branch("flagfit", &flagfit, "flagfit/I");
    APDtrees[i]->Branch("pn0", &pn0, "pn0/D");
    APDtrees[i]->Branch("pn1", &pn1, "pn1/D");

    APDtrees[i]->SetBranchAddress("event", &event);
    APDtrees[i]->SetBranchAddress("color", &color);
    APDtrees[i]->SetBranchAddress("iphi", &iphi);
    APDtrees[i]->SetBranchAddress("ieta", &ieta);
    APDtrees[i]->SetBranchAddress("side", &side);
    APDtrees[i]->SetBranchAddress("dccID", &dccID);
    APDtrees[i]->SetBranchAddress("towerID", &towerID);
    APDtrees[i]->SetBranchAddress("channelID", &channelID);
    APDtrees[i]->SetBranchAddress("apdAmpl", &apdAmpl);
    APDtrees[i]->SetBranchAddress("apdTime", &apdTime);
    if (_saveallevents)
      APDtrees[i]->SetBranchAddress("adc", adc);
    APDtrees[i]->SetBranchAddress("flagfit", &flagfit);
    APDtrees[i]->SetBranchAddress("pn0", &pn0);
    APDtrees[i]->SetBranchAddress("pn1", &pn1);
  }

  for (unsigned int iref = 0; iref < nRefChan; iref++) {
    for (unsigned int imod = 0; imod < nMod; imod++) {
      int jmod = modules[imod];

      stringstream nameref;
      nameref << "refAPDTree" << imod << "_" << iref;

      RefAPDtrees[iref][jmod] = new TTree(nameref.str().c_str(), nameref.str().c_str());

      RefAPDtrees[iref][jmod]->Branch("eventref", &eventref, "eventref/I");
      RefAPDtrees[iref][jmod]->Branch("colorref", &colorref, "colorref/I");
      if (iref == 0)
        RefAPDtrees[iref][jmod]->Branch("apdAmplA", &apdAmplA, "apdAmplA/D");
      if (iref == 1)
        RefAPDtrees[iref][jmod]->Branch("apdAmplB", &apdAmplB, "apdAmplB/D");

      RefAPDtrees[iref][jmod]->SetBranchAddress("eventref", &eventref);
      RefAPDtrees[iref][jmod]->SetBranchAddress("colorref", &colorref);
      if (iref == 0)
        RefAPDtrees[iref][jmod]->SetBranchAddress("apdAmplA", &apdAmplA);
      if (iref == 1)
        RefAPDtrees[iref][jmod]->SetBranchAddress("apdAmplB", &apdAmplB);
    }
  }

  assert(colors.size() <= nColor);
  unsigned int nCol = colors.size();

  // Declare PN stuff
  //===================

  for (unsigned int iM = 0; iM < nMod; iM++) {
    unsigned int iMod = modules[iM] - 1;
    for (unsigned int ich = 0; ich < nPNPerMod; ich++) {
      for (unsigned int icol = 0; icol < nCol; icol++) {
        PNFirstAnal[iMod][ich][icol] = new TPN(ich);
        PNAnal[iMod][ich][icol] = new TPN(ich);
      }
    }
  }

  // Declare function for APD ampl fit
  //===================================

  PulseFitWithShape* psfit = new PulseFitWithShape();

  for (unsigned int iCry = 0; iCry < nCrys; iCry++) {
    for (unsigned int icol = 0; icol < nCol; icol++) {
      // Declare APD stuff
      //===================

      APDFirstAnal[iCry][icol] = new TAPD();
      IsThereDataADC[iCry][icol] = 1;
      stringstream cut;
      cut << "color==" << colors.at(icol);
      if (ADCtrees[iCry]->GetEntries(cut.str().c_str()) < 10)
        IsThereDataADC[iCry][icol] = 0;
    }

    unsigned int iMod = iModule[iCry] - 1;
    assert(iMod <= nMod);

    if (isSPRFine)
      psfit->init(_nsamples, _firstsample, _lastsample, _niter, nSamplesShapes, shapesVec, _noise);

    // Loop on events
    //================

    Long64_t nbytes = 0, nb = 0;
    for (Long64_t jentry = 0; jentry < ADCtrees[iCry]->GetEntriesFast(); jentry++) {  // Loop on events
      nb = ADCtrees[iCry]->GetEntry(jentry);
      nbytes += nb;

      flagfit = 1;
      apdAmpl = 0.0;
      apdTime = 0.0;
      ieta = eta;
      iphi = phi;

      // Get back color

      unsigned int iCol = 0;
      for (unsigned int i = 0; i < nCol; i++) {
        if (color == colors[i]) {
          iCol = i;
          i = colors.size();
        }
      }

      // Amplitude calculation

      APDPulse->setPulse(adc);
      adcNoPed = APDPulse->getAdcWithoutPedestal();

      if (isSPRFine && APDPulse->isPulseOK()) {
        psfit->doFit(&adcNoPed[0]);
        apdAmpl = psfit->getAmpl();
        apdTime = psfit->getTime();

      } else {
        apdAmpl = 0;
        apdTime = 0;
        flagfit = 0;
      }

      if (_debug >= 1)
        cout << "-- debug test -- endJob -- apdAmpl:" << apdAmpl << " apdTime:" << apdTime << endl;

      double pnmean;
      if (pn0 < 10 && pn1 > 10) {
        pnmean = pn1;
      } else if (pn1 < 10 && pn0 > 10) {
        pnmean = pn0;
      } else
        pnmean = 0.5 * (pn0 + pn1);

      if (_debug >= 1)
        cout << "-- debug test -- endJob -- pnMean:" << pnmean << endl;

      // Fill PN stuff
      //===============

      if (firstChanMod[iMod] == iCry && IsThereDataADC[iCry][iCol] == 1) {
        for (unsigned int ichan = 0; ichan < nPNPerMod; ichan++) {
          PNFirstAnal[iMod][ichan][iCol]->addEntry(pnmean, pn0, pn1);
        }
      }

      // Fill APD stuff
      //================

      if (apdAmpl != 0.0)
        APDFirstAnal[iCry][iCol]->addEntry(apdAmpl, pnmean, pn0, pn1, apdTime);
      if (_debug >= 1)
        cout << "-- debug test -- endJob -- filling APDTree" << endl;

      APDtrees[iCry]->Fill();

      // Fill reference trees
      //=====================

      if (apdRefMap[0][iMod + 1] == iCry || apdRefMap[1][iMod + 1] == iCry) {
        apdAmplA = 0.0;
        apdAmplB = 0.0;
        eventref = event;
        colorref = color;

        for (unsigned int ir = 0; ir < nRefChan; ir++) {
          if (_debug >= 1)
            cout << "-- debug test -- ir:" << ir << " tt:" << towerID << " refmap:" << apdRefMap[ir][iMod + 1]
                 << " iCry:" << iCry << endl;

          if (apdRefMap[ir][iMod + 1] == iCry) {
            if (_debug >= 1)
              cout << "-- debug test -- cut passed " << endl;
            if (ir == 0)
              apdAmplA = apdAmpl;
            else if (ir == 1)
              apdAmplB = apdAmpl;
            if (_debug >= 1)
              cout << "-- debug test -- apdAmplA=" << apdAmplA << endl;
            if (_debug >= 1)
              cout << "-- debug test -- apdAmplB=" << apdAmplB << endl;
            if (_debug >= 1)
              cout << "-- debug test -- color=" << color << ", event:" << event << ", ir:" << ir
                   << " tt-1:" << towerID - 1 << endl;

            RefAPDtrees[ir][iMod + 1]->Fill();

            if (_debug >= 1)
              cout << "-- debug test -- tree filled" << event << endl;
          }
        }
      }
    }
  }

  delete psfit;

  ADCFile->Close();

  if (_debug == 1)
    cout << "-- debug test -- endJob -- after apdAmpl Loop" << endl;

  // Remove temporary file
  //=======================

  stringstream del;
  del << "rm " << ADCfile;
  system(del.str().c_str());

  // Create output trees
  //=====================

  resFile = new TFile(resfile.c_str(), "RECREATE");

  for (unsigned int iColor = 0; iColor < nCol; iColor++) {
    stringstream nametree;
    nametree << "APDCol" << colors.at(iColor);
    stringstream nametree2;
    nametree2 << "PNCol" << colors.at(iColor);

    restrees[iColor] = new TTree(nametree.str().c_str(), nametree.str().c_str());
    respntrees[iColor] = new TTree(nametree2.str().c_str(), nametree2.str().c_str());

    restrees[iColor]->Branch("iphi", &iphi, "iphi/I");
    restrees[iColor]->Branch("ieta", &ieta, "ieta/I");
    restrees[iColor]->Branch("side", &side, "side/I");
    restrees[iColor]->Branch("dccID", &dccID, "dccID/I");
    restrees[iColor]->Branch("moduleID", &moduleID, "moduleID/I");
    restrees[iColor]->Branch("towerID", &towerID, "towerID/I");
    restrees[iColor]->Branch("channelID", &channelID, "channelID/I");
    restrees[iColor]->Branch("APD", &APD, "APD[6]/D");
    restrees[iColor]->Branch("Time", &Time, "Time[6]/D");
    restrees[iColor]->Branch("APDoPN", &APDoPN, "APDoPN[6]/D");
    restrees[iColor]->Branch("APDoPNA", &APDoPNA, "APDoPNA[6]/D");
    restrees[iColor]->Branch("APDoPNB", &APDoPNB, "APDoPNB[6]/D");
    restrees[iColor]->Branch("APDoAPDA", &APDoAPDA, "APDoAPDA[6]/D");
    restrees[iColor]->Branch("APDoAPDB", &APDoAPDB, "APDoAPDB[6]/D");
    restrees[iColor]->Branch("ShapeCor", &ShapeCor, "ShapeCor/D");
    restrees[iColor]->Branch("flag", &flag, "flag/I");

    respntrees[iColor]->Branch("moduleID", &moduleID, "moduleID/I");
    respntrees[iColor]->Branch("pnID", &pnID, "pnID/I");
    respntrees[iColor]->Branch("PN", &PN, "PN[6]/D");
    respntrees[iColor]->Branch("PNoPN", &PNoPN, "PNoPN[6]/D");
    respntrees[iColor]->Branch("PNoPNA", &PNoPNA, "PNoPNA[6]/D");
    respntrees[iColor]->Branch("PNoPNB", &PNoPNB, "PNoPNB[6]/D");

    restrees[iColor]->SetBranchAddress("iphi", &iphi);
    restrees[iColor]->SetBranchAddress("ieta", &ieta);
    restrees[iColor]->SetBranchAddress("dccID", &dccID);
    restrees[iColor]->SetBranchAddress("moduleID", &moduleID);
    restrees[iColor]->SetBranchAddress("towerID", &towerID);
    restrees[iColor]->SetBranchAddress("channelID", &channelID);
    restrees[iColor]->SetBranchAddress("APD", APD);
    restrees[iColor]->SetBranchAddress("Time", Time);
    restrees[iColor]->SetBranchAddress("APDoPN", APDoPN);
    restrees[iColor]->SetBranchAddress("APDoPNA", APDoPNA);
    restrees[iColor]->SetBranchAddress("APDoPNB", APDoPNB);
    restrees[iColor]->SetBranchAddress("APDoAPDA", APDoAPDA);
    restrees[iColor]->SetBranchAddress("APDoAPDB", APDoAPDB);
    restrees[iColor]->SetBranchAddress("ShapeCor", &ShapeCor);
    restrees[iColor]->SetBranchAddress("flag", &flag);

    respntrees[iColor]->SetBranchAddress("moduleID", &moduleID);
    respntrees[iColor]->SetBranchAddress("pnID", &pnID);
    respntrees[iColor]->SetBranchAddress("PN", PN);
    respntrees[iColor]->SetBranchAddress("PNoPN", PNoPN);
    respntrees[iColor]->SetBranchAddress("PNoPNA", PNoPNA);
    respntrees[iColor]->SetBranchAddress("PNoPNB", PNoPNB);
  }

  // Set Cuts for PN stuff
  //=======================

  for (unsigned int iM = 0; iM < nMod; iM++) {
    unsigned int iMod = modules[iM] - 1;

    for (unsigned int ich = 0; ich < nPNPerMod; ich++) {
      for (unsigned int icol = 0; icol < nCol; icol++) {
        PNAnal[iMod][ich][icol]->setPNCut(PNFirstAnal[iMod][ich][icol]->getPN().at(0),
                                          PNFirstAnal[iMod][ich][icol]->getPN().at(1));
      }
    }
  }

  // Build ref trees indexes
  //========================

  for (unsigned int imod = 0; imod < nMod; imod++) {
    int jmod = modules[imod];
    if (RefAPDtrees[0][jmod]->GetEntries() != 0 && RefAPDtrees[1][jmod]->GetEntries() != 0) {
      RefAPDtrees[0][jmod]->BuildIndex("eventref");
      RefAPDtrees[1][jmod]->BuildIndex("eventref");
    }
  }

  // Final loop on crystals
  //=======================

  for (unsigned int iCry = 0; iCry < nCrys; iCry++) {
    unsigned int iMod = iModule[iCry] - 1;

    // Set cuts on APD stuff
    //=======================

    for (unsigned int iCol = 0; iCol < nCol; iCol++) {
      std::vector<double> lowcut;
      std::vector<double> highcut;
      double cutMin;
      double cutMax;

      cutMin = APDFirstAnal[iCry][iCol]->getAPD().at(0) - 2.0 * APDFirstAnal[iCry][iCol]->getAPD().at(1);
      if (cutMin < 0)
        cutMin = 0;
      cutMax = APDFirstAnal[iCry][iCol]->getAPD().at(0) + 2.0 * APDFirstAnal[iCry][iCol]->getAPD().at(1);

      lowcut.push_back(cutMin);
      highcut.push_back(cutMax);

      cutMin = APDFirstAnal[iCry][iCol]->getTime().at(0) - 2.0 * APDFirstAnal[iCry][iCol]->getTime().at(1);
      cutMax = APDFirstAnal[iCry][iCol]->getTime().at(0) + 2.0 * APDFirstAnal[iCry][iCol]->getTime().at(1);
      lowcut.push_back(cutMin);
      highcut.push_back(cutMax);

      APDAnal[iCry][iCol] = new TAPD();
      APDAnal[iCry][iCol]->setAPDCut(APDFirstAnal[iCry][iCol]->getAPD().at(0),
                                     APDFirstAnal[iCry][iCol]->getAPD().at(1));
      APDAnal[iCry][iCol]->setAPDoPNCut(APDFirstAnal[iCry][iCol]->getAPDoPN().at(0),
                                        APDFirstAnal[iCry][iCol]->getAPDoPN().at(1));
      APDAnal[iCry][iCol]->setAPDoPN0Cut(APDFirstAnal[iCry][iCol]->getAPDoPN0().at(0),
                                         APDFirstAnal[iCry][iCol]->getAPDoPN0().at(1));
      APDAnal[iCry][iCol]->setAPDoPN1Cut(APDFirstAnal[iCry][iCol]->getAPDoPN1().at(0),
                                         APDFirstAnal[iCry][iCol]->getAPDoPN1().at(1));
      APDAnal[iCry][iCol]->setTimeCut(APDFirstAnal[iCry][iCol]->getTime().at(0),
                                      APDFirstAnal[iCry][iCol]->getTime().at(1));

      APDAnal[iCry][iCol]->set2DAPDoAPD0Cut(lowcut, highcut);
      APDAnal[iCry][iCol]->set2DAPDoAPD1Cut(lowcut, highcut);
    }

    // Final loop on events
    //=======================

    Long64_t nbytes = 0, nb = 0;
    for (Long64_t jentry = 0; jentry < APDtrees[iCry]->GetEntriesFast(); jentry++) {
      nb = APDtrees[iCry]->GetEntry(jentry);
      nbytes += nb;

      double pnmean;
      if (pn0 < 10 && pn1 > 10) {
        pnmean = pn1;
      } else if (pn1 < 10 && pn0 > 10) {
        pnmean = pn0;
      } else
        pnmean = 0.5 * (pn0 + pn1);

      // Get back color
      //===============

      unsigned int iCol = 0;
      for (unsigned int i = 0; i < nCol; i++) {
        if (color == colors[i]) {
          iCol = i;
          i = colors.size();
        }
      }

      // Fill PN stuff
      //===============

      if (firstChanMod[iMod] == iCry && IsThereDataADC[iCry][iCol] == 1) {
        for (unsigned int ichan = 0; ichan < nPNPerMod; ichan++) {
          PNAnal[iMod][ichan][iCol]->addEntry(pnmean, pn0, pn1);
        }
      }

      // Get ref amplitudes
      //===================

      if (_debug >= 1)
        cout << "-- debug test -- LastLoop event:" << event << " apdAmpl:" << apdAmpl << endl;
      apdAmplA = 0.0;
      apdAmplB = 0.0;

      for (unsigned int iRef = 0; iRef < nRefChan; iRef++) {
        RefAPDtrees[iRef][iMod + 1]->GetEntryWithIndex(event);
      }

      if (_debug == 1)
        cout << "-- debug test -- LastLoop apdAmplA:" << apdAmplA << " apdAmplB:" << apdAmplB << ", event:" << event
             << ", eventref:" << eventref << endl;

      // Fill APD stuff
      //===============

      APDAnal[iCry][iCol]->addEntry(apdAmpl, pnmean, pn0, pn1, apdTime, apdAmplA, apdAmplB);
    }

    moduleID = iMod + 1;

    if (moduleID >= 20)
      moduleID -= 2;  // Trick to fix endcap specificity

    // Get final results for APD
    //===========================

    for (unsigned int iColor = 0; iColor < nCol; iColor++) {
      std::vector<double> apdvec = APDAnal[iCry][iColor]->getAPD();
      std::vector<double> apdpnvec = APDAnal[iCry][iColor]->getAPDoPN();
      std::vector<double> apdpn0vec = APDAnal[iCry][iColor]->getAPDoPN0();
      std::vector<double> apdpn1vec = APDAnal[iCry][iColor]->getAPDoPN1();
      std::vector<double> timevec = APDAnal[iCry][iColor]->getTime();
      std::vector<double> apdapd0vec = APDAnal[iCry][iColor]->getAPDoAPD0();
      std::vector<double> apdapd1vec = APDAnal[iCry][iColor]->getAPDoAPD1();

      for (unsigned int i = 0; i < apdvec.size(); i++) {
        APD[i] = apdvec.at(i);
        APDoPN[i] = apdpnvec.at(i);
        APDoPNA[i] = apdpn0vec.at(i);
        APDoPNB[i] = apdpn1vec.at(i);
        APDoAPDA[i] = apdapd0vec.at(i);
        APDoAPDB[i] = apdapd1vec.at(i);
        Time[i] = timevec.at(i);
        ShapeCor = shapeCorrection;
      }

      // Fill APD results trees
      //========================

      iphi = iPhi[iCry];
      ieta = iEta[iCry];
      dccID = idccID[iCry];
      towerID = iTowerID[iCry];
      side = iside[iCry];
      channelID = iChannelID[iCry];

      if (!wasGainOK[iCry] || !wasTimingOK[iCry] || IsThereDataADC[iCry][iColor] == 0)
        flag = 1;
      else
        flag = 0;

      if (_debug >= 1)
        cout << "-- debug test -- endJob -- APD[0]" << APD[0] << " APDoPN[0] " << APDoPN[0] << " APDoAPDA[0] "
             << APDoAPDA[0] << " flag " << flag << endl;
      restrees[iColor]->Fill();
    }
  }

  // Get final results for PN
  //==========================

  for (unsigned int iM = 0; iM < nMod; iM++) {
    unsigned int iMod = modules[iM] - 1;

    side = iside[firstChanMod[iMod]];

    for (unsigned int ch = 0; ch < nPNPerMod; ch++) {
      pnID = ch;
      moduleID = iMod + 1;

      if (moduleID >= 20)
        moduleID -= 2;  // Trick to fix endcap specificity

      for (unsigned int iColor = 0; iColor < nCol; iColor++) {
        std::vector<double> pnvec = PNAnal[iMod][ch][iColor]->getPN();
        std::vector<double> pnopnvec = PNAnal[iMod][ch][iColor]->getPNoPN();
        std::vector<double> pnopn0vec = PNAnal[iMod][ch][iColor]->getPNoPN0();
        std::vector<double> pnopn1vec = PNAnal[iMod][ch][iColor]->getPNoPN1();

        for (unsigned int i = 0; i < pnvec.size(); i++) {
          PN[i] = pnvec.at(i);
          PNoPN[i] = pnopnvec.at(i);
          PNoPNA[i] = pnopn0vec.at(i);
          PNoPNB[i] = pnopn1vec.at(i);
        }

        if (_debug >= 1)
          cout << "-- debug test -- endJob -- filling pn results'tree: PN[0]:" << PN[0] << " iModule:" << iMod
               << " iColor:" << iColor << " ch:" << ch << endl;

        // Fill PN results trees
        //========================

        respntrees[iColor]->Fill();
      }
    }
  }

  // Remove temporary files
  //========================

  if (!_saveallevents) {
    APDFile->Close();
    stringstream del2;
    del2 << "rm " << APDfile;
    system(del2.str().c_str());

  } else {
    APDFile->cd();
    APDtrees[0]->Write();
    APDFile->Close();
    resFile->cd();
  }

  // Save results
  //===============

  for (unsigned int i = 0; i < nCol; i++) {
    restrees[i]->Write();
    respntrees[i]->Write();
  }

  resFile->Close();

  cout << "\t+=+    .................................................. done  +=+" << endl;
  cout << "\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << endl;
}

//========================================================================
bool EcalLaserAnalyzer2::getShapes() {
  //========================================================================

  // Get Pulse From Matacq Analysis:
  //================================

  bool IsMatacqOK = false;

  int doesMatFileExist = 0;
  int doesMatShapeExist = 0;
  FILE* test2;
  TProfile* laserShape = nullptr;
  test2 = fopen(matfile.c_str(), "r");
  if (test2)
    doesMatFileExist = 1;

  TFile* MatShapeFile;
  if (doesMatFileExist == 1) {
    MatShapeFile = new TFile(matfile.c_str());
    laserShape = (TProfile*)MatShapeFile->Get("shapeLaser");
    if (laserShape) {
      doesMatShapeExist = 1;
      double y = laserShape->Integral("w");
      if (y != 0)
        laserShape->Scale(1.0 / y);
    }
  } else {
    cout << " ERROR! Matacq shape file not found !" << endl;
  }
  if (doesMatShapeExist)
    IsMatacqOK = true;

  // Get SPR from the average elec shape in SM6:
  //============================================

  int doesElecFileExist = 0;
  FILE* test;
  test = fopen(elecfile_.c_str(), "r");
  if (test)
    doesElecFileExist = 1;

  TFile* ElecShapesFile;
  TH1D* elecShape = nullptr;

  if (doesElecFileExist == 1) {
    ElecShapesFile = new TFile(elecfile_.c_str());
    stringstream name;
    name << "MeanElecShape";
    elecShape = (TH1D*)ElecShapesFile->Get(name.str().c_str());
    if (elecShape && doesMatShapeExist == 1) {
      double x = elecShape->GetMaximum();
      if (x != 0)
        elecShape->Scale(1.0 / x);
      isSPRFine = true;
    } else {
      isSPRFine = false;
    }

  } else {
    cout << " ERROR! Elec shape file not found !" << endl;
  }

  if (IsMatacqOK) {
    ShapeFile = new TFile(shapefile.c_str(), "RECREATE");

    unsigned int nBins = int(laserShape->GetEntries());
    assert(nSamplesShapes == nBins);
    double elec_jj, laser_iiMinusjj;
    double sum_jj;

    if (isSPRFine == true) {
      unsigned int nBins2 = int(elecShape->GetNbinsX());

      if (nBins2 < nBins) {
        cout << "EcalLaserAnalyzer2::getShapes: wrong configuration of the shapes' number of bins" << std::endl;
        isSPRFine = false;
      }
      assert(nSamplesShapes == nBins2);

      stringstream name;
      name << "PulseShape";

      PulseShape = new TProfile(name.str().c_str(), name.str().c_str(), nBins, -0.5, double(nBins) - 0.5);

      // shift shapes to have max close to the real APD max

      for (int ii = 0; ii < 50; ii++) {
        shapes[ii] = 0.0;
        PulseShape->Fill(ii, 0.0);
      }

      for (unsigned int ii = 0; ii < nBins - 50; ii++) {
        sum_jj = 0.0;
        for (unsigned int jj = 0; jj < ii; jj++) {
          elec_jj = elecShape->GetBinContent(jj + 1);
          laser_iiMinusjj = laserShape->GetBinContent(ii - jj + 1);
          sum_jj += elec_jj * laser_iiMinusjj;
        }
        PulseShape->Fill(ii + 50, sum_jj);
        shapes[ii + 50] = sum_jj;
      }

      double scale = PulseShape->GetMaximum();
      shapeCorrection = scale;

      if (scale != 0) {
        PulseShape->Scale(1.0 / scale);
        for (unsigned int ii = 0; ii < nBins; ii++) {
          shapesVec.push_back(shapes[ii] / scale);
        }
      }

      if (_saveshapes)
        PulseShape->Write();
    }
  }
  ShapeFile->Close();

  if (!_saveshapes) {
    stringstream del;
    del << "rm " << shapefile;
    system(del.str().c_str());
  }

  return IsMatacqOK;
}

void EcalLaserAnalyzer2::setGeomEB(
    int etaG, int phiG, int module, int tower, int strip, int xtal, int apdRefTT, int channel, int lmr) {
  side = MEEBGeom::side(etaG, phiG);

  assert(module >= *min_element(modules.begin(), modules.end()) &&
         module <= *max_element(modules.begin(), modules.end()));

  eta = etaG;
  phi = phiG;
  channelID = 5 * (strip - 1) + xtal - 1;
  towerID = tower;

  vector<int> apdRefChan = ME::apdRefChannels(module, lmr);
  for (unsigned int iref = 0; iref < nRefChan; iref++) {
    if (channelID == apdRefChan[iref] && towerID == apdRefTT && apdRefMap[iref].count(module) == 0) {
      apdRefMap[iref][module] = channel;
    }
  }

  if (isFirstChanModFilled[module - 1] == 0) {
    firstChanMod[module - 1] = channel;
    isFirstChanModFilled[module - 1] = 1;
  }

  iEta[channel] = eta;
  iPhi[channel] = phi;
  iModule[channel] = module;
  iTowerID[channel] = towerID;
  iChannelID[channel] = channelID;
  idccID[channel] = dccID;
  iside[channel] = side;
}

void EcalLaserAnalyzer2::setGeomEE(
    int etaG, int phiG, int iX, int iY, int iZ, int module, int tower, int ch, int apdRefTT, int channel, int lmr) {
  side = MEEEGeom::side(iX, iY, iZ);

  assert(module >= *min_element(modules.begin(), modules.end()) &&
         module <= *max_element(modules.begin(), modules.end()));

  eta = etaG;
  phi = phiG;
  channelID = ch;
  towerID = tower;

  vector<int> apdRefChan = ME::apdRefChannels(module, lmr);
  for (unsigned int iref = 0; iref < nRefChan; iref++) {
    if (channelID == apdRefChan[iref] && towerID == apdRefTT && apdRefMap[iref].count(module) == 0) {
      apdRefMap[iref][module] = channel;
    }
  }

  if (isFirstChanModFilled[module - 1] == 0) {
    firstChanMod[module - 1] = channel;
    isFirstChanModFilled[module - 1] = 1;
  }

  iEta[channel] = eta;
  iPhi[channel] = phi;
  iModule[channel] = module;
  iTowerID[channel] = towerID;
  iChannelID[channel] = channelID;
  idccID[channel] = dccID;
  iside[channel] = side;
}
DEFINE_FWK_MODULE(EcalLaserAnalyzer2);