HiJetAlgos/plugins/HiFJRhoProducer.cc

0001 // Original Author:  Marta Verweij
0002 //         Created:  Thu, 16 Jul 2015 10:57:12 GMT
0003
0004 // system include files
0005 #include <memory>
0006 #include <string>
0007 #include <vector>
0008
0009 // user include files
0010 #include "DataFormats/Common/interface/Handle.h"
0011 #include "DataFormats/Common/interface/View.h"
0012 #include "DataFormats/JetReco/interface/Jet.h"
0013 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0014 #include "DataFormats/PatCandidates/interface/Jet.h"
0015 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0016 #include "FWCore/Framework/interface/Event.h"
0017 #include "FWCore/Framework/interface/global/EDProducer.h"
0018 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0019 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0020 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0021 #include "FWCore/Utilities/interface/StreamID.h"
0022
0023 // class declaration
0024 class HiFJRhoProducer : public edm::global::EDProducer<> {
0025 public:
0026   explicit HiFJRhoProducer(const edm::ParameterSet&);
0027   ~HiFJRhoProducer() override = default;
0028
0029   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0030
0031 private:
0032   void produce(edm::StreamID, edm::Event&, edm::EventSetup const&) const override;
0033
0034   double calcMedian(std::vector<double>& v) const;
0035   double calcMd(reco::Jet const& jet) const;
0036   bool isPackedCandidate(reco::Candidate const* candidate) const;
0037
0038   // members
0039   const edm::EDGetTokenT<edm::View<reco::Jet>> jetsToken_;  // input kt jet source
0040   const unsigned int nExcl_;                                // number of leading jets to exclude
0041   const unsigned int nExcl2_;                               // number of leading jets to exclude in 2nd eta region
0042   const double etaMaxExcl_;                                 // max eta for jets to exclude
0043   const double ptMinExcl_;                                  // min pt for excluded jets
0044   const double etaMaxExcl2_;                                // max eta for jets to exclude in 2nd eta region
0045   const double ptMinExcl2_;                                 // min pt for excluded jets in 2nd eta region
0046   const std::vector<double> etaRanges_;                     // eta boundaries for rho calculation regions
0047   mutable std::once_flag checkJetCand_;  // check if using packed candidates and cache the information
0048   mutable bool usingPackedCand_ = false;
0049 };
0050
0051 using namespace reco;
0052
0053 // constructor
0054 HiFJRhoProducer::HiFJRhoProducer(const edm::ParameterSet& iConfig)
0055     : jetsToken_(consumes<edm::View<reco::Jet>>(iConfig.getParameter<edm::InputTag>("jetSource"))),
0056       nExcl_(iConfig.getParameter<int>("nExcl")),
0057       nExcl2_(iConfig.getParameter<int>("nExcl2")),
0058       etaMaxExcl_(iConfig.getParameter<double>("etaMaxExcl")),
0059       ptMinExcl_(iConfig.getParameter<double>("ptMinExcl")),
0060       etaMaxExcl2_(iConfig.getParameter<double>("etaMaxExcl2")),
0061       ptMinExcl2_(iConfig.getParameter<double>("ptMinExcl2")),
0062       etaRanges_(iConfig.getParameter<std::vector<double>>("etaRanges")) {
0063   // register your products
0064   produces<std::vector<double>>("mapEtaEdges");
0065   produces<std::vector<double>>("mapToRho");
0066   produces<std::vector<double>>("mapToRhoM");
0067   produces<std::vector<double>>("ptJets");
0068   produces<std::vector<double>>("areaJets");
0069   produces<std::vector<double>>("etaJets");
0070 }
0071
0072 // method called for each event to produce the data
0073 void HiFJRhoProducer::produce(edm::StreamID, edm::Event& iEvent, edm::EventSetup const& iSetup) const {
0074   // get the inputs jets
0075   edm::Handle<edm::View<reco::Jet>> jets;
0076   iEvent.getByToken(jetsToken_, jets);
0077
0078   int neta = static_cast<int>(etaRanges_.size());
0079   auto mapEtaRangesOut = std::make_unique<std::vector<double>>(neta, -999.);
0080
0081   for (int ieta = 0; ieta < neta; ieta++) {
0082     mapEtaRangesOut->at(ieta) = etaRanges_[ieta];
0083   }
0084   auto mapToRhoOut = std::make_unique<std::vector<double>>(neta - 1, 1e-6);
0085   auto mapToRhoMOut = std::make_unique<std::vector<double>>(neta - 1, 1e-6);
0086
0087   int njets = static_cast<int>(jets->size());
0088
0089   auto ptJetsOut = std::make_unique<std::vector<double>>();
0090   ptJetsOut->reserve(njets);
0091   auto areaJetsOut = std::make_unique<std::vector<double>>();
0092   areaJetsOut->reserve(njets);
0093   auto etaJetsOut = std::make_unique<std::vector<double>>();
0094   etaJetsOut->reserve(njets);
0095
0096   std::vector<double> rhoVec;
0097   rhoVec.reserve(njets);
0098   std::vector<double> rhomVec;
0099   rhomVec.reserve(njets);
0100
0101   int nacc = 0;
0102   unsigned int njetsEx = 0;
0103   unsigned int njetsEx2 = 0;
0104   for (auto const& jet : *jets) {
0105     if (njetsEx < nExcl_ and fabs(jet.eta()) < etaMaxExcl_ and jet.pt() > ptMinExcl_) {
0106       ++njetsEx;
0107       continue;
0108     }
0109     if (njetsEx2 < nExcl2_ and fabs(jet.eta()) < etaMaxExcl2_ and fabs(jet.eta()) > etaMaxExcl_ and
0110         jet.pt() > ptMinExcl2_) {
0111       ++njetsEx2;
0112       continue;
0113     }
0114     float pt = jet.pt();
0115     float area = jet.jetArea();
0116     float eta = jet.eta();
0117
0118     if (eta < mapEtaRangesOut->at(0) || eta > mapEtaRangesOut->at(neta - 1))
0119       continue;
0120     if (area > 0.) {
0121       rhoVec.push_back(pt / area);
0122       rhomVec.push_back(calcMd(jet) / area);
0123       ptJetsOut->push_back(pt);
0124       areaJetsOut->push_back(area);
0125       etaJetsOut->push_back(eta);
0126       ++nacc;
0127     }
0128   }
0129
0130   // calculate rho and rhom in eta ranges
0131   const double radius = 0.2;  // distance kt clusters needs to be from edge
0132   for (int ieta = 0; ieta < (neta - 1); ++ieta) {
0133     std::vector<double> rhoVecCur;
0134     rhoVecCur.reserve(nacc);
0135     std::vector<double> rhomVecCur;
0136     rhomVecCur.reserve(nacc);
0137
0138     double etaMin = mapEtaRangesOut->at(ieta) + radius;
0139     double etaMax = mapEtaRangesOut->at(ieta + 1) - radius;
0140
0141     for (int i = 0; i < nacc; ++i) {
0142       if ((*etaJetsOut)[i] >= etaMin and (*etaJetsOut)[i] < etaMax) {
0143         rhoVecCur.push_back(rhoVec[i]);
0144         rhomVecCur.push_back(rhomVec[i]);
0145       }  // eta selection
0146     }  // accepted jet loop
0147
0148     if (not rhoVecCur.empty()) {
0149       mapToRhoOut->at(ieta) = calcMedian(rhoVecCur);
0150       mapToRhoMOut->at(ieta) = calcMedian(rhomVecCur);
0151     }
0152   }  // eta ranges
0153
0154   iEvent.put(std::move(mapEtaRangesOut), "mapEtaEdges");
0155   iEvent.put(std::move(mapToRhoOut), "mapToRho");
0156   iEvent.put(std::move(mapToRhoMOut), "mapToRhoM");
0157   iEvent.put(std::move(ptJetsOut), "ptJets");
0158   iEvent.put(std::move(areaJetsOut), "areaJets");
0159   iEvent.put(std::move(etaJetsOut), "etaJets");
0160 }
0161
0162 double HiFJRhoProducer::calcMd(reco::Jet const& jet) const {
0163   // compute md as defined in http://arxiv.org/pdf/1211.2811.pdf
0164
0165   // loop over the jet constituents
0166   double sum = 0.;
0167   for (auto const* daughter : jet.getJetConstituentsQuick()) {
0168     if (isPackedCandidate(daughter)) {  // packed candidate situation
0169       auto part = static_cast<pat::PackedCandidate const*>(daughter);
0170       sum += sqrt(part->mass() * part->mass() + part->pt() * part->pt()) - part->pt();
0171     } else {
0172       auto part = static_cast<reco::PFCandidate const*>(daughter);
0173       sum += sqrt(part->mass() * part->mass() + part->pt() * part->pt()) - part->pt();
0174     }
0175   }
0176
0177   return sum;
0178 }
0179
0180 bool HiFJRhoProducer::isPackedCandidate(const reco::Candidate* candidate) const {
0181   // check if using packed candidates on the first call and cache the information
0182   std::call_once(checkJetCand_, [&]() {
0183     if (typeid(pat::PackedCandidate) == typeid(*candidate))
0184       usingPackedCand_ = true;
0185     else if (typeid(reco::PFCandidate) == typeid(*candidate))
0186       usingPackedCand_ = false;
0187     else
0188       throw cms::Exception("WrongJetCollection", "Jet constituents are not particle flow candidates");
0189   });
0190   return usingPackedCand_;
0191 }
0192
0193 void HiFJRhoProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0194   edm::ParameterSetDescription desc;
0195   desc.add<edm::InputTag>("jetSource", edm::InputTag("kt4PFJets"));
0196   desc.add<int>("nExcl", 2);
0197   desc.add<double>("etaMaxExcl", 2.);
0198   desc.add<double>("ptMinExcl", 20.);
0199   desc.add<int>("nExcl2", 2);
0200   desc.add<double>("etaMaxExcl2", 2.);
0201   desc.add<double>("ptMinExcl2", 20.);
0202   desc.add<std::vector<double>>("etaRanges", {});
0203   descriptions.add("hiFJRhoProducer", desc);
0204 }
0205
0206 //--------- method to calculate median ------------------
0207 double HiFJRhoProducer::calcMedian(std::vector<double>& v) const {
0208   // post-condition: After returning, the elements in v may be reordered and the resulting order is implementation defined.
0209   // works for even and odd collections
0210   if (v.empty()) {
0211     return 0.0;
0212   }
0213   auto n = v.size() / 2;
0214   std::nth_element(v.begin(), v.begin() + n, v.end());
0215   auto med = v[n];
0216   if (!(v.size() & 1)) {  // if the set size is even
0217     auto max_it = std::max_element(v.begin(), v.begin() + n);
0218     med = (*max_it + med) / 2.0;
0219   }
0220   return med;
0221 }
0222
0223 // define this as a plug-in
0224 #include "FWCore/Framework/interface/MakerMacros.h"
0225 DEFINE_FWK_MODULE(HiFJRhoProducer);