IsolatedParticles/src/FindDistCone.cc

0001 #include "Calibration/IsolatedParticles/interface/CaloConstants.h"
0002 #include "Calibration/IsolatedParticles/interface/FindDistCone.h"
0003 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005
0006 namespace spr {
0007
0008   // Cone clustering core
0009   double getDistInPlaneTrackDir(const GlobalPoint& caloPoint,
0010                                 const GlobalVector& caloVector,
0011                                 const GlobalPoint& rechitPoint,
0012                                 bool debug) {
0013     const GlobalVector caloIntersectVector(caloPoint.x(), caloPoint.y(),
0014                                            caloPoint.z());  //p
0015
0016     const GlobalVector caloUnitVector = caloVector.unit();
0017     const GlobalVector rechitVector(rechitPoint.x(), rechitPoint.y(), rechitPoint.z());
0018     const GlobalVector rechitUnitVector = rechitVector.unit();
0019     double dotprod_denominator = caloUnitVector.dot(rechitUnitVector);
0020     double dotprod_numerator = caloUnitVector.dot(caloIntersectVector);
0021     double rechitdist = dotprod_numerator / dotprod_denominator;
0022     const GlobalVector effectiveRechitVector = rechitdist * rechitUnitVector;
0023     const GlobalPoint effectiveRechitPoint(
0024         effectiveRechitVector.x(), effectiveRechitVector.y(), effectiveRechitVector.z());
0025     GlobalVector distance_vector = effectiveRechitPoint - caloPoint;
0026     if (debug) {
0027       edm::LogVerbatim("IsoTrack") << "getDistInPlaneTrackDir: point " << caloPoint << " dirn " << caloVector
0028                                    << " numerator " << dotprod_numerator << " denominator " << dotprod_denominator
0029                                    << " distance " << distance_vector.mag();
0030     }
0031     if (dotprod_denominator > 0. && dotprod_numerator > 0.) {
0032       return distance_vector.mag();
0033     } else {
0034       return 999999.;
0035     }
0036   }
0037
0038   // Not used, but here for reference
0039   double getDistInCMatEcal(double eta1, double phi1, double eta2, double phi2, bool debug) {
0040     double dR, Rec;
0041     if (fabs(eta1) < spr::etaBEEcal)
0042       Rec = spr::rFrontEB;
0043     else
0044       Rec = spr::zFrontEE;
0045     double ce1 = cosh(eta1);
0046     double ce2 = cosh(eta2);
0047     double te1 = tanh(eta1);
0048     double te2 = tanh(eta2);
0049
0050     double z = cos(phi1 - phi2) / ce1 / ce2 + te1 * te2;
0051     if (z != 0)
0052       dR = fabs(Rec * ce1 * sqrt(1. / z / z - 1.));
0053     else
0054       dR = 999999.;
0055     if (debug)
0056       edm::LogVerbatim("IsoTrack") << "getDistInCMatEcal: between (" << eta1 << ", " << phi1 << ") and (" << eta2
0057                                    << ", " << phi2 << " is " << dR;
0058     return dR;
0059   }
0060
0061   // Not used, but here for reference
0062   double getDistInCMatHcal(double eta1, double phi1, double eta2, double phi2, bool debug) {
0063     // Radii and eta from Geometry/HcalCommonData/data/hcalendcapalgo.xml
0064     // and Geometry/HcalCommonData/data/hcalbarrelalgo.xml
0065
0066     double dR, Rec;
0067     if (fabs(eta1) < spr::etaBEHcal)
0068       Rec = spr::rFrontHB;
0069     else
0070       Rec = spr::zFrontHE;
0071     double ce1 = cosh(eta1);
0072     double ce2 = cosh(eta2);
0073     double te1 = tanh(eta1);
0074     double te2 = tanh(eta2);
0075
0076     double z = cos(phi1 - phi2) / ce1 / ce2 + te1 * te2;
0077     if (z != 0)
0078       dR = fabs(Rec * ce1 * sqrt(1. / z / z - 1.));
0079     else
0080       dR = 999999.;
0081     if (debug)
0082       edm::LogVerbatim("IsoTrack") << "getDistInCMatHcal: between (" << eta1 << ", " << phi1 << ") and (" << eta2
0083                                    << ", " << phi2 << " is " << dR;
0084     return dR;
0085   }
0086
0087   void getEtaPhi(HBHERecHitCollection::const_iterator hit,
0088                  std::vector<int>& RH_ieta,
0089                  std::vector<int>& RH_iphi,
0090                  std::vector<double>& RH_ene,
0091                  bool) {
0092     RH_ieta.push_back(hit->id().ieta());
0093     RH_iphi.push_back(hit->id().iphi());
0094     RH_ene.push_back(hit->energy());
0095   }
0096
0097   void getEtaPhi(edm::PCaloHitContainer::const_iterator hit,
0098                  std::vector<int>& RH_ieta,
0099                  std::vector<int>& RH_iphi,
0100                  std::vector<double>& RH_ene,
0101                  bool) {
0102     // SimHit function not yet implemented.
0103     RH_ieta.push_back(-9);
0104     RH_iphi.push_back(-9);
0105     RH_ene.push_back(-9.);
0106   }
0107
0108   void getEtaPhi(EcalRecHitCollection::const_iterator hit,
0109                  std::vector<int>& RH_ieta,
0110                  std::vector<int>& RH_iphi,
0111                  std::vector<double>& RH_ene,
0112                  bool) {
0113     // Ecal function not yet implemented.
0114     RH_ieta.push_back(-9);
0115     RH_iphi.push_back(-9);
0116     RH_ene.push_back(-9.);
0117   }
0118
0119   void getEtaPhi(HBHERecHitCollection::const_iterator hit, int& ieta, int& iphi, bool) {
0120     ieta = hit->id().ieta();
0121     iphi = hit->id().iphi();
0122   }
0123
0124   void getEtaPhi(edm::PCaloHitContainer::const_iterator hit, int& ieta, int& iphi, bool) {
0125     DetId id = DetId(hit->id());
0126     if (id.det() == DetId::Hcal) {
0127       ieta = ((HcalDetId)(hit->id())).ieta();
0128       iphi = ((HcalDetId)(hit->id())).iphi();
0129     } else if (id.det() == DetId::Ecal && id.subdetId() == EcalBarrel) {
0130       ieta = ((EBDetId)(id)).ieta();
0131       iphi = ((EBDetId)(id)).iphi();
0132     } else {
0133       ieta = 999;
0134       iphi = 999;
0135     }
0136   }
0137
0138   void getEtaPhi(EcalRecHitCollection::const_iterator hit, int& ieta, int& iphi, bool) {
0139     DetId id = hit->id();
0140     if (id.subdetId() == EcalBarrel) {
0141       ieta = ((EBDetId)(id)).ieta();
0142       iphi = ((EBDetId)(id)).iphi();
0143     } else {
0144       ieta = 999;
0145       iphi = 999;
0146     }
0147   }
0148
0149   double getEnergy(HBHERecHitCollection::const_iterator hit, int useRaw, bool) {
0150     double energy = ((useRaw == 1) ? hit->eraw() : ((useRaw == 2) ? hit->eaux() : hit->energy()));
0151     return energy;
0152   }
0153
0154   double getEnergy(EcalRecHitCollection::const_iterator hit, int, bool) { return hit->energy(); }
0155
0156   double getEnergy(edm::PCaloHitContainer::const_iterator hit, int, bool) {
0157     // This will not yet handle Ecal CaloHits!!
0158     double samplingWeight = 1.;
0159     // Hard coded sampling weights from JFH analysis of iso tracks
0160     // Sept 2009.
0161     DetId id = hit->id();
0162     if (id.det() == DetId::Hcal) {
0163       HcalDetId detId(id);
0164       if (detId.subdet() == HcalBarrel)
0165         samplingWeight = 114.1;
0166       else if (detId.subdet() == HcalEndcap)
0167         samplingWeight = 167.3;
0168       else {
0169         // ONLY protection against summing HO, HF simhits
0170         return 0.;
0171       }
0172     }
0173
0174     return samplingWeight * hit->energy();
0175   }
0176
0177   GlobalPoint getGpos(const CaloGeometry* geo, HBHERecHitCollection::const_iterator hit, bool) {
0178     DetId detId(hit->id());
0179     return (static_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(detId)))->getPosition(detId);
0180   }
0181
0182   GlobalPoint getGpos(const CaloGeometry* geo, edm::PCaloHitContainer::const_iterator hit, bool) {
0183     DetId detId(hit->id());
0184     GlobalPoint point = (detId.det() == DetId::Hcal)
0185                             ? (static_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(detId)))->getPosition(detId)
0186                             : geo->getPosition(detId);
0187     return point;
0188   }
0189
0190   GlobalPoint getGpos(const CaloGeometry* geo, EcalRecHitCollection::const_iterator hit, bool) {
0191     // Not tested for EcalRecHits!!
0192     if (hit->id().subdetId() == EcalEndcap) {
0193       EEDetId EEid = EEDetId(hit->id());
0194       return geo->getPosition(EEid);
0195     } else {  // (hit->id().subdetId() == EcalBarrel)
0196       EBDetId EBid = EBDetId(hit->id());
0197       return geo->getPosition(EBid);
0198     }
0199   }
0200
0201   double getRawEnergy(HBHERecHitCollection::const_iterator hit, int useRaw) {
0202     double energy = ((useRaw == 1) ? hit->eraw() : ((useRaw == 2) ? hit->eaux() : hit->energy()));
0203     return energy;
0204   }
0205
0206   double getRawEnergy(EcalRecHitCollection::const_iterator hit, int) { return hit->energy(); }
0207
0208   double getRawEnergy(edm::PCaloHitContainer::const_iterator hit, int) { return hit->energy(); }
0209 }  // namespace spr