CMS 3D CMS Logo
HiPuRhoProducer.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: RecoHI/HiJetAlgos/plugins/HiPuRhoProducer
4 // Class: HiPuRhoProducer
5 //
11 //
12 // Original Author: Chris McGinn - in Style of HiFJRhoProducer
13 // Created: Mon, 29 May 2017
14 //
15 //
16 #include "DataFormats/CaloTowers/interface/CaloTower.h"
17 #include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
18 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
19 #include "DataFormats/Candidate/interface/Candidate.h"
20 #include "DataFormats/Common/interface/Handle.h"
21 #include "DataFormats/Common/interface/View.h"
22 #include "DataFormats/DetId/interface/DetId.h"
23 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
24 #include "DataFormats/Math/interface/LorentzVector.h"
25 #include "DataFormats/Math/interface/deltaR.h"
26 #include "FWCore/Framework/interface/ESHandle.h"
27 #include "FWCore/Framework/interface/Event.h"
28 #include "FWCore/Framework/interface/MakerMacros.h"
29 #include "FWCore/MessageLogger/interface/MessageLogger.h"
30 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
31 #include "FWCore/ParameterSet/interface/ParameterSet.h"
32 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
33 #include "FWCore/ParameterSet/interface/ParameterSetDescriptionFiller.h"
34 #include "FWCore/Utilities/interface/Exception.h"
35 #include "FWCore/Utilities/interface/isFinite.h"
36 #include "Geometry/Records/interface/CaloGeometryRecord.h"
37 #include "DataFormats/Candidate/interface/CandidateFwd.h"
38 #include "DataFormats/Common/interface/Ptr.h"
39 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
40 #include "FWCore/Framework/interface/stream/EDProducer.h"
41 #include "FWCore/Framework/interface/EventSetup.h"
42 #include "FWCore/Utilities/interface/EDGetToken.h"
43 #include "FWCore/Utilities/interface/InputTag.h"
44 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
45 #include "RecoHI/HiJetAlgos/plugins/HITowerHelper.h"
46 
47 #include "fastjet/ClusterSequence.hh"
48 #include "fastjet/JetDefinition.hh"
49 #include "fastjet/PseudoJet.hh"
50 
51 #include <cmath>
52 #include <cstdint>
53 #include <cstdlib>
54 #include <algorithm>
55 #include <limits>
56 #include <memory>
57 #include <sstream>
58 #include <string>
59 #include <utility>
60 #include <map>
61 #include <vector>
62 
63 struct EtaPhiTower {
64  int ieta, iphi;
65  float eta, phi;
66 };
67 
68 class HiPuRhoProducer : public edm::stream::EDProducer<> {
69 public:
70  explicit HiPuRhoProducer(const edm::ParameterSet&);
71 
72  using ClusterSequencePtr = std::shared_ptr<fastjet::ClusterSequence>;
73  using JetDefPtr = std::shared_ptr<fastjet::JetDefinition>;
74 
75  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
76  virtual void setupGeometryMap(edm::Event& iEvent, const edm::EventSetup& iSetup);
77  virtual void calculatePedestal(std::vector<fastjet::PseudoJet> const& coll);
78  virtual void subtractPedestal(std::vector<fastjet::PseudoJet>& coll);
79  virtual void calculateOrphanInput(std::vector<fastjet::PseudoJet>& orphanInput);
80  virtual void putRho(edm::Event& iEvent, const edm::EventSetup& iSetup);
81 
82 private:
83  void produce(edm::Event&, const edm::EventSetup&) override;
84 
85  // This checks if the tower is anomalous (if a calo tower).
86  virtual void inputTowers();
87 
88  bool postOrphan_;
89  bool setInitialValue_;
90 
91  const bool dropZeroTowers_;
92  const int medianWindowWidth_;
93  const double minimumTowersFraction_;
94  const double nSigmaPU_; // number of sigma for pileup
95  const double puPtMin_;
96  const double radiusPU_; // pileup radius
97  const double rParam_; // the R parameter to use
98  const double towSigmaCut_;
99  const edm::EDGetTokenT<CaloTowerCollection> caloTowerToken_;
100  const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometryToken_;
101  const int initialValue_ = -99;
102  constexpr static int nEtaTow_ = 82;
103 
104  std::array<int, nEtaTow_> vngeom_;
105  std::array<int, nEtaTow_> vntow_;
106  std::array<float, nEtaTow_> vmean0_;
107  std::array<float, nEtaTow_> vrms0_;
108  std::array<float, nEtaTow_> vrho0_;
109  std::array<float, nEtaTow_> vmean1_;
110  std::array<float, nEtaTow_> vrms1_;
111  std::array<float, nEtaTow_> vrho1_;
112 
113  std::vector<double> etaEdgeLow_;
114  std::vector<double> etaEdgeHi_;
115  std::vector<double> etaEdges_;
116 
117  std::vector<double> rho_;
118  std::vector<double> rhoExtra_;
119  std::vector<double> rhoM_;
120  std::vector<int> nTow_;
121 
122  std::vector<double> towExcludePt_;
123  std::vector<double> towExcludePhi_;
124  std::vector<double> towExcludeEta_;
125 
126  std::vector<const CaloTower*> inputs_;
127  ClusterSequencePtr fjClusterSeq_; // fastjet cluster sequence
128  JetDefPtr fjJetDefinition_; // fastjet jet definition
129 
130  std::vector<fastjet::PseudoJet> fjInputs_; // fastjet inputs
131  std::vector<fastjet::PseudoJet> fjJets_; // fastjet jets
132  std::vector<fastjet::PseudoJet> fjOriginalInputs_; // to back-up unsubtracted fastjet inputs
133 
134  CaloGeometry const* geo_ = nullptr; // geometry
135  std::vector<HcalDetId> allgeomid_; // all det ids in the geometry
136 
137  int ietamax_; // maximum eta in geometry
138  int ietamin_; // minimum eta in geometry
139  std::map<int, int> ntowersWithJets_; // number of towers with jets
140  std::map<int, int> geomtowers_; // map of geometry towers to det id
141  std::map<int, double> esigma_; // energy sigma
142  std::map<int, double> emean_; // energy mean
143  std::map<int, std::array<double, 4>> eTop4_; // energy mean
144 
145  typedef std::pair<double, double> EtaPhi;
146  std::vector<EtaPhiTower> towermap_;
147 };
148 
149 HiPuRhoProducer::HiPuRhoProducer(const edm::ParameterSet& iConfig)
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 }
171 
172 // ------------ method called to produce the data ------------
173 void HiPuRhoProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
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();
187  fjOriginalInputs_ = fjInputs_;
188  setInitialValue_ = true;
189  calculatePedestal(fjInputs_);
190  subtractPedestal(fjInputs_);
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 }
219 
220 void HiPuRhoProducer::inputTowers() {
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 }
234 
235 void HiPuRhoProducer::setupGeometryMap(edm::Event& iEvent, const edm::EventSetup& iSetup) {
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 }
274 
275 void HiPuRhoProducer::calculatePedestal(std::vector<fastjet::PseudoJet> const& coll) {
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 }
421 
422 void HiPuRhoProducer::subtractPedestal(std::vector<fastjet::PseudoJet>& coll) {
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 
437  math::XYZTLorentzVectorD towP4(
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 }
450 
451 void HiPuRhoProducer::calculateOrphanInput(std::vector<fastjet::PseudoJet>& orphanInput) {
452  LogDebug("PileUpSubtractor") << "The subtractor calculating orphan input...\n";
453 
454  fjInputs_ = fjOriginalInputs_;
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  (im.ieta - ntowersWithJets_[im.ieta]) > minimumTowersFraction_ * 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 }
501 
502 void HiPuRhoProducer::putRho(edm::Event& iEvent, const edm::EventSetup& iSetup) {
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 }
561 
562 // ------------ method fills 'descriptions' with the allowed parameters for the module ------------
563 void HiPuRhoProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
564  edm::ParameterSetDescription desc;
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 }
576 
577 //define this as a plug-in
578 DEFINE_FWK_MODULE(HiPuRhoProducer);
HiPuRhoProducer::fjInputs_
std::vector< fastjet::PseudoJet > fjInputs_
Definition: HiPuRhoProducer.cc:130
ConfigurationDescriptions.h
HiPuRhoProducer::towExcludePt_
std::vector< double > towExcludePt_
Definition: HiPuRhoProducer.cc:122
HiPuRhoProducer::towExcludeEta_
std::vector< double > towExcludeEta_
Definition: HiPuRhoProducer.cc:124
Handle.h
electrons_cff.bool
bool
Definition: electrons_cff.py:366
mps_fire.i
i
Definition: mps_fire.py:428
input
static const std::string input
Definition: EdmProvDump.cc:48
HiPuRhoProducer::emean_
std::map< int, double > emean_
Definition: HiPuRhoProducer.cc:142
MessageLogger.h
CaloTower.h
HiPuRhoProducer::nEtaTow_
constexpr static int nEtaTow_
Definition: HiPuRhoProducer.cc:102
HiPuRhoProducer
Definition: HiPuRhoProducer.cc:68
HiPuRhoProducer::fjJets_
std::vector< fastjet::PseudoJet > fjJets_
Definition: HiPuRhoProducer.cc:131
ESHandle.h
ESInputTag
edm::isNotFinite
constexpr bool isNotFinite(T x)
Definition: isFinite.h:9
HcalDetId::iphi
constexpr int iphi() const
get the cell iphi
Definition: HcalDetId.h:157
nt
int nt
Definition: AMPTWrapper.h:42
min
T min(T a, T b)
Definition: MathUtil.h:58
CaloGeometry::getPosition
GlobalPoint getPosition(const DetId &id) const
Get the position of a given detector id.
Definition: CaloGeometry.cc:50
edm::EDGetTokenT
Definition: EDGetToken.h:33
CaloGeometryRecord
Definition: CaloGeometryRecord.h:30
edm
HLT enums.
Definition: AlignableModifier.h:19
EtaPhiTower::ieta
int ieta
Definition: HiPuRhoProducer.cc:64
HiPuRhoProducer::vngeom_
std::array< int, nEtaTow_ > vngeom_
Definition: HiPuRhoProducer.cc:104
math::XYZTLorentzVectorD
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > XYZTLorentzVectorD
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:14
DetId::Hcal
Definition: DetId.h:28
EtaPhiTower::eta
float eta
Definition: HiPuRhoProducer.cc:65
HLT_FULL_cff.InputTag
InputTag
Definition: HLT_FULL_cff.py:89301
edm::ParameterSetDescription
Definition: ParameterSetDescription.h:52
HiPuRhoProducer::ntowersWithJets_
std::map< int, int > ntowersWithJets_
Definition: HiPuRhoProducer.cc:139
HiPuRhoProducer::inputs_
std::vector< const CaloTower * > inputs_
Definition: HiPuRhoProducer.cc:126
EDProducer.h
edm::SortedCollection< CaloTower >
Validation_hcalonly_cfi.sign
sign
Definition: Validation_hcalonly_cfi.py:32
reco::LeafCandidate::pt
double pt() const final
transverse momentum
Definition: LeafCandidate.h:146
HiPuRhoProducer::nSigmaPU_
const double nSigmaPU_
Definition: HiPuRhoProducer.cc:94
HiPuRhoProducer::towSigmaCut_
const double towSigmaCut_
Definition: HiPuRhoProducer.cc:98
HiPuRhoProducer::fjClusterSeq_
ClusterSequencePtr fjClusterSeq_
Definition: HiPuRhoProducer.cc:127
spr::find
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:19
HiPuRhoProducer::subtractPedestal
virtual void subtractPedestal(std::vector< fastjet::PseudoJet > &coll)
Definition: HiPuRhoProducer.cc:422
geometryDiff.epsilon
int epsilon
Definition: geometryDiff.py:26
HiPuRhoProducer::rhoM_
std::vector< double > rhoM_
Definition: HiPuRhoProducer.cc:119
HiPuRhoProducer::calculateOrphanInput
virtual void calculateOrphanInput(std::vector< fastjet::PseudoJet > &orphanInput)
Definition: HiPuRhoProducer.cc:451
HiPuRhoProducer::setInitialValue_
bool setInitialValue_
Definition: HiPuRhoProducer.cc:89
HiPuRhoProducer::JetDefPtr
std::shared_ptr< fastjet::JetDefinition > JetDefPtr
Definition: HiPuRhoProducer.cc:73
CandidateFwd.h
MakerMacros.h
CaloGeometry
Definition: CaloGeometry.h:21
HiPuRhoProducer::vrho1_
std::array< float, nEtaTow_ > vrho1_
Definition: HiPuRhoProducer.cc:111
HiPuRhoProducer::etaEdgeLow_
std::vector< double > etaEdgeLow_
Definition: HiPuRhoProducer.cc:113
submitPVValidationJobs.gt
list gt
Definition: submitPVValidationJobs.py:663
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
HiPuRhoProducer::fjJetDefinition_
JetDefPtr fjJetDefinition_
Definition: HiPuRhoProducer.cc:128
edm::ConfigurationDescriptions::add
void add(std::string const &label, ParameterSetDescription const &psetDescription)
Definition: ConfigurationDescriptions.cc:57
reco::LeafCandidate::py
double py() const final
y coordinate of momentum vector
Definition: LeafCandidate.h:142
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
hi::etaedge
constexpr std::array< double, 42 > etaedge
Definition: HITowerHelper.h:5
HiPuRhoProducer::initialValue_
const int initialValue_
Definition: HiPuRhoProducer.cc:101
CaloTower::ieta
int ieta() const
Definition: CaloTower.h:200
HiPuRhoProducer::fillDescriptions
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
Definition: HiPuRhoProducer.cc:563
HiPuRhoProducer::produce
void produce(edm::Event &, const edm::EventSetup &) override
Definition: HiPuRhoProducer.cc:173
HiPuRhoProducer::dropZeroTowers_
const bool dropZeroTowers_
Definition: HiPuRhoProducer.cc:91
HiPuRhoProducer::towExcludePhi_
std::vector< double > towExcludePhi_
Definition: HiPuRhoProducer.cc:123
ParameterSetDescription.h
EtaPhiTower
Definition: HiPuRhoProducer.cc:63
EDGetToken.h
HiPuRhoProducer::rho_
std::vector< double > rho_
Definition: HiPuRhoProducer.cc:117
CaloGeometryRecord.h
edm::ConfigurationDescriptions
Definition: ConfigurationDescriptions.h:28
HiPuRhoProducer::towermap_
std::vector< EtaPhiTower > towermap_
Definition: HiPuRhoProducer.cc:146
HiPuRhoProducer::putRho
virtual void putRho(edm::Event &iEvent, const edm::EventSetup &iSetup)
Definition: HiPuRhoProducer.cc:502
HiPuRhoProducer::minimumTowersFraction_
const double minimumTowersFraction_
Definition: HiPuRhoProducer.cc:93
HcalDetId::ieta
constexpr int ieta() const
get the cell ieta
Definition: HcalDetId.h:155
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:233
edm::ParameterSet
Definition: ParameterSet.h:47
EtaPhiTower::iphi
int iphi
Definition: HiPuRhoProducer.cc:64
Event.h
reco::LeafCandidate::eta
double eta() const final
momentum pseudorapidity
Definition: LeafCandidate.h:152
HiPuRhoProducer::HiPuRhoProducer
HiPuRhoProducer(const edm::ParameterSet &)
Definition: HiPuRhoProducer.cc:149
deltaR.h
HiPuRhoProducer::vrms1_
std::array< float, nEtaTow_ > vrms1_
Definition: HiPuRhoProducer.cc:110
HiPuRhoProducer::vrho0_
std::array< float, nEtaTow_ > vrho0_
Definition: HiPuRhoProducer.cc:108
HiPuRhoProducer::fjOriginalInputs_
std::vector< fastjet::PseudoJet > fjOriginalInputs_
Definition: HiPuRhoProducer.cc:132
jetUpdater_cfi.sort
sort
Definition: jetUpdater_cfi.py:29
HcalDetId.h
reco::deltaR2
constexpr auto deltaR2(const T1 &t1, const T2 &t2) -> decltype(t1.eta())
Definition: deltaR.h:16
StorageManager_cfg.e1
e1
Definition: StorageManager_cfg.py:16
PV3DBase::eta
T eta() const
Definition: PV3DBase.h:73
eventshapeDQM_cfi.order
order
Definition: eventshapeDQM_cfi.py:8
LorentzVector.h
Ptr.h
HcalDetId
Definition: HcalDetId.h:12
createfilelist.int
int
Definition: createfilelist.py:10
iEvent
int iEvent
Definition: GenABIO.cc:224
CaloTower
Definition: CaloTower.h:26
HiPuRhoProducer::eTop4_
std::map< int, std::array< double, 4 > > eTop4_
Definition: HiPuRhoProducer.cc:143
HiPuRhoProducer::postOrphan_
bool postOrphan_
Definition: HiPuRhoProducer.cc:88
M_PI
#define M_PI
Definition: BXVectorInputProducer.cc:49
CaloTower::et
double et(double vtxZ) const
Definition: CaloTower.h:150
HiPuRhoProducer::rhoExtra_
std::vector< double > rhoExtra_
Definition: HiPuRhoProducer.cc:118
edm::stream::EDProducer
Definition: EDProducer.h:36
CaloTowerDetId.h
edm::EventSetup
Definition: EventSetup.h:58
EtaPhiTower::phi
float phi
Definition: HiPuRhoProducer.cc:65
edm::ESGetToken< CaloGeometry, CaloGeometryRecord >
InputTag.h
HiPuRhoProducer::vntow_
std::array< int, nEtaTow_ > vntow_
Definition: HiPuRhoProducer.cc:105
edm::EventSetup::getData
bool getData(T &iHolder) const
Definition: EventSetup.h:127
HiPuRhoProducer::inputTowers
virtual void inputTowers()
Definition: HiPuRhoProducer.cc:220
HiPuRhoProducer::ietamax_
int ietamax_
Definition: HiPuRhoProducer.cc:137
CaloTowerCollection.h
HiPuRhoProducer::ClusterSequencePtr
std::shared_ptr< fastjet::ClusterSequence > ClusterSequencePtr
Definition: HiPuRhoProducer.cc:72
HiPuRhoProducer::EtaPhi
std::pair< double, double > EtaPhi
Definition: HiPuRhoProducer.cc:145
submitPVResolutionJobs.desc
string desc
Definition: submitPVResolutionJobs.py:251
HiPuRhoProducer::allgeomid_
std::vector< HcalDetId > allgeomid_
Definition: HiPuRhoProducer.cc:135
eostools.move
def move(src, dest)
Definition: eostools.py:511
CaloGeometry::getValidDetIds
std::vector< DetId > getValidDetIds() const
Get the list of all valid detector ids.
Definition: CaloGeometry.cc:75
HiPuRhoProducer::radiusPU_
const double radiusPU_
Definition: HiPuRhoProducer.cc:96
HiPuRhoProducer::caloTowerToken_
const edm::EDGetTokenT< CaloTowerCollection > caloTowerToken_
Definition: HiPuRhoProducer.cc:99
reco::LeafCandidate::phi
double phi() const final
momentum azimuthal angle
Definition: LeafCandidate.h:148
isFinite.h
HiPuRhoProducer::etaEdges_
std::vector< double > etaEdges_
Definition: HiPuRhoProducer.cc:115
DetId.h
HiPuRhoProducer::vmean1_
std::array< float, nEtaTow_ > vmean1_
Definition: HiPuRhoProducer.cc:109
CaloGeometry.h
HiPuRhoProducer::puPtMin_
const double puPtMin_
Definition: HiPuRhoProducer.cc:95
HiPuRhoProducer::caloGeometryToken_
const edm::ESGetToken< CaloGeometry, CaloGeometryRecord > caloGeometryToken_
Definition: HiPuRhoProducer.cc:100
HiPuRhoProducer::vmean0_
std::array< float, nEtaTow_ > vmean0_
Definition: HiPuRhoProducer.cc:106
EventSetup.h
HiPuRhoProducer::nTow_
std::vector< int > nTow_
Definition: HiPuRhoProducer.cc:120
Exception.h
CaloTower::iphi
int iphi() const
Definition: CaloTower.h:202
reco::LeafCandidate::energy
double energy() const final
energy
Definition: LeafCandidate.h:125
AlignmentPI::index
index
Definition: AlignmentPayloadInspectorHelper.h:46
HiPuRhoProducer::setupGeometryMap
virtual void setupGeometryMap(edm::Event &iEvent, const edm::EventSetup &iSetup)
Definition: HiPuRhoProducer.cc:235
Candidate.h
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
HiPuRhoProducer::vrms0_
std::array< float, nEtaTow_ > vrms0_
Definition: HiPuRhoProducer.cc:107
View.h
HiPuRhoProducer::medianWindowWidth_
const int medianWindowWidth_
Definition: HiPuRhoProducer.cc:92
ParameterSet.h
ParameterSetDescriptionFiller.h
HITowerHelper.h
HiPuRhoProducer::geo_
CaloGeometry const * geo_
Definition: HiPuRhoProducer.cc:134
DeDxTools::esConsumes
ESGetTokenH3DDVariant esConsumes(std::string const &Reccord, edm::ConsumesCollector &)
Definition: DeDxTools.cc:283
HiPuRhoProducer::geomtowers_
std::map< int, int > geomtowers_
Definition: HiPuRhoProducer.cc:140
edm::Event
Definition: Event.h:73
HiPuRhoProducer::etaEdgeHi_
std::vector< double > etaEdgeHi_
Definition: HiPuRhoProducer.cc:114
reco::LeafCandidate::px
double px() const final
x coordinate of momentum vector
Definition: LeafCandidate.h:140
HiPuRhoProducer::esigma_
std::map< int, double > esigma_
Definition: HiPuRhoProducer.cc:141
GlobalPoint.h
reco::LeafCandidate::pz
double pz() const final
z coordinate of momentum vector
Definition: LeafCandidate.h:144
edm::InputTag
Definition: InputTag.h:15
PV3DBase::phi
Geom::Phi< T > phi() const
Definition: PV3DBase.h:66
HiPuRhoProducer::rParam_
const double rParam_
Definition: HiPuRhoProducer.cc:97
findQualityFiles.size
size
Write out results.
Definition: findQualityFiles.py:443
HiPuRhoProducer::ietamin_
int ietamin_
Definition: HiPuRhoProducer.cc:138
HiPuRhoProducer::calculatePedestal
virtual void calculatePedestal(std::vector< fastjet::PseudoJet > const &coll)
Definition: HiPuRhoProducer.cc:275
Generated for CMSSW Reference Manual by doxygen 1.8.16