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/Framework/src/WorkerMaker.h"
30 #include "FWCore/MessageLogger/interface/MessageLogger.h"
31 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
32 #include "FWCore/ParameterSet/interface/ParameterSet.h"
33 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
34 #include "FWCore/ParameterSet/interface/ParameterSetDescriptionFiller.h"
35 #include "FWCore/Utilities/interface/Exception.h"
36 #include "FWCore/Utilities/interface/isFinite.h"
37 #include "Geometry/Records/interface/CaloGeometryRecord.h"
38 #include "DataFormats/Candidate/interface/CandidateFwd.h"
39 #include "DataFormats/Common/interface/Ptr.h"
40 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
41 #include "FWCore/Framework/interface/stream/EDProducer.h"
42 #include "FWCore/Framework/interface/EventSetup.h"
43 #include "FWCore/Utilities/interface/EDGetToken.h"
44 #include "FWCore/Utilities/interface/InputTag.h"
45 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
46 #include "RecoHI/HiJetAlgos/plugins/HITowerHelper.h"
47 
48 #include "fastjet/ClusterSequence.hh"
49 #include "fastjet/JetDefinition.hh"
50 #include "fastjet/PseudoJet.hh"
51 
52 #include <cmath>
53 #include <cstdint>
54 #include <cstdlib>
55 #include <algorithm>
56 #include <limits>
57 #include <memory>
58 #include <sstream>
59 #include <string>
60 #include <utility>
61 #include <map>
62 #include <vector>
63 
64 struct EtaPhiTower {
65  int ieta, iphi;
66  float eta, phi;
67 };
68 
69 class HiPuRhoProducer : public edm::stream::EDProducer<> {
70 public:
71  explicit HiPuRhoProducer(const edm::ParameterSet&);
72 
73  using ClusterSequencePtr = std::shared_ptr<fastjet::ClusterSequence>;
74  using JetDefPtr = std::shared_ptr<fastjet::JetDefinition>;
75 
76  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
77  virtual void setupGeometryMap(edm::Event& iEvent, const edm::EventSetup& iSetup);
78  virtual void calculatePedestal(std::vector<fastjet::PseudoJet> const& coll);
79  virtual void subtractPedestal(std::vector<fastjet::PseudoJet>& coll);
80  virtual void calculateOrphanInput(std::vector<fastjet::PseudoJet>& orphanInput);
81  virtual void putRho(edm::Event& iEvent, const edm::EventSetup& iSetup);
82 
83 private:
84  void produce(edm::Event&, const edm::EventSetup&) override;
85 
86  // This checks if the tower is anomalous (if a calo tower).
87  virtual void inputTowers();
88 
89  bool postOrphan_;
90  bool setInitialValue_;
91 
92  const bool dropZeroTowers_;
93  const int medianWindowWidth_;
94  const double minimumTowersFraction_;
95  const double nSigmaPU_; // number of sigma for pileup
96  const double puPtMin_;
97  const double radiusPU_; // pileup radius
98  const double rParam_; // the R parameter to use
99  const double towSigmaCut_;
100  const edm::EDGetTokenT<CaloTowerCollection> caloTowerToken_;
101  const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometryToken_;
102  const int initialValue_ = -99;
103  constexpr static int nEtaTow_ = 82;
104 
105  std::array<int, nEtaTow_> vngeom_;
106  std::array<int, nEtaTow_> vntow_;
107  std::array<float, nEtaTow_> vmean0_;
108  std::array<float, nEtaTow_> vrms0_;
109  std::array<float, nEtaTow_> vrho0_;
110  std::array<float, nEtaTow_> vmean1_;
111  std::array<float, nEtaTow_> vrms1_;
112  std::array<float, nEtaTow_> vrho1_;
113 
114  std::vector<double> etaEdgeLow_;
115  std::vector<double> etaEdgeHi_;
116  std::vector<double> etaEdges_;
117 
118  std::vector<double> rho_;
119  std::vector<double> rhoExtra_;
120  std::vector<double> rhoM_;
121  std::vector<int> nTow_;
122 
123  std::vector<double> towExcludePt_;
124  std::vector<double> towExcludePhi_;
125  std::vector<double> towExcludeEta_;
126 
127  std::vector<const CaloTower*> inputs_;
128  ClusterSequencePtr fjClusterSeq_; // fastjet cluster sequence
129  JetDefPtr fjJetDefinition_; // fastjet jet definition
130 
131  std::vector<fastjet::PseudoJet> fjInputs_; // fastjet inputs
132  std::vector<fastjet::PseudoJet> fjJets_; // fastjet jets
133  std::vector<fastjet::PseudoJet> fjOriginalInputs_; // to back-up unsubtracted fastjet inputs
134 
135  CaloGeometry const* geo_ = nullptr; // geometry
136  std::vector<HcalDetId> allgeomid_; // all det ids in the geometry
137 
138  int ietamax_; // maximum eta in geometry
139  int ietamin_; // minimum eta in geometry
140  std::map<int, int> ntowersWithJets_; // number of towers with jets
141  std::map<int, int> geomtowers_; // map of geometry towers to det id
142  std::map<int, double> esigma_; // energy sigma
143  std::map<int, double> emean_; // energy mean
144  std::map<int, std::array<double, 4>> eTop4_; // energy mean
145 
146  typedef std::pair<double, double> EtaPhi;
147  std::vector<EtaPhiTower> towermap_;
148 };
149 
150 HiPuRhoProducer::HiPuRhoProducer(const edm::ParameterSet& iConfig)
151  : dropZeroTowers_(iConfig.getParameter<bool>("dropZeroTowers")),
152  medianWindowWidth_(iConfig.getParameter<int>("medianWindowWidth")),
153  minimumTowersFraction_(iConfig.getParameter<double>("minimumTowersFraction")),
154  nSigmaPU_(iConfig.getParameter<double>("nSigmaPU")),
155  puPtMin_(iConfig.getParameter<double>("puPtMin")),
156  radiusPU_(iConfig.getParameter<double>("radiusPU")),
157  rParam_(iConfig.getParameter<double>("rParam")),
158  towSigmaCut_(iConfig.getParameter<double>("towSigmaCut")),
159  caloTowerToken_(consumes<CaloTowerCollection>(iConfig.getParameter<edm::InputTag>("src"))),
160  caloGeometryToken_(esConsumes<CaloGeometry, CaloGeometryRecord>(edm::ESInputTag{})) {
161  //register your products
162  produces<std::vector<double>>("mapEtaEdges");
163  produces<std::vector<double>>("mapToRho");
164  produces<std::vector<double>>("mapToRhoMedian");
165  produces<std::vector<double>>("mapToRhoExtra");
166  produces<std::vector<double>>("mapToRhoM");
167  produces<std::vector<int>>("mapToNTow");
168  produces<std::vector<double>>("mapToTowExcludePt");
169  produces<std::vector<double>>("mapToTowExcludePhi");
170  produces<std::vector<double>>("mapToTowExcludeEta");
171 }
172 
173 // ------------ method called to produce the data ------------
174 void HiPuRhoProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
175  setupGeometryMap(iEvent, iSetup);
176 
177  for (int i = ietamin_; i < ietamax_ + 1; i++) {
178  ntowersWithJets_[i] = 0;
179  }
180 
181  auto const& inputView = iEvent.get(caloTowerToken_);
182  inputs_.reserve(inputView.size());
183  for (auto const& input : inputView)
184  inputs_.push_back(&input);
185 
186  fjInputs_.reserve(inputs_.size());
187  inputTowers();
188  fjOriginalInputs_ = fjInputs_;
189  setInitialValue_ = true;
190  calculatePedestal(fjInputs_);
191  subtractPedestal(fjInputs_);
192 
193  fjJetDefinition_ = std::make_shared<fastjet::JetDefinition>(fastjet::antikt_algorithm, rParam_);
194  fjClusterSeq_ = std::make_shared<fastjet::ClusterSequence>(fjInputs_, *fjJetDefinition_);
195  fjJets_ = fastjet::sorted_by_pt(fjClusterSeq_->inclusive_jets(puPtMin_));
196 
197  etaEdgeLow_.clear();
198  etaEdgeHi_.clear();
199  etaEdges_.clear();
200 
201  rho_.clear();
202  rhoExtra_.clear();
203  rhoM_.clear();
204  nTow_.clear();
205 
206  towExcludePt_.clear();
207  towExcludePhi_.clear();
208  towExcludeEta_.clear();
209 
210  setInitialValue_ = false;
211  std::vector<fastjet::PseudoJet> orphanInput;
212  calculateOrphanInput(orphanInput);
213  calculatePedestal(orphanInput);
214  putRho(iEvent, iSetup);
215 
216  inputs_.clear();
217  fjInputs_.clear();
218  fjJets_.clear();
219 }
220 
221 void HiPuRhoProducer::inputTowers() {
222  int index = -1;
223  for (auto const& input : inputs_) {
224  index++;
225 
226  if (edm::isNotFinite(input->pt()))
227  continue;
228  if (input->pt() < 100 * std::numeric_limits<double>::epsilon())
229  continue;
230 
231  fjInputs_.push_back(fastjet::PseudoJet(input->px(), input->py(), input->pz(), input->energy()));
232  fjInputs_.back().set_user_index(index);
233  }
234 }
235 
236 void HiPuRhoProducer::setupGeometryMap(edm::Event& iEvent, const edm::EventSetup& iSetup) {
237  LogDebug("PileUpSubtractor") << "The subtractor setting up geometry...\n";
238  const auto& pG = iSetup.getData(caloGeometryToken_);
239  geo_ = &pG;
240  std::vector<DetId> alldid = geo_->getValidDetIds();
241  int ietaold = -10000;
242  ietamax_ = -10000;
243  ietamin_ = 10000;
244  towermap_.clear();
245 
246  for (auto const& did : alldid) {
247  if (did.det() == DetId::Hcal) {
248  HcalDetId hid = HcalDetId(did);
249  allgeomid_.push_back(did);
250  towermap_.push_back({hid.ieta(), hid.iphi(), geo_->getPosition(did).eta(), geo_->getPosition(did).phi()});
251  if (hid.ieta() != ietaold) {
252  ietaold = hid.ieta();
253  geomtowers_[hid.ieta()] = 1;
254  if (hid.ieta() > ietamax_)
255  ietamax_ = hid.ieta();
256  if (hid.ieta() < ietamin_)
257  ietamin_ = hid.ieta();
258  } else {
259  geomtowers_[hid.ieta()]++;
260  }
261  }
262  }
263 
264  for (auto const& gt : geomtowers_) {
265  int it = gt.first;
266  int vi = it - 1;
267 
268  if (it < 0)
269  vi = nEtaTow_ + it;
270 
271  if (vi < nEtaTow_)
272  vngeom_[vi] = gt.second;
273  }
274 }
275 
276 void HiPuRhoProducer::calculatePedestal(std::vector<fastjet::PseudoJet> const& coll) {
277  LogDebug("PileUpSubtractor") << "The subtractor calculating pedestals...\n";
278 
279  std::map<int, double> emean2;
280  std::map<int, int> ntowers;
281 
282  int ietaold = -10000;
283  // Initial values for emean_, emean2, esigma_, ntowers
284 
285  for (int vi = 0; vi < nEtaTow_; ++vi) {
286  vntow_[vi] = initialValue_;
287  vmean1_[vi] = initialValue_;
288  vrms1_[vi] = initialValue_;
289  vrho1_[vi] = initialValue_;
290 
291  if (setInitialValue_) {
292  vmean0_[vi] = initialValue_;
293  vrms0_[vi] = initialValue_;
294  vrho0_[vi] = initialValue_;
295  }
296  }
297 
298  for (int i = ietamin_; i < ietamax_ + 1; i++) {
299  emean_[i] = 0.;
300  emean2[i] = 0.;
301  esigma_[i] = 0.;
302  ntowers[i] = 0;
303 
304  eTop4_[i] = {{0., 0., 0., 0.}};
305  }
306 
307  for (auto const& input_object : coll) {
308  const CaloTower* originalTower = inputs_[input_object.user_index()];
309  double original_Et = originalTower->et();
310  int ieta0 = originalTower->ieta();
311 
312  if (original_Et > eTop4_[ieta0][0]) {
313  eTop4_[ieta0][3] = eTop4_[ieta0][2];
314  eTop4_[ieta0][2] = eTop4_[ieta0][1];
315  eTop4_[ieta0][1] = eTop4_[ieta0][0];
316  eTop4_[ieta0][0] = original_Et;
317  } else if (original_Et > eTop4_[ieta0][1]) {
318  eTop4_[ieta0][3] = eTop4_[ieta0][2];
319  eTop4_[ieta0][2] = eTop4_[ieta0][1];
320  eTop4_[ieta0][1] = original_Et;
321  } else if (original_Et > eTop4_[ieta0][2]) {
322  eTop4_[ieta0][3] = eTop4_[ieta0][2];
323  eTop4_[ieta0][2] = original_Et;
324  } else if (original_Et > eTop4_[ieta0][3]) {
325  eTop4_[ieta0][3] = original_Et;
326  }
327 
328  emean_[ieta0] = emean_[ieta0] + original_Et;
329  emean2[ieta0] = emean2[ieta0] + original_Et * original_Et;
330  if (ieta0 - ietaold != 0) {
331  ntowers[ieta0] = 1;
332  ietaold = ieta0;
333  } else {
334  ntowers[ieta0]++;
335  }
336  }
337 
338  for (auto const& gt : geomtowers_) {
339  int it = gt.first;
340 
341  int vi = it - 1;
342 
343  if (it < 0)
344  vi = nEtaTow_ + it;
345 
346  double e1 = emean_[it];
347  double e2 = emean2[it];
348  int nt = gt.second - ntowersWithJets_[it];
349 
350  if (vi < nEtaTow_) {
351  vntow_[vi] = nt;
352  }
353 
354  LogDebug("PileUpSubtractor") << " ieta: " << it << " number of towers: " << nt << " e1: " << e1 << " e2: " << e2
355  << "\n";
356 
357  if (nt > 0) {
358  if (postOrphan_) {
359  if (nt > (int)minimumTowersFraction_ * (gt.second)) {
360  emean_[it] = e1 / (double)nt;
361  double eee = e2 / (double)nt - e1 * e1 / (double)(nt * nt);
362  if (eee < 0.)
363  eee = 0.;
364  esigma_[it] = nSigmaPU_ * sqrt(eee);
365 
366  uint32_t numToCheck = std::min(int(eTop4_[it].size()), nt - (int)minimumTowersFraction_ * (gt.second));
367 
368  for (unsigned int lI = 0; lI < numToCheck; ++lI) {
369  if (eTop4_[it][lI] >= emean_[it] + towSigmaCut_ * esigma_[it] && towSigmaCut_ > 0) {
370  e1 -= eTop4_[it][lI];
371  nt -= 1;
372  } else
373  break;
374  }
375 
376  if (e1 < .000000001)
377  e1 = 0;
378  }
379  }
380 
381  emean_[it] = e1 / (double)nt;
382  double eee = e2 / nt - e1 * e1 / (nt * nt);
383  if (eee < 0.)
384  eee = 0.;
385  esigma_[it] = nSigmaPU_ * sqrt(eee);
386 
387  double etaWidth = std::abs(hi::etaedge[abs(it)] - hi::etaedge[abs(it) - 1]);
388 
389  int sign = (it < 0) ? -1 : 1;
390  auto absIt = std::abs(it);
391  if (it < 0) {
392  etaEdgeLow_.push_back(sign * hi::etaedge[absIt]);
393  etaEdgeHi_.push_back(sign * hi::etaedge[absIt - 1]);
394  } else {
395  etaEdgeHi_.push_back(sign * hi::etaedge[absIt]);
396  etaEdgeLow_.push_back(sign * hi::etaedge[absIt - 1]);
397  }
398 
399  if (vi < nEtaTow_) {
400  vmean1_[vi] = emean_[it];
401  vrho1_[vi] = emean_[it] / (etaWidth * (2. * M_PI / (double)vngeom_[vi]));
402  rho_.push_back(vrho1_[vi]);
403  rhoM_.push_back(0);
404  vrms1_[vi] = esigma_[it];
405  if (vngeom_[vi] == vntow_[vi]) {
406  vmean0_[vi] = emean_[it];
407  vrho0_[vi] = emean_[it] / (etaWidth * (2. * M_PI / (double)vngeom_[vi]));
408  vrms0_[vi] = esigma_[it];
409  }
410  rhoExtra_.push_back(vrho0_[vi]);
411  nTow_.push_back(vntow_[vi]);
412  }
413  } else {
414  emean_[it] = 0.;
415  esigma_[it] = 0.;
416  }
417  LogDebug("PileUpSubtractor") << " ieta: " << it << " Pedestals: " << emean_[it] << " " << esigma_[it] << "\n";
418  }
419 
420  postOrphan_ = false;
421 }
422 
423 void HiPuRhoProducer::subtractPedestal(std::vector<fastjet::PseudoJet>& coll) {
424  LogDebug("PileUpSubtractor") << "The subtractor subtracting pedestals...\n";
425 
426  std::vector<fastjet::PseudoJet> newcoll;
427  for (auto& input_object : coll) {
428  int index = input_object.user_index();
429  CaloTower const* itow = inputs_[index];
430  int it = itow->ieta();
431 
432  double original_Et = itow->et();
433  double etnew = original_Et - emean_.at(it) - esigma_.at(it);
434  float mScale = etnew / input_object.Et();
435  if (etnew < 0.)
436  mScale = 0.;
437 
438  math::XYZTLorentzVectorD towP4(
439  input_object.px() * mScale, input_object.py() * mScale, input_object.pz() * mScale, input_object.e() * mScale);
440 
441  input_object.reset_momentum(towP4.px(), towP4.py(), towP4.pz(), towP4.energy());
442  input_object.set_user_index(index);
443 
444  if (etnew > 0. && dropZeroTowers_)
445  newcoll.push_back(input_object);
446  }
447 
448  if (dropZeroTowers_)
449  coll = newcoll;
450 }
451 
452 void HiPuRhoProducer::calculateOrphanInput(std::vector<fastjet::PseudoJet>& orphanInput) {
453  LogDebug("PileUpSubtractor") << "The subtractor calculating orphan input...\n";
454 
455  fjInputs_ = fjOriginalInputs_;
456 
457  std::vector<int> jettowers; // vector of towers indexed by "user_index"
458  std::map<std::pair<int, int>, int> excludedTowers; // map of excluded ieta, iphi values
459  for (auto const& pseudojetTMP : fjJets_) {
460  EtaPhi jet_etaphi(pseudojetTMP.eta(), pseudojetTMP.phi());
461 
462  if (pseudojetTMP.perp() < puPtMin_)
463  continue;
464 
465  for (auto const& im : towermap_) {
466  double dr2 = reco::deltaR2(im.eta, im.phi, jet_etaphi.first, jet_etaphi.second);
467  if (dr2 < radiusPU_ * radiusPU_ && !excludedTowers[std::pair(im.ieta, im.iphi)] &&
468  (im.ieta - ntowersWithJets_[im.ieta]) > minimumTowersFraction_ * im.ieta) {
469  ntowersWithJets_[im.ieta]++;
470  excludedTowers[std::pair(im.ieta, im.iphi)] = 1;
471  }
472  }
473 
474  for (auto const& input : fjInputs_) {
475  int index = input.user_index();
476  const CaloTower* originalTower = inputs_[index];
477  int ie = originalTower->ieta();
478  int ip = originalTower->iphi();
479 
480  if (excludedTowers[std::pair<int, int>(ie, ip)]) {
481  jettowers.push_back(index);
482  }
483  }
484  }
485  // Create a new collections from the towers not included in jets
486  for (auto const& input : fjInputs_) {
487  int index = input.user_index();
488  const CaloTower* originalTower = inputs_[index];
489  auto itjet = std::find(jettowers.begin(), jettowers.end(), index);
490  if (itjet == jettowers.end()) {
491  orphanInput.emplace_back(originalTower->px(), originalTower->py(), originalTower->pz(), originalTower->energy());
492  orphanInput.back().set_user_index(index);
493  } else {
494  towExcludePt_.push_back(originalTower->pt());
495  towExcludePhi_.push_back(originalTower->phi());
496  towExcludeEta_.push_back(originalTower->eta());
497  }
498  }
499 
500  postOrphan_ = true;
501 }
502 
503 void HiPuRhoProducer::putRho(edm::Event& iEvent, const edm::EventSetup& iSetup) {
504  std::size_t size = etaEdgeLow_.size();
505 
506  std::vector<std::pair<std::size_t, double>> order;
507  for (std::size_t i = 0; i < size; ++i) {
508  order.emplace_back(i, etaEdgeLow_[i]);
509  }
510 
511  std::sort(
512  order.begin(), order.end(), [](auto const& pair0, auto const& pair1) { return pair0.second < pair1.second; });
513 
514  std::vector<double> sortedEtaEdgeLow(size);
515  std::vector<double> sortedEtaEdgeHigh(size);
516 
517  auto mapToRhoOut = std::make_unique<std::vector<double>>(size);
518  auto mapToRhoExtraOut = std::make_unique<std::vector<double>>(size);
519  auto mapToRhoMOut = std::make_unique<std::vector<double>>(size);
520  auto mapToNTowOut = std::make_unique<std::vector<int>>(size);
521 
522  for (std::size_t i = 0; i < size; ++i) {
523  auto const& index = order[i].first;
524 
525  sortedEtaEdgeLow[i] = etaEdgeLow_[index];
526  sortedEtaEdgeHigh[i] = etaEdgeHi_[index];
527 
528  (*mapToRhoOut)[i] = rho_[index];
529  (*mapToRhoExtraOut)[i] = rhoExtra_[index];
530  (*mapToRhoMOut)[i] = rhoM_[index];
531  (*mapToNTowOut)[i] = nTow_[index];
532  }
533 
534  auto mapToRhoMedianOut = std::make_unique<std::vector<double>>(size);
535 
536  for (uint32_t index = medianWindowWidth_; index < size - medianWindowWidth_; ++index) {
537  auto centre = mapToRhoOut->begin() + index;
538  std::vector<float> rhoRange(centre - medianWindowWidth_, centre + medianWindowWidth_);
539  std::nth_element(rhoRange.begin(), rhoRange.begin() + medianWindowWidth_, rhoRange.end());
540  (*mapToRhoMedianOut)[index] = rhoRange[medianWindowWidth_];
541  }
542 
543  auto mapEtaRangesOut = std::make_unique<std::vector<double>>();
544 
545  mapEtaRangesOut->assign(sortedEtaEdgeLow.begin(), sortedEtaEdgeLow.end());
546  mapEtaRangesOut->push_back(sortedEtaEdgeHigh.back());
547 
548  auto mapToTowExcludePtOut = std::make_unique<std::vector<double>>(std::move(towExcludePt_));
549  auto mapToTowExcludePhiOut = std::make_unique<std::vector<double>>(std::move(towExcludePhi_));
550  auto mapToTowExcludeEtaOut = std::make_unique<std::vector<double>>(std::move(towExcludeEta_));
551 
552  iEvent.put(std::move(mapEtaRangesOut), "mapEtaEdges");
553  iEvent.put(std::move(mapToRhoOut), "mapToRho");
554  iEvent.put(std::move(mapToRhoMedianOut), "mapToRhoMedian");
555  iEvent.put(std::move(mapToRhoExtraOut), "mapToRhoExtra");
556  iEvent.put(std::move(mapToRhoMOut), "mapToRhoM");
557  iEvent.put(std::move(mapToNTowOut), "mapToNTow");
558  iEvent.put(std::move(mapToTowExcludePtOut), "mapToTowExcludePt");
559  iEvent.put(std::move(mapToTowExcludePhiOut), "mapToTowExcludePhi");
560  iEvent.put(std::move(mapToTowExcludeEtaOut), "mapToTowExcludeEta");
561 }
562 
563 // ------------ method fills 'descriptions' with the allowed parameters for the module ------------
564 void HiPuRhoProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
565  edm::ParameterSetDescription desc;
566  desc.add<edm::InputTag>("src", edm::InputTag("PFTowers"));
567  desc.add<int>("medianWindowWidth", 2);
568  desc.add<double>("towSigmaCut", 5.0);
569  desc.add<double>("puPtMin", 15.0);
570  desc.add<double>("rParam", 0.3);
571  desc.add<double>("nSigmaPU", 1.0);
572  desc.add<double>("radiusPU", 0.5);
573  desc.add<double>("minimumTowersFraction", 0.7);
574  desc.add<bool>("dropZeroTowers", true);
575  descriptions.add("hiPuRhoProducer", desc);
576 }
577 
578 //define this as a plug-in
579 DEFINE_FWK_MODULE(HiPuRhoProducer);
HiPuRhoProducer::fjInputs_
std::vector< fastjet::PseudoJet > fjInputs_
Definition: HiPuRhoProducer.cc:131
ConfigurationDescriptions.h
HiPuRhoProducer::towExcludePt_
std::vector< double > towExcludePt_
Definition: HiPuRhoProducer.cc:123
HiPuRhoProducer::towExcludeEta_
std::vector< double > towExcludeEta_
Definition: HiPuRhoProducer.cc:125
Handle.h
electrons_cff.bool
bool
Definition: electrons_cff.py:393
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:143
MessageLogger.h
CaloTower.h
HiPuRhoProducer::nEtaTow_
constexpr static int nEtaTow_
Definition: HiPuRhoProducer.cc:103
HiPuRhoProducer
Definition: HiPuRhoProducer.cc:69
HiPuRhoProducer::fjJets_
std::vector< fastjet::PseudoJet > fjJets_
Definition: HiPuRhoProducer.cc:132
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:65
HiPuRhoProducer::vngeom_
std::array< int, nEtaTow_ > vngeom_
Definition: HiPuRhoProducer.cc:105
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:66
HLT_FULL_cff.InputTag
InputTag
Definition: HLT_FULL_cff.py:89287
edm::ParameterSetDescription
Definition: ParameterSetDescription.h:52
HiPuRhoProducer::ntowersWithJets_
std::map< int, int > ntowersWithJets_
Definition: HiPuRhoProducer.cc:140
HiPuRhoProducer::inputs_
std::vector< const CaloTower * > inputs_
Definition: HiPuRhoProducer.cc:127
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:95
HiPuRhoProducer::towSigmaCut_
const double towSigmaCut_
Definition: HiPuRhoProducer.cc:99
HiPuRhoProducer::fjClusterSeq_
ClusterSequencePtr fjClusterSeq_
Definition: HiPuRhoProducer.cc:128
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:423
WorkerMaker.h
geometryDiff.epsilon
int epsilon
Definition: geometryDiff.py:26
HiPuRhoProducer::rhoM_
std::vector< double > rhoM_
Definition: HiPuRhoProducer.cc:120
HiPuRhoProducer::calculateOrphanInput
virtual void calculateOrphanInput(std::vector< fastjet::PseudoJet > &orphanInput)
Definition: HiPuRhoProducer.cc:452
HiPuRhoProducer::setInitialValue_
bool setInitialValue_
Definition: HiPuRhoProducer.cc:90
HiPuRhoProducer::JetDefPtr
std::shared_ptr< fastjet::JetDefinition > JetDefPtr
Definition: HiPuRhoProducer.cc:74
CandidateFwd.h
MakerMacros.h
CaloGeometry
Definition: CaloGeometry.h:21
HiPuRhoProducer::vrho1_
std::array< float, nEtaTow_ > vrho1_
Definition: HiPuRhoProducer.cc:112
HiPuRhoProducer::etaEdgeLow_
std::vector< double > etaEdgeLow_
Definition: HiPuRhoProducer.cc:114
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:129
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:102
CaloTower::ieta
int ieta() const
Definition: CaloTower.h:200
HiPuRhoProducer::fillDescriptions
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
Definition: HiPuRhoProducer.cc:564
HiPuRhoProducer::produce
void produce(edm::Event &, const edm::EventSetup &) override
Definition: HiPuRhoProducer.cc:174
HiPuRhoProducer::dropZeroTowers_
const bool dropZeroTowers_
Definition: HiPuRhoProducer.cc:92
HiPuRhoProducer::towExcludePhi_
std::vector< double > towExcludePhi_
Definition: HiPuRhoProducer.cc:124
ParameterSetDescription.h
EtaPhiTower
Definition: HiPuRhoProducer.cc:64
EDGetToken.h
HiPuRhoProducer::rho_
std::vector< double > rho_
Definition: HiPuRhoProducer.cc:118
CaloGeometryRecord.h
edm::ConfigurationDescriptions
Definition: ConfigurationDescriptions.h:28
HiPuRhoProducer::towermap_
std::vector< EtaPhiTower > towermap_
Definition: HiPuRhoProducer.cc:147
HiPuRhoProducer::putRho
virtual void putRho(edm::Event &iEvent, const edm::EventSetup &iSetup)
Definition: HiPuRhoProducer.cc:503
HiPuRhoProducer::minimumTowersFraction_
const double minimumTowersFraction_
Definition: HiPuRhoProducer.cc:94
HcalDetId::ieta
constexpr int ieta() const
get the cell ieta
Definition: HcalDetId.h:155
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:223
edm::ParameterSet
Definition: ParameterSet.h:47
EtaPhiTower::iphi
int iphi
Definition: HiPuRhoProducer.cc:65
Event.h
reco::LeafCandidate::eta
double eta() const final
momentum pseudorapidity
Definition: LeafCandidate.h:152
HiPuRhoProducer::HiPuRhoProducer
HiPuRhoProducer(const edm::ParameterSet &)
Definition: HiPuRhoProducer.cc:150
deltaR.h
HiPuRhoProducer::vrms1_
std::array< float, nEtaTow_ > vrms1_
Definition: HiPuRhoProducer.cc:111
HiPuRhoProducer::vrho0_
std::array< float, nEtaTow_ > vrho0_
Definition: HiPuRhoProducer.cc:109
HiPuRhoProducer::fjOriginalInputs_
std::vector< fastjet::PseudoJet > fjOriginalInputs_
Definition: HiPuRhoProducer.cc:133
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:144
HiPuRhoProducer::postOrphan_
bool postOrphan_
Definition: HiPuRhoProducer.cc:89
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:119
edm::stream::EDProducer
Definition: EDProducer.h:38
CaloTowerDetId.h
edm::EventSetup
Definition: EventSetup.h:57
EtaPhiTower::phi
float phi
Definition: HiPuRhoProducer.cc:66
edm::ESGetToken< CaloGeometry, CaloGeometryRecord >
InputTag.h
HiPuRhoProducer::vntow_
std::array< int, nEtaTow_ > vntow_
Definition: HiPuRhoProducer.cc:106
edm::EventSetup::getData
bool getData(T &iHolder) const
Definition: EventSetup.h:120
HiPuRhoProducer::inputTowers
virtual void inputTowers()
Definition: HiPuRhoProducer.cc:221
HiPuRhoProducer::ietamax_
int ietamax_
Definition: HiPuRhoProducer.cc:138
CaloTowerCollection.h
HiPuRhoProducer::ClusterSequencePtr
std::shared_ptr< fastjet::ClusterSequence > ClusterSequencePtr
Definition: HiPuRhoProducer.cc:73
HiPuRhoProducer::EtaPhi
std::pair< double, double > EtaPhi
Definition: HiPuRhoProducer.cc:146
submitPVResolutionJobs.desc
string desc
Definition: submitPVResolutionJobs.py:251
HiPuRhoProducer::allgeomid_
std::vector< HcalDetId > allgeomid_
Definition: HiPuRhoProducer.cc:136
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:97
HiPuRhoProducer::caloTowerToken_
const edm::EDGetTokenT< CaloTowerCollection > caloTowerToken_
Definition: HiPuRhoProducer.cc:100
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:116
DetId.h
HiPuRhoProducer::vmean1_
std::array< float, nEtaTow_ > vmean1_
Definition: HiPuRhoProducer.cc:110
CaloGeometry.h
HiPuRhoProducer::puPtMin_
const double puPtMin_
Definition: HiPuRhoProducer.cc:96
HiPuRhoProducer::caloGeometryToken_
const edm::ESGetToken< CaloGeometry, CaloGeometryRecord > caloGeometryToken_
Definition: HiPuRhoProducer.cc:101
HiPuRhoProducer::vmean0_
std::array< float, nEtaTow_ > vmean0_
Definition: HiPuRhoProducer.cc:107
EventSetup.h
HiPuRhoProducer::nTow_
std::vector< int > nTow_
Definition: HiPuRhoProducer.cc:121
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:236
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:108
View.h
HiPuRhoProducer::medianWindowWidth_
const int medianWindowWidth_
Definition: HiPuRhoProducer.cc:93
ParameterSet.h
ParameterSetDescriptionFiller.h
HITowerHelper.h
HiPuRhoProducer::geo_
CaloGeometry const * geo_
Definition: HiPuRhoProducer.cc:135
HiPuRhoProducer::geomtowers_
std::map< int, int > geomtowers_
Definition: HiPuRhoProducer.cc:141
edm::Event
Definition: Event.h:73
HiPuRhoProducer::etaEdgeHi_
std::vector< double > etaEdgeHi_
Definition: HiPuRhoProducer.cc:115
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:142
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:98
findQualityFiles.size
size
Write out results.
Definition: findQualityFiles.py:443
HiPuRhoProducer::ietamin_
int ietamin_
Definition: HiPuRhoProducer.cc:139
HiPuRhoProducer::calculatePedestal
virtual void calculatePedestal(std::vector< fastjet::PseudoJet > const &coll)
Definition: HiPuRhoProducer.cc:276
Generated for CMSSW Reference Manual by doxygen 1.8.16