CandUtils/src/EventShapeVariables.cc

0001 #include "PhysicsTools/CandUtils/interface/EventShapeVariables.h"
0002
0003 #include "TMath.h"
0004
0005 /// constructor from reco::Candidates
0006 EventShapeVariables::EventShapeVariables(const edm::View<reco::Candidate>& inputVectors) : eigenVectors_(3, 3) {
0007   inputVectors_.reserve(inputVectors.size());
0008   for (const auto& vec : inputVectors) {
0009     inputVectors_.push_back(math::XYZVector(vec.px(), vec.py(), vec.pz()));
0010   }
0011   //default values
0012   set_r(2.);
0013   setFWmax(10);
0014 }
0015
0016 /// constructor from XYZ coordinates
0017 EventShapeVariables::EventShapeVariables(const std::vector<math::XYZVector>& inputVectors)
0018     : inputVectors_(inputVectors), eigenVectors_(3, 3) {
0019   //default values
0020   set_r(2.);
0021   setFWmax(10);
0022 }
0023
0024 /// constructor from rho eta phi coordinates
0025 EventShapeVariables::EventShapeVariables(const std::vector<math::RhoEtaPhiVector>& inputVectors) : eigenVectors_(3, 3) {
0026   inputVectors_.reserve(inputVectors.size());
0027   for (const auto& vec : inputVectors) {
0028     inputVectors_.push_back(math::XYZVector(vec.x(), vec.y(), vec.z()));
0029   }
0030   //default values
0031   set_r(2.);
0032   setFWmax(10);
0033 }
0034
0035 /// constructor from r theta phi coordinates
0036 EventShapeVariables::EventShapeVariables(const std::vector<math::RThetaPhiVector>& inputVectors) : eigenVectors_(3, 3) {
0037   inputVectors_.reserve(inputVectors.size());
0038   for (const auto& vec : inputVectors) {
0039     inputVectors_.push_back(math::XYZVector(vec.x(), vec.y(), vec.z()));
0040   }
0041   //default values
0042   set_r(2.);
0043   setFWmax(10);
0044 }
0045
0046 /// the return value is 1 for spherical events and 0 for events linear in r-phi. This function
0047 /// needs the number of steps to determine how fine the granularity of the algorithm in phi
0048 /// should be
0049 double EventShapeVariables::isotropy(const unsigned int& numberOfSteps) const {
0050   const double deltaPhi = 2 * TMath::Pi() / numberOfSteps;
0051   double phi = 0, eIn = -1., eOut = -1.;
0052   for (unsigned int i = 0; i < numberOfSteps; ++i) {
0053     phi += deltaPhi;
0054     double sum = 0;
0055     double cosphi = TMath::Cos(phi);
0056     double sinphi = TMath::Sin(phi);
0057     for (const auto& vec : inputVectors_) {
0058       // sum over inner product of unit vectors and momenta
0059       sum += TMath::Abs(cosphi * vec.x() + sinphi * vec.y());
0060     }
0061     if (eOut < 0. || sum < eOut)
0062       eOut = sum;
0063     if (eIn < 0. || sum > eIn)
0064       eIn = sum;
0065   }
0066   return (eIn - eOut) / eIn;
0067 }
0068
0069 /// the return value is 1 for spherical and 0 linear events in r-phi. This function needs the
0070 /// number of steps to determine how fine the granularity of the algorithm in phi should be
0071 double EventShapeVariables::circularity(const unsigned int& numberOfSteps) const {
0072   const double deltaPhi = 2 * TMath::Pi() / numberOfSteps;
0073   double circularity = -1, phi = 0, area = 0;
0074   for (const auto& vec : inputVectors_) {
0075     area += TMath::Sqrt(vec.x() * vec.x() + vec.y() * vec.y());
0076   }
0077   for (unsigned int i = 0; i < numberOfSteps; ++i) {
0078     phi += deltaPhi;
0079     double sum = 0, tmp = 0.;
0080     double cosphi = TMath::Cos(phi);
0081     double sinphi = TMath::Sin(phi);
0082     for (const auto& vec : inputVectors_) {
0083       sum += TMath::Abs(cosphi * vec.x() + sinphi * vec.y());
0084     }
0085     tmp = TMath::Pi() / 2 * sum / area;
0086     if (circularity < 0 || tmp < circularity) {
0087       circularity = tmp;
0088     }
0089   }
0090   return circularity;
0091 }
0092
0093 /// set exponent for computation of momentum tensor and related products
0094 void EventShapeVariables::set_r(double r) {
0095   r_ = r;
0096   /// invalidate previous cached computations
0097   tensors_computed_ = false;
0098   eigenValues_ = std::vector<double>(3, 0);
0099   eigenValuesNoNorm_ = std::vector<double>(3, 0);
0100 }
0101
0102 /// helper function to fill the 3 dimensional momentum tensor from the inputVectors where needed
0103 /// also fill the 3 dimensional vectors of eigen-values and eigen-vectors;
0104 /// the largest (smallest) eigen-value is stored at index position 0 (2)
0105 void EventShapeVariables::compTensorsAndVectors() {
0106   if (tensors_computed_)
0107     return;
0108
0109   if (inputVectors_.size() < 2) {
0110     tensors_computed_ = true;
0111     return;
0112   }
0113
0114   TMatrixDSym momentumTensor(3);
0115   momentumTensor.Zero();
0116
0117   // fill momentumTensor from inputVectors
0118   double norm = 0.;
0119   for (const auto& vec : inputVectors_) {
0120     double p2 = vec.Dot(vec);
0121     double pR = (r_ == 2.) ? p2 : TMath::Power(p2, 0.5 * r_);
0122     norm += pR;
0123     double pRminus2 = (r_ == 2.) ? 1. : TMath::Power(p2, 0.5 * r_ - 1.);
0124     momentumTensor(0, 0) += pRminus2 * vec.x() * vec.x();
0125     momentumTensor(0, 1) += pRminus2 * vec.x() * vec.y();
0126     momentumTensor(0, 2) += pRminus2 * vec.x() * vec.z();
0127     momentumTensor(1, 0) += pRminus2 * vec.y() * vec.x();
0128     momentumTensor(1, 1) += pRminus2 * vec.y() * vec.y();
0129     momentumTensor(1, 2) += pRminus2 * vec.y() * vec.z();
0130     momentumTensor(2, 0) += pRminus2 * vec.z() * vec.x();
0131     momentumTensor(2, 1) += pRminus2 * vec.z() * vec.y();
0132     momentumTensor(2, 2) += pRminus2 * vec.z() * vec.z();
0133   }
0134
0135   if (momentumTensor.IsSymmetric() && (momentumTensor.NonZeros() != 0)) {
0136     momentumTensor.EigenVectors(eigenValuesNoNormTmp_);
0137   }
0138   eigenValuesNoNorm_[0] = eigenValuesNoNormTmp_(0);
0139   eigenValuesNoNorm_[1] = eigenValuesNoNormTmp_(1);
0140   eigenValuesNoNorm_[2] = eigenValuesNoNormTmp_(2);
0141
0142   // momentumTensor normalized to determinant 1
0143   momentumTensor *= (1. / norm);
0144
0145   // now get eigens
0146   if (momentumTensor.IsSymmetric() && (momentumTensor.NonZeros() != 0)) {
0147     eigenVectors_ = momentumTensor.EigenVectors(eigenValuesTmp_);
0148   }
0149   eigenValues_[0] = eigenValuesTmp_(0);
0150   eigenValues_[1] = eigenValuesTmp_(1);
0151   eigenValues_[2] = eigenValuesTmp_(2);
0152
0153   tensors_computed_ = true;
0154 }
0155
0156 /// 1.5*(q1+q2) where q0>=q1>=q2>=0 are the eigenvalues of the momentum tensor sum{p_j[a]*p_j[b]}/sum{p_j**2}
0157 /// normalized to 1. Return values are 1 for spherical, 3/4 for plane and 0 for linear events
0158 double EventShapeVariables::sphericity() {
0159   if (!tensors_computed_)
0160     compTensorsAndVectors();
0161   return 1.5 * (eigenValues_[1] + eigenValues_[2]);
0162 }
0163
0164 /// 1.5*q2 where q0>=q1>=q2>=0 are the eigenvalues of the momentum tensor sum{p_j[a]*p_j[b]}/sum{p_j**2}
0165 /// normalized to 1. Return values are 0.5 for spherical and 0 for plane and linear events
0166 double EventShapeVariables::aplanarity() {
0167   if (!tensors_computed_)
0168     compTensorsAndVectors();
0169   return 1.5 * eigenValues_[2];
0170 }
0171
0172 /// 3.*(q0*q1+q0*q2+q1*q2) where q0>=q1>=q2>=0 are the eigenvalues of the momentum tensor sum{p_j[a]*p_j[b]}/sum{p_j**2}
0173 /// normalized to 1. Return value is between 0 and 1
0174 /// and measures the 3-jet structure of the event (C vanishes for a "perfect" 2-jet event)
0175 double EventShapeVariables::C() {
0176   if (!tensors_computed_)
0177     compTensorsAndVectors();
0178   return 3. *
0179          (eigenValues_[0] * eigenValues_[1] + eigenValues_[0] * eigenValues_[2] + eigenValues_[1] * eigenValues_[2]);
0180 }
0181
0182 /// 27.*(q0*q1*q2) where q0>=q1>=q2>=0 are the eigenvalues of the momemtum tensor sum{p_j[a]*p_j[b]}/sum{p_j**2}
0183 /// normalized to 1. Return value is between 0 and 1
0184 /// and measures the 4-jet structure of the event (D vanishes for a planar event)
0185 double EventShapeVariables::D() {
0186   if (!tensors_computed_)
0187     compTensorsAndVectors();
0188   return 27. * eigenValues_[0] * eigenValues_[1] * eigenValues_[2];
0189 }
0190
0191 //========================================================================================================
0192
0193 /// set number of Fox-Wolfram moments to compute
0194 void EventShapeVariables::setFWmax(unsigned m) {
0195   fwmom_maxl_ = m;
0196   fwmom_computed_ = false;
0197   fwmom_ = std::vector<double>(fwmom_maxl_, 0.);
0198 }
0199
0200 double EventShapeVariables::getFWmoment(unsigned l) {
0201   if (l > fwmom_maxl_)
0202     return 0.;
0203
0204   if (!fwmom_computed_)
0205     computeFWmoments();
0206
0207   return fwmom_[l];
0208
0209 }  // getFWmoment
0210
0211 const std::vector<double>& EventShapeVariables::getFWmoments() {
0212   if (!fwmom_computed_)
0213     computeFWmoments();
0214
0215   return fwmom_;
0216 }
0217
0218 void EventShapeVariables::computeFWmoments() {
0219   if (fwmom_computed_)
0220     return;
0221
0222   double esum_total(0.);
0223   for (unsigned int i = 0; i < inputVectors_.size(); i++) {
0224     esum_total += inputVectors_[i].R();
0225   }  // i
0226   double esum_total_sq = esum_total * esum_total;
0227
0228   for (unsigned int i = 0; i < inputVectors_.size(); i++) {
0229     double p_i = inputVectors_[i].R();
0230     if (p_i <= 0)
0231       continue;
0232
0233     for (unsigned int j = 0; j <= i; j++) {
0234       double p_j = inputVectors_[j].R();
0235       if (p_j <= 0)
0236         continue;
0237
0238       /// reduce computation by exploiting symmetry:
0239       /// all off-diagonal elements appear twice in the sum
0240       int symmetry_factor = 2;
0241       if (j == i)
0242         symmetry_factor = 1;
0243       double p_ij = p_i * p_j;
0244       double cosTheta = inputVectors_[i].Dot(inputVectors_[j]) / (p_ij);
0245       double pi_pj_over_etot2 = p_ij / esum_total_sq;
0246
0247       /// compute higher legendre polynomials recursively
0248       /// need to keep track of two previous values
0249       double Pn1 = 0;
0250       double Pn2 = 0;
0251       for (unsigned n = 0; n < fwmom_maxl_; n++) {
0252         /// initial cases
0253         if (n == 0) {
0254           Pn2 = pi_pj_over_etot2;
0255           fwmom_[0] += Pn2 * symmetry_factor;
0256         } else if (n == 1) {
0257           Pn1 = pi_pj_over_etot2 * cosTheta;
0258           fwmom_[1] += Pn1 * symmetry_factor;
0259         } else {
0260           double Pn = ((2 * n - 1) * cosTheta * Pn1 - (n - 1) * Pn2) / n;
0261           fwmom_[n] += Pn * symmetry_factor;
0262           /// store new value
0263           Pn2 = Pn1;
0264           Pn1 = Pn;
0265         }
0266       }
0267
0268     }  // j
0269   }    // i
0270
0271   fwmom_computed_ = true;
0272
0273 }  // computeFWmoments
0274
0275 //========================================================================================================