File indexing completed on 2023-03-17 10:51:04
0001 #include <cppunit/extensions/HelperMacros.h>
0002 #include <algorithm>
0003 #include <iterator>
0004 #include <iostream>
0005 #include <iomanip>
0006
0007 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0008
0009 class testPackedCandidate : public CppUnit::TestFixture {
0010 CPPUNIT_TEST_SUITE(testPackedCandidate);
0011
0012 CPPUNIT_TEST(testDefaultConstructor);
0013 CPPUNIT_TEST(testCopyConstructor);
0014 CPPUNIT_TEST(testPackUnpack);
0015 CPPUNIT_TEST(testSimulateReadFromRoot);
0016 CPPUNIT_TEST(testPackUnpackTime);
0017 CPPUNIT_TEST(testQualityFlags);
0018
0019 CPPUNIT_TEST_SUITE_END();
0020
0021 public:
0022 void setUp() {}
0023 void tearDown() {}
0024
0025 void testDefaultConstructor();
0026 void testCopyConstructor();
0027 void testPackUnpack();
0028 void testSimulateReadFromRoot();
0029
0030 void testPackUnpackTime();
0031 void testQualityFlags();
0032
0033 private:
0034 };
0035
0036 CPPUNIT_TEST_SUITE_REGISTRATION(testPackedCandidate);
0037
0038 void testPackedCandidate::testDefaultConstructor() {
0039 pat::PackedCandidate pc;
0040
0041 CPPUNIT_ASSERT(pc.polarP4() == pat::PackedCandidate::PolarLorentzVector(0, 0, 0, 0));
0042 CPPUNIT_ASSERT(pc.p4() == pat::PackedCandidate::LorentzVector(0, 0, 0, 0));
0043 CPPUNIT_ASSERT(pc.vertex() == pat::PackedCandidate::Point(0, 0, 0));
0044 }
0045
0046 static bool tolerance(double iLHS, double iRHS, double fraction) {
0047 return std::abs(iLHS - iRHS) <= fraction * std::abs(iLHS + iRHS) / 2.;
0048 }
0049
0050 void testPackedCandidate::testCopyConstructor() {
0051 pat::PackedCandidate::LorentzVector lv(1., 0.5, 0., std::sqrt(1. + 0.25 + 0.120 * 0.120));
0052 pat::PackedCandidate::PolarLorentzVector plv(lv.Pt(), lv.Eta(), lv.Phi(), lv.M());
0053
0054 pat::PackedCandidate::Point v(0.01, 0.02, 0.);
0055
0056
0057 pat::PackedCandidate pc(lv, v, 1., 1., 1., 11, reco::VertexRefProd(), reco::VertexRef().key());
0058
0059
0060
0061
0062
0063
0064 pat::PackedCandidate copy_pc(pc);
0065
0066
0067
0068
0069
0070 CPPUNIT_ASSERT(©_pc.polarP4() != &pc.polarP4());
0071 CPPUNIT_ASSERT(©_pc.p4() != &pc.p4());
0072 CPPUNIT_ASSERT(©_pc.vertex() != &pc.vertex());
0073
0074 CPPUNIT_ASSERT(tolerance(pc.vertex().X(), copy_pc.vertex().X(), 0.001));
0075 CPPUNIT_ASSERT(tolerance(pc.vertex().Y(), copy_pc.vertex().Y(), 0.001));
0076 CPPUNIT_ASSERT(tolerance(pc.vertex().Z(), copy_pc.vertex().Z(), 0.01));
0077 CPPUNIT_ASSERT(tolerance(pc.dxy(), copy_pc.dxy(), 0.001));
0078 CPPUNIT_ASSERT(tolerance(pc.dzAssociatedPV(), copy_pc.dzAssociatedPV(), 0.01));
0079 CPPUNIT_ASSERT(tolerance(pc.phiAtVtx(), copy_pc.phiAtVtx(), 0.001));
0080 CPPUNIT_ASSERT(tolerance(pc.etaAtVtx(), copy_pc.etaAtVtx(), 0.001));
0081 CPPUNIT_ASSERT(tolerance(pc.ptTrk(), copy_pc.ptTrk(), 0.001));
0082
0083 {
0084
0085 pat::PackedCandidate def;
0086 pat::PackedCandidate copy_def(def);
0087
0088 CPPUNIT_ASSERT(copy_def.polarP4() == pat::PackedCandidate::PolarLorentzVector(0, 0, 0, 0));
0089 CPPUNIT_ASSERT(copy_def.p4() == pat::PackedCandidate::LorentzVector(0, 0, 0, 0));
0090 CPPUNIT_ASSERT(copy_def.vertex() == pat::PackedCandidate::Point(0, 0, 0));
0091 }
0092 }
0093
0094 void testPackedCandidate::testPackUnpack() {
0095 pat::PackedCandidate::LorentzVector lv(1., 1., 0., std::sqrt(2. + 0.120 * 0.120));
0096 pat::PackedCandidate::PolarLorentzVector plv(lv.Pt(), lv.Eta(), lv.Phi(), lv.M());
0097
0098 pat::PackedCandidate::Point v(-0.005, 0.005, 0.1);
0099 float trkPt = plv.Pt() + 0.5;
0100 float trkEta = plv.Eta() - 0.1;
0101 float trkPhi = -3. / 4. * 3.1416;
0102
0103
0104 pat::PackedCandidate pc(lv, v, trkPt, trkEta, trkPhi, 11, reco::VertexRefProd(), reco::VertexRef().key());
0105
0106 pc.pack(true);
0107 pc.packVtx(true);
0108
0109 CPPUNIT_ASSERT(tolerance(pc.polarP4().Pt(), plv.Pt(), 0.001));
0110 CPPUNIT_ASSERT(tolerance(pc.polarP4().Eta(), plv.Eta(), 0.001));
0111 CPPUNIT_ASSERT(tolerance(pc.polarP4().Phi(), plv.Phi(), 0.001));
0112 CPPUNIT_ASSERT(tolerance(pc.polarP4().M(), plv.M(), 0.001));
0113 CPPUNIT_ASSERT(tolerance(pc.p4().X(), lv.X(), 0.001));
0114 CPPUNIT_ASSERT(tolerance(pc.p4().Y(), lv.Y(), 0.001));
0115 CPPUNIT_ASSERT(tolerance(pc.p4().Z(), lv.Z(), 0.001));
0116 CPPUNIT_ASSERT(tolerance(pc.p4().E(), lv.E(), 0.001));
0117 CPPUNIT_ASSERT(tolerance(pc.vertex().X(), v.X(), 0.001));
0118 CPPUNIT_ASSERT(tolerance(pc.vertex().Y(), v.Y(), 0.001));
0119 CPPUNIT_ASSERT(tolerance(pc.vertex().Z(), v.Z(), 0.01));
0120
0121
0122
0123
0124 CPPUNIT_ASSERT(tolerance(pc.ptTrk(), trkPt, 0.001));
0125 CPPUNIT_ASSERT(tolerance(pc.etaAtVtx(), trkEta, 0.001));
0126 CPPUNIT_ASSERT(tolerance(pc.phiAtVtx(), trkPhi, 0.001));
0127
0128 pc.setP4(plv * 2);
0129 CPPUNIT_ASSERT(tolerance(pc.ptTrk(), trkPt, 0.001));
0130 CPPUNIT_ASSERT(tolerance(pc.etaAtVtx(), trkEta, 0.001));
0131 CPPUNIT_ASSERT(tolerance(pc.phiAtVtx(), trkPhi, 0.001));
0132 pc.setP4(lv * 3);
0133 CPPUNIT_ASSERT(tolerance(pc.ptTrk(), trkPt, 0.001));
0134 CPPUNIT_ASSERT(tolerance(pc.etaAtVtx(), trkEta, 0.001));
0135 CPPUNIT_ASSERT(tolerance(pc.phiAtVtx(), trkPhi, 0.001));
0136 pc.setPz(pc.p4().Pz() + 3.3);
0137 CPPUNIT_ASSERT(tolerance(pc.ptTrk(), trkPt, 0.005));
0138 CPPUNIT_ASSERT(tolerance(pc.etaAtVtx(), trkEta, 0.005));
0139 CPPUNIT_ASSERT(tolerance(pc.phiAtVtx(), trkPhi, 0.005));
0140 }
0141
0142 void testPackedCandidate::testSimulateReadFromRoot() {
0143 pat::PackedCandidate::LorentzVector lv(1., 1., 0., std::sqrt(2. + 0.120 * 0.120));
0144 pat::PackedCandidate::PolarLorentzVector plv(lv.Pt(), lv.Eta(), lv.Phi(), lv.M());
0145 pat::PackedCandidate::Point v(-0.005, 0.005, 0.1);
0146
0147 float trkPt = plv.Pt() + 0.5;
0148 float trkEta = plv.Eta() - 0.1;
0149 float trkPhi = -3. / 4. * 3.1416;
0150
0151
0152 pat::PackedCandidate pc(lv, v, trkPt, trkEta, trkPhi, 11, reco::VertexRefProd(), reco::VertexRef().key());
0153
0154
0155
0156
0157
0158
0159
0160 delete pc.p4_.exchange(nullptr);
0161 delete pc.p4c_.exchange(nullptr);
0162 delete pc.vertex_.exchange(nullptr);
0163 delete pc.track_.exchange(nullptr);
0164 delete pc.m_.exchange(nullptr);
0165
0166 CPPUNIT_ASSERT(tolerance(pc.polarP4().Pt(), plv.Pt(), 0.001));
0167 CPPUNIT_ASSERT(tolerance(pc.polarP4().Eta(), plv.Eta(), 0.001));
0168 CPPUNIT_ASSERT(tolerance(pc.polarP4().Phi(), plv.Phi(), 0.001));
0169 CPPUNIT_ASSERT(tolerance(pc.polarP4().M(), plv.M(), 0.001));
0170 CPPUNIT_ASSERT(tolerance(pc.p4().X(), lv.X(), 0.001));
0171 CPPUNIT_ASSERT(tolerance(pc.p4().Y(), lv.Y(), 0.001));
0172 CPPUNIT_ASSERT(tolerance(pc.p4().Z(), lv.Z(), 0.001));
0173 CPPUNIT_ASSERT(tolerance(pc.p4().E(), lv.E(), 0.001));
0174 CPPUNIT_ASSERT(tolerance(pc.vertex().X(), v.X(), 0.001));
0175 CPPUNIT_ASSERT(tolerance(pc.vertex().Y(), v.Y(), 0.001));
0176 CPPUNIT_ASSERT(tolerance(pc.vertex().Z(), v.Z(), 0.01));
0177
0178
0179
0180
0181 CPPUNIT_ASSERT(tolerance(pc.ptTrk(), trkPt, 0.001));
0182 CPPUNIT_ASSERT(tolerance(pc.etaAtVtx(), trkEta, 0.001));
0183 CPPUNIT_ASSERT(tolerance(pc.phiAtVtx(), trkPhi, 0.001));
0184 }
0185
0186 void testPackedCandidate::testPackUnpackTime() {
0187 bool debug =
0188 false;
0189
0190 if (debug)
0191 std::cout << std::endl;
0192 if (debug)
0193 std::cout << "Minimum time error: " << pat::PackedCandidate::unpackTimeError(1) << std::endl;
0194 if (debug)
0195 std::cout << "Maximum time error: " << pat::PackedCandidate::unpackTimeError(255) << std::endl;
0196 float avgres = 0;
0197 int navg = 0;
0198 for (int i = 2; i < 255; ++i) {
0199 float unp = pat::PackedCandidate::unpackTimeError(i);
0200 float res = 0.5 * (pat::PackedCandidate::unpackTimeError(i + 1) - pat::PackedCandidate::unpackTimeError(i - 1));
0201 avgres += (res / unp);
0202 navg++;
0203 if (debug)
0204 std::cout << " i = " << i << " unp = " << unp << " quant error = " << res
0205 << " packed = " << int(pat::PackedCandidate::packTimeError(unp + 0.3 * res)) << " and "
0206 << int(pat::PackedCandidate::packTimeError(unp - 0.3 * res)) << std::endl;
0207 CPPUNIT_ASSERT(pat::PackedCandidate::packTimeError(unp + 0.3 * res) == i);
0208 CPPUNIT_ASSERT(pat::PackedCandidate::packTimeError(unp - 0.3 * res) == i);
0209 }
0210 if (debug)
0211 std::cout << "Average rel uncertainty: " << (avgres / navg) << std::endl;
0212 if (debug)
0213 std::cout << std::endl;
0214
0215 if (debug)
0216 std::cout << "Zero time standalone (pos): " << pat::PackedCandidate::unpackTimeNoError(0) << std::endl;
0217 if (debug)
0218 std::cout << "Minimum time standalone (pos): " << pat::PackedCandidate::unpackTimeNoError(+1) << std::endl;
0219 if (debug)
0220 std::cout << "Minimum time standalone (neg): " << pat::PackedCandidate::unpackTimeNoError(-1) << std::endl;
0221 if (debug)
0222 std::cout << "Maximum time standalone, 8 bits (pos): " << pat::PackedCandidate::unpackTimeNoError(+255)
0223 << std::endl;
0224 if (debug)
0225 std::cout << "Maximum time standalone, 8 bits (neg): " << pat::PackedCandidate::unpackTimeNoError(-255)
0226 << std::endl;
0227 if (debug)
0228 std::cout << "Maximum time standalone, 10 bits (pos): " << pat::PackedCandidate::unpackTimeNoError(+1023)
0229 << std::endl;
0230 if (debug)
0231 std::cout << "Maximum time standalone, 10 bits (neg): " << pat::PackedCandidate::unpackTimeNoError(-1023)
0232 << std::endl;
0233 if (debug)
0234 std::cout << "Maximum time standalone, 11 bits (pos): " << pat::PackedCandidate::unpackTimeNoError(+2047)
0235 << std::endl;
0236 if (debug)
0237 std::cout << "Maximum time standalone, 11 bits (neg): " << pat::PackedCandidate::unpackTimeNoError(-2047)
0238 << std::endl;
0239 avgres = 0;
0240 navg = 0;
0241 for (int i = 2; i < 2040; i *= 1.5) {
0242 float unp = pat::PackedCandidate::unpackTimeNoError(i);
0243 float res = 0.5 * (pat::PackedCandidate::unpackTimeNoError(i + 1) - pat::PackedCandidate::unpackTimeNoError(i - 1));
0244 avgres += (res / unp);
0245 navg++;
0246 if (debug)
0247 std::cout << " i = +" << i << " unp = +" << unp << " quant error = " << res
0248 << " packed = " << int(pat::PackedCandidate::packTimeNoError(unp + 0.3 * res)) << " and +"
0249 << int(pat::PackedCandidate::packTimeNoError(unp - 0.3 * res)) << std::endl;
0250 CPPUNIT_ASSERT(pat::PackedCandidate::packTimeNoError(unp + 0.3 * res) == i);
0251 CPPUNIT_ASSERT(pat::PackedCandidate::packTimeNoError(unp - 0.3 * res) == i);
0252 unp = pat::PackedCandidate::unpackTimeNoError(-i);
0253 res = 0.5 * (pat::PackedCandidate::unpackTimeNoError(-i + 1) - pat::PackedCandidate::unpackTimeNoError(-i - 1));
0254 avgres += std::abs(res / unp);
0255 navg++;
0256 if (debug)
0257 std::cout << " i = " << -i << " unp = " << unp << " quant error = " << res
0258 << " packed = " << int(pat::PackedCandidate::packTimeNoError(unp + 0.3 * res)) << " and "
0259 << int(pat::PackedCandidate::packTimeNoError(unp - 0.3 * res)) << std::endl;
0260 CPPUNIT_ASSERT(pat::PackedCandidate::packTimeNoError(unp + 0.3 * res) == -i);
0261 CPPUNIT_ASSERT(pat::PackedCandidate::packTimeNoError(unp - 0.3 * res) == -i);
0262 }
0263 if (debug)
0264 std::cout << "Average rel uncertainty: " << (avgres / navg) << std::endl;
0265 if (debug)
0266 std::cout << std::endl;
0267
0268 for (float aTimeErr = 2.0e-3; aTimeErr <= 1000e-3; aTimeErr *= std::sqrt(5.f)) {
0269 uint8_t packedTimeErr = pat::PackedCandidate::packTimeError(aTimeErr);
0270 float unpackedTimeErr = pat::PackedCandidate::unpackTimeError(packedTimeErr);
0271 if (debug)
0272 std::cout << "For a timeError of " << aTimeErr << " ns (uint8: " << unsigned(packedTimeErr) << ", unpack "
0273 << unpackedTimeErr << ")" << std::endl;
0274 if (debug)
0275 std::cout << "Minimum time (pos): " << pat::PackedCandidate::unpackTimeWithError(+1, packedTimeErr) << std::endl;
0276 if (debug)
0277 std::cout << "Minimum time (neg): " << pat::PackedCandidate::unpackTimeWithError(-1, packedTimeErr) << std::endl;
0278 if (debug)
0279 std::cout << "Maximum time 8 bits (pos): " << pat::PackedCandidate::unpackTimeWithError(+254, packedTimeErr)
0280 << std::endl;
0281 if (debug)
0282 std::cout << "Maximum time 8 bits (neg): " << pat::PackedCandidate::unpackTimeWithError(-254, packedTimeErr)
0283 << std::endl;
0284 if (debug)
0285 std::cout << "Maximum time 10 bits (pos): " << pat::PackedCandidate::unpackTimeWithError(+1022, packedTimeErr)
0286 << std::endl;
0287 if (debug)
0288 std::cout << "Maximum time 10 bits (neg): " << pat::PackedCandidate::unpackTimeWithError(-1022, packedTimeErr)
0289 << std::endl;
0290 if (debug)
0291 std::cout << "Maximum time 12 bits (pos): " << pat::PackedCandidate::unpackTimeWithError(+4094, packedTimeErr)
0292 << std::endl;
0293 if (debug)
0294 std::cout << "Maximum time 12 bits (neg): " << pat::PackedCandidate::unpackTimeWithError(-4094, packedTimeErr)
0295 << std::endl;
0296 avgres = 0;
0297 navg = 0;
0298 for (int i = 2; i < 4096; i *= 3) {
0299 float unp = pat::PackedCandidate::unpackTimeWithError(i, packedTimeErr);
0300 float res = 0.5 * (pat::PackedCandidate::unpackTimeWithError(i + 2, packedTimeErr) -
0301 pat::PackedCandidate::unpackTimeWithError(i - 2, packedTimeErr));
0302 avgres += (res);
0303 navg++;
0304 if (debug)
0305 std::cout << " i = +" << i << " unp = +" << unp << " quant error = " << res << " packed = +"
0306 << int(pat::PackedCandidate::packTimeWithError(unp + 0.2 * res, unpackedTimeErr)) << " and +"
0307 << int(pat::PackedCandidate::packTimeWithError(unp - 0.2 * res, unpackedTimeErr)) << std::endl;
0308 CPPUNIT_ASSERT(pat::PackedCandidate::packTimeWithError(unp + 0.2 * res, unpackedTimeErr) == i);
0309 CPPUNIT_ASSERT(pat::PackedCandidate::packTimeWithError(unp - 0.2 * res, unpackedTimeErr) == i);
0310 unp = pat::PackedCandidate::unpackTimeWithError(-i, packedTimeErr);
0311 res = 0.5 * (pat::PackedCandidate::unpackTimeWithError(-i + 2, packedTimeErr) -
0312 pat::PackedCandidate::unpackTimeWithError(-i - 2, packedTimeErr));
0313 avgres += std::abs(res);
0314 navg++;
0315 if (debug)
0316 std::cout << " i = " << -i << " unp = " << unp << " quant error = " << res
0317 << " packed = " << int(pat::PackedCandidate::packTimeWithError(unp + 0.2 * res, unpackedTimeErr))
0318 << " and " << int(pat::PackedCandidate::packTimeWithError(unp - 0.2 * res, unpackedTimeErr))
0319 << std::endl;
0320 CPPUNIT_ASSERT(pat::PackedCandidate::packTimeWithError(unp + 0.2 * res, unpackedTimeErr) == -i);
0321 CPPUNIT_ASSERT(pat::PackedCandidate::packTimeWithError(unp - 0.2 * res, unpackedTimeErr) == -i);
0322 }
0323 if (debug)
0324 std::cout << "Average abs uncertainty: " << (avgres / navg) << std::endl;
0325 if (debug)
0326 std::cout << "Now testing overflows: " << std::endl;
0327 avgres = 0;
0328 navg = 0;
0329 for (float aTime = aTimeErr; aTime <= 200; aTime *= std::sqrt(5.f)) {
0330 int i = pat::PackedCandidate::packTimeWithError(aTime, unpackedTimeErr);
0331 float res = 0.5 * std::abs(pat::PackedCandidate::unpackTimeWithError(i + 2, packedTimeErr) -
0332 pat::PackedCandidate::unpackTimeWithError(i - 2, packedTimeErr));
0333 float unp = pat::PackedCandidate::unpackTimeWithError(i, packedTimeErr);
0334 if (debug)
0335 std::cout << " t = +" << aTime << " i = +" << i << " quant error = " << res << " unpacked = +" << unp
0336 << " diff/res = " << (aTime - unp) / res << std::endl;
0337 CPPUNIT_ASSERT(std::abs(unp - aTime) < res);
0338 avgres = std::max(avgres, std::abs(res / unp));
0339 i = pat::PackedCandidate::packTimeWithError(-aTime, unpackedTimeErr);
0340 res = 0.5 * std::abs(pat::PackedCandidate::unpackTimeWithError(i + 2, packedTimeErr) -
0341 pat::PackedCandidate::unpackTimeWithError(i - 2, packedTimeErr));
0342 unp = pat::PackedCandidate::unpackTimeWithError(i, packedTimeErr);
0343 if (debug)
0344 std::cout << " t = " << -aTime << " i = " << i << " quant error = " << res << " unpacked = " << unp
0345 << " diff/res = " << (-aTime - unp) / res << std::endl;
0346 CPPUNIT_ASSERT(std::abs(unp + aTime) < res);
0347 avgres = std::max(avgres, std::abs(res / unp));
0348 }
0349 if (debug)
0350 std::cout << "Worst rel uncertainty: " << (avgres) << std::endl;
0351 if (debug)
0352 std::cout << std::endl;
0353 }
0354 if (debug)
0355 std::cout << std::endl;
0356 }
0357
0358 void testPackedCandidate::testQualityFlags() {
0359 const std::vector<pat::PackedCandidate::PVAssociationQuality> pvAssocVals = {
0360 pat::PackedCandidate::NotReconstructedPrimary,
0361 pat::PackedCandidate::OtherDeltaZ,
0362 pat::PackedCandidate::CompatibilityBTag,
0363 pat::PackedCandidate::CompatibilityDz,
0364 pat::PackedCandidate::UsedInFitLoose,
0365 pat::PackedCandidate::UsedInFitTight};
0366 const std::vector<pat::PackedCandidate::LostInnerHits> lostHitsVals = {
0367 pat::PackedCandidate::validHitInFirstPixelBarrelLayer,
0368 pat::PackedCandidate::noLostInnerHits,
0369 pat::PackedCandidate::oneLostInnerHit,
0370 pat::PackedCandidate::moreLostInnerHits};
0371 const std::vector<bool> trackQualVals = {false, true};
0372 const std::vector<bool> glbMuonVals = {false, true};
0373 const std::vector<bool> staMuonVals = {false, true};
0374 const std::vector<bool> goodEGVals = {false, true};
0375
0376 pat::PackedCandidate cand;
0377
0378 for (auto pvAssoc : pvAssocVals) {
0379 for (auto lostHits : lostHitsVals) {
0380 for (auto trackQual : trackQualVals) {
0381 for (auto glbMuon : glbMuonVals) {
0382 for (auto staMuon : staMuonVals) {
0383 for (auto goodEGamma : goodEGVals) {
0384 cand.setMuonID(staMuon, glbMuon);
0385 cand.setGoodEgamma(goodEGamma);
0386 cand.setAssociationQuality(pvAssoc);
0387 cand.setLostInnerHits(lostHits);
0388 cand.setTrackHighPurity(trackQual);
0389
0390 CPPUNIT_ASSERT(lostHits == cand.lostInnerHits());
0391 CPPUNIT_ASSERT(goodEGamma == cand.isGoodEgamma());
0392 CPPUNIT_ASSERT(staMuon == cand.isStandAloneMuon());
0393 CPPUNIT_ASSERT(glbMuon == cand.isGlobalMuon());
0394 CPPUNIT_ASSERT(trackQual == cand.trackHighPurity());
0395 CPPUNIT_ASSERT(pvAssoc == cand.pvAssociationQuality());
0396 }
0397 }
0398 }
0399 }
0400 }
0401 }
0402 }