Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:04:58 

0001 #include "DataFormats/PatCandidates/interface/ResolutionHelper.h"
0002 #include "FWCore/Utilities/interface/Exception.h"
0003 #include <cmath>
0004 #include <iostream>
0005 
0006 using namespace std;  // faster sqrt, sqrt, pow(x,2)....
0007 
0008 void pat::helper::ResolutionHelper::rescaleForKinFitter(pat::CandKinResolution::Parametrization parametrization,
0009                                                         AlgebraicSymMatrix44 &covariance,
0010                                                         const math::XYZTLorentzVector &initialP4) {
0011   double inv;
0012   switch (parametrization) {
0013     case pat::CandKinResolution::Cart:
0014     case pat::CandKinResolution::Spher:
0015       // for us parameter[3] = mass, for KinFitter parameter[3] = mass/initialP4.mass();
0016       inv = 1.0 / initialP4.mass();
0017       for (int i = 0; i < 4; i++) {
0018         covariance(3, i) *= inv;
0019       }
0020       covariance(3, 3) *= inv;
0021       break;
0022     case pat::CandKinResolution::ESpher:
0023       // for us parameter[3] = energy, for KinFitter parameter[3] = energy/initialP4.energy();
0024       inv = 1.0 / initialP4.energy();
0025       for (int i = 0; i < 4; i++) {
0026         covariance(3, i) *= inv;
0027       }
0028       covariance(3, 3) *= inv;
0029       break;
0030     default:;  // nothing to do
0031   }
0032 }
0033 
0034 double pat::helper::ResolutionHelper::getResolP(pat::CandKinResolution::Parametrization parametrization,
0035                                                 const AlgebraicSymMatrix44 &covariance,
0036                                                 const pat::CandKinResolution::LorentzVector &p4) {
0037   switch (parametrization) {
0038     // ==== CARTESIAN ====
0039     case pat::CandKinResolution::Cart:
0040     case pat::CandKinResolution::ECart:
0041     case pat::CandKinResolution::MCCart: {
0042       // p2 = px^2 + py^2 + pz^2
0043       // dp/dp_i = p_i/p ==> it is a unit vector!
0044       AlgebraicVector3 derivs(p4.X(), p4.Y(), p4.Z());
0045       derivs.Unit();
0046       return sqrt(ROOT::Math::Similarity(derivs, covariance.Sub<AlgebraicSymMatrix33>(0, 0)));
0047     }
0048     // ==== SPHERICAL (P)  ====
0049     case pat::CandKinResolution::Spher:
0050     case pat::CandKinResolution::ESpher:
0051     case pat::CandKinResolution::MCSpher:
0052       return sqrt(covariance(0, 0));
0053     // ==== SPHERICAL (1/P)  ====
0054     case pat::CandKinResolution::MCPInvSpher:
0055       return sqrt(covariance(0, 0)) * (p4.P2());
0056     // ==== HEP CYLINDRICAL (with Pt = Et, P = E) ====
0057     case pat::CandKinResolution::EtThetaPhi:
0058     case pat::CandKinResolution::EtEtaPhi:
0059       return getResolE(parametrization, covariance, p4);
0060     // ==== OTHER ====
0061     case pat::CandKinResolution::Invalid:
0062       throw cms::Exception("Invalid parametrization") << parametrization;
0063     default:
0064       throw cms::Exception("Not Implemented")
0065           << "getResolP not yet implemented for parametrization " << parametrization;
0066   }
0067 }
0068 double pat::helper::ResolutionHelper::getResolPt(pat::CandKinResolution::Parametrization parametrization,
0069                                                  const AlgebraicSymMatrix44 &covariance,
0070                                                  const pat::CandKinResolution::LorentzVector &p4) {
0071   switch (parametrization) {
0072     // ==== CARTESIAN ====
0073     case pat::CandKinResolution::Cart:
0074     case pat::CandKinResolution::ECart:
0075     case pat::CandKinResolution::MCCart: {
0076       double pti2 = 1.0 / (p4.Perp2());
0077       return sqrt((covariance(0, 0) * p4.Px() * p4.Px() + covariance(1, 1) * p4.Py() * p4.Py() +
0078                    2 * covariance(0, 1) * p4.Px() * p4.Py()) *
0079                   pti2);
0080     }
0081     // ==== SPHERICAL (P, Theta)  ====
0082     case pat::CandKinResolution::Spher:
0083     case pat::CandKinResolution::ESpher:
0084     case pat::CandKinResolution::MCSpher: {
0085       // pt = p * sin(theta)
0086       double a = sin(p4.Theta());
0087       double b = p4.P() * cos(p4.Theta());
0088       return sqrt(a * a * covariance(0, 0) + b * b * covariance(1, 1) + 2 * a * b * covariance(0, 1));
0089     }
0090     // ==== SPHERICAL (1/P)  ====
0091     case pat::CandKinResolution::MCPInvSpher: {
0092       // pt = (1/pi) * sin(theta)
0093       double p = p4.P();
0094       double a = -(p * p) * sin(p4.Theta());
0095       double b = p * cos(p4.Theta());
0096       return sqrt(a * a * covariance(0, 0) + b * b * covariance(1, 1) + 2 * a * b * covariance(0, 1));
0097     }
0098     // ==== HEP CYLINDRICAL (with Pt = Et) ====
0099     case pat::CandKinResolution::EtThetaPhi:
0100     case pat::CandKinResolution::EtEtaPhi:
0101       return sqrt(covariance(0, 0));
0102     case pat::CandKinResolution::Invalid:
0103       throw cms::Exception("Invalid parametrization") << parametrization;
0104     default:
0105       throw cms::Exception("Not Implemented")
0106           << "getResolPt not yet implemented for parametrization " << parametrization;
0107   }
0108 }
0109 
0110 double pat::helper::ResolutionHelper::getResolPInv(pat::CandKinResolution::Parametrization parametrization,
0111                                                    const AlgebraicSymMatrix44 &covariance,
0112                                                    const pat::CandKinResolution::LorentzVector &p4) {
0113   switch (parametrization) {
0114     // ==== SPHERICAL (P)  ====
0115     case pat::CandKinResolution::Spher:
0116     case pat::CandKinResolution::ESpher:
0117     case pat::CandKinResolution::MCSpher:
0118       return 1.0 / p4.P2() * sqrt(covariance(0, 0));
0119     // ==== SPHERICAL (1/P)  ====
0120     case pat::CandKinResolution::MCPInvSpher:
0121       return sqrt(covariance(0, 0));
0122     // ==== OTHER ====
0123     case pat::CandKinResolution::Cart:
0124     case pat::CandKinResolution::ECart:
0125     case pat::CandKinResolution::MCCart:
0126     case pat::CandKinResolution::EtThetaPhi:
0127     case pat::CandKinResolution::EtEtaPhi:
0128       return 1.0 / p4.P2() * getResolP(parametrization, covariance, p4);
0129     case pat::CandKinResolution::Invalid:
0130       throw cms::Exception("Invalid parametrization") << parametrization;
0131     default:
0132       throw cms::Exception("Not Implemented")
0133           << "getResolPInv not yet implemented for parametrization " << parametrization;
0134   }
0135 }
0136 
0137 double pat::helper::ResolutionHelper::getResolPx(pat::CandKinResolution::Parametrization parametrization,
0138                                                  const AlgebraicSymMatrix44 &covariance,
0139                                                  const pat::CandKinResolution::LorentzVector &p4) {
0140   switch (parametrization) {
0141     // ==== CARTESIAN ====
0142     case pat::CandKinResolution::Cart:
0143     case pat::CandKinResolution::ECart:
0144     case pat::CandKinResolution::MCCart:
0145       return sqrt(covariance(0, 0));
0146     // ==== SPHERICAL (P)  ====
0147     case pat::CandKinResolution::Spher:
0148     case pat::CandKinResolution::ESpher:
0149     case pat::CandKinResolution::MCSpher:
0150     case pat::CandKinResolution::MCPInvSpher: {
0151       // Px = P * sin(theta) * cos(phi)
0152       double p = p4.P();
0153       AlgebraicVector3 derivs;
0154       derivs[0] = sin(p4.Theta()) * cos(p4.Phi());  // now let's hope gcc does common subexpr optimiz.
0155       if (parametrization == pat::CandKinResolution::MCPInvSpher) {
0156         derivs[0] *= -(p * p);
0157       }
0158       derivs[1] = p * cos(p4.Theta()) * cos(p4.Phi());
0159       derivs[2] = p * sin(p4.Theta()) * -sin(p4.Phi());
0160       return sqrt(ROOT::Math::Similarity(derivs, covariance.Sub<AlgebraicSymMatrix33>(0, 0)));
0161     }
0162     // ==== HEP CYLINDRICAL (with Pt = Et) ====
0163     case pat::CandKinResolution::EtThetaPhi:
0164     case pat::CandKinResolution::EtEtaPhi: {
0165       // Px = Pt * cos(phi)
0166       double a = cos(p4.Phi());
0167       double b = -p4.Pt() * sin(p4.Phi());
0168       return sqrt(a * a * covariance(0, 0) + 2 * a * b * covariance(2, 0) + b * b * covariance(2, 2));
0169     }
0170     // ==== OTHERS ====
0171     case pat::CandKinResolution::Invalid:
0172       throw cms::Exception("Invalid parametrization") << parametrization;
0173     default:
0174       throw cms::Exception("Not Implemented")
0175           << "getResolPx not yet implemented for parametrization " << parametrization;
0176   }
0177 }
0178 double pat::helper::ResolutionHelper::getResolPy(pat::CandKinResolution::Parametrization parametrization,
0179                                                  const AlgebraicSymMatrix44 &covariance,
0180                                                  const pat::CandKinResolution::LorentzVector &p4) {
0181   switch (parametrization) {
0182     // ==== CARTESIAN ====
0183     case pat::CandKinResolution::Cart:
0184     case pat::CandKinResolution::ECart:
0185     case pat::CandKinResolution::MCCart:
0186       return sqrt(covariance(1, 1));
0187     // ==== SPHERICAL (P)  ====
0188     case pat::CandKinResolution::Spher:
0189     case pat::CandKinResolution::ESpher:
0190     case pat::CandKinResolution::MCSpher:
0191     case pat::CandKinResolution::MCPInvSpher: {
0192       // Py = P * sin(theta) * sin(phi)
0193       double p = p4.P();
0194       AlgebraicVector3 derivs;
0195       derivs[0] = sin(p4.Theta()) * sin(p4.Phi());  // now let's hope gcc does common subexpr optimiz.
0196       if (parametrization == pat::CandKinResolution::MCPInvSpher) {
0197         derivs[0] *= -(p * p);
0198       }
0199       derivs[1] = p * cos(p4.Theta()) * sin(p4.Phi());
0200       derivs[2] = p * sin(p4.Theta()) * cos(p4.Phi());
0201       return sqrt(ROOT::Math::Similarity(derivs, covariance.Sub<AlgebraicSymMatrix33>(0, 0)));
0202     }
0203     // ==== HEP CYLINDRICAL (with Pt = Et) ====
0204     case pat::CandKinResolution::EtThetaPhi:
0205     case pat::CandKinResolution::EtEtaPhi: {
0206       // Py = Pt * sin(phi)
0207       double a = sin(p4.Phi());
0208       double b = p4.Pt() * cos(p4.Phi());
0209       return sqrt(a * a * covariance(0, 0) + 2 * a * b * covariance(2, 0) + b * b * covariance(2, 2));
0210     }
0211     // ==== OTHERS ====
0212     case pat::CandKinResolution::Invalid:
0213       throw cms::Exception("Invalid parametrization") << parametrization;
0214     default:
0215       throw cms::Exception("Not Implemented")
0216           << "getResolPy not yet implemented for parametrization " << parametrization;
0217   }
0218 }
0219 double pat::helper::ResolutionHelper::getResolPz(pat::CandKinResolution::Parametrization parametrization,
0220                                                  const AlgebraicSymMatrix44 &covariance,
0221                                                  const pat::CandKinResolution::LorentzVector &p4) {
0222   switch (parametrization) {
0223     // ==== CARTESIAN ====
0224     case pat::CandKinResolution::Cart:
0225     case pat::CandKinResolution::ECart:
0226     case pat::CandKinResolution::MCCart:
0227       return sqrt(covariance(2, 2));
0228     // ==== SPHERICAL (P)  ====
0229     case pat::CandKinResolution::Spher:
0230     case pat::CandKinResolution::ESpher:
0231     case pat::CandKinResolution::MCSpher: {
0232       // Pz = P * cos(theta)
0233       double a = cos(p4.Theta());
0234       double b = -p4.P() * sin(p4.Theta());
0235       return sqrt(a * a * covariance(0, 0) + 2 * a * b * covariance(1, 0) + b * b * covariance(1, 1));
0236     }
0237     case pat::CandKinResolution::MCPInvSpher: {
0238       // Pz = P * cos(theta)
0239       double p = p4.P();
0240       double a = -p * p * cos(p4.Theta());
0241       double b = -p * sin(p4.Theta());
0242       return sqrt(a * a * covariance(0, 0) + 2 * a * b * covariance(1, 0) + b * b * covariance(1, 1));
0243     }
0244     // ==== HEP CYLINDRICAL (with Pt = Et) ====
0245     case pat::CandKinResolution::EtThetaPhi: {
0246       // Pz = Pt * ctg(theta)    d ctg(x) = -1/sin^2(x)
0247       double s = sin(p4.Theta()), c = cos(p4.Theta());
0248       double a = c / s;
0249       double b = -p4.Pt() / (s * s);
0250       return sqrt(a * a * covariance(0, 0) + 2 * a * b * covariance(1, 0) + b * b * covariance(1, 1));
0251     }
0252     case pat::CandKinResolution::EtEtaPhi: {
0253       // Pz = Pt * sinh(eta)
0254       double a = sinh(p4.Eta());
0255       double b = p4.Et() * cosh(p4.Eta());
0256       return sqrt(a * a * covariance(0, 0) + 2 * a * b * covariance(1, 0) + b * b * covariance(1, 1));
0257     }
0258     // ==== OTHERS ====
0259     case pat::CandKinResolution::Invalid:
0260       throw cms::Exception("Invalid parametrization") << parametrization;
0261     default:
0262       throw cms::Exception("Not Implemented")
0263           << "getResolPz not yet implemented for parametrization " << parametrization;
0264   }
0265 }
0266 
0267 double pat::helper::ResolutionHelper::getResolE(pat::CandKinResolution::Parametrization parametrization,
0268                                                 const AlgebraicSymMatrix44 &covariance,
0269                                                 const pat::CandKinResolution::LorentzVector &p4) {
0270   switch (parametrization) {
0271     // ======= ENERGY BASED ==========
0272     case pat::CandKinResolution::ECart:
0273     case pat::CandKinResolution::ESpher:
0274       return sqrt(covariance(3, 3));
0275       // ======= ET BASED ==========
0276     case pat::CandKinResolution::EtThetaPhi: {
0277       // E = Et/Sin(theta)
0278       double a = 1.0 / sin(p4.Theta());               // dE/dEt
0279       double b = -a * a * p4.Et() * cos(p4.Theta());  // dE/dTh
0280       return sqrt(a * a * covariance(0, 0) + b * b * covariance(1, 1) + 2 * a * b * covariance(0, 1));
0281     }
0282     case pat::CandKinResolution::EtEtaPhi: {
0283       // E = Et/Sin(Theta(eta))
0284       double th = p4.Theta();
0285       double a = 1.0 / sin(th);          // dE/dEt
0286       double b = a * p4.Et() * cos(th);  // dE/dEta: dTh/dEta = - 1.0/sin(theta) = - dE/dEt
0287       return sqrt(a * a * covariance(0, 0) + b * b * covariance(1, 1) + 2 * a * b * covariance(0, 1));
0288     }
0289       // ======= MASS BASED ==========
0290     case pat::CandKinResolution::Cart: {
0291       AlgebraicVector4 xoE(p4.X(), p4.Y(), p4.Z(), p4.M());
0292       xoE *= 1 / p4.E();
0293       return sqrt(ROOT::Math::Similarity(xoE, covariance));
0294     }
0295     case pat::CandKinResolution::MCCart: {
0296       AlgebraicVector4 xoE(p4.X(), p4.Y(), p4.Z(), 0);
0297       xoE *= 1 / p4.E();
0298       return sqrt(ROOT::Math::Similarity(xoE, covariance));
0299     }
0300     case pat::CandKinResolution::Spher: {
0301       // E = sqrt(P^2 + m^2)
0302       double einv = 1.0 / p4.E();
0303       double a = p4.P() * einv;  // dE/dP
0304       double b = p4.M() * einv;  // dE/dm
0305       return sqrt(a * a * covariance(0, 0) + b * b * covariance(3, 3) + 2 * a * b * covariance(0, 3));
0306     }
0307     case pat::CandKinResolution::MCSpher: {
0308       // E = sqrt(P^2 + m^2); |dE/dP| = |P/E| = P/E
0309       return p4.P() / p4.E() * sqrt(covariance(0, 0));
0310     }
0311     case pat::CandKinResolution::MCPInvSpher:  //
0312     {
0313       // E = sqrt(P^2 + m^2); |dE/d(1/P)| = P^2 |dE/dP| = P^3/E
0314       double p = p4.P();
0315       return p * p * p / p4.E() * sqrt(covariance(0, 0));
0316     }
0317       // ======= OTHER ==========
0318     case pat::CandKinResolution::Invalid:
0319       throw cms::Exception("Invalid parametrization") << parametrization;
0320     default:
0321       throw cms::Exception("Not Implemented")
0322           << "getResolE not yet implemented for parametrization " << parametrization;
0323   }
0324 }
0325 
0326 double pat::helper::ResolutionHelper::getResolEt(pat::CandKinResolution::Parametrization parametrization,
0327                                                  const AlgebraicSymMatrix44 &covariance,
0328                                                  const pat::CandKinResolution::LorentzVector &p4) {
0329   switch (parametrization) {
0330     // ======= ENERGY BASED ==========
0331     case pat::CandKinResolution::ECart: {
0332       // Et^2 = E^2 * (Pt^2/P^2)
0333       double pt2 = p4.Perp2();
0334       double pz2 = ROOT::Math::Square(p4.Pz()), p2 = pt2 + pz2;
0335       double e2OverP4 = ROOT::Math::Square(p4.E() / p2);
0336       AlgebraicVector4 derivs(p4.Px(), p4.Py(), p4.Pz(), p4.E());
0337       derivs *= (1.0 / p4.Et());
0338       derivs[0] *= pz2 * e2OverP4;
0339       derivs[1] *= pz2 * e2OverP4;
0340       derivs[2] *= -pt2 * e2OverP4;
0341       derivs[3] *= pt2 / p2;
0342       return sqrt(ROOT::Math::Similarity(derivs, covariance));
0343     }
0344     case pat::CandKinResolution::ESpher: {
0345       // Et = E * Sin(Theta)
0346       double st = sin(p4.Theta()), ct = cos(p4.Theta());
0347       return sqrt(st * st * covariance(3, 3) + ROOT::Math::Square(ct * p4.E()) * covariance(1, 1) +
0348                   2 * st * ct * p4.E() * covariance(1, 3));
0349     }
0350       // ======= ET BASED ==========
0351     case pat::CandKinResolution::EtThetaPhi:
0352     case pat::CandKinResolution::EtEtaPhi:
0353       return sqrt(covariance(0, 0));
0354       // ======= MASS BASED ==========
0355     case pat::CandKinResolution::Cart:
0356     case pat::CandKinResolution::MCCart: {
0357       // Et^2 = E^2 Sin^2(th) = (p^2 + m^2) * (pt^2) / p^2
0358       double pt2 = p4.Perp2();
0359       double p2 = pt2 + ROOT::Math::Square(p4.Pz());
0360       double e2 = p2 + p4.M2();
0361       double s2 = pt2 / p2, pi2 = 1.0 / p2;
0362       double et = sqrt(e2 * s2);
0363       AlgebraicVector4 derivs(p4.Px(), p4.Py(), p4.Pz(), p4.M());
0364       derivs *= 1.0 / et;
0365       derivs[0] *= (s2 + e2 * pi2 * (1.0 - pt2 * pi2));  /// dEt/dPx * Et
0366       derivs[1] *= (s2 + e2 * pi2 * (1.0 - pt2 * pi2));  /// dEt/dPx * Et
0367       derivs[2] *= (s2 - e2 * pt2 * pi2 * pi2);          /// dEt/dPx * Et
0368       if (parametrization == pat::CandKinResolution::Cart) {
0369         derivs[3] *= s2;
0370         return sqrt(ROOT::Math::Similarity(derivs, covariance));
0371       } else {
0372         derivs[3] = 0;
0373         return sqrt(ROOT::Math::Similarity(derivs, covariance));  // test if Sub<33> is faster
0374       }
0375     }
0376     case pat::CandKinResolution::Spher: {
0377       // Et = E sin(theta); dE/dM = M/E, dE/dP = P/E
0378       double s = sin(p4.Theta()), c = cos(p4.Theta());
0379       double e = p4.E();
0380       AlgebraicVector4 derivs(p4.P() / e * s, e * c, 0, p4.M() / e * s);
0381       return sqrt(ROOT::Math::Similarity(derivs, covariance));
0382     }
0383     case pat::CandKinResolution::MCSpher: {
0384       // Et = E sin(theta); dE/dP = P/E
0385       double s = sin(p4.Theta()), c = cos(p4.Theta());
0386       double e = p4.E();
0387       double a = p4.P() * s / e;
0388       double b = e * c;
0389       return sqrt(a * a * covariance(0, 0) + b * b * covariance(1, 1) + 2 * a * b * covariance(0, 1));
0390     }
0391     case pat::CandKinResolution::MCPInvSpher: {
0392       // Et = E sin(theta); dE/dP = P/E -> dE/d(1/P) = - P^2 dE/dP = - P^3 / E
0393       double s = sin(p4.Theta()), c = cos(p4.Theta());
0394       double p = p4.P(), e = p4.E();
0395       double a = (-p * p * p / e) * s;
0396       double b = e * c;
0397       return sqrt(a * a * covariance(0, 0) + b * b * covariance(1, 1) + 2 * a * b * covariance(0, 1));
0398     }
0399       // ======= OTHER ==========
0400     case pat::CandKinResolution::Invalid:
0401       throw cms::Exception("Invalid parametrization") << parametrization;
0402     default:
0403       throw cms::Exception("Not Implemented")
0404           << "getResolEt not yet implemented for parametrization " << parametrization;
0405   }
0406 }
0407 
0408 double pat::helper::ResolutionHelper::getResolM(pat::CandKinResolution::Parametrization parametrization,
0409                                                 const AlgebraicSymMatrix44 &covariance,
0410                                                 const pat::CandKinResolution::LorentzVector &p4) {
0411   switch (parametrization) {
0412     // ====== MASS CONSTRAINED =====
0413     case pat::CandKinResolution::MCSpher:
0414     case pat::CandKinResolution::MCPInvSpher:
0415     case pat::CandKinResolution::MCCart:
0416     case pat::CandKinResolution::EtThetaPhi:
0417     case pat::CandKinResolution::EtEtaPhi:
0418       return 0;
0419     // ======= MASS BASED ==========
0420     case pat::CandKinResolution::Cart:
0421     case pat::CandKinResolution::Spher:
0422       return sqrt(covariance(3, 3));
0423     // ======= ENERGY BASED ==========
0424     case pat::CandKinResolution::ESpher: {  // M^2 = E^2 - P^2
0425       double dMdE = p4.E() / p4.M(), dMdP = -p4.P() / p4.M();
0426       return sqrt(dMdP * dMdP * covariance(0, 0) + 2 * dMdP * dMdE * covariance(0, 3) + dMdE * dMdE * covariance(3, 3));
0427     }
0428     case pat::CandKinResolution::ECart: {  // M^2 = E^2 - sum_i P_i^2
0429       AlgebraicVector4 derivs(-p4.Px(), -p4.Py(), -p4.Pz(), p4.E());
0430       derivs *= 1.0 / p4.M();
0431       return sqrt(ROOT::Math::Similarity(derivs, covariance));
0432     }
0433       throw cms::Exception("Not Implemented")
0434           << "getResolM not yet implemented for parametrization " << parametrization;
0435     case pat::CandKinResolution::Invalid:
0436       throw cms::Exception("Invalid parametrization") << parametrization;
0437     default:
0438       throw cms::Exception("Not Implemented")
0439           << "getResolM not yet implemented for parametrization " << parametrization;
0440   }
0441 }
0442 
0443 inline double DetaDtheta(double theta) {
0444   // y  = -ln(tg(x/2)) =>
0445   // y' = - 1/tg(x/2) * 1/(cos(x/2))^2 * 1/2 = - 1 / (2 * sin(x/2) * cos(x/2)) = -1/sin(x)
0446   return -1.0 / sin(theta);
0447 }
0448 inline double DthetaDeta(double eta) {
0449   // y = 2 atan(exp(-x))
0450   // y' = 2 * 1/(1+exp^2) * exp(-x) * (-1) = - 2 * exp/(1+exp^2) = - 2 / (exp + 1/exp)
0451   double e = exp(-eta);
0452   return -2.0 / (e + 1.0 / e);
0453 }
0454 
0455 double pat::helper::ResolutionHelper::getResolEta(pat::CandKinResolution::Parametrization parametrization,
0456                                                   const AlgebraicSymMatrix44 &covariance,
0457                                                   const pat::CandKinResolution::LorentzVector &p4) {
0458   switch (parametrization) {
0459     case pat::CandKinResolution::Cart:
0460     case pat::CandKinResolution::ECart:
0461     case pat::CandKinResolution::MCCart:
0462       // dEta = dTheta * dEta/dTheta
0463       return abs(DetaDtheta(p4.Theta())) * getResolTheta(parametrization, covariance, p4);
0464     case pat::CandKinResolution::ESpher:       //  all the ones which have
0465     case pat::CandKinResolution::MCPInvSpher:  //  theta as parameter 1
0466     case pat::CandKinResolution::MCSpher:
0467     case pat::CandKinResolution::EtThetaPhi:
0468     case pat::CandKinResolution::Spher:
0469       return sqrt(covariance(1, 1)) * abs(DetaDtheta(p4.Theta()));
0470     case pat::CandKinResolution::EtEtaPhi:  // as simple as that
0471       return sqrt(covariance(1, 1));
0472     case pat::CandKinResolution::Invalid:
0473       throw cms::Exception("Invalid parametrization") << parametrization;
0474     default:
0475       throw cms::Exception("Not Implemented")
0476           << "getResolEta not yet implemented for parametrization " << parametrization;
0477   }
0478 }
0479 double pat::helper::ResolutionHelper::getResolTheta(pat::CandKinResolution::Parametrization parametrization,
0480                                                     const AlgebraicSymMatrix44 &covariance,
0481                                                     const pat::CandKinResolution::LorentzVector &p4) {
0482   switch (parametrization) {
0483     case pat::CandKinResolution::Cart:
0484     case pat::CandKinResolution::ECart:
0485     case pat::CandKinResolution::MCCart: {
0486       // theta = acos( pz / p )        ; d acos(x) = - 1 / sqrt( 1 - x*x) dx = - p/pt dx
0487       double pt2 = p4.Perp2();
0488       double p = p4.P(), pi = 1.0 / p, pi3 = pi * pi * pi;
0489       double dacos = -p / sqrt(pt2);
0490       AlgebraicVector3 derivs;
0491       derivs[0] = -p4.Px() * p4.Pz() * dacos * pi3;
0492       derivs[1] = -p4.Py() * p4.Pz() * dacos * pi3;
0493       derivs[2] = pt2 * dacos * pi3;
0494       return sqrt(ROOT::Math::Similarity(derivs, covariance.Sub<AlgebraicSymMatrix33>(0, 0)));
0495     }
0496     case pat::CandKinResolution::ESpher:       //  all the ones which have
0497     case pat::CandKinResolution::MCPInvSpher:  //  theta as parameter 1
0498     case pat::CandKinResolution::MCSpher:
0499     case pat::CandKinResolution::EtThetaPhi:
0500     case pat::CandKinResolution::Spher:
0501       return sqrt(covariance(1, 1));
0502     case pat::CandKinResolution::EtEtaPhi:
0503       return sqrt(covariance(1, 1)) * abs(DthetaDeta(p4.Eta()));
0504     case pat::CandKinResolution::Invalid:
0505       throw cms::Exception("Invalid parametrization") << parametrization;
0506     default:
0507       throw cms::Exception("Not Implemented")
0508           << "getResolTheta not yet implemented for parametrization " << parametrization;
0509   }
0510 }
0511 double pat::helper::ResolutionHelper::getResolPhi(pat::CandKinResolution::Parametrization parametrization,
0512                                                   const AlgebraicSymMatrix44 &covariance,
0513                                                   const pat::CandKinResolution::LorentzVector &p4) {
0514   double pt2 = p4.Perp2();
0515   switch (parametrization) {
0516     case pat::CandKinResolution::Cart:
0517     case pat::CandKinResolution::ECart:
0518     case pat::CandKinResolution::MCCart:
0519       return sqrt(ROOT::Math::Square(p4.Px()) * covariance(1, 1) + ROOT::Math::Square(p4.Py()) * covariance(0, 0) +
0520                   -2 * p4.Px() * p4.Py() * covariance(0, 1)) /
0521              pt2;
0522     case pat::CandKinResolution::ESpher:       //  all the ones which have
0523     case pat::CandKinResolution::MCPInvSpher:  //  phi as parameter 2
0524     case pat::CandKinResolution::MCSpher:
0525     case pat::CandKinResolution::EtThetaPhi:
0526     case pat::CandKinResolution::Spher:
0527     case pat::CandKinResolution::EtEtaPhi:
0528       return sqrt(covariance(2, 2));
0529     case pat::CandKinResolution::Invalid:
0530       throw cms::Exception("Invalid parametrization") << parametrization;
0531     default:
0532       throw cms::Exception("Not Implemented")
0533           << "getResolPhi not yet implemented for parametrization " << parametrization;
0534   }
0535 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!