CalibratedElectronProducer.cc

Go to the documentation of this file.

// This file is imported from:

// -*- C++ -*-
//
// Package:    EgammaElectronProducers
// Class:      CalibratedElectronProducer
//
#include "EgammaAnalysis/ElectronTools/plugins/CalibratedElectronProducer.h"

#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
#include "FWCore/ParameterSet/interface/FileInPath.h"

#include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
#include "DataFormats/EgammaReco/interface/ElectronSeed.h"
#include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
#include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
#include "Geometry/CaloEventSetup/interface/CaloTopologyRecord.h"
#include "Geometry/Records/interface/CaloGeometryRecord.h"
#include "EgammaAnalysis/ElectronTools/interface/SuperClusterHelper.h"

#include <iostream>

using namespace edm ;
using namespace std ;
using namespace reco ;

CalibratedElectronProducer::CalibratedElectronProducer( const edm::ParameterSet & cfg )
{
    inputElectrons_ = cfg.getParameter<edm::InputTag>("inputElectronsTag");
    
    nameEnergyReg_      = cfg.getParameter<edm::InputTag>("nameEnergyReg");
    nameEnergyErrorReg_ = cfg.getParameter<edm::InputTag>("nameEnergyErrorReg");
    
    recHitCollectionEB_ = cfg.getParameter<edm::InputTag>("recHitCollectionEB");
    recHitCollectionEE_ = cfg.getParameter<edm::InputTag>("recHitCollectionEE");
    
    
    nameNewEnergyReg_      = cfg.getParameter<std::string>("nameNewEnergyReg");
    nameNewEnergyErrorReg_ = cfg.getParameter<std::string>("nameNewEnergyErrorReg");
    newElectronName_ = cfg.getParameter<std::string>("outputGsfElectronCollectionLabel");
    
    
    dataset = cfg.getParameter<std::string>("inputDataset");
    isMC = cfg.getParameter<bool>("isMC");
    updateEnergyError = cfg.getParameter<bool>("updateEnergyError");
    lumiRatio = cfg.getParameter<double>("lumiRatio");
    correctionsType = cfg.getParameter<int>("correctionsType");
    applyLinearityCorrection = cfg.getParameter<bool>("applyLinearityCorrection");
    combinationType = cfg.getParameter<int>("combinationType");
    verbose = cfg.getParameter<bool>("verbose");
    synchronization = cfg.getParameter<bool>("synchronization");
    combinationRegressionInputPath = cfg.getParameter<std::string>("combinationRegressionInputPath");
    scaleCorrectionsInputPath = cfg.getParameter<std::string>("scaleCorrectionsInputPath");
    linCorrectionsInputPath   = cfg.getParameter<std::string>("linearityCorrectionsInputPath");
  
    //basic checks
    if ( isMC && ( dataset != "Summer11" && dataset != "Fall11" 
        && dataset!= "Summer12" && dataset != "Summer12_DR53X_HCP2012"
        && dataset != "Summer12_LegacyPaper" ) )
    { 
        throw cms::Exception("CalibratedgsfElectronProducer|ConfigError") << "Unknown MC dataset"; 
    }
    if ( !isMC && ( dataset != "Prompt" && dataset != "ReReco" 
        && dataset != "Jan16ReReco" && dataset != "ICHEP2012" 
        && dataset != "Moriond2013" && dataset != "22Jan2013ReReco" ) )
    { 
        throw cms::Exception("CalibratedgsfElectronProducer|ConfigError") << "Unknown Data dataset"; 
    }

    // Linearity correction only applied on combined momentum obtain with regression combination
    if(combinationType!=3 && applyLinearityCorrection)
    {
        std::cout << "[CalibratedElectronProducer] "
            << "Warning: you chose combinationType!=3 and applyLinearityCorrection=True. Linearity corrections are only applied on top of combination 3." << std::endl;
    }
    
    std::cout << "[CalibratedGsfElectronProducer] Correcting scale for dataset " << dataset << std::endl;

    //initializations
    std::string pathToDataCorr;
    switch (correctionsType)
    {
        case 0:
            break;
        case 1: 
            if ( verbose ) 
            {
                std::cout << "You choose regression 1 scale corrections" << std::endl;
            }
            break;
        case 2: 
            if ( verbose ) 
            {
                std::cout << "You choose regression 2 scale corrections." << std::endl;
            }
            break;
        case 3: 
            throw cms::Exception("CalibratedgsfElectronProducer|ConfigError")
                << "You choose standard non-regression ecal energy scale corrections. They are not implemented yet."; 
            break;
        default: 
            throw cms::Exception("CalibratedgsfElectronProducer|ConfigError")
                << "Unknown correctionsType !!!" ;
    }
  
    theEnCorrector = new ElectronEnergyCalibrator
        (
            edm::FileInPath(scaleCorrectionsInputPath.c_str()).fullPath().c_str(), 
            edm::FileInPath(linCorrectionsInputPath.c_str()).fullPath().c_str(),
            dataset, 
            correctionsType, 
            applyLinearityCorrection,
            lumiRatio, 
            isMC, 
            updateEnergyError, 
            verbose, 
            synchronization
        );
    
    if ( verbose ) 
    {
        std::cout<<"[CalibratedGsfElectronProducer] "
        << "ElectronEnergyCalibrator object is created" << std::endl;
    }

    myEpCombinationTool = new EpCombinationTool();
    myEpCombinationTool->init
        (
            edm::FileInPath(combinationRegressionInputPath.c_str()).fullPath().c_str(),
            "CombinationWeight"
        );

    myCombinator = new ElectronEPcombinator();
    
    if ( verbose ) 
    {
        std::cout << "[CalibratedGsfElectronProducer] "
        << "Combination tools are created and initialized" << std::endl;
    }
    
    produces<edm::ValueMap<double> >(nameNewEnergyReg_);
    produces<edm::ValueMap<double> >(nameNewEnergyErrorReg_);
    produces<GsfElectronCollection> (newElectronName_);
    geomInitialized_ = false;
}
 
CalibratedElectronProducer::~CalibratedElectronProducer()
{}

void CalibratedElectronProducer::produce( edm::Event & event, const edm::EventSetup & setup )
{
    if (!geomInitialized_) 
    {
        edm::ESHandle<CaloTopology> theCaloTopology;
        setup.get<CaloTopologyRecord>().get(theCaloTopology);
        ecalTopology_ = & (*theCaloTopology);
        
        edm::ESHandle<CaloGeometry> theCaloGeometry;
        setup.get<CaloGeometryRecord>().get(theCaloGeometry); 
        caloGeometry_ = & (*theCaloGeometry);
        geomInitialized_ = true;
    }

    // Read GsfElectrons
    edm::Handle<reco::GsfElectronCollection>  oldElectronsH ;
    event.getByLabel(inputElectrons_,oldElectronsH) ;
   
    // Read RecHits
    edm::Handle< EcalRecHitCollection > pEBRecHits;
    edm::Handle< EcalRecHitCollection > pEERecHits;
    event.getByLabel( recHitCollectionEB_, pEBRecHits );
    event.getByLabel( recHitCollectionEE_, pEERecHits );

    // ReadValueMaps
    edm::Handle<edm::ValueMap<double> > valMapEnergyH;
    event.getByLabel(nameEnergyReg_,valMapEnergyH);
    edm::Handle<edm::ValueMap<double> > valMapEnergyErrorH;
    event.getByLabel(nameEnergyErrorReg_,valMapEnergyErrorH);
  
    // Prepare output collections
    std::auto_ptr<GsfElectronCollection> electrons( new reco::GsfElectronCollection ) ;
    // Fillers for ValueMaps:
    std::auto_ptr<edm::ValueMap<double> > regrNewEnergyMap(new edm::ValueMap<double>() );
    edm::ValueMap<double>::Filler energyFiller(*regrNewEnergyMap);
    
    std::auto_ptr<edm::ValueMap<double> > regrNewEnergyErrorMap(new edm::ValueMap<double>() );
    edm::ValueMap<double>::Filler energyErrorFiller(*regrNewEnergyErrorMap);
    
    // first clone the initial collection
    unsigned nElectrons = oldElectronsH->size();
    for( unsigned iele = 0; iele < nElectrons; ++iele )
    {    
      electrons->push_back((*oldElectronsH)[iele]);
    }

    std::vector<double> regressionValues;
    std::vector<double> regressionErrorValues;
    regressionValues.reserve(nElectrons);
    regressionErrorValues.reserve(nElectrons);

    if ( correctionsType != 0 )
    {
        for ( unsigned iele = 0; iele < nElectrons ; ++iele) 
        {
            reco::GsfElectron & ele  ( (*electrons)[iele]);
            reco::GsfElectronRef elecRef(oldElectronsH,iele);
            double regressionEnergy = (*valMapEnergyH)[elecRef];
            double regressionEnergyError = (*valMapEnergyErrorH)[elecRef];
            
            regressionValues.push_back(regressionEnergy);
            regressionErrorValues.push_back(regressionEnergyError);

            //    r9 
            const EcalRecHitCollection * recHits=0;
            if( ele.isEB() ) 
            {
                recHits = pEBRecHits.product();
            } else recHits = pEERecHits.product();

            SuperClusterHelper mySCHelper( &(ele), recHits, ecalTopology_, caloGeometry_ );

            int elClass = -1;
            int run = event.run(); 
            
            float r9 = mySCHelper.r9(); 
            double correctedEcalEnergy = ele.correctedEcalEnergy();
            double correctedEcalEnergyError = ele.correctedEcalEnergyError();
            double trackMomentum = ele.trackMomentumAtVtx().R();
            double trackMomentumError = ele.trackMomentumError();
            double combinedMomentum = ele.p();
            double combinedMomentumError = ele.p4Error(ele.candidateP4Kind());
            // FIXME : p4Error not filled for pure tracker electrons
            // Recompute it using the parametrization implemented in 
            // RecoEgamma/EgammaElectronAlgos/src/ElectronEnergyCorrector.cc::simpleParameterizationUncertainty() 
            if( !ele.ecalDrivenSeed() )
            {
                double error = 999. ;
                double momentum = (combinedMomentum<15. ? 15. : combinedMomentum);
                if ( ele.isEB() )
                {
                    float parEB[3] = { 5.24e-02,  2.01e-01, 1.00e-02};
                    error = momentum * sqrt( pow(parEB[0]/sqrt(momentum),2) + pow(parEB[1]/momentum,2) + pow(parEB[2],2) );
                }
                else if ( ele.isEE() )
                {
                    float parEE[3] = { 1.46e-01, 9.21e-01, 1.94e-03} ;
                    error = momentum * sqrt( pow(parEE[0]/sqrt(momentum),2) + pow(parEE[1]/momentum,2) + pow(parEE[2],2) );
                }
                combinedMomentumError = error;
            }
    
            if (ele.classification() == reco::GsfElectron::GOLDEN) {elClass = 0;}
            if (ele.classification() == reco::GsfElectron::BIGBREM) {elClass = 1;}
            if (ele.classification() == reco::GsfElectron::BADTRACK) {elClass = 2;}
            if (ele.classification() == reco::GsfElectron::SHOWERING) {elClass = 3;}
            if (ele.classification() == reco::GsfElectron::GAP) {elClass = 4;}

            SimpleElectron mySimpleElectron
                (
                    run, 
                    elClass, 
                    r9, 
                    correctedEcalEnergy, 
                    correctedEcalEnergyError, 
                    trackMomentum, 
                    trackMomentumError, 
                    regressionEnergy, 
                    regressionEnergyError, 
                    combinedMomentum, 
                    combinedMomentumError, 
                    ele.superCluster()->eta(), 
                    ele.isEB(), 
                    isMC, 
                    ele.ecalDriven(), 
                    ele.trackerDrivenSeed()
                );

            // energy calibration for ecalDriven electrons
            if ( ele.core()->ecalDrivenSeed() || correctionsType==2 || combinationType==3 ) 
            {
                theEnCorrector->calibrate(mySimpleElectron);
        
                // E-p combination  
        
                switch ( combinationType )
                {
                    case 0: 
                        if ( verbose ) 
                        {
                            std::cout << "[CalibratedGsfElectronProducer] " 
                            << "You choose not to combine." << std::endl;
                        }
                        break;
                    case 1: 
                        if ( verbose ) 
                        {
                            std::cout << "[CalibratedGsfElectronProducer] " 
                            << "You choose corrected regression energy for standard combination" << std::endl;
                        }
                        myCombinator->setCombinationMode(1);
                        myCombinator->combine(mySimpleElectron);
                        break;
                    case 2: 
                        if ( verbose ) 
                        {
                            std::cout << "[CalibratedGsfElectronProducer] " 
                            << "You choose uncorrected regression energy for standard combination" << std::endl;
                        }
                        myCombinator->setCombinationMode(2);
                        myCombinator->combine(mySimpleElectron);
                        break;
                    case 3: 
                        if ( verbose ) 
                        {
                            std::cout << "[CalibratedGsfElectronProducer] "
                            << "You choose regression combination." << std::endl;
                        }
                        myEpCombinationTool->combine(mySimpleElectron);
                        theEnCorrector->correctLinearity(mySimpleElectron);
                        break;
                    default: 
                        throw cms::Exception("CalibratedgsfElectronProducer|ConfigError")
                            << "Unknown combination Type !!!" ;
                }

                math::XYZTLorentzVector oldMomentum = ele.p4() ;
                math::XYZTLorentzVector newMomentum_ ;
                newMomentum_ = math::XYZTLorentzVector
                 ( oldMomentum.x()*mySimpleElectron.getCombinedMomentum()/oldMomentum.t(),
                   oldMomentum.y()*mySimpleElectron.getCombinedMomentum()/oldMomentum.t(),
                   oldMomentum.z()*mySimpleElectron.getCombinedMomentum()/oldMomentum.t(),
                   mySimpleElectron.getCombinedMomentum() ) ;

                ele.correctMomentum
                    (
                        newMomentum_,
                        mySimpleElectron.getTrackerMomentumError(),
                        mySimpleElectron.getCombinedMomentumError()
                    );

                if ( verbose ) 
                {
                    std::cout << "[CalibratedGsfElectronProducer] Combined momentum after saving "
                        << ele.p4().t() << std::endl;
                }
            }// end of if (ele.core()->ecalDrivenSeed()) 
        }// end of loop on electrons 
    } else 
    {
        if ( verbose ) 
        {
            std::cout << "[CalibratedGsfElectronProducer] "
            << "You choose not to correct. Uncorrected Regression Energy is taken." << std::endl;
        }
    }

    // Save the electrons
    const edm::OrphanHandle<reco::GsfElectronCollection> gsfNewElectronHandle = event.put(electrons, newElectronName_) ;
    energyFiller.insert(gsfNewElectronHandle,regressionValues.begin(),regressionValues.end());
    energyFiller.fill();
    energyErrorFiller.insert(gsfNewElectronHandle,regressionErrorValues.begin(),regressionErrorValues.end());
    energyErrorFiller.fill();
    
    event.put(regrNewEnergyMap,nameNewEnergyReg_);
    event.put(regrNewEnergyErrorMap,nameNewEnergyErrorReg_);
}

#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/interface/ESProducer.h"
#include "FWCore/Framework/interface/ModuleFactory.h"
DEFINE_FWK_MODULE(CalibratedElectronProducer);