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

std::vector< HcalDetIdallgeomid_
 
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 145 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 149 of file HiPuRhoProducer.cc.

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

451  {
452  LogDebug("PileUpSubtractor") << "The subtractor calculating orphan input...\n";
453 
455 
456  std::vector<int> jettowers; // vector of towers indexed by "user_index"
457  std::map<std::pair<int, int>, int> excludedTowers; // map of excluded ieta, iphi values
458  for (auto const& pseudojetTMP : fjJets_) {
459  EtaPhi jet_etaphi(pseudojetTMP.eta(), pseudojetTMP.phi());
460 
461  if (pseudojetTMP.perp() < puPtMin_)
462  continue;
463 
464  for (auto const& im : towermap_) {
465  double dr2 = reco::deltaR2(im.eta, im.phi, jet_etaphi.first, jet_etaphi.second);
466  if (dr2 < radiusPU_ * radiusPU_ && !excludedTowers[std::pair(im.ieta, im.iphi)] &&
467  (geomtowers_[im.ieta] - ntowersWithJets_[im.ieta]) > minimumTowersFraction_ * geomtowers_[im.ieta]) {
468  ntowersWithJets_[im.ieta]++;
469  excludedTowers[std::pair(im.ieta, im.iphi)] = 1;
470  }
471  }
472 
473  for (auto const& input : fjInputs_) {
474  int index = input.user_index();
475  const CaloTower* originalTower = inputs_[index];
476  int ie = originalTower->ieta();
477  int ip = originalTower->iphi();
478 
479  if (excludedTowers[std::pair<int, int>(ie, ip)]) {
480  jettowers.push_back(index);
481  }
482  }
483  }
484  // Create a new collections from the towers not included in jets
485  for (auto const& input : fjInputs_) {
486  int index = input.user_index();
487  const CaloTower* originalTower = inputs_[index];
488  auto itjet = std::find(jettowers.begin(), jettowers.end(), index);
489  if (itjet == jettowers.end()) {
490  orphanInput.emplace_back(originalTower->px(), originalTower->py(), originalTower->pz(), originalTower->energy());
491  orphanInput.back().set_user_index(index);
492  } else {
493  towExcludePt_.push_back(originalTower->pt());
494  towExcludePhi_.push_back(originalTower->phi());
495  towExcludeEta_.push_back(originalTower->eta());
496  }
497  }
498 
499  postOrphan_ = true;
500 }
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 275 of file HiPuRhoProducer.cc.

References funct::abs(), StorageManager_cfg::e1, emean_, esigma_, CaloTower::et(), hi::etaedge, etaEdgeHi_, etaEdgeLow_, photons_cff::etaWidth, eTop4_, geomtowers_, submitPVValidationJobs::gt, mps_fire::i, CaloTower::ieta(), ietamax_, ietamin_, initialValue_, inputs_, 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().

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

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

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

◆ inputTowers()

void HiPuRhoProducer::inputTowers ( )
privatevirtual

Definition at line 220 of file HiPuRhoProducer.cc.

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

Referenced by produce().

220  {
221  int index = -1;
222  for (auto const& input : inputs_) {
223  index++;
224 
225  if (edm::isNotFinite(input->pt()))
226  continue;
227  if (input->pt() < 100 * std::numeric_limits<double>::epsilon())
228  continue;
229 
230  fjInputs_.push_back(fastjet::PseudoJet(input->px(), input->py(), input->pz(), input->energy()));
231  fjInputs_.back().set_user_index(index);
232  }
233 }
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 173 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_.

173  {
174  setupGeometryMap(iEvent, iSetup);
175 
176  for (int i = ietamin_; i < ietamax_ + 1; i++) {
177  ntowersWithJets_[i] = 0;
178  }
179 
180  auto const& inputView = iEvent.get(caloTowerToken_);
181  inputs_.reserve(inputView.size());
182  for (auto const& input : inputView)
183  inputs_.push_back(&input);
184 
185  fjInputs_.reserve(inputs_.size());
186  inputTowers();
188  setInitialValue_ = true;
191 
192  fjJetDefinition_ = std::make_shared<fastjet::JetDefinition>(fastjet::antikt_algorithm, rParam_);
193  fjClusterSeq_ = std::make_shared<fastjet::ClusterSequence>(fjInputs_, *fjJetDefinition_);
194  fjJets_ = fastjet::sorted_by_pt(fjClusterSeq_->inclusive_jets(puPtMin_));
195 
196  etaEdgeLow_.clear();
197  etaEdgeHi_.clear();
198  etaEdges_.clear();
199 
200  rho_.clear();
201  rhoExtra_.clear();
202  rhoM_.clear();
203  nTow_.clear();
204 
205  towExcludePt_.clear();
206  towExcludePhi_.clear();
207  towExcludeEta_.clear();
208 
209  setInitialValue_ = false;
210  std::vector<fastjet::PseudoJet> orphanInput;
211  calculateOrphanInput(orphanInput);
212  calculatePedestal(orphanInput);
213  putRho(iEvent, iSetup);
214 
215  inputs_.clear();
216  fjInputs_.clear();
217  fjJets_.clear();
218 }
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 502 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().

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

References allgeomid_, 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(), LogDebug, nEtaTow_, PV3DBase< T, PVType, FrameType >::phi(), towermap_, and vngeom_.

Referenced by produce().

235  {
236  LogDebug("PileUpSubtractor") << "The subtractor setting up geometry...\n";
237  const auto& pG = iSetup.getData(caloGeometryToken_);
238  geo_ = &pG;
239  std::vector<DetId> alldid = geo_->getValidDetIds();
240  int ietaold = -10000;
241  ietamax_ = -10000;
242  ietamin_ = 10000;
243  towermap_.clear();
244 
245  for (auto const& did : alldid) {
246  if (did.det() == DetId::Hcal) {
247  HcalDetId hid = HcalDetId(did);
248  allgeomid_.push_back(did);
249  towermap_.push_back({hid.ieta(), hid.iphi(), geo_->getPosition(did).eta(), geo_->getPosition(did).phi()});
250  if (hid.ieta() != ietaold) {
251  ietaold = hid.ieta();
252  geomtowers_[hid.ieta()] = 1;
253  if (hid.ieta() > ietamax_)
254  ietamax_ = hid.ieta();
255  if (hid.ieta() < ietamin_)
256  ietamin_ = hid.ieta();
257  } else {
258  geomtowers_[hid.ieta()]++;
259  }
260  }
261  }
262 
263  for (auto const& gt : geomtowers_) {
264  int it = gt.first;
265  int vi = it - 1;
266 
267  if (it < 0)
268  vi = nEtaTow_ + it;
269 
270  if (vi < nEtaTow_)
271  vngeom_[vi] = gt.second;
272  }
273 }
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< HcalDetId > allgeomid_
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 422 of file HiPuRhoProducer.cc.

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

Referenced by produce().

422  {
423  LogDebug("PileUpSubtractor") << "The subtractor subtracting pedestals...\n";
424 
425  std::vector<fastjet::PseudoJet> newcoll;
426  for (auto& input_object : coll) {
427  int index = input_object.user_index();
428  CaloTower const* itow = inputs_[index];
429  int it = itow->ieta();
430 
431  double original_Et = itow->et();
432  double etnew = original_Et - emean_.at(it) - esigma_.at(it);
433  float mScale = etnew / input_object.Et();
434  if (etnew < 0.)
435  mScale = 0.;
436 
438  input_object.px() * mScale, input_object.py() * mScale, input_object.pz() * mScale, input_object.e() * mScale);
439 
440  input_object.reset_momentum(towP4.px(), towP4.py(), towP4.pz(), towP4.energy());
441  input_object.set_user_index(index);
442 
443  if (etnew > 0. && dropZeroTowers_)
444  newcoll.push_back(input_object);
445  }
446 
447  if (dropZeroTowers_)
448  coll = newcoll;
449 }
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

◆ allgeomid_

std::vector<HcalDetId> HiPuRhoProducer::allgeomid_
private

Definition at line 135 of file HiPuRhoProducer.cc.

Referenced by setupGeometryMap().

◆ 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 142 of file HiPuRhoProducer.cc.

Referenced by calculatePedestal(), and subtractPedestal().

◆ esigma_

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

Definition at line 141 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 143 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 140 of file HiPuRhoProducer.cc.

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

◆ ietamax_

int HiPuRhoProducer::ietamax_
private

Definition at line 137 of file HiPuRhoProducer.cc.

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

◆ ietamin_

int HiPuRhoProducer::ietamin_
private

Definition at line 138 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 139 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 146 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().