Reference Manual

Go to the documentation of this file.
00001 /* 
00002  *  \class EcalMatacqAnalyzer
00003  *
00004  *  $Date: 2012/02/09 10:07:37 $
00005  *  primary author: Gautier Hamel De Monchenault - CEA/Saclay
00006  *  author: Julie Malcles - CEA/Saclay
00007  */
00008 
00009 
00010 #include <TFile.h>
00011 #include <TTree.h>
00012 #include <TProfile.h>
00013 #include <TChain.h>
00014 #include <vector>
00015 
00016 #include <CalibCalorimetry/EcalLaserAnalyzer/plugins/EcalMatacqAnalyzer.h>
00017 
00018 #include <sstream>
00019 #include <iostream>
00020 #include <iomanip>
00021 
00022 #include <FWCore/MessageLogger/interface/MessageLogger.h>
00023 #include <FWCore/Utilities/interface/Exception.h>
00024 
00025 #include <DataFormats/EcalDetId/interface/EcalDetIdCollections.h>
00026 #include <DataFormats/EcalRawData/interface/EcalRawDataCollections.h>
00027 
00028 #include <FWCore/Framework/interface/Event.h>
00029 #include <FWCore/Framework/interface/MakerMacros.h>
00030 #include <FWCore/ParameterSet/interface/ParameterSet.h>
00031 #include <DataFormats/EcalDigi/interface/EcalDigiCollections.h>
00032 
00033 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TMatacq.h>
00034 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TMom.h>
00035 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TMTQ.h>
00036 
00037 using namespace std;
00038 
00039 //========================================================================
00040 EcalMatacqAnalyzer::EcalMatacqAnalyzer(const edm::ParameterSet& iConfig)
00041  //========================================================================
00042   :
00043 
00044 iEvent(0),
00045 
00046 // framework parameters with default values
00047 
00048 _presample(     iConfig.getUntrackedParameter< double       >( "nPresamples",    6.7 ) ),
00049 _nsamplesaftmax(iConfig.getUntrackedParameter< unsigned int >( "nSamplesAftMax",  80 ) ),
00050 _noiseCut(      iConfig.getUntrackedParameter< unsigned int >( "noiseCut",         7 ) ),
00051 _parabnbefmax(  iConfig.getUntrackedParameter< unsigned int >( "paraBeforeMax",    8 ) ),
00052 _parabnaftmax(  iConfig.getUntrackedParameter< unsigned int >( "paraAfterMax",     7 ) ),
00053 _thres(         iConfig.getUntrackedParameter< unsigned int >( "threshold",       10 ) ),
00054 _lowlev(        iConfig.getUntrackedParameter< unsigned int >( "lowLevel",        20 ) ),
00055 _highlev(       iConfig.getUntrackedParameter< unsigned int >( "highLevel",       80 ) ),
00056 _nevlasers(     iConfig.getUntrackedParameter< unsigned int >( "nEventLaser",    600 ) ),
00057 _timebefmax(    iConfig.getUntrackedParameter< unsigned int >( "timeBefMax",     100 ) ),
00058 _timeaftmax(    iConfig.getUntrackedParameter< unsigned int >( "timeAftMax",     250 ) ),
00059 _cutwindow(     iConfig.getUntrackedParameter< double       >( "cutWindow",      0.1 ) ),
00060 _nsamplesshape( iConfig.getUntrackedParameter< unsigned int >( "nSamplesShape",  250 ) ),
00061 _presampleshape(iConfig.getUntrackedParameter< unsigned int >( "nPresamplesShape",50 ) ),
00062 _slide(         iConfig.getUntrackedParameter< unsigned int >( "nSlide",         100 ) ),
00063 _fedid(         iConfig.getUntrackedParameter< int          >( "fedID",         -999 ) ),
00064 _debug(         iConfig.getUntrackedParameter< int          >( "debug",           0  ) ),             
00065 nSides(NSIDES), lightside(0), runType(-1), runNum(0), 
00066 event(0), color(-1), maxsamp(0), nsamples(0), tt(0)
00067   
00068   //========================================================================
00069 {
00070   
00071 
00072   //now do what ever initialization is needed
00073   
00074   resdir_                  = iConfig.getUntrackedParameter<std::string>("resDir");
00075   
00076   digiCollection_          = iConfig.getParameter<std::string>("digiCollection");
00077   digiProducer_            = iConfig.getParameter<std::string>("digiProducer");
00078   
00079   eventHeaderCollection_   = iConfig.getParameter<std::string>("eventHeaderCollection");
00080   eventHeaderProducer_     = iConfig.getParameter<std::string>("eventHeaderProducer");
00081   
00082 }
00083 
00084 //========================================================================
00085 EcalMatacqAnalyzer::~EcalMatacqAnalyzer(){
00086 //========================================================================
00087 
00088 // do anything here that needs to be done at desctruction time
00089 // (e.g. close files, deallocate resources etc.)
00090 
00091 }
00092 
00093 
00094 
00095 //========================================================================
00096 void EcalMatacqAnalyzer::beginJob() {
00097 //========================================================================
00098 
00099 // Define temporary file name
00100 
00101   sampfile=resdir_;
00102   sampfile+="/TmpTreeMatacqAnalyzer.root";
00103   
00104   sampFile = new TFile(sampfile.c_str(),"RECREATE");
00105 
00106 
00107   // declaration of the tree to fill
00108   
00109   tree = new TTree("MatacqTree","MatacqTree");
00110 
00111 
00112     //List of branches
00113 
00114     tree->Branch( "event",       &event,        "event/I"        );
00115     tree->Branch( "color",       &color ,       "color/I"        );
00116     tree->Branch( "matacq",      &matacq ,      "matacq[2560]/D" );
00117     tree->Branch( "nsamples",    &nsamples ,    "nsamples/I"     );
00118     tree->Branch( "maxsamp",     &maxsamp ,     "maxsamp/I"      );
00119     tree->Branch( "tt",          &tt ,          "tt/D"           );
00120     tree->Branch( "lightside",   &lightside ,   "lightside/I"    );
00121     
00122     tree->SetBranchAddress( "event",       &event       );
00123     tree->SetBranchAddress( "color",       &color       );
00124     tree->SetBranchAddress( "matacq",      matacq       ); 
00125     tree->SetBranchAddress( "nsamples",    &nsamples    );
00126     tree->SetBranchAddress( "maxsamp",     &maxsamp     );
00127     tree->SetBranchAddress( "tt",          &tt          );
00128     tree->SetBranchAddress( "lightside",   &lightside   );
00129 
00130 
00131     
00132     // Define output results files' names
00133     
00134     std::stringstream namefile;
00135     namefile << resdir_ <<"/MATACQ.root";      
00136     outfile=namefile.str();
00137     
00138     
00139     // Laser events counter
00140     laserEvents=0;
00141     isThereMatacq=false;      
00142 
00143 }
00144 
00145 
00146 //========================================================================
00147 void EcalMatacqAnalyzer:: analyze( const edm::Event & e, const  edm::EventSetup& c){
00148 //========================================================================
00149 
00150   ++iEvent;
00151   
00152   if (_debug==1 )std::cout << "-- debug test -- Entering Analyze -- event= "<<iEvent<< std::endl; 
00153 
00154   // retrieving MATACQ :
00155   edm::Handle<EcalMatacqDigiCollection> pmatacqDigi;
00156   const EcalMatacqDigiCollection* matacqDigi=0;
00157   try {
00158     e.getByLabel(digiProducer_,digiCollection_, pmatacqDigi); 
00159     matacqDigi=pmatacqDigi.product();
00160     if (_debug==1 )std::cout << "-- debug test -- Matacq Digis Found -- "<< std::endl; 
00161     
00162   }catch ( std::exception& ex ) {
00163     std::cerr << "Error! can't get the product EcalMatacqDigi producer:" << digiProducer_.c_str()<<" collection:"<<digiCollection_.c_str() << std::endl;
00164     if (_debug==1 )std::cout << "-- debug test -- No Matacq Digis Found -- "<< std::endl; 
00165     return;
00166   }
00167   
00168   // retrieving DCC header
00169 
00170   edm::Handle<EcalRawDataCollection> pDCCHeader;
00171   const  EcalRawDataCollection* DCCHeader=0;
00172   try {
00173     e.getByLabel(eventHeaderProducer_,eventHeaderCollection_, pDCCHeader);
00174     DCCHeader=pDCCHeader.product();
00175   }catch ( std::exception& ex ) {
00176     std::cerr << "Error! can't get the product EcalRawData producer:" << eventHeaderProducer_.c_str()<<" collection:"<<eventHeaderCollection_.c_str() << std::endl;
00177     return;
00178   }
00179 
00180   // ====================================
00181   // Decode Basic DCCHeader Information 
00182   // ====================================
00183   
00184   if (_debug==1)  std::cout <<"-- debug test -- Before header -- "<< std::endl; 
00185   
00186   for ( EcalRawDataCollection::const_iterator headerItr= DCCHeader->begin();headerItr != DCCHeader->end(); 
00187         ++headerItr ) {
00188     
00189     EcalDCCHeaderBlock::EcalDCCEventSettings settings = headerItr->getEventSettings(); 
00190     color = (int) settings.wavelength;
00191     if( color<0 ) return;
00192 
00193     // Get run type and run number 
00194 
00195     int fed = headerItr->fedId();  
00196     
00197     if(fed!=_fedid && _fedid!=-999) continue; 
00198     
00199     runType=headerItr->getRunType();
00200     runNum=headerItr->getRunNumber();
00201     event=headerItr->getLV1();
00202 
00203     if (_debug==1)  std::cout <<"-- debug test -- runtype:"<<runType<<" event:"<<event <<" runNum:"<<runNum <<std::endl; 
00204 
00205     dccID=headerItr->getDccInTCCCommand();
00206     fedID=headerItr->fedId();  
00207     lightside=headerItr->getRtHalf();
00208 
00209     //assert (lightside<2 && lightside>=0);
00210 
00211     if( lightside!=1 && lightside!=0 ) {
00212       std::cout << "Unexpected lightside: "<< lightside<<" for event "<<iEvent << std::endl;
00213       return;
00214     }
00215     if (_debug==1) {
00216       std::cout <<"-- debug test -- Inside header before fed cut -- color="<<color<< ", dcc="<<dccID<<", fed="<< fedID<<",  lightside="<<lightside<<", runType="<<runType<< std::endl;
00217     }
00218 
00219     // take event only if the fed corresponds to the DCC in TCC
00220     if( 600+dccID != fedID ) continue;
00221     
00222     if (_debug==1) {
00223       std::cout <<"-- debug test -- Inside header after fed cut -- color="<<color<< ", dcc="<<dccID<<", fed="<< fedID<<",  lightside="<<lightside<<", runType="<<runType<< std::endl;
00224     }
00225 
00226     // Cut on runType
00227 
00228     if ( runType!=EcalDCCHeaderBlock::LASER_STD && runType!=EcalDCCHeaderBlock::LASER_GAP && runType!=EcalDCCHeaderBlock::LASER_POWER_SCAN && runType!=EcalDCCHeaderBlock::LASER_DELAY_SCAN ) return; 
00229     
00230     
00231     std::vector<int>::iterator iter= find( colors.begin(), colors.end(), color );
00232     if( iter==colors.end() ){
00233       colors.push_back( color );
00234       std::cout <<" new color found "<< color<<" "<< colors.size()<< std::endl;
00235     }
00236     
00237   }
00238   
00239   if (_debug==1) std::cout <<"-- debug test -- Before digis -- Event:"<<iEvent<< std::endl; 
00240 
00241   // Count laser events
00242   laserEvents++;
00243   
00244 
00245 
00246 // ===========================
00247 // Decode Matacq Information
00248 // ===========================
00249 
00250   int iCh=0;
00251   double max=0;
00252 
00253   for(EcalMatacqDigiCollection::const_iterator it = matacqDigi->begin(); it!=matacqDigi->end(); ++it){ // Loop on matacq channel 
00254     
00255     // 
00256     const EcalMatacqDigi& digis = *it;
00257     
00258     //if(digis.size()==0 || iCh>=N_channels) continue; 
00259     if (_debug==1) {
00260       std::cout <<"-- debug test -- Inside digis -- digi size="<<digis.size()<< std::endl;
00261     }
00262 
00263     if(digis.size()== 0 ) continue; 
00264     else isThereMatacq=true;
00265 
00266     max=0;
00267     maxsamp=0;
00268     nsamples=digis.size();
00269     tt=digis.tTrig();
00270 
00271     for(int i=0; i<digis.size(); ++i){ // Loop on matacq samples      
00272       matacq[i]=-digis.adcCount(i);
00273       if(matacq[i]>max) {
00274         max=matacq[i];
00275         maxsamp=i;
00276       }
00277     }  
00278     if (_debug==1) {
00279       std::cout <<"-- debug test -- Inside digis -- nsamples="<<nsamples<< ", max="<<max<< std::endl;
00280     }
00281     
00282     iCh++; 
00283   }
00284   
00285   if (_debug==1)std::cout <<"-- debug test -- After digis -- Event: "<<iEvent<< std::endl; 
00286   tree->Fill();
00287   
00288 } // analyze
00289 
00290 
00291 //========================================================================
00292 void EcalMatacqAnalyzer::endJob() 
00293 {
00294   
00295   // Don't do anything if there is no events
00296   if( !isThereMatacq ) {
00297     
00298     std::cout << "\n\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << std::endl;
00299     std::cout <<   "\t+=+     WARNING! NO MATACQ        +=+" << std::endl;
00300     std::cout <<   "\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << std::endl;
00301     
00302 
00303     // Remove temporary file    
00304     FILE *test; 
00305     test = fopen(sampfile.c_str(),"r");
00306     if (test){
00307       std::stringstream del2;
00308       del2 << "rm " <<sampfile;
00309       system(del2.str().c_str()); 
00310     }
00311     return;
00312   }
00313   
00314   assert( colors.size()<= nColor );
00315   unsigned int nCol=colors.size();
00316 
00317   for(unsigned int iCol=0;iCol<nCol;iCol++){
00318     for(unsigned int iSide=0;iSide<nSide;iSide++){
00319       MTQ[iCol][iSide]=new TMTQ();
00320     }
00321   }
00322   
00323   outFile = new TFile(outfile.c_str(),"RECREATE");
00324     
00325   TProfile *shapeMat = new TProfile("shapeLaser","shapeLaser",_nsamplesshape,-0.5,double(_nsamplesshape)-0.5);
00326   TProfile *shapeMatTmp = new TProfile("shapeLaserTmp","shapeLaserTmp",_timeaftmax+_timebefmax,-0.5,double(_timeaftmax+_timebefmax)-0.5);
00327 
00328   std::cout << "\n\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << std::endl;
00329   std::cout <<   "\t+=+     Analyzing MATACQ data     +=+" << std::endl;
00330   
00331   
00332   //
00333   // create output ntuple
00334   //
00335 
00336   mtqShape = new TTree("MatacqShape","MatacqShape");
00337 
00338   // list of branches 
00339   // keep Patrice's notations
00340 
00341   mtqShape->Branch( "event",       &event,       "event/I"       );
00342   mtqShape->Branch( "color",       &color,       "color/I"       );
00343   mtqShape->Branch( "status",      &status,      "status/I"      );
00344   mtqShape->Branch( "peak",        &peak ,       "peak/D"        );
00345   mtqShape->Branch( "sigma",       &sigma ,      "sigma/D"       );
00346   mtqShape->Branch( "fit",         &fit ,        "fit/D"         );
00347   mtqShape->Branch( "ampl",        &ampl ,       "ampl/D"        );
00348   mtqShape->Branch( "trise",       &trise ,      "trise/D"       );  
00349   mtqShape->Branch( "fwhm",        &fwhm ,       "fwhm/D"        );    
00350   mtqShape->Branch( "fw20",        &fw20 ,       "fw20/D"        );  
00351   mtqShape->Branch( "fw80",        &fw80 ,       "fw80/D"        );  
00352   mtqShape->Branch( "ped",         &ped ,        "ped/D"         );   
00353   mtqShape->Branch( "pedsig",      &pedsig ,     "pedsig/D"      );   
00354   mtqShape->Branch( "ttrig",       &ttrig ,      "ttrig/D"       );  
00355   mtqShape->Branch( "sliding",     &sliding ,    "sliding/D"     );  
00356 
00357   mtqShape->SetBranchAddress( "event",       &event       );
00358   mtqShape->SetBranchAddress( "color",       &color       );
00359   mtqShape->SetBranchAddress( "status",      &status      ); 
00360   mtqShape->SetBranchAddress( "peak",        &peak        ); 
00361   mtqShape->SetBranchAddress( "sigma",       &sigma       ); 
00362   mtqShape->SetBranchAddress( "fit",         &fit         ); 
00363   mtqShape->SetBranchAddress( "ampl",        &ampl        ); 
00364   mtqShape->SetBranchAddress( "fwhm",        &fwhm        ); 
00365   mtqShape->SetBranchAddress( "fw20",        &fw20        ); 
00366   mtqShape->SetBranchAddress( "fw80",        &fw80        ); 
00367   mtqShape->SetBranchAddress( "trise",       &trise       ); 
00368   mtqShape->SetBranchAddress( "ped",         &ped         ); 
00369   mtqShape->SetBranchAddress( "pedsig",      &pedsig      ); 
00370   mtqShape->SetBranchAddress( "ttrig",       &ttrig       ); 
00371   mtqShape->SetBranchAddress( "sliding",     &sliding     ); 
00372 
00373   
00374   unsigned int endsample;
00375   unsigned int presample;
00376   
00377   
00378   
00379   // loop over the entries of the tree
00380   TChain* fChain = (TChain*)tree;
00381   Long64_t nentries = fChain->GetEntriesFast();
00382   Long64_t nbytes = 0, nb = 0;
00383   
00384   for( Long64_t jentry=0; jentry<nentries; jentry++ ) 
00385     {
00386       // load the event
00387       Long64_t ientry = fChain->LoadTree(jentry);
00388       if (ientry < 0) break;
00389       nb = fChain->GetEntry( jentry );   nbytes += nb;      
00390 
00391       status   = 0;
00392       peak    = -1;
00393       sigma   =  0;
00394       fit     = -1;
00395       ampl    = -1;
00396       trise   = -1; 
00397       ttrig   = tt; 
00398       fwhm    =  0;
00399       fw20    =  0;
00400       fw80    =  0;
00401       ped     =  0;
00402       pedsig  =  0;
00403       sliding =  0;
00404 
00405 
00406       if (_debug==1)std::cout <<"-- debug test -- inside loop 1  -- jentry:"<<jentry<<" over nentries="<<nentries<< std::endl; 
00407 
00408       // create the object for Matacq data analysis
00409 
00410       endsample = maxsamp+_nsamplesaftmax;
00411       presample=int(_presample*nsamples/100.);
00412       TMatacq* mtq = new TMatacq( nsamples, presample, endsample,
00413                                   _noiseCut, _parabnbefmax, _parabnaftmax,
00414                                   _thres, _lowlev, _highlev,
00415                                   _nevlasers , _slide);
00416       
00417       if (_debug==1)std::cout <<"-- debug test -- inside loop 2  -- "<< std::endl; 
00418 
00419       // analyse the Matacq data
00420       if( mtq->rawPulseAnalysis( nsamples, &matacq[0] )==0 ) 
00421         {
00422           status = 1;
00423           ped =  mtq->getBaseLine();
00424           pedsig =  mtq->getsigBaseLine();
00425 
00426           if (_debug==1)std::cout <<"-- debug test -- inside loop 3  -- ped:"<<ped << std::endl; 
00427           if( mtq->findPeak()==0 ) 
00428             {
00429               peak = mtq->getTimpeak();
00430               sigma = mtq->getsigTimpeak();
00431             }
00432           if (_debug==1)std::cout <<"-- debug test -- inside loop 4  -- peak:"<<peak<< std::endl; 
00433           if( mtq->doFit()==0 ) 
00434             {
00435               fit  = mtq->getTimax();
00436               ampl = mtq->getAmpl(); 
00437               fwhm = mtq->getFwhm();
00438               fw20 = mtq->getWidth20();
00439               fw80 = mtq->getWidth80();
00440               sliding = mtq->getSlide();
00441             }
00442           if (_debug==1)std::cout <<"-- debug test -- inside loop 4  -- ampl:"<<ampl<< std::endl; 
00443           if( mtq->compute_trise()==0 ) 
00444             {
00445               trise = mtq->getTrise();
00446             }
00447           if (_debug==1)std::cout <<"-- debug test -- inside loop 5  -- trise:"<<trise<< std::endl; 
00448         }
00449       
00450       if (_debug==1)std::cout <<"-- debug test -- inside loop 6  -- status:"<<status<< std::endl; 
00451     
00452       if( status == 1 && mtq->findPeak()==0 ){
00453         
00454         int firstS=int(peak-double(_timebefmax));
00455         int lastS=int(peak+double(_timeaftmax));
00456         
00457 
00458         // Fill histo if there are enough samples
00459         if (_debug==1)std::cout <<"-- debug test -- inside loop 7  -- firstS:"<<firstS<<", nsamples:"<< nsamples<< std::endl;
00460 
00461         if(firstS>=0 && lastS<=nsamples){
00462           
00463           for (int i=firstS;i<lastS;i++){
00464             shapeMatTmp->Fill(double(i)-firstS,matacq[i]);
00465           }
00466           
00467 
00468         }else{  // else extrapolate 
00469           
00470           
00471           int firstSBis;
00472           
00473           if(firstS<0){ // fill first bins with 0
00474 
00475             double thisped;
00476             thisped=(matacq[0]+matacq[1]+matacq[2]+matacq[4]+matacq[5])/5.0;
00477 
00478             for(int i=firstS;i<0;i++){
00479               shapeMatTmp->Fill(double(i)-firstS,thisped);
00480               
00481             }
00482             firstSBis=0; 
00483             
00484           }else{
00485             firstSBis=firstS;
00486           }
00487           
00488           if(lastS>nsamples){
00489 
00490             for(int i=firstSBis;i<int(nsamples);i++){
00491               shapeMatTmp->Fill(double(i)-firstS,matacq[i]);
00492             }
00493             
00494             //extrapolate with expo tail
00495 
00496             double expb=0.998;
00497             double matacqval=expb*matacq[nsamples-1];
00498             
00499             for(int i=nsamples;i<lastS;i++){
00500               shapeMatTmp->Fill(double(i)-firstS,matacqval);
00501               matacqval*=expb;
00502             }
00503             
00504           }else{    
00505             for (int i=firstSBis;i<lastS;i++){
00506               shapeMatTmp->Fill(double(i)-firstS,matacq[i]);
00507             }       
00508           }
00509         }
00510         
00511       }
00512       if (_debug==1)std::cout <<"-- debug test -- inside loop 8"<< std::endl;
00513 
00514       // get back color
00515  
00516       int iCol=nCol;
00517       for(unsigned int i=0;i<nCol;i++){
00518         if(color==colors[i]) {
00519           iCol=i;
00520           i=nCol;
00521         }
00522       }
00523       if (_debug==1)std::cout <<"-- debug test -- inside loop 8bis color:"<<color<<" iCol:"<<iCol<<" nCol:"<< nCol<< std::endl;
00524       
00525       // fill TMTQ 
00526 
00527       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);
00528       
00529       // fill the output tree
00530       
00531       if (_debug==1)std::cout <<"-- debug test -- inside loop 9"<< std::endl;
00532       mtqShape->Fill();
00533       
00534       // clean up
00535       delete mtq;
00536     }
00537   
00538   if (_debug==1)std::cout <<"-- debug test -- after loop "<< std::endl;
00539   sampFile->Close();
00540   
00541   double Peak[6], Sigma[6], Fit[6], Ampl[6], Trise[6], Fwhm[6], Fw20[6], Fw80[6], Ped[6], Pedsig[6], Sliding[6];
00542   int  Side;
00543   
00544   for (unsigned int iColor=0;iColor<nCol;iColor++){
00545     
00546     std::stringstream nametree;
00547     nametree <<"MatacqCol"<<colors[iColor];
00548     meanTree[iColor]= new TTree(nametree.str().c_str(),nametree.str().c_str());
00549     meanTree[iColor]->Branch( "side",        &Side ,       "Side/I"           );
00550     meanTree[iColor]->Branch( "peak",        &Peak ,       "Peak[6]/D"        );
00551     meanTree[iColor]->Branch( "sigma",       &Sigma ,      "Sigma[6]/D"       );
00552     meanTree[iColor]->Branch( "fit",         &Fit ,        "Fit[6]/D"         );
00553     meanTree[iColor]->Branch( "ampl",        &Ampl ,       "Ampl[6]/D"        );
00554     meanTree[iColor]->Branch( "trise",       &Trise ,      "Trise[6]/D"       );  
00555     meanTree[iColor]->Branch( "fwhm",        &Fwhm ,       "Fwhm[6]/D"        );    
00556     meanTree[iColor]->Branch( "fw20",        &Fw20 ,       "Fw20[6]/D"        );  
00557     meanTree[iColor]->Branch( "fw80",        &Fw80 ,       "Fw80[6]/D"        );  
00558     meanTree[iColor]->Branch( "ped",         &Ped ,        "Ped[6]/D"         );   
00559     meanTree[iColor]->Branch( "pedsig",      &Pedsig ,     "Pedsig[6]/D"      );   
00560     meanTree[iColor]->Branch( "sliding",     &Sliding ,    "Sliding[6]/D"     );  
00561     
00562     meanTree[iColor]->SetBranchAddress( "side",        &Side       );
00563     meanTree[iColor]->SetBranchAddress( "peak",        Peak        ); 
00564     meanTree[iColor]->SetBranchAddress( "sigma",       Sigma       ); 
00565     meanTree[iColor]->SetBranchAddress( "fit",         Fit         ); 
00566     meanTree[iColor]->SetBranchAddress( "ampl",        Ampl        ); 
00567     meanTree[iColor]->SetBranchAddress( "fwhm",        Fwhm        ); 
00568     meanTree[iColor]->SetBranchAddress( "fw20",        Fw20        ); 
00569     meanTree[iColor]->SetBranchAddress( "fw80",        Fw80        ); 
00570     meanTree[iColor]->SetBranchAddress( "trise",       Trise       ); 
00571     meanTree[iColor]->SetBranchAddress( "ped",         Ped         ); 
00572     meanTree[iColor]->SetBranchAddress( "pedsig",      Pedsig      ); 
00573     meanTree[iColor]->SetBranchAddress( "sliding",     Sliding     );
00574     
00575   }
00576 
00577   for(unsigned int iCol=0;iCol<nCol;iCol++){
00578     for(unsigned int iSide=0;iSide<nSides;iSide++){
00579       
00580       Side=iSide;
00581       std::vector<double> val[TMTQ::nOutVar];
00582       
00583       for(int iVar=0;iVar<TMTQ::nOutVar;iVar++){
00584         val[iVar] = MTQ[iCol][iSide]->get(iVar);
00585 
00586         for(unsigned int i=0;i<val[iVar].size();i++){
00587           
00588           switch (iVar){
00589             
00590           case TMTQ::iPeak: Peak[i]=val[iVar].at(i);
00591           case TMTQ::iSigma: Sigma[i]=val[iVar].at(i);
00592           case TMTQ::iFit: Fit[i]=val[iVar].at(i);
00593           case TMTQ::iAmpl: Ampl[i]=val[iVar].at(i);
00594           case TMTQ::iFwhm: Fwhm[i]=val[iVar].at(i);
00595           case TMTQ::iFw20: Fw20[i]=val[iVar].at(i);
00596           case TMTQ::iFw80: Fw80[i]=val[iVar].at(i);
00597           case TMTQ::iTrise: Trise[i]=val[iVar].at(i);
00598           case TMTQ::iPed: Ped[i]=val[iVar].at(i);
00599           case TMTQ::iPedsig: Pedsig[i]=val[iVar].at(i);
00600           case TMTQ::iSlide: Sliding[i]=val[iVar].at(i);
00601           }
00602         }
00603       }
00604       meanTree[iCol]->Fill();
00605       if (_debug==1)std::cout <<"-- debug test -- inside final loop  "<< std::endl;
00606     }
00607   }
00608   
00609   // Calculate maximum with pol 2
00610   
00611   int im = shapeMatTmp->GetMaximumBin();
00612   double q1=shapeMatTmp->GetBinContent(im-1);
00613   double q2=shapeMatTmp->GetBinContent(im);
00614   double q3=shapeMatTmp->GetBinContent(im+1);
00615 
00616   double a2=(q3+q1)/2.0-q2;
00617   double a1=q2-q1+a2*(1-2*im);
00618   double a0=q2-a1*im-a2*im*im;
00619   
00620   double am=a0-a1*a1/(4*a2);
00621   
00622 
00623   // Compute pedestal 
00624 
00625   double bl=0;
00626   for (unsigned int i=1; i<_presampleshape+1;i++){ 
00627     bl+=shapeMatTmp->GetBinContent(i);
00628   }
00629   bl/=_presampleshape;
00630 
00631   // Compute and save laser shape
00632 
00633   if (_debug==1)std::cout <<"-- debug test -- computing shape  "<< std::endl;
00634 
00635   int firstBin=0;
00636   double height=0.0;
00637   
00638   for (unsigned int i=_timebefmax; i>_presampleshape;i--){ 
00639     height=shapeMatTmp->GetBinContent(i)-bl;
00640     
00641     if(height<(am-bl)*_cutwindow){
00642       firstBin=i;
00643       i=_presampleshape;
00644     }
00645   }
00646   
00647   unsigned int lastBin=firstBin+_nsamplesshape;
00648   
00649   for(unsigned int i=firstBin;i<lastBin;i++){
00650     shapeMat->Fill(i-firstBin,shapeMatTmp->GetBinContent(i)-bl);
00651   }
00652   
00653   mtqShape->Write();
00654   for (unsigned int iColor=0;iColor<nCol;iColor++){
00655     meanTree[iColor]->Write();
00656   }
00657   if (_debug==1)std::cout <<"-- debug test -- writing  "<< std::endl;
00658   shapeMat->Write();
00659   
00660   // close the output file
00661   outFile->Close();
00662   
00663   // Remove temporary file    
00664   FILE *test; 
00665   test = fopen(sampfile.c_str(),"r");
00666   if (test){
00667   std::stringstream del2;
00668   del2 << "rm " <<sampfile;
00669   system(del2.str().c_str()); 
00670   }
00671 
00672   std::cout <<   "\t+=+    .................... done  +=+" << std::endl;
00673   std::cout <<   "\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << std::endl;
00674  
00675 }
00676 
00677 DEFINE_FWK_MODULE(EcalMatacqAnalyzer);
00678
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/CalibCalorimetry/EcalLaserAnalyzer/plugins/EcalMatacqAnalyzer.cc
