EcalMatacqAnalyzer.cc

Go to the documentation of this file.

/* 
 *  \class EcalMatacqAnalyzer
 *
 *  primary author: Gautier Hamel De Monchenault - CEA/Saclay
 *  author: Julie Malcles - CEA/Saclay
 */
 
#include <TFile.h>
#include <TTree.h>
#include <TProfile.h>
#include <TChain.h>
#include <vector>
 
#include <CalibCalorimetry/EcalLaserAnalyzer/plugins/EcalMatacqAnalyzer.h>
 
#include <sstream>
#include <iostream>
#include <iomanip>
 
#include <FWCore/MessageLogger/interface/MessageLogger.h>
#include <FWCore/Utilities/interface/Exception.h>
 
#include <DataFormats/EcalDetId/interface/EcalDetIdCollections.h>
#include <DataFormats/EcalRawData/interface/EcalRawDataCollections.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 <CalibCalorimetry/EcalLaserAnalyzer/interface/TMatacq.h>
#include <CalibCalorimetry/EcalLaserAnalyzer/interface/TMom.h>
#include <CalibCalorimetry/EcalLaserAnalyzer/interface/TMTQ.h>
 
using namespace std;
 
//========================================================================
EcalMatacqAnalyzer::EcalMatacqAnalyzer(const edm::ParameterSet& iConfig)
    //========================================================================
    :
 
      iEvent(0),
 
      // framework parameters with default values
 
      _presample(iConfig.getUntrackedParameter<double>("nPresamples", 6.7)),
      _nsamplesaftmax(iConfig.getUntrackedParameter<unsigned int>("nSamplesAftMax", 80)),
      _noiseCut(iConfig.getUntrackedParameter<unsigned int>("noiseCut", 7)),
      _parabnbefmax(iConfig.getUntrackedParameter<unsigned int>("paraBeforeMax", 8)),
      _parabnaftmax(iConfig.getUntrackedParameter<unsigned int>("paraAfterMax", 7)),
      _thres(iConfig.getUntrackedParameter<unsigned int>("threshold", 10)),
      _lowlev(iConfig.getUntrackedParameter<unsigned int>("lowLevel", 20)),
      _highlev(iConfig.getUntrackedParameter<unsigned int>("highLevel", 80)),
      _nevlasers(iConfig.getUntrackedParameter<unsigned int>("nEventLaser", 600)),
      _timebefmax(iConfig.getUntrackedParameter<unsigned int>("timeBefMax", 100)),
      _timeaftmax(iConfig.getUntrackedParameter<unsigned int>("timeAftMax", 250)),
      _cutwindow(iConfig.getUntrackedParameter<double>("cutWindow", 0.1)),
      _nsamplesshape(iConfig.getUntrackedParameter<unsigned int>("nSamplesShape", 250)),
      _presampleshape(iConfig.getUntrackedParameter<unsigned int>("nPresamplesShape", 50)),
      _slide(iConfig.getUntrackedParameter<unsigned int>("nSlide", 100)),
      _fedid(iConfig.getUntrackedParameter<int>("fedID", -999)),
      _debug(iConfig.getUntrackedParameter<int>("debug", 0)),
      nSides(NSIDES),
      lightside(0),
      runType(-1),
      runNum(0),
      event(0),
      color(-1),
      maxsamp(0),
      nsamples(0),
      tt(0)
 
//========================================================================
{
  //now do what ever initialization is needed
 
  resdir_ = iConfig.getUntrackedParameter<std::string>("resDir");
 
  digiCollection_ = iConfig.getParameter<std::string>("digiCollection");
  digiProducer_ = iConfig.getParameter<std::string>("digiProducer");
 
  eventHeaderCollection_ = iConfig.getParameter<std::string>("eventHeaderCollection");
  eventHeaderProducer_ = iConfig.getParameter<std::string>("eventHeaderProducer");
}
 
//========================================================================
EcalMatacqAnalyzer::~EcalMatacqAnalyzer() {
  //========================================================================
 
  // do anything here that needs to be done at desctruction time
  // (e.g. close files, deallocate resources etc.)
}
 
//========================================================================
void EcalMatacqAnalyzer::beginJob() {
  //========================================================================
 
  // Define temporary file name
 
  sampfile = resdir_;
  sampfile += "/TmpTreeMatacqAnalyzer.root";
 
  sampFile = new TFile(sampfile.c_str(), "RECREATE");
 
  // declaration of the tree to fill
 
  tree = new TTree("MatacqTree", "MatacqTree");
 
  //List of branches
 
  tree->Branch("event", &event, "event/I");
  tree->Branch("color", &color, "color/I");
  tree->Branch("matacq", &matacq, "matacq[2560]/D");
  tree->Branch("nsamples", &nsamples, "nsamples/I");
  tree->Branch("maxsamp", &maxsamp, "maxsamp/I");
  tree->Branch("tt", &tt, "tt/D");
  tree->Branch("lightside", &lightside, "lightside/I");
 
  tree->SetBranchAddress("event", &event);
  tree->SetBranchAddress("color", &color);
  tree->SetBranchAddress("matacq", matacq);
  tree->SetBranchAddress("nsamples", &nsamples);
  tree->SetBranchAddress("maxsamp", &maxsamp);
  tree->SetBranchAddress("tt", &tt);
  tree->SetBranchAddress("lightside", &lightside);
 
  // Define output results files' names
 
  std::stringstream namefile;
  namefile << resdir_ << "/MATACQ.root";
  outfile = namefile.str();
 
  // Laser events counter
  laserEvents = 0;
  isThereMatacq = false;
}
 
//========================================================================
void EcalMatacqAnalyzer::analyze(const edm::Event& e, const edm::EventSetup& c) {
  //========================================================================
 
  ++iEvent;
 
  if (_debug == 1)
    std::cout << "-- debug test -- Entering Analyze -- event= " << iEvent << std::endl;
 
  // retrieving MATACQ :
  edm::Handle<EcalMatacqDigiCollection> pmatacqDigi;
  const EcalMatacqDigiCollection* matacqDigi = nullptr;
  try {
    e.getByLabel(digiProducer_, digiCollection_, pmatacqDigi);
    matacqDigi = pmatacqDigi.product();
    if (_debug == 1)
      std::cout << "-- debug test -- Matacq Digis Found -- " << std::endl;
 
  } catch (std::exception& ex) {
    std::cerr << "Error! can't get the product EcalMatacqDigi producer:" << digiProducer_.c_str()
              << " collection:" << digiCollection_.c_str() << std::endl;
    if (_debug == 1)
      std::cout << "-- debug test -- No Matacq Digis Found -- " << std::endl;
    return;
  }
 
  // 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 EcalRawData producer:" << eventHeaderProducer_.c_str()
              << " collection:" << eventHeaderCollection_.c_str() << std::endl;
    return;
  }
 
  // ====================================
  // Decode Basic DCCHeader Information
  // ====================================
 
  if (_debug == 1)
    std::cout << "-- debug test -- Before header -- " << std::endl;
 
  for (EcalRawDataCollection::const_iterator headerItr = DCCHeader->begin(); headerItr != DCCHeader->end();
       ++headerItr) {
    EcalDCCHeaderBlock::EcalDCCEventSettings settings = headerItr->getEventSettings();
    color = (int)settings.wavelength;
    if (color < 0)
      return;
 
    // 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();
 
    if (_debug == 1)
      std::cout << "-- debug test -- runtype:" << runType << " event:" << event << " runNum:" << runNum << std::endl;
 
    dccID = headerItr->getDccInTCCCommand();
    fedID = headerItr->fedId();
    lightside = headerItr->getRtHalf();
 
    //assert (lightside<2 && lightside>=0);
 
    if (lightside != 1 && lightside != 0) {
      std::cout << "Unexpected lightside: " << lightside << " for event " << iEvent << std::endl;
      return;
    }
    if (_debug == 1) {
      std::cout << "-- debug test -- Inside header before fed cut -- color=" << color << ", dcc=" << dccID
                << ", fed=" << fedID << ",  lightside=" << lightside << ", runType=" << runType << std::endl;
    }
 
    // take event only if the fed corresponds to the DCC in TCC
    if (600 + dccID != fedID)
      continue;
 
    if (_debug == 1) {
      std::cout << "-- debug test -- Inside header after fed cut -- color=" << color << ", dcc=" << dccID
                << ", fed=" << fedID << ",  lightside=" << lightside << ", runType=" << runType << std::endl;
    }
 
    // Cut on runType
 
    if (runType != EcalDCCHeaderBlock::LASER_STD && runType != EcalDCCHeaderBlock::LASER_GAP &&
        runType != EcalDCCHeaderBlock::LASER_POWER_SCAN && runType != EcalDCCHeaderBlock::LASER_DELAY_SCAN)
      return;
 
    std::vector<int>::iterator iter = find(colors.begin(), colors.end(), color);
    if (iter == colors.end()) {
      colors.push_back(color);
      std::cout << " new color found " << color << " " << colors.size() << std::endl;
    }
  }
 
  if (_debug == 1)
    std::cout << "-- debug test -- Before digis -- Event:" << iEvent << std::endl;
 
  // Count laser events
  laserEvents++;
 
  // ===========================
  // Decode Matacq Information
  // ===========================
 
  int iCh = 0;
  double max = 0;
 
  for (EcalMatacqDigiCollection::const_iterator it = matacqDigi->begin(); it != matacqDigi->end();
       ++it) {  // Loop on matacq channel
 
    //
    const EcalMatacqDigi& digis = *it;
 
    //if(digis.size()==0 || iCh>=N_channels) continue;
    if (_debug == 1) {
      std::cout << "-- debug test -- Inside digis -- digi size=" << digis.size() << std::endl;
    }
 
    if (digis.size() == 0)
      continue;
    else
      isThereMatacq = true;
 
    max = 0;
    maxsamp = 0;
    nsamples = digis.size();
    tt = digis.tTrig();
 
    for (int i = 0; i < digis.size(); ++i) {  // Loop on matacq samples
      matacq[i] = -digis.adcCount(i);
      if (matacq[i] > max) {
        max = matacq[i];
        maxsamp = i;
      }
    }
    if (_debug == 1) {
      std::cout << "-- debug test -- Inside digis -- nsamples=" << nsamples << ", max=" << max << std::endl;
    }
 
    iCh++;
  }
 
  if (_debug == 1)
    std::cout << "-- debug test -- After digis -- Event: " << iEvent << std::endl;
  tree->Fill();
 
}  // analyze
 
//========================================================================
void EcalMatacqAnalyzer::endJob() {
  // Don't do anything if there is no events
  if (!isThereMatacq) {
    std::cout << "\n\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << std::endl;
    std::cout << "\t+=+     WARNING! NO MATACQ        +=+" << std::endl;
    std::cout << "\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << std::endl;
 
    // Remove temporary file
    FILE* test;
    test = fopen(sampfile.c_str(), "r");
    if (test) {
      std::stringstream del2;
      del2 << "rm " << sampfile;
      system(del2.str().c_str());
    }
    return;
  }
 
  assert(colors.size() <= nColor);
  unsigned int nCol = colors.size();
 
  for (unsigned int iCol = 0; iCol < nCol; iCol++) {
    for (unsigned int iSide = 0; iSide < nSide; iSide++) {
      MTQ[iCol][iSide] = new TMTQ();
    }
  }
 
  outFile = new TFile(outfile.c_str(), "RECREATE");
 
  TProfile* shapeMat = new TProfile("shapeLaser", "shapeLaser", _nsamplesshape, -0.5, double(_nsamplesshape) - 0.5);
  TProfile* shapeMatTmp = new TProfile(
      "shapeLaserTmp", "shapeLaserTmp", _timeaftmax + _timebefmax, -0.5, double(_timeaftmax + _timebefmax) - 0.5);
 
  std::cout << "\n\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << std::endl;
  std::cout << "\t+=+     Analyzing MATACQ data     +=+" << std::endl;
 
  //
  // create output ntuple
  //
 
  mtqShape = new TTree("MatacqShape", "MatacqShape");
 
  // list of branches
  // keep Patrice's notations
 
  mtqShape->Branch("event", &event, "event/I");
  mtqShape->Branch("color", &color, "color/I");
  mtqShape->Branch("status", &status, "status/I");
  mtqShape->Branch("peak", &peak, "peak/D");
  mtqShape->Branch("sigma", &sigma, "sigma/D");
  mtqShape->Branch("fit", &fit, "fit/D");
  mtqShape->Branch("ampl", &ampl, "ampl/D");
  mtqShape->Branch("trise", &trise, "trise/D");
  mtqShape->Branch("fwhm", &fwhm, "fwhm/D");
  mtqShape->Branch("fw20", &fw20, "fw20/D");
  mtqShape->Branch("fw80", &fw80, "fw80/D");
  mtqShape->Branch("ped", &ped, "ped/D");
  mtqShape->Branch("pedsig", &pedsig, "pedsig/D");
  mtqShape->Branch("ttrig", &ttrig, "ttrig/D");
  mtqShape->Branch("sliding", &sliding, "sliding/D");
 
  mtqShape->SetBranchAddress("event", &event);
  mtqShape->SetBranchAddress("color", &color);
  mtqShape->SetBranchAddress("status", &status);
  mtqShape->SetBranchAddress("peak", &peak);
  mtqShape->SetBranchAddress("sigma", &sigma);
  mtqShape->SetBranchAddress("fit", &fit);
  mtqShape->SetBranchAddress("ampl", &ampl);
  mtqShape->SetBranchAddress("fwhm", &fwhm);
  mtqShape->SetBranchAddress("fw20", &fw20);
  mtqShape->SetBranchAddress("fw80", &fw80);
  mtqShape->SetBranchAddress("trise", &trise);
  mtqShape->SetBranchAddress("ped", &ped);
  mtqShape->SetBranchAddress("pedsig", &pedsig);
  mtqShape->SetBranchAddress("ttrig", &ttrig);
  mtqShape->SetBranchAddress("sliding", &sliding);
 
  unsigned int endsample;
  unsigned int presample;
 
  // loop over the entries of the tree
  TChain* fChain = (TChain*)tree;
  Long64_t nentries = fChain->GetEntriesFast();
  Long64_t nbytes = 0, nb = 0;
 
  for (Long64_t jentry = 0; jentry < nentries; jentry++) {
    // load the event
    Long64_t ientry = fChain->LoadTree(jentry);
    if (ientry < 0)
      break;
    nb = fChain->GetEntry(jentry);
    nbytes += nb;
 
    status = 0;
    peak = -1;
    sigma = 0;
    fit = -1;
    ampl = -1;
    trise = -1;
    ttrig = tt;
    fwhm = 0;
    fw20 = 0;
    fw80 = 0;
    ped = 0;
    pedsig = 0;
    sliding = 0;
 
    if (_debug == 1)
      std::cout << "-- debug test -- inside loop 1  -- jentry:" << jentry << " over nentries=" << nentries << std::endl;
 
    // create the object for Matacq data analysis
 
    endsample = maxsamp + _nsamplesaftmax;
    presample = int(_presample * nsamples / 100.);
    TMatacq* mtq = new TMatacq(nsamples,
                               presample,
                               endsample,
                               _noiseCut,
                               _parabnbefmax,
                               _parabnaftmax,
                               _thres,
                               _lowlev,
                               _highlev,
                               _nevlasers,
                               _slide);
 
    if (_debug == 1)
      std::cout << "-- debug test -- inside loop 2  -- " << std::endl;
 
    // analyse the Matacq data
    if (mtq->rawPulseAnalysis(nsamples, &matacq[0]) == 0) {
      status = 1;
      ped = mtq->getBaseLine();
      pedsig = mtq->getsigBaseLine();
 
      if (_debug == 1)
        std::cout << "-- debug test -- inside loop 3  -- ped:" << ped << std::endl;
      if (mtq->findPeak() == 0) {
        peak = mtq->getTimpeak();
        sigma = mtq->getsigTimpeak();
      }
      if (_debug == 1)
        std::cout << "-- debug test -- inside loop 4  -- peak:" << peak << std::endl;
      if (mtq->doFit() == 0) {
        fit = mtq->getTimax();
        ampl = mtq->getAmpl();
        fwhm = mtq->getFwhm();
        fw20 = mtq->getWidth20();
        fw80 = mtq->getWidth80();
        sliding = mtq->getSlide();
      }
      if (_debug == 1)
        std::cout << "-- debug test -- inside loop 4  -- ampl:" << ampl << std::endl;
      if (mtq->compute_trise() == 0) {
        trise = mtq->getTrise();
      }
      if (_debug == 1)
        std::cout << "-- debug test -- inside loop 5  -- trise:" << trise << std::endl;
    }
 
    if (_debug == 1)
      std::cout << "-- debug test -- inside loop 6  -- status:" << status << std::endl;
 
    if (status == 1 && mtq->findPeak() == 0) {
      int firstS = int(peak - double(_timebefmax));
      int lastS = int(peak + double(_timeaftmax));
 
      // Fill histo if there are enough samples
      if (_debug == 1)
        std::cout << "-- debug test -- inside loop 7  -- firstS:" << firstS << ", nsamples:" << nsamples << std::endl;
 
      if (firstS >= 0 && lastS <= nsamples) {
        for (int i = firstS; i < lastS; i++) {
          shapeMatTmp->Fill(double(i) - firstS, matacq[i]);
        }
 
      } else {  // else extrapolate
 
        int firstSBis;
 
        if (firstS < 0) {  // fill first bins with 0
 
          double thisped;
          thisped = (matacq[0] + matacq[1] + matacq[2] + matacq[4] + matacq[5]) / 5.0;
 
          for (int i = firstS; i < 0; i++) {
            shapeMatTmp->Fill(double(i) - firstS, thisped);
          }
          firstSBis = 0;
 
        } else {
          firstSBis = firstS;
        }
 
        if (lastS > nsamples) {
          for (int i = firstSBis; i < int(nsamples); i++) {
            shapeMatTmp->Fill(double(i) - firstS, matacq[i]);
          }
 
          //extrapolate with expo tail
 
          double expb = 0.998;
          double matacqval = expb * matacq[nsamples - 1];
 
          for (int i = nsamples; i < lastS; i++) {
            shapeMatTmp->Fill(double(i) - firstS, matacqval);
            matacqval *= expb;
          }
 
        } else {
          for (int i = firstSBis; i < lastS; i++) {
            shapeMatTmp->Fill(double(i) - firstS, matacq[i]);
          }
        }
      }
    }
    if (_debug == 1)
      std::cout << "-- debug test -- inside loop 8" << std::endl;
 
    // get back color
 
    int iCol = nCol;
    for (unsigned int i = 0; i < nCol; i++) {
      if (color == colors[i]) {
        iCol = i;
        i = nCol;
      }
    }
    if (_debug == 1)
      std::cout << "-- debug test -- inside loop 8bis color:" << color << " iCol:" << iCol << " nCol:" << nCol
                << std::endl;
 
    // fill TMTQ
 
    if (status == 1 && mtq->findPeak() == 0 && mtq->doFit() == 0 && mtq->compute_trise() == 0)
      MTQ[iCol][lightside]->addEntry(peak, sigma, fit, ampl, trise, fwhm, fw20, fw80, ped, pedsig, sliding);
 
    // fill the output tree
 
    if (_debug == 1)
      std::cout << "-- debug test -- inside loop 9" << std::endl;
    mtqShape->Fill();
 
    // clean up
    delete mtq;
  }
 
  if (_debug == 1)
    std::cout << "-- debug test -- after loop " << std::endl;
  sampFile->Close();
 
  double Peak[6], Sigma[6], Fit[6], Ampl[6], Trise[6], Fwhm[6], Fw20[6], Fw80[6], Ped[6], Pedsig[6], Sliding[6];
  int Side;
 
  for (unsigned int iColor = 0; iColor < nCol; iColor++) {
    std::stringstream nametree;
    nametree << "MatacqCol" << colors[iColor];
    meanTree[iColor] = new TTree(nametree.str().c_str(), nametree.str().c_str());
    meanTree[iColor]->Branch("side", &Side, "Side/I");
    meanTree[iColor]->Branch("peak", &Peak, "Peak[6]/D");
    meanTree[iColor]->Branch("sigma", &Sigma, "Sigma[6]/D");
    meanTree[iColor]->Branch("fit", &Fit, "Fit[6]/D");
    meanTree[iColor]->Branch("ampl", &Ampl, "Ampl[6]/D");
    meanTree[iColor]->Branch("trise", &Trise, "Trise[6]/D");
    meanTree[iColor]->Branch("fwhm", &Fwhm, "Fwhm[6]/D");
    meanTree[iColor]->Branch("fw20", &Fw20, "Fw20[6]/D");
    meanTree[iColor]->Branch("fw80", &Fw80, "Fw80[6]/D");
    meanTree[iColor]->Branch("ped", &Ped, "Ped[6]/D");
    meanTree[iColor]->Branch("pedsig", &Pedsig, "Pedsig[6]/D");
    meanTree[iColor]->Branch("sliding", &Sliding, "Sliding[6]/D");
 
    meanTree[iColor]->SetBranchAddress("side", &Side);
    meanTree[iColor]->SetBranchAddress("peak", Peak);
    meanTree[iColor]->SetBranchAddress("sigma", Sigma);
    meanTree[iColor]->SetBranchAddress("fit", Fit);
    meanTree[iColor]->SetBranchAddress("ampl", Ampl);
    meanTree[iColor]->SetBranchAddress("fwhm", Fwhm);
    meanTree[iColor]->SetBranchAddress("fw20", Fw20);
    meanTree[iColor]->SetBranchAddress("fw80", Fw80);
    meanTree[iColor]->SetBranchAddress("trise", Trise);
    meanTree[iColor]->SetBranchAddress("ped", Ped);
    meanTree[iColor]->SetBranchAddress("pedsig", Pedsig);
    meanTree[iColor]->SetBranchAddress("sliding", Sliding);
  }
 
  for (unsigned int iCol = 0; iCol < nCol; iCol++) {
    for (unsigned int iSide = 0; iSide < nSides; iSide++) {
      Side = iSide;
      std::vector<double> val[TMTQ::nOutVar];
 
      for (int iVar = 0; iVar < TMTQ::nOutVar; iVar++) {
        val[iVar] = MTQ[iCol][iSide]->get(iVar);
 
        for (unsigned int i = 0; i < val[iVar].size(); i++) {
          switch (iVar) {
            case TMTQ::iPeak:
              Peak[i] = val[iVar].at(i);
              break;
            case TMTQ::iSigma:
              Sigma[i] = val[iVar].at(i);
              break;
            case TMTQ::iFit:
              Fit[i] = val[iVar].at(i);
              break;
            case TMTQ::iAmpl:
              Ampl[i] = val[iVar].at(i);
              break;
            case TMTQ::iFwhm:
              Fwhm[i] = val[iVar].at(i);
              break;
            case TMTQ::iFw20:
              Fw20[i] = val[iVar].at(i);
              break;
            case TMTQ::iFw80:
              Fw80[i] = val[iVar].at(i);
              break;
            case TMTQ::iTrise:
              Trise[i] = val[iVar].at(i);
              break;
            case TMTQ::iPed:
              Ped[i] = val[iVar].at(i);
              break;
            case TMTQ::iPedsig:
              Pedsig[i] = val[iVar].at(i);
              break;
            case TMTQ::iSlide:
              Sliding[i] = val[iVar].at(i);
              break;
          }
        }
      }
      meanTree[iCol]->Fill();
      if (_debug == 1)
        std::cout << "-- debug test -- inside final loop  " << std::endl;
    }
  }
 
  // Calculate maximum with pol 2
 
  int im = shapeMatTmp->GetMaximumBin();
  double q1 = shapeMatTmp->GetBinContent(im - 1);
  double q2 = shapeMatTmp->GetBinContent(im);
  double q3 = shapeMatTmp->GetBinContent(im + 1);
 
  double a2 = (q3 + q1) / 2.0 - q2;
  double a1 = q2 - q1 + a2 * (1 - 2 * im);
  double a0 = q2 - a1 * im - a2 * im * im;
 
  double am = a0 - a1 * a1 / (4 * a2);
 
  // Compute pedestal
 
  double bl = 0;
  for (unsigned int i = 1; i < _presampleshape + 1; i++) {
    bl += shapeMatTmp->GetBinContent(i);
  }
  bl /= _presampleshape;
 
  // Compute and save laser shape
 
  if (_debug == 1)
    std::cout << "-- debug test -- computing shape  " << std::endl;
 
  int firstBin = 0;
  double height = 0.0;
 
  for (unsigned int i = _timebefmax; i > _presampleshape; i--) {
    height = shapeMatTmp->GetBinContent(i) - bl;
 
    if (height < (am - bl) * _cutwindow) {
      firstBin = i;
      i = _presampleshape;
    }
  }
 
  unsigned int lastBin = firstBin + _nsamplesshape;
 
  for (unsigned int i = firstBin; i < lastBin; i++) {
    shapeMat->Fill(i - firstBin, shapeMatTmp->GetBinContent(i) - bl);
  }
 
  mtqShape->Write();
  for (unsigned int iColor = 0; iColor < nCol; iColor++) {
    meanTree[iColor]->Write();
  }
  if (_debug == 1)
    std::cout << "-- debug test -- writing  " << std::endl;
  shapeMat->Write();
 
  // close the output file
  outFile->Close();
 
  // Remove temporary file
  FILE* test;
  test = fopen(sampfile.c_str(), "r");
  if (test) {
    std::stringstream del2;
    del2 << "rm " << sampfile;
    system(del2.str().c_str());
  }
 
  std::cout << "\t+=+    .................... done  +=+" << std::endl;
  std::cout << "\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << std::endl;
}
 
DEFINE_FWK_MODULE(EcalMatacqAnalyzer);