CMS 3D CMS Logo

List of all members | Public Types | Public Member Functions | Static Public Member Functions | Private Types | Private Member Functions | Private Attributes | Static Private Attributes
HiPuRhoProducer Class Reference

#include <RecoHI/HiJetAlgos/plugins/HiPuRhoProducer.cc>

Inheritance diagram for HiPuRhoProducer:
edm::stream::EDProducer<>

Public Types

using ClusterSequencePtr = std::shared_ptr< fastjet::ClusterSequence >
 
using JetDefPtr = std::shared_ptr< fastjet::JetDefinition >
 
- Public Types inherited from edm::stream::EDProducer<>
using CacheTypes = CacheContexts< T... >
 
using GlobalCache = typename CacheTypes::GlobalCache
 
using HasAbility = AbilityChecker< T... >
 
using InputProcessBlockCache = typename CacheTypes::InputProcessBlockCache
 
using LuminosityBlockCache = typename CacheTypes::LuminosityBlockCache
 
using LuminosityBlockContext = LuminosityBlockContextT< LuminosityBlockCache, RunCache, GlobalCache >
 
using LuminosityBlockSummaryCache = typename CacheTypes::LuminosityBlockSummaryCache
 
using RunCache = typename CacheTypes::RunCache
 
using RunContext = RunContextT< RunCache, GlobalCache >
 
using RunSummaryCache = typename CacheTypes::RunSummaryCache
 

Public Member Functions

virtual void calculateOrphanInput (std::vector< fastjet::PseudoJet > &orphanInput)
 
virtual void calculatePedestal (std::vector< fastjet::PseudoJet > const &coll)
 
 HiPuRhoProducer (const edm::ParameterSet &)
 
virtual void putRho (edm::Event &iEvent, const edm::EventSetup &iSetup)
 
virtual void setupGeometryMap (edm::Event &iEvent, const edm::EventSetup &iSetup)
 
virtual void subtractPedestal (std::vector< fastjet::PseudoJet > &coll)
 
- Public Member Functions inherited from edm::stream::EDProducer<>
 EDProducer ()=default
 
 EDProducer (const EDProducer &)=delete
 
bool hasAbilityToProduceInBeginLumis () const final
 
bool hasAbilityToProduceInBeginProcessBlocks () const final
 
bool hasAbilityToProduceInBeginRuns () const final
 
bool hasAbilityToProduceInEndLumis () const final
 
bool hasAbilityToProduceInEndProcessBlocks () const final
 
bool hasAbilityToProduceInEndRuns () const final
 
const EDProduceroperator= (const EDProducer &)=delete
 

Static Public Member Functions

static void fillDescriptions (edm::ConfigurationDescriptions &descriptions)
 

Private Types

typedef std::pair< double, double > EtaPhi
 

Private Member Functions

virtual void inputTowers ()
 
void produce (edm::Event &, const edm::EventSetup &) override
 

Private Attributes

const edm::ESGetToken< CaloGeometry, CaloGeometryRecordcaloGeometryToken_
 
const edm::EDGetTokenT< CaloTowerCollectioncaloTowerToken_
 
const bool dropZeroTowers_
 
std::map< int, double > emean_
 
std::map< int, double > esigma_
 
std::vector< double > etaEdgeHi_
 
std::vector< double > etaEdgeLow_
 
std::vector< double > etaEdges_
 
std::map< int, std::array< double, 4 > > eTop4_
 
ClusterSequencePtr fjClusterSeq_
 
std::vector< fastjet::PseudoJet > fjInputs_
 
JetDefPtr fjJetDefinition_
 
std::vector< fastjet::PseudoJet > fjJets_
 
std::vector< fastjet::PseudoJet > fjOriginalInputs_
 
CaloGeometry const * geo_ = nullptr
 
std::map< int, int > geomtowers_
 
int ietamax_
 
int ietamin_
 
const int initialValue_ = -99
 
std::vector< const CaloTower * > inputs_
 
const int medianWindowWidth_
 
const double minimumTowersFraction_
 
const double nSigmaPU_
 
std::vector< int > nTow_
 
std::map< int, int > ntowersWithJets_
 
bool postOrphan_ = false
 
const double puPtMin_
 
const double radiusPU_
 
std::vector< double > rho_
 
std::vector< double > rhoExtra_
 
std::vector< double > rhoM_
 
const double rParam_
 
bool setInitialValue_
 
std::vector< EtaPhiTowertowermap_
 
std::vector< double > towExcludeEta_
 
std::vector< double > towExcludePhi_
 
std::vector< double > towExcludePt_
 
const double towSigmaCut_
 
std::array< float, nEtaTow_vmean0_
 
std::array< float, nEtaTow_vmean1_
 
std::array< int, nEtaTow_vngeom_
 
std::array< int, nEtaTow_vntow_
 
std::array< float, nEtaTow_vrho0_
 
std::array< float, nEtaTow_vrho1_
 
std::array< float, nEtaTow_vrms0_
 
std::array< float, nEtaTow_vrms1_
 

Static Private Attributes

static constexpr int nEtaTow_ = 82
 

Detailed Description

Description: Producer to dump Pu-jet style rho into event content Implementation: Just see MultipleAlgoIterator - re-implenting for use in CS jet with sigma subtraction and zeroing

Definition at line 68 of file HiPuRhoProducer.cc.

Member Typedef Documentation

◆ ClusterSequencePtr

using HiPuRhoProducer::ClusterSequencePtr = std::shared_ptr<fastjet::ClusterSequence>

Definition at line 72 of file HiPuRhoProducer.cc.

◆ EtaPhi

typedef std::pair<double, double> HiPuRhoProducer::EtaPhi
private

Definition at line 144 of file HiPuRhoProducer.cc.

◆ JetDefPtr

using HiPuRhoProducer::JetDefPtr = std::shared_ptr<fastjet::JetDefinition>

Definition at line 73 of file HiPuRhoProducer.cc.

Constructor & Destructor Documentation

◆ HiPuRhoProducer()

HiPuRhoProducer::HiPuRhoProducer ( const edm::ParameterSet iConfig)
explicit

Definition at line 148 of file HiPuRhoProducer.cc.

149  : dropZeroTowers_(iConfig.getParameter<bool>("dropZeroTowers")),
150  medianWindowWidth_(iConfig.getParameter<int>("medianWindowWidth")),
151  minimumTowersFraction_(iConfig.getParameter<double>("minimumTowersFraction")),
152  nSigmaPU_(iConfig.getParameter<double>("nSigmaPU")),
153  puPtMin_(iConfig.getParameter<double>("puPtMin")),
154  radiusPU_(iConfig.getParameter<double>("radiusPU")),
155  rParam_(iConfig.getParameter<double>("rParam")),
156  towSigmaCut_(iConfig.getParameter<double>("towSigmaCut")),
157  caloTowerToken_(consumes<CaloTowerCollection>(iConfig.getParameter<edm::InputTag>("src"))),
158  caloGeometryToken_(esConsumes<CaloGeometry, CaloGeometryRecord>(edm::ESInputTag{})) {
159  //register your products
160  produces<std::vector<double>>("mapEtaEdges");
161  produces<std::vector<double>>("mapToRho");
162  produces<std::vector<double>>("mapToRhoMedian");
163  produces<std::vector<double>>("mapToRhoExtra");
164  produces<std::vector<double>>("mapToRhoM");
165  produces<std::vector<int>>("mapToNTow");
166  produces<std::vector<double>>("mapToTowExcludePt");
167  produces<std::vector<double>>("mapToTowExcludePhi");
168  produces<std::vector<double>>("mapToTowExcludeEta");
169 }
T getParameter(std::string const &) const
Definition: ParameterSet.h:307
const double radiusPU_
const edm::ESGetToken< CaloGeometry, CaloGeometryRecord > caloGeometryToken_
const double nSigmaPU_
const edm::EDGetTokenT< CaloTowerCollection > caloTowerToken_
const double rParam_
const double towSigmaCut_
const double puPtMin_
const bool dropZeroTowers_
const int medianWindowWidth_
const double minimumTowersFraction_

Member Function Documentation

◆ calculateOrphanInput()

void HiPuRhoProducer::calculateOrphanInput ( std::vector< fastjet::PseudoJet > &  orphanInput)
virtual

Definition at line 449 of file HiPuRhoProducer.cc.

References reco::deltaR2(), reco::LeafCandidate::energy(), reco::LeafCandidate::eta(), spr::find(), fjInputs_, fjJets_, fjOriginalInputs_, geomtowers_, CaloTower::ieta(), input, inputs_, CaloTower::iphi(), LogDebug, minimumTowersFraction_, ntowersWithJets_, reco::LeafCandidate::phi(), postOrphan_, reco::LeafCandidate::pt(), puPtMin_, reco::LeafCandidate::px(), reco::LeafCandidate::py(), reco::LeafCandidate::pz(), radiusPU_, towermap_, towExcludeEta_, towExcludePhi_, and towExcludePt_.

Referenced by produce().

449  {
450  LogDebug("PileUpSubtractor") << "The subtractor calculating orphan input...\n";
451 
453 
454  std::vector<int> jettowers; // vector of towers indexed by "user_index"
455  std::map<std::pair<int, int>, int> excludedTowers; // map of excluded ieta, iphi values
456  for (auto const& pseudojetTMP : fjJets_) {
457  EtaPhi jet_etaphi(pseudojetTMP.eta(), pseudojetTMP.phi());
458 
459  if (pseudojetTMP.perp() < puPtMin_)
460  continue;
461 
462  for (auto const& im : towermap_) {
463  double dr2 = reco::deltaR2(im.eta, im.phi, jet_etaphi.first, jet_etaphi.second);
464  if (dr2 < radiusPU_ * radiusPU_ && !excludedTowers[std::pair(im.ieta, im.iphi)] &&
465  (geomtowers_[im.ieta] - ntowersWithJets_[im.ieta]) > minimumTowersFraction_ * geomtowers_[im.ieta]) {
466  ntowersWithJets_[im.ieta]++;
467  excludedTowers[std::pair(im.ieta, im.iphi)] = 1;
468  }
469  }
470 
471  for (auto const& input : fjInputs_) {
472  int index = input.user_index();
473  const CaloTower* originalTower = inputs_[index];
474  int ie = originalTower->ieta();
475  int ip = originalTower->iphi();
476 
477  if (excludedTowers[std::pair<int, int>(ie, ip)]) {
478  jettowers.push_back(index);
479  }
480  }
481  }
482  // Create a new collections from the towers not included in jets
483  for (auto const& input : fjInputs_) {
484  int index = input.user_index();
485  const CaloTower* originalTower = inputs_[index];
486  auto itjet = std::find(jettowers.begin(), jettowers.end(), index);
487  if (itjet == jettowers.end()) {
488  orphanInput.emplace_back(originalTower->px(), originalTower->py(), originalTower->pz(), originalTower->energy());
489  orphanInput.back().set_user_index(index);
490  } else {
491  towExcludePt_.push_back(originalTower->pt());
492  towExcludePhi_.push_back(originalTower->phi());
493  towExcludeEta_.push_back(originalTower->eta());
494  }
495  }
496 
497  postOrphan_ = true;
498 }
std::vector< double > towExcludePt_
const double radiusPU_
double pz() const final
z coordinate of momentum vector
double pt() const final
transverse momentum
std::vector< EtaPhiTower > towermap_
std::vector< double > towExcludeEta_
std::vector< fastjet::PseudoJet > fjInputs_
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:19
std::vector< double > towExcludePhi_
static std::string const input
Definition: EdmProvDump.cc:50
double px() const final
x coordinate of momentum vector
std::vector< const CaloTower * > inputs_
double py() const final
y coordinate of momentum vector
std::pair< double, double > EtaPhi
const double puPtMin_
constexpr auto deltaR2(const T1 &t1, const T2 &t2) -> decltype(t1.eta())
Definition: deltaR.h:16
int iphi() const
Definition: CaloTower.h:202
int ieta() const
Definition: CaloTower.h:200
std::vector< fastjet::PseudoJet > fjJets_
std::map< int, int > geomtowers_
std::vector< fastjet::PseudoJet > fjOriginalInputs_
double phi() const final
momentum azimuthal angle
std::map< int, int > ntowersWithJets_
const double minimumTowersFraction_
#define LogDebug(id)
double energy() const final
energy
double eta() const final
momentum pseudorapidity

◆ calculatePedestal()

void HiPuRhoProducer::calculatePedestal ( std::vector< fastjet::PseudoJet > const &  coll)
virtual

Definition at line 273 of file HiPuRhoProducer.cc.

References funct::abs(), StorageManager_cfg::e1, emean_, esigma_, CaloTower::et(), hi::etaedge, etaEdgeHi_, etaEdgeLow_, electrons_cff::etaWidth, eTop4_, geomtowers_, submitPVValidationJobs::gt, mps_fire::i, CaloTower::ieta(), ietamax_, ietamin_, initialValue_, inputs_, ALPAKA_ACCELERATOR_NAMESPACE::vertexFinder::it, LogDebug, M_PI, SiStripPI::min, minimumTowersFraction_, nEtaTow_, nSigmaPU_, nt, nTow_, ntowersWithJets_, postOrphan_, rho_, rhoExtra_, rhoM_, setInitialValue_, Validation_hcalonly_cfi::sign, findQualityFiles::size, mathSSE::sqrt(), towSigmaCut_, vmean0_, vmean1_, vngeom_, vntow_, vrho0_, vrho1_, vrms0_, and vrms1_.

Referenced by produce().

273  {
274  LogDebug("PileUpSubtractor") << "The subtractor calculating pedestals...\n";
275 
276  std::map<int, double> emean2;
277  std::map<int, int> ntowers;
278 
279  int ietaold = -10000;
280  // Initial values for emean_, emean2, esigma_, ntowers
281 
282  for (int vi = 0; vi < nEtaTow_; ++vi) {
283  vntow_[vi] = initialValue_;
284  vmean1_[vi] = initialValue_;
285  vrms1_[vi] = initialValue_;
286  vrho1_[vi] = initialValue_;
287 
288  if (setInitialValue_) {
289  vmean0_[vi] = initialValue_;
290  vrms0_[vi] = initialValue_;
291  vrho0_[vi] = initialValue_;
292  }
293  }
294 
295  for (int i = ietamin_; i < ietamax_ + 1; i++) {
296  emean_[i] = 0.;
297  emean2[i] = 0.;
298  esigma_[i] = 0.;
299  ntowers[i] = 0;
300 
301  eTop4_[i] = {{0., 0., 0., 0.}};
302  }
303 
304  for (auto const& input_object : coll) {
305  const CaloTower* originalTower = inputs_[input_object.user_index()];
306  double original_Et = originalTower->et();
307  int ieta0 = originalTower->ieta();
308 
309  if (original_Et > eTop4_[ieta0][0]) {
310  eTop4_[ieta0][3] = eTop4_[ieta0][2];
311  eTop4_[ieta0][2] = eTop4_[ieta0][1];
312  eTop4_[ieta0][1] = eTop4_[ieta0][0];
313  eTop4_[ieta0][0] = original_Et;
314  } else if (original_Et > eTop4_[ieta0][1]) {
315  eTop4_[ieta0][3] = eTop4_[ieta0][2];
316  eTop4_[ieta0][2] = eTop4_[ieta0][1];
317  eTop4_[ieta0][1] = original_Et;
318  } else if (original_Et > eTop4_[ieta0][2]) {
319  eTop4_[ieta0][3] = eTop4_[ieta0][2];
320  eTop4_[ieta0][2] = original_Et;
321  } else if (original_Et > eTop4_[ieta0][3]) {
322  eTop4_[ieta0][3] = original_Et;
323  }
324 
325  emean_[ieta0] = emean_[ieta0] + original_Et;
326  emean2[ieta0] = emean2[ieta0] + original_Et * original_Et;
327  if (ieta0 - ietaold != 0) {
328  ntowers[ieta0] = 1;
329  ietaold = ieta0;
330  } else {
331  ntowers[ieta0]++;
332  }
333  }
334 
335  for (auto const& gt : geomtowers_) {
336  int it = gt.first;
337 
338  int vi = it - 1;
339 
340  if (it < 0)
341  vi = nEtaTow_ + it;
342 
343  double e1 = emean_[it];
344  double e2 = emean2[it];
345  int nt = gt.second - ntowersWithJets_[it];
346 
347  if (vi < nEtaTow_) {
348  vntow_[vi] = nt;
349  }
350 
351  LogDebug("PileUpSubtractor") << " ieta: " << it << " number of towers: " << nt << " e1: " << e1 << " e2: " << e2
352  << "\n";
353 
354  if (nt > 0) {
355  if (postOrphan_) {
356  if (nt > (int)minimumTowersFraction_ * (gt.second)) {
357  emean_[it] = e1 / (double)nt;
358  double eee = e2 / (double)nt - e1 * e1 / (double)(nt * nt);
359  if (eee < 0.)
360  eee = 0.;
361  esigma_[it] = nSigmaPU_ * sqrt(eee);
362 
363  uint32_t numToCheck = std::min(int(eTop4_[it].size()), nt - (int)minimumTowersFraction_ * (gt.second));
364 
365  for (unsigned int lI = 0; lI < numToCheck; ++lI) {
366  if (eTop4_[it][lI] >= emean_[it] + towSigmaCut_ * esigma_[it] && towSigmaCut_ > 0) {
367  e1 -= eTop4_[it][lI];
368  nt -= 1;
369  } else
370  break;
371  }
372 
373  if (e1 < .000000001)
374  e1 = 0;
375  }
376  }
377 
378  emean_[it] = e1 / (double)nt;
379  double eee = e2 / nt - e1 * e1 / (nt * nt);
380  if (eee < 0.)
381  eee = 0.;
382  esigma_[it] = nSigmaPU_ * sqrt(eee);
383 
384  double etaWidth = std::abs(hi::etaedge[abs(it)] - hi::etaedge[abs(it) - 1]);
385 
386  int sign = (it < 0) ? -1 : 1;
387  auto absIt = std::abs(it);
388  if (it < 0) {
389  etaEdgeLow_.push_back(sign * hi::etaedge[absIt]);
390  etaEdgeHi_.push_back(sign * hi::etaedge[absIt - 1]);
391  } else {
392  etaEdgeHi_.push_back(sign * hi::etaedge[absIt]);
393  etaEdgeLow_.push_back(sign * hi::etaedge[absIt - 1]);
394  }
395 
396  if (vi < nEtaTow_) {
397  vmean1_[vi] = emean_[it];
398  vrho1_[vi] = emean_[it] / (etaWidth * (2. * M_PI / (double)vngeom_[vi]));
399  rho_.push_back(vrho1_[vi]);
400  rhoM_.push_back(0);
401  vrms1_[vi] = esigma_[it];
402  if (vngeom_[vi] == vntow_[vi]) {
403  vmean0_[vi] = emean_[it];
404  vrho0_[vi] = emean_[it] / (etaWidth * (2. * M_PI / (double)vngeom_[vi]));
405  vrms0_[vi] = esigma_[it];
406  }
407  rhoExtra_.push_back(vrho0_[vi]);
408  nTow_.push_back(vntow_[vi]);
409  }
410  } else {
411  emean_[it] = 0.;
412  esigma_[it] = 0.;
413  }
414  LogDebug("PileUpSubtractor") << " ieta: " << it << " Pedestals: " << emean_[it] << " " << esigma_[it] << "\n";
415  }
416 
417  postOrphan_ = false;
418 }
size
Write out results.
const double nSigmaPU_
std::vector< double > rhoExtra_
std::array< float, nEtaTow_ > vrms0_
std::array< float, nEtaTow_ > vrho1_
std::vector< double > etaEdgeHi_
std::vector< double > rho_
std::map< int, std::array< double, 4 > > eTop4_
std::array< float, nEtaTow_ > vmean0_
std::array< float, nEtaTow_ > vmean1_
static constexpr int nEtaTow_
double et(double vtxZ) const
Definition: CaloTower.h:150
std::vector< double > etaEdgeLow_
std::array< float, nEtaTow_ > vrms1_
std::vector< const CaloTower * > inputs_
const double towSigmaCut_
T sqrt(T t)
Definition: SSEVec.h:19
std::array< int, nEtaTow_ > vngeom_
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
int nt
Definition: AMPTWrapper.h:42
#define M_PI
const int initialValue_
std::array< int, nEtaTow_ > vntow_
int ieta() const
Definition: CaloTower.h:200
std::map< int, int > geomtowers_
std::map< int, double > emean_
std::vector< int > nTow_
std::vector< double > rhoM_
std::map< int, double > esigma_
std::map< int, int > ntowersWithJets_
constexpr std::array< double, 42 > etaedge
Definition: HITowerHelper.h:5
const double minimumTowersFraction_
std::array< float, nEtaTow_ > vrho0_
#define LogDebug(id)

◆ fillDescriptions()

void HiPuRhoProducer::fillDescriptions ( edm::ConfigurationDescriptions descriptions)
static

Definition at line 561 of file HiPuRhoProducer.cc.

References edm::ConfigurationDescriptions::add(), submitPVResolutionJobs::desc, and ProducerED_cfi::InputTag.

561  {
563  desc.add<edm::InputTag>("src", edm::InputTag("PFTowers"));
564  desc.add<int>("medianWindowWidth", 2);
565  desc.add<double>("towSigmaCut", 5.0);
566  desc.add<double>("puPtMin", 15.0);
567  desc.add<double>("rParam", 0.3);
568  desc.add<double>("nSigmaPU", 1.0);
569  desc.add<double>("radiusPU", 0.5);
570  desc.add<double>("minimumTowersFraction", 0.7);
571  desc.add<bool>("dropZeroTowers", true);
572  descriptions.add("hiPuRhoProducer", desc);
573 }
void add(std::string const &label, ParameterSetDescription const &psetDescription)

◆ inputTowers()

void HiPuRhoProducer::inputTowers ( )
privatevirtual

Definition at line 219 of file HiPuRhoProducer.cc.

References geometryDiff::epsilon, fjInputs_, input, inputs_, and edm::isNotFinite().

Referenced by produce().

219  {
220  int index = -1;
221  for (auto const& input : inputs_) {
222  index++;
223 
224  if (edm::isNotFinite(input->pt()))
225  continue;
226  if (input->pt() < 100 * std::numeric_limits<double>::epsilon())
227  continue;
228 
229  fjInputs_.push_back(fastjet::PseudoJet(input->px(), input->py(), input->pz(), input->energy()));
230  fjInputs_.back().set_user_index(index);
231  }
232 }
constexpr bool isNotFinite(T x)
Definition: isFinite.h:9
std::vector< fastjet::PseudoJet > fjInputs_
static std::string const input
Definition: EdmProvDump.cc:50
std::vector< const CaloTower * > inputs_

◆ produce()

void HiPuRhoProducer::produce ( edm::Event iEvent,
const edm::EventSetup iSetup 
)
overrideprivate

Definition at line 172 of file HiPuRhoProducer.cc.

References calculateOrphanInput(), calculatePedestal(), caloTowerToken_, etaEdgeHi_, etaEdgeLow_, etaEdges_, fjClusterSeq_, fjInputs_, fjJetDefinition_, fjJets_, fjOriginalInputs_, mps_fire::i, ietamax_, ietamin_, iEvent, input, inputs_, inputTowers(), nTow_, ntowersWithJets_, puPtMin_, putRho(), rho_, rhoExtra_, rhoM_, rParam_, setInitialValue_, setupGeometryMap(), subtractPedestal(), towExcludeEta_, towExcludePhi_, and towExcludePt_.

172  {
173  setupGeometryMap(iEvent, iSetup);
174 
175  for (int i = ietamin_; i < ietamax_ + 1; i++) {
176  ntowersWithJets_[i] = 0;
177  }
178 
179  auto const& inputView = iEvent.get(caloTowerToken_);
180  inputs_.reserve(inputView.size());
181  for (auto const& input : inputView)
182  inputs_.push_back(&input);
183 
184  fjInputs_.reserve(inputs_.size());
185  inputTowers();
187  setInitialValue_ = true;
190 
191  fjJetDefinition_ = std::make_shared<fastjet::JetDefinition>(fastjet::antikt_algorithm, rParam_);
192  fjClusterSeq_ = std::make_shared<fastjet::ClusterSequence>(fjInputs_, *fjJetDefinition_);
193  fjJets_ = fastjet::sorted_by_pt(fjClusterSeq_->inclusive_jets(puPtMin_));
194 
195  etaEdgeLow_.clear();
196  etaEdgeHi_.clear();
197  etaEdges_.clear();
198 
199  rho_.clear();
200  rhoExtra_.clear();
201  rhoM_.clear();
202  nTow_.clear();
203 
204  towExcludePt_.clear();
205  towExcludePhi_.clear();
206  towExcludeEta_.clear();
207 
208  setInitialValue_ = false;
209  std::vector<fastjet::PseudoJet> orphanInput;
210  calculateOrphanInput(orphanInput);
211  calculatePedestal(orphanInput);
212  putRho(iEvent, iSetup);
213 
214  inputs_.clear();
215  fjInputs_.clear();
216  fjJets_.clear();
217 }
std::vector< double > towExcludePt_
std::vector< double > etaEdges_
std::vector< double > rhoExtra_
virtual void inputTowers()
std::vector< double > towExcludeEta_
const edm::EDGetTokenT< CaloTowerCollection > caloTowerToken_
const double rParam_
std::vector< double > etaEdgeHi_
std::vector< fastjet::PseudoJet > fjInputs_
std::vector< double > rho_
std::vector< double > towExcludePhi_
static std::string const input
Definition: EdmProvDump.cc:50
std::vector< double > etaEdgeLow_
virtual void subtractPedestal(std::vector< fastjet::PseudoJet > &coll)
int iEvent
Definition: GenABIO.cc:224
std::vector< const CaloTower * > inputs_
virtual void calculateOrphanInput(std::vector< fastjet::PseudoJet > &orphanInput)
virtual void setupGeometryMap(edm::Event &iEvent, const edm::EventSetup &iSetup)
virtual void putRho(edm::Event &iEvent, const edm::EventSetup &iSetup)
const double puPtMin_
ClusterSequencePtr fjClusterSeq_
std::vector< fastjet::PseudoJet > fjJets_
std::vector< int > nTow_
std::vector< fastjet::PseudoJet > fjOriginalInputs_
std::vector< double > rhoM_
JetDefPtr fjJetDefinition_
std::map< int, int > ntowersWithJets_
virtual void calculatePedestal(std::vector< fastjet::PseudoJet > const &coll)

◆ putRho()

void HiPuRhoProducer::putRho ( edm::Event iEvent,
const edm::EventSetup iSetup 
)
virtual

Definition at line 500 of file HiPuRhoProducer.cc.

References etaEdgeHi_, etaEdgeLow_, mps_fire::i, iEvent, medianWindowWidth_, eostools::move(), nTow_, eventshapeDQM_cfi::order, rho_, rhoExtra_, rhoM_, findQualityFiles::size, jetUpdater_cfi::sort, towExcludeEta_, towExcludePhi_, and towExcludePt_.

Referenced by produce().

500  {
501  std::size_t size = etaEdgeLow_.size();
502 
503  std::vector<std::pair<std::size_t, double>> order;
504  for (std::size_t i = 0; i < size; ++i) {
505  order.emplace_back(i, etaEdgeLow_[i]);
506  }
507 
508  std::sort(
509  order.begin(), order.end(), [](auto const& pair0, auto const& pair1) { return pair0.second < pair1.second; });
510 
511  std::vector<double> sortedEtaEdgeLow(size);
512  std::vector<double> sortedEtaEdgeHigh(size);
513 
514  auto mapToRhoOut = std::make_unique<std::vector<double>>(size);
515  auto mapToRhoExtraOut = std::make_unique<std::vector<double>>(size);
516  auto mapToRhoMOut = std::make_unique<std::vector<double>>(size);
517  auto mapToNTowOut = std::make_unique<std::vector<int>>(size);
518 
519  for (std::size_t i = 0; i < size; ++i) {
520  auto const& index = order[i].first;
521 
522  sortedEtaEdgeLow[i] = etaEdgeLow_[index];
523  sortedEtaEdgeHigh[i] = etaEdgeHi_[index];
524 
525  (*mapToRhoOut)[i] = rho_[index];
526  (*mapToRhoExtraOut)[i] = rhoExtra_[index];
527  (*mapToRhoMOut)[i] = rhoM_[index];
528  (*mapToNTowOut)[i] = nTow_[index];
529  }
530 
531  auto mapToRhoMedianOut = std::make_unique<std::vector<double>>(size);
532 
533  for (uint32_t index = medianWindowWidth_; index < size - medianWindowWidth_; ++index) {
534  auto centre = mapToRhoOut->begin() + index;
535  std::vector<float> rhoRange(centre - medianWindowWidth_, centre + medianWindowWidth_);
536  std::nth_element(rhoRange.begin(), rhoRange.begin() + medianWindowWidth_, rhoRange.end());
537  (*mapToRhoMedianOut)[index] = rhoRange[medianWindowWidth_];
538  }
539 
540  auto mapEtaRangesOut = std::make_unique<std::vector<double>>();
541 
542  mapEtaRangesOut->assign(sortedEtaEdgeLow.begin(), sortedEtaEdgeLow.end());
543  mapEtaRangesOut->push_back(sortedEtaEdgeHigh.back());
544 
545  auto mapToTowExcludePtOut = std::make_unique<std::vector<double>>(std::move(towExcludePt_));
546  auto mapToTowExcludePhiOut = std::make_unique<std::vector<double>>(std::move(towExcludePhi_));
547  auto mapToTowExcludeEtaOut = std::make_unique<std::vector<double>>(std::move(towExcludeEta_));
548 
549  iEvent.put(std::move(mapEtaRangesOut), "mapEtaEdges");
550  iEvent.put(std::move(mapToRhoOut), "mapToRho");
551  iEvent.put(std::move(mapToRhoMedianOut), "mapToRhoMedian");
552  iEvent.put(std::move(mapToRhoExtraOut), "mapToRhoExtra");
553  iEvent.put(std::move(mapToRhoMOut), "mapToRhoM");
554  iEvent.put(std::move(mapToNTowOut), "mapToNTow");
555  iEvent.put(std::move(mapToTowExcludePtOut), "mapToTowExcludePt");
556  iEvent.put(std::move(mapToTowExcludePhiOut), "mapToTowExcludePhi");
557  iEvent.put(std::move(mapToTowExcludeEtaOut), "mapToTowExcludeEta");
558 }
size
Write out results.
std::vector< double > towExcludePt_
std::vector< double > rhoExtra_
std::vector< double > towExcludeEta_
std::vector< double > etaEdgeHi_
std::vector< double > rho_
std::vector< double > towExcludePhi_
std::vector< double > etaEdgeLow_
int iEvent
Definition: GenABIO.cc:224
std::vector< int > nTow_
std::vector< double > rhoM_
const int medianWindowWidth_
def move(src, dest)
Definition: eostools.py:511

◆ setupGeometryMap()

void HiPuRhoProducer::setupGeometryMap ( edm::Event iEvent,
const edm::EventSetup iSetup 
)
virtual

Definition at line 234 of file HiPuRhoProducer.cc.

References caloGeometryToken_, PV3DBase< T, PVType, FrameType >::eta(), geo_, geomtowers_, edm::EventSetup::getData(), CaloGeometry::getPosition(), CaloGeometry::getValidDetIds(), submitPVValidationJobs::gt, DetId::Hcal, HcalDetId::ieta(), ietamax_, ietamin_, HcalDetId::iphi(), ALPAKA_ACCELERATOR_NAMESPACE::vertexFinder::it, LogDebug, nEtaTow_, PV3DBase< T, PVType, FrameType >::phi(), towermap_, and vngeom_.

Referenced by produce().

234  {
235  LogDebug("PileUpSubtractor") << "The subtractor setting up geometry...\n";
236  const auto& pG = iSetup.getData(caloGeometryToken_);
237  geo_ = &pG;
238  std::vector<DetId> alldid = geo_->getValidDetIds();
239  int ietaold = -10000;
240  ietamax_ = -10000;
241  ietamin_ = 10000;
242  towermap_.clear();
243 
244  for (auto const& did : alldid) {
245  if (did.det() == DetId::Hcal) {
246  HcalDetId hid = HcalDetId(did);
247  towermap_.push_back({hid.ieta(), hid.iphi(), geo_->getPosition(did).eta(), geo_->getPosition(did).phi()});
248  if (hid.ieta() != ietaold) {
249  ietaold = hid.ieta();
250  geomtowers_[hid.ieta()] = 1;
251  if (hid.ieta() > ietamax_)
252  ietamax_ = hid.ieta();
253  if (hid.ieta() < ietamin_)
254  ietamin_ = hid.ieta();
255  } else {
256  geomtowers_[hid.ieta()]++;
257  }
258  }
259  }
260 
261  for (auto const& gt : geomtowers_) {
262  int it = gt.first;
263  int vi = it - 1;
264 
265  if (it < 0)
266  vi = nEtaTow_ + it;
267 
268  if (vi < nEtaTow_)
269  vngeom_[vi] = gt.second;
270  }
271 }
const edm::ESGetToken< CaloGeometry, CaloGeometryRecord > caloGeometryToken_
T const & getData(const ESGetToken< T, R > &iToken) const noexcept(false)
Definition: EventSetup.h:119
std::vector< EtaPhiTower > towermap_
Geom::Phi< T > phi() const
Definition: PV3DBase.h:66
T eta() const
Definition: PV3DBase.h:73
static constexpr int nEtaTow_
GlobalPoint getPosition(const DetId &id) const
Get the position of a given detector id.
Definition: CaloGeometry.cc:50
constexpr int ieta() const
get the cell ieta
Definition: HcalDetId.h:155
std::array< int, nEtaTow_ > vngeom_
CaloGeometry const * geo_
std::vector< DetId > getValidDetIds() const
Get the list of all valid detector ids.
Definition: CaloGeometry.cc:75
std::map< int, int > geomtowers_
constexpr int iphi() const
get the cell iphi
Definition: HcalDetId.h:157
#define LogDebug(id)

◆ subtractPedestal()

void HiPuRhoProducer::subtractPedestal ( std::vector< fastjet::PseudoJet > &  coll)
virtual

Definition at line 420 of file HiPuRhoProducer.cc.

References dropZeroTowers_, emean_, esigma_, CaloTower::et(), CaloTower::ieta(), inputs_, ALPAKA_ACCELERATOR_NAMESPACE::vertexFinder::it, and LogDebug.

Referenced by produce().

420  {
421  LogDebug("PileUpSubtractor") << "The subtractor subtracting pedestals...\n";
422 
423  std::vector<fastjet::PseudoJet> newcoll;
424  for (auto& input_object : coll) {
425  int index = input_object.user_index();
426  CaloTower const* itow = inputs_[index];
427  int it = itow->ieta();
428 
429  double original_Et = itow->et();
430  double etnew = original_Et - emean_.at(it) - esigma_.at(it);
431  float mScale = etnew / input_object.Et();
432  if (etnew < 0.)
433  mScale = 0.;
434 
436  input_object.px() * mScale, input_object.py() * mScale, input_object.pz() * mScale, input_object.e() * mScale);
437 
438  input_object.reset_momentum(towP4.px(), towP4.py(), towP4.pz(), towP4.energy());
439  input_object.set_user_index(index);
440 
441  if (etnew > 0. && dropZeroTowers_)
442  newcoll.push_back(input_object);
443  }
444 
445  if (dropZeroTowers_)
446  coll = newcoll;
447 }
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > XYZTLorentzVectorD
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:14
double et(double vtxZ) const
Definition: CaloTower.h:150
std::vector< const CaloTower * > inputs_
int ieta() const
Definition: CaloTower.h:200
const bool dropZeroTowers_
std::map< int, double > emean_
std::map< int, double > esigma_
#define LogDebug(id)

Member Data Documentation

◆ caloGeometryToken_

const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> HiPuRhoProducer::caloGeometryToken_
private

Definition at line 100 of file HiPuRhoProducer.cc.

Referenced by setupGeometryMap().

◆ caloTowerToken_

const edm::EDGetTokenT<CaloTowerCollection> HiPuRhoProducer::caloTowerToken_
private

Definition at line 99 of file HiPuRhoProducer.cc.

Referenced by produce().

◆ dropZeroTowers_

const bool HiPuRhoProducer::dropZeroTowers_
private

Definition at line 91 of file HiPuRhoProducer.cc.

Referenced by subtractPedestal().

◆ emean_

std::map<int, double> HiPuRhoProducer::emean_
private

Definition at line 141 of file HiPuRhoProducer.cc.

Referenced by calculatePedestal(), and subtractPedestal().

◆ esigma_

std::map<int, double> HiPuRhoProducer::esigma_
private

Definition at line 140 of file HiPuRhoProducer.cc.

Referenced by calculatePedestal(), and subtractPedestal().

◆ etaEdgeHi_

std::vector<double> HiPuRhoProducer::etaEdgeHi_
private

Definition at line 114 of file HiPuRhoProducer.cc.

Referenced by calculatePedestal(), produce(), and putRho().

◆ etaEdgeLow_

std::vector<double> HiPuRhoProducer::etaEdgeLow_
private

Definition at line 113 of file HiPuRhoProducer.cc.

Referenced by calculatePedestal(), produce(), and putRho().

◆ etaEdges_

std::vector<double> HiPuRhoProducer::etaEdges_
private

Definition at line 115 of file HiPuRhoProducer.cc.

Referenced by produce().

◆ eTop4_

std::map<int, std::array<double, 4> > HiPuRhoProducer::eTop4_
private

Definition at line 142 of file HiPuRhoProducer.cc.

Referenced by calculatePedestal().

◆ fjClusterSeq_

ClusterSequencePtr HiPuRhoProducer::fjClusterSeq_
private

Definition at line 127 of file HiPuRhoProducer.cc.

Referenced by produce().

◆ fjInputs_

std::vector<fastjet::PseudoJet> HiPuRhoProducer::fjInputs_
private

Definition at line 130 of file HiPuRhoProducer.cc.

Referenced by calculateOrphanInput(), inputTowers(), and produce().

◆ fjJetDefinition_

JetDefPtr HiPuRhoProducer::fjJetDefinition_
private

Definition at line 128 of file HiPuRhoProducer.cc.

Referenced by produce().

◆ fjJets_

std::vector<fastjet::PseudoJet> HiPuRhoProducer::fjJets_
private

Definition at line 131 of file HiPuRhoProducer.cc.

Referenced by calculateOrphanInput(), and produce().

◆ fjOriginalInputs_

std::vector<fastjet::PseudoJet> HiPuRhoProducer::fjOriginalInputs_
private

Definition at line 132 of file HiPuRhoProducer.cc.

Referenced by calculateOrphanInput(), and produce().

◆ geo_

CaloGeometry const* HiPuRhoProducer::geo_ = nullptr
private

Definition at line 134 of file HiPuRhoProducer.cc.

Referenced by setupGeometryMap().

◆ geomtowers_

std::map<int, int> HiPuRhoProducer::geomtowers_
private

Definition at line 139 of file HiPuRhoProducer.cc.

Referenced by calculateOrphanInput(), calculatePedestal(), and setupGeometryMap().

◆ ietamax_

int HiPuRhoProducer::ietamax_
private

Definition at line 136 of file HiPuRhoProducer.cc.

Referenced by calculatePedestal(), produce(), and setupGeometryMap().

◆ ietamin_

int HiPuRhoProducer::ietamin_
private

Definition at line 137 of file HiPuRhoProducer.cc.

Referenced by calculatePedestal(), produce(), and setupGeometryMap().

◆ initialValue_

const int HiPuRhoProducer::initialValue_ = -99
private

Definition at line 101 of file HiPuRhoProducer.cc.

Referenced by calculatePedestal().

◆ inputs_

std::vector<const CaloTower*> HiPuRhoProducer::inputs_
private

◆ medianWindowWidth_

const int HiPuRhoProducer::medianWindowWidth_
private

Definition at line 92 of file HiPuRhoProducer.cc.

Referenced by putRho().

◆ minimumTowersFraction_

const double HiPuRhoProducer::minimumTowersFraction_
private

Definition at line 93 of file HiPuRhoProducer.cc.

Referenced by calculateOrphanInput(), and calculatePedestal().

◆ nEtaTow_

constexpr int HiPuRhoProducer::nEtaTow_ = 82
staticprivate

Definition at line 102 of file HiPuRhoProducer.cc.

Referenced by calculatePedestal(), and setupGeometryMap().

◆ nSigmaPU_

const double HiPuRhoProducer::nSigmaPU_
private

Definition at line 94 of file HiPuRhoProducer.cc.

Referenced by calculatePedestal().

◆ nTow_

std::vector<int> HiPuRhoProducer::nTow_
private

Definition at line 120 of file HiPuRhoProducer.cc.

Referenced by calculatePedestal(), produce(), and putRho().

◆ ntowersWithJets_

std::map<int, int> HiPuRhoProducer::ntowersWithJets_
private

Definition at line 138 of file HiPuRhoProducer.cc.

Referenced by calculateOrphanInput(), calculatePedestal(), and produce().

◆ postOrphan_

bool HiPuRhoProducer::postOrphan_ = false
private

Definition at line 88 of file HiPuRhoProducer.cc.

Referenced by calculateOrphanInput(), and calculatePedestal().

◆ puPtMin_

const double HiPuRhoProducer::puPtMin_
private

Definition at line 95 of file HiPuRhoProducer.cc.

Referenced by calculateOrphanInput(), and produce().

◆ radiusPU_

const double HiPuRhoProducer::radiusPU_
private

Definition at line 96 of file HiPuRhoProducer.cc.

Referenced by calculateOrphanInput().

◆ rho_

std::vector<double> HiPuRhoProducer::rho_
private

Definition at line 117 of file HiPuRhoProducer.cc.

Referenced by calculatePedestal(), produce(), and putRho().

◆ rhoExtra_

std::vector<double> HiPuRhoProducer::rhoExtra_
private

Definition at line 118 of file HiPuRhoProducer.cc.

Referenced by calculatePedestal(), produce(), and putRho().

◆ rhoM_

std::vector<double> HiPuRhoProducer::rhoM_
private

Definition at line 119 of file HiPuRhoProducer.cc.

Referenced by calculatePedestal(), produce(), and putRho().

◆ rParam_

const double HiPuRhoProducer::rParam_
private

Definition at line 97 of file HiPuRhoProducer.cc.

Referenced by produce().

◆ setInitialValue_

bool HiPuRhoProducer::setInitialValue_
private

Definition at line 89 of file HiPuRhoProducer.cc.

Referenced by calculatePedestal(), and produce().

◆ towermap_

std::vector<EtaPhiTower> HiPuRhoProducer::towermap_
private

Definition at line 145 of file HiPuRhoProducer.cc.

Referenced by calculateOrphanInput(), and setupGeometryMap().

◆ towExcludeEta_

std::vector<double> HiPuRhoProducer::towExcludeEta_
private

Definition at line 124 of file HiPuRhoProducer.cc.

Referenced by calculateOrphanInput(), produce(), and putRho().

◆ towExcludePhi_

std::vector<double> HiPuRhoProducer::towExcludePhi_
private

Definition at line 123 of file HiPuRhoProducer.cc.

Referenced by calculateOrphanInput(), produce(), and putRho().

◆ towExcludePt_

std::vector<double> HiPuRhoProducer::towExcludePt_
private

Definition at line 122 of file HiPuRhoProducer.cc.

Referenced by calculateOrphanInput(), produce(), and putRho().

◆ towSigmaCut_

const double HiPuRhoProducer::towSigmaCut_
private

Definition at line 98 of file HiPuRhoProducer.cc.

Referenced by calculatePedestal().

◆ vmean0_

std::array<float, nEtaTow_> HiPuRhoProducer::vmean0_
private

Definition at line 106 of file HiPuRhoProducer.cc.

Referenced by calculatePedestal().

◆ vmean1_

std::array<float, nEtaTow_> HiPuRhoProducer::vmean1_
private

Definition at line 109 of file HiPuRhoProducer.cc.

Referenced by calculatePedestal().

◆ vngeom_

std::array<int, nEtaTow_> HiPuRhoProducer::vngeom_
private

Definition at line 104 of file HiPuRhoProducer.cc.

Referenced by calculatePedestal(), and setupGeometryMap().

◆ vntow_

std::array<int, nEtaTow_> HiPuRhoProducer::vntow_
private

Definition at line 105 of file HiPuRhoProducer.cc.

Referenced by calculatePedestal().

◆ vrho0_

std::array<float, nEtaTow_> HiPuRhoProducer::vrho0_
private

Definition at line 108 of file HiPuRhoProducer.cc.

Referenced by calculatePedestal().

◆ vrho1_

std::array<float, nEtaTow_> HiPuRhoProducer::vrho1_
private

Definition at line 111 of file HiPuRhoProducer.cc.

Referenced by calculatePedestal().

◆ vrms0_

std::array<float, nEtaTow_> HiPuRhoProducer::vrms0_
private

Definition at line 107 of file HiPuRhoProducer.cc.

Referenced by calculatePedestal().

◆ vrms1_

std::array<float, nEtaTow_> HiPuRhoProducer::vrms1_
private

Definition at line 110 of file HiPuRhoProducer.cc.

Referenced by calculatePedestal().