Back to home page

Project CMSSW displayed by LXR
 
 

    

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 // JetIDHelper needs a much more detailed description that the one in HcalTopology,
0013 // so to be consistent, all needed constants are hardwired in JetIDHelper.cc itself
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 // defining stuff useful only for this class (not needed in header)
0025 namespace reco {
0026 
0027   namespace helper {
0028 
0029     // select2nd exists only in some std and boost implementations, so let's control our own fate
0030     // and it can't be a non-static member function.
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   }  // namespace helper
0040 }  // namespace reco
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   // 1) jet ID variable derived from existing fractions
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   // 2) tower based variables
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   // counts
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();  // ignores negative energies from HF!
0149 
0150   // energy fractions
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   // 3) cell based variables
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     // counts
0188     n90Hits_ = nCarrying(0.9, energies);
0189 
0190     // energy fractions
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   // slightly odd loop structure avoids round-off problems when runningE never catches up with totalE
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   // slightly odd loop structure avoids round-off problems when runningE never catches up with totalE
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   // the jet only contains DetIds, so first read recHit collection
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();  // cell == recHit
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 {  // HBHE
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;  // Depth 3 at the HE forward edge is split over tower 28 & 29, and jet reco. assigns half each
0385 
0386           int iHPD = 100 * region;
0387           int iRBX = 100 * region + ((hitIPhi + 1) % 72) / 4;  // 71,72,1,2 are in the same RBX module
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);  // the upper iPhi indices in the HE wedge
0396             // remap the iPhi so it fits the one in the inner HE regions, change in needed in two cases:
0397             // 1) in the upper iPhis of the module, the even iEtas belong to the higher iPhi
0398             // 2) in the loewr iPhis of the module, the odd  iEtas belong to the higher iPhi
0399             if (upperIPhi != oddnessIEta)
0400               ++hitIPhi;
0401             // note that hitIPhi could not be 72 before, so it's still in the legal range [1,72]
0402           }
0403           iHPD += hitIPhi;
0404 
0405           // book the energies
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         }  // if HBHE
0419         if (hitE == 0)
0420           edm::LogWarning("UnexpectedEventContents") << "HCal hitE==0? (or unknown subdetector?)";
0421 
0422       }  // if HCAL
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     }    // loop on cells
0452   }      // loop on towers
0453 
0454   /* Disabling check until HO is accounted for in EMF. Check was used in CMSSW_2, where HE was excluded.
0455   double expHcalE = jet.energy() * (1-jet.emEnergyFraction());
0456   if( totHcalE + expHcalE > 0 && 
0457       TMath::Abs( totHcalE - expHcalE ) > 0.01 && 
0458       ( totHcalE - expHcalE ) / ( totHcalE + expHcalE ) > 0.0001 ) {
0459     edm::LogWarning("CodeAssumptionsViolated")<<"failed to account for all Hcal energies"
0460                           <<totHcalE<<"!="<<expHcalE;
0461                           } */
0462   // concatenate Hcal and Ecal lists
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   // sort the energies
0471   // std::sort( Hcal_energies.begin(), Hcal_energies.end(), greater<double>() );
0472   // std::sort( Ecal_energies.begin(), Ecal_energies.end(), greater<double>() );
0473 
0474   // put the energy sums (the 2nd entry in each pair of the maps) into the output vectors and sort them
0475   std::transform(
0476       HPD_energy_map.begin(), HPD_energy_map.end(), std::inserter(HPD_energies, HPD_energies.end()), select2nd);
0477   //          std::select2nd<std::map<int,double>::value_type>());
0478   // std::sort( HPD_energies.begin(), HPD_energies.end(), greater<double>() );
0479   std::transform(
0480       RBX_energy_map.begin(), RBX_energy_map.end(), std::inserter(RBX_energies, RBX_energies.end()), select2nd);
0481   //          std::select2nd<std::map<int,double>::value_type>());
0482   // std::sort( RBX_energies.begin(), RBX_energies.end(), greater<double>() );
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   // prepare sub detector energies, then turn them into fractions
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();  // cell == recHit
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       // totHcalE += E_had;
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;  // can't assign this tower to a single readout component
0590 
0591         int iHPD = 100 * reg;
0592         int iRBX = 100 * reg + ((iPhi + 1) % 72) / 4;  // 71,72,1,2 are in the same RBX module
0593 
0594         if ((reg == HEneg || reg == HEpos) && std::abs(iEta) >= 21) {  // at low-granularity edge of HE
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);  // the upper iPhi indices in the HE wedge
0601           // remap the iPhi so it fits the one in the inner HE regions, change in needed in two cases:
0602           // 1) in the upper iPhis of the module, the even IEtas belong to the higher iPhi
0603           // 2) in the loewr iPhis of the module, the odd  IEtas belong to the higher iPhi
0604           if (upperIPhi != boolOddnessIEta)
0605             ++iPhi;
0606           // note that iPhi could not be 72 before, so it's still in the legal range [1,72]
0607         }  // if at low-granularity edge of HE
0608         iHPD += iPhi;
0609 
0610         // book the energies
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       }  // HBHE
0616     }    // E_had > 0
0617   }      // loop on towers
0618 
0619   // sort the subtowers
0620   // std::sort( Hcal_subtowers.begin(), Hcal_subtowers.end(), subtower_has_greater_E );
0621   // std::sort( Ecal_subtowers.begin(), Ecal_subtowers.end(), subtower_has_greater_E );
0622 
0623   // put the energy sums (the 2nd entry in each pair of the maps) into the output vectors and sort them
0624   std::transform(
0625       HPD_energy_map.begin(), HPD_energy_map.end(), std::inserter(HPD_energies, HPD_energies.end()), select2nd);
0626   //          std::select2nd<std::map<int,double>::value_type>());
0627   // std::sort( HPD_energies.begin(), HPD_energies.end(), greater<double>() );
0628   std::transform(
0629       RBX_energy_map.begin(), RBX_energy_map.end(), std::inserter(RBX_energies, RBX_energies.end()), select2nd);
0630   //          std::select2nd<std::map<int,double>::value_type>());
0631   // std::sort( RBX_energies.begin(), RBX_energies.end(), greater<double>() );
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 // private helper functions to figure out detector & readout geometry as needed
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;  // observed that: hits are at depths 1 & 2; 1 goes with the even pattern
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;  // can't figure it out without RecHits
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;  // both HB and HE cells belong to these towers
0672   if (iEta == 29 || iEta == -29)
0673     return unknown_region;  // both HE and HF cells belong to these towers
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 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!