EgammaTools/plugins/EGEtScaleSysModifier.cc

0001 #include "CommonTools/CandAlgos/interface/ModifyObjectValueBase.h"
0002 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0003 #include "RecoEgamma/EgammaTools/interface/EpCombinationTool.h"
0004 #include "RecoEgamma/EgammaTools/interface/EGEnergySysIndex.h"
0005
0006 #include <vdt/vdtMath.h>
0007
0008 //in the legacy re-reco we came across an Et scale issue, specifically an inflection at 45 GeV
0009 //it is problematic to address in the normal scale and smearing code therefore
0010 //this patch modifies the e/gamma object to "patch in" this additional systematic
0011 //Questions:
0012 //1 ) Why dont we add this to the overall scale systematic?
0013 //Well we could but the whole point of this is to get a proper template which can be evolved
0014 //correctly. The issue is not that the current systematic doesnt cover this systematic on a
0015 //per bin basis, it does, the problem is the sign flips at 45 GeV so its fine if you use
0016 //only eles/phos > 45 or <45 GeV but the template cant simultaneously cover the entire spectrum
0017 //and you'll get a nasty constrained fit.
0018 //And if you care about these sort of things (and you probably should), you shouldnt be using
0019 //the overall systematic anyways but its seperate parts
0020 //
0021 //2 ) Why dont you do it inside the standard class
0022 //We could but we're rather hoping to solve this issue more cleanly in the future
0023 //but we would hold up the reminiAOD too long to do so so this is just to get something
0024 //which should be fine for 90% of analyses in the miniAOD. So this is a temporary fix
0025 //which we hope to rip out soon (if you're reading this in 2024 because we're still doing it
0026 //this way, all I can say is sorry, we thought it would go away soon! )
0027
0028 class EGEtScaleSysModifier : public ModifyObjectValueBase {
0029 public:
0030   EGEtScaleSysModifier(const edm::ParameterSet& conf, edm::ConsumesCollector&);
0031   ~EGEtScaleSysModifier() override {}
0032
0033   void setEvent(const edm::Event&) final;
0034   void setEventContent(const edm::EventSetup&) final;
0035
0036   void modifyObject(pat::Electron& ele) const final;
0037   void modifyObject(pat::Photon& pho) const final;
0038
0039 private:
0040   std::pair<float, float> calCombinedMom(reco::GsfElectron& ele, const float scale, const float smear) const;
0041   void setEcalEnergy(reco::GsfElectron& ele, const float scale, const float smear) const;
0042
0043   class UncertFuncBase {
0044   public:
0045     UncertFuncBase() {}
0046     virtual ~UncertFuncBase() {}
0047     virtual float val(const float et) const = 0;
0048   };
0049
0050   //defines two uncertaintes, one Et<X and one Et>Y
0051   //for X<Et<Y, it linearly extrapolates betwen the two values
0052   class UncertFuncV1 : public UncertFuncBase {
0053   public:
0054     UncertFuncV1(const edm::ParameterSet& conf)
0055         : lowEt_(conf.getParameter<double>("lowEt")),
0056           highEt_(conf.getParameter<double>("highEt")),
0057           lowEtUncert_(conf.getParameter<double>("lowEtUncert")),
0058           highEtUncert_(conf.getParameter<double>("highEtUncert")),
0059           dEt_(highEt_ - lowEt_),
0060           dUncert_(highEtUncert_ - lowEtUncert_) {
0061       if (highEt_ <= lowEt_)
0062         throw cms::Exception("ConfigError") << " highEt " << highEt_ << " is not higher than lowEt " << lowEt_;
0063     }
0064     ~UncertFuncV1() override {}
0065
0066     float val(const float et) const override {
0067       if (et <= lowEt_)
0068         return lowEtUncert_;
0069       else if (et >= highEt_)
0070         return highEtUncert_;
0071       else {
0072         return (et - lowEt_) * dUncert_ / dEt_ + lowEtUncert_;
0073       }
0074     }
0075
0076   private:
0077     float lowEt_;
0078     float highEt_;
0079     float lowEtUncert_;
0080     float highEtUncert_;
0081     float dEt_;
0082     float dUncert_;
0083   };
0084
0085   EpCombinationTool epCombTool_;
0086   std::unique_ptr<UncertFuncBase> uncertFunc_;
0087 };
0088
0089 EGEtScaleSysModifier::EGEtScaleSysModifier(const edm::ParameterSet& conf, edm::ConsumesCollector& cc)
0090     : ModifyObjectValueBase(conf), epCombTool_{conf.getParameterSet("epCombConfig"), std::move(cc)} {
0091   const edm::ParameterSet& funcPSet = conf.getParameterSet("uncertFunc");
0092   const std::string& funcName = funcPSet.getParameter<std::string>("name");
0093   if (funcName == "UncertFuncV1") {
0094     uncertFunc_ = std::make_unique<UncertFuncV1>(funcPSet);
0095   } else {
0096     throw cms::Exception("ConfigError") << "Error constructing EGEtScaleSysModifier, function name " << funcName
0097                                         << " not valid";
0098   }
0099 }
0100
0101 void EGEtScaleSysModifier::setEvent(const edm::Event& iEvent) {}
0102
0103 void EGEtScaleSysModifier::setEventContent(const edm::EventSetup& iSetup) { epCombTool_.setEventContent(iSetup); }
0104
0105 void EGEtScaleSysModifier::modifyObject(pat::Electron& ele) const {
0106   auto getVal = [](const pat::Electron& ele, EGEnergySysIndex::Index valIndex) {
0107     return ele.userFloat(EGEnergySysIndex::name(valIndex));
0108   };
0109   //so ele.energy() may be either pre or post corrections, we have no idea
0110   //so we explicity access the pre and post correction ecal energies
0111   //we need the pre corrected to properly do the e/p combination
0112   //we need the post corrected to get et uncertainty
0113   const float ecalEnergyPostCorr = getVal(ele, EGEnergySysIndex::kEcalPostCorr);
0114   const float ecalEnergyPreCorr = getVal(ele, EGEnergySysIndex::kEcalPreCorr);
0115   const float ecalEnergyErrPreCorr = getVal(ele, EGEnergySysIndex::kEcalErrPreCorr);
0116
0117   //the et cut is in terms of ecal et using the track angle and post corr ecal energy
0118   const float etUncert = uncertFunc_->val(ele.et() / ele.energy() * ecalEnergyPostCorr);
0119   const float smear = getVal(ele, EGEnergySysIndex::kSmearValue);
0120   const float corr = getVal(ele, EGEnergySysIndex::kScaleValue);
0121
0122   //get the values we have to reset back to
0123   const float oldEcalEnergy = ele.ecalEnergy();
0124   const float oldEcalEnergyErr = ele.ecalEnergyError();
0125   const auto oldP4 = ele.p4();
0126   const float oldP4Err = ele.p4Error(reco::GsfElectron::P4_COMBINATION);
0127   const float oldTrkMomErr = ele.trackMomentumError();
0128
0129   ele.setCorrectedEcalEnergy(ecalEnergyPreCorr);
0130   ele.setCorrectedEcalEnergyError(ecalEnergyErrPreCorr);
0131
0132   const float energyEtUncertUp = calCombinedMom(ele, corr + etUncert, smear).first;
0133   const float energyEtUncertDn = calCombinedMom(ele, corr - etUncert, smear).first;
0134
0135   //reset it back to how it was
0136   ele.setCorrectedEcalEnergy(oldEcalEnergy);
0137   ele.setCorrectedEcalEnergyError(oldEcalEnergyErr);
0138   ele.correctMomentum(oldP4, oldTrkMomErr, oldP4Err);
0139
0140   ele.addUserFloat("energyScaleEtUp", energyEtUncertUp);
0141   ele.addUserFloat("energyScaleEtDown", energyEtUncertDn);
0142 }
0143
0144 void EGEtScaleSysModifier::modifyObject(pat::Photon& pho) const {
0145   auto getVal = [](const pat::Photon& pho, EGEnergySysIndex::Index valIndex) {
0146     return pho.userFloat(EGEnergySysIndex::name(valIndex));
0147   };
0148   //so pho.energy() may be either pre or post corrections, we have no idea
0149   //so we explicity access the pre and post correction ecal energies
0150   //post corr for the et value for the systematic, pre corr to apply them
0151   const float ecalEnergyPostCorr = getVal(pho, EGEnergySysIndex::kEcalPostCorr);
0152   const float ecalEnergyPreCorr = getVal(pho, EGEnergySysIndex::kEcalPreCorr);
0153
0154   //the et cut is in terms of post corr ecal energy
0155   const float etUncert = uncertFunc_->val(pho.et() / pho.energy() * ecalEnergyPostCorr);
0156   const float corr = getVal(pho, EGEnergySysIndex::kScaleValue);
0157
0158   const float energyEtUncertUp = ecalEnergyPreCorr * (corr + etUncert);
0159   const float energyEtUncertDn = ecalEnergyPreCorr * (corr - etUncert);
0160
0161   pho.addUserFloat("energyScaleEtUp", energyEtUncertUp);
0162   pho.addUserFloat("energyScaleEtDown", energyEtUncertDn);
0163 }
0164
0165 std::pair<float, float> EGEtScaleSysModifier::calCombinedMom(reco::GsfElectron& ele,
0166                                                              const float scale,
0167                                                              const float smear) const {
0168   const float oldEcalEnergy = ele.ecalEnergy();
0169   const float oldEcalEnergyErr = ele.ecalEnergyError();
0170   const auto oldP4 = ele.p4();
0171   const float oldP4Err = ele.p4Error(reco::GsfElectron::P4_COMBINATION);
0172   const float oldTrkMomErr = ele.trackMomentumError();
0173
0174   setEcalEnergy(ele, scale, smear);
0175   const auto& combinedMomentum = epCombTool_.combine(ele);
0176   ele.setCorrectedEcalEnergy(oldEcalEnergy);
0177   ele.setCorrectedEcalEnergyError(oldEcalEnergyErr);
0178   ele.correctMomentum(oldP4, oldTrkMomErr, oldP4Err);
0179
0180   return combinedMomentum;
0181 }
0182
0183 void EGEtScaleSysModifier::setEcalEnergy(reco::GsfElectron& ele, const float scale, const float smear) const {
0184   const float oldEcalEnergy = ele.ecalEnergy();
0185   const float oldEcalEnergyErr = ele.ecalEnergyError();
0186   ele.setCorrectedEcalEnergy(oldEcalEnergy * scale);
0187   ele.setCorrectedEcalEnergyError(std::hypot(oldEcalEnergyErr * scale, oldEcalEnergy * smear * scale));
0188 }
0189
0190 DEFINE_EDM_PLUGIN(ModifyObjectValueFactory, EGEtScaleSysModifier, "EGEtScaleSysModifier");