JetProducers/plugins/HFJetShowerShape.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    HFJetShowerShape/HFJetShowerShape
0004 // Class:      HFJetShowerShape
0005 //
0006 /**\class HFJetShowerShape HFJetShowerShape.cc HFJetShowerShape/HFJetShowerShape/plugins/HFJetShowerShape.cc
0007
0008  Description: [one line class summary]
0009
0010  Implementation:
0011      [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  Laurent Thomas
0015 //         Created:  Tue, 25 Aug 2020 09:17:42 GMT
0016 //
0017 //
0018
0019 // system include files
0020 #include <memory>
0021
0022 // user include files
0023 #include "FWCore/Framework/interface/Frameworkfwd.h"
0024 #include "FWCore/Framework/interface/stream/EDProducer.h"
0025
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/Framework/interface/MakerMacros.h"
0028
0029 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0030 #include "FWCore/Utilities/interface/StreamID.h"
0031
0032 #include "DataFormats/JetReco/interface/PFJet.h"
0033 #include "DataFormats/VertexReco/interface/Vertex.h"
0034
0035 #include "DataFormats/Math/interface/deltaPhi.h"
0036
0037 #include <iostream>
0038
0039 //
0040 // class declaration
0041 //
0042
0043 class HFJetShowerShape : public edm::stream::EDProducer<> {
0044 public:
0045   explicit HFJetShowerShape(const edm::ParameterSet&);
0046
0047   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0048
0049 private:
0050   void produce(edm::Event&, const edm::EventSetup&) override;
0051   template <typename T>
0052   void putInEvent(const std::string&, const edm::Handle<edm::View<reco::Jet>>&, std::vector<T>, edm::Event&) const;
0053
0054   const edm::EDGetTokenT<edm::View<reco::Jet>> jets_token_;
0055   const edm::EDGetTokenT<std::vector<reco::Vertex>> vertices_token_;
0056
0057   //Jet pt/eta thresholds
0058   const double jetPtThreshold_, jetEtaThreshold_;
0059   //HF geometry
0060   const double hfTowerEtaWidth_, hfTowerPhiWidth_;
0061   //Variables for PU subtraction
0062   const double vertexRecoEffcy_, offsetPerPU_, jetReferenceRadius_;
0063   //Pt thresholds for showershape variable calculation
0064   const double stripPtThreshold_, widthPtThreshold_;
0065 };
0066
0067 //
0068 // constants, enums and typedefs
0069 //
0070
0071 //
0072 // static data member definitions
0073 //
0074
0075 //
0076 // constructors and destructor
0077 //
0078 HFJetShowerShape::HFJetShowerShape(const edm::ParameterSet& iConfig)
0079     : jets_token_(consumes<edm::View<reco::Jet>>(iConfig.getParameter<edm::InputTag>("jets"))),
0080       vertices_token_(consumes<std::vector<reco::Vertex>>(iConfig.getParameter<edm::InputTag>("vertices"))),
0081       jetPtThreshold_(iConfig.getParameter<double>("jetPtThreshold")),
0082       jetEtaThreshold_(iConfig.getParameter<double>("jetEtaThreshold")),
0083       hfTowerEtaWidth_(iConfig.getParameter<double>("hfTowerEtaWidth")),
0084       hfTowerPhiWidth_(iConfig.getParameter<double>("hfTowerPhiWidth")),
0085       vertexRecoEffcy_(iConfig.getParameter<double>("vertexRecoEffcy")),
0086       offsetPerPU_(iConfig.getParameter<double>("offsetPerPU")),
0087       jetReferenceRadius_(iConfig.getParameter<double>("jetReferenceRadius")),
0088       stripPtThreshold_(iConfig.getParameter<double>("stripPtThreshold")),
0089       widthPtThreshold_(iConfig.getParameter<double>("widthPtThreshold")) {
0090   produces<edm::ValueMap<float>>("sigmaEtaEta");
0091   produces<edm::ValueMap<float>>("sigmaPhiPhi");
0092   produces<edm::ValueMap<int>>("centralEtaStripSize");
0093   produces<edm::ValueMap<int>>("adjacentEtaStripsSize");
0094 }
0095
0096 //
0097 // member functions
0098 //
0099
0100 // ------------ method called to produce the data  ------------
0101 void HFJetShowerShape::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0102   using namespace edm;
0103
0104   //Jets
0105   auto theJets = iEvent.getHandle(jets_token_);
0106
0107   //Vertices
0108   int nPV = iEvent.get(vertices_token_).size();
0109
0110   //Products
0111   std::vector<float> v_sigmaEtaEta, v_sigmaPhiPhi;
0112   v_sigmaEtaEta.reserve(theJets->size());
0113   v_sigmaPhiPhi.reserve(theJets->size());
0114   std::vector<int> v_size_CentralEtaStrip, v_size_AdjacentEtaStrips;
0115   v_size_CentralEtaStrip.reserve(theJets->size());
0116   v_size_AdjacentEtaStrips.reserve(theJets->size());
0117
0118   //Et offset for HF PF candidates
0119   double puoffset = offsetPerPU_ / (M_PI * jetReferenceRadius_ * jetReferenceRadius_) * nPV / vertexRecoEffcy_ *
0120                     (hfTowerEtaWidth_ * hfTowerPhiWidth_);
0121
0122   for (auto const& jet : *theJets) {
0123     double pt_jet = jet.pt();
0124     double eta_jet = jet.eta();
0125
0126     //If central or low pt jets, fill with dummy variables
0127     if (pt_jet <= jetPtThreshold_ || std::abs(eta_jet) <= jetEtaThreshold_) {
0128       v_sigmaEtaEta.push_back(-1.);
0129       v_sigmaPhiPhi.push_back(-1.);
0130       v_size_CentralEtaStrip.push_back(0);
0131       v_size_AdjacentEtaStrips.push_back(0);
0132     } else {
0133       //First loop over PF candidates to compute some global variables needed for shower shape calculations
0134       double sumptPFcands = 0.;
0135
0136       for (unsigned i = 0; i < jet.numberOfSourceCandidatePtrs(); ++i) {
0137         const reco::Candidate* icand = jet.sourceCandidatePtr(i).get();
0138         //Only look at pdgId =1,2 (HF PF cands)
0139         if (std::abs(icand->pdgId()) > 2)
0140           continue;
0141         double pt_PUsub = icand->pt() - puoffset;
0142         if (pt_PUsub < widthPtThreshold_)
0143           continue;
0144         sumptPFcands += pt_PUsub;
0145       }
0146
0147       //Second loop to compute the various shower shape variables
0148       int size_CentralEtaStrip(0.), size_AdjacentEtaStrips(0.);
0149       double sigmaEtaEtaSq(0.), sigmaPhiPhiSq(0.);
0150       double sumweightsPFcands = 0;
0151       for (unsigned i = 0; i < jet.numberOfSourceCandidatePtrs(); ++i) {
0152         const reco::Candidate* icand = jet.sourceCandidatePtr(i).get();
0153         //Only look at pdgId =1,2 (HF PF cands)
0154         if (std::abs(icand->pdgId()) > 2)
0155           continue;
0156
0157         double deta = std::abs(icand->eta() - jet.eta());
0158         double dphi = std::abs(deltaPhi(icand->phi(), jet.phi()));
0159         double pt_PUsub = icand->pt() - puoffset;
0160
0161         //This is simply the size of the central eta strip and the adjacent strips
0162         if (pt_PUsub >= stripPtThreshold_) {
0163           if (dphi <= hfTowerPhiWidth_ * 0.5)
0164             size_CentralEtaStrip++;
0165           else if (dphi <= hfTowerPhiWidth_ * 1.5)
0166             size_AdjacentEtaStrips++;
0167         }
0168
0169         //Now computing sigmaEtaEta/PhiPhi
0170         if (pt_PUsub >= widthPtThreshold_ && sumptPFcands > 0) {
0171           double weight = pt_PUsub / sumptPFcands;
0172           sigmaEtaEtaSq += deta * deta * weight;
0173           sigmaPhiPhiSq += dphi * dphi * weight;
0174           sumweightsPFcands += weight;
0175         }
0176       }
0177
0178       v_size_CentralEtaStrip.push_back(size_CentralEtaStrip);
0179       v_size_AdjacentEtaStrips.push_back(size_AdjacentEtaStrips);
0180
0181       if (sumweightsPFcands > 0 && sigmaEtaEtaSq > 0 && sigmaPhiPhiSq > 0) {
0182         v_sigmaEtaEta.push_back(sqrt(sigmaEtaEtaSq / sumweightsPFcands));
0183         v_sigmaPhiPhi.push_back(sqrt(sigmaPhiPhiSq / sumweightsPFcands));
0184       } else {
0185         v_sigmaEtaEta.push_back(-1.);
0186         v_sigmaPhiPhi.push_back(-1.);
0187       }
0188
0189     }  //End loop over jets
0190   }
0191
0192   putInEvent("sigmaEtaEta", theJets, v_sigmaEtaEta, iEvent);
0193   putInEvent("sigmaPhiPhi", theJets, v_sigmaPhiPhi, iEvent);
0194   putInEvent("centralEtaStripSize", theJets, v_size_CentralEtaStrip, iEvent);
0195   putInEvent("adjacentEtaStripsSize", theJets, v_size_AdjacentEtaStrips, iEvent);
0196 }
0197
0198 /// Function to put product into event
0199 template <typename T>
0200 void HFJetShowerShape::putInEvent(const std::string& name,
0201                                   const edm::Handle<edm::View<reco::Jet>>& jets,
0202                                   std::vector<T> product,
0203                                   edm::Event& iEvent) const {
0204   auto out = std::make_unique<edm::ValueMap<T>>();
0205   typename edm::ValueMap<T>::Filler filler(*out);
0206   filler.insert(jets, product.begin(), product.end());
0207   filler.fill();
0208   iEvent.put(std::move(out), name);
0209 }
0210
0211 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0212 void HFJetShowerShape::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0213   edm::ParameterSetDescription desc;
0214   desc.add<edm::InputTag>("jets", edm::InputTag("ak4PFJetsCHS"));
0215   desc.add<edm::InputTag>("vertices", edm::InputTag("offlineSlimmedPrimaryVertices"));
0216   desc.add<double>("jetPtThreshold", 25.);
0217   desc.add<double>("jetEtaThreshold", 2.9);
0218   desc.add<double>("hfTowerEtaWidth", 0.175);
0219   desc.add<double>("hfTowerPhiWidth", 0.175);
0220   desc.add<double>("vertexRecoEffcy", 0.7);
0221   desc.add<double>("offsetPerPU", 0.4);
0222   desc.add<double>("jetReferenceRadius", 0.4);
0223   desc.add<double>("stripPtThreshold", 10.);
0224   desc.add<double>("widthPtThreshold", 3.);
0225   descriptions.add("hfJetShowerShape", desc);
0226 }
0227
0228 //define this as a plug-in
0229 DEFINE_FWK_MODULE(HFJetShowerShape);