Reference Manual

Go to the documentation of this file.
00001 /* 
00002  *  $id: EcalTestPulseAnalyzer.cc,v 1.4 2007/08/30 17:06:41 ganzhur Exp $
00003  *
00004  *  \class EcalTestPulseAnalyzer
00005  *
00006  *  $Date: 2012/02/09 10:07:37 $
00007  *  primary author: Julie Malcles - CEA/Saclay
00008  *  author: Gautier Hamel De Monchenault - CEA/Saclay
00009  */
00010 #include <TFile.h>
00011 #include <TTree.h>
00012 
00013 using namespace std;
00014 #include "EcalTestPulseAnalyzer.h"
00015 
00016 #include <sstream>
00017 #include <iostream>
00018 #include <iomanip>
00019 
00020 
00021 
00022 #include <FWCore/MessageLogger/interface/MessageLogger.h>
00023 #include <FWCore/Utilities/interface/Exception.h>
00024 
00025 #include <FWCore/Framework/interface/Event.h>
00026 #include <FWCore/Framework/interface/MakerMacros.h>
00027 #include <FWCore/Framework/interface/ESHandle.h>
00028 #include <FWCore/ParameterSet/interface/ParameterSet.h>
00029 #include <FWCore/Framework/interface/EventSetup.h>
00030 
00031 #include <DataFormats/EcalDigi/interface/EcalDigiCollections.h>
00032 #include <DataFormats/EcalDetId/interface/EcalDetIdCollections.h>
00033 #include <DataFormats/EcalRawData/interface/EcalRawDataCollections.h>
00034 
00035 #include <Geometry/EcalMapping/interface/EcalElectronicsMapping.h>
00036 #include <Geometry/EcalMapping/interface/EcalMappingRcd.h>
00037 #include <DataFormats/EcalDetId/interface/EcalElectronicsId.h>
00038 
00039 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TPNFit.h>
00040 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TSFit.h>  
00041 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TMom.h>  
00042 
00043 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/ME.h>
00044 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/MEGeom.h>
00045 
00046 using namespace std;
00047 
00048 
00049 
00050 //========================================================================
00051 EcalTestPulseAnalyzer::EcalTestPulseAnalyzer(const edm::ParameterSet& iConfig)
00052   //========================================================================
00053   :
00054 iEvent(0),
00055 
00056 // framework parameters with default values
00057 
00058 _nsamples(      iConfig.getUntrackedParameter< unsigned int >( "nSamples",       10 ) ),
00059 _presample(     iConfig.getUntrackedParameter< unsigned int >( "nPresamples",     3 ) ),
00060 _firstsample(   iConfig.getUntrackedParameter< unsigned int >( "firstSample",     1 ) ),
00061 _lastsample(    iConfig.getUntrackedParameter< unsigned int >( "lastSample",      2 ) ), 
00062 _samplemin(     iConfig.getUntrackedParameter< unsigned int >( "sampleMin",       3 ) ),
00063 _samplemax(     iConfig.getUntrackedParameter< unsigned int >( "sampleMax",       9 ) ),
00064 _nsamplesPN(    iConfig.getUntrackedParameter< unsigned int >( "nSamplesPN",     50 ) ), 
00065 _presamplePN(   iConfig.getUntrackedParameter< unsigned int >( "nPresamplesPN",   6 ) ),
00066 _firstsamplePN( iConfig.getUntrackedParameter< unsigned int >( "firstSamplePN",   7 ) ),
00067 _lastsamplePN(  iConfig.getUntrackedParameter< unsigned int >( "lastSamplePN",    8 ) ),
00068 _niter(         iConfig.getUntrackedParameter< unsigned int >( "nIter",           3 ) ),
00069 _chi2max(       iConfig.getUntrackedParameter< double       >( "chi2Max",      10.0 ) ),
00070 _timeofmax(     iConfig.getUntrackedParameter< double       >( "timeOfMax",     4.5 ) ),
00071 _ecalPart(      iConfig.getUntrackedParameter< std::string  >( "ecalPart",     "EB" ) ),
00072 _fedid(         iConfig.getUntrackedParameter< int          >( "fedID",        -999 ) ),
00073 nCrys(                                                                         NCRYSEB),
00074 nTT(                                                                             NTTEB),
00075 nMod(                                                                           NMODEB),
00076 nGainPN(                                                                       NGAINPN),
00077 nGainAPD(                                                                     NGAINAPD),
00078 towerID(-1), channelID(-1),runType(-1), runNum(0), fedID(-1), dccID(-1), side(-1), iZ(1),
00079 phi(-1), eta(-1), event(0), apdAmpl(0),apdTime(0),pnAmpl(0),
00080 pnID(-1), moduleID(-1), channelIteratorEE(0)
00081 
00082 
00083   //========================================================================
00084 
00085 {
00086 
00087   //now do what ever initialization is needed
00088 
00089   resdir_                = iConfig.getUntrackedParameter<std::string>("resDir");
00090 
00091   digiCollection_        = iConfig.getParameter<std::string>("digiCollection");
00092   digiPNCollection_      = iConfig.getParameter<std::string>("digiPNCollection");
00093   digiProducer_          = iConfig.getParameter<std::string>("digiProducer");
00094 
00095   eventHeaderCollection_ = iConfig.getParameter<std::string>("eventHeaderCollection");
00096   eventHeaderProducer_   = iConfig.getParameter<std::string>("eventHeaderProducer");
00097 
00098 
00099   // Define geometrical constants 
00100 
00101   if (_ecalPart == "EB") {
00102     nCrys     = NCRYSEB;
00103     nTT       = NTTEB;
00104   } else {
00105     nCrys     = NCRYSEE;
00106     nTT       = NTTEE;
00107   }
00108   
00109   iZ        = 1;
00110   if(  _fedid <= 609 ) iZ=-1;
00111   
00112   dccMEM    = ME::memFromDcc(_fedid);
00113   modules   = ME::lmmodFromDcc(_fedid);
00114   nMod      = modules.size(); 
00115 
00116 }
00117 
00118 
00119 //========================================================================
00120 EcalTestPulseAnalyzer::~EcalTestPulseAnalyzer(){
00121   //========================================================================
00122 
00123   // do anything here that needs to be done at desctruction time
00124   // (e.g. close files, deallocate resources etc.)
00125 
00126 }
00127 
00128 
00129 
00130 //========================================================================
00131 void EcalTestPulseAnalyzer::beginJob() {
00132   //========================================================================
00133 
00134   // Define temporary file
00135 
00136   rootfile=resdir_;
00137   rootfile+="/TmpTreeTestPulseAnalyzer.root";
00138 
00139   outFile = new TFile(rootfile.c_str(),"RECREATE");
00140   
00141   for (unsigned int i=0;i<nCrys;i++){
00142     
00143     std::stringstream name;
00144     name << "Tree" <<i;
00145     
00146     trees[i]= new TTree(name.str().c_str(),name.str().c_str());
00147     
00148     //List of branches
00149 
00150     trees[i]->Branch( "iphi",      &phi,         "phi/I"         );
00151     trees[i]->Branch( "ieta",      &eta,         "eta/I"         );
00152     trees[i]->Branch( "side",      &side,        "side/I"        );
00153     trees[i]->Branch( "dccID",     &dccID,       "dccID/I"       );
00154     trees[i]->Branch( "towerID",   &towerID,     "towerID/I"     );
00155     trees[i]->Branch( "channelID", &channelID,   "channelID/I"   );
00156     trees[i]->Branch( "event",     &event,       "event/I"       );
00157     trees[i]->Branch( "apdGain",   &apdGain,     "apdGain/I"     );
00158     trees[i]->Branch( "pnGain",    &pnGain,      "pnGain/I"      );
00159     trees[i]->Branch( "apdAmpl",   &apdAmpl,     "apdAmpl/D"     );
00160     trees[i]->Branch( "pnAmpl0",   &pnAmpl0,     "pnAmpl0/D"     ); 
00161     trees[i]->Branch( "pnAmpl1",   &pnAmpl1,     "pnAmpl1/D"     ); 
00162   
00163     trees[i]->SetBranchAddress( "ieta",        &eta         );  
00164     trees[i]->SetBranchAddress( "iphi",        &phi         ); 
00165     trees[i]->SetBranchAddress( "side",        &side        ); 
00166     trees[i]->SetBranchAddress( "dccID",       &dccID       ); 
00167     trees[i]->SetBranchAddress( "towerID",     &towerID     ); 
00168     trees[i]->SetBranchAddress( "channelID",   &channelID   ); 
00169     trees[i]->SetBranchAddress( "event",       &event       );
00170     trees[i]->SetBranchAddress( "apdGain",     &apdGain     );
00171     trees[i]->SetBranchAddress( "pnGain",      &pnGain      );
00172     trees[i]->SetBranchAddress( "apdAmpl",     &apdAmpl     );
00173     trees[i]->SetBranchAddress( "pnAmpl0",     &pnAmpl0     );
00174     trees[i]->SetBranchAddress( "pnAmpl1",     &pnAmpl1     );
00175 
00176   }  
00177   
00178   // Initializations 
00179 
00180   for(unsigned int j=0;j<nCrys;j++){
00181     iEta[j]=-1;
00182     iPhi[j]=-1;
00183     iModule[j]=10;
00184     iTowerID[j]=-1;
00185     iChannelID[j]=-1;
00186     idccID[j]=-1;
00187     iside[j]=-1;
00188   }
00189 
00190   for(unsigned int j=0;j<nMod;j++){
00191     firstChanMod[j]=0;
00192     isFirstChanModFilled[j]=0;
00193   }
00194 
00195   // Define output results file name
00196   
00197   std::stringstream namefile;
00198   namefile << resdir_ <<"/APDPN_TESTPULSE.root";      
00199   resfile=namefile.str();
00200 
00201   // TP events counter
00202   TPEvents=0;
00203   
00204 }
00205 
00206 
00207 //========================================================================
00208 void EcalTestPulseAnalyzer:: analyze( const edm::Event & e, const  edm::EventSetup& c){
00209   //========================================================================
00210 
00211   ++iEvent;
00212 
00213   // Retrieve DCC header
00214   edm::Handle<EcalRawDataCollection> pDCCHeader;
00215   const  EcalRawDataCollection* DCCHeader=0;
00216   try {
00217     e.getByLabel(eventHeaderProducer_,eventHeaderCollection_, pDCCHeader);
00218     DCCHeader=pDCCHeader.product();
00219   }catch ( std::exception& ex ) {
00220     std::cerr << "Error! can't get the product  retrieving DCC header" << eventHeaderCollection_.c_str() << std::endl;
00221   }
00222  
00223 
00224   // retrieving crystal EB data from Event
00225   edm::Handle<EBDigiCollection>  pEBDigi;
00226   const  EBDigiCollection* EBDigi=0;
00227 
00228   // retrieving crystal EE data from Event
00229   edm::Handle<EEDigiCollection>  pEEDigi;
00230   const  EEDigiCollection* EEDigi=0;
00231 
00232   if (_ecalPart == "EB") {
00233     try {
00234       e.getByLabel(digiProducer_,digiCollection_, pEBDigi); 
00235       EBDigi=pEBDigi.product(); 
00236     }catch ( std::exception& ex ) {
00237       std::cerr << "Error! can't get the product retrieving EB crystal data " << digiCollection_.c_str() << std::endl;
00238     } 
00239   } else {
00240     try {
00241       e.getByLabel(digiProducer_,digiCollection_, pEEDigi); 
00242       EEDigi=pEEDigi.product(); 
00243     }catch ( std::exception& ex ) {
00244       std::cerr << "Error! can't get the product retrieving EE crystal data " << digiCollection_.c_str() << std::endl;
00245     } 
00246   }
00247 
00248 
00249   // Retrieve crystal PN diodes from Event
00250   edm::Handle<EcalPnDiodeDigiCollection>  pPNDigi;
00251   const  EcalPnDiodeDigiCollection* PNDigi=0;
00252   try {
00253     e.getByLabel(digiProducer_,digiPNCollection_, pPNDigi);
00254     PNDigi=pPNDigi.product(); 
00255   }catch ( std::exception& ex ) {
00256     std::cerr << "Error! can't get the product " << digiCollection_.c_str() << std::endl;
00257   }
00258   
00259 
00260   // retrieving electronics mapping
00261   edm::ESHandle< EcalElectronicsMapping > ecalmapping;
00262   const EcalElectronicsMapping* TheMapping=0; 
00263   try{
00264     c.get< EcalMappingRcd >().get(ecalmapping);
00265     TheMapping = ecalmapping.product();
00266   }catch ( std::exception& ex ) {
00267     std::cerr << "Error! can't get the product EcalMappingRcd"<< std::endl;
00268   }
00269 
00270 
00271   // ====================================
00272   // Decode Basic DCCHeader Information 
00273   // ====================================
00274 
00275   for ( EcalRawDataCollection::const_iterator headerItr= DCCHeader->begin();headerItr != DCCHeader->end(); 
00276         ++headerItr ) {
00277 
00278     int fed = headerItr->fedId();  
00279     
00280     if(fed!=_fedid && _fedid!=-999) continue; 
00281     
00282     runType=headerItr->getRunType();
00283     runNum=headerItr->getRunNumber();
00284     event=headerItr->getLV1();
00285 
00286     dccID=headerItr->getDccInTCCCommand();
00287     fedID=headerItr->fedId();  
00288 
00289 
00290     if( 600+dccID != fedID ) continue;
00291 
00292     // Cut on runType
00293 
00294     if(runType!=EcalDCCHeaderBlock::TESTPULSE_MGPA && runType!=EcalDCCHeaderBlock::TESTPULSE_GAP 
00295       && runType!=EcalDCCHeaderBlock::TESTPULSE_SCAN_MEM  ) return;
00296      
00297   }
00298   
00299   // Cut on fedID
00300 
00301   if(fedID!=_fedid && _fedid!=-999) return; 
00302   
00303   // Count TP events
00304   TPEvents++;
00305 
00306   // ======================
00307   // Decode PN Information
00308   // ======================
00309   
00310   TPNFit * pnfit = new TPNFit();
00311   pnfit -> init(_nsamplesPN,_firstsamplePN,  _lastsamplePN);
00312   
00313   double chi2pn=0;
00314   double ypnrange[50];
00315   double dsum=0.;
00316   double dsum1=0.;
00317   double bl=0.;
00318   double val_max=0.;
00319   int samplemax=0;
00320   unsigned int k;
00321   int pnG[50];
00322   int pngain=0;
00323   
00324   std::map <int, std::vector<double> > allPNAmpl;
00325   std::map <int, std::vector<int> > allPNGain;
00326   
00327   for ( EcalPnDiodeDigiCollection::const_iterator pnItr = PNDigi->begin(); pnItr != PNDigi->end(); ++pnItr ) { // Loop on PNs 
00328     
00329     EcalPnDiodeDetId pnDetId = EcalPnDiodeDetId((*pnItr).id());
00330 
00331     bool isMemRelevant=false;
00332     for (unsigned int imem=0;imem<dccMEM.size();imem++){
00333       if(pnDetId.iDCCId() == dccMEM[imem]) {
00334         isMemRelevant=true;
00335       }
00336     }
00337     
00338     // skip mem dcc without relevant data
00339     if(!isMemRelevant) continue;
00340     
00341     for ( int samId=0; samId < (*pnItr).size() ; samId++ ) { // Loop on PN samples  
00342       pn[samId]=(*pnItr).sample(samId).adc(); 
00343       pnG[samId]=(*pnItr).sample(samId).gainId();
00344       
00345       if(pnG[samId]!=1) std::cout << "PN gain different from 1 for sample "<<samId<< std::endl;
00346       if (samId==0) pngain=pnG[samId];
00347       if (samId>0) pngain=TMath::Max(pnG[samId],pngain); 
00348     }
00349     
00350     for(dsum=0.,k=0;k<_presamplePN;k++) {
00351       dsum+=pn[k];                                                        
00352     }
00353     bl=dsum/((double) _presamplePN);
00354     
00355 
00356     for(val_max=0.,k=0;k<_nsamplesPN;k++) {
00357       ypnrange[k]=pn[k] - bl;
00358       
00359       if(ypnrange[k] > val_max) {
00360         val_max= ypnrange[k]; samplemax=k; 
00361       } 
00362     }
00363     
00364     chi2pn = pnfit -> doFit(samplemax,&ypnrange[0]); 
00365     
00366     if(chi2pn == 101 || chi2pn == 102 || chi2pn == 103) pnAmpl=0.;
00367     else pnAmpl= pnfit -> getAmpl();
00368 
00369     allPNAmpl[pnDetId.iDCCId()].push_back(pnAmpl);
00370     allPNGain[pnDetId.iDCCId()].push_back(pngain); 
00371     
00372   }
00373     
00374   // ===========================
00375   // Decode EBDigis Information
00376   // ===========================
00377   
00378   TSFit * pstpfit = new TSFit(_nsamples,650);
00379   pstpfit -> set_params(_nsamples, _niter, _presample, _samplemin, _samplemax, _timeofmax ,  _chi2max, _firstsample,  _lastsample);
00380   pstpfit -> init_errmat(10.);
00381 
00382   double chi2=0;
00383   double yrange[10]; 
00384   int  adcgain=0;
00385   int  adcG[10];
00386 
00387   
00388   if (EBDigi){
00389     for ( EBDigiCollection::const_iterator digiItr= EBDigi->begin(); digiItr != EBDigi->end(); ++digiItr ) {  // Loop on EB crystals
00390 
00391       EBDetId id_crystal(digiItr->id()) ;
00392       EBDataFrame df( *digiItr );
00393 
00394       int etaG = id_crystal.ieta() ;  // global
00395       int phiG = id_crystal.iphi() ;  // global
00396 
00397       int etaL ;  // local
00398       int phiL ;  // local
00399       std::pair<int, int> LocalCoord=MEEBGeom::localCoord( etaG , phiG );
00400       
00401       etaL=LocalCoord.first ;
00402       phiL=LocalCoord.second ;
00403       
00404       eta = etaG;
00405       phi = phiG;
00406 
00407       side=MEEBGeom::side(etaG,phiG);
00408       
00409       EcalElectronicsId elecid_crystal = TheMapping->getElectronicsId(id_crystal);
00410 
00411       towerID=elecid_crystal.towerId(); 
00412       int strip=elecid_crystal.stripId();
00413       int xtal=elecid_crystal.xtalId(); 
00414       channelID=  5*(strip-1) + xtal-1; // FIXME
00415 
00416   
00417       int module= MEEBGeom::lmmod(etaG, phiG);
00418       int iMod = module-1;
00419 
00420       assert( module>=*min_element(modules.begin(),modules.end()) && module<=*max_element(modules.begin(),modules.end()) );
00421 
00422 
00423       std::pair<int,int> pnpair=MEEBGeom::pn(module);
00424       unsigned int MyPn0=pnpair.first;
00425       unsigned int MyPn1=pnpair.second;
00426                 
00427       unsigned int channel=MEEBGeom::electronic_channel( etaL, phiL );
00428       
00429       if(isFirstChanModFilled[iMod]==0) {
00430         firstChanMod[iMod]=channel;
00431         isFirstChanModFilled[iMod]=1;
00432       }
00433 
00434       iEta[channel]=eta;
00435       iPhi[channel]=phi;
00436       iModule[channel]= module ;
00437       iTowerID[channel]=towerID;
00438       iChannelID[channel]=channelID;
00439       idccID[channel]=dccID;
00440       iside[channel]=side;
00441 
00442   
00443   
00444       // get adc samples 
00445       //====================
00446       
00447       for (unsigned int i=0; i< (*digiItr).size() ; ++i ) { 
00448 
00449         EcalMGPASample samp_crystal(df.sample(i));
00450         adc[i]=samp_crystal.adc() ;  
00451         adcG[i]=samp_crystal.gainId();   
00452   
00453         if (i==0) adcgain=adcG[i];
00454         if (i>0) adcgain=TMath::Max(adcG[i],adcgain);      
00455       }
00456       // Remove pedestal
00457       //====================
00458       for(dsum=0.,dsum1=0.,k=0;k<_presample;k++) {
00459         dsum+=adc[k]; 
00460         if(k<_presample-1) dsum1+=adc[k];
00461       }
00462       
00463       bl=dsum/((double)_presample);
00464       
00465       for(val_max=0.,k=0;k<_nsamples;k++) {
00466         yrange[k]=adc[k] - bl;
00467         if(yrange[k] > val_max) {
00468           val_max= yrange[k]; samplemax=k;
00469         }
00470       } 
00471       
00472       apdGain=adcgain;
00473       
00474       if(allPNAmpl[dccMEM[0]].size()>MyPn0) pnAmpl0=allPNAmpl[dccMEM[0]][MyPn0];
00475       else pnAmpl0=0;
00476       if(allPNAmpl[dccMEM[0]].size()>MyPn1) pnAmpl1=allPNAmpl[dccMEM[0]][MyPn1];
00477       else pnAmpl1=0;
00478 
00479       if(allPNGain[dccMEM[0]].size()>MyPn0) pnGain=allPNGain[dccMEM[0]][MyPn0];
00480       else pnGain=0;
00481 
00482       // Perform the fit on apd samples
00483       //================================
00484 
00485       chi2 =  pstpfit -> fit_third_degree_polynomial(&yrange[0],ret_data);
00486       
00487       //Retrieve APD amplitude from fit
00488       //================================
00489 
00490       if( val_max > 100000. || chi2 < 0. || chi2 == 102 ) {
00491       
00492         apdAmpl=0;
00493         apdTime=0;
00494 
00495       }else{
00496       
00497         apdAmpl = ret_data[0];
00498         apdTime = ret_data[1];
00499         
00500       }
00501 
00502       trees[channel]->Fill();    
00503       
00504     }
00505     
00506   } else {
00507 
00508     for ( EEDigiCollection::const_iterator digiItr= EEDigi->begin(); digiItr != EEDigi->end(); ++digiItr ) {  // Loop on EE crystals
00509       
00510       EEDetId id_crystal(digiItr->id()) ;
00511       EEDataFrame df( *digiItr );
00512 
00513       phi = id_crystal.ix() ; 
00514       eta = id_crystal.iy() ; 
00515 
00516       int iX = (phi-1)/5+1;
00517       int iY = (eta-1)/5+1;
00518   
00519       side=MEEEGeom::side( iX, iY ,iZ); 
00520 
00521 
00522       // Recover the TT id and the electronic crystal numbering from EcalElectronicsMapping
00523 
00524       EcalElectronicsId elecid_crystal = TheMapping->getElectronicsId(id_crystal);
00525     
00526 
00527       towerID=elecid_crystal.towerId();
00528       channelID=elecid_crystal.channelId()-1;  
00529       
00530       int module=MEEEGeom::lmmod( iX, iY );
00531       if( module>=18 && side==1 ) module+=2;  // Trick to fix endcap specificity
00532       int iMod = module-1;
00533 
00534       assert( module>=*min_element(modules.begin(),modules.end()) && module<=*max_element(modules.begin(),modules.end()) );
00535 
00536       std::pair<int,int> pnpair=MEEEGeom::pn( module, _fedid ) ; 
00537 
00538       unsigned int MyPn0=pnpair.first;
00539       unsigned int MyPn1=pnpair.second;
00540 
00541       int hashedIndex=100000*eta+phi;
00542      
00543       if( channelMapEE.count(hashedIndex) == 0 ){
00544         channelMapEE[hashedIndex]=channelIteratorEE;
00545         channelIteratorEE++;
00546       }
00547       
00548       unsigned int channel=channelMapEE[hashedIndex];
00549 
00550       if(isFirstChanModFilled[iMod]==0) {
00551         firstChanMod[iMod]=channel;
00552         isFirstChanModFilled[iMod]=1;
00553       }
00554       
00555       iEta[channel]=eta;
00556       iPhi[channel]=phi;
00557       iModule[channel]= module ;
00558       iTowerID[channel]=towerID;
00559       iChannelID[channel]=channelID;
00560       idccID[channel]=dccID;
00561       iside[channel]=side;
00562 
00563       assert (channel < nCrys);
00564       
00565       // Get adc samples
00566       //====================
00567 
00568       for (unsigned int i=0; i< (*digiItr).size() ; ++i ) { 
00569 
00570         EcalMGPASample samp_crystal(df.sample(i));
00571         adc[i]=samp_crystal.adc() ;  
00572         adcG[i]=samp_crystal.gainId();      
00573 
00574         if (i==0) adcgain=adcG[i];
00575         if (i>0) adcgain=TMath::Max(adcG[i],adcgain); 
00576       }
00577       
00578 
00579       // Remove pedestal
00580       //====================
00581       for(dsum=0.,dsum1=0.,k=0;k<_presample;k++) {
00582         dsum+=adc[k]; 
00583         if(k<_presample-1) dsum1+=adc[k];
00584       }
00585       
00586       bl=dsum/((double)_presample);
00587       
00588       for(val_max=0.,k=0;k<_nsamples;k++) {
00589         yrange[k]=adc[k] - bl;
00590         if(yrange[k] > val_max) {
00591           val_max= yrange[k]; samplemax=k;
00592         }
00593       } 
00594       apdGain=adcgain;
00595       
00596       int dccMEMIndex=0;
00597       if(side==1) dccMEMIndex+=2; // Trick to fix endcap specificity
00598       
00599       if(allPNAmpl[dccMEM[dccMEMIndex]].size()>MyPn0) pnAmpl0=allPNAmpl[dccMEM[dccMEMIndex]][MyPn0];
00600       else pnAmpl0=0;
00601       if(allPNAmpl[dccMEM[dccMEMIndex+1]].size()>MyPn1) pnAmpl1=allPNAmpl[dccMEM[dccMEMIndex+1]][MyPn1];
00602       else pnAmpl1=0;
00603       
00604       if(allPNGain[dccMEM[dccMEMIndex]].size()>MyPn0) pnGain=allPNGain[dccMEM[dccMEMIndex]][MyPn0];
00605       else pnGain=0;
00606 
00607 
00608       // Perform the fit on apd samples
00609       //=================================
00610 
00611       chi2 =  pstpfit -> fit_third_degree_polynomial(&yrange[0],ret_data);
00612       
00613       //Retrieve APD amplitude from fit
00614       //=================================
00615       
00616       if( val_max > 100000. || chi2 < 0. || chi2 == 102 ) {
00617         
00618         apdAmpl=0;
00619         apdTime=0;
00620         
00621       }else{
00622       
00623         apdAmpl = ret_data[0];
00624         apdTime = ret_data[1];
00625         
00626       }
00627 
00628       trees[channel]->Fill();    
00629     }
00630  
00631   }
00632   
00633   
00634 } // end of analyze
00635 
00636 
00637 //========================================================================
00638 void EcalTestPulseAnalyzer::endJob() {
00639   //========================================================================
00640   
00641   // Don't do anything if there is no events
00642   if( TPEvents == 0 ){
00643     
00644     outFile->Close(); 
00645     
00646     // Remove temporary file
00647     
00648     std::stringstream del;
00649     del << "rm " <<rootfile;
00650     system(del.str().c_str()); 
00651     
00652     std::cout << " No TP Events "<< std::endl;
00653     return;
00654   }
00655 
00656   std::cout <<  "\n\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << std::endl;
00657   std::cout <<    "\t+=+     Analyzing test pulse data: getting APD, PN  +=+" << std::endl;
00658   
00659   
00660   // Create output ntuples:
00661   
00662   //std::cout<< "TP Test Name File "<< resfile.c_str() << std::endl;
00663 
00664   resFile = new TFile(resfile.c_str(),"RECREATE");
00665   
00666   restrees= new TTree("TPAPD","TPAPD");
00667   respntrees= new TTree("TPPN","TPPN");
00668 
00669   restrees->Branch( "iphi",        &iphi,        "iphi/I"           );
00670   restrees->Branch( "ieta",        &ieta,        "ieta/I"           );
00671   restrees->Branch( "dccID",       &dccID,       "dccID/I"          );
00672   restrees->Branch( "side",        &side,        "side/I"           );
00673   restrees->Branch( "towerID",     &towerID,     "towerID/I"        );
00674   restrees->Branch( "channelID",   &channelID,   "channelID/I"      );
00675   restrees->Branch( "moduleID",    &moduleID,    "moduleID/I"       );
00676   restrees->Branch( "flag",        &flag,        "flag/I"           );
00677   restrees->Branch( "gain",        &gain,        "gain/I"           );
00678   restrees->Branch( "APD",         &APD,         "APD[6]/D"         );
00679   
00680   respntrees->Branch( "pnID",      &pnID,      "pnID/I"         );
00681   respntrees->Branch( "moduleID",  &moduleID,  "moduleID/I"     );
00682   respntrees->Branch( "gain",      &gain,      "gain/I"         ); 
00683   respntrees->Branch( "PN",        &PN,        "PN[6]/D"        ); 
00684   
00685   restrees->SetBranchAddress( "iphi",        &iphi       );
00686   restrees->SetBranchAddress( "ieta",        &ieta       );
00687   restrees->SetBranchAddress( "dccID",       &dccID      );
00688   restrees->SetBranchAddress( "side",        &side       );
00689   restrees->SetBranchAddress( "towerID",     &towerID    );
00690   restrees->SetBranchAddress( "channelID",   &channelID  );
00691   restrees->SetBranchAddress( "moduleID",    &moduleID   ); 
00692   restrees->SetBranchAddress( "flag",        &flag       );  
00693   restrees->SetBranchAddress( "gain",        &gain       );  
00694   restrees->SetBranchAddress( "APD",         APD         );  
00695   
00696   respntrees->SetBranchAddress( "pnID",      &pnID     );
00697   respntrees->SetBranchAddress( "moduleID",  &moduleID );
00698   respntrees->SetBranchAddress( "gain",      &gain     ); 
00699   respntrees->SetBranchAddress( "PN",        PN        ); 
00700   
00701   
00702 
00703   TMom *APDAnal[1700][10];
00704   TMom *PNAnal[9][2][10];
00705 
00706 
00707   for (unsigned int iMod=0;iMod<nMod;iMod++){
00708     for (unsigned int ich=0;ich<2;ich++){
00709       for (unsigned int ig=0;ig<nGainPN;ig++){
00710         PNAnal[iMod][ich][ig]=new TMom();
00711       }
00712     }
00713   }
00714 
00715   for (unsigned int iCry=0;iCry<nCrys;iCry++){ // Loop on data trees (ie on cristals)
00716 
00717     for(unsigned int iG=0;iG<nGainAPD;iG++){
00718       APDAnal[iCry][iG]=new TMom();
00719     }
00720 
00721 
00722     
00723     // Define submodule and channel number inside the submodule (as Patrice)
00724     
00725     unsigned int iMod=iModule[iCry]-1;
00726   
00727     moduleID=iMod+1;    
00728     if( moduleID>=20 ) moduleID-=2;  // Trick to fix endcap specificity
00729     
00730     Long64_t nbytes = 0, nb = 0;
00731     for (Long64_t jentry=0; jentry< trees[iCry]->GetEntriesFast();jentry++) { 
00732       nb = trees[iCry]->GetEntry(jentry);   nbytes += nb; 
00733       
00734       // PN Means and RMS 
00735 
00736       if( firstChanMod[iMod]==iCry ){ 
00737         PNAnal[iMod][0][pnGain]->addEntry(pnAmpl0);
00738         PNAnal[iMod][1][pnGain]->addEntry(pnAmpl1);
00739       }
00740       
00741       // APD means and RMS
00742       
00743       APDAnal[iCry][apdGain]->addEntry(apdAmpl);
00744       
00745     }
00746     
00747     if (trees[iCry]->GetEntries()<10){
00748       flag=-1;
00749       for (int j=0;j<6;j++){
00750         APD[j]=0.0;
00751       }
00752     }
00753     else flag=1;
00754     
00755     iphi=iPhi[iCry];
00756     ieta=iEta[iCry];
00757     dccID=idccID[iCry];
00758     side=iside[iCry];
00759     towerID=iTowerID[iCry];
00760     channelID=iChannelID[iCry];
00761 
00762     for (unsigned int ig=0;ig<nGainAPD;ig++){
00763       
00764       APD[0]= APDAnal[iCry][ig]->getMean();
00765       APD[1]= APDAnal[iCry][ig]->getRMS();
00766       APD[2]= APDAnal[iCry][ig]->getM3();
00767       APD[3]= APDAnal[iCry][ig]->getNevt();
00768       APD[4]= APDAnal[iCry][ig]->getMin();
00769       APD[5]= APDAnal[iCry][ig]->getMax();
00770       gain=ig;
00771       
00772       // Fill APD tree
00773       
00774       restrees->Fill(); 
00775       
00776     }
00777   }
00778   
00779   // Get final results for PN and PN/PN
00780   
00781   for (unsigned int ig=0;ig<nGainPN;ig++){
00782     for (unsigned int iMod=0;iMod<nMod;iMod++){
00783       for (int ch=0;ch<2;ch++){
00784         
00785         pnID=ch;
00786         moduleID=iMod;
00787         if( moduleID>=20 ) moduleID-=2;  // Trick to fix endcap specificity
00788 
00789         PN[0]= PNAnal[iMod][ch][ig]->getMean();
00790         PN[1]= PNAnal[iMod][ch][ig]->getRMS();
00791         PN[2]= PNAnal[iMod][ch][ig]->getM3();
00792         PN[3]= PNAnal[iMod][ch][ig]->getNevt();
00793         PN[4]= PNAnal[iMod][ch][ig]->getMin();
00794         PN[5]= PNAnal[iMod][ch][ig]->getMax();
00795         gain=ig;
00796         
00797         // Fill PN tree
00798         respntrees->Fill();
00799         
00800       }
00801     }
00802   }
00803   
00804   outFile->Close(); 
00805 
00806   // Remove temporary file
00807 
00808   std::stringstream del;
00809   del << "rm " <<rootfile;
00810   system(del.str().c_str()); 
00811 
00812 
00813   // Save final results 
00814 
00815   restrees->Write();
00816   respntrees->Write();
00817   resFile->Close(); 
00818   
00819   std::cout <<    "\t+=+    ...................................... done  +=+" << std::endl;
00820   std::cout <<    "\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << std::endl;
00821 }
00822 
00823 
00824 DEFINE_FWK_MODULE(EcalTestPulseAnalyzer);
00825
/data/refman/pasoursint/CMSSW_5_3_3/src/CalibCalorimetry/EcalLaserAnalyzer/plugins/EcalTestPulseAnalyzer.cc
