HiJetAlgos/plugins/HiFJRhoProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    HiJetBackground/HiFJRhoProducer
0004 // Class:      HiFJRhoProducer
0005 //
0006 /**\class HiFJRhoProducer HiFJRhoProducer.cc HiJetBackground/HiFJRhoProducer/plugins/HiFJRhoProducer.cc
0007
0008  Description: [one line class summary]
0009
0010  Implementation:
0011      [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  Marta Verweij
0015 //         Created:  Thu, 16 Jul 2015 10:57:12 GMT
0016 //
0017 //
0018
0019 #include <memory>
0020 #include <string>
0021
0022 #include "RecoHI/HiJetAlgos/plugins/HiFJRhoProducer.h"
0023
0024 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0025 #include "FWCore/Framework/interface/ESHandle.h"
0026 #include "FWCore/Framework/interface/MakerMacros.h"
0027 #include "FWCore/Framework/interface/ConsumesCollector.h"
0028
0029 #include "DataFormats/Common/interface/View.h"
0030 #include "DataFormats/Common/interface/Handle.h"
0031 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0032 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0033 #include "DataFormats/PatCandidates/interface/Jet.h"
0034
0035 //using namespace edm;
0036 using namespace reco;
0037 //using namespace pat;
0038
0039 //
0040 // constants, enums and typedefs
0041 //
0042
0043 //
0044 // static data member definitions
0045 //
0046
0047 //
0048 // constructors and destructor
0049 //
0050 HiFJRhoProducer::HiFJRhoProducer(const edm::ParameterSet& iConfig)
0051     : src_(iConfig.getParameter<edm::InputTag>("jetSource")),
0052       nExcl_(iConfig.getParameter<int>("nExcl")),
0053       etaMaxExcl_(iConfig.getParameter<double>("etaMaxExcl")),
0054       ptMinExcl_(iConfig.getParameter<double>("ptMinExcl")),
0055       nExcl2_(iConfig.getParameter<int>("nExcl2")),
0056       etaMaxExcl2_(iConfig.getParameter<double>("etaMaxExcl2")),
0057       ptMinExcl2_(iConfig.getParameter<double>("ptMinExcl2")),
0058       checkJetCand(true),
0059       usingPackedCand(false) {
0060   jetsToken_ = consumes<edm::View<reco::Jet>>(src_);
0061
0062   //register your products
0063   produces<std::vector<double>>("mapEtaEdges");
0064   produces<std::vector<double>>("mapToRho");
0065   produces<std::vector<double>>("mapToRhoM");
0066   produces<std::vector<double>>("ptJets");
0067   produces<std::vector<double>>("areaJets");
0068   produces<std::vector<double>>("etaJets");
0069   etaRanges = iConfig.getParameter<std::vector<double>>("etaRanges");
0070 }
0071
0072 HiFJRhoProducer::~HiFJRhoProducer() {
0073   // do anything here that needs to be done at desctruction time
0074   // (e.g. close files, deallocate resources etc.)
0075 }
0076
0077 // ------------ method called to produce the data  ------------
0078 void HiFJRhoProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0079   // Get the vector of jets
0080   edm::Handle<edm::View<reco::Jet>> jets;
0081   iEvent.getByToken(jetsToken_, jets);
0082
0083   int neta = (int)etaRanges.size();
0084   auto mapEtaRangesOut = std::make_unique<std::vector<double>>(neta, -999.);
0085
0086   for (int ieta = 0; ieta < neta; ieta++) {
0087     mapEtaRangesOut->at(ieta) = etaRanges[ieta];
0088   }
0089   auto mapToRhoOut = std::make_unique<std::vector<double>>(neta - 1, 1e-6);
0090   auto mapToRhoMOut = std::make_unique<std::vector<double>>(neta - 1, 1e-6);
0091
0092   int njets = jets->size();
0093
0094   auto ptJetsOut = std::make_unique<std::vector<double>>(njets, 1e-6);
0095   auto areaJetsOut = std::make_unique<std::vector<double>>(njets, 1e-6);
0096   auto etaJetsOut = std::make_unique<std::vector<double>>(njets, 1e-6);
0097
0098   double rhoVec[999];
0099   double rhomVec[999];
0100   double etaVec[999];
0101
0102   // int neta = (int)mapEtaRangesOut->size();
0103   int nacc = 0;
0104   unsigned int njetsEx = 0;
0105   unsigned int njetsEx2 = 0;
0106   for (auto jet = jets->begin(); jet != jets->end(); ++jet) {
0107     if (njetsEx < nExcl_ && fabs(jet->eta()) < etaMaxExcl_ && jet->pt() > ptMinExcl_) {
0108       njetsEx++;
0109       continue;
0110     }
0111     if (njetsEx2 < nExcl2_ && fabs(jet->eta()) < etaMaxExcl2_ && fabs(jet->eta()) > etaMaxExcl_ &&
0112         jet->pt() > ptMinExcl2_) {
0113       njetsEx2++;
0114       continue;
0115     }
0116     float pt = jet->pt();
0117     float area = jet->jetArea();
0118     float eta = jet->eta();
0119
0120     if (eta < mapEtaRangesOut->at(0) || eta > mapEtaRangesOut->at(neta - 1))
0121       continue;
0122     if (area > 0.) {
0123       rhoVec[nacc] = pt / area;
0124       rhomVec[nacc] = calcMd(&*jet) / area;
0125       etaVec[nacc] = eta;
0126       ptJetsOut->at(nacc) = pt;
0127       areaJetsOut->at(nacc) = area;
0128       etaJetsOut->at(nacc) = eta;
0129       ++nacc;
0130     }
0131   }
0132
0133   ptJetsOut->resize(nacc);
0134   areaJetsOut->resize(nacc);
0135   etaJetsOut->resize(nacc);
0136   //calculate rho and rhom in eta ranges
0137   double radius = 0.2;  //distance kt clusters needs to be from edge
0138   for (int ieta = 0; ieta < (neta - 1); ieta++) {
0139     std::vector<double> rhoVecCur;
0140     std::vector<double> rhomVecCur;
0141
0142     double etaMin = mapEtaRangesOut->at(ieta) + radius;
0143     double etaMax = mapEtaRangesOut->at(ieta + 1) - radius;
0144
0145     int naccCur = 0;
0146     for (int i = 0; i < nacc; i++) {
0147       if (etaVec[i] >= etaMin && etaVec[i] < etaMax) {
0148         rhoVecCur.push_back(rhoVec[i]);
0149         rhomVecCur.push_back(rhomVec[i]);
0150         ++naccCur;
0151       }  //eta selection
0152     }    //accepted jet loop
0153
0154     if (naccCur > 0) {
0155       double rhoCur = calcMedian(rhoVecCur);
0156       double rhomCur = calcMedian(rhomVecCur);
0157       mapToRhoOut->at(ieta) = rhoCur;
0158       mapToRhoMOut->at(ieta) = rhomCur;
0159     }
0160   }  //eta ranges
0161
0162   iEvent.put(std::move(mapEtaRangesOut), "mapEtaEdges");
0163   iEvent.put(std::move(mapToRhoOut), "mapToRho");
0164   iEvent.put(std::move(mapToRhoMOut), "mapToRhoM");
0165
0166   iEvent.put(std::move(ptJetsOut), "ptJets");
0167   iEvent.put(std::move(areaJetsOut), "areaJets");
0168   iEvent.put(std::move(etaJetsOut), "etaJets");
0169
0170   return;
0171 }
0172
0173 // ------------ method called once each stream before processing any runs, lumis or events  ------------
0174 void HiFJRhoProducer::beginStream(edm::StreamID) {}
0175
0176 // ------------ method called once each stream after processing all runs, lumis and events  ------------
0177 void HiFJRhoProducer::endStream() {}
0178
0179 double HiFJRhoProducer::calcMd(const reco::Jet* jet) {
0180   //
0181   //get md as defined in http://arxiv.org/pdf/1211.2811.pdf
0182   //
0183
0184   //Loop over the jet constituents
0185   double sum = 0.;
0186   for (auto daughter : jet->getJetConstituentsQuick()) {
0187     if (isPackedCandidate(daughter)) {  //packed candidate situation
0188       auto part = static_cast<const pat::PackedCandidate*>(daughter);
0189       sum += sqrt(part->mass() * part->mass() + part->pt() * part->pt()) - part->pt();
0190     } else {
0191       auto part = static_cast<const reco::PFCandidate*>(daughter);
0192       sum += sqrt(part->mass() * part->mass() + part->pt() * part->pt()) - part->pt();
0193     }
0194   }
0195
0196   return sum;
0197 }
0198
0199 bool HiFJRhoProducer::isPackedCandidate(const reco::Candidate* candidate) {
0200   if (checkJetCand) {
0201     if (typeid(pat::PackedCandidate) == typeid(*candidate))
0202       usingPackedCand = true;
0203     else if (typeid(reco::PFCandidate) == typeid(*candidate))
0204       usingPackedCand = false;
0205     else
0206       throw cms::Exception("WrongJetCollection", "Jet constituents are not particle flow candidates");
0207     checkJetCand = false;
0208   }
0209   return usingPackedCand;
0210 }
0211
0212 void HiFJRhoProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0213   edm::ParameterSetDescription desc;
0214   desc.add<edm::InputTag>("jetSource", edm::InputTag("kt4PFJets"));
0215   desc.add<int>("nExcl", 2);
0216   desc.add<double>("etaMaxExcl", 2.);
0217   desc.add<double>("ptMinExcl", 20.);
0218   desc.add<int>("nExcl2", 2);
0219   desc.add<double>("etaMaxExcl2", 2.);
0220   desc.add<double>("ptMinExcl2", 20.);
0221   desc.add<std::vector<double>>("etaRanges", {});
0222   descriptions.add("hiFJRhoProducer", desc);
0223 }
0224
0225 //--------- method to calculate median ------------------
0226 double HiFJRhoProducer::calcMedian(std::vector<double>& v) {
0227   //post-condition: After returning, the elements in v may be reordered and the resulting order is implementation defined.
0228   //works for even and odd collections
0229   if (v.empty()) {
0230     return 0.0;
0231   }
0232   auto n = v.size() / 2;
0233   std::nth_element(v.begin(), v.begin() + n, v.end());
0234   auto med = v[n];
0235   if (!(v.size() & 1)) {  //If the set size is even
0236     auto max_it = std::max_element(v.begin(), v.begin() + n);
0237     med = (*max_it + med) / 2.0;
0238   }
0239   return med;
0240 }
0241
0242 //define this as a plug-in
0243 DEFINE_FWK_MODULE(HiFJRhoProducer);