Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:34:44

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     double rhoCurSum = 0.;
0147     double rhomCurSum = 0.;
0148     for (int i = 0; i < nacc; i++) {
0149       if (etaVec[i] >= etaMin && etaVec[i] < etaMax) {
0150         rhoVecCur.push_back(rhoVec[i]);
0151         rhomVecCur.push_back(rhomVec[i]);
0152 
0153         rhoCurSum += rhoVec[i];
0154         rhomCurSum += rhomVec[i];
0155         ++naccCur;
0156       }  //eta selection
0157     }    //accepted jet loop
0158 
0159     if (naccCur > 0) {
0160       double rhoCur = calcMedian(rhoVecCur);
0161       double rhomCur = calcMedian(rhomVecCur);
0162       mapToRhoOut->at(ieta) = rhoCur;
0163       mapToRhoMOut->at(ieta) = rhomCur;
0164     }
0165   }  //eta ranges
0166 
0167   iEvent.put(std::move(mapEtaRangesOut), "mapEtaEdges");
0168   iEvent.put(std::move(mapToRhoOut), "mapToRho");
0169   iEvent.put(std::move(mapToRhoMOut), "mapToRhoM");
0170 
0171   iEvent.put(std::move(ptJetsOut), "ptJets");
0172   iEvent.put(std::move(areaJetsOut), "areaJets");
0173   iEvent.put(std::move(etaJetsOut), "etaJets");
0174 
0175   return;
0176 }
0177 
0178 // ------------ method called once each stream before processing any runs, lumis or events  ------------
0179 void HiFJRhoProducer::beginStream(edm::StreamID) {}
0180 
0181 // ------------ method called once each stream after processing all runs, lumis and events  ------------
0182 void HiFJRhoProducer::endStream() {}
0183 
0184 double HiFJRhoProducer::calcMd(const reco::Jet* jet) {
0185   //
0186   //get md as defined in http://arxiv.org/pdf/1211.2811.pdf
0187   //
0188 
0189   //Loop over the jet constituents
0190   double sum = 0.;
0191   for (auto daughter : jet->getJetConstituentsQuick()) {
0192     if (isPackedCandidate(daughter)) {  //packed candidate situation
0193       auto part = static_cast<const pat::PackedCandidate*>(daughter);
0194       sum += sqrt(part->mass() * part->mass() + part->pt() * part->pt()) - part->pt();
0195     } else {
0196       auto part = static_cast<const reco::PFCandidate*>(daughter);
0197       sum += sqrt(part->mass() * part->mass() + part->pt() * part->pt()) - part->pt();
0198     }
0199   }
0200 
0201   return sum;
0202 }
0203 
0204 bool HiFJRhoProducer::isPackedCandidate(const reco::Candidate* candidate) {
0205   if (checkJetCand) {
0206     if (typeid(pat::PackedCandidate) == typeid(*candidate))
0207       usingPackedCand = true;
0208     else if (typeid(reco::PFCandidate) == typeid(*candidate))
0209       usingPackedCand = false;
0210     else
0211       throw cms::Exception("WrongJetCollection", "Jet constituents are not particle flow candidates");
0212     checkJetCand = false;
0213   }
0214   return usingPackedCand;
0215 }
0216 
0217 void HiFJRhoProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0218   edm::ParameterSetDescription desc;
0219   desc.add<edm::InputTag>("jetSource", edm::InputTag("kt4PFJets"));
0220   desc.add<int>("nExcl", 2);
0221   desc.add<double>("etaMaxExcl", 2.);
0222   desc.add<double>("ptMinExcl", 20.);
0223   desc.add<int>("nExcl2", 2);
0224   desc.add<double>("etaMaxExcl2", 2.);
0225   desc.add<double>("ptMinExcl2", 20.);
0226   desc.add<std::vector<double>>("etaRanges", {});
0227   descriptions.add("hiFJRhoProducer", desc);
0228 }
0229 
0230 //--------- method to calculate median ------------------
0231 double HiFJRhoProducer::calcMedian(std::vector<double>& v) {
0232   //post-condition: After returning, the elements in v may be reordered and the resulting order is implementation defined.
0233   //works for even and odd collections
0234   if (v.empty()) {
0235     return 0.0;
0236   }
0237   auto n = v.size() / 2;
0238   std::nth_element(v.begin(), v.begin() + n, v.end());
0239   auto med = v[n];
0240   if (!(v.size() & 1)) {  //If the set size is even
0241     auto max_it = std::max_element(v.begin(), v.begin() + n);
0242     med = (*max_it + med) / 2.0;
0243   }
0244   return med;
0245 }
0246 
0247 //define this as a plug-in
0248 DEFINE_FWK_MODULE(HiFJRhoProducer);