40 std::atomic<int> JetIDHelper::sanity_checks_left_{100};
77 approximatefHPD_ = -1.0;
78 approximatefRBX_ = -1.0;
80 fEB_ = fEE_ = fHB_ = fHE_ = fHO_ = fLong_ = fShort_ = -1.0;
81 fLS_ = fHFOOT_ = -1.0;
93 )->
setComment(
"If using RecHits to calculate the precise jet ID variables that need them, " 94 "their sources need to be specified");
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;
120 vector<subtower> subtowers, Ecal_subtowers, Hcal_subtowers, HO_subtowers;
121 vector<double> HPD_energies, RBX_energies;
123 classifyJetTowers( event, jet,
124 subtowers, Ecal_subtowers, Hcal_subtowers, HO_subtowers,
125 HPD_energies, RBX_energies, iDbg );
128 for (
unsigned int i=0;
i<subtowers.size(); ++
i)
cout<<
" "<<subtowers[
i].E<<
","<<subtowers[
i].Nhit;
130 for (
unsigned int i=0;
i<Ecal_subtowers.size(); ++
i)
cout<<
" "<<Ecal_subtowers[
i].E<<
","<<Ecal_subtowers[
i].Nhit;
132 for (
unsigned int i=0;
i<Hcal_subtowers.size(); ++
i)
cout<<
" "<<Hcal_subtowers[
i].E<<
","<<Hcal_subtowers[
i].Nhit;
134 for (
unsigned int i=0;
i<HO_subtowers.size(); ++
i)
cout<<
" "<<HO_subtowers[
i].E<<
","<<HO_subtowers[
i].Nhit;
136 for (
unsigned int i=0;
i<HPD_energies.size(); ++
i)
cout<<
" "<<HPD_energies[
i];
138 for (
unsigned int i=0;
i<RBX_energies.size(); ++
i)
cout<<
" "<<RBX_energies[
i];
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();
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;
156 if ( it_max_RBX_energy != RBX_energies.end() ) {
157 max_RBX_energy = *it_max_RBX_energy;
160 if( !HPD_energies.empty() ) approximatefHPD_ = max_HPD_energy / jet.
energy();
161 if( !RBX_energies.empty() ) approximatefRBX_ = max_HPD_energy / jet.
energy();
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 );
175 n90Hits_ = nCarrying( 0.9, energies );
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();
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 );
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;
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;
216 for(
unsigned int i = 0;
i < descending_energies.size(); ++
i ) totalE += descending_energies[
i ];
219 unsigned int NC = descending_energies.size();
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;
233 for(
unsigned int i = 0;
i < descending_towers.size(); ++
i ) totalE += descending_towers[
i ].E;
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;
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 )
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.;
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;
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;
279 if( iDbg > 9 )
cout<<
"In classifyJetComponents. # of towers found: "<<nTowers<<endl;
281 for(
int iTower = 0; iTower <
nTowers ; iTower++ ) {
286 if( iDbg )
cout<<
"tower #"<<iTower<<
" has "<<nCells<<
" cells. " 287 <<
"It's at iEta: "<<tower->
ieta()<<
", iPhi: "<<tower->
iphi()<<endl;
291 for(
int iCell = 0; iCell < nCells; ++iCell ) {
299 if (theRecHit == HORecHits->
end()) {
300 edm::LogWarning(
"UnexpectedEventContents")<<
"Can't find the HO recHit with ID: "<<HcalID;
303 hitE = theRecHit->energy();
304 HO_energies.push_back( hitE );
309 if( theRecHit == HFRecHits->
end() ) {
310 edm::LogWarning(
"UnexpectedEventContents")<<
"Can't find the HF recHit with ID: "<<HcalID;
313 hitE = theRecHit->energy();
315 <<
"hit #"<<iCell<<
" is HF , E: "<<hitE<<
" iEta: "<<theRecHit->id().ieta()
316 <<
", depth: "<<theRecHit->id().depth()<<
", iPhi: " 317 <<theRecHit->id().iphi();
319 if( HcalID.
depth() == 1 )
320 long_energies.push_back( hitE );
322 short_energies.push_back( hitE );
324 uint32_t
flags = theRecHit->flags();
331 if( iDbg>4 && flags )
cout<<
"flags: "<<flags
332 <<
" -> LS_bad_energy: "<<LS_bad_energy<<
", HF_OOT_energy: "<<HF_OOT_energy<<endl;
337 if( theRecHit == HBHERecHits->
end() ) {
338 edm::LogWarning(
"UnexpectedEventContents")<<
"Can't find the HBHE recHit with ID: "<<HcalID;
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()
349 int absIEta =
TMath::Abs( theRecHit->id().ieta() );
350 if( depth == 3 && (absIEta == 28 || absIEta == 29) ) {
353 int iHPD = 100 * region;
354 int iRBX = 100 * region + ((hitIPhi + 1) % 72) / 4;
357 if( (0
x1 & hitIPhi) == 0 ) {
358 edm::LogError(
"CodeAssumptionsViolated")<<
"Bug?! Jet ID code assumes no even iPhi recHits at HE edges";
361 bool oddnessIEta = HBHE_oddness( iEta, depth );
362 bool upperIPhi = (( hitIPhi%4 ) == 1 || ( hitIPhi%4 ) == 2);
366 if( upperIPhi != oddnessIEta ) ++hitIPhi;
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;
378 if( region == HBneg || region == HBpos )
379 HB_energies.push_back( hitE );
381 HE_energies.push_back( hitE );
384 if( hitE == 0 )
edm::LogWarning(
"UnexpectedEventContents")<<
"HCal hitE==0? (or unknown subdetector?)";
390 int EcalNum = cellIDs[iCell].subdetId();
393 EBDetId EcalID = cellIDs[iCell];
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];
403 <<
"Can't find the EE recHit with ID: "<<EcalID;
continue;}
404 hitE = theRecHit->energy();
405 EE_energies.push_back( hitE );
407 if( hitE == 0 )
edm::LogWarning(
"UnexpectedEventContents")<<
"ECal hitE==0? (or unknown subdetector?)";
408 if( iDbg > 6 )
cout<<
"EcalNum: "<<EcalNum<<
" hitE: "<<hitE<<endl;
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() );
435 std::inserter (HPD_energies, HPD_energies.end()),
select2nd );
439 std::inserter (RBX_energies, RBX_energies.end()),
select2nd );
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>() );
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>() );
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";
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,
487 subtowers.clear(); Ecal_subtowers.clear(); Hcal_subtowers.clear(); HO_subtowers.clear();
488 HPD_energies.clear(); RBX_energies.clear();
490 std::map< int, double > HPD_energy_map, RBX_energy_map;
494 if( iDbg > 9 )
cout<<
"classifyJetTowers started. # of towers found: "<<nTowers<<endl;
496 for(
int iTower = 0; iTower <
nTowers ; iTower++ ) {
500 int nEM = 0, nHad = 0, nHO = 0;
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;
506 for(
int iCell = 0; iCell < nCells; ++iCell ) {
522 if( E_em != 0 ) Ecal_subtowers.push_back(
subtower( E_em, nEM ) );
525 if( E_HO != 0 ) HO_subtowers.push_back(
subtower( E_HO, nHO ) );
529 Hcal_subtowers.push_back(
subtower( E_had, nHad ) );
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;
538 if( reg == HEneg || reg == HBneg || reg == HBpos || reg == HEpos ) {
539 int oddnessIEta = HBHE_oddness( iEta );
540 if( oddnessIEta < 0 )
break;
542 int iHPD = 100 * reg;
543 int iRBX = 100 * reg + ((iPhi + 1) % 72) / 4;
545 if(( reg == HEneg || reg == HEpos ) &&
std::abs( iEta ) >= 21 ) {
546 if( (0
x1 & iPhi) == 0 ) {
548 "Bug?! Jet ID code assumes no even iPhi recHits at HE edges";
551 bool boolOddnessIEta = oddnessIEta;
552 bool upperIPhi = (( iPhi%4 ) == 1 || ( iPhi%4 ) == 2);
556 if( upperIPhi != boolOddnessIEta ) ++iPhi;
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];
576 std::inserter (HPD_energies, HPD_energies.end()),
select2nd );
580 std::inserter (RBX_energies, RBX_energies.end()),
select2nd );
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() );
596 if( ae == 29 && depth == 1 ) ae += 1;
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;
612 if( ae == 29 )
return -1;
618 if( iEta == 16 || iEta == -16 )
return unknown_region;
619 if( iEta == 29 || iEta == -29 )
return unknown_region;
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;
float hadEnergyInHE() const
float emEnergyInEE() const
T getParameter(std::string const &) const
void setComment(std::string const &value)
bool subtower_has_greater_E(reco::helper::JetIDHelper::subtower i, reco::helper::JetIDHelper::subtower j)
Jets made from CaloTowers.
size_t constituentsSize() const
HcalSubdetector subdet() const
get the subdetector
std::vector< HORecHit >::const_iterator const_iterator
float emEnergyInHF() const
std::vector< Variable::Flags > flags
double pt() const final
transverse momentum
unsigned int nCarrying(double fraction, const std::vector< double > &descending_energies)
virtual std::vector< CaloTowerPtr > getCaloConstituents() const
get all constituents
bool hasNonPositiveE(reco::helper::JetIDHelper::subtower x)
float hadEnergyInHO() const
int depth() const
get the tower depth
Container::value_type value_type
double et() const final
transverse energy
float emEnergyInEB() const
ParameterDescriptionNode * ifValue(ParameterDescription< T > const &switchParameter, std::unique_ptr< ParameterDescriptionCases< T > > cases)
int ieta() const
get the cell ieta
double energy() const final
energy
Abs< T >::type abs(const T &t)
static double select2nd(std::map< int, double >::value_type const &pair)
void fillDescription(edm::ParameterSetDescription &iDesc)
const std::vector< DetId > & constituents() const
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)
const_iterator end() const
Region HBHE_region(int iEta, int depth)
int HBHE_oddness(int iEta, int depth)
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)
iterator find(key_type k)
float hadEnergyInHB() const
float hadEnergyInHF() const
double outerEnergy() const
void calculate(const edm::Event &event, const reco::CaloJet &jet, const int iDbg=0)
unsigned int hitsInNCarrying(double fraction, const std::vector< subtower > &descending_towers)