Back to home page

Project CMSSW displayed by LXR

 
 

    


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

0001 #include "DataFormats/PatCandidates/interface/ParametrizationHelper.h"
0002 #include <cmath>
0003 #include <Math/CylindricalEta3D.h>
0004 #include <Math/Polar3D.h>
0005 #include <Math/Cartesian3D.h>
0006 #include <Math/DisplacementVector3D.h>
0007 #include <Math/Functions.h>
0008 
0009 using namespace std;
0010 
0011 //pat::helper::ParametrizationHelper::POLAR_ZERO();
0012 //pat::helper::ParametrizationHelper::CARTESIAN_ZERO();
0013 
0014 pat::CandKinResolution::Parametrization pat::helper::ParametrizationHelper::fromString(const std::string &name) {
0015   using namespace std;
0016 
0017   typedef pat::CandKinResolution::Parametrization Parametrization;
0018 
0019   static map<string, Parametrization> const parMaps = {{"Cart", pat::CandKinResolution::Cart},
0020                                                        {"ECart", pat::CandKinResolution::ECart},
0021                                                        {"MCCart", pat::CandKinResolution::MCCart},
0022                                                        {"Spher", pat::CandKinResolution::Spher},
0023                                                        {"ESpher", pat::CandKinResolution::ESpher},
0024                                                        {"MCSpher", pat::CandKinResolution::MCSpher},
0025                                                        {"MCPInvSpher", pat::CandKinResolution::MCPInvSpher},
0026                                                        {"EtEtaPhi", pat::CandKinResolution::EtEtaPhi},
0027                                                        {"EtThetaPhi", pat::CandKinResolution::EtThetaPhi},
0028                                                        {"MomDev", pat::CandKinResolution::MomDev},
0029                                                        {"EMomDev", pat::CandKinResolution::EMomDev},
0030                                                        {"MCMomDev", pat::CandKinResolution::MCMomDev},
0031                                                        {"EScaledMomDev", pat::CandKinResolution::EScaledMomDev}};
0032   map<string, Parametrization>::const_iterator itP = parMaps.find(name);
0033   if (itP == parMaps.end()) {
0034     throw cms::Exception("StringResolutionProvider") << "Bad parametrization '" << name.c_str() << "'";
0035   }
0036   return itP->second;
0037 }
0038 
0039 const char *pat::helper::ParametrizationHelper::name(pat::CandKinResolution::Parametrization param) {
0040   switch (param) {
0041     case pat::CandKinResolution::Cart:
0042       return "Cart";
0043     case pat::CandKinResolution::ECart:
0044       return "ECart";
0045     case pat::CandKinResolution::MCCart:
0046       return "MCCart";
0047     case pat::CandKinResolution::Spher:
0048       return "Spher";
0049     case pat::CandKinResolution::ESpher:
0050       return "ESpher";
0051     case pat::CandKinResolution::MCSpher:
0052       return "MCSpher";
0053     case pat::CandKinResolution::MCPInvSpher:
0054       return "MCPInvSpher";
0055     case pat::CandKinResolution::EtEtaPhi:
0056       return "EtEtaPhi";
0057     case pat::CandKinResolution::EtThetaPhi:
0058       return "EtThetaPhi";
0059     case pat::CandKinResolution::MomDev:
0060       return "MomDev";
0061     case pat::CandKinResolution::EMomDev:
0062       return "EMomDev";
0063     case pat::CandKinResolution::MCMomDev:
0064       return "MCMomDev";
0065     case pat::CandKinResolution::EScaledMomDev:
0066       return "EScaledMomDev";
0067     case pat::CandKinResolution::Invalid:
0068       return "Invalid";
0069     default:
0070       return "UNKNOWN";
0071   }
0072 }
0073 
0074 math::PtEtaPhiMLorentzVector pat::helper::ParametrizationHelper::polarP4fromParameters(
0075     pat::CandKinResolution::Parametrization parametrization,
0076     const AlgebraicVector4 &parameters,
0077     const math::PtEtaPhiMLorentzVector &initialP4) {
0078   math::PtEtaPhiMLorentzVector ret;
0079   ROOT::Math::CylindricalEta3D<double> converter;
0080   double m2;
0081   switch (parametrization) {
0082     // ======= CARTESIAN ==========
0083     case pat::CandKinResolution::Cart:
0084       ret = math::XYZTLorentzVector(parameters[0],
0085                                     parameters[1],
0086                                     parameters[2],
0087                                     sqrt(parameters[0] * parameters[0] + parameters[1] * parameters[1] +
0088                                          parameters[2] * parameters[2] + parameters[3] * parameters[3]));
0089       ret.SetM(parameters[3]);  // to be sure about roundoffs
0090       break;
0091     case pat::CandKinResolution::MCCart:
0092       ret = math::XYZTLorentzVector(parameters[0],
0093                                     parameters[1],
0094                                     parameters[2],
0095                                     sqrt(parameters[0] * parameters[0] + parameters[1] * parameters[1] +
0096                                          parameters[2] * parameters[2] + initialP4.mass() * initialP4.mass()));
0097       ret.SetM(initialP4.mass());  // to be sure about roundoffs
0098       break;
0099     case pat::CandKinResolution::ECart:
0100       ret = math::XYZTLorentzVector(parameters[0], parameters[1], parameters[2], parameters[3]);
0101       break;
0102     // ======= SPHERICAL ==========
0103     case pat::CandKinResolution::Spher:
0104       converter = ROOT::Math::Polar3D<double>(parameters[0], parameters[1], 0);
0105       ret.SetCoordinates(converter.Rho(), converter.Eta(), parameters[2], parameters[3]);
0106       break;
0107     case pat::CandKinResolution::MCSpher:  // same as above
0108       converter = ROOT::Math::Polar3D<double>(parameters[0], parameters[1], 0);
0109       ret.SetCoordinates(converter.Rho(), converter.Eta(), parameters[2], initialP4.mass());
0110       break;
0111     case pat::CandKinResolution::ESpher:  //
0112       converter = ROOT::Math::Polar3D<double>(parameters[0], parameters[1], 0);
0113       m2 = -parameters[0] * parameters[0] + parameters[3] * parameters[3];
0114       ret.SetCoordinates(converter.Rho(), converter.Eta(), parameters[2], (m2 > 0 ? sqrt(m2) : 0.0));
0115       break;
0116     case pat::CandKinResolution::MCPInvSpher:  //
0117       converter = ROOT::Math::Polar3D<double>(1.0 / parameters[0], parameters[1], 0);
0118       ret.SetCoordinates(converter.Rho(), converter.Eta(), parameters[2], initialP4.mass());
0119       break;
0120     // ======= HEP CYLINDRICAL ==========
0121     case pat::CandKinResolution::EtThetaPhi:
0122       converter = ROOT::Math::Polar3D<double>(1.0, parameters[1], 0);
0123       ret.SetCoordinates(parameters[0], converter.Eta(), parameters[2], 0);
0124       break;
0125     case pat::CandKinResolution::EtEtaPhi:  // as simple as that
0126       ret.SetCoordinates(parameters[0], parameters[1], parameters[2], 0);
0127       break;
0128     // ======= MomentumDeviates ==========
0129     case pat::CandKinResolution::MomDev:
0130     case pat::CandKinResolution::EMomDev:
0131     case pat::CandKinResolution::MCMomDev:
0132     case pat::CandKinResolution::EScaledMomDev: {
0133       ROOT::Math::DisplacementVector3D<ROOT::Math::Cartesian3D<double> > p = initialP4.Vect(), uz(0, 0, 1), uph, uth;
0134       uph = uz.Cross(p).Unit();
0135       uth = p.Cross(uph).Unit();
0136       p *= parameters[0];
0137       p += uth * parameters[1] + uph * parameters[2];
0138       if (parametrization == pat::CandKinResolution::MomDev) {
0139         ret.SetCoordinates(p.Rho(), p.Eta(), p.Phi(), initialP4.mass() * parameters[3]);
0140       } else if (parametrization == pat::CandKinResolution::EMomDev) {
0141         double m2 = ROOT::Math::Square(parameters[3] * initialP4.energy()) - p.Mag2();
0142         ret.SetCoordinates(p.Rho(), p.Eta(), p.Phi(), (m2 > 0 ? sqrt(m2) : 0.0));
0143       } else if (parametrization == pat::CandKinResolution::EScaledMomDev) {
0144         double m2 = ROOT::Math::Square(p.R() * initialP4.E() / initialP4.P()) - p.Mag2();
0145         ret.SetCoordinates(p.Rho(), p.Eta(), p.Phi(), (m2 > 0 ? sqrt(m2) : 0.0));
0146       } else if (parametrization == pat::CandKinResolution::MCMomDev) {
0147         ret.SetCoordinates(p.Rho(), p.Eta(), p.Phi(), initialP4.mass());
0148       }
0149       break;
0150     }
0151     // ======= OTHER ==========
0152     case pat::CandKinResolution::Invalid:
0153       throw cms::Exception("Invalid parametrization") << parametrization;
0154     default:
0155       throw cms::Exception("Not Implemented")
0156           << "getResolEta not yet implemented for parametrization " << parametrization;
0157   }
0158   return ret;
0159 }
0160 
0161 math::XYZTLorentzVector pat::helper::ParametrizationHelper::p4fromParameters(
0162     pat::CandKinResolution::Parametrization parametrization,
0163     const AlgebraicVector4 &parameters,
0164     const math::XYZTLorentzVector &initialP4) {
0165   math::XYZTLorentzVector ret;
0166   switch (parametrization) {
0167     // ======= CARTESIAN ==========
0168     case pat::CandKinResolution::Cart:
0169       ret.SetCoordinates(parameters[0],
0170                          parameters[1],
0171                          parameters[2],
0172                          sqrt(parameters[0] * parameters[0] + parameters[1] * parameters[1] +
0173                               parameters[2] * parameters[2] + parameters[3] * parameters[3]));
0174       break;
0175     case pat::CandKinResolution::MCCart:  // same as above
0176       ret.SetCoordinates(parameters[0],
0177                          parameters[1],
0178                          parameters[2],
0179                          sqrt(parameters[0] * parameters[0] + parameters[1] * parameters[1] +
0180                               parameters[2] * parameters[2] + initialP4.mass() * initialP4.mass()));
0181       break;
0182     case pat::CandKinResolution::ECart:
0183       ret.SetCoordinates(parameters[0], parameters[1], parameters[2], parameters[3]);
0184       break;
0185     // ======= MomentumDeviates ==========
0186     case pat::CandKinResolution::MomDev:
0187     case pat::CandKinResolution::EMomDev:
0188     case pat::CandKinResolution::MCMomDev:
0189     case pat::CandKinResolution::EScaledMomDev: {
0190       ROOT::Math::DisplacementVector3D<ROOT::Math::Cartesian3D<double> > p = initialP4.Vect(), uz(0, 0, 1), uph, uth;
0191       uph = uz.Cross(p).Unit();
0192       uth = p.Cross(uph).Unit();
0193       p *= parameters[0];
0194       p += uth * parameters[1] + uph * parameters[2];
0195       if (parametrization == pat::CandKinResolution::MomDev) {
0196         ret.SetCoordinates(p.X(), p.Y(), p.Z(), sqrt(p.Mag2() + ROOT::Math::Square(initialP4.mass() * parameters[3])));
0197       } else if (parametrization == pat::CandKinResolution::EMomDev) {
0198         ret.SetCoordinates(p.X(), p.Y(), p.Z(), parameters[3] * initialP4.energy());
0199       } else if (parametrization == pat::CandKinResolution::EMomDev) {
0200         ret.SetCoordinates(p.X(), p.Y(), p.Z(), p.R() * initialP4.E() / initialP4.P());
0201       } else {
0202         ret.SetCoordinates(p.X(), p.Y(), p.Z(), sqrt(p.Mag2() + initialP4.mass() * initialP4.mass()));
0203       }
0204       break;
0205     }
0206     // ======= ALL OTHERS ==========
0207     default:
0208       ret = polarP4fromParameters(parametrization, parameters, math::PtEtaPhiMLorentzVector(initialP4));
0209   }
0210   return ret;
0211 }
0212 
0213 template <typename T>
0214 void pat::helper::ParametrizationHelper::setParametersFromAnyVector(
0215     pat::CandKinResolution::Parametrization parametrization, AlgebraicVector4 &ret, const T &p4) {
0216   switch (parametrization) {
0217     // ======= CARTESIAN ==========
0218     case pat::CandKinResolution::Cart:
0219       ret[0] = p4.px();
0220       ret[1] = p4.py();
0221       ret[2] = p4.pz();
0222       ret[3] = p4.mass();
0223       break;
0224     case pat::CandKinResolution::MCCart:
0225       ret[0] = p4.px();
0226       ret[1] = p4.py();
0227       ret[2] = p4.pz();
0228       ret[3] = p4.mass();
0229       break;
0230     case pat::CandKinResolution::ECart:
0231       ret[0] = p4.px();
0232       ret[1] = p4.py();
0233       ret[2] = p4.pz();
0234       ret[3] = p4.energy();
0235       break;
0236       // ======= SPHERICAL ==========
0237     case pat::CandKinResolution::Spher:
0238       ret[0] = p4.P();
0239       ret[1] = p4.theta();
0240       ret[2] = p4.phi();
0241       ret[3] = p4.mass();
0242       break;
0243     case pat::CandKinResolution::MCSpher:
0244       ret[0] = p4.P();
0245       ret[1] = p4.theta();
0246       ret[2] = p4.phi();
0247       ret[3] = p4.mass();
0248       break;
0249     case pat::CandKinResolution::ESpher:
0250       ret[0] = p4.P();
0251       ret[1] = p4.theta();
0252       ret[2] = p4.phi();
0253       ret[3] = p4.energy();
0254       break;
0255     case pat::CandKinResolution::MCPInvSpher:
0256       ret[0] = 1.0 / p4.P();
0257       ret[1] = p4.theta();
0258       ret[2] = p4.phi();
0259       ret[3] = p4.mass();
0260       break;
0261     // ======= HEP CYLINDRICAL ==========
0262     case pat::CandKinResolution::EtThetaPhi:
0263       ret[0] = p4.pt();
0264       ret[1] = p4.theta();
0265       ret[2] = p4.phi();
0266       ret[3] = 0;
0267       break;
0268     case pat::CandKinResolution::EtEtaPhi:
0269       ret[0] = p4.pt();
0270       ret[1] = p4.eta();
0271       ret[2] = p4.phi();
0272       ret[3] = 0;
0273       break;
0274     // ======= DEVIATES ==========
0275     case pat::CandKinResolution::MomDev:
0276     case pat::CandKinResolution::EMomDev:
0277       ret[0] = 1.0;
0278       ret[1] = 0.0;
0279       ret[2] = 0.0;
0280       ret[3] = 1.0;
0281       break;
0282     case pat::CandKinResolution::MCMomDev:
0283       ret[0] = 1.0;
0284       ret[1] = 0.0;
0285       ret[2] = 0.0;
0286       ret[3] = p4.mass();
0287       break;
0288     case pat::CandKinResolution::EScaledMomDev:
0289       ret[0] = 1.0;
0290       ret[1] = 0.0;
0291       ret[2] = 0.0;
0292       ret[3] = p4.E() / p4.P();
0293       break;
0294     // ======= OTHER ==========
0295     case pat::CandKinResolution::Invalid:
0296       throw cms::Exception("Invalid parametrization") << parametrization;
0297     default:
0298       throw cms::Exception("Not Implemented")
0299           << "getResolEta not yet implemented for parametrization " << parametrization;
0300   }
0301 }
0302 
0303 void pat::helper::ParametrizationHelper::setParametersFromP4(pat::CandKinResolution::Parametrization parametrization,
0304                                                              AlgebraicVector4 &v,
0305                                                              const math::PtEtaPhiMLorentzVector &p4) {
0306   setParametersFromAnyVector(parametrization, v, p4);
0307 }
0308 
0309 void pat::helper::ParametrizationHelper::setParametersFromP4(pat::CandKinResolution::Parametrization parametrization,
0310                                                              AlgebraicVector4 &v,
0311                                                              const math::XYZTLorentzVector &p4) {
0312   setParametersFromAnyVector(parametrization, v, p4);
0313 }
0314 
0315 AlgebraicVector4 pat::helper::ParametrizationHelper::parametersFromP4(
0316     pat::CandKinResolution::Parametrization parametrization, const math::PtEtaPhiMLorentzVector &p4) {
0317   AlgebraicVector4 ret;
0318   setParametersFromP4(parametrization, ret, p4);
0319   return ret;
0320 }
0321 
0322 AlgebraicVector4 pat::helper::ParametrizationHelper::parametersFromP4(
0323     pat::CandKinResolution::Parametrization parametrization, const math::XYZTLorentzVector &p4) {
0324   AlgebraicVector4 ret;
0325   setParametersFromP4(parametrization, ret, p4);
0326   return ret;
0327 }
0328 
0329 AlgebraicVector4 pat::helper::ParametrizationHelper::diffToParameters(
0330     pat::CandKinResolution::Parametrization parametrization,
0331     const math::PtEtaPhiMLorentzVector &p4ini,
0332     const math::PtEtaPhiMLorentzVector &p4fin) {
0333   AlgebraicVector4 ret;
0334   switch (parametrization) {
0335     case pat::CandKinResolution::Cart:
0336     case pat::CandKinResolution::ECart:
0337     case pat::CandKinResolution::MCCart:
0338       ret = parametersFromP4(parametrization, p4fin) - parametersFromP4(parametrization, p4ini);
0339       break;
0340     case pat::CandKinResolution::Spher:
0341     case pat::CandKinResolution::ESpher:
0342     case pat::CandKinResolution::MCSpher:
0343     case pat::CandKinResolution::MCPInvSpher:
0344     case pat::CandKinResolution::EtThetaPhi:
0345     case pat::CandKinResolution::EtEtaPhi:
0346       ret = parametersFromP4(parametrization, p4fin) - parametersFromP4(parametrization, p4ini);
0347       while (ret[2] > +M_PI)
0348         ret[2] -= (2 * M_PI);
0349       while (ret[2] < -M_PI)
0350         ret[2] += (2 * M_PI);
0351       break;
0352     case pat::CandKinResolution::MCMomDev:
0353     case pat::CandKinResolution::MomDev:
0354     case pat::CandKinResolution::EMomDev:
0355     case pat::CandKinResolution::EScaledMomDev:
0356       return diffToParameters(parametrization, math::XYZTLorentzVector(p4ini), math::XYZTLorentzVector(p4fin));
0357     case pat::CandKinResolution::Invalid:
0358       throw cms::Exception("Invalid parametrization") << parametrization;
0359     default:
0360       throw cms::Exception("Not Implemented")
0361           << "diffToParameters not yet implemented for parametrization " << parametrization;
0362   }
0363   return ret;
0364 }
0365 
0366 AlgebraicVector4 pat::helper::ParametrizationHelper::diffToParameters(
0367     pat::CandKinResolution::Parametrization parametrization,
0368     const math::XYZTLorentzVector &p4ini,
0369     const math::XYZTLorentzVector &p4fin) {
0370   AlgebraicVector4 ret;
0371   switch (parametrization) {
0372     case pat::CandKinResolution::Cart:
0373     case pat::CandKinResolution::ECart:
0374     case pat::CandKinResolution::MCCart:
0375     case pat::CandKinResolution::Spher:
0376       ret = parametersFromP4(parametrization, p4fin) - parametersFromP4(parametrization, p4ini);
0377       break;
0378     case pat::CandKinResolution::ESpher:
0379     case pat::CandKinResolution::MCSpher:
0380     case pat::CandKinResolution::MCPInvSpher:
0381     case pat::CandKinResolution::EtThetaPhi:
0382     case pat::CandKinResolution::EtEtaPhi:
0383       ret = parametersFromP4(parametrization, p4fin) - parametersFromP4(parametrization, p4ini);
0384       while (ret[2] > +M_PI)
0385         ret[2] -= (2 * M_PI);
0386       while (ret[2] < -M_PI)
0387         ret[2] += (2 * M_PI);
0388       break;
0389     case pat::CandKinResolution::MCMomDev:
0390     case pat::CandKinResolution::MomDev:
0391     case pat::CandKinResolution::EMomDev:
0392     case pat::CandKinResolution::EScaledMomDev: {
0393       typedef ROOT::Math::DisplacementVector3D<ROOT::Math::Cartesian3D<double> > V3Cart;
0394       V3Cart p1 = p4ini.Vect(), p2 = p4fin.Vect();
0395       V3Cart ur = p1.Unit();
0396       V3Cart uz(0, 0, 1);
0397       V3Cart uph = uz.Cross(ur).Unit();
0398       V3Cart uth = ur.Cross(uph).Unit();
0399       ret[0] = p2.Dot(ur) / p1.R() - 1.0;
0400       ret[1] = (p2 - p1).Dot(uth);
0401       ret[2] = (p2 - p1).Dot(uph);
0402       if (parametrization == pat::CandKinResolution::MomDev) {
0403         ret[3] = p4fin.mass() / p4ini.mass() - 1.0;
0404       } else if (parametrization == pat::CandKinResolution::EMomDev) {
0405         ret[3] = p4fin.energy() / p4ini.energy() - 1.0;
0406       }
0407     } break;
0408     case pat::CandKinResolution::Invalid:
0409       throw cms::Exception("Invalid parametrization") << parametrization;
0410     default:
0411       throw cms::Exception("Not Implemented")
0412           << "diffToParameters not yet implemented for parametrization " << parametrization;
0413   }
0414   return ret;
0415 }
0416 
0417 bool pat::helper::ParametrizationHelper::isAlwaysMassless(pat::CandKinResolution::Parametrization parametrization) {
0418   switch (parametrization) {
0419     case pat::CandKinResolution::Cart:
0420     case pat::CandKinResolution::ECart:
0421     case pat::CandKinResolution::MCCart:
0422     case pat::CandKinResolution::Spher:
0423     case pat::CandKinResolution::ESpher:
0424     case pat::CandKinResolution::MCSpher:
0425     case pat::CandKinResolution::MCPInvSpher:
0426     case pat::CandKinResolution::MCMomDev:
0427     case pat::CandKinResolution::MomDev:
0428     case pat::CandKinResolution::EMomDev:
0429     case pat::CandKinResolution::EScaledMomDev:
0430       return false;
0431     case pat::CandKinResolution::EtThetaPhi:
0432     case pat::CandKinResolution::EtEtaPhi:
0433       return true;
0434     case pat::CandKinResolution::Invalid:
0435       throw cms::Exception("Invalid parametrization") << parametrization;
0436     default:
0437       throw cms::Exception("Not Implemented")
0438           << "isAlwaysMassless not yet implemented for parametrization " << parametrization;
0439   }
0440 }
0441 
0442 bool pat::helper::ParametrizationHelper::isAlwaysMassive(pat::CandKinResolution::Parametrization parametrization) {
0443   switch (parametrization) {
0444     case pat::CandKinResolution::Cart:
0445     case pat::CandKinResolution::ECart:
0446     case pat::CandKinResolution::MCCart:
0447     case pat::CandKinResolution::Spher:
0448     case pat::CandKinResolution::ESpher:
0449     case pat::CandKinResolution::MCSpher:
0450     case pat::CandKinResolution::MCPInvSpher:
0451     case pat::CandKinResolution::MCMomDev:
0452     case pat::CandKinResolution::EMomDev:
0453     case pat::CandKinResolution::EScaledMomDev:
0454     case pat::CandKinResolution::EtThetaPhi:
0455     case pat::CandKinResolution::EtEtaPhi:
0456       return false;
0457     case pat::CandKinResolution::MomDev:
0458       return true;
0459     case pat::CandKinResolution::Invalid:
0460       throw cms::Exception("Invalid parametrization") << parametrization;
0461     default:
0462       throw cms::Exception("Not Implemented")
0463           << "isAlwaysMassless not yet implemented for parametrization " << parametrization;
0464   }
0465 }
0466 
0467 bool pat::helper::ParametrizationHelper::isMassConstrained(pat::CandKinResolution::Parametrization parametrization) {
0468   switch (parametrization) {
0469     case pat::CandKinResolution::Cart:
0470     case pat::CandKinResolution::ECart:
0471     case pat::CandKinResolution::Spher:
0472     case pat::CandKinResolution::ESpher:
0473     case pat::CandKinResolution::EMomDev:
0474     case pat::CandKinResolution::EScaledMomDev:
0475     case pat::CandKinResolution::EtThetaPhi:
0476     case pat::CandKinResolution::EtEtaPhi:
0477     case pat::CandKinResolution::MomDev:
0478       return false;
0479     case pat::CandKinResolution::MCCart:
0480     case pat::CandKinResolution::MCSpher:
0481     case pat::CandKinResolution::MCPInvSpher:
0482     case pat::CandKinResolution::MCMomDev:
0483       return true;
0484     case pat::CandKinResolution::Invalid:
0485       throw cms::Exception("Invalid parametrization") << parametrization;
0486     default:
0487       throw cms::Exception("Not Implemented")
0488           << "isAlwaysMassless not yet implemented for parametrization " << parametrization;
0489   }
0490 }
0491 
0492 bool pat::helper::ParametrizationHelper::isPhysical(pat::CandKinResolution::Parametrization parametrization,
0493                                                     const AlgebraicVector4 &parameters,
0494                                                     const math::PtEtaPhiMLorentzVector &initialP4) {
0495   switch (parametrization) {
0496     // ======= CARTESIAN ==========
0497     case pat::CandKinResolution::Cart:
0498       return parameters[3] >= 0;  // M >= 0
0499     case pat::CandKinResolution::MCCart:
0500       return true;
0501     case pat::CandKinResolution::ECart:
0502       return (parameters[0] * parameters[0] + parameters[1] * parameters[1] + parameters[2] * parameters[2] <=
0503               parameters[3] * parameters[3]);  // E >= P
0504     case pat::CandKinResolution::Spher:
0505       return (parameters[0] >= 0) &&   // P >= 0
0506              (parameters[3] >= 0) &&   // M >= 0
0507              (parameters[1] >= 0) &&   // theta >= 0
0508              (parameters[1] <= M_PI);  // theta <= pi
0509     case pat::CandKinResolution::MCSpher:
0510       return (parameters[0] >= 0) &&              // P >= 0
0511              (parameters[1] >= 0) &&              // theta >= 0
0512              (parameters[1] <= M_PI);             // theta <= pi
0513     case pat::CandKinResolution::ESpher:          //
0514       return (parameters[0] >= 0) &&              // P >= 0
0515              (parameters[3] >= parameters[0]) &&  // E >= P
0516              (parameters[1] >= 0) &&              // theta >= 0
0517              (parameters[1] <= M_PI);             // theta <= PI
0518     case pat::CandKinResolution::MCPInvSpher:     //
0519       return (parameters[0] > 0) &&               // 1/P > 0
0520              (parameters[1] >= 0) &&              // theta >= 0
0521              (parameters[1] <= M_PI);             // theta <= pi
0522     // ======= HEP CYLINDRICAL ==========
0523     case pat::CandKinResolution::EtThetaPhi:
0524       return (parameters[0] > 0) &&         // Et >= 0
0525              (parameters[1] >= 0) &&        // theta >= 0
0526              (parameters[1] <= M_PI);       // theta <= pi
0527     case pat::CandKinResolution::EtEtaPhi:  // as simple as that
0528       return (parameters[0] > 0);           // Et >= 0
0529     // ======= MomentumDeviates ==========
0530     case pat::CandKinResolution::MomDev:
0531       return (parameters[0] >= 0) &&  // P >= 0
0532              (parameters[3] >= 0);    // m/M0 >= 0
0533     case pat::CandKinResolution::MCMomDev:
0534       return (parameters[0] >= 0);  // P >= 0
0535     case pat::CandKinResolution::EMomDev:
0536     case pat::CandKinResolution::EScaledMomDev: {
0537       if (parameters[0] <= 0)
0538         return false;
0539       ROOT::Math::DisplacementVector3D<ROOT::Math::Cartesian3D<double> > p = initialP4.Vect(), uz(0, 0, 1), uph, uth;
0540       uph = uz.Cross(p).Unit();
0541       uth = p.Cross(uph).Unit();
0542       p *= parameters[0];
0543       p += uth * parameters[1] + uph * parameters[2];
0544       if (parametrization == pat::CandKinResolution::EMomDev) {
0545         if (parameters[3] < 0)
0546           return false;
0547         double m2 = ROOT::Math::Square(parameters[3] * initialP4.energy()) - p.Mag2();
0548         if (m2 < 0)
0549           return false;
0550       } else if (parametrization == pat::CandKinResolution::EScaledMomDev) {
0551         if (parameters[3] < 0)
0552           return false;
0553         double m2 = ROOT::Math::Square(p.R() * initialP4.E() / initialP4.P()) - p.Mag2();
0554         if (m2 < 0)
0555           return false;
0556       }
0557       return true;
0558     }
0559     // ======= OTHER ==========
0560     case pat::CandKinResolution::Invalid:
0561       throw cms::Exception("Invalid parametrization") << parametrization;
0562     default:
0563       throw cms::Exception("Not Implemented")
0564           << "getResolEta not yet implemented for parametrization " << parametrization;
0565   }
0566 }