File indexing completed on 2023-03-17 11:18:35
0001 #include "RecoJets/JetProducers/interface/JetIDHelper.h"
0002
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004
0005 #include "FWCore/Framework/interface/Frameworkfwd.h"
0006 #include "FWCore/Framework/interface/Event.h"
0007 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0008 #include "FWCore/Utilities/interface/InputTag.h"
0009
0010 #include "DataFormats/METReco/interface/HcalCaloFlagLabels.h"
0011 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0012
0013
0014 #include "Geometry/Records/interface/HcalRecNumberingRecord.h"
0015 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0016
0017 #include "TMath.h"
0018 #include <vector>
0019 #include <numeric>
0020 #include <iostream>
0021
0022 using namespace std;
0023
0024
0025 namespace reco {
0026
0027 namespace helper {
0028
0029
0030
0031 static double select2nd(std::map<int, double>::value_type const &pair) { return pair.second; }
0032
0033 bool hasNonPositiveE(reco::helper::JetIDHelper::subtower x) { return x.E <= 0; }
0034
0035 bool subtower_has_greater_E(reco::helper::JetIDHelper::subtower i, reco::helper::JetIDHelper::subtower j) {
0036 return i.E > j.E;
0037 }
0038 std::atomic<int> JetIDHelper::sanity_checks_left_{100};
0039 }
0040 }
0041
0042 reco::helper::JetIDHelper::JetIDHelper(edm::ParameterSet const &pset, edm::ConsumesCollector &&iC) {
0043 useRecHits_ = pset.getParameter<bool>("useRecHits");
0044 if (useRecHits_) {
0045 hbheRecHitsColl_ = pset.getParameter<edm::InputTag>("hbheRecHitsColl");
0046 hoRecHitsColl_ = pset.getParameter<edm::InputTag>("hoRecHitsColl");
0047 hfRecHitsColl_ = pset.getParameter<edm::InputTag>("hfRecHitsColl");
0048 ebRecHitsColl_ = pset.getParameter<edm::InputTag>("ebRecHitsColl");
0049 eeRecHitsColl_ = pset.getParameter<edm::InputTag>("eeRecHitsColl");
0050 }
0051 initValues();
0052
0053 input_HBHERecHits_token_ = iC.consumes<HBHERecHitCollection>(hbheRecHitsColl_);
0054 input_HORecHits_token_ = iC.consumes<HORecHitCollection>(hoRecHitsColl_);
0055 input_HFRecHits_token_ = iC.consumes<HFRecHitCollection>(hfRecHitsColl_);
0056 input_EBRecHits_token_ = iC.consumes<EBRecHitCollection>(ebRecHitsColl_);
0057 input_EERecHits_token_ = iC.consumes<EERecHitCollection>(eeRecHitsColl_);
0058
0059 hcal_topo_token_ = iC.esConsumes();
0060 }
0061
0062 void reco::helper::JetIDHelper::initValues() {
0063 fHPD_ = -1.0;
0064 fRBX_ = -1.0;
0065 n90Hits_ = -1;
0066 fSubDetector1_ = -1.0;
0067 fSubDetector2_ = -1.0;
0068 fSubDetector3_ = -1.0;
0069 fSubDetector4_ = -1.0;
0070 restrictedEMF_ = -1.0;
0071 nHCALTowers_ = -1;
0072 nECALTowers_ = -1;
0073 approximatefHPD_ = -1.0;
0074 approximatefRBX_ = -1.0;
0075 hitsInN90_ = -1;
0076 fEB_ = fEE_ = fHB_ = fHE_ = fHO_ = fLong_ = fShort_ = -1.0;
0077 fLS_ = fHFOOT_ = -1.0;
0078 }
0079
0080 void reco::helper::JetIDHelper::fillDescription(edm::ParameterSetDescription &iDesc) {
0081 iDesc
0082 .ifValue(edm::ParameterDescription<bool>("useRecHits", true, true),
0083 true >> (edm::ParameterDescription<edm::InputTag>("hbheRecHitsColl", edm::InputTag(), true) and
0084 edm::ParameterDescription<edm::InputTag>("hoRecHitsColl", edm::InputTag(), true) and
0085 edm::ParameterDescription<edm::InputTag>("hfRecHitsColl", edm::InputTag(), true) and
0086 edm::ParameterDescription<edm::InputTag>("ebRecHitsColl", edm::InputTag(), true) and
0087 edm::ParameterDescription<edm::InputTag>("eeRecHitsColl", edm::InputTag(), true)))
0088 ->setComment(
0089 "If using RecHits to calculate the precise jet ID variables that need them, "
0090 "their sources need to be specified");
0091 }
0092
0093 void reco::helper::JetIDHelper::calculate(const edm::Event &event,
0094 const edm::EventSetup &setup,
0095 const reco::CaloJet &jet,
0096 const int iDbg) {
0097 initValues();
0098
0099
0100
0101
0102
0103 double E_EM = TMath::Max(float(0.), jet.emEnergyInHF()) + TMath::Max(float(0.), jet.emEnergyInEB()) +
0104 TMath::Max(float(0.), jet.emEnergyInEE());
0105 double E_Had = TMath::Max(float(0.), jet.hadEnergyInHB()) + TMath::Max(float(0.), jet.hadEnergyInHE()) +
0106 TMath::Max(float(0.), jet.hadEnergyInHO()) + TMath::Max(float(0.), jet.hadEnergyInHF());
0107 if (E_Had + E_EM > 0)
0108 restrictedEMF_ = E_EM / (E_EM + E_Had);
0109 if (iDbg > 1)
0110 cout << "jet pT: " << jet.pt() << ", eT: " << jet.et() << ", E: " << jet.energy() << " rEMF: " << restrictedEMF_
0111 << endl;
0112
0113
0114
0115
0116 vector<subtower> subtowers, Ecal_subtowers, Hcal_subtowers, HO_subtowers;
0117 vector<double> HPD_energies, RBX_energies;
0118
0119 classifyJetTowers(
0120 event, jet, subtowers, Ecal_subtowers, Hcal_subtowers, HO_subtowers, HPD_energies, RBX_energies, iDbg);
0121 if (iDbg > 1) {
0122 cout << "E:";
0123 for (unsigned int i = 0; i < subtowers.size(); ++i)
0124 cout << " " << subtowers[i].E << "," << subtowers[i].Nhit;
0125 cout << "\nECal_E:";
0126 for (unsigned int i = 0; i < Ecal_subtowers.size(); ++i)
0127 cout << " " << Ecal_subtowers[i].E << "," << Ecal_subtowers[i].Nhit;
0128 cout << "\nHCal_E:";
0129 for (unsigned int i = 0; i < Hcal_subtowers.size(); ++i)
0130 cout << " " << Hcal_subtowers[i].E << "," << Hcal_subtowers[i].Nhit;
0131 cout << "\nHO_E:";
0132 for (unsigned int i = 0; i < HO_subtowers.size(); ++i)
0133 cout << " " << HO_subtowers[i].E << "," << HO_subtowers[i].Nhit;
0134 cout << "\nHPD_E:";
0135 for (unsigned int i = 0; i < HPD_energies.size(); ++i)
0136 cout << " " << HPD_energies[i];
0137 cout << "\nRBX_E:";
0138 for (unsigned int i = 0; i < RBX_energies.size(); ++i)
0139 cout << " " << RBX_energies[i];
0140 cout << endl;
0141 }
0142
0143
0144 hitsInN90_ = hitsInNCarrying(0.9, subtowers);
0145 nHCALTowers_ = Hcal_subtowers.size();
0146 vector<subtower>::const_iterator it;
0147 it = find_if(Ecal_subtowers.begin(), Ecal_subtowers.end(), hasNonPositiveE);
0148 nECALTowers_ = it - Ecal_subtowers.begin();
0149
0150
0151 double max_HPD_energy = 0., max_RBX_energy = 0.;
0152 std::vector<double>::const_iterator it_max_HPD_energy = std::max_element(HPD_energies.begin(), HPD_energies.end());
0153 std::vector<double>::const_iterator it_max_RBX_energy = std::max_element(RBX_energies.begin(), RBX_energies.end());
0154 if (it_max_HPD_energy != HPD_energies.end()) {
0155 max_HPD_energy = *it_max_HPD_energy;
0156 }
0157 if (it_max_RBX_energy != RBX_energies.end()) {
0158 max_RBX_energy = *it_max_RBX_energy;
0159 }
0160 if (jet.energy() > 0) {
0161 if (!HPD_energies.empty())
0162 approximatefHPD_ = max_HPD_energy / jet.energy();
0163 if (!RBX_energies.empty())
0164 approximatefRBX_ = max_HPD_energy / jet.energy();
0165 }
0166
0167
0168
0169
0170 if (useRecHits_) {
0171 vector<double> energies, subdet_energies, Ecal_energies, Hcal_energies, HO_energies;
0172 double LS_bad_energy, HF_OOT_energy;
0173 classifyJetComponents(event,
0174 setup,
0175 jet,
0176 energies,
0177 subdet_energies,
0178 Ecal_energies,
0179 Hcal_energies,
0180 HO_energies,
0181 HPD_energies,
0182 RBX_energies,
0183 LS_bad_energy,
0184 HF_OOT_energy,
0185 iDbg);
0186
0187
0188 n90Hits_ = nCarrying(0.9, energies);
0189
0190
0191 if (jet.energy() > 0) {
0192 if (!HPD_energies.empty())
0193 fHPD_ = max_HPD_energy / jet.energy();
0194 if (!RBX_energies.empty())
0195 fRBX_ = max_RBX_energy / jet.energy();
0196 if (!subdet_energies.empty())
0197 fSubDetector1_ = subdet_energies.at(0) / jet.energy();
0198 if (subdet_energies.size() > 1)
0199 fSubDetector2_ = subdet_energies.at(1) / jet.energy();
0200 if (subdet_energies.size() > 2)
0201 fSubDetector3_ = subdet_energies.at(2) / jet.energy();
0202 if (subdet_energies.size() > 3)
0203 fSubDetector4_ = subdet_energies.at(3) / jet.energy();
0204 fLS_ = LS_bad_energy / jet.energy();
0205 fHFOOT_ = HF_OOT_energy / jet.energy();
0206
0207 if (iDbg > 0 && sanity_checks_left_.load(std::memory_order_acquire) > 0) {
0208 --sanity_checks_left_;
0209 double EH_sum = accumulate(Ecal_energies.begin(), Ecal_energies.end(), 0.);
0210 EH_sum = accumulate(Hcal_energies.begin(), Hcal_energies.end(), EH_sum);
0211 double EHO_sum = accumulate(HO_energies.begin(), HO_energies.end(), EH_sum);
0212 if (jet.energy() > 0.001 && abs(EH_sum / jet.energy() - 1) > 0.01 && abs(EHO_sum / jet.energy() - 1) > 0.01)
0213 edm::LogWarning("BadInput") << "Jet energy (" << jet.energy()
0214 << ") does not match the total energy in the recHits (" << EHO_sum
0215 << ", or without HO: " << EH_sum << ") . Are these the right recHits? "
0216 << "Were jet energy corrections mistakenly applied before jet ID? A bug?";
0217 if (iDbg > 1)
0218 cout << "Sanity check - E: " << jet.energy() << " =? EH_sum: " << EH_sum << " / EHO_sum: " << EHO_sum << endl;
0219 }
0220 }
0221
0222 if (iDbg > 1) {
0223 cout << "DBG - fHPD: " << fHPD_ << ", fRBX: " << fRBX_ << ", nh90: " << n90Hits_ << ", fLS: " << fLS_
0224 << ", fHFOOT: " << fHFOOT_ << endl;
0225 cout << " -~fHPD: " << approximatefHPD_ << ", ~fRBX: " << approximatefRBX_ << ", hits in n90: " << hitsInN90_
0226 << endl;
0227 cout << " - nHCALTowers: " << nHCALTowers_ << ", nECALTowers: " << nECALTowers_
0228 << "; subdet fractions: " << fSubDetector1_ << ", " << fSubDetector2_ << ", " << fSubDetector3_ << ", "
0229 << fSubDetector4_ << endl;
0230 }
0231 }
0232 }
0233
0234 unsigned int reco::helper::JetIDHelper::nCarrying(double fraction, const std::vector<double> &descending_energies) {
0235 double totalE = 0;
0236 for (unsigned int i = 0; i < descending_energies.size(); ++i)
0237 totalE += descending_energies[i];
0238
0239 double runningE = 0;
0240 unsigned int NC = descending_energies.size();
0241
0242
0243 for (unsigned int i = descending_energies.size(); i > 0; --i) {
0244 runningE += descending_energies[i - 1];
0245 if (runningE < (1 - fraction) * totalE)
0246 NC = i - 1;
0247 }
0248 return NC;
0249 }
0250
0251 unsigned int reco::helper::JetIDHelper::hitsInNCarrying(double fraction,
0252 const std::vector<subtower> &descending_towers) {
0253 double totalE = 0;
0254 for (unsigned int i = 0; i < descending_towers.size(); ++i)
0255 totalE += descending_towers[i].E;
0256
0257 double runningE = 0;
0258 unsigned int NH = 0;
0259
0260
0261 for (unsigned int i = descending_towers.size(); i > 0; --i) {
0262 runningE += descending_towers[i - 1].E;
0263 if (runningE >= (1 - fraction) * totalE)
0264 NH += descending_towers[i - 1].Nhit;
0265 }
0266 return NH;
0267 }
0268
0269 void reco::helper::JetIDHelper::classifyJetComponents(const edm::Event &event,
0270 const edm::EventSetup &setup,
0271 const reco::CaloJet &jet,
0272 vector<double> &energies,
0273 vector<double> &subdet_energies,
0274 vector<double> &Ecal_energies,
0275 vector<double> &Hcal_energies,
0276 vector<double> &HO_energies,
0277 vector<double> &HPD_energies,
0278 vector<double> &RBX_energies,
0279 double &LS_bad_energy,
0280 double &HF_OOT_energy,
0281 const int iDbg) {
0282 energies.clear();
0283 subdet_energies.clear();
0284 Ecal_energies.clear();
0285 Hcal_energies.clear();
0286 HO_energies.clear();
0287 HPD_energies.clear();
0288 RBX_energies.clear();
0289 LS_bad_energy = HF_OOT_energy = 0.;
0290
0291 edm::ESHandle<HcalTopology> theHcalTopology = setup.getHandle(hcal_topo_token_);
0292
0293 std::map<int, double> HPD_energy_map, RBX_energy_map;
0294 vector<double> EB_energies, EE_energies, HB_energies, HE_energies, short_energies, long_energies;
0295 edm::Handle<HBHERecHitCollection> HBHERecHits;
0296 edm::Handle<HORecHitCollection> HORecHits;
0297 edm::Handle<HFRecHitCollection> HFRecHits;
0298 edm::Handle<EBRecHitCollection> EBRecHits;
0299 edm::Handle<EERecHitCollection> EERecHits;
0300
0301 event.getByToken(input_HBHERecHits_token_, HBHERecHits);
0302 event.getByToken(input_HORecHits_token_, HORecHits);
0303 event.getByToken(input_HFRecHits_token_, HFRecHits);
0304 event.getByToken(input_EBRecHits_token_, EBRecHits);
0305 event.getByToken(input_EERecHits_token_, EERecHits);
0306 if (iDbg > 2)
0307 cout << "# of rechits found - HBHE: " << HBHERecHits->size() << ", HO: " << HORecHits->size()
0308 << ", HF: " << HFRecHits->size() << ", EB: " << EBRecHits->size() << ", EE: " << EERecHits->size() << endl;
0309
0310 vector<CaloTowerPtr> towers = jet.getCaloConstituents();
0311 int nTowers = towers.size();
0312 if (iDbg > 9)
0313 cout << "In classifyJetComponents. # of towers found: " << nTowers << endl;
0314
0315 for (int iTower = 0; iTower < nTowers; iTower++) {
0316 CaloTowerPtr &tower = towers[iTower];
0317
0318 int nCells = tower->constituentsSize();
0319 if (iDbg)
0320 cout << "tower #" << iTower << " has " << nCells << " cells. "
0321 << "It's at iEta: " << tower->ieta() << ", iPhi: " << tower->iphi() << endl;
0322
0323 const vector<DetId> &cellIDs = tower->constituents();
0324
0325 for (int iCell = 0; iCell < nCells; ++iCell) {
0326 DetId::Detector detNum = cellIDs[iCell].det();
0327 if (detNum == DetId::Hcal) {
0328 HcalDetId HcalID = cellIDs[iCell];
0329 HcalSubdetector HcalNum = HcalID.subdet();
0330 double hitE = 0;
0331 if (HcalNum == HcalOuter) {
0332 HORecHitCollection::const_iterator theRecHit = HORecHits->find(HcalID);
0333 if (theRecHit == HORecHits->end()) {
0334 edm::LogWarning("UnexpectedEventContents") << "Can't find the HO recHit with ID: " << HcalID;
0335 continue;
0336 }
0337 hitE = theRecHit->energy();
0338 HO_energies.push_back(hitE);
0339
0340 } else if (HcalNum == HcalForward) {
0341 HFRecHitCollection::const_iterator theRecHit = HFRecHits->find(HcalID);
0342 if (theRecHit == HFRecHits->end()) {
0343 edm::LogWarning("UnexpectedEventContents") << "Can't find the HF recHit with ID: " << HcalID;
0344 continue;
0345 }
0346 hitE = theRecHit->energy();
0347 if (iDbg > 4)
0348 cout << "hit #" << iCell << " is HF , E: " << hitE << " iEta: " << theRecHit->id().ieta()
0349 << ", depth: " << theRecHit->id().depth() << ", iPhi: " << theRecHit->id().iphi();
0350
0351 if (HcalID.depth() == 1)
0352 long_energies.push_back(hitE);
0353 else
0354 short_energies.push_back(hitE);
0355
0356 uint32_t flags = theRecHit->flags();
0357 if (flags & (1 << HcalCaloFlagLabels::HFLongShort))
0358 LS_bad_energy += hitE;
0359 if (flags & ((1 << HcalCaloFlagLabels::HFDigiTime) | (3 << HcalCaloFlagLabels::HFTimingTrustBits) |
0360 (1 << HcalCaloFlagLabels::TimingSubtractedBit) | (1 << HcalCaloFlagLabels::TimingAddedBit) |
0361 (1 << HcalCaloFlagLabels::TimingErrorBit)))
0362 HF_OOT_energy += hitE;
0363 if (iDbg > 4 && flags)
0364 cout << "flags: " << flags << " -> LS_bad_energy: " << LS_bad_energy << ", HF_OOT_energy: " << HF_OOT_energy
0365 << endl;
0366
0367 } else {
0368
0369 HBHERecHitCollection::const_iterator theRecHit = HBHERecHits->find(HcalID);
0370 if (theRecHit == HBHERecHits->end()) {
0371 edm::LogWarning("UnexpectedEventContents") << "Can't find the HBHE recHit with ID: " << HcalID;
0372 continue;
0373 }
0374 hitE = theRecHit->energy();
0375 int iEta = theRecHit->id().ieta();
0376 int depth = theRecHit->id().depth();
0377 Region region = HBHE_region(theRecHit->id().rawId());
0378 int hitIPhi = theRecHit->id().iphi();
0379 if (iDbg > 3)
0380 cout << "hit #" << iCell << " is HBHE, E: " << hitE << " iEta: " << iEta << ", depth: " << depth
0381 << ", iPhi: " << theRecHit->id().iphi() << " -> " << region;
0382
0383 if (theHcalTopology->mergedDepth29(theRecHit->id()))
0384 hitE /= 2;
0385
0386 int iHPD = 100 * region;
0387 int iRBX = 100 * region + ((hitIPhi + 1) % 72) / 4;
0388
0389 if (std::abs(iEta) >= 21) {
0390 if ((0x1 & hitIPhi) == 0) {
0391 edm::LogError("CodeAssumptionsViolated") << "Bug?! Jet ID code assumes no even iPhi recHits at HE edges";
0392 return;
0393 }
0394 bool oddnessIEta = HBHE_oddness(iEta, depth);
0395 bool upperIPhi = ((hitIPhi % 4) == 1 || (hitIPhi % 4) == 2);
0396
0397
0398
0399 if (upperIPhi != oddnessIEta)
0400 ++hitIPhi;
0401
0402 }
0403 iHPD += hitIPhi;
0404
0405
0406 HPD_energy_map[iHPD] += hitE;
0407 RBX_energy_map[iRBX] += hitE;
0408 if (iDbg > 5)
0409 cout << " --> H[" << iHPD << "]=" << HPD_energy_map[iHPD] << ", R[" << iRBX << "]=" << RBX_energy_map[iRBX];
0410 if (iDbg > 1)
0411 cout << endl;
0412
0413 if (region == HBneg || region == HBpos)
0414 HB_energies.push_back(hitE);
0415 else
0416 HE_energies.push_back(hitE);
0417
0418 }
0419 if (hitE == 0)
0420 edm::LogWarning("UnexpectedEventContents") << "HCal hitE==0? (or unknown subdetector?)";
0421
0422 }
0423
0424 else if (detNum == DetId::Ecal) {
0425 int EcalNum = cellIDs[iCell].subdetId();
0426 double hitE = 0;
0427 if (EcalNum == 1) {
0428 EBDetId EcalID = cellIDs[iCell];
0429 EBRecHitCollection::const_iterator theRecHit = EBRecHits->find(EcalID);
0430 if (theRecHit == EBRecHits->end()) {
0431 edm::LogWarning("UnexpectedEventContents") << "Can't find the EB recHit with ID: " << EcalID;
0432 continue;
0433 }
0434 hitE = theRecHit->energy();
0435 EB_energies.push_back(hitE);
0436 } else if (EcalNum == 2) {
0437 EEDetId EcalID = cellIDs[iCell];
0438 EERecHitCollection::const_iterator theRecHit = EERecHits->find(EcalID);
0439 if (theRecHit == EERecHits->end()) {
0440 edm::LogWarning("UnexpectedEventContents") << "Can't find the EE recHit with ID: " << EcalID;
0441 continue;
0442 }
0443 hitE = theRecHit->energy();
0444 EE_energies.push_back(hitE);
0445 }
0446 if (hitE == 0)
0447 edm::LogWarning("UnexpectedEventContents") << "ECal hitE==0? (or unknown subdetector?)";
0448 if (iDbg > 6)
0449 cout << "EcalNum: " << EcalNum << " hitE: " << hitE << endl;
0450 }
0451 }
0452 }
0453
0454
0455
0456
0457
0458
0459
0460
0461
0462
0463 Hcal_energies.insert(Hcal_energies.end(), HB_energies.begin(), HB_energies.end());
0464 Hcal_energies.insert(Hcal_energies.end(), HE_energies.begin(), HE_energies.end());
0465 Hcal_energies.insert(Hcal_energies.end(), short_energies.begin(), short_energies.end());
0466 Hcal_energies.insert(Hcal_energies.end(), long_energies.begin(), long_energies.end());
0467 Ecal_energies.insert(Ecal_energies.end(), EB_energies.begin(), EB_energies.end());
0468 Ecal_energies.insert(Ecal_energies.end(), EE_energies.begin(), EE_energies.end());
0469
0470
0471
0472
0473
0474
0475 std::transform(
0476 HPD_energy_map.begin(), HPD_energy_map.end(), std::inserter(HPD_energies, HPD_energies.end()), select2nd);
0477
0478
0479 std::transform(
0480 RBX_energy_map.begin(), RBX_energy_map.end(), std::inserter(RBX_energies, RBX_energies.end()), select2nd);
0481
0482
0483
0484 energies.insert(energies.end(), Hcal_energies.begin(), Hcal_energies.end());
0485 energies.insert(energies.end(), Ecal_energies.begin(), Ecal_energies.end());
0486 energies.insert(energies.end(), HO_energies.begin(), HO_energies.end());
0487 std::sort(energies.begin(), energies.end(), greater<double>());
0488
0489
0490 fEB_ = std::accumulate(EB_energies.begin(), EB_energies.end(), 0.);
0491 fEE_ = std::accumulate(EE_energies.begin(), EE_energies.end(), 0.);
0492 fHB_ = std::accumulate(HB_energies.begin(), HB_energies.end(), 0.);
0493 fHE_ = std::accumulate(HE_energies.begin(), HE_energies.end(), 0.);
0494 fHO_ = std::accumulate(HO_energies.begin(), HO_energies.end(), 0.);
0495 fShort_ = std::accumulate(short_energies.begin(), short_energies.end(), 0.);
0496 fLong_ = std::accumulate(long_energies.begin(), long_energies.end(), 0.);
0497 subdet_energies.push_back(fEB_);
0498 subdet_energies.push_back(fEE_);
0499 subdet_energies.push_back(fHB_);
0500 subdet_energies.push_back(fHE_);
0501 subdet_energies.push_back(fHO_);
0502 subdet_energies.push_back(fShort_);
0503 subdet_energies.push_back(fLong_);
0504 std::sort(subdet_energies.begin(), subdet_energies.end(), greater<double>());
0505 if (jet.energy() > 0) {
0506 fEB_ /= jet.energy();
0507 fEE_ /= jet.energy();
0508 fHB_ /= jet.energy();
0509 fHE_ /= jet.energy();
0510 fHO_ /= jet.energy();
0511 fShort_ /= jet.energy();
0512 fLong_ /= jet.energy();
0513 } else {
0514 if (fEB_ > 0 || fEE_ > 0 || fHB_ > 0 || fHE_ > 0 || fHO_ > 0 || fShort_ > 0 || fLong_ > 0)
0515 edm::LogError("UnexpectedEventContents") << "Jet ID Helper found energy in subdetectors and jet E <= 0";
0516 }
0517 }
0518
0519 void reco::helper::JetIDHelper::classifyJetTowers(const edm::Event &event,
0520 const reco::CaloJet &jet,
0521 vector<subtower> &subtowers,
0522 vector<subtower> &Ecal_subtowers,
0523 vector<subtower> &Hcal_subtowers,
0524 vector<subtower> &HO_subtowers,
0525 vector<double> &HPD_energies,
0526 vector<double> &RBX_energies,
0527 const int iDbg) {
0528 subtowers.clear();
0529 Ecal_subtowers.clear();
0530 Hcal_subtowers.clear();
0531 HO_subtowers.clear();
0532 HPD_energies.clear();
0533 RBX_energies.clear();
0534
0535 std::map<int, double> HPD_energy_map, RBX_energy_map;
0536
0537 vector<CaloTowerPtr> towers = jet.getCaloConstituents();
0538 int nTowers = towers.size();
0539 if (iDbg > 9)
0540 cout << "classifyJetTowers started. # of towers found: " << nTowers << endl;
0541
0542 for (int iTower = 0; iTower < nTowers; iTower++) {
0543 CaloTowerPtr &tower = towers[iTower];
0544
0545 int nEM = 0, nHad = 0, nHO = 0;
0546 const vector<DetId> &cellIDs = tower->constituents();
0547 int nCells = cellIDs.size();
0548 if (iDbg)
0549 cout << "tower #" << iTower << " has " << nCells << " cells. "
0550 << "It's at iEta: " << tower->ieta() << ", iPhi: " << tower->iphi() << endl;
0551
0552 for (int iCell = 0; iCell < nCells; ++iCell) {
0553 DetId::Detector detNum = cellIDs[iCell].det();
0554 if (detNum == DetId::Hcal) {
0555 HcalDetId HcalID = cellIDs[iCell];
0556 HcalSubdetector HcalNum = HcalID.subdet();
0557 if (HcalNum == HcalOuter) {
0558 ++nHO;
0559 } else {
0560 ++nHad;
0561 }
0562 } else if (detNum == DetId::Ecal) {
0563 ++nEM;
0564 }
0565 }
0566
0567 double E_em = tower->emEnergy();
0568 if (E_em != 0)
0569 Ecal_subtowers.push_back(subtower(E_em, nEM));
0570
0571 double E_HO = tower->outerEnergy();
0572 if (E_HO != 0)
0573 HO_subtowers.push_back(subtower(E_HO, nHO));
0574
0575 double E_had = tower->hadEnergy();
0576 if (E_had != 0) {
0577 Hcal_subtowers.push_back(subtower(E_had, nHad));
0578
0579
0580 int iEta = tower->ieta();
0581 Region reg = region(iEta);
0582 int iPhi = tower->iphi();
0583 if (iDbg > 3)
0584 cout << "tower has E_had: " << E_had << " iEta: " << iEta << ", iPhi: " << iPhi << " -> " << reg;
0585
0586 if (reg == HEneg || reg == HBneg || reg == HBpos || reg == HEpos) {
0587 int oddnessIEta = HBHE_oddness(iEta);
0588 if (oddnessIEta < 0)
0589 break;
0590
0591 int iHPD = 100 * reg;
0592 int iRBX = 100 * reg + ((iPhi + 1) % 72) / 4;
0593
0594 if ((reg == HEneg || reg == HEpos) && std::abs(iEta) >= 21) {
0595 if ((0x1 & iPhi) == 0) {
0596 edm::LogError("CodeAssumptionsViolated") << "Bug?! Jet ID code assumes no even iPhi recHits at HE edges";
0597 return;
0598 }
0599 bool boolOddnessIEta = oddnessIEta;
0600 bool upperIPhi = ((iPhi % 4) == 1 || (iPhi % 4) == 2);
0601
0602
0603
0604 if (upperIPhi != boolOddnessIEta)
0605 ++iPhi;
0606
0607 }
0608 iHPD += iPhi;
0609
0610
0611 HPD_energy_map[iHPD] += E_had;
0612 RBX_energy_map[iRBX] += E_had;
0613 if (iDbg > 5)
0614 cout << " --> H[" << iHPD << "]=" << HPD_energy_map[iHPD] << ", R[" << iRBX << "]=" << RBX_energy_map[iRBX];
0615 }
0616 }
0617 }
0618
0619
0620
0621
0622
0623
0624 std::transform(
0625 HPD_energy_map.begin(), HPD_energy_map.end(), std::inserter(HPD_energies, HPD_energies.end()), select2nd);
0626
0627
0628 std::transform(
0629 RBX_energy_map.begin(), RBX_energy_map.end(), std::inserter(RBX_energies, RBX_energies.end()), select2nd);
0630
0631
0632
0633 subtowers.insert(subtowers.end(), Hcal_subtowers.begin(), Hcal_subtowers.end());
0634 subtowers.insert(subtowers.end(), Ecal_subtowers.begin(), Ecal_subtowers.end());
0635 subtowers.insert(subtowers.end(), HO_subtowers.begin(), HO_subtowers.end());
0636 std::sort(subtowers.begin(), subtowers.end(), subtower_has_greater_E);
0637 }
0638
0639
0640
0641
0642 int reco::helper::JetIDHelper::HBHE_oddness(int iEta, int depth) {
0643 int ae = TMath::Abs(iEta);
0644 if (ae == 29 && depth == 1)
0645 ae += 1;
0646 return ae & 0x1;
0647 }
0648
0649 reco::helper::JetIDHelper::Region reco::helper::JetIDHelper::HBHE_region(uint32_t rawid) {
0650 HcalDetId id(rawid);
0651 if (id.subdet() == HcalEndcap) {
0652 if (id.ieta() < 0)
0653 return HEneg;
0654 else
0655 return HEpos;
0656 }
0657 if (id.subdet() == HcalBarrel && id.ieta() < 0)
0658 return HBneg;
0659 return HBpos;
0660 }
0661
0662 int reco::helper::JetIDHelper::HBHE_oddness(int iEta) {
0663 int ae = TMath::Abs(iEta);
0664 if (ae == 29)
0665 return -1;
0666 return ae & 0x1;
0667 }
0668
0669 reco::helper::JetIDHelper::Region reco::helper::JetIDHelper::region(int iEta) {
0670 if (iEta == 16 || iEta == -16)
0671 return unknown_region;
0672 if (iEta == 29 || iEta == -29)
0673 return unknown_region;
0674 if (iEta <= -30)
0675 return HFneg;
0676 if (iEta >= 30)
0677 return HFpos;
0678 if (iEta <= -17)
0679 return HEneg;
0680 if (iEta >= 17)
0681 return HEpos;
0682 if (iEta < 0)
0683 return HBneg;
0684 return HBpos;
0685 }