JetProducers/plugins/BetaStarVarProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    PhysicsTools/NanoAOD
0004 // Class:      BetaStarVarProducer
0005 //
0006 /**\class BetaStarVarProducer BetaStarVarProducer.cc PhysicsTools/NanoAOD/plugins/BetaStarVarProducer.cc
0007
0008  Description: This produces value maps to store CHS-related variables for JERC.
0009               This includes the charged hadrons associated to CHS jets,
0010               and those associated to PU that are within the CHS jet.
0011
0012  Implementation:
0013      This uses a ValueMap producer functionality, and
0014      loops over the input candidates (usually PackedCandidates)
0015      that are associated to each jet, counting the candidates associated
0016      to the PV and those not.
0017 */
0018 //
0019 // Original Author:  Sal Rappoccio
0020 // Further editing:  Hannu Siikonen
0021 //
0022 //
0023
0024 #include <memory>
0025
0026 #include "FWCore/Framework/interface/Frameworkfwd.h"
0027 #include "FWCore/Framework/interface/global/EDProducer.h"
0028
0029 #include "FWCore/Framework/interface/Event.h"
0030 #include "FWCore/Framework/interface/MakerMacros.h"
0031
0032 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0033 #include "FWCore/Utilities/interface/StreamID.h"
0034
0035 #include "DataFormats/PatCandidates/interface/Jet.h"
0036 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0037 #include "DataFormats/VertexReco/interface/Vertex.h"
0038
0039 #include "DataFormats/Common/interface/View.h"
0040
0041 template <typename T>
0042 class BetaStarVarProducer : public edm::global::EDProducer<> {
0043 public:
0044   explicit BetaStarVarProducer(const edm::ParameterSet &iConfig)
0045       : srcJet_(consumes<edm::View<pat::Jet>>(iConfig.getParameter<edm::InputTag>("srcJet"))),
0046         srcPF_(consumes<edm::View<T>>(iConfig.getParameter<edm::InputTag>("srcPF"))),
0047         maxDR_(iConfig.getParameter<double>("maxDR")) {
0048     produces<edm::ValueMap<float>>("chargedFromPV0EnergyFraction");
0049     produces<edm::ValueMap<float>>("chargedFromPV1EnergyFraction");
0050     produces<edm::ValueMap<float>>("chargedFromPV2EnergyFraction");
0051     produces<edm::ValueMap<float>>("chargedFromPV3EnergyFraction");
0052   }
0053   ~BetaStarVarProducer() override {}
0054
0055   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0056
0057 private:
0058   void produce(edm::StreamID, edm::Event &, const edm::EventSetup &) const override;
0059
0060   void fillValueMaps(std::string valueName,
0061                      const std::vector<float> &values,
0062                      edm::Handle<edm::View<pat::Jet>> &srcJet,
0063                      edm::Event &iEvent) const;
0064   std::tuple<float, float, float, float> calculateCHSEnergies(edm::Ptr<pat::Jet> const &jet, double chefrompv0) const;
0065
0066   edm::EDGetTokenT<edm::View<pat::Jet>> srcJet_;
0067   edm::EDGetTokenT<edm::View<T>> srcPF_;
0068   double maxDR_;
0069 };
0070
0071 template <typename T>
0072 void BetaStarVarProducer<T>::produce(edm::StreamID streamID, edm::Event &iEvent, const edm::EventSetup &iSetup) const {
0073   edm::Handle<edm::View<pat::Jet>> srcJet;
0074   iEvent.getByToken(srcJet_, srcJet);
0075   edm::Handle<edm::View<T>> srcPF;
0076   iEvent.getByToken(srcPF_, srcPF);
0077
0078   unsigned int nJet = srcJet->size();
0079
0080   // Create an injective mapping from jets to charged PF candidates with fromPV==0 within the jet cone.
0081   // These are the PF candidates supposedly removed by the CHS method - i.e. removed charged pileup.
0082   std::vector<double> jet2pue(nJet, 0);
0083   for (const auto &cand : *srcPF) {
0084     if (cand.charge() == 0 or cand.fromPV() != 0)
0085       continue;
0086     // Looking for the closest match for this charged frompv==0 candidate
0087     float dRMin = 999.;
0088     int bestjet = -1, jetidx = 0;
0089     for (const auto &jet : *srcJet) {
0090       float dR = reco::deltaR(jet, cand);
0091       if (dR < dRMin) {
0092         dRMin = dR;
0093         bestjet = jetidx;
0094       }
0095       ++jetidx;
0096     }
0097     // If the candidate is closer than the jet radius to the jet axis, this is a PU particle of interest for the selected jet
0098     if (dRMin < maxDR_)
0099       jet2pue[bestjet] += cand.energy();
0100   }
0101
0102   std::vector<float> chargedFromPV0EnergyFraction(nJet, -1);
0103   std::vector<float> chargedFromPV1EnergyFraction(nJet, -1);
0104   std::vector<float> chargedFromPV2EnergyFraction(nJet, -1);
0105   std::vector<float> chargedFromPV3EnergyFraction(nJet, -1);
0106
0107   for (unsigned int ij = 0; ij < nJet; ++ij) {
0108     auto vals = calculateCHSEnergies(srcJet->ptrAt(ij), jet2pue[ij]);
0109     chargedFromPV0EnergyFraction[ij] = std::get<0>(vals);
0110     chargedFromPV1EnergyFraction[ij] = std::get<1>(vals);
0111     chargedFromPV2EnergyFraction[ij] = std::get<2>(vals);
0112     chargedFromPV3EnergyFraction[ij] = std::get<3>(vals);
0113   }
0114
0115   fillValueMaps("chargedFromPV0EnergyFraction", chargedFromPV0EnergyFraction, srcJet, iEvent);
0116   fillValueMaps("chargedFromPV1EnergyFraction", chargedFromPV1EnergyFraction, srcJet, iEvent);
0117   fillValueMaps("chargedFromPV2EnergyFraction", chargedFromPV2EnergyFraction, srcJet, iEvent);
0118   fillValueMaps("chargedFromPV3EnergyFraction", chargedFromPV3EnergyFraction, srcJet, iEvent);
0119 }
0120
0121 template <typename T>
0122 void BetaStarVarProducer<T>::fillValueMaps(std::string valueName,
0123                                            const std::vector<float> &values,
0124                                            edm::Handle<edm::View<pat::Jet>> &srcJet,
0125                                            edm::Event &iEvent) const {
0126   std::unique_ptr<edm::ValueMap<float>> valuesPtr(new edm::ValueMap<float>());
0127   edm::ValueMap<float>::Filler valuesFiller(*valuesPtr);
0128   valuesFiller.insert(srcJet, values.begin(), values.end());
0129   valuesFiller.fill();
0130   iEvent.put(std::move(valuesPtr), valueName);
0131 }
0132
0133 template <typename T>
0134 std::tuple<float, float, float, float> BetaStarVarProducer<T>::calculateCHSEnergies(edm::Ptr<pat::Jet> const &ijet,
0135                                                                                     double chefrompv0) const {
0136   // Keep track of energy (present in the jet) for PU stuff
0137   double chefrompv1 = 0.0;
0138   double chefrompv2 = 0.0;
0139   double chefrompv3 = 0.0;
0140
0141   // Loop through the PF candidates within the jet.
0142   // Store the sum of their energy, and their indices.
0143   for (auto pidx = 0u; pidx < ijet->numberOfDaughters(); ++pidx) {
0144     auto *dtr = dynamic_cast<const pat::PackedCandidate *>(ijet->daughter(pidx));
0145     if (dtr->charge() == 0)
0146       continue;
0147     if (dtr->fromPV() == 1)
0148       chefrompv1 += dtr->energy();
0149     else if (dtr->fromPV() == 2)
0150       chefrompv2 += dtr->energy();
0151     else if (dtr->fromPV() == 3)
0152       chefrompv3 += dtr->energy();
0153   }
0154
0155   // Now get the fractions relative to the raw jet.
0156   auto rawP4 = ijet->correctedP4(0);
0157   double chffpv0 = chefrompv0 / rawP4.energy();
0158   double chffpv1 = chefrompv1 / rawP4.energy();
0159   double chffpv2 = chefrompv2 / rawP4.energy();
0160   double chffpv3 = chefrompv3 / rawP4.energy();
0161
0162   return std::tuple<float, float, float, float>(chffpv0, chffpv1, chffpv2, chffpv3);
0163 }
0164
0165 template <typename T>
0166 void BetaStarVarProducer<T>::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0167   edm::ParameterSetDescription desc;
0168   desc.add<edm::InputTag>("srcJet")->setComment("jet input collection");
0169   desc.add<edm::InputTag>("srcPF")->setComment("PF candidate input collection");
0170   desc.add<double>("maxDR")->setComment("Maximum DR to consider for jet-->pf cand association");
0171   std::string modname("BetaStarVarProducer");
0172   descriptions.add(modname, desc);
0173 }
0174
0175 typedef BetaStarVarProducer<pat::PackedCandidate> BetaStarPackedCandidateVarProducer;
0176
0177 DEFINE_FWK_MODULE(BetaStarPackedCandidateVarProducer);