HcalCalibAlgos/interface/CommonUsefulStuff.h

0001 #ifndef Calibration_HcalCalibALgos_CommonUsefulStuff_h
0002 #define Calibration_HcalCalibALgos_CommonUsefulStuff_h
0003
0004 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0005 #include "DataFormats/DetId/interface/DetId.h"
0006 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0007 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0008 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0009
0010 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0011 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0012
0013 #include "FWCore/Utilities/interface/Exception.h"
0014
0015 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0016 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0017 #include "CommonTools/UtilAlgos/interface/DeltaR.h"
0018
0019 /* getDist functions by Jim:
0020 /uscms/home/jhirsch/IsoTracks_314/CMSSW_3_1_4/src/JetMETCorrections/IsolatedParticles/interface/FindCaloHit.icc
0021 */
0022
0023 struct MaxHit_struct {
0024   int iphihitm;
0025   int ietahitm;
0026   int depthhit;
0027   float hitenergy;
0028   float dr;
0029   GlobalPoint posMax;
0030   MaxHit_struct() : iphihitm(0), ietahitm(0), depthhit(0), hitenergy(-100), dr(0) {}
0031 };
0032
0033 inline double getDistInPlaneSimple(const GlobalPoint caloPoint, const GlobalPoint rechitPoint) {
0034   // Simplified version of getDistInPlane
0035   // Assume track direction is origin -> point of hcal intersection
0036
0037   const GlobalVector caloIntersectVector(caloPoint.x(), caloPoint.y(), caloPoint.z());
0038
0039   const GlobalVector caloIntersectUnitVector = caloIntersectVector.unit();
0040
0041   const GlobalVector rechitVector(rechitPoint.x(), rechitPoint.y(), rechitPoint.z());
0042
0043   const GlobalVector rechitUnitVector = rechitVector.unit();
0044   double dotprod = caloIntersectUnitVector.dot(rechitUnitVector);
0045   double rechitdist = caloIntersectVector.mag() / dotprod;
0046
0047   const GlobalVector effectiveRechitVector = rechitdist * rechitUnitVector;
0048   const GlobalPoint effectiveRechitPoint(
0049       effectiveRechitVector.x(), effectiveRechitVector.y(), effectiveRechitVector.z());
0050
0051   GlobalVector distance_vector = effectiveRechitPoint - caloPoint;
0052
0053   if (dotprod > 0.) {
0054     return distance_vector.mag();
0055   } else {
0056     return 999999.;
0057   }
0058 }
0059
0060 inline double getDistInPlaneTrackDir(const GlobalPoint caloPoint,
0061                                      const GlobalVector caloVector,
0062                                      const GlobalPoint rechitPoint) {
0063   // Simplified version of getDistInPlane : no cone "within" Hcal, but
0064   // don't assume track direction is origin -> point of hcal
0065   // intersection.
0066   const GlobalVector caloIntersectVector(caloPoint.x(), caloPoint.y(),
0067                                          caloPoint.z());  //p
0068
0069   const GlobalVector caloUnitVector = caloVector.unit();
0070   const GlobalVector rechitVector(rechitPoint.x(), rechitPoint.y(), rechitPoint.z());
0071   const GlobalVector rechitUnitVector = rechitVector.unit();
0072   double dotprod_denominator = caloUnitVector.dot(rechitUnitVector);
0073   double dotprod_numerator = caloUnitVector.dot(caloIntersectVector);
0074   double rechitdist = dotprod_numerator / dotprod_denominator;
0075   const GlobalVector effectiveRechitVector = rechitdist * rechitUnitVector;
0076   const GlobalPoint effectiveRechitPoint(
0077       effectiveRechitVector.x(), effectiveRechitVector.y(), effectiveRechitVector.z());
0078   GlobalVector distance_vector = effectiveRechitPoint - caloPoint;
0079   if (dotprod_denominator > 0. && dotprod_numerator > 0.) {
0080     return distance_vector.mag();
0081   } else {
0082     return 999999.;
0083   }
0084 }
0085
0086 inline double getDistInPlane(const GlobalVector trackDirection,
0087                              const GlobalPoint caloPoint,
0088                              const GlobalPoint rechitPoint,
0089                              double coneHeight) {
0090   // The iso track candidate hits the Calo (Ecal or Hcal) at "caloPoint"
0091   // with direction "trackDirection".
0092
0093   // "rechitPoint" is the position of the rechit.  We only care about
0094   // the direction of the rechit.
0095
0096   // Consider the rechitPoint as characterized by angles theta and phi
0097   // wrt the origin which points at the calo cell of the rechit.  In
0098   // some sense the distance along the line theta/phi is arbitrary. A
0099   // simplified choice might be to put the rechit at the surface of
0100   // the Hcal.  Our approach is to see whether/where this line
0101   // intersects a cone of height "coneHeight" with vertex at caloPoint
0102   // aligned with trackDirection.
0103   // To that end, this function returns the distance between the
0104   // center of the base of the cone and the intersection of the rechit
0105   // line and the plane that contains the base of the cone.  This
0106   // distance is compared with the radius of the cone outside this
0107   // function.
0108
0109   // Make vector of length cone height along track direction
0110   const GlobalVector heightVector = trackDirection * coneHeight;
0111
0112   // Make vector from origin to point where iso track intersects the
0113   // calorimeter.
0114   const GlobalVector caloIntersectVector(caloPoint.x(), caloPoint.y(), caloPoint.z());
0115
0116   // Make vector from origin to point in center of base of cone
0117   const GlobalVector coneBaseVector = heightVector + caloIntersectVector;
0118
0119   // Make point in center of base of cone
0120   const GlobalPoint coneBasePoint(coneBaseVector.x(), coneBaseVector.y(), coneBaseVector.z());
0121
0122   // Make unit vector pointing at rechit.
0123   const GlobalVector rechitVector(rechitPoint.x(), rechitPoint.y(), rechitPoint.z());
0124   const GlobalVector rechitDirection = rechitVector.unit();
0125
0126   // Find distance "r" along "rechit line" (with angles theta2 and
0127   // phi2) where line intersects plane defined by base of cone.
0128
0129   // Definition plane of that contains base of cone:
0130   // trackDirection.x() (x - coneBaseVector.x()) +
0131   // trackDirection.y() (y - coneBaseVector.y()) +
0132   // trackDirection.z() (z - coneBaseVector.z()) = 0
0133
0134   // Point P_{rh} where rechit line intersects plane:
0135   // (rechitdist sin(theta2) cos(phi2),
0136   //  rechitdist sin(theta2) cos(phi2),
0137   //  rechitdist cos(theta2))
0138
0139   // Substitute P_{rh} into equation for plane and solve for rechitdist.
0140   // rechitDist turns out to be the ratio of dot products:
0141
0142   double rechitdist = trackDirection.dot(coneBaseVector) / trackDirection.dot(rechitDirection);
0143
0144   // Now find distance between point at center of cone base and point
0145   // where the "rechit line" intersects the plane defined by the base
0146   // of the cone; i.e. the effectiveRecHitPoint.
0147   const GlobalVector effectiveRechitVector = rechitdist * rechitDirection;
0148   const GlobalPoint effectiveRechitPoint(
0149       effectiveRechitVector.x(), effectiveRechitVector.y(), effectiveRechitVector.z());
0150
0151   GlobalVector distance_vector = effectiveRechitPoint - coneBasePoint;
0152   return distance_vector.mag();
0153 }
0154
0155 /*  Function to calculate Ecal Energy in Cone (given in cm) */
0156 inline double ecalEnergyInCone(const GlobalPoint center,
0157                                double radius,
0158                                const EcalRecHitCollection ecalCol,
0159                                const CaloGeometry* geo) {
0160   double eECALcone = 0;
0161   std::vector<int> usedHitsEcal;
0162   usedHitsEcal.clear();
0163
0164   for (std::vector<EcalRecHit>::const_iterator ehit = ecalCol.begin(); ehit != ecalCol.end(); ehit++) {
0165     //This is a precaution for the case when hitCollection contains duplicats.
0166     bool hitIsUsed = false;
0167     int hitHashedIndex = -10000;
0168     if (ehit->id().subdetId() == EcalBarrel) {
0169       EBDetId did(ehit->id());
0170       hitHashedIndex = did.hashedIndex();
0171     }
0172     if (ehit->id().subdetId() == EcalEndcap) {
0173       EEDetId did(ehit->id());
0174       hitHashedIndex = did.hashedIndex();
0175     }
0176     for (uint32_t i = 0; i < usedHitsEcal.size(); i++) {
0177       if (usedHitsEcal[i] == hitHashedIndex)
0178         hitIsUsed = true;
0179     }
0180     if (hitIsUsed)
0181       continue;
0182     usedHitsEcal.push_back(hitHashedIndex);
0183     // -----------------------------------------------
0184
0185     const GlobalPoint& pos = geo->getPosition((*ehit).detid());
0186
0187     if (getDistInPlaneSimple(center, pos) < radius) {
0188       eECALcone += ehit->energy();
0189     }
0190   }
0191   return eECALcone;
0192 }
0193
0194 /*  This is another version of ecalEnergy calculation using the getDistInPlaneTrackDir()  */
0195 inline double ecalEnergyInCone(const GlobalVector trackMom,
0196                                const GlobalPoint center,
0197                                double radius,
0198                                const EcalRecHitCollection ecalCol,
0199                                const CaloGeometry* geo) {
0200   double eECALcone = 0;
0201   std::vector<int> usedHitsEcal;
0202   usedHitsEcal.clear();
0203   for (std::vector<EcalRecHit>::const_iterator ehit = ecalCol.begin(); ehit != ecalCol.end(); ehit++) {
0204     //This is a precaution for the case when hitCollection contains duplicats.
0205     bool hitIsUsed = false;
0206     int hitHashedIndex = -10000;
0207     if (ehit->id().subdetId() == EcalBarrel) {
0208       EBDetId did(ehit->id());
0209       hitHashedIndex = did.hashedIndex();
0210     }
0211     if (ehit->id().subdetId() == EcalEndcap) {
0212       EEDetId did(ehit->id());
0213       hitHashedIndex = did.hashedIndex();
0214     }
0215     for (uint32_t i = 0; i < usedHitsEcal.size(); i++) {
0216       if (usedHitsEcal[i] == hitHashedIndex)
0217         hitIsUsed = true;
0218     }
0219     if (hitIsUsed)
0220       continue;
0221     usedHitsEcal.push_back(hitHashedIndex);
0222     // -----------------------------------------------
0223
0224     const GlobalPoint& pos = geo->getPosition((*ehit).detid());
0225
0226     if (getDistInPlaneTrackDir(center, trackMom, pos) < radius) {
0227       eECALcone += ehit->energy();
0228     }
0229   }
0230   return eECALcone;
0231 }
0232
0233 #endif