HcalHaloAlgo.cc

Go to the documentation of this file.

#include "RecoMET/METAlgorithms/interface/HcalHaloAlgo.h"
#include <map>
 
/*
  [class]:  HcalHaloAlgo
  [authors]: R. Remington, The University of Florida
  [description]: See HcalHaloAlgo.h
  [date]: October 15, 2009
*/
namespace {
  constexpr float c_cm_per_ns = 29.9792458;
  constexpr float zseparation_HBHE = 380.;
};  // namespace
 
using namespace reco;
 
#include <iomanip>
bool CompareTime(const HBHERecHit* x, const HBHERecHit* y) { return x->time() < y->time(); }
bool CompareTowers(const CaloTower* x, const CaloTower* y) {
  return x->iphi() * 1000 + x->ieta() < y->iphi() * 1000 + y->ieta();
}
 
HcalHaloAlgo::HcalHaloAlgo() : geo_(nullptr), hgeo_(nullptr) {
  HBRecHitEnergyThreshold = 0.;
  HERecHitEnergyThreshold = 0.;
  SumEnergyThreshold = 0.;
  NHitsThreshold = 0;
}
 
HcalHaloData HcalHaloAlgo::Calculate(const CaloGeometry& TheCaloGeometry,
                                     edm::Handle<HBHERecHitCollection>& TheHBHERecHits,
                                     edm::Handle<EBRecHitCollection>& TheEBRecHits,
                                     edm::Handle<EERecHitCollection>& TheEERecHits,
                                     const edm::EventSetup& TheSetup) {
  edm::Handle<CaloTowerCollection> TheCaloTowers;
  return Calculate(TheCaloGeometry, TheHBHERecHits, TheCaloTowers, TheEBRecHits, TheEERecHits, TheSetup);
}
 
HcalHaloData HcalHaloAlgo::Calculate(const CaloGeometry& TheCaloGeometry,
                                     edm::Handle<HBHERecHitCollection>& TheHBHERecHits,
                                     edm::Handle<CaloTowerCollection>& TheCaloTowers,
                                     edm::Handle<EBRecHitCollection>& TheEBRecHits,
                                     edm::Handle<EERecHitCollection>& TheEERecHits,
                                     const edm::EventSetup& TheSetup) {
  HcalHaloData TheHcalHaloData;
  // ieta overlap geometrically w/ HB
  const int iEtaOverlap = 22;
  const int nPhiMax = 73;
  // Store Energy sum of rechits as a function of iPhi (iPhi goes from 1 to 72)
  float SumE[nPhiMax];
  // Store Number of rechits as a function of iPhi
  int NumHits[nPhiMax];
  // Store minimum time of rechit as a function of iPhi
  float MinTimeHits[nPhiMax];
  // Store maximum time of rechit as a function of iPhi
  float MaxTimeHits[nPhiMax];
  for (unsigned int i = 0; i < nPhiMax; i++) {
    SumE[i] = 0;
    NumHits[i] = 0;
    MinTimeHits[i] = 0.;
    MaxTimeHits[i] = 0.;
  }
 
  for (const auto& hit : (*TheHBHERecHits)) {
    HcalDetId id = HcalDetId(hit.id());
    switch (id.subdet()) {
      case HcalBarrel:
        if (hit.energy() < HBRecHitEnergyThreshold)
          continue;
        break;
      case HcalEndcap:
        if (hit.energy() < HERecHitEnergyThreshold)
          continue;
        break;
      default:
        continue;
    }
 
    int iEta = id.ieta();
    int iPhi = id.iphi();
    if (iPhi < nPhiMax && std::abs(iEta) <= iEtaOverlap) {
      SumE[iPhi] += hit.energy();
      NumHits[iPhi]++;
 
      float time = hit.time();
      MinTimeHits[iPhi] = time < MinTimeHits[iPhi] ? time : MinTimeHits[iPhi];
      MaxTimeHits[iPhi] = time > MaxTimeHits[iPhi] ? time : MaxTimeHits[iPhi];
    }
  }
 
  for (int iPhi = 1; iPhi < nPhiMax; iPhi++) {
    if (SumE[iPhi] >= SumEnergyThreshold && NumHits[iPhi] > NHitsThreshold) {
      // Build PhiWedge and store to HcalHaloData if energy or #hits pass thresholds
      PhiWedge wedge(SumE[iPhi], iPhi, NumHits[iPhi], MinTimeHits[iPhi], MaxTimeHits[iPhi]);
 
      // Loop over rechits again to calculate direction based on timing info
      std::vector<const HBHERecHit*> Hits;
      for (const auto& hit : (*TheHBHERecHits)) {
        HcalDetId id = HcalDetId(hit.id());
        if (id.iphi() != iPhi)
          continue;
        if (std::abs(id.ieta()) > iEtaOverlap)
          continue;  // has to overlap geometrically w/ HB
        switch (id.subdet()) {
          case HcalBarrel:
            if (hit.energy() < HBRecHitEnergyThreshold)
              continue;
            break;
          case HcalEndcap:
            if (hit.energy() < HERecHitEnergyThreshold)
              continue;
            break;
          default:
            continue;
        }
        Hits.push_back(&(hit));
      }
 
      std::sort(Hits.begin(), Hits.end(), CompareTime);
      float MinusToPlus = 0.;
      float PlusToMinus = 0.;
      for (unsigned int i = 0; i < Hits.size(); i++) {
        HcalDetId id_i = HcalDetId(Hits[i]->id());
        int ieta_i = id_i.ieta();
        for (unsigned int j = (i + 1); j < Hits.size(); j++) {
          HcalDetId id_j = HcalDetId(Hits[j]->id());
          int ieta_j = id_j.ieta();
          if (ieta_i > ieta_j)
            PlusToMinus += std::abs(ieta_i - ieta_j);
          else
            MinusToPlus += std::abs(ieta_i - ieta_j);
        }
      }
      float PlusZOriginConfidence = (PlusToMinus + MinusToPlus) ? PlusToMinus / (PlusToMinus + MinusToPlus) : -1.;
      wedge.SetPlusZOriginConfidence(PlusZOriginConfidence);
      TheHcalHaloData.GetPhiWedges().push_back(wedge);
    }
  }
 
  // Don't use HF.
  int maxAbsIEta = 29;
 
  std::map<int, float> iPhiHadEtMap;
  std::vector<const CaloTower*> sortedCaloTowers;
  for (const auto& tower : (*TheCaloTowers)) {
    if (std::abs(tower.ieta()) > maxAbsIEta)
      continue;
 
    int iPhi = tower.iphi();
    if (!iPhiHadEtMap.count(iPhi))
      iPhiHadEtMap[iPhi] = 0.0;
    iPhiHadEtMap[iPhi] += tower.hadEt();
 
    if (tower.numProblematicHcalCells() > 0)
      sortedCaloTowers.push_back(&(tower));
  }
 
  // Sort towers such that lowest iphi and ieta are first, highest last, and towers
  // with same iphi value are consecutive. Then we can do everything else in one loop.
  std::sort(sortedCaloTowers.begin(), sortedCaloTowers.end(), CompareTowers);
 
  HaloTowerStrip strip;
 
  int prevIEta = -99, prevIPhi = -99;
  float prevHadEt = 0.;
  float prevEmEt = 0.;
  std::pair<uint8_t, CaloTowerDetId> prevPair, towerPair;
  bool wasContiguous = true;
 
  // Loop through and store a vector of pairs (problematicCells, DetId) for each contiguous strip we find
  for (unsigned int i = 0; i < sortedCaloTowers.size(); i++) {
    const CaloTower* tower = sortedCaloTowers[i];
 
    towerPair = std::make_pair((uint8_t)tower->numProblematicHcalCells(), tower->id());
 
    bool newIPhi = tower->iphi() != prevIPhi;
    bool isContiguous = tower->ieta() == 1 ? tower->ieta() - 2 == prevIEta : tower->ieta() - 1 == prevIEta;
 
    isContiguous = isContiguous || (tower->ieta() == -maxAbsIEta);
    if (newIPhi)
      isContiguous = false;
 
    if (!wasContiguous && isContiguous) {
      strip.cellTowerIds.push_back(prevPair);
      strip.cellTowerIds.push_back(towerPair);
      strip.hadEt += prevHadEt + tower->hadEt();
      strip.emEt += prevEmEt + tower->emEt();
    }
 
    if (wasContiguous && isContiguous) {
      strip.cellTowerIds.push_back(towerPair);
      strip.hadEt += tower->hadEt();
      strip.emEt += tower->emEt();
    }
 
    if ((wasContiguous && !isContiguous) || i == sortedCaloTowers.size() - 1) {  //ended the strip, so flush it
 
      if (strip.cellTowerIds.size() > 3) {
        int iPhi = strip.cellTowerIds.at(0).second.iphi();
        int iPhiLower = (iPhi == 1) ? 72 : iPhi - 1;
        int iPhiUpper = (iPhi == 72) ? 1 : iPhi + 1;
 
        float energyRatio = 0.0;
        if (iPhiHadEtMap.count(iPhiLower))
          energyRatio += iPhiHadEtMap[iPhiLower];
        if (iPhiHadEtMap.count(iPhiUpper))
          energyRatio += iPhiHadEtMap[iPhiUpper];
        iPhiHadEtMap[iPhi] = std::max(iPhiHadEtMap[iPhi], 0.001F);
 
        energyRatio /= iPhiHadEtMap[iPhi];
        strip.energyRatio = energyRatio;
 
        TheHcalHaloData.getProblematicStrips().push_back(strip);
      }
      strip = HaloTowerStrip();
    }
 
    wasContiguous = isContiguous;
    prevPair = towerPair;
    prevEmEt = tower->emEt();
    prevIPhi = tower->iphi();
    prevIEta = tower->ieta();
    prevHadEt = tower->hadEt();
  }
 
  edm::ESHandle<CaloGeometry> pGeo;
  TheSetup.get<CaloGeometryRecord>().get(pGeo);
  geo_ = pGeo.product();
  hgeo_ = dynamic_cast<const HcalGeometry*>(geo_->getSubdetectorGeometry(DetId::Hcal, 1));
 
  //Halo cluster building:
  //Various clusters are built, depending on the subdetector.
  //In barrel, one looks for deposits narrow in phi.
  //In endcaps, one looks for localized deposits (dr condition in EE where r =sqrt(dphi*dphi+deta*deta)
  //E/H condition is also applied.
  //The halo cluster building step targets a large efficiency (ideally >99%) for beam halo deposits.
  //These clusters are used as input for the halo pattern finding methods in HcalHaloAlgo and for the CSC-calo matching methods in GlobalHaloAlgo.
 
  //Et threshold hardcoded for now. Might one to get it from config
 
  std::vector<HaloClusterCandidateHCAL> haloclustercands_HB;
  haloclustercands_HB = GetHaloClusterCandidateHB(TheEBRecHits, TheHBHERecHits, 5);
 
  std::vector<HaloClusterCandidateHCAL> haloclustercands_HE;
  haloclustercands_HE = GetHaloClusterCandidateHE(TheEERecHits, TheHBHERecHits, 10);
 
  TheHcalHaloData.setHaloClusterCandidatesHB(haloclustercands_HB);
  TheHcalHaloData.setHaloClusterCandidatesHE(haloclustercands_HE);
 
  return TheHcalHaloData;
}
 
std::vector<HaloClusterCandidateHCAL> HcalHaloAlgo::GetHaloClusterCandidateHB(
    edm::Handle<EcalRecHitCollection>& ecalrechitcoll,
    edm::Handle<HBHERecHitCollection>& hbherechitcoll,
    float et_thresh_seedrh) {
  std::vector<HaloClusterCandidateHCAL> TheHaloClusterCandsHB;
 
  reco::Vertex::Point vtx(0, 0, 0);
 
  for (size_t ihit = 0; ihit < hbherechitcoll->size(); ++ihit) {
    HaloClusterCandidateHCAL clustercand;
 
    const HBHERecHit& rechit = (*hbherechitcoll)[ihit];
    math::XYZPoint rhpos = getPosition(rechit.id(), vtx);
    //Et condition
    double rhet = rechit.energy() * sqrt(rhpos.perp2() / rhpos.mag2());
    if (rhet < et_thresh_seedrh)
      continue;
    if (std::abs(rhpos.z()) > zseparation_HBHE)
      continue;
    double eta = rhpos.eta();
    double phi = rhpos.phi();
 
    bool isiso = true;
    double etcluster(0);
    int nbtowerssameeta(0);
    double timediscriminatorITBH(0), timediscriminatorOTBH(0);
    double etstrip_phiseedplus1(0), etstrip_phiseedminus1(0);
 
    //Building the cluster
    edm::RefVector<HBHERecHitCollection> bhrhcandidates;
    for (size_t jhit = 0; jhit < hbherechitcoll->size(); ++jhit) {
      const HBHERecHit& rechitj = (*hbherechitcoll)[jhit];
      HBHERecHitRef rhRef(hbherechitcoll, jhit);
      math::XYZPoint rhposj = getPosition(rechitj.id(), vtx);
      double rhetj = rechitj.energy() * sqrt(rhposj.perp2() / rhposj.mag2());
      if (rhetj < 2)
        continue;
      if (std::abs(rhposj.z()) > zseparation_HBHE)
        continue;
      double etaj = rhposj.eta();
      double phij = rhposj.phi();
      double deta = eta - etaj;
      double dphi = deltaPhi(phi, phij);
      if (std::abs(deta) > 0.4)
        continue;  //This means +/-4 towers in eta
      if (std::abs(dphi) > 0.2)
        continue;  //This means +/-2 towers in phi
      if (std::abs(dphi) > 0.1 && std::abs(deta) < 0.2) {
        isiso = false;
        break;
      }  //The strip should be isolated
      if (std::abs(dphi) > 0.1)
        continue;
      if (std::abs(dphi) < 0.05)
        nbtowerssameeta++;
      if (dphi > 0.05)
        etstrip_phiseedplus1 += rhetj;
      if (dphi < -0.05)
        etstrip_phiseedminus1 += rhetj;
 
      etcluster += rhetj;
      //Timing discriminator
      //We assign a weight to the rechit defined as:
      //Log10(Et)*f(T,R,Z)
      //where f(T,R,Z) is the separation curve between halo-like and IP-like times.
      //The time difference between a deposit from a outgoing IT halo and a deposit coming from a particle emitted at the IP is given by:
      //dt= ( - sqrt(R^2+z^2) + |z| )/c
      // For OT beam halo, the time difference is:
      //dt= ( 25 + sqrt(R^2+z^2) + |z| )/c
      //only consider the central part of HB as things get hard at large z.
      //The best fitted value for R leads to 240 cm (IT) and 330 cm (OT)
      double rhtj = rechitj.time();
      timediscriminatorITBH +=
          std::log10(rhetj) *
          (rhtj + 0.5 * (sqrt(240. * 240. + rhposj.z() * rhposj.z()) - std::abs(rhposj.z())) / c_cm_per_ns);
      if (std::abs(rhposj.z()) < 300)
        timediscriminatorOTBH +=
            std::log10(rhetj) *
            (rhtj - 0.5 * (25 - (sqrt(330. * 330. + rhposj.z() * rhposj.z()) + std::abs(rhposj.z())) / c_cm_per_ns));
      bhrhcandidates.push_back(rhRef);
    }
    //Isolation conditions
    if (!isiso)
      continue;
    if (etstrip_phiseedplus1 / etcluster > 0.2 && etstrip_phiseedminus1 / etcluster > 0.2)
      continue;
 
    //Calculate E/H
    double eoh(0);
    for (size_t jhit = 0; jhit < ecalrechitcoll->size(); ++jhit) {
      const EcalRecHit& rechitj = (*ecalrechitcoll)[jhit];
      math::XYZPoint rhposj = getPosition(rechitj.id(), vtx);
      double rhetj = rechitj.energy() * sqrt(rhposj.perp2() / rhposj.mag2());
      if (rhetj < 2)
        continue;
      double etaj = rhposj.eta();
      double phij = rhposj.phi();
      if (std::abs(eta - etaj) > 0.2)
        continue;
      if (std::abs(deltaPhi(phi, phij)) > 0.2)
        continue;
      eoh += rhetj / etcluster;
    }
    //E/H condition
    if (eoh > 0.1)
      continue;
 
    clustercand.setClusterEt(etcluster);
    clustercand.setSeedEt(rhet);
    clustercand.setSeedEta(eta);
    clustercand.setSeedPhi(phi);
    clustercand.setSeedZ(rhpos.Z());
    clustercand.setSeedR(sqrt(rhpos.perp2()));
    clustercand.setSeedTime(rechit.time());
    clustercand.setEoverH(eoh);
    clustercand.setNbTowersInEta(nbtowerssameeta);
    clustercand.setEtStripPhiSeedPlus1(etstrip_phiseedplus1);
    clustercand.setEtStripPhiSeedMinus1(etstrip_phiseedminus1);
    clustercand.setTimeDiscriminatorITBH(timediscriminatorITBH);
    clustercand.setTimeDiscriminatorOTBH(timediscriminatorOTBH);
    clustercand.setBeamHaloRecHitsCandidates(bhrhcandidates);
 
    bool isbeamhalofrompattern = HBClusterShapeandTimeStudy(clustercand, false);
    clustercand.setIsHaloFromPattern(isbeamhalofrompattern);
    bool isbeamhalofrompattern_hlt = HBClusterShapeandTimeStudy(clustercand, true);
    clustercand.setIsHaloFromPattern_HLT(isbeamhalofrompattern_hlt);
 
    TheHaloClusterCandsHB.push_back(clustercand);
  }
 
  return TheHaloClusterCandsHB;
}
 
std::vector<HaloClusterCandidateHCAL> HcalHaloAlgo::GetHaloClusterCandidateHE(
    edm::Handle<EcalRecHitCollection>& ecalrechitcoll,
    edm::Handle<HBHERecHitCollection>& hbherechitcoll,
    float et_thresh_seedrh) {
  std::vector<HaloClusterCandidateHCAL> TheHaloClusterCandsHE;
 
  reco::Vertex::Point vtx(0, 0, 0);
 
  for (size_t ihit = 0; ihit < hbherechitcoll->size(); ++ihit) {
    HaloClusterCandidateHCAL clustercand;
 
    const HBHERecHit& rechit = (*hbherechitcoll)[ihit];
    math::XYZPoint rhpos = getPosition(rechit.id(), vtx);
    //Et condition
    double rhet = rechit.energy() * sqrt(rhpos.perp2() / rhpos.mag2());
    if (rhet < et_thresh_seedrh)
      continue;
    if (std::abs(rhpos.z()) < zseparation_HBHE)
      continue;
    double eta = rhpos.eta();
    double phi = rhpos.phi();
    double rhr = sqrt(rhpos.perp2());
    bool isiso = true;
    double etcluster(0), hdepth1(0);
    int clustersize(0);
    double etstrip_phiseedplus1(0), etstrip_phiseedminus1(0);
 
    //Building the cluster
    edm::RefVector<HBHERecHitCollection> bhrhcandidates;
    for (size_t jhit = 0; jhit < hbherechitcoll->size(); ++jhit) {
      const HBHERecHit& rechitj = (*hbherechitcoll)[jhit];
      HBHERecHitRef rhRef(hbherechitcoll, jhit);
      math::XYZPoint rhposj = getPosition(rechitj.id(), vtx);
      double rhetj = rechitj.energy() * sqrt(rhposj.perp2() / rhposj.mag2());
      if (rhetj < 2)
        continue;
      if (std::abs(rhposj.z()) < zseparation_HBHE)
        continue;
      if (rhpos.z() * rhposj.z() < 0)
        continue;
      double phij = rhposj.phi();
      double dphi = deltaPhi(phi, phij);
      if (std::abs(dphi) > 0.4)
        continue;
      double rhrj = sqrt(rhposj.perp2());
      if (std::abs(rhr - rhrj) > 50)
        continue;
      if (std::abs(dphi) > 0.2 || std::abs(rhr - rhrj) > 20) {
        isiso = false;
        break;
      }  //The deposit should be isolated
      if (dphi > 0.05)
        etstrip_phiseedplus1 += rhetj;
      if (dphi < -0.05)
        etstrip_phiseedminus1 += rhetj;
      clustersize++;
      etcluster += rhetj;
      if (std::abs(rhposj.z()) < 405)
        hdepth1 += rhetj;
      //No timing condition for now in HE
      bhrhcandidates.push_back(rhRef);
    }
    //Isolation conditions
    if (!isiso)
      continue;
    if (etstrip_phiseedplus1 / etcluster > 0.1 && etstrip_phiseedminus1 / etcluster > 0.1)
      continue;
 
    //Calculate E/H
    double eoh(0);
    for (size_t jhit = 0; jhit < ecalrechitcoll->size(); ++jhit) {
      const EcalRecHit& rechitj = (*ecalrechitcoll)[jhit];
      math::XYZPoint rhposj = getPosition(rechitj.id(), vtx);
      double rhetj = rechitj.energy() * sqrt(rhposj.perp2() / rhposj.mag2());
      if (rhetj < 2)
        continue;
      if (rhpos.z() * rhposj.z() < 0)
        continue;
      double etaj = rhposj.eta();
      double phij = rhposj.phi();
      double dr = sqrt((eta - etaj) * (eta - etaj) + deltaPhi(phi, phij) * deltaPhi(phi, phij));
      if (dr > 0.3)
        continue;
 
      eoh += rhetj / etcluster;
    }
    //E/H condition
    if (eoh > 0.1)
      continue;
 
    clustercand.setClusterEt(etcluster);
    clustercand.setSeedEt(rhet);
    clustercand.setSeedEta(eta);
    clustercand.setSeedPhi(phi);
    clustercand.setSeedZ(rhpos.Z());
    clustercand.setSeedR(sqrt(rhpos.perp2()));
    clustercand.setSeedTime(rechit.time());
    clustercand.setEoverH(eoh);
    clustercand.setH1overH123(hdepth1 / etcluster);
    clustercand.setClusterSize(clustersize);
    clustercand.setEtStripPhiSeedPlus1(etstrip_phiseedplus1);
    clustercand.setEtStripPhiSeedMinus1(etstrip_phiseedminus1);
    clustercand.setTimeDiscriminator(0);
    clustercand.setBeamHaloRecHitsCandidates(bhrhcandidates);
 
    bool isbeamhalofrompattern = HEClusterShapeandTimeStudy(clustercand, false);
    clustercand.setIsHaloFromPattern(isbeamhalofrompattern);
    bool isbeamhalofrompattern_hlt = HEClusterShapeandTimeStudy(clustercand, true);
    clustercand.setIsHaloFromPattern_HLT(isbeamhalofrompattern_hlt);
 
    TheHaloClusterCandsHE.push_back(clustercand);
  }
 
  return TheHaloClusterCandsHE;
}
 
bool HcalHaloAlgo::HBClusterShapeandTimeStudy(HaloClusterCandidateHCAL hcand, bool ishlt) {
  //Conditions on the central strip size in eta.
  //For low size, extra conditions on seed et, isolation and cluster timing
  //Here we target both IT and OT beam halo. Two separate discriminators were built for the two cases.
 
  if (hcand.getSeedEt() < 10)
    return false;
 
  if (hcand.getNbTowersInEta() < 3)
    return false;
  //Isolation criteria for very short eta strips
  if (hcand.getNbTowersInEta() == 3 && (hcand.getEtStripPhiSeedPlus1() > 0.1 || hcand.getEtStripPhiSeedMinus1() > 0.1))
    return false;
  if (hcand.getNbTowersInEta() <= 5 && (hcand.getEtStripPhiSeedPlus1() > 0.1 && hcand.getEtStripPhiSeedMinus1() > 0.1))
    return false;
 
  //Timing conditions for short eta strips
  if (hcand.getNbTowersInEta() == 3 && hcand.getTimeDiscriminatorITBH() >= 0.)
    return false;
  if (hcand.getNbTowersInEta() <= 6 && hcand.getTimeDiscriminatorITBH() >= 5. && hcand.getTimeDiscriminatorOTBH() < 0.)
    return false;
 
  //For HLT, only use conditions without timing
  if (ishlt && hcand.getNbTowersInEta() < 7)
    return false;
 
  hcand.setIsHaloFromPattern(true);
 
  return true;
}
 
bool HcalHaloAlgo::HEClusterShapeandTimeStudy(HaloClusterCandidateHCAL hcand, bool ishlt) {
  //Conditions on H1/H123 to spot halo interacting only in one HCAL layer.
  //For R> about 170cm, HE has only one layer and this condition cannot be applied
  //Note that for R>170 cm, the halo is in CSC acceptance and will most likely be spotted by the CSC-calo matching method
  //A method to identify halos interacting in both H1 and H2/H3 at low R is still missing.
 
  if (hcand.getSeedEt() < 20)
    return false;
  if (hcand.getSeedR() > 170)
    return false;
 
  if (hcand.getH1overH123() > 0.02 && hcand.getH1overH123() < 0.98)
    return false;
 
  //This method is one of the ones with the highest fake rate: in JetHT dataset, it happens in around 0.1% of the cases that a low pt jet (pt= 20) leaves all of its energy in only one HCAL layer.
  //At HLT, one only cares about large deposits from BH that would lead to a MET/SinglePhoton trigger to be fired.
  //Rising the seed Et threshold at HLT has therefore little impact on the HLT performances but ensures that possible controversial events are still recorded.
  if (ishlt && hcand.getSeedEt() < 50)
    return false;
 
  hcand.setIsHaloFromPattern(true);
 
  return true;
}
 
math::XYZPoint HcalHaloAlgo::getPosition(const DetId& id, reco::Vertex::Point vtx) {
  const GlobalPoint pos = ((id.det() == DetId::Hcal) ? hgeo_->getPosition(id) : GlobalPoint(geo_->getPosition(id)));
  math::XYZPoint posV(pos.x() - vtx.x(), pos.y() - vtx.y(), pos.z() - vtx.z());
  return posV;
}