#include <EcalMatacqAnalyzer.h>

Inheritance diagram for EcalMatacqAnalyzer:

Public Types
enum	VarCol { iBlue, iRed, nColor }
enum	VarSide { iSide0, iSide1, nSide }
Public Member Functions
virtual void	analyze (const edm::Event &e, const edm::EventSetup &c)
virtual void	beginJob ()
	EcalMatacqAnalyzer (const edm::ParameterSet &iConfig)
virtual void	endJob ()
	~EcalMatacqAnalyzer ()
Private Attributes
double	_cutwindow
int	_debug
int	_fedid
unsigned int	_highlev
unsigned int	_lowlev
unsigned int	_nevlasers
unsigned int	_noiseCut
unsigned int	_nsamplesaftmax
unsigned int	_nsamplesshape
unsigned int	_parabnaftmax
unsigned int	_parabnbefmax
double	_presample
unsigned int	_presampleshape
unsigned int	_slide
unsigned int	_thres
unsigned int	_timeaftmax
unsigned int	_timebefmax
double	ampl
int	color
std::vector< int >	colors
int	dccID
std::string	digiCollection_
std::string	digiProducer_
int	event
std::string	eventHeaderCollection_
std::string	eventHeaderProducer_
int	fedID
double	fit
double	fw20
double	fw80
double	fwhm
int	iEvent
bool	isThereMatacq
int	laserEvents
int	lightside
double	matacq [N_samples]
int	maxsamp
TTree *	meanTree [nColor]
TMTQ *	MTQ [nColor][nSide]
TTree *	mtqShape
int	nsamples
unsigned int	nSides
TFile *	outFile
std::string	outfile
double	peak
double	ped
double	pedsig
std::string	resdir_
int	runNum
int	runType
std::string	sampfile
TFile *	sampFile
double	sigma
double	sliding
int	status
TTree *	tree
double	trise
double	tt
double	ttrig

Detailed Description

Definition at line 14 of file EcalMatacqAnalyzer.h.

Member Enumeration Documentation

enum EcalMatacqAnalyzer::VarCol

Enumerator:

iBlue
iRed
nColor

Definition at line 26 of file EcalMatacqAnalyzer.h.

{ iBlue, iRed, nColor };

enum EcalMatacqAnalyzer::VarSide

Enumerator:

iSide0
iSide1
nSide

Definition at line 27 of file EcalMatacqAnalyzer.h.

{ iSide0, iSide1, nSide };

Constructor & Destructor Documentation

EcalMatacqAnalyzer::EcalMatacqAnalyzer ( const edm::ParameterSet & iConfig ) [explicit]

Definition at line 40 of file EcalMatacqAnalyzer.cc.

References digiCollection_, digiProducer_, eventHeaderCollection_, eventHeaderProducer_, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), and resdir_.

  :

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 ( )

Definition at line 85 of file EcalMatacqAnalyzer.cc.

                                       {
//========================================================================

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

}

Member Function Documentation

void EcalMatacqAnalyzer::analyze	(	const edm::Event &	e,
		const edm::EventSetup &	c
	)		`[virtual]`

Implements edm::EDAnalyzer.

Definition at line 147 of file EcalMatacqAnalyzer.cc.

References _debug, _fedid, EcalMatacqDigi::adcCount(), edm::SortedCollection< T, SORT >::begin(), benchmark_cfg::cerr, color, colors, gather_cfg::cout, dccID, digiCollection_, digiProducer_, edm::SortedCollection< T, SORT >::end(), eventHeaderCollection_, eventHeaderProducer_, exception, fedID, spr::find(), edm::Event::getByLabel(), i, iEvent, isThereMatacq, EcalDCCHeaderBlock::LASER_DELAY_SCAN, EcalDCCHeaderBlock::LASER_GAP, EcalDCCHeaderBlock::LASER_POWER_SCAN, EcalDCCHeaderBlock::LASER_STD, laserEvents, lightside, matacq, max(), maxsamp, nsamples, edm::Handle< T >::product(), runNum, runType, EcalMatacqDigi::size(), tree, tt, EcalMatacqDigi::tTrig(), and EcalDCCHeaderBlock::EcalDCCEventSettings::wavelength.

                                                                             {
//========================================================================

  ++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::beginJob ( void ) [virtual]

Reimplemented from edm::EDAnalyzer.

Definition at line 96 of file EcalMatacqAnalyzer.cc.

References color, event, isThereMatacq, laserEvents, lightside, matacq, maxsamp, prof2calltree::namefile, nsamples, outfile, resdir_, sampfile, sampFile, tree, and tt.

                                  {
//========================================================================

// 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::endJob ( void ) [virtual]

Reimplemented from edm::EDAnalyzer.

Definition at line 292 of file EcalMatacqAnalyzer.cc.

References _cutwindow, _debug, _highlev, _lowlev, _nevlasers, _noiseCut, _nsamplesaftmax, _nsamplesshape, _parabnaftmax, _parabnbefmax, _presample, _presampleshape, _slide, _thres, _timeaftmax, _timebefmax, revisionDML::addEntry(), ampl, color, colors, TMatacq::compute_trise(), gather_cfg::cout, TMatacq::doFit(), event, TMatacq::findPeak(), fit, fw20, fw80, fwhm, TMTQ::get(), TMatacq::getAmpl(), TMatacq::getBaseLine(), TMatacq::getFwhm(), TMatacq::getsigBaseLine(), TMatacq::getsigTimpeak(), TMatacq::getSlide(), TMatacq::getTimax(), TMatacq::getTimpeak(), TMatacq::getTrise(), TMatacq::getWidth20(), TMatacq::getWidth80(), i, TMTQ::iAmpl, TMTQ::iFit, TMTQ::iFw20, TMTQ::iFw80, TMTQ::iFwhm, TMTQ::iPeak, TMTQ::iPed, TMTQ::iPedsig, TMTQ::iSigma, TMTQ::iSlide, isThereMatacq, TMTQ::iTrise, lightside, matacq, maxsamp, meanTree, MTQ, mtqShape, nColor, TMTQ::nOutVar, nsamples, nSide, nSides, outfile, outFile, peak, ped, pedsig, q1, q2, TMatacq::rawPulseAnalysis(), sampfile, sampFile, FitTarget::Sigma, sigma, sliding, status, MultipleCompare::test, tree, trise, tt, and ttrig.

{
  
  // 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 tm;
  if(a2!=0) tm=-a1/(2.0*a2);
  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;
 
}