FindDistCone.cc

Go to the documentation of this file.

#include "Calibration/IsolatedParticles/interface/CaloConstants.h"
#include "Calibration/IsolatedParticles/interface/FindDistCone.h"
#include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
 
#include <iostream>
 
//#define EDM_ML_DEBUG
 
namespace spr {
 
  // Cone clustering core
  double getDistInPlaneTrackDir(const GlobalPoint& caloPoint,
                                const GlobalVector& caloVector,
                                const GlobalPoint& rechitPoint,
                                bool
#ifdef EDM_ML_DEBUG
                                    debug
#endif
  ) {
 
    const GlobalVector caloIntersectVector(caloPoint.x(), caloPoint.y(),
                                           caloPoint.z());  //p
 
    const GlobalVector caloUnitVector = caloVector.unit();
    const GlobalVector rechitVector(rechitPoint.x(), rechitPoint.y(), rechitPoint.z());
    const GlobalVector rechitUnitVector = rechitVector.unit();
    double dotprod_denominator = caloUnitVector.dot(rechitUnitVector);
    double dotprod_numerator = caloUnitVector.dot(caloIntersectVector);
    double rechitdist = dotprod_numerator / dotprod_denominator;
    const GlobalVector effectiveRechitVector = rechitdist * rechitUnitVector;
    const GlobalPoint effectiveRechitPoint(
        effectiveRechitVector.x(), effectiveRechitVector.y(), effectiveRechitVector.z());
    GlobalVector distance_vector = effectiveRechitPoint - caloPoint;
#ifdef EDM_ML_DEBUG
    if (debug) {
      std::cout << "getDistInPlaneTrackDir: point " << caloPoint << " dirn " << caloVector << " numerator "
                << dotprod_numerator << " denominator " << dotprod_denominator << " distance " << distance_vector.mag()
                << std::endl;
    }
#endif
    if (dotprod_denominator > 0. && dotprod_numerator > 0.) {
      return distance_vector.mag();
    } else {
      return 999999.;
    }
  }
 
  // Not used, but here for reference
  double getDistInCMatEcal(double eta1,
                           double phi1,
                           double eta2,
                           double phi2,
                           bool
#ifdef EDM_ML_DEBUG
                               debug
#endif
  ) {
 
    double dR, Rec;
    if (fabs(eta1) < spr::etaBEEcal)
      Rec = spr::rFrontEB;
    else
      Rec = spr::zFrontEE;
    double ce1 = cosh(eta1);
    double ce2 = cosh(eta2);
    double te1 = tanh(eta1);
    double te2 = tanh(eta2);
 
    double z = cos(phi1 - phi2) / ce1 / ce2 + te1 * te2;
    if (z != 0)
      dR = fabs(Rec * ce1 * sqrt(1. / z / z - 1.));
    else
      dR = 999999.;
#ifdef EDM_ML_DEBUG
    if (debug)
      std::cout << "getDistInCMatEcal: between (" << eta1 << ", " << phi1 << ") and (" << eta2 << ", " << phi2 << " is "
                << dR << std::endl;
#endif
    return dR;
  }
 
  // Not used, but here for reference
  double getDistInCMatHcal(double eta1,
                           double phi1,
                           double eta2,
                           double phi2,
                           bool
#ifdef EDM_ML_DEBUG
                               debug
#endif
  ) {
 
    // Radii and eta from Geometry/HcalCommonData/data/hcalendcapalgo.xml
    // and Geometry/HcalCommonData/data/hcalbarrelalgo.xml
 
    double dR, Rec;
    if (fabs(eta1) < spr::etaBEHcal)
      Rec = spr::rFrontHB;
    else
      Rec = spr::zFrontHE;
    double ce1 = cosh(eta1);
    double ce2 = cosh(eta2);
    double te1 = tanh(eta1);
    double te2 = tanh(eta2);
 
    double z = cos(phi1 - phi2) / ce1 / ce2 + te1 * te2;
    if (z != 0)
      dR = fabs(Rec * ce1 * sqrt(1. / z / z - 1.));
    else
      dR = 999999.;
    return dR;
#ifdef EDM_ML_DEBUG
    if (debug)
      std::cout << "getDistInCMatHcal: between (" << eta1 << ", " << phi1 << ") and (" << eta2 << ", " << phi2 << " is "
                << dR << std::endl;
#endif
  }
 
  void getEtaPhi(HBHERecHitCollection::const_iterator hit,
                 std::vector<int>& RH_ieta,
                 std::vector<int>& RH_iphi,
                 std::vector<double>& RH_ene,
                 bool) {
    RH_ieta.push_back(hit->id().ieta());
    RH_iphi.push_back(hit->id().iphi());
    RH_ene.push_back(hit->energy());
  }
 
  void getEtaPhi(edm::PCaloHitContainer::const_iterator hit,
                 std::vector<int>& RH_ieta,
                 std::vector<int>& RH_iphi,
                 std::vector<double>& RH_ene,
                 bool) {
    // SimHit function not yet implemented.
    RH_ieta.push_back(-9);
    RH_iphi.push_back(-9);
    RH_ene.push_back(-9.);
  }
 
  void getEtaPhi(EcalRecHitCollection::const_iterator hit,
                 std::vector<int>& RH_ieta,
                 std::vector<int>& RH_iphi,
                 std::vector<double>& RH_ene,
                 bool) {
    // Ecal function not yet implemented.
    RH_ieta.push_back(-9);
    RH_iphi.push_back(-9);
    RH_ene.push_back(-9.);
  }
 
  void getEtaPhi(HBHERecHitCollection::const_iterator hit, int& ieta, int& iphi, bool) {
    ieta = hit->id().ieta();
    iphi = hit->id().iphi();
  }
 
  void getEtaPhi(edm::PCaloHitContainer::const_iterator hit, int& ieta, int& iphi, bool) {
    DetId id = DetId(hit->id());
    if (id.det() == DetId::Hcal) {
      ieta = ((HcalDetId)(hit->id())).ieta();
      iphi = ((HcalDetId)(hit->id())).iphi();
    } else if (id.det() == DetId::Ecal && id.subdetId() == EcalBarrel) {
      ieta = ((EBDetId)(id)).ieta();
      iphi = ((EBDetId)(id)).iphi();
    } else {
      ieta = 999;
      iphi = 999;
    }
  }
 
  void getEtaPhi(EcalRecHitCollection::const_iterator hit, int& ieta, int& iphi, bool) {
    DetId id = hit->id();
    if (id.subdetId() == EcalBarrel) {
      ieta = ((EBDetId)(id)).ieta();
      iphi = ((EBDetId)(id)).iphi();
    } else {
      ieta = 999;
      iphi = 999;
    }
  }
 
  double getEnergy(HBHERecHitCollection::const_iterator hit, int useRaw, bool) {
    double energy = ((useRaw == 1) ? hit->eraw() : ((useRaw == 2) ? hit->eaux() : hit->energy()));
    return energy;
  }
 
  double getEnergy(EcalRecHitCollection::const_iterator hit, int, bool) { return hit->energy(); }
 
  double getEnergy(edm::PCaloHitContainer::const_iterator hit, int, bool) {
    // This will not yet handle Ecal CaloHits!!
    double samplingWeight = 1.;
    // Hard coded sampling weights from JFH analysis of iso tracks
    // Sept 2009.
    DetId id = hit->id();
    if (id.det() == DetId::Hcal) {
      HcalDetId detId(id);
      if (detId.subdet() == HcalBarrel)
        samplingWeight = 114.1;
      else if (detId.subdet() == HcalEndcap)
        samplingWeight = 167.3;
      else {
        // ONLY protection against summing HO, HF simhits
        return 0.;
      }
    }
 
    return samplingWeight * hit->energy();
  }
 
  GlobalPoint getGpos(const CaloGeometry* geo, HBHERecHitCollection::const_iterator hit, bool) {
    DetId detId(hit->id());
    return (static_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(detId)))->getPosition(detId);
  }
 
  GlobalPoint getGpos(const CaloGeometry* geo, edm::PCaloHitContainer::const_iterator hit, bool) {
    DetId detId(hit->id());
    GlobalPoint point = (detId.det() == DetId::Hcal)
                            ? (static_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(detId)))->getPosition(detId)
                            : geo->getPosition(detId);
    return point;
  }
 
  GlobalPoint getGpos(const CaloGeometry* geo, EcalRecHitCollection::const_iterator hit, bool) {
    // Not tested for EcalRecHits!!
    if (hit->id().subdetId() == EcalEndcap) {
      EEDetId EEid = EEDetId(hit->id());
      return geo->getPosition(EEid);
    } else {  // (hit->id().subdetId() == EcalBarrel)
      EBDetId EBid = EBDetId(hit->id());
      return geo->getPosition(EBid);
    }
  }
 
  double getRawEnergy(HBHERecHitCollection::const_iterator hit, int useRaw) {
    double energy = ((useRaw == 1) ? hit->eraw() : ((useRaw == 2) ? hit->eaux() : hit->energy()));
    return energy;
  }
 
  double getRawEnergy(EcalRecHitCollection::const_iterator hit, int) { return hit->energy(); }
 
  double getRawEnergy(edm::PCaloHitContainer::const_iterator hit, int) { return hit->energy(); }
}  // namespace spr