Reference Manual

Go to the documentation of this file.
00001 // This file is imported from:
00002 
00003 // -*- C++ -*-
00004 //
00005 // Package:    EgammaElectronProducers
00006 // Class:      CalibratedElectronProducer
00007 //
00016 #include "EgammaAnalysis/ElectronTools/plugins/CalibratedElectronProducer.h"
00017 
00018 #include "FWCore/Framework/interface/Frameworkfwd.h"
00019 #include "FWCore/Framework/interface/MakerMacros.h"
00020 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00021 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00022 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
00023 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
00024 #include "FWCore/ParameterSet/interface/FileInPath.h"
00025 
00026 #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
00027 #include "DataFormats/EgammaReco/interface/ElectronSeed.h"
00028 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
00029 #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
00030 #include "Geometry/CaloEventSetup/interface/CaloTopologyRecord.h"
00031 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00032 #include "EgammaAnalysis/ElectronTools/interface/SuperClusterHelper.h"
00033 
00034 #include <iostream>
00035 
00036 using namespace edm ;
00037 using namespace std ;
00038 using namespace reco ;
00039 
00040 CalibratedElectronProducer::CalibratedElectronProducer( const edm::ParameterSet & cfg )
00041 {
00042     inputElectrons_ = cfg.getParameter<edm::InputTag>("inputElectronsTag");
00043     
00044     nameEnergyReg_      = cfg.getParameter<edm::InputTag>("nameEnergyReg");
00045     nameEnergyErrorReg_ = cfg.getParameter<edm::InputTag>("nameEnergyErrorReg");
00046     
00047     recHitCollectionEB_ = cfg.getParameter<edm::InputTag>("recHitCollectionEB");
00048     recHitCollectionEE_ = cfg.getParameter<edm::InputTag>("recHitCollectionEE");
00049     
00050     
00051     nameNewEnergyReg_      = cfg.getParameter<std::string>("nameNewEnergyReg");
00052     nameNewEnergyErrorReg_ = cfg.getParameter<std::string>("nameNewEnergyErrorReg");
00053     newElectronName_ = cfg.getParameter<std::string>("outputGsfElectronCollectionLabel");
00054     
00055     
00056     dataset = cfg.getParameter<std::string>("inputDataset");
00057     isMC = cfg.getParameter<bool>("isMC");
00058     updateEnergyError = cfg.getParameter<bool>("updateEnergyError");
00059     lumiRatio = cfg.getParameter<double>("lumiRatio");
00060     correctionsType = cfg.getParameter<int>("correctionsType");
00061     applyLinearityCorrection = cfg.getParameter<bool>("applyLinearityCorrection");
00062     combinationType = cfg.getParameter<int>("combinationType");
00063     verbose = cfg.getParameter<bool>("verbose");
00064     synchronization = cfg.getParameter<bool>("synchronization");
00065     combinationRegressionInputPath = cfg.getParameter<std::string>("combinationRegressionInputPath");
00066     scaleCorrectionsInputPath = cfg.getParameter<std::string>("scaleCorrectionsInputPath");
00067     linCorrectionsInputPath   = cfg.getParameter<std::string>("linearityCorrectionsInputPath");
00068   
00069     //basic checks
00070     if ( isMC && ( dataset != "Summer11" && dataset != "Fall11" 
00071         && dataset!= "Summer12" && dataset != "Summer12_DR53X_HCP2012"
00072         && dataset != "Summer12_LegacyPaper" ) )
00073     { 
00074         throw cms::Exception("CalibratedgsfElectronProducer|ConfigError") << "Unknown MC dataset"; 
00075     }
00076     if ( !isMC && ( dataset != "Prompt" && dataset != "ReReco" 
00077         && dataset != "Jan16ReReco" && dataset != "ICHEP2012" 
00078         && dataset != "Moriond2013" && dataset != "22Jan2013ReReco" ) )
00079     { 
00080         throw cms::Exception("CalibratedgsfElectronProducer|ConfigError") << "Unknown Data dataset"; 
00081     }
00082 
00083     // Linearity correction only applied on combined momentum obtain with regression combination
00084     if(combinationType!=3 && applyLinearityCorrection)
00085     {
00086         std::cout << "[CalibratedElectronProducer] "
00087             << "Warning: you chose combinationType!=3 and applyLinearityCorrection=True. Linearity corrections are only applied on top of combination 3." << std::endl;
00088     }
00089     
00090     std::cout << "[CalibratedGsfElectronProducer] Correcting scale for dataset " << dataset << std::endl;
00091 
00092     //initializations
00093     std::string pathToDataCorr;
00094     switch (correctionsType)
00095     {
00096         case 0:
00097             break;
00098         case 1: 
00099                 if ( verbose ) 
00100             {
00101                 std::cout << "You choose regression 1 scale corrections" << std::endl;
00102             }
00103             break;
00104         case 2: 
00105                 if ( verbose ) 
00106             {
00107                 std::cout << "You choose regression 2 scale corrections." << std::endl;
00108             }
00109             break;
00110         case 3: 
00111             throw cms::Exception("CalibratedgsfElectronProducer|ConfigError")
00112                 << "You choose standard non-regression ecal energy scale corrections. They are not implemented yet."; 
00113                 break;
00114         default: 
00115             throw cms::Exception("CalibratedgsfElectronProducer|ConfigError")
00116                 << "Unknown correctionsType !!!" ;
00117     }
00118   
00119     theEnCorrector = new ElectronEnergyCalibrator
00120         (
00121             edm::FileInPath(scaleCorrectionsInputPath.c_str()).fullPath().c_str(), 
00122             edm::FileInPath(linCorrectionsInputPath.c_str()).fullPath().c_str(),
00123             dataset, 
00124             correctionsType, 
00125             applyLinearityCorrection,
00126             lumiRatio, 
00127             isMC, 
00128             updateEnergyError, 
00129             verbose, 
00130             synchronization
00131         );
00132     
00133     if ( verbose ) 
00134     {
00135         std::cout<<"[CalibratedGsfElectronProducer] "
00136         << "ElectronEnergyCalibrator object is created" << std::endl;
00137     }
00138 
00139     myEpCombinationTool = new EpCombinationTool();
00140     myEpCombinationTool->init
00141         (
00142             edm::FileInPath(combinationRegressionInputPath.c_str()).fullPath().c_str(),
00143             "CombinationWeight"
00144         );
00145 
00146     myCombinator = new ElectronEPcombinator();
00147     
00148     if ( verbose ) 
00149     {
00150         std::cout << "[CalibratedGsfElectronProducer] "
00151         << "Combination tools are created and initialized" << std::endl;
00152     }
00153     
00154     produces<edm::ValueMap<double> >(nameNewEnergyReg_);
00155     produces<edm::ValueMap<double> >(nameNewEnergyErrorReg_);
00156     produces<GsfElectronCollection> (newElectronName_);
00157     geomInitialized_ = false;
00158 }
00159  
00160 CalibratedElectronProducer::~CalibratedElectronProducer()
00161 {}
00162 
00163 void CalibratedElectronProducer::produce( edm::Event & event, const edm::EventSetup & setup )
00164 {
00165     if (!geomInitialized_) 
00166     {
00167         edm::ESHandle<CaloTopology> theCaloTopology;
00168         setup.get<CaloTopologyRecord>().get(theCaloTopology);
00169         ecalTopology_ = & (*theCaloTopology);
00170         
00171         edm::ESHandle<CaloGeometry> theCaloGeometry;
00172         setup.get<CaloGeometryRecord>().get(theCaloGeometry); 
00173         caloGeometry_ = & (*theCaloGeometry);
00174         geomInitialized_ = true;
00175     }
00176 
00177     // Read GsfElectrons
00178     edm::Handle<reco::GsfElectronCollection>  oldElectronsH ;
00179     event.getByLabel(inputElectrons_,oldElectronsH) ;
00180    
00181     // Read RecHits
00182     edm::Handle< EcalRecHitCollection > pEBRecHits;
00183     edm::Handle< EcalRecHitCollection > pEERecHits;
00184     event.getByLabel( recHitCollectionEB_, pEBRecHits );
00185     event.getByLabel( recHitCollectionEE_, pEERecHits );
00186 
00187     // ReadValueMaps
00188     edm::Handle<edm::ValueMap<double> > valMapEnergyH;
00189     event.getByLabel(nameEnergyReg_,valMapEnergyH);
00190     edm::Handle<edm::ValueMap<double> > valMapEnergyErrorH;
00191     event.getByLabel(nameEnergyErrorReg_,valMapEnergyErrorH);
00192   
00193     // Prepare output collections
00194     std::auto_ptr<GsfElectronCollection> electrons( new reco::GsfElectronCollection ) ;
00195     // Fillers for ValueMaps:
00196     std::auto_ptr<edm::ValueMap<double> > regrNewEnergyMap(new edm::ValueMap<double>() );
00197     edm::ValueMap<double>::Filler energyFiller(*regrNewEnergyMap);
00198     
00199     std::auto_ptr<edm::ValueMap<double> > regrNewEnergyErrorMap(new edm::ValueMap<double>() );
00200     edm::ValueMap<double>::Filler energyErrorFiller(*regrNewEnergyErrorMap);
00201     
00202     // first clone the initial collection
00203     unsigned nElectrons = oldElectronsH->size();
00204     for( unsigned iele = 0; iele < nElectrons; ++iele )
00205     {    
00206       electrons->push_back((*oldElectronsH)[iele]);
00207     }
00208 
00209     std::vector<double> regressionValues;
00210     std::vector<double> regressionErrorValues;
00211     regressionValues.reserve(nElectrons);
00212     regressionErrorValues.reserve(nElectrons);
00213 
00214     if ( correctionsType != 0 )
00215     {
00216         for ( unsigned iele = 0; iele < nElectrons ; ++iele) 
00217         {
00218             reco::GsfElectron & ele  ( (*electrons)[iele]);
00219             reco::GsfElectronRef elecRef(oldElectronsH,iele);
00220             double regressionEnergy = (*valMapEnergyH)[elecRef];
00221             double regressionEnergyError = (*valMapEnergyErrorH)[elecRef];
00222             
00223             regressionValues.push_back(regressionEnergy);
00224             regressionErrorValues.push_back(regressionEnergyError);
00225 
00226             //    r9 
00227             const EcalRecHitCollection * recHits=0;
00228             if( ele.isEB() ) 
00229             {
00230                 recHits = pEBRecHits.product();
00231             } else recHits = pEERecHits.product();
00232 
00233             SuperClusterHelper mySCHelper( &(ele), recHits, ecalTopology_, caloGeometry_ );
00234 
00235             int elClass = -1;
00236             int run = event.run(); 
00237             
00238             float r9 = mySCHelper.r9(); 
00239             double correctedEcalEnergy = ele.correctedEcalEnergy();
00240             double correctedEcalEnergyError = ele.correctedEcalEnergyError();
00241             double trackMomentum = ele.trackMomentumAtVtx().R();
00242             double trackMomentumError = ele.trackMomentumError();
00243             double combinedMomentum = ele.p();
00244             double combinedMomentumError = ele.p4Error(ele.candidateP4Kind());
00245             // FIXME : p4Error not filled for pure tracker electrons
00246             // Recompute it using the parametrization implemented in 
00247             // RecoEgamma/EgammaElectronAlgos/src/ElectronEnergyCorrector.cc::simpleParameterizationUncertainty() 
00248             if( !ele.ecalDrivenSeed() )
00249             {
00250                 double error = 999. ;
00251                 double momentum = (combinedMomentum<15. ? 15. : combinedMomentum);
00252                 if ( ele.isEB() )
00253                 {
00254                     float parEB[3] = { 5.24e-02,  2.01e-01, 1.00e-02};
00255                     error = momentum * sqrt( pow(parEB[0]/sqrt(momentum),2) + pow(parEB[1]/momentum,2) + pow(parEB[2],2) );
00256                 }
00257                 else if ( ele.isEE() )
00258                 {
00259                     float parEE[3] = { 1.46e-01, 9.21e-01, 1.94e-03} ;
00260                     error = momentum * sqrt( pow(parEE[0]/sqrt(momentum),2) + pow(parEE[1]/momentum,2) + pow(parEE[2],2) );
00261                 }
00262                 combinedMomentumError = error;
00263             }
00264     
00265             if (ele.classification() == reco::GsfElectron::GOLDEN) {elClass = 0;}
00266             if (ele.classification() == reco::GsfElectron::BIGBREM) {elClass = 1;}
00267             if (ele.classification() == reco::GsfElectron::BADTRACK) {elClass = 2;}
00268             if (ele.classification() == reco::GsfElectron::SHOWERING) {elClass = 3;}
00269             if (ele.classification() == reco::GsfElectron::GAP) {elClass = 4;}
00270 
00271             SimpleElectron mySimpleElectron
00272                 (
00273                     run, 
00274                     elClass, 
00275                     r9, 
00276                     correctedEcalEnergy, 
00277                     correctedEcalEnergyError, 
00278                     trackMomentum, 
00279                     trackMomentumError, 
00280                     regressionEnergy, 
00281                     regressionEnergyError, 
00282                     combinedMomentum, 
00283                     combinedMomentumError, 
00284                     ele.superCluster()->eta(), 
00285                     ele.isEB(), 
00286                     isMC, 
00287                     ele.ecalDriven(), 
00288                     ele.trackerDrivenSeed()
00289                 );
00290 
00291             // energy calibration for ecalDriven electrons
00292             if ( ele.core()->ecalDrivenSeed() || correctionsType==2 || combinationType==3 ) 
00293             {
00294                 theEnCorrector->calibrate(mySimpleElectron);
00295         
00296                 // E-p combination  
00297         
00298                 switch ( combinationType )
00299                 {
00300                     case 0: 
00301                         if ( verbose ) 
00302                         {
00303                             std::cout << "[CalibratedGsfElectronProducer] " 
00304                             << "You choose not to combine." << std::endl;
00305                         }
00306                         break;
00307                     case 1: 
00308                         if ( verbose ) 
00309                         {
00310                             std::cout << "[CalibratedGsfElectronProducer] " 
00311                             << "You choose corrected regression energy for standard combination" << std::endl;
00312                         }
00313                         myCombinator->setCombinationMode(1);
00314                         myCombinator->combine(mySimpleElectron);
00315                         break;
00316                     case 2: 
00317                         if ( verbose ) 
00318                         {
00319                             std::cout << "[CalibratedGsfElectronProducer] " 
00320                             << "You choose uncorrected regression energy for standard combination" << std::endl;
00321                         }
00322                         myCombinator->setCombinationMode(2);
00323                         myCombinator->combine(mySimpleElectron);
00324                         break;
00325                     case 3: 
00326                         if ( verbose ) 
00327                         {
00328                             std::cout << "[CalibratedGsfElectronProducer] "
00329                             << "You choose regression combination." << std::endl;
00330                         }
00331                         myEpCombinationTool->combine(mySimpleElectron);
00332                         theEnCorrector->correctLinearity(mySimpleElectron);
00333                         break;
00334                     default: 
00335                                 throw cms::Exception("CalibratedgsfElectronProducer|ConfigError")
00336                             << "Unknown combination Type !!!" ;
00337                 }
00338 
00339                 math::XYZTLorentzVector oldMomentum = ele.p4() ;
00340                 math::XYZTLorentzVector newMomentum_ ;
00341                 newMomentum_ = math::XYZTLorentzVector
00342                  ( oldMomentum.x()*mySimpleElectron.getCombinedMomentum()/oldMomentum.t(),
00343                    oldMomentum.y()*mySimpleElectron.getCombinedMomentum()/oldMomentum.t(),
00344                    oldMomentum.z()*mySimpleElectron.getCombinedMomentum()/oldMomentum.t(),
00345                    mySimpleElectron.getCombinedMomentum() ) ;
00346 
00347                 ele.correctMomentum
00348                     (
00349                         newMomentum_,
00350                         mySimpleElectron.getTrackerMomentumError(),
00351                         mySimpleElectron.getCombinedMomentumError()
00352                     );
00353 
00354                 if ( verbose ) 
00355                 {
00356                     std::cout << "[CalibratedGsfElectronProducer] Combined momentum after saving "
00357                         << ele.p4().t() << std::endl;
00358                 }
00359             }// end of if (ele.core()->ecalDrivenSeed()) 
00360         }// end of loop on electrons 
00361     } else 
00362     {
00363         if ( verbose ) 
00364         {
00365             std::cout << "[CalibratedGsfElectronProducer] "
00366             << "You choose not to correct. Uncorrected Regression Energy is taken." << std::endl;
00367         }
00368     }
00369 
00370     // Save the electrons
00371     const edm::OrphanHandle<reco::GsfElectronCollection> gsfNewElectronHandle = event.put(electrons, newElectronName_) ;
00372     energyFiller.insert(gsfNewElectronHandle,regressionValues.begin(),regressionValues.end());
00373     energyFiller.fill();
00374     energyErrorFiller.insert(gsfNewElectronHandle,regressionErrorValues.begin(),regressionErrorValues.end());
00375     energyErrorFiller.fill();
00376     
00377     event.put(regrNewEnergyMap,nameNewEnergyReg_);
00378     event.put(regrNewEnergyErrorMap,nameNewEnergyErrorReg_);
00379 }
00380 
00381 #include "FWCore/Framework/interface/MakerMacros.h"
00382 #include "FWCore/Framework/interface/ESProducer.h"
00383 #include "FWCore/Framework/interface/ModuleFactory.h"
00384 DEFINE_FWK_MODULE(CalibratedElectronProducer);
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_14/src/EgammaAnalysis/ElectronTools/plugins/CalibratedElectronProducer.cc
