Reference Manual

Go to the documentation of this file.
00001 #include "DataFormats/PatCandidates/interface/ParametrizationHelper.h"
00002 #include <cmath>
00003 #include <Math/CylindricalEta3D.h>
00004 #include <Math/Polar3D.h>
00005 #include <Math/Cartesian3D.h>
00006 #include <Math/DisplacementVector3D.h>
00007 #include <Math/Functions.h>
00008 
00009 using namespace std; 
00010 
00011 //pat::helper::ParametrizationHelper::POLAR_ZERO();
00012 //pat::helper::ParametrizationHelper::CARTESIAN_ZERO();
00013 
00014 pat::CandKinResolution::Parametrization 
00015 pat::helper::ParametrizationHelper::fromString(const std::string &name) {
00016     using namespace std;
00017 
00018     typedef pat::CandKinResolution::Parametrization Parametrization;
00019 
00020     static map<string,Parametrization> parMaps;
00021     if (parMaps.empty()) {
00022         parMaps["Cart"]          = pat::CandKinResolution::Cart;  
00023         parMaps["ECart"]         = pat::CandKinResolution::ECart;  
00024         parMaps["MCCart"]        = pat::CandKinResolution::MCCart;  
00025         parMaps["Spher"]         = pat::CandKinResolution::Spher;  
00026         parMaps["ESpher"]        = pat::CandKinResolution::ESpher;  
00027         parMaps["MCSpher"]       = pat::CandKinResolution::MCSpher;  
00028         parMaps["MCPInvSpher"]   = pat::CandKinResolution::MCPInvSpher;  
00029         parMaps["EtEtaPhi"]      = pat::CandKinResolution::EtEtaPhi;  
00030         parMaps["EtThetaPhi"]    = pat::CandKinResolution::EtThetaPhi;
00031         parMaps["MomDev"]        = pat::CandKinResolution::MomDev;
00032         parMaps["EMomDev"]       = pat::CandKinResolution::EMomDev;
00033         parMaps["MCMomDev"]      = pat::CandKinResolution::MCMomDev;
00034         parMaps["EScaledMomDev"] = pat::CandKinResolution::EScaledMomDev;
00035     }
00036     map<string,Parametrization>::const_iterator itP = parMaps.find(name);
00037     if (itP == parMaps.end()) {
00038         throw cms::Exception("StringResolutionProvider") << "Bad parametrization '" << name.c_str() << "'";
00039     }
00040     return itP->second;
00041 }
00042 
00043 const char * 
00044 pat::helper::ParametrizationHelper::name(pat::CandKinResolution::Parametrization param) {
00045     switch (param) {
00046         case pat::CandKinResolution::Cart:          return "Cart";  
00047         case pat::CandKinResolution::ECart:         return "ECart";  
00048         case pat::CandKinResolution::MCCart:        return "MCCart";  
00049         case pat::CandKinResolution::Spher:         return "Spher";  
00050         case pat::CandKinResolution::ESpher:        return "ESpher";  
00051         case pat::CandKinResolution::MCSpher:       return "MCSpher";  
00052         case pat::CandKinResolution::MCPInvSpher:   return "MCPInvSpher";  
00053         case pat::CandKinResolution::EtEtaPhi:      return "EtEtaPhi";  
00054         case pat::CandKinResolution::EtThetaPhi:    return "EtThetaPhi";
00055         case pat::CandKinResolution::MomDev:        return "MomDev";
00056         case pat::CandKinResolution::EMomDev:       return "EMomDev";
00057         case pat::CandKinResolution::MCMomDev:      return "MCMomDev";
00058         case pat::CandKinResolution::EScaledMomDev: return "EScaledMomDev";
00059         case pat::CandKinResolution::Invalid:       return "Invalid";
00060         default: return "UNKNOWN";
00061     }
00062 }
00063  
00064 math::PtEtaPhiMLorentzVector 
00065 pat::helper::ParametrizationHelper::polarP4fromParameters(pat::CandKinResolution::Parametrization parametrization, 
00066             const AlgebraicVector4 &parameters,
00067             const math::PtEtaPhiMLorentzVector &initialP4) {
00068     math::PtEtaPhiMLorentzVector  ret;
00069     ROOT::Math::CylindricalEta3D<double> converter;
00070     double m2;
00071     switch (parametrization) {
00072         // ======= CARTESIAN ==========
00073         case pat::CandKinResolution::Cart:
00074             ret = math::XYZTLorentzVector(parameters[0], parameters[1], parameters[2], 
00075                     sqrt(parameters[0]*parameters[0] + 
00076                         parameters[1]*parameters[1] + 
00077                         parameters[2]*parameters[2] + 
00078                         parameters[3]*parameters[3])  );
00079             ret.SetM(parameters[3]); // to be sure about roundoffs
00080             break;
00081         case pat::CandKinResolution::MCCart:
00082             ret = math::XYZTLorentzVector(parameters[0], parameters[1], parameters[2], 
00083                     sqrt(parameters[0]*parameters[0] + 
00084                         parameters[1]*parameters[1] + 
00085                         parameters[2]*parameters[2] + 
00086                         initialP4.mass()*initialP4.mass())  );
00087             ret.SetM(initialP4.mass()); // to be sure about roundoffs
00088             break;
00089         case pat::CandKinResolution::ECart:    
00090             ret = math::XYZTLorentzVector(parameters[0], parameters[1], parameters[2], parameters[3]);
00091             break;
00092         // ======= SPHERICAL ==========
00093         case pat::CandKinResolution::Spher:  
00094             converter = ROOT::Math::Polar3D<double>(parameters[0], parameters[1], 0); 
00095             ret.SetCoordinates(converter.Rho(),converter.Eta(),parameters[2],parameters[3]);
00096             break;
00097         case pat::CandKinResolution::MCSpher:    // same as above
00098             converter = ROOT::Math::Polar3D<double>(parameters[0], parameters[1], 0); 
00099             ret.SetCoordinates(converter.Rho(),converter.Eta(),parameters[2],initialP4.mass());
00100             break;
00101         case pat::CandKinResolution::ESpher:        //  
00102             converter = ROOT::Math::Polar3D<double>(parameters[0], parameters[1], 0); 
00103             m2 = - parameters[0]*parameters[0] + parameters[3]*parameters[3];
00104             ret.SetCoordinates(converter.Rho(),converter.Eta(),parameters[2],(m2 > 0 ? sqrt(m2) : 0.0));
00105             break;
00106         case pat::CandKinResolution::MCPInvSpher:   //  
00107             converter = ROOT::Math::Polar3D<double>(1.0/parameters[0], parameters[1], 0); 
00108             ret.SetCoordinates(converter.Rho(),converter.Eta(),parameters[2],initialP4.mass());
00109             break;
00110         // ======= HEP CYLINDRICAL ==========
00111         case pat::CandKinResolution::EtThetaPhi: 
00112             converter = ROOT::Math::Polar3D<double>(1.0, parameters[1], 0); 
00113             ret.SetCoordinates(parameters[0],converter.Eta(),parameters[2],0);
00114             break;
00115         case pat::CandKinResolution::EtEtaPhi:      // as simple as that
00116             ret.SetCoordinates(parameters[0], parameters[1], parameters[2], 0);
00117             break;
00118         // ======= MomentumDeviates ==========
00119         case pat::CandKinResolution::MomDev:
00120         case pat::CandKinResolution::EMomDev:
00121         case pat::CandKinResolution::MCMomDev:
00122         case pat::CandKinResolution::EScaledMomDev:
00123             {
00124                 ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D<double> > p = initialP4.Vect(), uz(0,0,1), uph, uth;
00125                 uph = uz.Cross(p).Unit();
00126                 uth = p.Cross(uph).Unit();
00127                 p *= parameters[0]; 
00128                 p += uth * parameters[1] + uph * parameters[2];
00129                 if (parametrization == pat::CandKinResolution::MomDev) {
00130                     ret.SetCoordinates(p.Rho(), p.Eta(), p.Phi(), initialP4.mass() * parameters[3]);
00131                 } else if (parametrization == pat::CandKinResolution::EMomDev) {
00132                     double m2 = ROOT::Math::Square(parameters[3] * initialP4.energy()) - p.Mag2();
00133                     ret.SetCoordinates(p.Rho(), p.Eta(), p.Phi(), (m2 > 0 ? sqrt(m2) : 0.0));
00134                 } else if (parametrization == pat::CandKinResolution::EScaledMomDev) {
00135                     double m2 = ROOT::Math::Square(p.R()*initialP4.E()/initialP4.P()) - p.Mag2();
00136                     ret.SetCoordinates(p.Rho(), p.Eta(), p.Phi(), (m2 > 0 ? sqrt(m2) : 0.0));
00137                 } else if (parametrization == pat::CandKinResolution::MCMomDev) {
00138                     ret.SetCoordinates(p.Rho(), p.Eta(), p.Phi(), initialP4.mass());
00139                 }
00140                 break;
00141             }
00142         // ======= OTHER ==========
00143         case pat::CandKinResolution::Invalid:
00144             throw cms::Exception("Invalid parametrization") << parametrization;
00145         default:
00146             throw cms::Exception("Not Implemented") << "getResolEta not yet implemented for parametrization " << parametrization ;
00147     }
00148     return ret;
00149 }
00150 
00151 math::XYZTLorentzVector 
00152 pat::helper::ParametrizationHelper::p4fromParameters(pat::CandKinResolution::Parametrization parametrization, 
00153             const AlgebraicVector4 &parameters,
00154             const math::XYZTLorentzVector &initialP4) {
00155     math::XYZTLorentzVector ret;
00156     switch (parametrization) {
00157         // ======= CARTESIAN ==========
00158         case pat::CandKinResolution::Cart:
00159             ret.SetCoordinates(parameters[0], parameters[1], parameters[2], 
00160                     sqrt(parameters[0]*parameters[0] + 
00161                         parameters[1]*parameters[1] + 
00162                         parameters[2]*parameters[2] + 
00163                         parameters[3]*parameters[3])  );
00164             break;
00165         case pat::CandKinResolution::MCCart: // same as above
00166             ret.SetCoordinates(parameters[0], parameters[1], parameters[2], 
00167                     sqrt(parameters[0]*parameters[0] + 
00168                         parameters[1]*parameters[1] + 
00169                         parameters[2]*parameters[2] + 
00170                         initialP4.mass()*initialP4.mass())  );
00171             break;
00172         case pat::CandKinResolution::ECart:    
00173             ret.SetCoordinates(parameters[0], parameters[1], parameters[2], parameters[3]);
00174             break;
00175         // ======= MomentumDeviates ==========
00176         case pat::CandKinResolution::MomDev:
00177         case pat::CandKinResolution::EMomDev:
00178         case pat::CandKinResolution::MCMomDev:
00179         case pat::CandKinResolution::EScaledMomDev:
00180             {
00181                 ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D<double> > p = initialP4.Vect(), uz(0,0,1), uph, uth;
00182                 uph = uz.Cross(p).Unit();
00183                 uth = p.Cross(uph).Unit();
00184                 p *= parameters[0]; 
00185                 p += uth * parameters[1] + uph * parameters[2];
00186                 if (parametrization == pat::CandKinResolution::MomDev) {
00187                     ret.SetCoordinates(p.X(), p.Y(), p.Z(), sqrt(p.Mag2() + ROOT::Math::Square(initialP4.mass() * parameters[3])) );
00188                 } else if (parametrization == pat::CandKinResolution::EMomDev) {
00189                     ret.SetCoordinates(p.X(), p.Y(), p.Z(), parameters[3] * initialP4.energy());
00190                 } else if (parametrization == pat::CandKinResolution::EMomDev) {
00191                     ret.SetCoordinates(p.X(), p.Y(), p.Z(), p.R() * initialP4.E()/initialP4.P());
00192                 } else {
00193                     ret.SetCoordinates(p.X(), p.Y(), p.Z(), sqrt(p.Mag2() + initialP4.mass()*initialP4.mass()));
00194                 }
00195                 break;
00196             }
00197         // ======= ALL OTHERS ==========
00198         default:
00199             ret = polarP4fromParameters(parametrization, parameters, math::PtEtaPhiMLorentzVector(initialP4));
00200     }
00201     return ret;
00202 }
00203 
00204 template <typename T>
00205 void
00206 pat::helper::ParametrizationHelper::setParametersFromAnyVector(pat::CandKinResolution::Parametrization parametrization, 
00207             AlgebraicVector4 &ret, 
00208             const T &p4) {
00209     switch (parametrization) {
00210         // ======= CARTESIAN ==========
00211         case pat::CandKinResolution::Cart:
00212             ret[0] = p4.px(); ret[1] = p4.py(); ret[2] = p4.pz(); ret[3] = p4.mass();
00213             break;
00214         case pat::CandKinResolution::MCCart: 
00215             ret[0] = p4.px(); ret[1] = p4.py(); ret[2] = p4.pz(); ret[3] = p4.mass();        
00216             break;
00217         case pat::CandKinResolution::ECart:    
00218             ret[0] = p4.px(); ret[1] = p4.py(); ret[2] = p4.pz(); ret[3] = p4.energy();        
00219             break;
00220         // ======= SPHERICAL ==========
00221        case pat::CandKinResolution::Spher:  
00222             ret[0] = p4.P(); ret[1] = p4.theta(); ret[2] = p4.phi(); ret[3] = p4.mass();        
00223             break;
00224        case pat::CandKinResolution::MCSpher: 
00225             ret[0] = p4.P(); ret[1] = p4.theta(); ret[2] = p4.phi(); ret[3] = p4.mass();
00226             break;
00227        case pat::CandKinResolution::ESpher:
00228             ret[0] = p4.P(); ret[1] = p4.theta(); ret[2] = p4.phi(); ret[3] = p4.energy();        
00229             break;
00230         case pat::CandKinResolution::MCPInvSpher:  
00231             ret[0] = 1.0/p4.P(); ret[1] = p4.theta(); ret[2] = p4.phi(); ret[3] = p4.mass();
00232             break;
00233         // ======= HEP CYLINDRICAL ==========
00234         case pat::CandKinResolution::EtThetaPhi: 
00235             ret[0] = p4.pt(); ret[1] = p4.theta(); ret[2] = p4.phi(); ret[3] = 0;
00236             break;
00237         case pat::CandKinResolution::EtEtaPhi:
00238             ret[0] = p4.pt(); ret[1] = p4.eta(); ret[2] = p4.phi();   ret[3] = 0;
00239             break;
00240         // ======= DEVIATES ==========
00241         case pat::CandKinResolution::MomDev:
00242         case pat::CandKinResolution::EMomDev:
00243              ret[0] = 1.0; ret[1] = 0.0; ret[2] = 0.0; ret[3] = 1.0;
00244              break;
00245         case pat::CandKinResolution::MCMomDev:
00246              ret[0] = 1.0; ret[1] = 0.0; ret[2] = 0.0; ret[3] = p4.mass();
00247              break;
00248         case pat::CandKinResolution::EScaledMomDev:
00249              ret[0] = 1.0; ret[1] = 0.0; ret[2] = 0.0; ret[3] = p4.E()/p4.P();
00250              break;
00251         // ======= OTHER ==========
00252         case pat::CandKinResolution::Invalid:
00253             throw cms::Exception("Invalid parametrization") << parametrization;
00254         default:
00255             throw cms::Exception("Not Implemented") << "getResolEta not yet implemented for parametrization " << parametrization ;
00256     }
00257 }
00258 
00259 void 
00260 pat::helper::ParametrizationHelper::setParametersFromP4(pat::CandKinResolution::Parametrization parametrization, 
00261         AlgebraicVector4 &v, const math::PtEtaPhiMLorentzVector &p4) {
00262     setParametersFromAnyVector(parametrization, v, p4);
00263 }
00264 
00265 void 
00266 pat::helper::ParametrizationHelper::setParametersFromP4(pat::CandKinResolution::Parametrization parametrization, 
00267         AlgebraicVector4 &v, const math::XYZTLorentzVector &p4) {
00268     setParametersFromAnyVector(parametrization, v, p4);
00269 }
00270 
00271 AlgebraicVector4 
00272 pat::helper::ParametrizationHelper::parametersFromP4(pat::CandKinResolution::Parametrization parametrization, const math::PtEtaPhiMLorentzVector &p4) {
00273     AlgebraicVector4 ret;
00274     setParametersFromP4(parametrization, ret, p4);
00275     return ret;
00276 }
00277 
00278 AlgebraicVector4 
00279 pat::helper::ParametrizationHelper::parametersFromP4(pat::CandKinResolution::Parametrization parametrization, const math::XYZTLorentzVector &p4) {
00280     AlgebraicVector4 ret;
00281     setParametersFromP4(parametrization, ret, p4);
00282     return ret;
00283 }
00284 
00285 AlgebraicVector4 
00286 pat::helper::ParametrizationHelper::diffToParameters(pat::CandKinResolution::Parametrization parametrization, 
00287                 const math::PtEtaPhiMLorentzVector &p4ini, const math::PtEtaPhiMLorentzVector &p4fin) 
00288 {
00289     AlgebraicVector4 ret;
00290     switch (parametrization) {
00291         case pat::CandKinResolution::Cart:
00292         case pat::CandKinResolution::ECart:    
00293         case pat::CandKinResolution::MCCart:
00294             ret = parametersFromP4(parametrization,p4fin) - parametersFromP4(parametrization,p4ini);
00295             break;
00296         case pat::CandKinResolution::Spher: 
00297         case pat::CandKinResolution::ESpher:
00298         case pat::CandKinResolution::MCSpher:
00299         case pat::CandKinResolution::MCPInvSpher:
00300         case pat::CandKinResolution::EtThetaPhi: 
00301         case pat::CandKinResolution::EtEtaPhi:
00302             ret = parametersFromP4(parametrization,p4fin) - parametersFromP4(parametrization,p4ini);
00303             while(ret[2] > +M_PI) ret[2] -= (2*M_PI);
00304             while(ret[2] < -M_PI) ret[2] += (2*M_PI);
00305             break;
00306         case pat::CandKinResolution::MCMomDev:
00307         case pat::CandKinResolution::MomDev:
00308         case pat::CandKinResolution::EMomDev:
00309         case pat::CandKinResolution::EScaledMomDev:
00310             return diffToParameters(parametrization,
00311                                     math::XYZTLorentzVector(p4ini),
00312                                     math::XYZTLorentzVector(p4fin));
00313         case pat::CandKinResolution::Invalid:
00314             throw cms::Exception("Invalid parametrization") << parametrization;
00315         default:
00316             throw cms::Exception("Not Implemented") << "diffToParameters not yet implemented for parametrization " << parametrization ;
00317     }
00318     return ret;
00319 }
00320 
00321 
00322 AlgebraicVector4 
00323 pat::helper::ParametrizationHelper::diffToParameters(pat::CandKinResolution::Parametrization parametrization, 
00324                 const math::XYZTLorentzVector &p4ini, const math::XYZTLorentzVector &p4fin) 
00325 {
00326     AlgebraicVector4 ret;
00327     switch (parametrization) {
00328         case pat::CandKinResolution::Cart:
00329         case pat::CandKinResolution::ECart:    
00330         case pat::CandKinResolution::MCCart:
00331         case pat::CandKinResolution::Spher: 
00332             ret = parametersFromP4(parametrization,p4fin) - parametersFromP4(parametrization,p4ini);
00333             break;
00334         case pat::CandKinResolution::ESpher:
00335         case pat::CandKinResolution::MCSpher:
00336         case pat::CandKinResolution::MCPInvSpher:
00337         case pat::CandKinResolution::EtThetaPhi: 
00338         case pat::CandKinResolution::EtEtaPhi:
00339             ret = parametersFromP4(parametrization,p4fin) - parametersFromP4(parametrization,p4ini);
00340             while(ret[2] > +M_PI) ret[2] -= (2*M_PI);
00341             while(ret[2] < -M_PI) ret[2] += (2*M_PI);
00342             break;
00343         case pat::CandKinResolution::MCMomDev:
00344         case pat::CandKinResolution::MomDev:
00345         case pat::CandKinResolution::EMomDev:
00346         case pat::CandKinResolution::EScaledMomDev:
00347             {
00348                 typedef ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D<double> > V3Cart;
00349                 V3Cart p1 = p4ini.Vect(), p2 = p4fin.Vect();
00350                 V3Cart ur = p1.Unit(); 
00351                 V3Cart uz(0,0,1);
00352                 V3Cart uph = uz.Cross(ur).Unit();
00353                 V3Cart uth = ur.Cross(uph).Unit();
00354                 ret[0] = p2.Dot(ur)/p1.R() - 1.0;
00355                 ret[1] = (p2 - p1).Dot(uth);
00356                 ret[2] = (p2 - p1).Dot(uph);
00357                 if (parametrization == pat::CandKinResolution::MomDev) {
00358                     ret[3] = p4fin.mass()/p4ini.mass() - 1.0;
00359                 } else if (parametrization == pat::CandKinResolution::EMomDev) {
00360                     ret[3] = p4fin.energy()/p4ini.energy() - 1.0;
00361                 }
00362             }
00363             break;
00364         case pat::CandKinResolution::Invalid:
00365             throw cms::Exception("Invalid parametrization") << parametrization;
00366         default:
00367             throw cms::Exception("Not Implemented") << "diffToParameters not yet implemented for parametrization " << parametrization ;
00368     }
00369     return ret;
00370 }
00371 
00372 bool
00373 pat::helper::ParametrizationHelper::isAlwaysMassless(pat::CandKinResolution::Parametrization parametrization) 
00374 {
00375     switch (parametrization) {
00376         case pat::CandKinResolution::Cart:
00377         case pat::CandKinResolution::ECart:    
00378         case pat::CandKinResolution::MCCart:
00379         case pat::CandKinResolution::Spher: 
00380         case pat::CandKinResolution::ESpher:
00381         case pat::CandKinResolution::MCSpher:
00382         case pat::CandKinResolution::MCPInvSpher:
00383         case pat::CandKinResolution::MCMomDev:
00384         case pat::CandKinResolution::MomDev:
00385         case pat::CandKinResolution::EMomDev:
00386         case pat::CandKinResolution::EScaledMomDev:
00387             return false;
00388         case pat::CandKinResolution::EtThetaPhi: 
00389         case pat::CandKinResolution::EtEtaPhi:
00390             return true;
00391         case pat::CandKinResolution::Invalid:
00392             throw cms::Exception("Invalid parametrization") << parametrization;
00393         default:
00394             throw cms::Exception("Not Implemented") << "isAlwaysMassless not yet implemented for parametrization " << parametrization ;
00395     }
00396 }
00397 
00398 bool
00399 pat::helper::ParametrizationHelper::isAlwaysMassive(pat::CandKinResolution::Parametrization parametrization) 
00400 {
00401     switch (parametrization) {
00402         case pat::CandKinResolution::Cart:
00403         case pat::CandKinResolution::ECart:    
00404         case pat::CandKinResolution::MCCart:
00405         case pat::CandKinResolution::Spher: 
00406         case pat::CandKinResolution::ESpher:
00407         case pat::CandKinResolution::MCSpher:
00408         case pat::CandKinResolution::MCPInvSpher:
00409         case pat::CandKinResolution::MCMomDev:
00410         case pat::CandKinResolution::EMomDev:
00411         case pat::CandKinResolution::EScaledMomDev:
00412         case pat::CandKinResolution::EtThetaPhi: 
00413         case pat::CandKinResolution::EtEtaPhi:
00414             return false;
00415         case pat::CandKinResolution::MomDev:
00416             return true;
00417         case pat::CandKinResolution::Invalid:
00418             throw cms::Exception("Invalid parametrization") << parametrization;
00419         default:
00420             throw cms::Exception("Not Implemented") << "isAlwaysMassless not yet implemented for parametrization " << parametrization ;
00421     }
00422 }
00423 
00424 bool
00425 pat::helper::ParametrizationHelper::isMassConstrained(pat::CandKinResolution::Parametrization parametrization) 
00426 {
00427     switch (parametrization) {
00428         case pat::CandKinResolution::Cart:
00429         case pat::CandKinResolution::ECart:    
00430         case pat::CandKinResolution::Spher: 
00431         case pat::CandKinResolution::ESpher:
00432         case pat::CandKinResolution::EMomDev:
00433         case pat::CandKinResolution::EScaledMomDev:
00434         case pat::CandKinResolution::EtThetaPhi: 
00435         case pat::CandKinResolution::EtEtaPhi:
00436         case pat::CandKinResolution::MomDev:
00437             return false;
00438         case pat::CandKinResolution::MCCart:
00439         case pat::CandKinResolution::MCSpher:
00440         case pat::CandKinResolution::MCPInvSpher:
00441         case pat::CandKinResolution::MCMomDev:
00442             return true;
00443         case pat::CandKinResolution::Invalid:
00444             throw cms::Exception("Invalid parametrization") << parametrization;
00445         default:
00446             throw cms::Exception("Not Implemented") << "isAlwaysMassless not yet implemented for parametrization " << parametrization ;
00447     }
00448 }
00449 
00450 bool
00451 pat::helper::ParametrizationHelper::isPhysical(pat::CandKinResolution::Parametrization parametrization, 
00452             const AlgebraicVector4 &parameters,
00453             const math::PtEtaPhiMLorentzVector &initialP4) {
00454     switch (parametrization) {
00455         // ======= CARTESIAN ==========
00456         case pat::CandKinResolution::Cart:
00457             return parameters[3] >= 0; // M >= 0
00458         case pat::CandKinResolution::MCCart:
00459             return true;
00460         case pat::CandKinResolution::ECart: 
00461             return (parameters[0]*parameters[0] + parameters[1]*parameters[1] + parameters[2]*parameters[2] <= parameters[3]*parameters[3]); // E >= P
00462         case pat::CandKinResolution::Spher:  
00463             return (parameters[0] >= 0   ) && // P >= 0
00464                    (parameters[3] >= 0   ) && // M >= 0
00465                    (parameters[1] >= 0   ) && // theta >= 0
00466                    (parameters[1] <= M_PI);   // theta <= pi
00467         case pat::CandKinResolution::MCSpher:    
00468             return (parameters[0] >= 0   ) && // P >= 0
00469                    (parameters[1] >= 0   ) && // theta >= 0
00470                    (parameters[1] <= M_PI);   // theta <= pi
00471         case pat::CandKinResolution::ESpher:        //  
00472             return (parameters[0] >= 0            ) && // P >= 0
00473                    (parameters[3] >= parameters[0]) && // E >= P
00474                    (parameters[1] >= 0            ) && // theta >= 0
00475                    (parameters[1] <= M_PI         ) ;  // theta <= PI
00476        case pat::CandKinResolution::MCPInvSpher:   //  
00477             return (parameters[0] >  0   ) && // 1/P > 0
00478                    (parameters[1] >= 0   ) && // theta >= 0
00479                    (parameters[1] <= M_PI);   // theta <= pi
00480         // ======= HEP CYLINDRICAL ==========
00481         case pat::CandKinResolution::EtThetaPhi: 
00482             return (parameters[0] >  0   ) && // Et >= 0
00483                    (parameters[1] >= 0   ) && // theta >= 0
00484                    (parameters[1] <= M_PI);   // theta <= pi
00485         case pat::CandKinResolution::EtEtaPhi:      // as simple as that
00486             return (parameters[0] >  0); // Et >= 0
00487         // ======= MomentumDeviates ==========
00488         case pat::CandKinResolution::MomDev:
00489             return (parameters[0] >= 0) && // P >= 0
00490                    (parameters[3] >= 0);   // m/M0 >= 0
00491         case pat::CandKinResolution::MCMomDev:
00492             return (parameters[0] >= 0); // P >= 0
00493         case pat::CandKinResolution::EMomDev:
00494         case pat::CandKinResolution::EScaledMomDev:
00495             {
00496                 if (parameters[0] <= 0) return false;
00497                 ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D<double> > p = initialP4.Vect(), uz(0,0,1), uph, uth;
00498                 uph = uz.Cross(p).Unit();
00499                 uth = p.Cross(uph).Unit();
00500                 p *= parameters[0]; 
00501                 p += uth * parameters[1] + uph * parameters[2];
00502                 if (parametrization == pat::CandKinResolution::EMomDev) {
00503                     if (parameters[3] < 0) return false;
00504                     double m2 = ROOT::Math::Square(parameters[3] * initialP4.energy()) - p.Mag2();
00505                     if (m2 < 0) return false;
00506                 } else if (parametrization == pat::CandKinResolution::EScaledMomDev) {
00507                     if (parameters[3] < 0) return false;
00508                     double m2 = ROOT::Math::Square(p.R()*initialP4.E()/initialP4.P()) - p.Mag2();
00509                     if (m2 < 0) return false;
00510                 }
00511                 return true;
00512             }
00513         // ======= OTHER ==========
00514         case pat::CandKinResolution::Invalid:
00515             throw cms::Exception("Invalid parametrization") << parametrization;
00516         default:
00517             throw cms::Exception("Not Implemented") << "getResolEta not yet implemented for parametrization " << parametrization ;
00518     }
00519 }
00520 
00521
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/DataFormats/PatCandidates/src/ParametrizationHelper.cc
