PatUtils/plugins/ShiftedParticleProducer.cc

0001 /** \class ShiftedParticleProducer
0002  *
0003  * Vary energy of electrons/muons/tau-jets by +/- 1 standard deviation,
0004  * in order to estimate resulting uncertainty on MET
0005  *
0006  * NOTE: energy scale uncertainties need to be specified in python config
0007  *
0008  * \author Matthieu marionneau ETH
0009  *
0010  *
0011  *
0012  */
0013
0014 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
0015 #include "DataFormats/Candidate/interface/Candidate.h"
0016 #include "DataFormats/Candidate/interface/LeafCandidate.h"
0017 #include "DataFormats/Common/interface/ValueMap.h"
0018 #include "FWCore/Framework/interface/Event.h"
0019 #include "FWCore/Framework/interface/EventSetup.h"
0020 #include "FWCore/Framework/interface/stream/EDProducer.h"
0021 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0022 #include "FWCore/Utilities/interface/InputTag.h"
0023 #include "FWCore/Utilities/interface/isFinite.h"
0024
0025 #include <TF2.h>
0026
0027 #include <string>
0028 #include <vector>
0029
0030 class ShiftedParticleProducer : public edm::stream::EDProducer<> {
0031   typedef edm::View<reco::Candidate> CandidateView;
0032
0033 public:
0034   explicit ShiftedParticleProducer(const edm::ParameterSet& cfg);
0035   ~ShiftedParticleProducer() override;
0036
0037 private:
0038   void produce(edm::Event& evt, const edm::EventSetup& es) override;
0039
0040   double getUncShift(const reco::Candidate& originalParticle);
0041
0042   std::string moduleLabel_;
0043
0044   edm::EDGetTokenT<CandidateView> srcToken_;
0045   edm::EDGetTokenT<edm::ValueMap<float>> weightsToken_;
0046
0047   struct binningEntryType {
0048     binningEntryType(std::string uncertainty, std::string moduleLabel)
0049         : binSelection_(nullptr), binUncertainty_(uncertainty), energyDep_(false) {
0050       binUncFormula_ =
0051           std::make_unique<TF2>(std::string(moduleLabel).append("_uncFormula").c_str(), binUncertainty_.c_str());
0052     }
0053     binningEntryType(const edm::ParameterSet& cfg, std::string moduleLabel)
0054         : binSelection_(new StringCutObjectSelector<reco::Candidate>(cfg.getParameter<std::string>("binSelection"))),
0055           binUncertainty_(cfg.getParameter<std::string>("binUncertainty")),
0056           energyDep_(false) {
0057       binUncFormula_ =
0058           std::make_unique<TF2>(std::string(moduleLabel).append("_uncFormula").c_str(), binUncertainty_.c_str());
0059       if (cfg.exists("energyDependency")) {
0060         energyDep_ = cfg.getParameter<bool>("energyDependency");
0061       }
0062     }
0063     ~binningEntryType() {}
0064     std::unique_ptr<StringCutObjectSelector<reco::Candidate>> binSelection_;
0065     //double binUncertainty_;
0066     std::string binUncertainty_;
0067     std::unique_ptr<TF2> binUncFormula_;
0068     bool energyDep_;
0069   };
0070   std::vector<binningEntryType*> binning_;
0071
0072   double shiftBy_;  // set to +1.0/-1.0 for up/down variation of energy scale
0073 };
0074
0075 ShiftedParticleProducer::ShiftedParticleProducer(const edm::ParameterSet& cfg) {
0076   moduleLabel_ = cfg.getParameter<std::string>("@module_label");
0077   srcToken_ = consumes<CandidateView>(cfg.getParameter<edm::InputTag>("src"));
0078   shiftBy_ = cfg.getParameter<double>("shiftBy");
0079   edm::InputTag srcWeights = cfg.getParameter<edm::InputTag>("srcWeights");
0080   if (!srcWeights.label().empty())
0081     weightsToken_ = consumes<edm::ValueMap<float>>(srcWeights);
0082
0083   if (cfg.exists("binning")) {
0084     typedef std::vector<edm::ParameterSet> vParameterSet;
0085     vParameterSet cfgBinning = cfg.getParameter<vParameterSet>("binning");
0086     for (vParameterSet::const_iterator cfgBinningEntry = cfgBinning.begin(); cfgBinningEntry != cfgBinning.end();
0087          ++cfgBinningEntry) {
0088       binning_.push_back(new binningEntryType(*cfgBinningEntry, moduleLabel_));
0089     }
0090   } else {
0091     std::string uncertainty = cfg.getParameter<std::string>("uncertainty");
0092     binning_.push_back(new binningEntryType(uncertainty, moduleLabel_));
0093   }
0094
0095   produces<reco::CandidateCollection>();
0096 }
0097
0098 ShiftedParticleProducer::~ShiftedParticleProducer() {
0099   for (std::vector<binningEntryType*>::const_iterator it = binning_.begin(); it != binning_.end(); ++it) {
0100     delete (*it);
0101   }
0102 }
0103
0104 void ShiftedParticleProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
0105   edm::Handle<CandidateView> originalParticles;
0106   evt.getByToken(srcToken_, originalParticles);
0107
0108   edm::Handle<edm::ValueMap<float>> weights;
0109   if (!weightsToken_.isUninitialized())
0110     evt.getByToken(weightsToken_, weights);
0111
0112   auto shiftedParticles = std::make_unique<reco::CandidateCollection>();
0113
0114   for (unsigned i = 0; i < originalParticles->size(); ++i) {
0115     float weight = 1.0;
0116     if (!weightsToken_.isUninitialized()) {
0117       edm::Ptr<reco::Candidate> particlePtr = originalParticles->ptrAt(i);
0118       while (!weights->contains(particlePtr.id()) && (particlePtr->numberOfSourceCandidatePtrs() > 0))
0119         particlePtr = particlePtr->sourceCandidatePtr(0);
0120       weight = (*weights)[particlePtr];
0121     }
0122     const reco::Candidate& originalParticle = originalParticles->at(i);
0123     reco::LeafCandidate weightedParticle(originalParticle);
0124     weightedParticle.setP4(originalParticle.p4() * weight);
0125     double uncertainty = getUncShift(weightedParticle);
0126     double shift = shiftBy_ * uncertainty;
0127
0128     reco::Candidate::LorentzVector shiftedParticleP4 = originalParticle.p4();
0129     //leave 0*nan = 0
0130     if ((weight > 0) && (!(edm::isNotFinite(shift) && shiftedParticleP4.mag2() == 0)))
0131       shiftedParticleP4 *= (1. + shift);
0132
0133     std::unique_ptr<reco::Candidate> shiftedParticle = std::make_unique<reco::LeafCandidate>(originalParticle);
0134     shiftedParticle->setP4(shiftedParticleP4);
0135
0136     shiftedParticles->push_back(shiftedParticle.release());
0137   }
0138
0139   evt.put(std::move(shiftedParticles));
0140 }
0141
0142 double ShiftedParticleProducer::getUncShift(const reco::Candidate& originalParticle) {
0143   double valx = 0;
0144   double valy = 0;
0145   for (std::vector<binningEntryType*>::iterator binningEntry = binning_.begin(); binningEntry != binning_.end();
0146        ++binningEntry) {
0147     if ((!(*binningEntry)->binSelection_) || (*(*binningEntry)->binSelection_)(originalParticle)) {
0148       if ((*binningEntry)->energyDep_)
0149         valx = originalParticle.energy();
0150       else
0151         valx = originalParticle.pt();
0152
0153       valy = originalParticle.eta();
0154       return (*binningEntry)->binUncFormula_->Eval(valx, valy);
0155     }
0156   }
0157   return 0;
0158 }
0159
0160 #include "FWCore/Framework/interface/MakerMacros.h"
0161 DEFINE_FWK_MODULE(ShiftedParticleProducer);