JetIDHelper.cc

Go to the documentation of this file.

#include "RecoJets/JetProducers/interface/JetIDHelper.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/EDProducer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Utilities/interface/InputTag.h"

#include "DataFormats/METReco/interface/HcalCaloFlagLabels.h"
#include "DataFormats/HcalDetId/interface/HcalDetId.h"
// JetIDHelper needs a much more detailed description that the one in HcalTopology,
// so to be consistent, all needed constants are hardwired in JetIDHelper.cc itself
#include "Geometry/Records/interface/HcalRecNumberingRecord.h"
#include "Geometry/CaloTopology/interface/HcalTopology.h"

#include "TMath.h"
#include <vector>
#include <numeric>
#include <iostream>

using namespace std;

// defining stuff useful only for this class (not needed in header)
namespace reco {

  namespace helper {

    // select2nd exists only in some std and boost implementations, so let's control our own fate
    // and it can't be a non-static member function.
    static double select2nd(std::map<int, double>::value_type const &pair) { return pair.second; }

    bool hasNonPositiveE(reco::helper::JetIDHelper::subtower x) { return x.E <= 0; }

    bool subtower_has_greater_E(reco::helper::JetIDHelper::subtower i, reco::helper::JetIDHelper::subtower j) {
      return i.E > j.E;
    }
    std::atomic<int> JetIDHelper::sanity_checks_left_{100};
  }  // namespace helper
}  // namespace reco

reco::helper::JetIDHelper::JetIDHelper(edm::ParameterSet const &pset, edm::ConsumesCollector &&iC) {
  useRecHits_ = pset.getParameter<bool>("useRecHits");
  if (useRecHits_) {
    hbheRecHitsColl_ = pset.getParameter<edm::InputTag>("hbheRecHitsColl");
    hoRecHitsColl_ = pset.getParameter<edm::InputTag>("hoRecHitsColl");
    hfRecHitsColl_ = pset.getParameter<edm::InputTag>("hfRecHitsColl");
    ebRecHitsColl_ = pset.getParameter<edm::InputTag>("ebRecHitsColl");
    eeRecHitsColl_ = pset.getParameter<edm::InputTag>("eeRecHitsColl");
  }
  initValues();

  input_HBHERecHits_token_ = iC.consumes<HBHERecHitCollection>(hbheRecHitsColl_);
  input_HORecHits_token_ = iC.consumes<HORecHitCollection>(hoRecHitsColl_);
  input_HFRecHits_token_ = iC.consumes<HFRecHitCollection>(hfRecHitsColl_);
  input_EBRecHits_token_ = iC.consumes<EBRecHitCollection>(ebRecHitsColl_);
  input_EERecHits_token_ = iC.consumes<EERecHitCollection>(eeRecHitsColl_);
}

void reco::helper::JetIDHelper::initValues() {
  fHPD_ = -1.0;
  fRBX_ = -1.0;
  n90Hits_ = -1;
  fSubDetector1_ = -1.0;
  fSubDetector2_ = -1.0;
  fSubDetector3_ = -1.0;
  fSubDetector4_ = -1.0;
  restrictedEMF_ = -1.0;
  nHCALTowers_ = -1;
  nECALTowers_ = -1;
  approximatefHPD_ = -1.0;
  approximatefRBX_ = -1.0;
  hitsInN90_ = -1;
  fEB_ = fEE_ = fHB_ = fHE_ = fHO_ = fLong_ = fShort_ = -1.0;
  fLS_ = fHFOOT_ = -1.0;
}

void reco::helper::JetIDHelper::fillDescription(edm::ParameterSetDescription &iDesc) {
  iDesc
      .ifValue(edm::ParameterDescription<bool>("useRecHits", true, true),
               true >> (edm::ParameterDescription<edm::InputTag>("hbheRecHitsColl", edm::InputTag(), true) and
                        edm::ParameterDescription<edm::InputTag>("hoRecHitsColl", edm::InputTag(), true) and
                        edm::ParameterDescription<edm::InputTag>("hfRecHitsColl", edm::InputTag(), true) and
                        edm::ParameterDescription<edm::InputTag>("ebRecHitsColl", edm::InputTag(), true) and
                        edm::ParameterDescription<edm::InputTag>("eeRecHitsColl", edm::InputTag(), true)))
      ->setComment(
          "If using RecHits to calculate the precise jet ID variables that need them, "
          "their sources need to be specified");
}

void reco::helper::JetIDHelper::calculate(const edm::Event &event,
                                          const edm::EventSetup &setup,
                                          const reco::CaloJet &jet,
                                          const int iDbg) {
  initValues();

  // --------------------------------------------------
  // 1) jet ID variable derived from existing fractions
  // --------------------------------------------------

  double E_EM = TMath::Max(float(0.), jet.emEnergyInHF()) + TMath::Max(float(0.), jet.emEnergyInEB()) +
                TMath::Max(float(0.), jet.emEnergyInEE());
  double E_Had = TMath::Max(float(0.), jet.hadEnergyInHB()) + TMath::Max(float(0.), jet.hadEnergyInHE()) +
                 TMath::Max(float(0.), jet.hadEnergyInHO()) + TMath::Max(float(0.), jet.hadEnergyInHF());
  if (E_Had + E_EM > 0)
    restrictedEMF_ = E_EM / (E_EM + E_Had);
  if (iDbg > 1)
    cout << "jet pT: " << jet.pt() << ", eT: " << jet.et() << ", E: " << jet.energy() << " rEMF: " << restrictedEMF_
         << endl;

  // ------------------------
  // 2) tower based variables
  // ------------------------
  vector<subtower> subtowers, Ecal_subtowers, Hcal_subtowers, HO_subtowers;
  vector<double> HPD_energies, RBX_energies;

  classifyJetTowers(
      event, jet, subtowers, Ecal_subtowers, Hcal_subtowers, HO_subtowers, HPD_energies, RBX_energies, iDbg);
  if (iDbg > 1) {
    cout << "E:";
    for (unsigned int i = 0; i < subtowers.size(); ++i)
      cout << " " << subtowers[i].E << "," << subtowers[i].Nhit;
    cout << "\nECal_E:";
    for (unsigned int i = 0; i < Ecal_subtowers.size(); ++i)
      cout << " " << Ecal_subtowers[i].E << "," << Ecal_subtowers[i].Nhit;
    cout << "\nHCal_E:";
    for (unsigned int i = 0; i < Hcal_subtowers.size(); ++i)
      cout << " " << Hcal_subtowers[i].E << "," << Hcal_subtowers[i].Nhit;
    cout << "\nHO_E:";
    for (unsigned int i = 0; i < HO_subtowers.size(); ++i)
      cout << " " << HO_subtowers[i].E << "," << HO_subtowers[i].Nhit;
    cout << "\nHPD_E:";
    for (unsigned int i = 0; i < HPD_energies.size(); ++i)
      cout << " " << HPD_energies[i];
    cout << "\nRBX_E:";
    for (unsigned int i = 0; i < RBX_energies.size(); ++i)
      cout << " " << RBX_energies[i];
    cout << endl;
  }

  // counts
  hitsInN90_ = hitsInNCarrying(0.9, subtowers);
  nHCALTowers_ = Hcal_subtowers.size();
  vector<subtower>::const_iterator it;
  it = find_if(Ecal_subtowers.begin(), Ecal_subtowers.end(), hasNonPositiveE);
  nECALTowers_ = it - Ecal_subtowers.begin();  // ignores negative energies from HF!

  // energy fractions
  double max_HPD_energy = 0., max_RBX_energy = 0.;
  std::vector<double>::const_iterator it_max_HPD_energy = std::max_element(HPD_energies.begin(), HPD_energies.end());
  std::vector<double>::const_iterator it_max_RBX_energy = std::max_element(RBX_energies.begin(), RBX_energies.end());
  if (it_max_HPD_energy != HPD_energies.end()) {
    max_HPD_energy = *it_max_HPD_energy;
  }
  if (it_max_RBX_energy != RBX_energies.end()) {
    max_RBX_energy = *it_max_RBX_energy;
  }
  if (jet.energy() > 0) {
    if (!HPD_energies.empty())
      approximatefHPD_ = max_HPD_energy / jet.energy();
    if (!RBX_energies.empty())
      approximatefRBX_ = max_HPD_energy / jet.energy();
  }

  // -----------------------
  // 3) cell based variables
  // -----------------------
  if (useRecHits_) {
    vector<double> energies, subdet_energies, Ecal_energies, Hcal_energies, HO_energies;
    double LS_bad_energy, HF_OOT_energy;
    classifyJetComponents(event,
                          setup,
                          jet,
                          energies,
                          subdet_energies,
                          Ecal_energies,
                          Hcal_energies,
                          HO_energies,
                          HPD_energies,
                          RBX_energies,
                          LS_bad_energy,
                          HF_OOT_energy,
                          iDbg);

    // counts
    n90Hits_ = nCarrying(0.9, energies);

    // energy fractions
    if (jet.energy() > 0) {
      if (!HPD_energies.empty())
        fHPD_ = max_HPD_energy / jet.energy();
      if (!RBX_energies.empty())
        fRBX_ = max_RBX_energy / jet.energy();
      if (!subdet_energies.empty())
        fSubDetector1_ = subdet_energies.at(0) / jet.energy();
      if (subdet_energies.size() > 1)
        fSubDetector2_ = subdet_energies.at(1) / jet.energy();
      if (subdet_energies.size() > 2)
        fSubDetector3_ = subdet_energies.at(2) / jet.energy();
      if (subdet_energies.size() > 3)
        fSubDetector4_ = subdet_energies.at(3) / jet.energy();
      fLS_ = LS_bad_energy / jet.energy();
      fHFOOT_ = HF_OOT_energy / jet.energy();

      if (iDbg > 0 && sanity_checks_left_.load(std::memory_order_acquire) > 0) {
        --sanity_checks_left_;
        double EH_sum = accumulate(Ecal_energies.begin(), Ecal_energies.end(), 0.);
        EH_sum = accumulate(Hcal_energies.begin(), Hcal_energies.end(), EH_sum);
        double EHO_sum = accumulate(HO_energies.begin(), HO_energies.end(), EH_sum);
        if (jet.energy() > 0.001 && abs(EH_sum / jet.energy() - 1) > 0.01 && abs(EHO_sum / jet.energy() - 1) > 0.01)
          edm::LogWarning("BadInput") << "Jet energy (" << jet.energy()
                                      << ") does not match the total energy in the recHits (" << EHO_sum
                                      << ", or without HO: " << EH_sum << ") . Are these the right recHits? "
                                      << "Were jet energy corrections mistakenly applied before jet ID? A bug?";
        if (iDbg > 1)
          cout << "Sanity check - E: " << jet.energy() << " =? EH_sum: " << EH_sum << " / EHO_sum: " << EHO_sum << endl;
      }
    }

    if (iDbg > 1) {
      cout << "DBG - fHPD: " << fHPD_ << ", fRBX: " << fRBX_ << ", nh90: " << n90Hits_ << ", fLS: " << fLS_
           << ", fHFOOT: " << fHFOOT_ << endl;
      cout << "    -~fHPD: " << approximatefHPD_ << ", ~fRBX: " << approximatefRBX_ << ", hits in n90: " << hitsInN90_
           << endl;
      cout << "    - nHCALTowers: " << nHCALTowers_ << ", nECALTowers: " << nECALTowers_
           << "; subdet fractions: " << fSubDetector1_ << ", " << fSubDetector2_ << ", " << fSubDetector3_ << ", "
           << fSubDetector4_ << endl;
    }
  }
}

unsigned int reco::helper::JetIDHelper::nCarrying(double fraction, const std::vector<double> &descending_energies) {
  double totalE = 0;
  for (unsigned int i = 0; i < descending_energies.size(); ++i)
    totalE += descending_energies[i];

  double runningE = 0;
  unsigned int NC = descending_energies.size();

  // slightly odd loop structure avoids round-off problems when runningE never catches up with totalE
  for (unsigned int i = descending_energies.size(); i > 0; --i) {
    runningE += descending_energies[i - 1];
    if (runningE < (1 - fraction) * totalE)
      NC = i - 1;
  }
  return NC;
}

unsigned int reco::helper::JetIDHelper::hitsInNCarrying(double fraction,
                                                        const std::vector<subtower> &descending_towers) {
  double totalE = 0;
  for (unsigned int i = 0; i < descending_towers.size(); ++i)
    totalE += descending_towers[i].E;

  double runningE = 0;
  unsigned int NH = 0;

  // slightly odd loop structure avoids round-off problems when runningE never catches up with totalE
  for (unsigned int i = descending_towers.size(); i > 0; --i) {
    runningE += descending_towers[i - 1].E;
    if (runningE >= (1 - fraction) * totalE)
      NH += descending_towers[i - 1].Nhit;
  }
  return NH;
}

void reco::helper::JetIDHelper::classifyJetComponents(const edm::Event &event,
                                                      const edm::EventSetup &setup,
                                                      const reco::CaloJet &jet,
                                                      vector<double> &energies,
                                                      vector<double> &subdet_energies,
                                                      vector<double> &Ecal_energies,
                                                      vector<double> &Hcal_energies,
                                                      vector<double> &HO_energies,
                                                      vector<double> &HPD_energies,
                                                      vector<double> &RBX_energies,
                                                      double &LS_bad_energy,
                                                      double &HF_OOT_energy,
                                                      const int iDbg) {
  energies.clear();
  subdet_energies.clear();
  Ecal_energies.clear();
  Hcal_energies.clear();
  HO_energies.clear();
  HPD_energies.clear();
  RBX_energies.clear();
  LS_bad_energy = HF_OOT_energy = 0.;

  edm::ESHandle<HcalTopology> theHcalTopology;
  setup.get<HcalRecNumberingRecord>().get(theHcalTopology);

  std::map<int, double> HPD_energy_map, RBX_energy_map;
  vector<double> EB_energies, EE_energies, HB_energies, HE_energies, short_energies, long_energies;
  edm::Handle<HBHERecHitCollection> HBHERecHits;
  edm::Handle<HORecHitCollection> HORecHits;
  edm::Handle<HFRecHitCollection> HFRecHits;
  edm::Handle<EBRecHitCollection> EBRecHits;
  edm::Handle<EERecHitCollection> EERecHits;
  // the jet only contains DetIds, so first read recHit collection
  event.getByToken(input_HBHERecHits_token_, HBHERecHits);
  event.getByToken(input_HORecHits_token_, HORecHits);
  event.getByToken(input_HFRecHits_token_, HFRecHits);
  event.getByToken(input_EBRecHits_token_, EBRecHits);
  event.getByToken(input_EERecHits_token_, EERecHits);
  if (iDbg > 2)
    cout << "# of rechits found - HBHE: " << HBHERecHits->size() << ", HO: " << HORecHits->size()
         << ", HF: " << HFRecHits->size() << ", EB: " << EBRecHits->size() << ", EE: " << EERecHits->size() << endl;

  vector<CaloTowerPtr> towers = jet.getCaloConstituents();
  int nTowers = towers.size();
  if (iDbg > 9)
    cout << "In classifyJetComponents. # of towers found: " << nTowers << endl;

  for (int iTower = 0; iTower < nTowers; iTower++) {
    CaloTowerPtr &tower = towers[iTower];

    int nCells = tower->constituentsSize();
    if (iDbg)
      cout << "tower #" << iTower << " has " << nCells << " cells. "
           << "It's at iEta: " << tower->ieta() << ", iPhi: " << tower->iphi() << endl;

    const vector<DetId> &cellIDs = tower->constituents();  // cell == recHit

    for (int iCell = 0; iCell < nCells; ++iCell) {
      DetId::Detector detNum = cellIDs[iCell].det();
      if (detNum == DetId::Hcal) {
        HcalDetId HcalID = cellIDs[iCell];
        HcalSubdetector HcalNum = HcalID.subdet();
        double hitE = 0;
        if (HcalNum == HcalOuter) {
          HORecHitCollection::const_iterator theRecHit = HORecHits->find(HcalID);
          if (theRecHit == HORecHits->end()) {
            edm::LogWarning("UnexpectedEventContents") << "Can't find the HO recHit with ID: " << HcalID;
            continue;
          }
          hitE = theRecHit->energy();
          HO_energies.push_back(hitE);

        } else if (HcalNum == HcalForward) {
          HFRecHitCollection::const_iterator theRecHit = HFRecHits->find(HcalID);
          if (theRecHit == HFRecHits->end()) {
            edm::LogWarning("UnexpectedEventContents") << "Can't find the HF recHit with ID: " << HcalID;
            continue;
          }
          hitE = theRecHit->energy();
          if (iDbg > 4)
            cout << "hit #" << iCell << " is  HF , E: " << hitE << " iEta: " << theRecHit->id().ieta()
                 << ", depth: " << theRecHit->id().depth() << ", iPhi: " << theRecHit->id().iphi();

          if (HcalID.depth() == 1)
            long_energies.push_back(hitE);
          else
            short_energies.push_back(hitE);

          uint32_t flags = theRecHit->flags();
          if (flags & (1 << HcalCaloFlagLabels::HFLongShort))
            LS_bad_energy += hitE;
          if (flags & ((1 << HcalCaloFlagLabels::HFDigiTime) | (3 << HcalCaloFlagLabels::HFTimingTrustBits) |
                       (1 << HcalCaloFlagLabels::TimingSubtractedBit) | (1 << HcalCaloFlagLabels::TimingAddedBit) |
                       (1 << HcalCaloFlagLabels::TimingErrorBit)))
            HF_OOT_energy += hitE;
          if (iDbg > 4 && flags)
            cout << "flags: " << flags << " -> LS_bad_energy: " << LS_bad_energy << ", HF_OOT_energy: " << HF_OOT_energy
                 << endl;

        } else {  // HBHE

          HBHERecHitCollection::const_iterator theRecHit = HBHERecHits->find(HcalID);
          if (theRecHit == HBHERecHits->end()) {
            edm::LogWarning("UnexpectedEventContents") << "Can't find the HBHE recHit with ID: " << HcalID;
            continue;
          }
          hitE = theRecHit->energy();
          int iEta = theRecHit->id().ieta();
          int depth = theRecHit->id().depth();
          Region region = HBHE_region(theRecHit->id().rawId());
          int hitIPhi = theRecHit->id().iphi();
          if (iDbg > 3)
            cout << "hit #" << iCell << " is HBHE, E: " << hitE << " iEta: " << iEta << ", depth: " << depth
                 << ", iPhi: " << theRecHit->id().iphi() << " -> " << region;

          if (theHcalTopology->mergedDepth29(theRecHit->id()))
            hitE /= 2;  // Depth 3 at the HE forward edge is split over tower 28 & 29, and jet reco. assigns half each

          int iHPD = 100 * region;
          int iRBX = 100 * region + ((hitIPhi + 1) % 72) / 4;  // 71,72,1,2 are in the same RBX module

          if (std::abs(iEta) >= 21) {
            if ((0x1 & hitIPhi) == 0) {
              edm::LogError("CodeAssumptionsViolated") << "Bug?! Jet ID code assumes no even iPhi recHits at HE edges";
              return;
            }
            bool oddnessIEta = HBHE_oddness(iEta, depth);
            bool upperIPhi = ((hitIPhi % 4) == 1 || (hitIPhi % 4) == 2);  // the upper iPhi indices in the HE wedge
            // remap the iPhi so it fits the one in the inner HE regions, change in needed in two cases:
            // 1) in the upper iPhis of the module, the even iEtas belong to the higher iPhi
            // 2) in the loewr iPhis of the module, the odd  iEtas belong to the higher iPhi
            if (upperIPhi != oddnessIEta)
              ++hitIPhi;
            // note that hitIPhi could not be 72 before, so it's still in the legal range [1,72]
          }
          iHPD += hitIPhi;

          // book the energies
          HPD_energy_map[iHPD] += hitE;
          RBX_energy_map[iRBX] += hitE;
          if (iDbg > 5)
            cout << " --> H[" << iHPD << "]=" << HPD_energy_map[iHPD] << ", R[" << iRBX << "]=" << RBX_energy_map[iRBX];
          if (iDbg > 1)
            cout << endl;

          if (region == HBneg || region == HBpos)
            HB_energies.push_back(hitE);
          else
            HE_energies.push_back(hitE);

        }  // if HBHE
        if (hitE == 0)
          edm::LogWarning("UnexpectedEventContents") << "HCal hitE==0? (or unknown subdetector?)";

      }  // if HCAL

      else if (detNum == DetId::Ecal) {
        int EcalNum = cellIDs[iCell].subdetId();
        double hitE = 0;
        if (EcalNum == 1) {
          EBDetId EcalID = cellIDs[iCell];
          EBRecHitCollection::const_iterator theRecHit = EBRecHits->find(EcalID);
          if (theRecHit == EBRecHits->end()) {
            edm::LogWarning("UnexpectedEventContents") << "Can't find the EB recHit with ID: " << EcalID;
            continue;
          }
          hitE = theRecHit->energy();
          EB_energies.push_back(hitE);
        } else if (EcalNum == 2) {
          EEDetId EcalID = cellIDs[iCell];
          EERecHitCollection::const_iterator theRecHit = EERecHits->find(EcalID);
          if (theRecHit == EERecHits->end()) {
            edm::LogWarning("UnexpectedEventContents") << "Can't find the EE recHit with ID: " << EcalID;
            continue;
          }
          hitE = theRecHit->energy();
          EE_energies.push_back(hitE);
        }
        if (hitE == 0)
          edm::LogWarning("UnexpectedEventContents") << "ECal hitE==0? (or unknown subdetector?)";
        if (iDbg > 6)
          cout << "EcalNum: " << EcalNum << " hitE: " << hitE << endl;
      }  //
    }    // loop on cells
  }      // loop on towers

  /* Disabling check until HO is accounted for in EMF. Check was used in CMSSW_2, where HE was excluded.
  double expHcalE = jet.energy() * (1-jet.emEnergyFraction());
  if( totHcalE + expHcalE > 0 && 
      TMath::Abs( totHcalE - expHcalE ) > 0.01 && 
      ( totHcalE - expHcalE ) / ( totHcalE + expHcalE ) > 0.0001 ) {
    edm::LogWarning("CodeAssumptionsViolated")<<"failed to account for all Hcal energies"
                          <<totHcalE<<"!="<<expHcalE;
                          } */
  // concatenate Hcal and Ecal lists
  Hcal_energies.insert(Hcal_energies.end(), HB_energies.begin(), HB_energies.end());
  Hcal_energies.insert(Hcal_energies.end(), HE_energies.begin(), HE_energies.end());
  Hcal_energies.insert(Hcal_energies.end(), short_energies.begin(), short_energies.end());
  Hcal_energies.insert(Hcal_energies.end(), long_energies.begin(), long_energies.end());
  Ecal_energies.insert(Ecal_energies.end(), EB_energies.begin(), EB_energies.end());
  Ecal_energies.insert(Ecal_energies.end(), EE_energies.begin(), EE_energies.end());

  // sort the energies
  // std::sort( Hcal_energies.begin(), Hcal_energies.end(), greater<double>() );
  // std::sort( Ecal_energies.begin(), Ecal_energies.end(), greater<double>() );

  // put the energy sums (the 2nd entry in each pair of the maps) into the output vectors and sort them
  std::transform(
      HPD_energy_map.begin(), HPD_energy_map.end(), std::inserter(HPD_energies, HPD_energies.end()), select2nd);
  //          std::select2nd<std::map<int,double>::value_type>());
  // std::sort( HPD_energies.begin(), HPD_energies.end(), greater<double>() );
  std::transform(
      RBX_energy_map.begin(), RBX_energy_map.end(), std::inserter(RBX_energies, RBX_energies.end()), select2nd);
  //          std::select2nd<std::map<int,double>::value_type>());
  // std::sort( RBX_energies.begin(), RBX_energies.end(), greater<double>() );

  energies.insert(energies.end(), Hcal_energies.begin(), Hcal_energies.end());
  energies.insert(energies.end(), Ecal_energies.begin(), Ecal_energies.end());
  energies.insert(energies.end(), HO_energies.begin(), HO_energies.end());
  std::sort(energies.begin(), energies.end(), greater<double>());

  // prepare sub detector energies, then turn them into fractions
  fEB_ = std::accumulate(EB_energies.begin(), EB_energies.end(), 0.);
  fEE_ = std::accumulate(EE_energies.begin(), EE_energies.end(), 0.);
  fHB_ = std::accumulate(HB_energies.begin(), HB_energies.end(), 0.);
  fHE_ = std::accumulate(HE_energies.begin(), HE_energies.end(), 0.);
  fHO_ = std::accumulate(HO_energies.begin(), HO_energies.end(), 0.);
  fShort_ = std::accumulate(short_energies.begin(), short_energies.end(), 0.);
  fLong_ = std::accumulate(long_energies.begin(), long_energies.end(), 0.);
  subdet_energies.push_back(fEB_);
  subdet_energies.push_back(fEE_);
  subdet_energies.push_back(fHB_);
  subdet_energies.push_back(fHE_);
  subdet_energies.push_back(fHO_);
  subdet_energies.push_back(fShort_);
  subdet_energies.push_back(fLong_);
  std::sort(subdet_energies.begin(), subdet_energies.end(), greater<double>());
  if (jet.energy() > 0) {
    fEB_ /= jet.energy();
    fEE_ /= jet.energy();
    fHB_ /= jet.energy();
    fHE_ /= jet.energy();
    fHO_ /= jet.energy();
    fShort_ /= jet.energy();
    fLong_ /= jet.energy();
  } else {
    if (fEB_ > 0 || fEE_ > 0 || fHB_ > 0 || fHE_ > 0 || fHO_ > 0 || fShort_ > 0 || fLong_ > 0)
      edm::LogError("UnexpectedEventContents") << "Jet ID Helper found energy in subdetectors and jet E <= 0";
  }
}

void reco::helper::JetIDHelper::classifyJetTowers(const edm::Event &event,
                                                  const reco::CaloJet &jet,
                                                  vector<subtower> &subtowers,
                                                  vector<subtower> &Ecal_subtowers,
                                                  vector<subtower> &Hcal_subtowers,
                                                  vector<subtower> &HO_subtowers,
                                                  vector<double> &HPD_energies,
                                                  vector<double> &RBX_energies,
                                                  const int iDbg) {
  subtowers.clear();
  Ecal_subtowers.clear();
  Hcal_subtowers.clear();
  HO_subtowers.clear();
  HPD_energies.clear();
  RBX_energies.clear();

  std::map<int, double> HPD_energy_map, RBX_energy_map;

  vector<CaloTowerPtr> towers = jet.getCaloConstituents();
  int nTowers = towers.size();
  if (iDbg > 9)
    cout << "classifyJetTowers started. # of towers found: " << nTowers << endl;

  for (int iTower = 0; iTower < nTowers; iTower++) {
    CaloTowerPtr &tower = towers[iTower];

    int nEM = 0, nHad = 0, nHO = 0;
    const vector<DetId> &cellIDs = tower->constituents();  // cell == recHit
    int nCells = cellIDs.size();
    if (iDbg)
      cout << "tower #" << iTower << " has " << nCells << " cells. "
           << "It's at iEta: " << tower->ieta() << ", iPhi: " << tower->iphi() << endl;

    for (int iCell = 0; iCell < nCells; ++iCell) {
      DetId::Detector detNum = cellIDs[iCell].det();
      if (detNum == DetId::Hcal) {
        HcalDetId HcalID = cellIDs[iCell];
        HcalSubdetector HcalNum = HcalID.subdet();
        if (HcalNum == HcalOuter) {
          ++nHO;
        } else {
          ++nHad;
        }
      } else if (detNum == DetId::Ecal) {
        ++nEM;
      }
    }

    double E_em = tower->emEnergy();
    if (E_em != 0)
      Ecal_subtowers.push_back(subtower(E_em, nEM));

    double E_HO = tower->outerEnergy();
    if (E_HO != 0)
      HO_subtowers.push_back(subtower(E_HO, nHO));

    double E_had = tower->hadEnergy();
    if (E_had != 0) {
      Hcal_subtowers.push_back(subtower(E_had, nHad));
      // totHcalE += E_had;

      int iEta = tower->ieta();
      Region reg = region(iEta);
      int iPhi = tower->iphi();
      if (iDbg > 3)
        cout << "tower has E_had: " << E_had << " iEta: " << iEta << ", iPhi: " << iPhi << " -> " << reg;

      if (reg == HEneg || reg == HBneg || reg == HBpos || reg == HEpos) {
        int oddnessIEta = HBHE_oddness(iEta);
        if (oddnessIEta < 0)
          break;  // can't assign this tower to a single readout component

        int iHPD = 100 * reg;
        int iRBX = 100 * reg + ((iPhi + 1) % 72) / 4;  // 71,72,1,2 are in the same RBX module

        if ((reg == HEneg || reg == HEpos) && std::abs(iEta) >= 21) {  // at low-granularity edge of HE
          if ((0x1 & iPhi) == 0) {
            edm::LogError("CodeAssumptionsViolated") << "Bug?! Jet ID code assumes no even iPhi recHits at HE edges";
            return;
          }
          bool boolOddnessIEta = oddnessIEta;
          bool upperIPhi = ((iPhi % 4) == 1 || (iPhi % 4) == 2);  // the upper iPhi indices in the HE wedge
          // remap the iPhi so it fits the one in the inner HE regions, change in needed in two cases:
          // 1) in the upper iPhis of the module, the even IEtas belong to the higher iPhi
          // 2) in the loewr iPhis of the module, the odd  IEtas belong to the higher iPhi
          if (upperIPhi != boolOddnessIEta)
            ++iPhi;
          // note that iPhi could not be 72 before, so it's still in the legal range [1,72]
        }  // if at low-granularity edge of HE
        iHPD += iPhi;

        // book the energies
        HPD_energy_map[iHPD] += E_had;
        RBX_energy_map[iRBX] += E_had;
        if (iDbg > 5)
          cout << " --> H[" << iHPD << "]=" << HPD_energy_map[iHPD] << ", R[" << iRBX << "]=" << RBX_energy_map[iRBX];
      }  // HBHE
    }    // E_had > 0
  }      // loop on towers

  // sort the subtowers
  // std::sort( Hcal_subtowers.begin(), Hcal_subtowers.end(), subtower_has_greater_E );
  // std::sort( Ecal_subtowers.begin(), Ecal_subtowers.end(), subtower_has_greater_E );

  // put the energy sums (the 2nd entry in each pair of the maps) into the output vectors and sort them
  std::transform(
      HPD_energy_map.begin(), HPD_energy_map.end(), std::inserter(HPD_energies, HPD_energies.end()), select2nd);
  //          std::select2nd<std::map<int,double>::value_type>());
  // std::sort( HPD_energies.begin(), HPD_energies.end(), greater<double>() );
  std::transform(
      RBX_energy_map.begin(), RBX_energy_map.end(), std::inserter(RBX_energies, RBX_energies.end()), select2nd);
  //          std::select2nd<std::map<int,double>::value_type>());
  // std::sort( RBX_energies.begin(), RBX_energies.end(), greater<double>() );

  subtowers.insert(subtowers.end(), Hcal_subtowers.begin(), Hcal_subtowers.end());
  subtowers.insert(subtowers.end(), Ecal_subtowers.begin(), Ecal_subtowers.end());
  subtowers.insert(subtowers.end(), HO_subtowers.begin(), HO_subtowers.end());
  std::sort(subtowers.begin(), subtowers.end(), subtower_has_greater_E);
}

// ---------------------------------------------------------------------------------------------------------
// private helper functions to figure out detector & readout geometry as needed

int reco::helper::JetIDHelper::HBHE_oddness(int iEta, int depth) {
  int ae = TMath::Abs(iEta);
  if (ae == 29 && depth == 1)
    ae += 1;  // observed that: hits are at depths 1 & 2; 1 goes with the even pattern
  return ae & 0x1;
}

reco::helper::JetIDHelper::Region reco::helper::JetIDHelper::HBHE_region(uint32_t rawid) {
  HcalDetId id(rawid);
  if (id.subdet() == HcalEndcap) {
    if (id.ieta() < 0)
      return HEneg;
    else
      return HEpos;
  }
  if (id.subdet() == HcalBarrel && id.ieta() < 0)
    return HBneg;
  return HBpos;
}

int reco::helper::JetIDHelper::HBHE_oddness(int iEta) {
  int ae = TMath::Abs(iEta);
  if (ae == 29)
    return -1;  // can't figure it out without RecHits
  return ae & 0x1;
}

reco::helper::JetIDHelper::Region reco::helper::JetIDHelper::region(int iEta) {
  if (iEta == 16 || iEta == -16)
    return unknown_region;  // both HB and HE cells belong to these towers
  if (iEta == 29 || iEta == -29)
    return unknown_region;  // both HE and HF cells belong to these towers
  if (iEta <= -30)
    return HFneg;
  if (iEta >= 30)
    return HFpos;
  if (iEta <= -17)
    return HEneg;
  if (iEta >= 17)
    return HEpos;
  if (iEta < 0)
    return HBneg;
  return HBpos;
}