EGEtScaleSysModifier.cc

Go to the documentation of this file.

#include "CommonTools/CandAlgos/interface/ModifyObjectValueBase.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "RecoEgamma/EgammaTools/interface/EpCombinationTool.h"
#include "RecoEgamma/EgammaTools/interface/EGEnergySysIndex.h"

#include <vdt/vdtMath.h>

//in the legacy re-reco we came across an Et scale issue, specifically an inflection at 45 GeV
//it is problematic to address in the normal scale and smearing code therefore
//this patch modifies the e/gamma object to "patch in" this additional systematic
//Questions:
//1 ) Why dont we add this to the overall scale systematic?
//Well we could but the whole point of this is to get a proper template which can be evolved
//correctly. The issue is not that the current systematic doesnt cover this systematic on a
//per bin basis, it does, the problem is the sign flips at 45 GeV so its fine if you use
//only eles/phos > 45 or <45 GeV but the template cant simultaneously cover the entire spectrum
//and you'll get a nasty constrained fit.
//And if you care about these sort of things (and you probably should), you shouldnt be using
//the overall systematic anyways but its seperate parts
//
//2 ) Why dont you do it inside the standard class
//We could but we're rather hoping to solve this issue more cleanly in the future
//but we would hold up the reminiAOD too long to do so so this is just to get something
//which should be fine for 90% of analyses in the miniAOD. So this is a temporary fix
//which we hope to rip out soon (if you're reading this in 2024 because we're still doing it
//this way, all I can say is sorry, we thought it would go away soon! )

class EGEtScaleSysModifier : public ModifyObjectValueBase {
public:
  EGEtScaleSysModifier(const edm::ParameterSet& conf, edm::ConsumesCollector&);
  ~EGEtScaleSysModifier() override {}

  void setEvent(const edm::Event&) final;
  void setEventContent(const edm::EventSetup&) final;

  void modifyObject(pat::Electron& ele) const final;
  void modifyObject(pat::Photon& pho) const final;

private:
  std::pair<float, float> calCombinedMom(reco::GsfElectron& ele, const float scale, const float smear) const;
  void setEcalEnergy(reco::GsfElectron& ele, const float scale, const float smear) const;

  class UncertFuncBase {
  public:
    UncertFuncBase() {}
    virtual ~UncertFuncBase() {}
    virtual float val(const float et) const = 0;
  };

  //defines two uncertaintes, one Et<X and one Et>Y
  //for X<Et<Y, it linearly extrapolates betwen the two values
  class UncertFuncV1 : public UncertFuncBase {
  public:
    UncertFuncV1(const edm::ParameterSet& conf)
        : lowEt_(conf.getParameter<double>("lowEt")),
          highEt_(conf.getParameter<double>("highEt")),
          lowEtUncert_(conf.getParameter<double>("lowEtUncert")),
          highEtUncert_(conf.getParameter<double>("highEtUncert")),
          dEt_(highEt_ - lowEt_),
          dUncert_(highEtUncert_ - lowEtUncert_) {
      if (highEt_ <= lowEt_)
        throw cms::Exception("ConfigError") << " highEt " << highEt_ << " is not higher than lowEt " << lowEt_;
    }
    ~UncertFuncV1() override {}

    float val(const float et) const override {
      if (et <= lowEt_)
        return lowEtUncert_;
      else if (et >= highEt_)
        return highEtUncert_;
      else {
        return (et - lowEt_) * dUncert_ / dEt_ + lowEtUncert_;
      }
    }

  private:
    float lowEt_;
    float highEt_;
    float lowEtUncert_;
    float highEtUncert_;
    float dEt_;
    float dUncert_;
  };

  EpCombinationTool epCombTool_;
  std::unique_ptr<UncertFuncBase> uncertFunc_;
};

EGEtScaleSysModifier::EGEtScaleSysModifier(const edm::ParameterSet& conf, edm::ConsumesCollector&)
    : ModifyObjectValueBase(conf), epCombTool_(conf.getParameter<edm::ParameterSet>("epCombConfig")) {
  const edm::ParameterSet funcPSet = conf.getParameter<edm::ParameterSet>("uncertFunc");
  const std::string& funcName = funcPSet.getParameter<std::string>("name");
  if (funcName == "UncertFuncV1") {
    uncertFunc_ = std::make_unique<UncertFuncV1>(funcPSet);
  } else {
    throw cms::Exception("ConfigError") << "Error constructing EGEtScaleSysModifier, function name " << funcName
                                        << " not valid";
  }
}

void EGEtScaleSysModifier::setEvent(const edm::Event& iEvent) {}

void EGEtScaleSysModifier::setEventContent(const edm::EventSetup& iSetup) { epCombTool_.setEventContent(iSetup); }

void EGEtScaleSysModifier::modifyObject(pat::Electron& ele) const {
  auto getVal = [](const pat::Electron& ele, EGEnergySysIndex::Index valIndex) {
    return ele.userFloat(EGEnergySysIndex::name(valIndex));
  };
  //so ele.energy() may be either pre or post corrections, we have no idea
  //so we explicity access the pre and post correction ecal energies
  //we need the pre corrected to properly do the e/p combination
  //we need the post corrected to get et uncertainty
  const float ecalEnergyPostCorr = getVal(ele, EGEnergySysIndex::kEcalPostCorr);
  const float ecalEnergyPreCorr = getVal(ele, EGEnergySysIndex::kEcalPreCorr);
  const float ecalEnergyErrPreCorr = getVal(ele, EGEnergySysIndex::kEcalErrPreCorr);

  //the et cut is in terms of ecal et using the track angle and post corr ecal energy
  const float etUncert = uncertFunc_->val(ele.et() / ele.energy() * ecalEnergyPostCorr);
  const float smear = getVal(ele, EGEnergySysIndex::kSmearValue);
  const float corr = getVal(ele, EGEnergySysIndex::kScaleValue);

  //get the values we have to reset back to
  const float oldEcalEnergy = ele.ecalEnergy();
  const float oldEcalEnergyErr = ele.ecalEnergyError();
  const auto oldP4 = ele.p4();
  const float oldP4Err = ele.p4Error(reco::GsfElectron::P4_COMBINATION);
  const float oldTrkMomErr = ele.trackMomentumError();

  ele.setCorrectedEcalEnergy(ecalEnergyPreCorr);
  ele.setCorrectedEcalEnergyError(ecalEnergyErrPreCorr);

  const float energyEtUncertUp = calCombinedMom(ele, corr + etUncert, smear).first;
  const float energyEtUncertDn = calCombinedMom(ele, corr - etUncert, smear).first;

  //reset it back to how it was
  ele.setCorrectedEcalEnergy(oldEcalEnergy);
  ele.setCorrectedEcalEnergyError(oldEcalEnergyErr);
  ele.correctMomentum(oldP4, oldTrkMomErr, oldP4Err);

  ele.addUserFloat("energyScaleEtUp", energyEtUncertUp);
  ele.addUserFloat("energyScaleEtDown", energyEtUncertDn);
}

void EGEtScaleSysModifier::modifyObject(pat::Photon& pho) const {
  auto getVal = [](const pat::Photon& pho, EGEnergySysIndex::Index valIndex) {
    return pho.userFloat(EGEnergySysIndex::name(valIndex));
  };
  //so pho.energy() may be either pre or post corrections, we have no idea
  //so we explicity access the pre and post correction ecal energies
  //post corr for the et value for the systematic, pre corr to apply them
  const float ecalEnergyPostCorr = getVal(pho, EGEnergySysIndex::kEcalPostCorr);
  const float ecalEnergyPreCorr = getVal(pho, EGEnergySysIndex::kEcalPreCorr);

  //the et cut is in terms of post corr ecal energy
  const float etUncert = uncertFunc_->val(pho.et() / pho.energy() * ecalEnergyPostCorr);
  const float corr = getVal(pho, EGEnergySysIndex::kScaleValue);

  const float energyEtUncertUp = ecalEnergyPreCorr * (corr + etUncert);
  const float energyEtUncertDn = ecalEnergyPreCorr * (corr - etUncert);

  pho.addUserFloat("energyScaleEtUp", energyEtUncertUp);
  pho.addUserFloat("energyScaleEtDown", energyEtUncertDn);
}

std::pair<float, float> EGEtScaleSysModifier::calCombinedMom(reco::GsfElectron& ele,
                                                             const float scale,
                                                             const float smear) const {
  const float oldEcalEnergy = ele.ecalEnergy();
  const float oldEcalEnergyErr = ele.ecalEnergyError();
  const auto oldP4 = ele.p4();
  const float oldP4Err = ele.p4Error(reco::GsfElectron::P4_COMBINATION);
  const float oldTrkMomErr = ele.trackMomentumError();

  setEcalEnergy(ele, scale, smear);
  const auto& combinedMomentum = epCombTool_.combine(ele);
  ele.setCorrectedEcalEnergy(oldEcalEnergy);
  ele.setCorrectedEcalEnergyError(oldEcalEnergyErr);
  ele.correctMomentum(oldP4, oldTrkMomErr, oldP4Err);

  return combinedMomentum;
}

void EGEtScaleSysModifier::setEcalEnergy(reco::GsfElectron& ele, const float scale, const float smear) const {
  const float oldEcalEnergy = ele.ecalEnergy();
  const float oldEcalEnergyErr = ele.ecalEnergyError();
  ele.setCorrectedEcalEnergy(oldEcalEnergy * scale);
  ele.setCorrectedEcalEnergyError(std::hypot(oldEcalEnergyErr * scale, oldEcalEnergy * smear * scale));
}

DEFINE_EDM_PLUGIN(ModifyObjectValueFactory, EGEtScaleSysModifier, "EGEtScaleSysModifier");