CMS 3D CMS Logo

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