File indexing completed on 2024-09-07 04:37:17
0001 #include "PhysicsTools/CandUtils/interface/EventShapeVariables.h"
0002
0003 #include "TMath.h"
0004
0005
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
0012 set_r(2.);
0013 setFWmax(10);
0014 }
0015
0016
0017 EventShapeVariables::EventShapeVariables(const std::vector<math::XYZVector>& inputVectors)
0018 : inputVectors_(inputVectors), eigenVectors_(3, 3) {
0019
0020 set_r(2.);
0021 setFWmax(10);
0022 }
0023
0024
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
0031 set_r(2.);
0032 setFWmax(10);
0033 }
0034
0035
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
0042 set_r(2.);
0043 setFWmax(10);
0044 }
0045
0046
0047
0048
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
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
0070
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
0094 void EventShapeVariables::set_r(double r) {
0095 r_ = r;
0096
0097 tensors_computed_ = false;
0098 eigenValues_ = std::vector<double>(3, 0);
0099 eigenValuesNoNorm_ = std::vector<double>(3, 0);
0100 }
0101
0102
0103
0104
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
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
0143 momentumTensor *= (1. / norm);
0144
0145
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
0157
0158 double EventShapeVariables::sphericity() {
0159 if (!tensors_computed_)
0160 compTensorsAndVectors();
0161 return 1.5 * (eigenValues_[1] + eigenValues_[2]);
0162 }
0163
0164
0165
0166 double EventShapeVariables::aplanarity() {
0167 if (!tensors_computed_)
0168 compTensorsAndVectors();
0169 return 1.5 * eigenValues_[2];
0170 }
0171
0172
0173
0174
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
0183
0184
0185 double EventShapeVariables::D() {
0186 if (!tensors_computed_)
0187 compTensorsAndVectors();
0188 return 27. * eigenValues_[0] * eigenValues_[1] * eigenValues_[2];
0189 }
0190
0191
0192
0193
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 }
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 }
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
0239
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
0248
0249 double Pn1 = 0;
0250 double Pn2 = 0;
0251 for (unsigned n = 0; n < fwmom_maxl_; n++) {
0252
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
0263 Pn2 = Pn1;
0264 Pn1 = Pn;
0265 }
0266 }
0267
0268 }
0269 }
0270
0271 fwmom_computed_ = true;
0272
0273 }
0274
0275