JetReco/src/Jet.cc

0001 // Jet.cc
0002 // Fedor Ratnikov, UMd
0003
0004 #include <sstream>
0005 #include <cmath>
0006
0007 #include "DataFormats/Math/interface/deltaR.h"
0008 #include "DataFormats/Math/interface/deltaPhi.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010
0011 //Own header file
0012 #include "DataFormats/JetReco/interface/Jet.h"
0013
0014 using namespace reco;
0015
0016 namespace {
0017   // approximate simple CALO geometry
0018   // abstract baseclass for geometry.
0019
0020   class CaloPoint {
0021   public:
0022     static const double depth;  // one for all relative depth of the reference point between ECAL begin and HCAL end
0023     static const double R_BARREL;
0024     static const double R_BARREL2;
0025     static const double Z_ENDCAP;   // 1/2(EEz+HEz)
0026     static const double R_FORWARD;  // eta=3
0027     static const double R_FORWARD2;
0028     static const double Z_FORWARD;
0029     static const double Z_BIG;
0030   };
0031
0032   // one for all relative depth of the reference point between ECAL begin and HCAL end
0033   const double CaloPoint::depth = 0.1;
0034   const double CaloPoint::R_BARREL = (1. - depth) * 143. + depth * 407.;
0035   const double CaloPoint::R_BARREL2 = R_BARREL * R_BARREL;
0036   const double CaloPoint::Z_ENDCAP = (1. - depth) * 320. + depth * 568.;                         // 1/2(EEz+HEz)
0037   const double CaloPoint::R_FORWARD = Z_ENDCAP / std::sqrt(std::cosh(3.) * std::cosh(3.) - 1.);  // eta=3
0038   const double CaloPoint::R_FORWARD2 = R_FORWARD * R_FORWARD;
0039   const double CaloPoint::Z_FORWARD = 1100. + depth * 165.;
0040   const double CaloPoint::Z_BIG = 1.e5;
0041
0042   //old zvertex only implementation:
0043   class CaloPointZ : private CaloPoint {
0044   public:
0045     CaloPointZ(double fZ, double fEta) {
0046       static const double ETA_MAX = 5.2;
0047
0048       if (fZ > Z_ENDCAP)
0049         fZ = Z_ENDCAP - 1.;
0050       if (fZ < -Z_ENDCAP)
0051         fZ = -Z_ENDCAP + 1;  // sanity check
0052
0053       double tanThetaAbs = std::sqrt(std::cosh(fEta) * std::cosh(fEta) - 1.);
0054       double tanTheta = fEta >= 0 ? tanThetaAbs : -tanThetaAbs;
0055
0056       double rEndcap = tanTheta == 0 ? 1.e10 : fEta > 0 ? (Z_ENDCAP - fZ) / tanTheta : (-Z_ENDCAP - fZ) / tanTheta;
0057       if (rEndcap > R_BARREL) {  // barrel
0058         mR = R_BARREL;
0059         mZ = fZ + R_BARREL * tanTheta;
0060       } else {
0061         double zRef = Z_BIG;  // very forward;
0062         if (rEndcap > R_FORWARD)
0063           zRef = Z_ENDCAP;  // endcap
0064         else if (std::fabs(fEta) < ETA_MAX)
0065           zRef = Z_FORWARD;  // forward
0066
0067         mZ = fEta > 0 ? zRef : -zRef;
0068         mR = std::fabs((mZ - fZ) / tanTheta);
0069       }
0070     }
0071     double etaReference(double fZ) {
0072       Jet::Point p(r(), 0., z() - fZ);
0073       return p.eta();
0074     }
0075
0076     double thetaReference(double fZ) {
0077       Jet::Point p(r(), 0., z() - fZ);
0078       return p.theta();
0079     }
0080
0081     double z() const { return mZ; }
0082     double r() const { return mR; }
0083
0084   private:
0085     CaloPointZ() {}
0086     double mZ;
0087     double mR;
0088   };
0089
0090   //new implementation to derive CaloPoint for free 3d vertex.
0091   //code provided thanks to Christophe Saout
0092   template <typename Point>
0093   class CaloPoint3D : private CaloPoint {
0094   public:
0095     template <typename Vector, typename Point2>
0096     CaloPoint3D(const Point2 &vertex, const Vector &dir) {
0097       // note: no sanity checks here, make sure vertex is inside the detector!
0098
0099       // check if positive or negative (or none) endcap should be tested
0100       int side = dir.z() < -1e-9 ? -1 : dir.z() > 1e-9 ? +1 : 0;
0101
0102       double dirR = dir.Rho();
0103
0104       // normalized direction in x-y plane
0105       double dirUnit[2] = {dir.x() / dirR, dir.y() / dirR};
0106
0107       // rotate the vertex into a coordinate system where dir lies along x
0108
0109       // vtxLong is the longitudinal coordinate of the vertex wrt/ dir
0110       double vtxLong = dirUnit[0] * vertex.x() + dirUnit[1] * vertex.y();
0111
0112       // tIP is the (signed) transverse impact parameter
0113       double tIP = dirUnit[0] * vertex.y() - dirUnit[1] * vertex.x();
0114
0115       // r and z coordinate
0116       double r, z;
0117
0118       if (side) {
0119         double slope = dirR / dir.z();
0120
0121         // check extrapolation to endcap
0122         r = vtxLong + slope * (side * Z_ENDCAP - vertex.z());
0123         double r2 = sqr(r) + sqr(tIP);
0124
0125         if (r2 < R_FORWARD2) {
0126           // we are in the forward calorimeter, recompute
0127           r = vtxLong + slope * (side * Z_FORWARD - vertex.z());
0128           z = side * Z_FORWARD;
0129         } else if (r2 < R_BARREL2) {
0130           // we are in the endcap
0131           z = side * Z_ENDCAP;
0132         } else {
0133           // we are in the barrel, do the intersection below
0134           side = 0;
0135         }
0136       }
0137
0138       if (!side) {
0139         // we are in the barrel
0140         double slope = dir.z() / dirR;
0141         r = std::sqrt(R_BARREL2 - sqr(tIP));
0142         z = vertex.z() + slope * (r - vtxLong);
0143       }
0144
0145       // rotate (r, tIP, z) back into original x-y coordinate system
0146       point = Point(dirUnit[0] * r - dirUnit[1] * tIP, dirUnit[1] * r + dirUnit[0] * tIP, z);
0147     }
0148
0149     const Point &caloPoint() const { return point; }
0150
0151   private:
0152     template <typename T>
0153     static inline T sqr(const T &value) {
0154       return value * value;
0155     }
0156
0157     Point point;
0158   };
0159
0160 }  // namespace
0161
0162 Jet::Jet(const LorentzVector &fP4, const Point &fVertex, const Constituents &fConstituents)
0163     : CompositePtrCandidate(0, fP4, fVertex), mJetArea(0), mPileupEnergy(0), mPassNumber(0), mIsWeighted(false) {
0164   for (unsigned i = 0; i < fConstituents.size(); i++) {
0165     addDaughter(fConstituents[i]);
0166   }
0167 }
0168
0169 Jet::Jet(const LorentzVector &fP4, const Point &fVertex)
0170     : CompositePtrCandidate(0, fP4, fVertex), mJetArea(0), mPileupEnergy(0), mPassNumber(0), mIsWeighted(false) {}
0171
0172 /// eta-phi statistics
0173 Jet::EtaPhiMoments Jet::etaPhiStatistics() const {
0174   std::vector<const Candidate *> towers = getJetConstituentsQuick();
0175   double phiRef = phi();
0176   double sumw = 0;
0177   double sumEta = 0;
0178   double sumEta2 = 0;
0179   double sumPhi = 0;
0180   double sumPhi2 = 0;
0181   double sumEtaPhi = 0;
0182   int i = towers.size();
0183   while (--i >= 0) {
0184     double eta = towers[i]->eta();
0185     double phi = deltaPhi(towers[i]->phi(), phiRef);
0186     double weight = towers[i]->et();
0187     sumw += weight;
0188     sumEta += eta * weight;
0189     sumEta2 += eta * eta * weight;
0190     sumPhi += phi * weight;
0191     sumPhi2 += phi * phi * weight;
0192     sumEtaPhi += eta * phi * weight;
0193   }
0194   Jet::EtaPhiMoments result;
0195   if (sumw > 0) {
0196     result.etaMean = sumEta / sumw;
0197     result.phiMean = deltaPhi(phiRef + sumPhi, 0.);
0198     result.etaEtaMoment = (sumEta2 - sumEta * sumEta / sumw) / sumw;
0199     result.phiPhiMoment = (sumPhi2 - sumPhi * sumPhi / sumw) / sumw;
0200     result.etaPhiMoment = (sumEtaPhi - sumEta * sumPhi / sumw) / sumw;
0201   } else {
0202     result.etaMean = 0;
0203     result.phiMean = 0;
0204     result.etaEtaMoment = 0;
0205     result.phiPhiMoment = 0;
0206     result.etaPhiMoment = 0;
0207   }
0208   return result;
0209 }
0210
0211 /// eta-eta second moment
0212 float Jet::etaetaMoment() const {
0213   std::vector<const Candidate *> towers = getJetConstituentsQuick();
0214   double sumw = 0;
0215   double sum = 0;
0216   double sum2 = 0;
0217   int i = towers.size();
0218   while (--i >= 0) {
0219     double value = towers[i]->eta();
0220     double weight = towers[i]->et();
0221     sumw += weight;
0222     sum += value * weight;
0223     sum2 += value * value * weight;
0224   }
0225   return sumw > 0 ? (sum2 - sum * sum / sumw) / sumw : 0;
0226 }
0227
0228 /// phi-phi second moment
0229 float Jet::phiphiMoment() const {
0230   double phiRef = phi();
0231   std::vector<const Candidate *> towers = getJetConstituentsQuick();
0232   double sumw = 0;
0233   double sum = 0;
0234   double sum2 = 0;
0235   int i = towers.size();
0236   while (--i >= 0) {
0237     double value = deltaPhi(towers[i]->phi(), phiRef);
0238     double weight = towers[i]->et();
0239     sumw += weight;
0240     sum += value * weight;
0241     sum2 += value * value * weight;
0242   }
0243   return sumw > 0 ? (sum2 - sum * sum / sumw) / sumw : 0;
0244 }
0245
0246 /// eta-phi second moment
0247 float Jet::etaphiMoment() const {
0248   double phiRef = phi();
0249   std::vector<const Candidate *> towers = getJetConstituentsQuick();
0250   double sumw = 0;
0251   double sumA = 0;
0252   double sumB = 0;
0253   double sumAB = 0;
0254   int i = towers.size();
0255   while (--i >= 0) {
0256     double valueA = towers[i]->eta();
0257     double valueB = deltaPhi(towers[i]->phi(), phiRef);
0258     double weight = towers[i]->et();
0259     sumw += weight;
0260     sumA += valueA * weight;
0261     sumB += valueB * weight;
0262     sumAB += valueA * valueB * weight;
0263   }
0264   return sumw > 0 ? (sumAB - sumA * sumB / sumw) / sumw : 0;
0265 }
0266
0267 /// et in annulus between rmin and rmax around jet direction
0268 float Jet::etInAnnulus(float fRmin, float fRmax) const {
0269   float result = 0;
0270   std::vector<const Candidate *> towers = getJetConstituentsQuick();
0271   int i = towers.size();
0272   while (--i >= 0) {
0273     double r = deltaR(*this, *(towers[i]));
0274     if (r >= fRmin && r < fRmax)
0275       result += towers[i]->et();
0276   }
0277   return result;
0278 }
0279
0280 /// return # of constituent carring fraction of energy. Assume ordered towers
0281 int Jet::nCarrying(float fFraction) const {
0282   std::vector<const Candidate *> towers = getJetConstituentsQuick();
0283   if (fFraction >= 1)
0284     return towers.size();
0285   double totalEt = 0;
0286   for (unsigned i = 0; i < towers.size(); ++i)
0287     totalEt += towers[i]->et();
0288   double fractionEnergy = totalEt * fFraction;
0289   unsigned result = 0;
0290   for (; result < towers.size(); ++result) {
0291     fractionEnergy -= towers[result]->et();
0292     if (fractionEnergy <= 0)
0293       return result + 1;
0294   }
0295   return 0;
0296 }
0297
0298 /// maximum distance from jet to constituent
0299 float Jet::maxDistance() const {
0300   float result = 0;
0301   std::vector<const Candidate *> towers = getJetConstituentsQuick();
0302   for (unsigned i = 0; i < towers.size(); ++i) {
0303     float dR = deltaR(*(towers[i]), *this);
0304     if (dR > result)
0305       result = dR;
0306   }
0307   return result;
0308 }
0309
0310 /// static function to convert detector eta to physics eta
0311 // kept for backwards compatibility, use detector/physicsP4 instead!
0312 float Jet::physicsEta(float fZVertex, float fDetectorEta) {
0313   CaloPointZ refPoint(0., fDetectorEta);
0314   return refPoint.etaReference(fZVertex);
0315 }
0316
0317 /// static function to convert physics eta to detector eta
0318 // kept for backwards compatibility, use detector/physicsP4 instead!
0319 float Jet::detectorEta(float fZVertex, float fPhysicsEta) {
0320   CaloPointZ refPoint(fZVertex, fPhysicsEta);
0321   return refPoint.etaReference(0.);
0322 }
0323
0324 Candidate::LorentzVector Jet::physicsP4(const Candidate::Point &newVertex,
0325                                         const Candidate &inParticle,
0326                                         const Candidate::Point &oldVertex) {
0327   CaloPoint3D<Point> caloPoint(oldVertex, inParticle.momentum());  // Jet position in Calo.
0328   Vector physicsDir = caloPoint.caloPoint() - newVertex;
0329   double p = inParticle.momentum().r();
0330   Vector p3 = p * physicsDir.unit();
0331   LorentzVector returnVector(p3.x(), p3.y(), p3.z(), inParticle.energy());
0332   return returnVector;
0333 }
0334
0335 Candidate::LorentzVector Jet::detectorP4(const Candidate::Point &vertex, const Candidate &inParticle) {
0336   CaloPoint3D<Point> caloPoint(vertex, inParticle.momentum());  // Jet position in Calo.
0337   static const Point np(0, 0, 0);
0338   Vector detectorDir = caloPoint.caloPoint() - np;
0339   double p = inParticle.momentum().r();
0340   Vector p3 = p * detectorDir.unit();
0341   LorentzVector returnVector(p3.x(), p3.y(), p3.z(), inParticle.energy());
0342   return returnVector;
0343 }
0344
0345 Jet::Constituents Jet::getJetConstituents() const {
0346   Jet::Constituents result;
0347   for (unsigned i = 0; i < numberOfDaughters(); i++) {
0348     result.push_back(daughterPtr(i));
0349   }
0350   return result;
0351 }
0352
0353 std::vector<const Candidate *> Jet::getJetConstituentsQuick() const {
0354   std::vector<const Candidate *> result;
0355   int nDaughters = numberOfDaughters();
0356   for (int i = 0; i < nDaughters; ++i) {
0357     if (!(this->sourceCandidatePtr(i).isNonnull() && this->sourceCandidatePtr(i).isAvailable())) {
0358       edm::LogInfo("JetConstituentPointer")
0359           << "Bad pointer to jet constituent found.  Check for possible dropped particle.";
0360       continue;
0361     }
0362     result.push_back(daughter(i));
0363   }
0364   return result;
0365 }
0366
0367 float Jet::constituentPtDistribution() const {
0368   Jet::Constituents constituents = this->getJetConstituents();
0369
0370   float sum_pt2 = 0.;
0371   float sum_pt = 0.;
0372
0373   for (unsigned iConst = 0; iConst < constituents.size(); ++iConst) {
0374     float pt = constituents[iConst]->p4().Pt();
0375     float pt2 = pt * pt;
0376
0377     sum_pt += pt;
0378     sum_pt2 += pt2;
0379
0380   }  //for constituents
0381
0382   float ptD_value = (sum_pt > 0.) ? sqrt(sum_pt2 / (sum_pt * sum_pt)) : 0.;
0383
0384   return ptD_value;
0385
0386 }  //constituentPtDistribution
0387
0388 float Jet::constituentEtaPhiSpread() const {
0389   Jet::Constituents constituents = this->getJetConstituents();
0390
0391   float sum_pt2 = 0.;
0392   float sum_pt2deltaR2 = 0.;
0393
0394   for (unsigned iConst = 0; iConst < constituents.size(); ++iConst) {
0395     LorentzVector thisConstituent = constituents[iConst]->p4();
0396
0397     float pt = thisConstituent.Pt();
0398     float pt2 = pt * pt;
0399     double dR = deltaR(*this, *(constituents[iConst]));
0400     float pt2deltaR2 = pt * pt * dR * dR;
0401
0402     sum_pt2 += pt2;
0403     sum_pt2deltaR2 += pt2deltaR2;
0404
0405   }  //for constituents
0406
0407   float rmsCand_value = (sum_pt2 > 0.) ? sum_pt2deltaR2 / sum_pt2 : 0.;
0408
0409   return rmsCand_value;
0410
0411 }  //constituentEtaPhiSpread
0412
0413 std::string Jet::print() const {
0414   std::ostringstream out;
0415   out << "Jet p/px/py/pz/pt: " << p() << '/' << px() << '/' << py() << '/' << pz() << '/' << pt() << std::endl
0416       << "    eta/phi: " << eta() << '/' << phi() << std::endl
0417       << "    # of constituents: " << nConstituents() << std::endl;
0418   out << "    Constituents:" << std::endl;
0419   for (unsigned index = 0; index < numberOfDaughters(); index++) {
0420     Constituent constituent = daughterPtr(index);  // deref
0421     if (constituent.isNonnull()) {
0422       out << "      #" << index << " p/pt/eta/phi: " << constituent->p() << '/' << constituent->pt() << '/'
0423           << constituent->eta() << '/' << constituent->phi() << "    productId/index: " << constituent.id() << '/'
0424           << constituent.key() << std::endl;
0425     } else {
0426       out << "      #" << index << " constituent is not available in the event" << std::endl;
0427     }
0428   }
0429   return out.str();
0430 }
0431
0432 void Jet::scaleEnergy(double fScale) { setP4(p4() * fScale); }
0433
0434 bool Jet::isJet() const { return true; }