EgammaTools/src/EpCombinationTool.cc

0001 #include "RecoEgamma/EgammaTools/interface/EpCombinationTool.h"
0002
0003 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0004 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0005 #include "FWCore/Framework/interface/EventSetup.h"
0006
0007 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0008 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0009 #include "FWCore/Framework/interface/ConsumesCollector.h"
0010
0011 #include <cmath>
0012 #include <vector>
0013 #include <vdt/vdtMath.h>
0014
0015 EpCombinationTool::EpCombinationTool(const edm::ParameterSet& iConfig, edm::ConsumesCollector&& cc)
0016     : ecalTrkEnergyRegress_(iConfig.getParameter<edm::ParameterSet>("ecalTrkRegressionConfig"), cc),
0017       ecalTrkEnergyRegressUncert_(iConfig.getParameter<edm::ParameterSet>("ecalTrkRegressionUncertConfig"), cc),
0018       maxEcalEnergyForComb_(iConfig.getParameter<double>("maxEcalEnergyForComb")),
0019       minEOverPForComb_(iConfig.getParameter<double>("minEOverPForComb")),
0020       maxEPDiffInSigmaForComb_(iConfig.getParameter<double>("maxEPDiffInSigmaForComb")),
0021       maxRelTrkMomErrForComb_(iConfig.getParameter<double>("maxRelTrkMomErrForComb")) {}
0022
0023 edm::ParameterSetDescription EpCombinationTool::makePSetDescription() {
0024   edm::ParameterSetDescription desc;
0025   desc.add<edm::ParameterSetDescription>("ecalTrkRegressionConfig", EgammaRegressionContainer::makePSetDescription());
0026   desc.add<edm::ParameterSetDescription>("ecalTrkRegressionUncertConfig",
0027                                          EgammaRegressionContainer::makePSetDescription());
0028   desc.add<double>("maxEcalEnergyForComb", 200.);
0029   desc.add<double>("minEOverPForComb", 0.025);
0030   desc.add<double>("maxEPDiffInSigmaForComb", 15.);
0031   desc.add<double>("maxRelTrkMomErrForComb", 10.);
0032   return desc;
0033 }
0034
0035 void EpCombinationTool::setEventContent(const edm::EventSetup& iSetup) {
0036   ecalTrkEnergyRegress_.setEventContent(iSetup);
0037   ecalTrkEnergyRegressUncert_.setEventContent(iSetup);
0038 }
0039
0040 std::pair<float, float> EpCombinationTool::combine(const reco::GsfElectron& ele) const {
0041   return combine(ele, ele.correctedEcalEnergyError());
0042 }
0043
0044 //when doing the E/p combination, its very important to ensure the ecalEnergyErr
0045 //that the regression is trained on is used, not the actual ecalEnergyErr of the electron
0046 //these differ when you correct the ecalEnergyErr by smearing value needed to get data/MC to agree
0047 std::pair<float, float> EpCombinationTool::combine(const reco::GsfElectron& ele, const float corrEcalEnergyErr) const {
0048   const float scRawEnergy = ele.superCluster()->rawEnergy();
0049   const float esEnergy = ele.superCluster()->preshowerEnergy();
0050
0051   const float corrEcalEnergy = ele.correctedEcalEnergy();
0052   const float ecalMean = ele.correctedEcalEnergy() / (scRawEnergy + esEnergy);
0053   const float ecalSigma = corrEcalEnergyErr / corrEcalEnergy;
0054
0055   auto gsfTrk = ele.gsfTrack();
0056
0057   const float trkP = gsfTrk->pMode();
0058   const float trkEta = gsfTrk->etaMode();
0059   const float trkPhi = gsfTrk->phiMode();
0060   const float trkPErr = std::abs(gsfTrk->qoverpModeError()) * trkP * trkP;
0061   const float eOverP = corrEcalEnergy / trkP;
0062   const float fbrem = ele.fbrem();
0063
0064   if (corrEcalEnergy < maxEcalEnergyForComb_ && eOverP > minEOverPForComb_ &&
0065       std::abs(corrEcalEnergy - trkP) <
0066           maxEPDiffInSigmaForComb_ * std::sqrt(trkPErr * trkPErr + corrEcalEnergyErr * corrEcalEnergyErr) &&
0067       trkPErr < maxRelTrkMomErrForComb_ * trkP) {
0068     std::array<float, 9> eval;
0069     eval[0] = corrEcalEnergy;
0070     eval[1] = ecalSigma / ecalMean;
0071     eval[2] = trkPErr / trkP;
0072     eval[3] = eOverP;
0073     eval[4] = ele.ecalDrivenSeed();
0074     eval[5] = ele.full5x5_showerShape().r9;
0075     eval[6] = fbrem;
0076     eval[7] = trkEta;
0077     eval[8] = trkPhi;
0078
0079     const float preCombinationEt = corrEcalEnergy / std::cosh(trkEta);
0080     float mean = ecalTrkEnergyRegress_(preCombinationEt, ele.isEB(), ele.nSaturatedXtals() != 0, eval.data());
0081     float sigma = ecalTrkEnergyRegressUncert_(preCombinationEt, ele.isEB(), ele.nSaturatedXtals() != 0, eval.data());
0082     // Final correction
0083     // A negative energy means that the correction went
0084     // outside the boundaries of the training. In this case uses raw.
0085     // The resolution estimation, on the other hand should be ok.
0086     if (mean < 0.)
0087       mean = 1.0;
0088
0089     //why this differs from the defination of corrEcalEnergyErr (it misses the mean) is not clear to me
0090     //still this is a direct port from EGExtraInfoModifierFromDB, potential bugs and all
0091     const float ecalSigmaTimesRawEnergy = ecalSigma * (scRawEnergy + esEnergy);
0092     const float rawCombEnergy =
0093         (corrEcalEnergy * trkPErr * trkPErr + trkP * ecalSigmaTimesRawEnergy * ecalSigmaTimesRawEnergy) /
0094         (trkPErr * trkPErr + ecalSigmaTimesRawEnergy * ecalSigmaTimesRawEnergy);
0095
0096     return std::make_pair(mean * rawCombEnergy, sigma * rawCombEnergy);
0097   } else {
0098     return std::make_pair(corrEcalEnergy, corrEcalEnergyErr);
0099   }
0100 }