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=0;
  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=0;
  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);
      case TMTQ::iSigma: Sigma[i]=val[iVar].at(i);
      case TMTQ::iFit: Fit[i]=val[iVar].at(i);
      case TMTQ::iAmpl: Ampl[i]=val[iVar].at(i);
      case TMTQ::iFwhm: Fwhm[i]=val[iVar].at(i);
      case TMTQ::iFw20: Fw20[i]=val[iVar].at(i);
      case TMTQ::iFw80: Fw80[i]=val[iVar].at(i);
      case TMTQ::iTrise: Trise[i]=val[iVar].at(i);
      case TMTQ::iPed: Ped[i]=val[iVar].at(i);
      case TMTQ::iPedsig: Pedsig[i]=val[iVar].at(i);
      case TMTQ::iSlide: Sliding[i]=val[iVar].at(i);
      }
    }
      }
      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);