Skip to content

Commit

Permalink
Merge pull request #45570 from fwyzard/cleanup_HiFJRhoProducer
Browse files Browse the repository at this point in the history 
Clean up `HiFJRhoProducer`
  • Loading branch information
@cmsbuild
cmsbuild authored Jul 27, 2024
2 parents 5c2d038 + 3ba8706 commit b08f87e
Show file tree
Hide file tree
Showing 2 changed files with 109 additions and 185 deletions.
236 changes: 109 additions & 127 deletions RecoHI/HiJetAlgos/plugins/HiFJRhoProducer.cc
View file
Original file line number Diff line number Diff line change
@@ -1,212 +1,193 @@
// -*- C++ -*-
//
// Package: HiJetBackground/HiFJRhoProducer
// Class: HiFJRhoProducer
//
/**\class HiFJRhoProducer HiFJRhoProducer.cc HiJetBackground/HiFJRhoProducer/plugins/HiFJRhoProducer.cc
Description: [one line class summary]
Implementation:
[Notes on implementation]
*/
//
// Original Author: Marta Verweij
// Created: Thu, 16 Jul 2015 10:57:12 GMT
//
//

// system include files
#include <memory>
#include <string>
#include <vector>

#include "RecoHI/HiJetAlgos/plugins/HiFJRhoProducer.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/interface/ConsumesCollector.h"

#include "DataFormats/Common/interface/View.h"
// user include files
#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/PatCandidates/interface/PackedCandidate.h"
#include "DataFormats/Common/interface/View.h"
#include "DataFormats/JetReco/interface/Jet.h"
#include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
#include "DataFormats/PatCandidates/interface/Jet.h"
#include "DataFormats/PatCandidates/interface/PackedCandidate.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/global/EDProducer.h"
#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
#include "FWCore/Utilities/interface/StreamID.h"

// class declaration
class HiFJRhoProducer : public edm::global::EDProducer<> {
public:
explicit HiFJRhoProducer(const edm::ParameterSet&);
~HiFJRhoProducer() override = default;

static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);

private:
void produce(edm::StreamID, edm::Event&, edm::EventSetup const&) const override;

double calcMedian(std::vector<double>& v) const;
double calcMd(reco::Jet const& jet) const;
bool isPackedCandidate(reco::Candidate const* candidate) const;

// members
const edm::EDGetTokenT<edm::View<reco::Jet>> jetsToken_; // input kt jet source
const unsigned int nExcl_; // number of leading jets to exclude
const unsigned int nExcl2_; // number of leading jets to exclude in 2nd eta region
const double etaMaxExcl_; // max eta for jets to exclude
const double ptMinExcl_; // min pt for excluded jets
const double etaMaxExcl2_; // max eta for jets to exclude in 2nd eta region
const double ptMinExcl2_; // min pt for excluded jets in 2nd eta region
const std::vector<double> etaRanges_; // eta boundaries for rho calculation regions
mutable std::once_flag checkJetCand_; // check if using packed candidates and cache the information
mutable bool usingPackedCand_ = false;
};

//using namespace edm;
using namespace reco;
//using namespace pat;

//
// constants, enums and typedefs
//

//
// static data member definitions
//

//
// constructors and destructor
//
// constructor
HiFJRhoProducer::HiFJRhoProducer(const edm::ParameterSet& iConfig)
: src_(iConfig.getParameter<edm::InputTag>("jetSource")),
: jetsToken_(consumes<edm::View<reco::Jet>>(iConfig.getParameter<edm::InputTag>("jetSource"))),
nExcl_(iConfig.getParameter<int>("nExcl")),
nExcl2_(iConfig.getParameter<int>("nExcl2")),
etaMaxExcl_(iConfig.getParameter<double>("etaMaxExcl")),
ptMinExcl_(iConfig.getParameter<double>("ptMinExcl")),
nExcl2_(iConfig.getParameter<int>("nExcl2")),
etaMaxExcl2_(iConfig.getParameter<double>("etaMaxExcl2")),
ptMinExcl2_(iConfig.getParameter<double>("ptMinExcl2")),
checkJetCand(true),
usingPackedCand(false) {
jetsToken_ = consumes<edm::View<reco::Jet>>(src_);

//register your products
etaRanges_(iConfig.getParameter<std::vector<double>>("etaRanges")) {
// register your products
produces<std::vector<double>>("mapEtaEdges");
produces<std::vector<double>>("mapToRho");
produces<std::vector<double>>("mapToRhoM");
produces<std::vector<double>>("ptJets");
produces<std::vector<double>>("areaJets");
produces<std::vector<double>>("etaJets");
etaRanges = iConfig.getParameter<std::vector<double>>("etaRanges");
}

HiFJRhoProducer::~HiFJRhoProducer() {
// do anything here that needs to be done at desctruction time
// (e.g. close files, deallocate resources etc.)
}

// ------------ method called to produce the data ------------
void HiFJRhoProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
// Get the vector of jets
// method called for each event to produce the data
void HiFJRhoProducer::produce(edm::StreamID, edm::Event& iEvent, edm::EventSetup const& iSetup) const {
// get the inputs jets
edm::Handle<edm::View<reco::Jet>> jets;
iEvent.getByToken(jetsToken_, jets);

int neta = (int)etaRanges.size();
int neta = static_cast<int>(etaRanges_.size());
auto mapEtaRangesOut = std::make_unique<std::vector<double>>(neta, -999.);

for (int ieta = 0; ieta < neta; ieta++) {
mapEtaRangesOut->at(ieta) = etaRanges[ieta];
mapEtaRangesOut->at(ieta) = etaRanges_[ieta];
}
auto mapToRhoOut = std::make_unique<std::vector<double>>(neta - 1, 1e-6);
auto mapToRhoMOut = std::make_unique<std::vector<double>>(neta - 1, 1e-6);

int njets = jets->size();
int njets = static_cast<int>(jets->size());

auto ptJetsOut = std::make_unique<std::vector<double>>(njets, 1e-6);
auto areaJetsOut = std::make_unique<std::vector<double>>(njets, 1e-6);
auto etaJetsOut = std::make_unique<std::vector<double>>(njets, 1e-6);
auto ptJetsOut = std::make_unique<std::vector<double>>();
ptJetsOut->reserve(njets);
auto areaJetsOut = std::make_unique<std::vector<double>>();
areaJetsOut->reserve(njets);
auto etaJetsOut = std::make_unique<std::vector<double>>();
etaJetsOut->reserve(njets);

double rhoVec[999];
double rhomVec[999];
double etaVec[999];
std::vector<double> rhoVec;
rhoVec.reserve(njets);
std::vector<double> rhomVec;
rhomVec.reserve(njets);

// int neta = (int)mapEtaRangesOut->size();
int nacc = 0;
unsigned int njetsEx = 0;
unsigned int njetsEx2 = 0;
for (auto jet = jets->begin(); jet != jets->end(); ++jet) {
if (njetsEx < nExcl_ && fabs(jet->eta()) < etaMaxExcl_ && jet->pt() > ptMinExcl_) {
njetsEx++;
for (auto const& jet : *jets) {
if (njetsEx < nExcl_ and fabs(jet.eta()) < etaMaxExcl_ and jet.pt() > ptMinExcl_) {
++njetsEx;
continue;
}
if (njetsEx2 < nExcl2_ && fabs(jet->eta()) < etaMaxExcl2_ && fabs(jet->eta()) > etaMaxExcl_ &&
jet->pt() > ptMinExcl2_) {
njetsEx2++;
if (njetsEx2 < nExcl2_ and fabs(jet.eta()) < etaMaxExcl2_ and fabs(jet.eta()) > etaMaxExcl_ and
jet.pt() > ptMinExcl2_) {
++njetsEx2;
continue;
}
float pt = jet->pt();
float area = jet->jetArea();
float eta = jet->eta();
float pt = jet.pt();
float area = jet.jetArea();
float eta = jet.eta();

if (eta < mapEtaRangesOut->at(0) || eta > mapEtaRangesOut->at(neta - 1))
continue;
if (area > 0.) {
rhoVec[nacc] = pt / area;
rhomVec[nacc] = calcMd(&*jet) / area;
etaVec[nacc] = eta;
ptJetsOut->at(nacc) = pt;
areaJetsOut->at(nacc) = area;
etaJetsOut->at(nacc) = eta;
rhoVec.push_back(pt / area);
rhomVec.push_back(calcMd(jet) / area);
ptJetsOut->push_back(pt);
areaJetsOut->push_back(area);
etaJetsOut->push_back(eta);
++nacc;
}
}

ptJetsOut->resize(nacc);
areaJetsOut->resize(nacc);
etaJetsOut->resize(nacc);
//calculate rho and rhom in eta ranges
double radius = 0.2; //distance kt clusters needs to be from edge
for (int ieta = 0; ieta < (neta - 1); ieta++) {
// calculate rho and rhom in eta ranges
const double radius = 0.2; // distance kt clusters needs to be from edge
for (int ieta = 0; ieta < (neta - 1); ++ieta) {
std::vector<double> rhoVecCur;
rhoVecCur.reserve(nacc);
std::vector<double> rhomVecCur;
rhomVecCur.reserve(nacc);

double etaMin = mapEtaRangesOut->at(ieta) + radius;
double etaMax = mapEtaRangesOut->at(ieta + 1) - radius;

int naccCur = 0;
for (int i = 0; i < nacc; i++) {
if (etaVec[i] >= etaMin && etaVec[i] < etaMax) {
for (int i = 0; i < nacc; ++i) {
if ((*etaJetsOut)[i] >= etaMin and (*etaJetsOut)[i] < etaMax) {
rhoVecCur.push_back(rhoVec[i]);
rhomVecCur.push_back(rhomVec[i]);
++naccCur;
} //eta selection
} //accepted jet loop

if (naccCur > 0) {
double rhoCur = calcMedian(rhoVecCur);
double rhomCur = calcMedian(rhomVecCur);
mapToRhoOut->at(ieta) = rhoCur;
mapToRhoMOut->at(ieta) = rhomCur;
} // eta selection
} // accepted jet loop

if (not rhoVecCur.empty()) {
mapToRhoOut->at(ieta) = calcMedian(rhoVecCur);
mapToRhoMOut->at(ieta) = calcMedian(rhomVecCur);
}
} //eta ranges
} // eta ranges

iEvent.put(std::move(mapEtaRangesOut), "mapEtaEdges");
iEvent.put(std::move(mapToRhoOut), "mapToRho");
iEvent.put(std::move(mapToRhoMOut), "mapToRhoM");

iEvent.put(std::move(ptJetsOut), "ptJets");
iEvent.put(std::move(areaJetsOut), "areaJets");
iEvent.put(std::move(etaJetsOut), "etaJets");

return;
}

// ------------ method called once each stream before processing any runs, lumis or events ------------
void HiFJRhoProducer::beginStream(edm::StreamID) {}

// ------------ method called once each stream after processing all runs, lumis and events ------------
void HiFJRhoProducer::endStream() {}

double HiFJRhoProducer::calcMd(const reco::Jet* jet) {
//
//get md as defined in http://arxiv.org/pdf/1211.2811.pdf
//
double HiFJRhoProducer::calcMd(reco::Jet const& jet) const {
// compute md as defined in http://arxiv.org/pdf/1211.2811.pdf

//Loop over the jet constituents
// loop over the jet constituents
double sum = 0.;
for (auto daughter : jet->getJetConstituentsQuick()) {
if (isPackedCandidate(daughter)) { //packed candidate situation
auto part = static_cast<const pat::PackedCandidate*>(daughter);
for (auto const* daughter : jet.getJetConstituentsQuick()) {
if (isPackedCandidate(daughter)) { // packed candidate situation
auto part = static_cast<pat::PackedCandidate const*>(daughter);
sum += sqrt(part->mass() * part->mass() + part->pt() * part->pt()) - part->pt();
} else {
auto part = static_cast<const reco::PFCandidate*>(daughter);
auto part = static_cast<reco::PFCandidate const*>(daughter);
sum += sqrt(part->mass() * part->mass() + part->pt() * part->pt()) - part->pt();
}
}

return sum;
}

bool HiFJRhoProducer::isPackedCandidate(const reco::Candidate* candidate) {
if (checkJetCand) {
bool HiFJRhoProducer::isPackedCandidate(const reco::Candidate* candidate) const {
// check if using packed candidates on the first call and cache the information
std::call_once(checkJetCand_, [&]() {
if (typeid(pat::PackedCandidate) == typeid(*candidate))
usingPackedCand = true;
usingPackedCand_ = true;
else if (typeid(reco::PFCandidate) == typeid(*candidate))
usingPackedCand = false;
usingPackedCand_ = false;
else
throw cms::Exception("WrongJetCollection", "Jet constituents are not particle flow candidates");
checkJetCand = false;
}
return usingPackedCand;
});
return usingPackedCand_;
}

void HiFJRhoProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
Expand All @@ -223,21 +204,22 @@ void HiFJRhoProducer::fillDescriptions(edm::ConfigurationDescriptions& descripti
}

//--------- method to calculate median ------------------
double HiFJRhoProducer::calcMedian(std::vector<double>& v) {
//post-condition: After returning, the elements in v may be reordered and the resulting order is implementation defined.
//works for even and odd collections
double HiFJRhoProducer::calcMedian(std::vector<double>& v) const {
// post-condition: After returning, the elements in v may be reordered and the resulting order is implementation defined.
// works for even and odd collections
if (v.empty()) {
return 0.0;
}
auto n = v.size() / 2;
std::nth_element(v.begin(), v.begin() + n, v.end());
auto med = v[n];
if (!(v.size() & 1)) { //If the set size is even
if (!(v.size() & 1)) { // if the set size is even
auto max_it = std::max_element(v.begin(), v.begin() + n);
med = (*max_it + med) / 2.0;
}
return med;
}

//define this as a plug-in
// define this as a plug-in
#include "FWCore/Framework/interface/MakerMacros.h"
DEFINE_FWK_MODULE(HiFJRhoProducer);
Loading

0 comments on commit b08f87e

Please sign in to comment.