CMS 3D CMS Logo
JetIDHelper.cc
Go to the documentation of this file.
1 #include "RecoJets/JetProducers/interface/JetIDHelper.h"
2 
3 #include "FWCore/MessageLogger/interface/MessageLogger.h"
4 
5 #include "FWCore/Framework/interface/Frameworkfwd.h"
6 #include "FWCore/Framework/interface/EDProducer.h"
7 #include "FWCore/Framework/interface/Event.h"
8 #include "FWCore/Framework/interface/MakerMacros.h"
9 #include "FWCore/ParameterSet/interface/ParameterSet.h"
10 #include "FWCore/Utilities/interface/InputTag.h"
11 
12 #include "DataFormats/METReco/interface/HcalCaloFlagLabels.h"
13 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
14 // JetIDHelper needs a much more detailed description that the one in HcalTopology,
15 // so to be consistent, all needed constants are hardwired in JetIDHelper.cc itself
16 #include "Geometry/Records/interface/HcalRecNumberingRecord.h"
17 #include "Geometry/CaloTopology/interface/HcalTopology.h"
18 
19 #include "TMath.h"
20 #include <vector>
21 #include <numeric>
22 #include <iostream>
23 
24 using namespace std;
25 
26 // defining stuff useful only for this class (not needed in header)
27 namespace reco {
28 
29  namespace helper {
30 
31  // select2nd exists only in some std and boost implementations, so let's control our own fate
32  // and it can't be a non-static member function.
33  static double select2nd(std::map<int, double>::value_type const &pair) { return pair.second; }
34 
35  bool hasNonPositiveE(reco::helper::JetIDHelper::subtower x) { return x.E <= 0; }
36 
37  bool subtower_has_greater_E(reco::helper::JetIDHelper::subtower i, reco::helper::JetIDHelper::subtower j) {
38  return i.E > j.E;
39  }
40  std::atomic<int> JetIDHelper::sanity_checks_left_{100};
41  } // namespace helper
42 } // namespace reco
43 
44 reco::helper::JetIDHelper::JetIDHelper(edm::ParameterSet const &pset, edm::ConsumesCollector &&iC) {
45  useRecHits_ = pset.getParameter<bool>("useRecHits");
46  if (useRecHits_) {
47  hbheRecHitsColl_ = pset.getParameter<edm::InputTag>("hbheRecHitsColl");
48  hoRecHitsColl_ = pset.getParameter<edm::InputTag>("hoRecHitsColl");
49  hfRecHitsColl_ = pset.getParameter<edm::InputTag>("hfRecHitsColl");
50  ebRecHitsColl_ = pset.getParameter<edm::InputTag>("ebRecHitsColl");
51  eeRecHitsColl_ = pset.getParameter<edm::InputTag>("eeRecHitsColl");
52  }
53  initValues();
54 
55  input_HBHERecHits_token_ = iC.consumes<HBHERecHitCollection>(hbheRecHitsColl_);
56  input_HORecHits_token_ = iC.consumes<HORecHitCollection>(hoRecHitsColl_);
57  input_HFRecHits_token_ = iC.consumes<HFRecHitCollection>(hfRecHitsColl_);
58  input_EBRecHits_token_ = iC.consumes<EBRecHitCollection>(ebRecHitsColl_);
59  input_EERecHits_token_ = iC.consumes<EERecHitCollection>(eeRecHitsColl_);
60 }
61 
62 void reco::helper::JetIDHelper::initValues() {
63  fHPD_ = -1.0;
64  fRBX_ = -1.0;
65  n90Hits_ = -1;
66  fSubDetector1_ = -1.0;
67  fSubDetector2_ = -1.0;
68  fSubDetector3_ = -1.0;
69  fSubDetector4_ = -1.0;
70  restrictedEMF_ = -1.0;
71  nHCALTowers_ = -1;
72  nECALTowers_ = -1;
73  approximatefHPD_ = -1.0;
74  approximatefRBX_ = -1.0;
75  hitsInN90_ = -1;
76  fEB_ = fEE_ = fHB_ = fHE_ = fHO_ = fLong_ = fShort_ = -1.0;
77  fLS_ = fHFOOT_ = -1.0;
78 }
79 
80 void reco::helper::JetIDHelper::fillDescription(edm::ParameterSetDescription &iDesc) {
81  iDesc
82  .ifValue(edm::ParameterDescription<bool>("useRecHits", true, true),
83  true >> (edm::ParameterDescription<edm::InputTag>("hbheRecHitsColl", edm::InputTag(), true) and
84  edm::ParameterDescription<edm::InputTag>("hoRecHitsColl", edm::InputTag(), true) and
85  edm::ParameterDescription<edm::InputTag>("hfRecHitsColl", edm::InputTag(), true) and
86  edm::ParameterDescription<edm::InputTag>("ebRecHitsColl", edm::InputTag(), true) and
87  edm::ParameterDescription<edm::InputTag>("eeRecHitsColl", edm::InputTag(), true)))
88  ->setComment(
89  "If using RecHits to calculate the precise jet ID variables that need them, "
90  "their sources need to be specified");
91 }
92 
93 void reco::helper::JetIDHelper::calculate(const edm::Event &event,
94  const edm::EventSetup &setup,
95  const reco::CaloJet &jet,
96  const int iDbg) {
97  initValues();
98 
99  // --------------------------------------------------
100  // 1) jet ID variable derived from existing fractions
101  // --------------------------------------------------
102 
103  double E_EM = TMath::Max(float(0.), jet.emEnergyInHF()) + TMath::Max(float(0.), jet.emEnergyInEB()) +
104  TMath::Max(float(0.), jet.emEnergyInEE());
105  double E_Had = TMath::Max(float(0.), jet.hadEnergyInHB()) + TMath::Max(float(0.), jet.hadEnergyInHE()) +
106  TMath::Max(float(0.), jet.hadEnergyInHO()) + TMath::Max(float(0.), jet.hadEnergyInHF());
107  if (E_Had + E_EM > 0)
108  restrictedEMF_ = E_EM / (E_EM + E_Had);
109  if (iDbg > 1)
110  cout << "jet pT: " << jet.pt() << ", eT: " << jet.et() << ", E: " << jet.energy() << " rEMF: " << restrictedEMF_
111  << endl;
112 
113  // ------------------------
114  // 2) tower based variables
115  // ------------------------
116  vector<subtower> subtowers, Ecal_subtowers, Hcal_subtowers, HO_subtowers;
117  vector<double> HPD_energies, RBX_energies;
118 
119  classifyJetTowers(
120  event, jet, subtowers, Ecal_subtowers, Hcal_subtowers, HO_subtowers, HPD_energies, RBX_energies, iDbg);
121  if (iDbg > 1) {
122  cout << "E:";
123  for (unsigned int i = 0; i < subtowers.size(); ++i)
124  cout << " " << subtowers[i].E << "," << subtowers[i].Nhit;
125  cout << "\nECal_E:";
126  for (unsigned int i = 0; i < Ecal_subtowers.size(); ++i)
127  cout << " " << Ecal_subtowers[i].E << "," << Ecal_subtowers[i].Nhit;
128  cout << "\nHCal_E:";
129  for (unsigned int i = 0; i < Hcal_subtowers.size(); ++i)
130  cout << " " << Hcal_subtowers[i].E << "," << Hcal_subtowers[i].Nhit;
131  cout << "\nHO_E:";
132  for (unsigned int i = 0; i < HO_subtowers.size(); ++i)
133  cout << " " << HO_subtowers[i].E << "," << HO_subtowers[i].Nhit;
134  cout << "\nHPD_E:";
135  for (unsigned int i = 0; i < HPD_energies.size(); ++i)
136  cout << " " << HPD_energies[i];
137  cout << "\nRBX_E:";
138  for (unsigned int i = 0; i < RBX_energies.size(); ++i)
139  cout << " " << RBX_energies[i];
140  cout << endl;
141  }
142 
143  // counts
144  hitsInN90_ = hitsInNCarrying(0.9, subtowers);
145  nHCALTowers_ = Hcal_subtowers.size();
146  vector<subtower>::const_iterator it;
147  it = find_if(Ecal_subtowers.begin(), Ecal_subtowers.end(), hasNonPositiveE);
148  nECALTowers_ = it - Ecal_subtowers.begin(); // ignores negative energies from HF!
149 
150  // energy fractions
151  double max_HPD_energy = 0., max_RBX_energy = 0.;
152  std::vector<double>::const_iterator it_max_HPD_energy = std::max_element(HPD_energies.begin(), HPD_energies.end());
153  std::vector<double>::const_iterator it_max_RBX_energy = std::max_element(RBX_energies.begin(), RBX_energies.end());
154  if (it_max_HPD_energy != HPD_energies.end()) {
155  max_HPD_energy = *it_max_HPD_energy;
156  }
157  if (it_max_RBX_energy != RBX_energies.end()) {
158  max_RBX_energy = *it_max_RBX_energy;
159  }
160  if (jet.energy() > 0) {
161  if (!HPD_energies.empty())
162  approximatefHPD_ = max_HPD_energy / jet.energy();
163  if (!RBX_energies.empty())
164  approximatefRBX_ = max_HPD_energy / jet.energy();
165  }
166 
167  // -----------------------
168  // 3) cell based variables
169  // -----------------------
170  if (useRecHits_) {
171  vector<double> energies, subdet_energies, Ecal_energies, Hcal_energies, HO_energies;
172  double LS_bad_energy, HF_OOT_energy;
173  classifyJetComponents(event,
174  setup,
175  jet,
176  energies,
177  subdet_energies,
178  Ecal_energies,
179  Hcal_energies,
180  HO_energies,
181  HPD_energies,
182  RBX_energies,
183  LS_bad_energy,
184  HF_OOT_energy,
185  iDbg);
186 
187  // counts
188  n90Hits_ = nCarrying(0.9, energies);
189 
190  // energy fractions
191  if (jet.energy() > 0) {
192  if (!HPD_energies.empty())
193  fHPD_ = max_HPD_energy / jet.energy();
194  if (!RBX_energies.empty())
195  fRBX_ = max_RBX_energy / jet.energy();
196  if (!subdet_energies.empty())
197  fSubDetector1_ = subdet_energies.at(0) / jet.energy();
198  if (subdet_energies.size() > 1)
199  fSubDetector2_ = subdet_energies.at(1) / jet.energy();
200  if (subdet_energies.size() > 2)
201  fSubDetector3_ = subdet_energies.at(2) / jet.energy();
202  if (subdet_energies.size() > 3)
203  fSubDetector4_ = subdet_energies.at(3) / jet.energy();
204  fLS_ = LS_bad_energy / jet.energy();
205  fHFOOT_ = HF_OOT_energy / jet.energy();
206 
207  if (iDbg > 0 && sanity_checks_left_.load(std::memory_order_acquire) > 0) {
208  --sanity_checks_left_;
209  double EH_sum = accumulate(Ecal_energies.begin(), Ecal_energies.end(), 0.);
210  EH_sum = accumulate(Hcal_energies.begin(), Hcal_energies.end(), EH_sum);
211  double EHO_sum = accumulate(HO_energies.begin(), HO_energies.end(), EH_sum);
212  if (jet.energy() > 0.001 && abs(EH_sum / jet.energy() - 1) > 0.01 && abs(EHO_sum / jet.energy() - 1) > 0.01)
213  edm::LogWarning("BadInput") << "Jet energy (" << jet.energy()
214  << ") does not match the total energy in the recHits (" << EHO_sum
215  << ", or without HO: " << EH_sum << ") . Are these the right recHits? "
216  << "Were jet energy corrections mistakenly applied before jet ID? A bug?";
217  if (iDbg > 1)
218  cout << "Sanity check - E: " << jet.energy() << " =? EH_sum: " << EH_sum << " / EHO_sum: " << EHO_sum << endl;
219  }
220  }
221 
222  if (iDbg > 1) {
223  cout << "DBG - fHPD: " << fHPD_ << ", fRBX: " << fRBX_ << ", nh90: " << n90Hits_ << ", fLS: " << fLS_
224  << ", fHFOOT: " << fHFOOT_ << endl;
225  cout << " -~fHPD: " << approximatefHPD_ << ", ~fRBX: " << approximatefRBX_ << ", hits in n90: " << hitsInN90_
226  << endl;
227  cout << " - nHCALTowers: " << nHCALTowers_ << ", nECALTowers: " << nECALTowers_
228  << "; subdet fractions: " << fSubDetector1_ << ", " << fSubDetector2_ << ", " << fSubDetector3_ << ", "
229  << fSubDetector4_ << endl;
230  }
231  }
232 }
233 
234 unsigned int reco::helper::JetIDHelper::nCarrying(double fraction, const std::vector<double> &descending_energies) {
235  double totalE = 0;
236  for (unsigned int i = 0; i < descending_energies.size(); ++i)
237  totalE += descending_energies[i];
238 
239  double runningE = 0;
240  unsigned int NC = descending_energies.size();
241 
242  // slightly odd loop structure avoids round-off problems when runningE never catches up with totalE
243  for (unsigned int i = descending_energies.size(); i > 0; --i) {
244  runningE += descending_energies[i - 1];
245  if (runningE < (1 - fraction) * totalE)
246  NC = i - 1;
247  }
248  return NC;
249 }
250 
251 unsigned int reco::helper::JetIDHelper::hitsInNCarrying(double fraction,
252  const std::vector<subtower> &descending_towers) {
253  double totalE = 0;
254  for (unsigned int i = 0; i < descending_towers.size(); ++i)
255  totalE += descending_towers[i].E;
256 
257  double runningE = 0;
258  unsigned int NH = 0;
259 
260  // slightly odd loop structure avoids round-off problems when runningE never catches up with totalE
261  for (unsigned int i = descending_towers.size(); i > 0; --i) {
262  runningE += descending_towers[i - 1].E;
263  if (runningE >= (1 - fraction) * totalE)
264  NH += descending_towers[i - 1].Nhit;
265  }
266  return NH;
267 }
268 
269 void reco::helper::JetIDHelper::classifyJetComponents(const edm::Event &event,
270  const edm::EventSetup &setup,
271  const reco::CaloJet &jet,
272  vector<double> &energies,
273  vector<double> &subdet_energies,
274  vector<double> &Ecal_energies,
275  vector<double> &Hcal_energies,
276  vector<double> &HO_energies,
277  vector<double> &HPD_energies,
278  vector<double> &RBX_energies,
279  double &LS_bad_energy,
280  double &HF_OOT_energy,
281  const int iDbg) {
282  energies.clear();
283  subdet_energies.clear();
284  Ecal_energies.clear();
285  Hcal_energies.clear();
286  HO_energies.clear();
287  HPD_energies.clear();
288  RBX_energies.clear();
289  LS_bad_energy = HF_OOT_energy = 0.;
290 
291  edm::ESHandle<HcalTopology> theHcalTopology;
292  setup.get<HcalRecNumberingRecord>().get(theHcalTopology);
293 
294  std::map<int, double> HPD_energy_map, RBX_energy_map;
295  vector<double> EB_energies, EE_energies, HB_energies, HE_energies, short_energies, long_energies;
296  edm::Handle<HBHERecHitCollection> HBHERecHits;
297  edm::Handle<HORecHitCollection> HORecHits;
298  edm::Handle<HFRecHitCollection> HFRecHits;
299  edm::Handle<EBRecHitCollection> EBRecHits;
300  edm::Handle<EERecHitCollection> EERecHits;
301  // the jet only contains DetIds, so first read recHit collection
302  event.getByToken(input_HBHERecHits_token_, HBHERecHits);
303  event.getByToken(input_HORecHits_token_, HORecHits);
304  event.getByToken(input_HFRecHits_token_, HFRecHits);
305  event.getByToken(input_EBRecHits_token_, EBRecHits);
306  event.getByToken(input_EERecHits_token_, EERecHits);
307  if (iDbg > 2)
308  cout << "# of rechits found - HBHE: " << HBHERecHits->size() << ", HO: " << HORecHits->size()
309  << ", HF: " << HFRecHits->size() << ", EB: " << EBRecHits->size() << ", EE: " << EERecHits->size() << endl;
310 
311  vector<CaloTowerPtr> towers = jet.getCaloConstituents();
312  int nTowers = towers.size();
313  if (iDbg > 9)
314  cout << "In classifyJetComponents. # of towers found: " << nTowers << endl;
315 
316  for (int iTower = 0; iTower < nTowers; iTower++) {
317  CaloTowerPtr &tower = towers[iTower];
318 
319  int nCells = tower->constituentsSize();
320  if (iDbg)
321  cout << "tower #" << iTower << " has " << nCells << " cells. "
322  << "It's at iEta: " << tower->ieta() << ", iPhi: " << tower->iphi() << endl;
323 
324  const vector<DetId> &cellIDs = tower->constituents(); // cell == recHit
325 
326  for (int iCell = 0; iCell < nCells; ++iCell) {
327  DetId::Detector detNum = cellIDs[iCell].det();
328  if (detNum == DetId::Hcal) {
329  HcalDetId HcalID = cellIDs[iCell];
330  HcalSubdetector HcalNum = HcalID.subdet();
331  double hitE = 0;
332  if (HcalNum == HcalOuter) {
333  HORecHitCollection::const_iterator theRecHit = HORecHits->find(HcalID);
334  if (theRecHit == HORecHits->end()) {
335  edm::LogWarning("UnexpectedEventContents") << "Can't find the HO recHit with ID: " << HcalID;
336  continue;
337  }
338  hitE = theRecHit->energy();
339  HO_energies.push_back(hitE);
340 
341  } else if (HcalNum == HcalForward) {
342  HFRecHitCollection::const_iterator theRecHit = HFRecHits->find(HcalID);
343  if (theRecHit == HFRecHits->end()) {
344  edm::LogWarning("UnexpectedEventContents") << "Can't find the HF recHit with ID: " << HcalID;
345  continue;
346  }
347  hitE = theRecHit->energy();
348  if (iDbg > 4)
349  cout << "hit #" << iCell << " is HF , E: " << hitE << " iEta: " << theRecHit->id().ieta()
350  << ", depth: " << theRecHit->id().depth() << ", iPhi: " << theRecHit->id().iphi();
351 
352  if (HcalID.depth() == 1)
353  long_energies.push_back(hitE);
354  else
355  short_energies.push_back(hitE);
356 
357  uint32_t flags = theRecHit->flags();
358  if (flags & (1 << HcalCaloFlagLabels::HFLongShort))
359  LS_bad_energy += hitE;
360  if (flags & ((1 << HcalCaloFlagLabels::HFDigiTime) | (3 << HcalCaloFlagLabels::HFTimingTrustBits) |
361  (1 << HcalCaloFlagLabels::TimingSubtractedBit) | (1 << HcalCaloFlagLabels::TimingAddedBit) |
362  (1 << HcalCaloFlagLabels::TimingErrorBit)))
363  HF_OOT_energy += hitE;
364  if (iDbg > 4 && flags)
365  cout << "flags: " << flags << " -> LS_bad_energy: " << LS_bad_energy << ", HF_OOT_energy: " << HF_OOT_energy
366  << endl;
367 
368  } else { // HBHE
369 
370  HBHERecHitCollection::const_iterator theRecHit = HBHERecHits->find(HcalID);
371  if (theRecHit == HBHERecHits->end()) {
372  edm::LogWarning("UnexpectedEventContents") << "Can't find the HBHE recHit with ID: " << HcalID;
373  continue;
374  }
375  hitE = theRecHit->energy();
376  int iEta = theRecHit->id().ieta();
377  int depth = theRecHit->id().depth();
378  Region region = HBHE_region(theRecHit->id().rawId());
379  int hitIPhi = theRecHit->id().iphi();
380  if (iDbg > 3)
381  cout << "hit #" << iCell << " is HBHE, E: " << hitE << " iEta: " << iEta << ", depth: " << depth
382  << ", iPhi: " << theRecHit->id().iphi() << " -> " << region;
383 
384  if (theHcalTopology->mergedDepth29(theRecHit->id()))
385  hitE /= 2; // Depth 3 at the HE forward edge is split over tower 28 & 29, and jet reco. assigns half each
386 
387  int iHPD = 100 * region;
388  int iRBX = 100 * region + ((hitIPhi + 1) % 72) / 4; // 71,72,1,2 are in the same RBX module
389 
390  if (std::abs(iEta) >= 21) {
391  if ((0x1 & hitIPhi) == 0) {
392  edm::LogError("CodeAssumptionsViolated") << "Bug?! Jet ID code assumes no even iPhi recHits at HE edges";
393  return;
394  }
395  bool oddnessIEta = HBHE_oddness(iEta, depth);
396  bool upperIPhi = ((hitIPhi % 4) == 1 || (hitIPhi % 4) == 2); // the upper iPhi indices in the HE wedge
397  // remap the iPhi so it fits the one in the inner HE regions, change in needed in two cases:
398  // 1) in the upper iPhis of the module, the even iEtas belong to the higher iPhi
399  // 2) in the loewr iPhis of the module, the odd iEtas belong to the higher iPhi
400  if (upperIPhi != oddnessIEta)
401  ++hitIPhi;
402  // note that hitIPhi could not be 72 before, so it's still in the legal range [1,72]
403  }
404  iHPD += hitIPhi;
405 
406  // book the energies
407  HPD_energy_map[iHPD] += hitE;
408  RBX_energy_map[iRBX] += hitE;
409  if (iDbg > 5)
410  cout << " --> H[" << iHPD << "]=" << HPD_energy_map[iHPD] << ", R[" << iRBX << "]=" << RBX_energy_map[iRBX];
411  if (iDbg > 1)
412  cout << endl;
413 
414  if (region == HBneg || region == HBpos)
415  HB_energies.push_back(hitE);
416  else
417  HE_energies.push_back(hitE);
418 
419  } // if HBHE
420  if (hitE == 0)
421  edm::LogWarning("UnexpectedEventContents") << "HCal hitE==0? (or unknown subdetector?)";
422 
423  } // if HCAL
424 
425  else if (detNum == DetId::Ecal) {
426  int EcalNum = cellIDs[iCell].subdetId();
427  double hitE = 0;
428  if (EcalNum == 1) {
429  EBDetId EcalID = cellIDs[iCell];
430  EBRecHitCollection::const_iterator theRecHit = EBRecHits->find(EcalID);
431  if (theRecHit == EBRecHits->end()) {
432  edm::LogWarning("UnexpectedEventContents") << "Can't find the EB recHit with ID: " << EcalID;
433  continue;
434  }
435  hitE = theRecHit->energy();
436  EB_energies.push_back(hitE);
437  } else if (EcalNum == 2) {
438  EEDetId EcalID = cellIDs[iCell];
439  EERecHitCollection::const_iterator theRecHit = EERecHits->find(EcalID);
440  if (theRecHit == EERecHits->end()) {
441  edm::LogWarning("UnexpectedEventContents") << "Can't find the EE recHit with ID: " << EcalID;
442  continue;
443  }
444  hitE = theRecHit->energy();
445  EE_energies.push_back(hitE);
446  }
447  if (hitE == 0)
448  edm::LogWarning("UnexpectedEventContents") << "ECal hitE==0? (or unknown subdetector?)";
449  if (iDbg > 6)
450  cout << "EcalNum: " << EcalNum << " hitE: " << hitE << endl;
451  } //
452  } // loop on cells
453  } // loop on towers
454 
455  /* Disabling check until HO is accounted for in EMF. Check was used in CMSSW_2, where HE was excluded.
456  double expHcalE = jet.energy() * (1-jet.emEnergyFraction());
457  if( totHcalE + expHcalE > 0 &&
458  TMath::Abs( totHcalE - expHcalE ) > 0.01 &&
459  ( totHcalE - expHcalE ) / ( totHcalE + expHcalE ) > 0.0001 ) {
460  edm::LogWarning("CodeAssumptionsViolated")<<"failed to account for all Hcal energies"
461  <<totHcalE<<"!="<<expHcalE;
462  } */
463  // concatenate Hcal and Ecal lists
464  Hcal_energies.insert(Hcal_energies.end(), HB_energies.begin(), HB_energies.end());
465  Hcal_energies.insert(Hcal_energies.end(), HE_energies.begin(), HE_energies.end());
466  Hcal_energies.insert(Hcal_energies.end(), short_energies.begin(), short_energies.end());
467  Hcal_energies.insert(Hcal_energies.end(), long_energies.begin(), long_energies.end());
468  Ecal_energies.insert(Ecal_energies.end(), EB_energies.begin(), EB_energies.end());
469  Ecal_energies.insert(Ecal_energies.end(), EE_energies.begin(), EE_energies.end());
470 
471  // sort the energies
472  // std::sort( Hcal_energies.begin(), Hcal_energies.end(), greater<double>() );
473  // std::sort( Ecal_energies.begin(), Ecal_energies.end(), greater<double>() );
474 
475  // put the energy sums (the 2nd entry in each pair of the maps) into the output vectors and sort them
476  std::transform(
477  HPD_energy_map.begin(), HPD_energy_map.end(), std::inserter(HPD_energies, HPD_energies.end()), select2nd);
478  // std::select2nd<std::map<int,double>::value_type>());
479  // std::sort( HPD_energies.begin(), HPD_energies.end(), greater<double>() );
480  std::transform(
481  RBX_energy_map.begin(), RBX_energy_map.end(), std::inserter(RBX_energies, RBX_energies.end()), select2nd);
482  // std::select2nd<std::map<int,double>::value_type>());
483  // std::sort( RBX_energies.begin(), RBX_energies.end(), greater<double>() );
484 
485  energies.insert(energies.end(), Hcal_energies.begin(), Hcal_energies.end());
486  energies.insert(energies.end(), Ecal_energies.begin(), Ecal_energies.end());
487  energies.insert(energies.end(), HO_energies.begin(), HO_energies.end());
488  std::sort(energies.begin(), energies.end(), greater<double>());
489 
490  // prepare sub detector energies, then turn them into fractions
491  fEB_ = std::accumulate(EB_energies.begin(), EB_energies.end(), 0.);
492  fEE_ = std::accumulate(EE_energies.begin(), EE_energies.end(), 0.);
493  fHB_ = std::accumulate(HB_energies.begin(), HB_energies.end(), 0.);
494  fHE_ = std::accumulate(HE_energies.begin(), HE_energies.end(), 0.);
495  fHO_ = std::accumulate(HO_energies.begin(), HO_energies.end(), 0.);
496  fShort_ = std::accumulate(short_energies.begin(), short_energies.end(), 0.);
497  fLong_ = std::accumulate(long_energies.begin(), long_energies.end(), 0.);
498  subdet_energies.push_back(fEB_);
499  subdet_energies.push_back(fEE_);
500  subdet_energies.push_back(fHB_);
501  subdet_energies.push_back(fHE_);
502  subdet_energies.push_back(fHO_);
503  subdet_energies.push_back(fShort_);
504  subdet_energies.push_back(fLong_);
505  std::sort(subdet_energies.begin(), subdet_energies.end(), greater<double>());
506  if (jet.energy() > 0) {
507  fEB_ /= jet.energy();
508  fEE_ /= jet.energy();
509  fHB_ /= jet.energy();
510  fHE_ /= jet.energy();
511  fHO_ /= jet.energy();
512  fShort_ /= jet.energy();
513  fLong_ /= jet.energy();
514  } else {
515  if (fEB_ > 0 || fEE_ > 0 || fHB_ > 0 || fHE_ > 0 || fHO_ > 0 || fShort_ > 0 || fLong_ > 0)
516  edm::LogError("UnexpectedEventContents") << "Jet ID Helper found energy in subdetectors and jet E <= 0";
517  }
518 }
519 
520 void reco::helper::JetIDHelper::classifyJetTowers(const edm::Event &event,
521  const reco::CaloJet &jet,
522  vector<subtower> &subtowers,
523  vector<subtower> &Ecal_subtowers,
524  vector<subtower> &Hcal_subtowers,
525  vector<subtower> &HO_subtowers,
526  vector<double> &HPD_energies,
527  vector<double> &RBX_energies,
528  const int iDbg) {
529  subtowers.clear();
530  Ecal_subtowers.clear();
531  Hcal_subtowers.clear();
532  HO_subtowers.clear();
533  HPD_energies.clear();
534  RBX_energies.clear();
535 
536  std::map<int, double> HPD_energy_map, RBX_energy_map;
537 
538  vector<CaloTowerPtr> towers = jet.getCaloConstituents();
539  int nTowers = towers.size();
540  if (iDbg > 9)
541  cout << "classifyJetTowers started. # of towers found: " << nTowers << endl;
542 
543  for (int iTower = 0; iTower < nTowers; iTower++) {
544  CaloTowerPtr &tower = towers[iTower];
545 
546  int nEM = 0, nHad = 0, nHO = 0;
547  const vector<DetId> &cellIDs = tower->constituents(); // cell == recHit
548  int nCells = cellIDs.size();
549  if (iDbg)
550  cout << "tower #" << iTower << " has " << nCells << " cells. "
551  << "It's at iEta: " << tower->ieta() << ", iPhi: " << tower->iphi() << endl;
552 
553  for (int iCell = 0; iCell < nCells; ++iCell) {
554  DetId::Detector detNum = cellIDs[iCell].det();
555  if (detNum == DetId::Hcal) {
556  HcalDetId HcalID = cellIDs[iCell];
557  HcalSubdetector HcalNum = HcalID.subdet();
558  if (HcalNum == HcalOuter) {
559  ++nHO;
560  } else {
561  ++nHad;
562  }
563  } else if (detNum == DetId::Ecal) {
564  ++nEM;
565  }
566  }
567 
568  double E_em = tower->emEnergy();
569  if (E_em != 0)
570  Ecal_subtowers.push_back(subtower(E_em, nEM));
571 
572  double E_HO = tower->outerEnergy();
573  if (E_HO != 0)
574  HO_subtowers.push_back(subtower(E_HO, nHO));
575 
576  double E_had = tower->hadEnergy();
577  if (E_had != 0) {
578  Hcal_subtowers.push_back(subtower(E_had, nHad));
579  // totHcalE += E_had;
580 
581  int iEta = tower->ieta();
582  Region reg = region(iEta);
583  int iPhi = tower->iphi();
584  if (iDbg > 3)
585  cout << "tower has E_had: " << E_had << " iEta: " << iEta << ", iPhi: " << iPhi << " -> " << reg;
586 
587  if (reg == HEneg || reg == HBneg || reg == HBpos || reg == HEpos) {
588  int oddnessIEta = HBHE_oddness(iEta);
589  if (oddnessIEta < 0)
590  break; // can't assign this tower to a single readout component
591 
592  int iHPD = 100 * reg;
593  int iRBX = 100 * reg + ((iPhi + 1) % 72) / 4; // 71,72,1,2 are in the same RBX module
594 
595  if ((reg == HEneg || reg == HEpos) && std::abs(iEta) >= 21) { // at low-granularity edge of HE
596  if ((0x1 & iPhi) == 0) {
597  edm::LogError("CodeAssumptionsViolated") << "Bug?! Jet ID code assumes no even iPhi recHits at HE edges";
598  return;
599  }
600  bool boolOddnessIEta = oddnessIEta;
601  bool upperIPhi = ((iPhi % 4) == 1 || (iPhi % 4) == 2); // the upper iPhi indices in the HE wedge
602  // remap the iPhi so it fits the one in the inner HE regions, change in needed in two cases:
603  // 1) in the upper iPhis of the module, the even IEtas belong to the higher iPhi
604  // 2) in the loewr iPhis of the module, the odd IEtas belong to the higher iPhi
605  if (upperIPhi != boolOddnessIEta)
606  ++iPhi;
607  // note that iPhi could not be 72 before, so it's still in the legal range [1,72]
608  } // if at low-granularity edge of HE
609  iHPD += iPhi;
610 
611  // book the energies
612  HPD_energy_map[iHPD] += E_had;
613  RBX_energy_map[iRBX] += E_had;
614  if (iDbg > 5)
615  cout << " --> H[" << iHPD << "]=" << HPD_energy_map[iHPD] << ", R[" << iRBX << "]=" << RBX_energy_map[iRBX];
616  } // HBHE
617  } // E_had > 0
618  } // loop on towers
619 
620  // sort the subtowers
621  // std::sort( Hcal_subtowers.begin(), Hcal_subtowers.end(), subtower_has_greater_E );
622  // std::sort( Ecal_subtowers.begin(), Ecal_subtowers.end(), subtower_has_greater_E );
623 
624  // put the energy sums (the 2nd entry in each pair of the maps) into the output vectors and sort them
625  std::transform(
626  HPD_energy_map.begin(), HPD_energy_map.end(), std::inserter(HPD_energies, HPD_energies.end()), select2nd);
627  // std::select2nd<std::map<int,double>::value_type>());
628  // std::sort( HPD_energies.begin(), HPD_energies.end(), greater<double>() );
629  std::transform(
630  RBX_energy_map.begin(), RBX_energy_map.end(), std::inserter(RBX_energies, RBX_energies.end()), select2nd);
631  // std::select2nd<std::map<int,double>::value_type>());
632  // std::sort( RBX_energies.begin(), RBX_energies.end(), greater<double>() );
633 
634  subtowers.insert(subtowers.end(), Hcal_subtowers.begin(), Hcal_subtowers.end());
635  subtowers.insert(subtowers.end(), Ecal_subtowers.begin(), Ecal_subtowers.end());
636  subtowers.insert(subtowers.end(), HO_subtowers.begin(), HO_subtowers.end());
637  std::sort(subtowers.begin(), subtowers.end(), subtower_has_greater_E);
638 }
639 
640 // ---------------------------------------------------------------------------------------------------------
641 // private helper functions to figure out detector & readout geometry as needed
642 
643 int reco::helper::JetIDHelper::HBHE_oddness(int iEta, int depth) {
644  int ae = TMath::Abs(iEta);
645  if (ae == 29 && depth == 1)
646  ae += 1; // observed that: hits are at depths 1 & 2; 1 goes with the even pattern
647  return ae & 0x1;
648 }
649 
650 reco::helper::JetIDHelper::Region reco::helper::JetIDHelper::HBHE_region(uint32_t rawid) {
651  HcalDetId id(rawid);
652  if (id.subdet() == HcalEndcap) {
653  if (id.ieta() < 0)
654  return HEneg;
655  else
656  return HEpos;
657  }
658  if (id.subdet() == HcalBarrel && id.ieta() < 0)
659  return HBneg;
660  return HBpos;
661 }
662 
663 int reco::helper::JetIDHelper::HBHE_oddness(int iEta) {
664  int ae = TMath::Abs(iEta);
665  if (ae == 29)
666  return -1; // can't figure it out without RecHits
667  return ae & 0x1;
668 }
669 
670 reco::helper::JetIDHelper::Region reco::helper::JetIDHelper::region(int iEta) {
671  if (iEta == 16 || iEta == -16)
672  return unknown_region; // both HB and HE cells belong to these towers
673  if (iEta == 29 || iEta == -29)
674  return unknown_region; // both HE and HF cells belong to these towers
675  if (iEta <= -30)
676  return HFneg;
677  if (iEta >= 30)
678  return HFpos;
679  if (iEta <= -17)
680  return HEneg;
681  if (iEta >= 17)
682  return HEpos;
683  if (iEta < 0)
684  return HBneg;
685  return HBpos;
686 }
HLT_FULL_cff.towers
towers
Definition: HLT_FULL_cff.py:36362
reco::CaloJet
Jets made from CaloTowers.
Definition: CaloJet.h:27
HcalCaloFlagLabels::TimingErrorBit
Definition: HcalCaloFlagLabels.h:52
HcalCaloFlagLabels::TimingSubtractedBit
Definition: HcalCaloFlagLabels.h:50
mps_fire.i
i
Definition: mps_fire.py:428
edm::SortedCollection::const_iterator
std::vector< T >::const_iterator const_iterator
Definition: SortedCollection.h:80
MessageLogger.h
EDProducer.h
HcalTopology::mergedDepth29
bool mergedDepth29(HcalDetId id) const
Definition: HcalTopology.h:111
HcalCaloFlagLabels::HFTimingTrustBits
Definition: HcalCaloFlagLabels.h:32
egHLT::errCodes::HFRecHits
Definition: EgHLTErrCodes.h:21
EBDetId
Definition: EBDetId.h:17
egHLT::errCodes::HBHERecHits
Definition: EgHLTErrCodes.h:20
gather_cfg.cout
cout
Definition: gather_cfg.py:144
reco::helper::JetIDHelper::calculate
void calculate(const edm::Event &event, const edm::EventSetup &setup, const reco::CaloJet &jet, const int iDbg=0)
Definition: JetIDHelper.cc:93
reco::helper::subtower_has_greater_E
bool subtower_has_greater_E(reco::helper::JetIDHelper::subtower i, reco::helper::JetIDHelper::subtower j)
Definition: JetIDHelper.cc:37
DetId::Hcal
Definition: DetId.h:28
edm::ParameterSetDescription
Definition: ParameterSetDescription.h:52
hgcalTowerProducer_cfi.tower
tower
Definition: hgcalTowerProducer_cfi.py:4
reco::helper::JetIDHelper::nCarrying
unsigned int nCarrying(double fraction, const std::vector< double > &descending_energies)
Definition: JetIDHelper.cc:234
reco::helper::JetIDHelper::region
Region region(int iEta)
Definition: JetIDHelper.cc:670
edm::SortedCollection
Definition: SortedCollection.h:49
edm::SortedCollection::size
size_type size() const
Definition: SortedCollection.h:215
HcalRecNumberingRecord.h
HcalDetId::depth
constexpr int depth() const
get the tower depth
Definition: HcalDetId.h:164
reco
fixed size matrix
Definition: AlignmentAlgorithmBase.h:45
HcalBarrel
Definition: HcalAssistant.h:33
JetIDHelper.h
reco::helper::hasNonPositiveE
bool hasNonPositiveE(reco::helper::JetIDHelper::subtower x)
Definition: JetIDHelper.cc:35
edm::Handle
Definition: AssociativeIterator.h:50
filterRecHits_cfi.EERecHits
EERecHits
Definition: filterRecHits_cfi.py:9
edm::LogWarning
Log< level::Warning, false > LogWarning
Definition: MessageLogger.h:122
singleTopDQM_cfi.setup
setup
Definition: singleTopDQM_cfi.py:37
MakerMacros.h
testProducerWithPsetDescEmpty_cfi.x1
x1
Definition: testProducerWithPsetDescEmpty_cfi.py:33
Abs
T Abs(T a)
Definition: MathUtil.h:49
HLT_FULL_cff.fraction
fraction
Definition: HLT_FULL_cff.py:52795
reco::helper::select2nd
static double select2nd(std::map< int, double >::value_type const &pair)
Definition: JetIDHelper.cc:33
HcalCaloFlagLabels::HFDigiTime
Definition: HcalCaloFlagLabels.h:38
edm::ESHandle< HcalTopology >
HcalOuter
Definition: HcalAssistant.h:35
HcalDetIdTransform::transform
unsigned transform(const HcalDetId &id, unsigned transformCode)
Definition: HcalDetIdTransform.cc:7
LEDCalibrationChannels.depth
depth
Definition: LEDCalibrationChannels.py:65
HcalRecNumberingRecord
Definition: HcalRecNumberingRecord.h:23
EEDetId
Definition: EEDetId.h:14
reco::helper::JetIDHelper::HBHE_region
Region HBHE_region(uint32_t)
Definition: JetIDHelper.cc:650
reco::helper::JetIDHelper::Region
Region
Definition: JetIDHelper.h:107
LEDCalibrationChannels.ieta
ieta
Definition: LEDCalibrationChannels.py:63
edm::ParameterSet
Definition: ParameterSet.h:47
HcalCaloFlagLabels.h
Event.h
reco::helper::JetIDHelper::subtower
Definition: JetIDHelper.h:68
reco::helper::JetIDHelper::fillDescription
void fillDescription(edm::ParameterSetDescription &iDesc)
Definition: JetIDHelper.cc:80
edm::SortedCollection::end
const_iterator end() const
Definition: SortedCollection.h:267
HcalDetId.h
helper
Definition: helper.py:1
edm::ParameterSetDescription::ifValue
ParameterDescriptionNode * ifValue(ParameterDescription< T > const &switchParameter, std::unique_ptr< ParameterDescriptionCases< T >> cases)
Definition: ParameterSetDescription.h:220
HcalDetId::subdet
constexpr HcalSubdetector subdet() const
get the subdetector
Definition: HcalDetId.h:138
HLT_FULL_cff.region
region
Definition: HLT_FULL_cff.py:88272
HcalDetId
Definition: HcalDetId.h:12
DetId::Detector
Detector
Definition: DetId.h:24
Max
T Max(T a, T b)
Definition: MathUtil.h:44
edm::EventSetup
Definition: EventSetup.h:57
reco::JetExtendedAssociation::value_type
Container::value_type value_type
Definition: JetExtendedAssociation.h:30
edm::LogError
Log< level::Error, false > LogError
Definition: MessageLogger.h:123
DetId::Ecal
Definition: DetId.h:27
get
#define get
reco::helper::JetIDHelper::classifyJetTowers
void classifyJetTowers(const edm::Event &event, const reco::CaloJet &jet, std::vector< subtower > &subtowers, std::vector< subtower > &Ecal_subtowers, std::vector< subtower > &Hcal_subtowers, std::vector< subtower > &HO_subtowers, std::vector< double > &HPD_energies, std::vector< double > &RBX_energies, const int iDbg=0)
Definition: JetIDHelper.cc:520
InputTag.h
edm::Ptr
Definition: AssociationVector.h:31
HcalSubdetector
HcalSubdetector
Definition: HcalAssistant.h:31
HcalForward
Definition: HcalAssistant.h:36
edm::SortedCollection::find
iterator find(key_type k)
Definition: SortedCollection.h:240
HcalTopology.h
std
Definition: JetResolutionObject.h:76
HcalEndcap
Definition: HcalAssistant.h:34
Frameworkfwd.h
triggerObjects_cff.id
id
Definition: triggerObjects_cff.py:31
metsig::jet
Definition: SignAlgoResolutions.h:47
HcalCaloFlagLabels::HFLongShort
Definition: HcalCaloFlagLabels.h:37
reco::helper::JetIDHelper::initValues
void initValues()
Definition: JetIDHelper.cc:62
L1TowerCalibrationProducer_cfi.iEta
iEta
Definition: L1TowerCalibrationProducer_cfi.py:60
reco::helper::JetIDHelper::HBHE_oddness
int HBHE_oddness(int iEta, int depth)
Definition: JetIDHelper.cc:643
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
ParameterSet.h
edm::ParameterDescriptionNode::setComment
void setComment(std::string const &value)
Definition: ParameterDescriptionNode.cc:106
reco::helper::JetIDHelper::hitsInNCarrying
unsigned int hitsInNCarrying(double fraction, const std::vector< subtower > &descending_towers)
Definition: JetIDHelper.cc:251
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
event
Definition: event.py:1
HLT_FULL_cff.flags
flags
Definition: HLT_FULL_cff.py:13150
edm::Event
Definition: Event.h:73
ecaldqm::nTowers
Definition: EcalDQMCommonUtils.h:115
reco::helper::JetIDHelper::classifyJetComponents
void classifyJetComponents(const edm::Event &event, const edm::EventSetup &setup, const reco::CaloJet &jet, std::vector< double > &energies, std::vector< double > &subdet_energies, std::vector< double > &Ecal_energies, std::vector< double > &Hcal_energies, std::vector< double > &HO_energies, std::vector< double > &HPD_energies, std::vector< double > &RBX_energies, double &LS_bad_energy, double &HF_OOT_energy, const int iDbg=0)
Definition: JetIDHelper.cc:269
HcalCaloFlagLabels::TimingAddedBit
Definition: HcalCaloFlagLabels.h:51
edm::ParameterDescription
Definition: ParameterDescription.h:110
reco::helper::JetIDHelper::JetIDHelper
JetIDHelper()
Definition: JetIDHelper.h:26
edm::InputTag
Definition: InputTag.h:15
edm::ConsumesCollector
Definition: ConsumesCollector.h:45
muonDTDigis_cfi.pset
pset
Definition: muonDTDigis_cfi.py:27
filterRecHits_cfi.EBRecHits
EBRecHits
Definition: filterRecHits_cfi.py:8
Generated for CMSSW Reference Manual by doxygen 1.8.16