File indexing completed on 2023-10-25 09:39:14
0001
0002
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
0012 #include "DataFormats/JetReco/interface/Jet.h"
0013
0014 using namespace reco;
0015
0016 namespace {
0017
0018
0019
0020 class CaloPoint {
0021 public:
0022 static const double depth;
0023 static const double R_BARREL;
0024 static const double R_BARREL2;
0025 static const double Z_ENDCAP;
0026 static const double R_FORWARD;
0027 static const double R_FORWARD2;
0028 static const double Z_FORWARD;
0029 static const double Z_BIG;
0030 };
0031
0032
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.;
0037 const double CaloPoint::R_FORWARD = Z_ENDCAP / std::sqrt(std::cosh(3.) * std::cosh(3.) - 1.);
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
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;
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) {
0058 mR = R_BARREL;
0059 mZ = fZ + R_BARREL * tanTheta;
0060 } else {
0061 double zRef = Z_BIG;
0062 if (rEndcap > R_FORWARD)
0063 zRef = Z_ENDCAP;
0064 else if (std::fabs(fEta) < ETA_MAX)
0065 zRef = Z_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
0091
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
0098
0099
0100 int side = dir.z() < -1e-9 ? -1 : dir.z() > 1e-9 ? +1 : 0;
0101
0102 double dirR = dir.Rho();
0103
0104
0105 double dirUnit[2] = {dir.x() / dirR, dir.y() / dirR};
0106
0107
0108
0109
0110 double vtxLong = dirUnit[0] * vertex.x() + dirUnit[1] * vertex.y();
0111
0112
0113 double tIP = dirUnit[0] * vertex.y() - dirUnit[1] * vertex.x();
0114
0115
0116 double r, z;
0117
0118 if (side) {
0119 double slope = dirR / dir.z();
0120
0121
0122 r = vtxLong + slope * (side * Z_ENDCAP - vertex.z());
0123 double r2 = sqr(r) + sqr(tIP);
0124
0125 if (r2 < R_FORWARD2) {
0126
0127 r = vtxLong + slope * (side * Z_FORWARD - vertex.z());
0128 z = side * Z_FORWARD;
0129 } else if (r2 < R_BARREL2) {
0130
0131 z = side * Z_ENDCAP;
0132 } else {
0133
0134 side = 0;
0135 }
0136 }
0137
0138 if (!side) {
0139
0140 double slope = dir.z() / dirR;
0141 r = std::sqrt(R_BARREL2 - sqr(tIP));
0142 z = vertex.z() + slope * (r - vtxLong);
0143 }
0144
0145
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 }
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
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
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
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
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
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
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
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
0311
0312 float Jet::physicsEta(float fZVertex, float fDetectorEta) {
0313 CaloPointZ refPoint(0., fDetectorEta);
0314 return refPoint.etaReference(fZVertex);
0315 }
0316
0317
0318
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());
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());
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 }
0381
0382 float ptD_value = (sum_pt > 0.) ? sqrt(sum_pt2 / (sum_pt * sum_pt)) : 0.;
0383
0384 return ptD_value;
0385
0386 }
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 }
0406
0407 float rmsCand_value = (sum_pt2 > 0.) ? sum_pt2deltaR2 / sum_pt2 : 0.;
0408
0409 return rmsCand_value;
0410
0411 }
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);
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; }