GlobalHaloAlgo.cc

Go to the documentation of this file.

#include "RecoMET/METAlgorithms/interface/GlobalHaloAlgo.h"
namespace {
  constexpr float c_cm_per_ns = 29.9792458;
};
/*
  [class]:  GlobalHaloAlgo
  [authors]: R. Remington, The University of Florida
  [description]: See GlobalHaloAlgo.h
  [date]: October 15, 2009
*/
using namespace std;
using namespace edm;
using namespace reco;

enum detectorregion { EB, EE, HB, HE };
int Phi_To_HcaliPhi(float phi) {
  phi = phi < 0 ? phi + 2. * TMath::Pi() : phi;
  float phi_degrees = phi * (360.) / (2. * TMath::Pi());
  int iPhi = (int)((phi_degrees / 5.) + 1.);

  return iPhi < 73 ? iPhi : 73;
}

int Phi_To_EcaliPhi(float phi) {
  phi = phi < 0 ? phi + 2. * TMath::Pi() : phi;
  float phi_degrees = phi * (360.) / (2. * TMath::Pi());
  int iPhi = (int)(phi_degrees + 1.);

  return iPhi < 361 ? iPhi : 360;
}

GlobalHaloAlgo::GlobalHaloAlgo() {
  // Defaults are "loose"
  Ecal_R_Min = 110.;  // Tight: 200.
  Ecal_R_Max = 330.;  // Tight: 250.
  Hcal_R_Min = 110.;  // Tight: 220.
  Hcal_R_Max = 490.;  // Tight: 350.
}

reco::GlobalHaloData GlobalHaloAlgo::Calculate(const CaloGeometry& TheCaloGeometry,
                                               const CSCGeometry& TheCSCGeometry,
                                               const reco::CaloMET& TheCaloMET,
                                               edm::Handle<edm::View<Candidate> >& TheCaloTowers,
                                               edm::Handle<CSCSegmentCollection>& TheCSCSegments,
                                               edm::Handle<CSCRecHit2DCollection>& TheCSCRecHits,
                                               edm::Handle<reco::MuonCollection>& TheMuons,
                                               const CSCHaloData& TheCSCHaloData,
                                               const EcalHaloData& TheEcalHaloData,
                                               const HcalHaloData& TheHcalHaloData,
                                               bool ishlt) {
  GlobalHaloData TheGlobalHaloData;
  float METOverSumEt = TheCaloMET.sumEt() ? TheCaloMET.pt() / TheCaloMET.sumEt() : 0;
  TheGlobalHaloData.SetMETOverSumEt(METOverSumEt);

  int EcalOverlapping_CSCRecHits[361] = {};
  int EcalOverlapping_CSCSegments[361] = {};
  int HcalOverlapping_CSCRecHits[73] = {};
  int HcalOverlapping_CSCSegments[73] = {};

  if (TheCSCSegments.isValid()) {
    for (CSCSegmentCollection::const_iterator iSegment = TheCSCSegments->begin(); iSegment != TheCSCSegments->end();
         iSegment++) {
      bool EcalOverlap[361];
      bool HcalOverlap[73];
      for (int i = 0; i < 361; i++) {
        EcalOverlap[i] = false;
        if (i < 73)
          HcalOverlap[i] = false;
      }

      std::vector<CSCRecHit2D> Hits = iSegment->specificRecHits();
      for (std::vector<CSCRecHit2D>::iterator iHit = Hits.begin(); iHit != Hits.end(); iHit++) {
        DetId TheDetUnitId(iHit->geographicalId());
        if (TheDetUnitId.det() != DetId::Muon)
          continue;
        if (TheDetUnitId.subdetId() != MuonSubdetId::CSC)
          continue;

        const GeomDetUnit* TheUnit = TheCSCGeometry.idToDetUnit(TheDetUnitId);
        LocalPoint TheLocalPosition = iHit->localPosition();
        const BoundPlane& TheSurface = TheUnit->surface();
        const GlobalPoint TheGlobalPosition = TheSurface.toGlobal(TheLocalPosition);

        int Hcal_iphi = Phi_To_HcaliPhi(TheGlobalPosition.phi());
        int Ecal_iphi = Phi_To_EcaliPhi(TheGlobalPosition.phi());
        float x = TheGlobalPosition.x();
        float y = TheGlobalPosition.y();

        float r = TMath::Sqrt(x * x + y * y);

        if (r < Ecal_R_Max && r > Ecal_R_Min)
          EcalOverlap[Ecal_iphi] = true;
        if (r < Hcal_R_Max && r > Hcal_R_Max)
          HcalOverlap[Hcal_iphi] = true;
      }
      for (int i = 0; i < 361; i++) {
        if (EcalOverlap[i])
          EcalOverlapping_CSCSegments[i]++;
        if (i < 73 && HcalOverlap[i])
          HcalOverlapping_CSCSegments[i]++;
      }
    }
  }
  if (TheCSCRecHits.isValid()) {
    for (CSCRecHit2DCollection::const_iterator iCSCRecHit = TheCSCRecHits->begin(); iCSCRecHit != TheCSCRecHits->end();
         iCSCRecHit++) {
      DetId TheDetUnitId(iCSCRecHit->geographicalId());
      if (TheDetUnitId.det() != DetId::Muon)
        continue;
      if (TheDetUnitId.subdetId() != MuonSubdetId::CSC)
        continue;

      const GeomDetUnit* TheUnit = TheCSCGeometry.idToDetUnit(TheDetUnitId);
      LocalPoint TheLocalPosition = iCSCRecHit->localPosition();
      const BoundPlane& TheSurface = TheUnit->surface();
      const GlobalPoint TheGlobalPosition = TheSurface.toGlobal(TheLocalPosition);

      int Hcaliphi = Phi_To_HcaliPhi(TheGlobalPosition.phi());
      int Ecaliphi = Phi_To_EcaliPhi(TheGlobalPosition.phi());
      float x = TheGlobalPosition.x();
      float y = TheGlobalPosition.y();

      float r = TMath::Sqrt(x * x + y * y);

      if (r < Ecal_R_Max && r > Ecal_R_Min)
        EcalOverlapping_CSCRecHits[Ecaliphi]++;
      if (r < Hcal_R_Max && r > Hcal_R_Max)
        HcalOverlapping_CSCRecHits[Hcaliphi]++;
    }
  }

  // In development....
  // Get Ecal Wedges
  std::vector<PhiWedge> EcalWedges = TheEcalHaloData.GetPhiWedges();

  // Get Hcal Wedges
  std::vector<PhiWedge> HcalWedges = TheHcalHaloData.GetPhiWedges();

  //Get Ref to CSC Tracks
  //edm::RefVector<reco::TrackCollection> TheCSCTracks = TheCSCHaloData.GetTracks();
  //for(unsigned int i = 0 ; i < TheCSCTracks.size() ; i++ )
  //edm::Ref<reco::TrackCollection> iTrack( TheCSCTracks, i );

  // Get global positions of central most rechit of CSC Halo tracks
  std::vector<GlobalPoint> TheGlobalPositions = TheCSCHaloData.GetCSCTrackImpactPositions();

  // Container to store Ecal/Hcal iPhi values matched to impact point of CSC tracks
  std::vector<int> vEcaliPhi, vHcaliPhi;

  // Keep track of number of calo pointing CSC halo tracks that do not match to Phi wedges
  int N_Unmatched_Tracks = 0;

  for (std::vector<GlobalPoint>::iterator Pos = TheGlobalPositions.begin(); Pos != TheGlobalPositions.end(); Pos++) {
    // Calculate global phi coordinate for central most rechit in the track
    float global_phi = Pos->phi();
    float global_r = TMath::Sqrt(Pos->x() * Pos->x() + Pos->y() * Pos->y());

    // Convert global phi to iPhi
    int global_EcaliPhi = Phi_To_EcaliPhi(global_phi);
    int global_HcaliPhi = Phi_To_HcaliPhi(global_phi);
    bool MATCHED = false;

    //Loop over Ecal Phi Wedges
    for (std::vector<PhiWedge>::iterator iWedge = EcalWedges.begin(); iWedge != EcalWedges.end(); iWedge++) {
      if ((TMath::Abs(global_EcaliPhi - iWedge->iPhi()) <= 5) && (global_r > Ecal_R_Min && global_r < Ecal_R_Max)) {
        bool StoreWedge = true;
        for (unsigned int i = 0; i < vEcaliPhi.size(); i++)
          if (vEcaliPhi[i] == iWedge->iPhi())
            StoreWedge = false;

        if (StoreWedge) {
          PhiWedge NewWedge(*iWedge);
          NewWedge.SetOverlappingCSCSegments(EcalOverlapping_CSCSegments[iWedge->iPhi()]);
          NewWedge.SetOverlappingCSCRecHits(EcalOverlapping_CSCRecHits[iWedge->iPhi()]);
          vEcaliPhi.push_back(iWedge->iPhi());
          TheGlobalHaloData.GetMatchedEcalPhiWedges().push_back(NewWedge);
        }
        MATCHED = true;
      }
    }
    //Loop over Hcal Phi Wedges
    for (std::vector<PhiWedge>::iterator iWedge = HcalWedges.begin(); iWedge != HcalWedges.end(); iWedge++) {
      if ((TMath::Abs(global_HcaliPhi - iWedge->iPhi()) <= 2) && (global_r > Hcal_R_Min && global_r < Hcal_R_Max)) {
        bool StoreWedge = true;
        for (unsigned int i = 0; i < vHcaliPhi.size(); i++)
          if (vHcaliPhi[i] == iWedge->iPhi())
            StoreWedge = false;

        if (StoreWedge) {
          vHcaliPhi.push_back(iWedge->iPhi());
          PhiWedge NewWedge(*iWedge);
          NewWedge.SetOverlappingCSCSegments(HcalOverlapping_CSCSegments[iWedge->iPhi()]);
          NewWedge.SetOverlappingCSCRecHits(HcalOverlapping_CSCRecHits[iWedge->iPhi()]);
          PhiWedge wedge(*iWedge);
          TheGlobalHaloData.GetMatchedHcalPhiWedges().push_back(NewWedge);
        }
        MATCHED = true;
      }
    }
    if (!MATCHED)
      N_Unmatched_Tracks++;
  }

  // Corrections to MEx, MEy
  float dMEx = 0.;
  float dMEy = 0.;
  // Loop over calotowers and correct the MET for the towers that lie in the trajectory of the CSC Halo Tracks
  for (edm::View<Candidate>::const_iterator iCandidate = TheCaloTowers->begin(); iCandidate != TheCaloTowers->end();
       iCandidate++) {
    const Candidate* c = &(*iCandidate);
    if (c) {
      const CaloTower* iTower = dynamic_cast<const CaloTower*>(c);
      if (iTower->et() < TowerEtThreshold)
        continue;
      if (abs(iTower->ieta()) > 24)
        continue;  // not in barrel/endcap
      int iphi = iTower->iphi();
      for (unsigned int x = 0; x < vEcaliPhi.size(); x++) {
        if (iphi == vEcaliPhi[x]) {
          dMEx += (TMath::Cos(iTower->phi()) * iTower->emEt());
          dMEy += (TMath::Sin(iTower->phi()) * iTower->emEt());
        }
      }
      for (unsigned int x = 0; x < vHcaliPhi.size(); x++) {
        if (iphi == vHcaliPhi[x]) {
          dMEx += (TMath::Cos(iTower->phi()) * iTower->hadEt());
          dMEy += (TMath::Sin(iTower->phi()) * iTower->hadEt());
        }
      }
    }
  }

  TheGlobalHaloData.SetMETCorrections(dMEx, dMEy);

  std::vector<HaloClusterCandidateECAL> hccandEB = TheEcalHaloData.getHaloClusterCandidatesEB();
  std::vector<HaloClusterCandidateECAL> hccandEE = TheEcalHaloData.getHaloClusterCandidatesEE();
  std::vector<HaloClusterCandidateHCAL> hccandHB = TheHcalHaloData.getHaloClusterCandidatesHB();
  std::vector<HaloClusterCandidateHCAL> hccandHE = TheHcalHaloData.getHaloClusterCandidatesHE();

  //CSC-calo matching
  bool ECALBmatched(false), ECALEmatched(false), HCALBmatched(false), HCALEmatched(false);

  if (TheCSCSegments.isValid()) {
    for (CSCSegmentCollection::const_iterator iSegment = TheCSCSegments->begin(); iSegment != TheCSCSegments->end();
         iSegment++) {
      CSCDetId iCscDetID = iSegment->cscDetId();
      bool Segment1IsGood = true;

      //avoid segments from collision muons
      if (TheMuons.isValid()) {
        for (reco::MuonCollection::const_iterator mu = TheMuons->begin(); mu != TheMuons->end() && (Segment1IsGood);
             mu++) {
          if (!mu->isTrackerMuon() && !mu->isGlobalMuon() && mu->isStandAloneMuon())
            continue;
          if (!mu->isGlobalMuon() && mu->isTrackerMuon() && mu->pt() < 3)
            continue;
          const std::vector<MuonChamberMatch> chambers = mu->matches();
          for (std::vector<MuonChamberMatch>::const_iterator kChamber = chambers.begin(); kChamber != chambers.end();
               kChamber++) {
            if (kChamber->detector() != MuonSubdetId::CSC)
              continue;
            for (std::vector<reco::MuonSegmentMatch>::const_iterator kSegment = kChamber->segmentMatches.begin();
                 kSegment != kChamber->segmentMatches.end();
                 kSegment++) {
              edm::Ref<CSCSegmentCollection> cscSegRef = kSegment->cscSegmentRef;
              CSCDetId kCscDetID = cscSegRef->cscDetId();

              if (kCscDetID == iCscDetID) {
                Segment1IsGood = false;
              }
            }
          }
        }
      }
      if (!Segment1IsGood)
        continue;

      // Get local direction vector; if direction runs parallel to beamline,
      // count this segment as beam halo candidate.
      LocalPoint iLocalPosition = iSegment->localPosition();
      LocalVector iLocalDirection = iSegment->localDirection();

      GlobalPoint iGlobalPosition = TheCSCGeometry.chamber(iCscDetID)->toGlobal(iLocalPosition);
      GlobalVector iGlobalDirection = TheCSCGeometry.chamber(iCscDetID)->toGlobal(iLocalDirection);

      float iTheta = iGlobalDirection.theta();
      if (iTheta > max_segment_theta && iTheta < TMath::Pi() - max_segment_theta)
        continue;

      float iPhi = iGlobalPosition.phi();
      float iR = sqrt(iGlobalPosition.perp2());
      float iZ = iGlobalPosition.z();
      float iT = iSegment->time();

      //CSC-calo matching:
      //Here, one checks if any halo cluster can be matched to a CSC segment.
      //The matching uses both geometric (dphi, dR) and timing information (dt).
      //The cut values depend on the subdetector considered (e.g. in HB, Rcalo-Rsegment is allowed to be very negative)

      bool ebmatched = SegmentMatchingEB(TheGlobalHaloData, hccandEB, iZ, iR, iT, iPhi, ishlt);
      bool eematched = SegmentMatchingEE(TheGlobalHaloData, hccandEE, iZ, iR, iT, iPhi, ishlt);
      bool hbmatched = SegmentMatchingHB(TheGlobalHaloData, hccandHB, iZ, iR, iT, iPhi, ishlt);
      bool hematched = SegmentMatchingHE(TheGlobalHaloData, hccandHE, iZ, iR, iT, iPhi, ishlt);

      ECALBmatched |= ebmatched;
      ECALEmatched |= eematched;
      HCALBmatched |= hbmatched;
      HCALEmatched |= hematched;
    }
  }

  TheGlobalHaloData.SetSegmentIsEBCaloMatched(ECALBmatched);
  TheGlobalHaloData.SetSegmentIsEECaloMatched(ECALEmatched);
  TheGlobalHaloData.SetSegmentIsHBCaloMatched(HCALBmatched);
  TheGlobalHaloData.SetSegmentIsHECaloMatched(HCALEmatched);

  //Now checking patterns from EcalHaloData and HcalHaloData:
  //Simply check whether any cluster has a halo pattern
  //In that case store the rhits in GlobalHaloData

  bool HaloPatternFoundInEB = false;
  for (auto& hcand : hccandEB) {
    if ((hcand.getIsHaloFromPattern() && !ishlt) || (hcand.getIsHaloFromPattern_HLT() && ishlt)) {
      HaloPatternFoundInEB = true;
      edm::RefVector<EcalRecHitCollection> bhrhcandidates = hcand.getBeamHaloRecHitsCandidates();
      AddtoBeamHaloEBEERechits(bhrhcandidates, TheGlobalHaloData, true);
    }
  }

  bool HaloPatternFoundInEE = false;
  for (auto& hcand : hccandEE) {
    if ((hcand.getIsHaloFromPattern() && !ishlt) || (hcand.getIsHaloFromPattern_HLT() && ishlt)) {
      HaloPatternFoundInEE = true;
      edm::RefVector<EcalRecHitCollection> bhrhcandidates = hcand.getBeamHaloRecHitsCandidates();
      AddtoBeamHaloEBEERechits(bhrhcandidates, TheGlobalHaloData, false);
    }
  }

  bool HaloPatternFoundInHB = false;
  for (auto& hcand : hccandHB) {
    if ((hcand.getIsHaloFromPattern() && !ishlt) || (hcand.getIsHaloFromPattern_HLT() && ishlt)) {
      HaloPatternFoundInHB = true;
      edm::RefVector<HBHERecHitCollection> bhrhcandidates = hcand.getBeamHaloRecHitsCandidates();
      AddtoBeamHaloHBHERechits(bhrhcandidates, TheGlobalHaloData);
    }
  }

  bool HaloPatternFoundInHE = false;
  for (auto& hcand : hccandHE) {
    if ((hcand.getIsHaloFromPattern() && !ishlt) || (hcand.getIsHaloFromPattern_HLT() && ishlt)) {
      HaloPatternFoundInHE = true;
      edm::RefVector<HBHERecHitCollection> bhrhcandidates = hcand.getBeamHaloRecHitsCandidates();
      AddtoBeamHaloHBHERechits(bhrhcandidates, TheGlobalHaloData);
    }
  }
  TheGlobalHaloData.SetHaloPatternFoundEB(HaloPatternFoundInEB);
  TheGlobalHaloData.SetHaloPatternFoundEE(HaloPatternFoundInEE);
  TheGlobalHaloData.SetHaloPatternFoundHB(HaloPatternFoundInHB);
  TheGlobalHaloData.SetHaloPatternFoundHE(HaloPatternFoundInHE);

  return TheGlobalHaloData;
}

bool GlobalHaloAlgo::SegmentMatchingEB(reco::GlobalHaloData& thehalodata,
                                       const std::vector<HaloClusterCandidateECAL>& haloclustercands,
                                       float iZ,
                                       float iR,
                                       float iT,
                                       float iPhi,
                                       bool ishlt) {
  bool rhmatchingfound = false;

  for (auto& hcand : haloclustercands) {
    if (!ApplyMatchingCuts(EB,
                           ishlt,
                           hcand.getSeedEt(),
                           iZ,
                           hcand.getSeedZ(),
                           iR,
                           hcand.getSeedR(),
                           iT,
                           hcand.getSeedTime(),
                           iPhi,
                           hcand.getSeedPhi()))
      continue;

    rhmatchingfound = true;

    edm::RefVector<EcalRecHitCollection> bhrhcandidates = hcand.getBeamHaloRecHitsCandidates();

    AddtoBeamHaloEBEERechits(bhrhcandidates, thehalodata, true);
  }

  return rhmatchingfound;
}

bool GlobalHaloAlgo::SegmentMatchingEE(reco::GlobalHaloData& thehalodata,
                                       const std::vector<HaloClusterCandidateECAL>& haloclustercands,
                                       float iZ,
                                       float iR,
                                       float iT,
                                       float iPhi,
                                       bool ishlt) {
  bool rhmatchingfound = false;

  for (auto& hcand : haloclustercands) {
    if (!ApplyMatchingCuts(EE,
                           ishlt,
                           hcand.getSeedEt(),
                           iZ,
                           hcand.getSeedZ(),
                           iR,
                           hcand.getSeedR(),
                           iT,
                           hcand.getSeedTime(),
                           iPhi,
                           hcand.getSeedPhi()))
      continue;

    rhmatchingfound = true;

    edm::RefVector<EcalRecHitCollection> bhrhcandidates = hcand.getBeamHaloRecHitsCandidates();

    AddtoBeamHaloEBEERechits(bhrhcandidates, thehalodata, false);
  }

  return rhmatchingfound;
}

bool GlobalHaloAlgo::SegmentMatchingHB(reco::GlobalHaloData& thehalodata,
                                       const std::vector<HaloClusterCandidateHCAL>& haloclustercands,
                                       float iZ,
                                       float iR,
                                       float iT,
                                       float iPhi,
                                       bool ishlt) {
  bool rhmatchingfound = false;

  for (auto& hcand : haloclustercands) {
    if (!ApplyMatchingCuts(HB,
                           ishlt,
                           hcand.getSeedEt(),
                           iZ,
                           hcand.getSeedZ(),
                           iR,
                           hcand.getSeedR(),
                           iT,
                           hcand.getSeedTime(),
                           iPhi,
                           hcand.getSeedPhi()))
      continue;

    rhmatchingfound = true;

    edm::RefVector<HBHERecHitCollection> bhrhcandidates = hcand.getBeamHaloRecHitsCandidates();

    AddtoBeamHaloHBHERechits(bhrhcandidates, thehalodata);
  }

  return rhmatchingfound;
}

bool GlobalHaloAlgo::SegmentMatchingHE(reco::GlobalHaloData& thehalodata,
                                       const std::vector<HaloClusterCandidateHCAL>& haloclustercands,
                                       float iZ,
                                       float iR,
                                       float iT,
                                       float iPhi,
                                       bool ishlt) {
  bool rhmatchingfound = false;

  for (auto& hcand : haloclustercands) {
    if (!ApplyMatchingCuts(HE,
                           ishlt,
                           hcand.getSeedEt(),
                           iZ,
                           hcand.getSeedZ(),
                           iR,
                           hcand.getSeedR(),
                           iT,
                           hcand.getSeedTime(),
                           iPhi,
                           hcand.getSeedPhi()))
      continue;

    rhmatchingfound = true;

    edm::RefVector<HBHERecHitCollection> bhrhcandidates = hcand.getBeamHaloRecHitsCandidates();

    AddtoBeamHaloHBHERechits(bhrhcandidates, thehalodata);
  }

  return rhmatchingfound;
}

bool GlobalHaloAlgo::ApplyMatchingCuts(int subdet,
                                       bool ishlt,
                                       double rhet,
                                       double segZ,
                                       double rhZ,
                                       double segR,
                                       double rhR,
                                       double segT,
                                       double rhT,
                                       double segPhi,
                                       double rhPhi) {
  //Std::Absolute time wrt BX
  double tBXrh = rhT + sqrt(rhR * rhR + rhZ * rhZ) / c_cm_per_ns;
  double tBXseg = segT + sqrt(segR * segR + segZ * segZ) / c_cm_per_ns;
  //Time at z=0, under beam halo hypothesis
  double tcorseg = tBXseg - std::abs(segZ) / c_cm_per_ns;       //Outgoing beam halo
  double tcorsegincbh = tBXseg + std::abs(segZ) / c_cm_per_ns;  //Ingoing beam halo
  double truedt[4] = {1000, 1000, 1000, 1000};
  //There are four types of segments associated to beam halo, test each hypothesis:
  //IT beam halo, ingoing track
  double twindow_seg = 15;
  if (std::abs(tcorsegincbh) < twindow_seg)
    truedt[0] = tBXrh - tBXseg - std::abs(rhZ - segZ) / c_cm_per_ns;
  //IT beam halo, outgoing track
  if (std::abs(tcorseg) < twindow_seg)
    truedt[1] = tBXseg - tBXrh - std::abs(rhZ - segZ) / c_cm_per_ns;
  //OT beam halo (from next BX), ingoing track
  if (tcorsegincbh > 25 - twindow_seg && std::abs(tcorsegincbh) < 25 + twindow_seg)
    truedt[2] = tBXrh - tBXseg - std::abs(rhZ - segZ) / c_cm_per_ns;
  //OT beam halo (from next BX), outgoing track
  if (tcorseg > 25 - twindow_seg && tcorseg < 25 + twindow_seg)
    truedt[3] = tBXseg - tBXrh - std::abs(rhZ - segZ) / c_cm_per_ns;

  if (subdet == EB) {
    if (rhet < et_thresh_rh_eb)
      return false;
    if (rhet < 20 && ishlt)
      return false;
    if (std::abs(deltaPhi(rhPhi, segPhi)) > dphi_thresh_segvsrh_eb)
      return false;
    if (rhR - segR < dr_lowthresh_segvsrh_eb)
      return false;
    if (rhR - segR > dr_highthresh_segvsrh_eb)
      return false;
    if (std::abs(truedt[0]) > dt_segvsrh_eb && std::abs(truedt[1]) > dt_segvsrh_eb &&
        std::abs(truedt[2]) > dt_segvsrh_eb && std::abs(truedt[3]) > dt_segvsrh_eb)
      return false;
    return true;
  }

  if (subdet == EE) {
    if (rhet < et_thresh_rh_ee)
      return false;
    if (rhet < 20 && ishlt)
      return false;
    if (std::abs(deltaPhi(rhPhi, segPhi)) > dphi_thresh_segvsrh_ee)
      return false;
    if (rhR - segR < dr_lowthresh_segvsrh_ee)
      return false;
    if (rhR - segR > dr_highthresh_segvsrh_ee)
      return false;
    if (std::abs(truedt[0]) > dt_segvsrh_ee && std::abs(truedt[1]) > dt_segvsrh_ee &&
        std::abs(truedt[2]) > dt_segvsrh_ee && std::abs(truedt[3]) > dt_segvsrh_ee)
      return false;
    return true;
  }

  if (subdet == HB) {
    if (rhet < et_thresh_rh_hb)
      return false;
    if (rhet < 20 && ishlt)
      return false;
    if (std::abs(deltaPhi(rhPhi, segPhi)) > dphi_thresh_segvsrh_hb)
      return false;
    if (rhR - segR < dr_lowthresh_segvsrh_hb)
      return false;
    if (rhR - segR > dr_highthresh_segvsrh_hb)
      return false;
    if (std::abs(truedt[0]) > dt_segvsrh_hb && std::abs(truedt[1]) > dt_segvsrh_hb &&
        std::abs(truedt[2]) > dt_segvsrh_hb && std::abs(truedt[3]) > dt_segvsrh_hb)
      return false;
    return true;
  }

  if (subdet == HE) {
    if (rhet < et_thresh_rh_he)
      return false;
    if (rhet < 20 && ishlt)
      return false;
    if (std::abs(deltaPhi(rhPhi, segPhi)) > dphi_thresh_segvsrh_he)
      return false;
    if (rhR - segR < dr_lowthresh_segvsrh_he)
      return false;
    if (rhR - segR > dr_highthresh_segvsrh_he)
      return false;
    if (std::abs(truedt[0]) > dt_segvsrh_he && std::abs(truedt[1]) > dt_segvsrh_he &&
        std::abs(truedt[2]) > dt_segvsrh_he && std::abs(truedt[3]) > dt_segvsrh_he)
      return false;
    return true;
  }

  return false;
}

void GlobalHaloAlgo::AddtoBeamHaloEBEERechits(edm::RefVector<EcalRecHitCollection>& bhtaggedrechits,
                                              reco::GlobalHaloData& thehalodata,
                                              bool isbarrel) {
  for (size_t ihit = 0; ihit < bhtaggedrechits.size(); ++ihit) {
    bool alreadyincl = false;
    edm::Ref<EcalRecHitCollection> rhRef(bhtaggedrechits[ihit]);
    edm::RefVector<EcalRecHitCollection> refrhcoll;
    if (isbarrel)
      refrhcoll = thehalodata.GetEBRechits();
    else
      refrhcoll = thehalodata.GetEERechits();
    for (size_t jhit = 0; jhit < refrhcoll.size(); jhit++) {
      edm::Ref<EcalRecHitCollection> rhitRef(refrhcoll[jhit]);
      if (rhitRef->detid() == rhRef->detid())
        alreadyincl = true;
      if (rhitRef->detid() == rhRef->detid())
        break;
    }
    if (!alreadyincl && isbarrel)
      thehalodata.GetEBRechits().push_back(rhRef);
    if (!alreadyincl && !isbarrel)
      thehalodata.GetEERechits().push_back(rhRef);
  }
}

void GlobalHaloAlgo::AddtoBeamHaloHBHERechits(edm::RefVector<HBHERecHitCollection>& bhtaggedrechits,
                                              reco::GlobalHaloData& thehalodata) {
  for (size_t ihit = 0; ihit < bhtaggedrechits.size(); ++ihit) {
    bool alreadyincl = false;
    edm::Ref<HBHERecHitCollection> rhRef(bhtaggedrechits[ihit]);
    edm::RefVector<HBHERecHitCollection> refrhcoll;
    refrhcoll = thehalodata.GetHBHERechits();
    for (size_t jhit = 0; jhit < refrhcoll.size(); jhit++) {
      edm::Ref<HBHERecHitCollection> rhitRef(refrhcoll[jhit]);
      if (rhitRef->detid() == rhRef->detid())
        alreadyincl = true;
      if (rhitRef->detid() == rhRef->detid())
        break;
    }
    if (!alreadyincl)
      thehalodata.GetHBHERechits().push_back(rhRef);
  }
}