PatCandidates/test/testKinParametrizations.cc

0001 #include <cppunit/extensions/HelperMacros.h>
0002 #include <algorithm>
0003 #include <iterator>
0004 #include <iostream>
0005
0006 #include "DataFormats/PatCandidates/interface/ParametrizationHelper.h"
0007 #include "DataFormats/Math/interface/deltaPhi.h"
0008
0009 class testKinParametrizations : public CppUnit::TestFixture {
0010   CPPUNIT_TEST_SUITE(testKinParametrizations);
0011
0012   CPPUNIT_TEST(testTrivialVec2Par_Cart);
0013   CPPUNIT_TEST(testTrivialVec2Par_ECart);
0014   CPPUNIT_TEST(testTrivialVec2Par_Spher);
0015   CPPUNIT_TEST(testTrivialVec2Par_ESpher);
0016   CPPUNIT_TEST(testTrivialVec2Par_MomDev);
0017   CPPUNIT_TEST(testTrivialVec2Par_EMomDev);
0018   CPPUNIT_TEST(testTrivialVec2Par_MCCart);
0019   CPPUNIT_TEST(testTrivialVec2Par_MCSpher);
0020   CPPUNIT_TEST(testTrivialVec2Par_MCPInvSpher);
0021   CPPUNIT_TEST(testTrivialVec2Par_EtEtaPhi);
0022   CPPUNIT_TEST(testTrivialVec2Par_EtThetaPhi);
0023   CPPUNIT_TEST(testTrivialVec2Par_MCMomDev);
0024   CPPUNIT_TEST(testTrivialVec2Par_EScaledMomDev);
0025
0026   CPPUNIT_TEST(testVecDiff2Par_Cart);
0027   CPPUNIT_TEST(testVecDiff2Par_ECart);
0028   CPPUNIT_TEST(testVecDiff2Par_Spher);
0029   CPPUNIT_TEST(testVecDiff2Par_ESpher);
0030   CPPUNIT_TEST(testVecDiff2Par_MomDev);
0031   CPPUNIT_TEST(testVecDiff2Par_EMomDev);
0032   CPPUNIT_TEST(testVecDiff2Par_MCCart);
0033   CPPUNIT_TEST(testVecDiff2Par_MCSpher);
0034   CPPUNIT_TEST(testVecDiff2Par_MCPInvSpher);
0035   CPPUNIT_TEST(testVecDiff2Par_EtEtaPhi);
0036   CPPUNIT_TEST(testVecDiff2Par_EtThetaPhi);
0037   CPPUNIT_TEST(testVecDiff2Par_MCMomDev);
0038   CPPUNIT_TEST(testVecDiff2Par_EScaledMomDev);
0039
0040   CPPUNIT_TEST(testVecVec2Diff_Cart);
0041   CPPUNIT_TEST(testVecVec2Diff_ECart);
0042   CPPUNIT_TEST(testVecVec2Diff_Spher);
0043   CPPUNIT_TEST(testVecVec2Diff_ESpher);
0044   CPPUNIT_TEST(testVecVec2Diff_MomDev);
0045   CPPUNIT_TEST(testVecVec2Diff_EMomDev);
0046   CPPUNIT_TEST(testVecVec2Diff_MCCart);
0047   CPPUNIT_TEST(testVecVec2Diff_MCSpher);
0048   CPPUNIT_TEST(testVecVec2Diff_MCPInvSpher);
0049   CPPUNIT_TEST(testVecVec2Diff_EtEtaPhi);
0050   CPPUNIT_TEST(testVecVec2Diff_EtThetaPhi);
0051   CPPUNIT_TEST(testVecVec2Diff_MCMomDev);
0052   CPPUNIT_TEST(testVecVec2Diff_EScaledMomDev);
0053
0054   CPPUNIT_TEST_SUITE_END();
0055
0056 public:
0057   void setUp() {}
0058   void tearDown() {}
0059
0060   void testTrivialVec2Par_Cart();
0061   void testTrivialVec2Par_ECart();
0062   void testTrivialVec2Par_Spher();
0063   void testTrivialVec2Par_ESpher();
0064   void testTrivialVec2Par_MomDev();
0065   void testTrivialVec2Par_EMomDev();
0066   void testTrivialVec2Par_MCCart();
0067   void testTrivialVec2Par_MCSpher();
0068   void testTrivialVec2Par_MCPInvSpher();
0069   void testTrivialVec2Par_EtEtaPhi();
0070   void testTrivialVec2Par_EtThetaPhi();
0071   void testTrivialVec2Par_MCMomDev();
0072   void testTrivialVec2Par_EScaledMomDev();
0073
0074   void testVecDiff2Par_Cart();
0075   void testVecDiff2Par_ECart();
0076   void testVecDiff2Par_Spher();
0077   void testVecDiff2Par_ESpher();
0078   void testVecDiff2Par_MomDev();
0079   void testVecDiff2Par_EMomDev();
0080   void testVecDiff2Par_MCCart();
0081   void testVecDiff2Par_MCSpher();
0082   void testVecDiff2Par_MCPInvSpher();
0083   void testVecDiff2Par_EtEtaPhi();
0084   void testVecDiff2Par_EtThetaPhi();
0085   void testVecDiff2Par_MCMomDev();
0086   void testVecDiff2Par_EScaledMomDev();
0087
0088   void testVecVec2Diff_Cart();
0089   void testVecVec2Diff_ECart();
0090   void testVecVec2Diff_Spher();
0091   void testVecVec2Diff_ESpher();
0092   void testVecVec2Diff_MomDev();
0093   void testVecVec2Diff_EMomDev();
0094   void testVecVec2Diff_MCCart();
0095   void testVecVec2Diff_MCSpher();
0096   void testVecVec2Diff_MCPInvSpher();
0097   void testVecVec2Diff_EtEtaPhi();
0098   void testVecVec2Diff_EtThetaPhi();
0099   void testVecVec2Diff_MCMomDev();
0100   void testVecVec2Diff_EScaledMomDev();
0101
0102   typedef math::XYZTLorentzVector P4C;
0103   typedef math::PtEtaPhiMLorentzVector P4P;
0104   typedef AlgebraicVector4 V4;
0105
0106 private:
0107   typedef pat::CandKinResolution::Parametrization Parametrization;
0108   bool testTrivialVec2Par(Parametrization p, int tries = 100000);
0109   bool testVecDiff2Par(Parametrization p, int tries = 100000);
0110   bool testVecVec2Diff(Parametrization p, int tries = 100000);
0111
0112   // Utilities
0113   static double r(double val = 1.0) { return rand() * val / double(RAND_MAX); }
0114   static double r(double from, double to) { return rand() * (to - from) / double(RAND_MAX) + from; }
0115   static P4P r4(double m = -1) {
0116     double mass = (m != -1 ? m : (r() > .3 ? r(.1, 30) : 0));
0117     return P4P(r(5, 25), r(-2.5, 5), r(-M_PI, M_PI), mass);
0118   }
0119   static P4P p4near(const P4P &p, double d = 0.2) {
0120     double mass = (p.mass() == 0 ? 0 : p.mass() * r(1 - d, 1 + d));
0121     return P4P(p.pt() * r(1 - d, 1 + d), p.eta() + r(-d, +d), p.phi() + r(-d, +d), mass);
0122   }
0123   static V4 v4near(const V4 &v4, double d = 0.2) {
0124     V4 ret;
0125     ret[0] = v4[0] * r(1 - d, 1 + d);
0126     ret[1] = v4[1] * r(1 - d, 1 + d);
0127     ret[2] = v4[2] * r(1 - d, 1 + d);
0128     ret[3] = v4[3] * r(1 - d, 1 + d);
0129     return ret;
0130   }
0131   static P4C r4c(double m = -1) { return P4C(r4()); }
0132   static bool d(double x, double y, double eps = 1e-6) {
0133     return std::abs(x - y) < eps * (1 + std::abs(x) + std::abs(y));
0134   }
0135   static bool d4(P4P &p1, P4P &p2, double eps = 1e-6) {
0136     return d(p1.E(), p2.E(), eps) && d(p1.Px(), p2.Px(), eps) && d(p1.Py(), p2.Py(), eps) && d(p1.Pz(), p2.Pz(), eps);
0137   }
0138   static bool d4(P4C &p1, P4C &p2, double eps = 1e-6) {
0139     return d(p1.Pt(), p2.Pt(), eps) && d(p1.Eta(), p2.Eta(), eps) && d(p1.M(), p2.M(), eps) &&
0140            (deltaPhi(p1.Phi(), p2.Py()) < eps);
0141   }
0142   static bool d4(const V4 &v1, const V4 &v2, double eps = 1e-6) {
0143     return d(v1[0], v2[0], eps) && d(v1[1], v2[1], eps) && d(v1[2], v2[2], eps) && d(v1[3], v2[3], eps);
0144   }
0145 };
0146
0147 CPPUNIT_TEST_SUITE_REGISTRATION(testKinParametrizations);
0148
0149 bool testKinParametrizations::testTrivialVec2Par(Parametrization par, int tries) {
0150   using namespace pat::helper::ParametrizationHelper;
0151   srand(37);
0152   for (int i = 0; i < tries; ++i) {
0153     P4P p1 = r4(isAlwaysMassless(par) ? 0 : (r() > .3 ? r(.1, 30) : 0));
0154     V4 v4 = parametersFromP4(par, p1);
0155     P4P p2 = polarP4fromParameters(par, v4, p1);
0156     if (!d4(p1, p2)) {
0157       std::cerr << "\nFailure: (try " << i << "):"
0158                 << "\n  p1 = " << p1 << "\n  p2 = " << p2 << "\n  v4 = " << v4 << std::endl;
0159       return false;
0160     }
0161   }
0162   return true;
0163 }
0164
0165 bool testKinParametrizations::testVecVec2Diff(Parametrization par, int tries) {
0166   using namespace pat::helper::ParametrizationHelper;
0167   srand(37);
0168   for (int i = 0; i < tries; ++i) {
0169     P4P p1 = r4(isAlwaysMassless(par) ? 0 : (isAlwaysMassive(par) ? r(.1, 30) : (r() > .3 ? r(.1, 30) : 0)));
0170     V4 v1 = parametersFromP4(par, p1), v2;
0171     int phystries = 0;
0172     do {
0173       v2 = v4near(v1);
0174       if (dimension(par) == 3)
0175         v2[3] = v1[3];  // keep constraint
0176       phystries++;
0177     } while (!isPhysical(par, v2, p1) && (phystries < 50));
0178     if (!isPhysical(par, v2, p1)) {
0179       std::cerr << "\nFailure (try " << i << "/" << phystries << "):"
0180                 << "\n  p1  = " << p1 << "\n  v1  = " << v1 << "\n  v2  = " << v2 << std::endl;
0181       return false;
0182     }
0183     V4 dv1 = v2 - v1;
0184     P4P p2 = polarP4fromParameters(par, v2, p1);
0185     V4 dv2 = diffToParameters(par, p1, p2);
0186     if (!d4(dv1, dv2)) {
0187       std::cerr << "\nFailure (try " << i << "):"
0188                 << "\n  p1  = " << p1 << "\n  v1  = " << v1 << "\n  dv1 = " << dv1 << "\n  v2  = " << v2
0189                 << "\n  p2  = " << p2 << "\n  dv2 = " << dv2 << std::endl;
0190       return false;
0191     }
0192   }
0193   return true;
0194 }
0195
0196 bool testKinParametrizations::testVecDiff2Par(Parametrization par, int tries) {
0197   using namespace pat::helper::ParametrizationHelper;
0198   srand(37);
0199   for (int i = 0; i < tries; ++i) {
0200     P4P p1 = r4(isAlwaysMassless(par) ? 0 : (isAlwaysMassive(par) ? r(.1, 30) : (r() > .3 ? r(.1, 30) : 0)));
0201     V4 v1 = parametersFromP4(par, p1);
0202     P4P p2 = p4near(p1);
0203     if (isMassConstrained(par)) {
0204       p2.SetM(p1.mass());
0205     } else if (par == pat::CandKinResolution::EScaledMomDev) {
0206       // This is more tricky to get
0207       p2 = P4C(p2.px(), p2.py(), p2.pz(), p1.E() / p1.P() * p2.P());
0208     }
0209     V4 dv = diffToParameters(par, p1, p2);
0210     V4 v2 = v1 + dv;
0211     P4P p3 = polarP4fromParameters(par, v2, p1);
0212     if (!d4(p3, p2)) {
0213       std::cerr << "\nFailure: (try " << i << "):"
0214                 << "\n  p1 = " << p1 << "\n  v1 = " << v1 << "\n  p2 = " << p2 << "\n  dv = " << dv << "\n  v2 = " << v2
0215                 << "\n  p3 = " << p3 << std::endl;
0216       return false;
0217     }
0218   }
0219   return true;
0220 }
0221
0222 void testKinParametrizations::testTrivialVec2Par_Cart() {
0223   CPPUNIT_ASSERT(testTrivialVec2Par(pat::CandKinResolution::Cart));
0224 }
0225 void testKinParametrizations::testTrivialVec2Par_ECart() {
0226   CPPUNIT_ASSERT(testTrivialVec2Par(pat::CandKinResolution::ECart));
0227 }
0228 void testKinParametrizations::testTrivialVec2Par_Spher() {
0229   CPPUNIT_ASSERT(testTrivialVec2Par(pat::CandKinResolution::Spher));
0230 }
0231 void testKinParametrizations::testTrivialVec2Par_ESpher() {
0232   CPPUNIT_ASSERT(testTrivialVec2Par(pat::CandKinResolution::ESpher));
0233 }
0234 void testKinParametrizations::testTrivialVec2Par_MomDev() {
0235   CPPUNIT_ASSERT(testTrivialVec2Par(pat::CandKinResolution::MomDev));
0236 }
0237 void testKinParametrizations::testTrivialVec2Par_EMomDev() {
0238   CPPUNIT_ASSERT(testTrivialVec2Par(pat::CandKinResolution::EMomDev));
0239 }
0240 void testKinParametrizations::testTrivialVec2Par_MCCart() {
0241   CPPUNIT_ASSERT(testTrivialVec2Par(pat::CandKinResolution::MCCart));
0242 }
0243 void testKinParametrizations::testTrivialVec2Par_MCSpher() {
0244   CPPUNIT_ASSERT(testTrivialVec2Par(pat::CandKinResolution::MCSpher));
0245 }
0246 void testKinParametrizations::testTrivialVec2Par_MCPInvSpher() {
0247   CPPUNIT_ASSERT(testTrivialVec2Par(pat::CandKinResolution::MCPInvSpher));
0248 }
0249 void testKinParametrizations::testTrivialVec2Par_EtEtaPhi() {
0250   CPPUNIT_ASSERT(testTrivialVec2Par(pat::CandKinResolution::EtEtaPhi));
0251 }
0252 void testKinParametrizations::testTrivialVec2Par_EtThetaPhi() {
0253   CPPUNIT_ASSERT(testTrivialVec2Par(pat::CandKinResolution::EtThetaPhi));
0254 }
0255 void testKinParametrizations::testTrivialVec2Par_MCMomDev() {
0256   CPPUNIT_ASSERT(testTrivialVec2Par(pat::CandKinResolution::MCMomDev));
0257 }
0258 void testKinParametrizations::testTrivialVec2Par_EScaledMomDev() {
0259   CPPUNIT_ASSERT(testTrivialVec2Par(pat::CandKinResolution::EScaledMomDev));
0260 }
0261
0262 void testKinParametrizations::testVecDiff2Par_Cart() { CPPUNIT_ASSERT(testVecDiff2Par(pat::CandKinResolution::Cart)); }
0263 void testKinParametrizations::testVecDiff2Par_ECart() {
0264   CPPUNIT_ASSERT(testVecDiff2Par(pat::CandKinResolution::ECart));
0265 }
0266 void testKinParametrizations::testVecDiff2Par_Spher() {
0267   CPPUNIT_ASSERT(testVecDiff2Par(pat::CandKinResolution::Spher));
0268 }
0269 void testKinParametrizations::testVecDiff2Par_ESpher() {
0270   CPPUNIT_ASSERT(testVecDiff2Par(pat::CandKinResolution::ESpher));
0271 }
0272 void testKinParametrizations::testVecDiff2Par_MCCart() {
0273   CPPUNIT_ASSERT(testVecDiff2Par(pat::CandKinResolution::MCCart));
0274 }
0275 void testKinParametrizations::testVecDiff2Par_MCSpher() {
0276   CPPUNIT_ASSERT(testVecDiff2Par(pat::CandKinResolution::MCSpher));
0277 }
0278 void testKinParametrizations::testVecDiff2Par_MCPInvSpher() {
0279   CPPUNIT_ASSERT(testVecDiff2Par(pat::CandKinResolution::MCPInvSpher));
0280 }
0281 void testKinParametrizations::testVecDiff2Par_EtEtaPhi() {
0282   CPPUNIT_ASSERT(testVecDiff2Par(pat::CandKinResolution::EtEtaPhi));
0283 }
0284 void testKinParametrizations::testVecDiff2Par_EtThetaPhi() {
0285   CPPUNIT_ASSERT(testVecDiff2Par(pat::CandKinResolution::EtThetaPhi));
0286 }
0287 void testKinParametrizations::testVecDiff2Par_MomDev() {
0288   CPPUNIT_ASSERT(testVecDiff2Par(pat::CandKinResolution::MomDev));
0289 }
0290 void testKinParametrizations::testVecDiff2Par_EMomDev() {
0291   CPPUNIT_ASSERT(testVecDiff2Par(pat::CandKinResolution::EMomDev));
0292 }
0293 void testKinParametrizations::testVecDiff2Par_MCMomDev() {
0294   CPPUNIT_ASSERT(testVecDiff2Par(pat::CandKinResolution::MCMomDev));
0295 }
0296 void testKinParametrizations::testVecDiff2Par_EScaledMomDev() {
0297   CPPUNIT_ASSERT(testVecDiff2Par(pat::CandKinResolution::EScaledMomDev));
0298 }
0299
0300 void testKinParametrizations::testVecVec2Diff_Cart() { CPPUNIT_ASSERT(testVecVec2Diff(pat::CandKinResolution::Cart)); }
0301 void testKinParametrizations::testVecVec2Diff_ECart() {
0302   CPPUNIT_ASSERT(testVecVec2Diff(pat::CandKinResolution::ECart));
0303 }
0304 void testKinParametrizations::testVecVec2Diff_Spher() {
0305   CPPUNIT_ASSERT(testVecVec2Diff(pat::CandKinResolution::Spher));
0306 }
0307 void testKinParametrizations::testVecVec2Diff_ESpher() {
0308   CPPUNIT_ASSERT(testVecVec2Diff(pat::CandKinResolution::ESpher));
0309 }
0310 void testKinParametrizations::testVecVec2Diff_MCCart() {
0311   CPPUNIT_ASSERT(testVecVec2Diff(pat::CandKinResolution::MCCart));
0312 }
0313 void testKinParametrizations::testVecVec2Diff_MCSpher() {
0314   CPPUNIT_ASSERT(testVecVec2Diff(pat::CandKinResolution::MCSpher));
0315 }
0316 void testKinParametrizations::testVecVec2Diff_MCPInvSpher() {
0317   CPPUNIT_ASSERT(testVecVec2Diff(pat::CandKinResolution::MCPInvSpher));
0318 }
0319 void testKinParametrizations::testVecVec2Diff_EtEtaPhi() {
0320   CPPUNIT_ASSERT(testVecVec2Diff(pat::CandKinResolution::EtEtaPhi));
0321 }
0322 void testKinParametrizations::testVecVec2Diff_EtThetaPhi() {
0323   CPPUNIT_ASSERT(testVecVec2Diff(pat::CandKinResolution::EtThetaPhi));
0324 }
0325 void testKinParametrizations::testVecVec2Diff_MomDev() {
0326   CPPUNIT_ASSERT(testVecVec2Diff(pat::CandKinResolution::MomDev));
0327 }
0328 void testKinParametrizations::testVecVec2Diff_EMomDev() {
0329   CPPUNIT_ASSERT(testVecVec2Diff(pat::CandKinResolution::EMomDev));
0330 }
0331 void testKinParametrizations::testVecVec2Diff_MCMomDev() {
0332   CPPUNIT_ASSERT(testVecVec2Diff(pat::CandKinResolution::MCMomDev));
0333 }
0334 void testKinParametrizations::testVecVec2Diff_EScaledMomDev() {
0335   CPPUNIT_ASSERT(testVecVec2Diff(pat::CandKinResolution::EScaledMomDev));
0336 }