PatCandidates/src/PackedCandidate.cc

0001 #include "DataFormats/Math/interface/deltaPhi.h"
0002 #include "DataFormats/Math/interface/libminifloat.h"
0003 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0004 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0005 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0006
0007 #include "DataFormats/Math/interface/liblogintpack.h"
0008 using namespace logintpack;
0009
0010 CovarianceParameterization pat::PackedCandidate::covarianceParameterization_;
0011 std::once_flag pat::PackedCandidate::covariance_load_flag;
0012
0013 void pat::PackedCandidate::pack(bool unpackAfterwards) {
0014   float unpackedPt = std::min<float>(p4_.load()->Pt(), MiniFloatConverter::max());
0015   packedPt_ = MiniFloatConverter::float32to16(unpackedPt);
0016   packedEta_ = int16_t(std::round(p4_.load()->Eta() / 6.0f * std::numeric_limits<int16_t>::max()));
0017   packedPhi_ = int16_t(std::round(p4_.load()->Phi() / 3.2f * std::numeric_limits<int16_t>::max()));
0018   packedM_ = MiniFloatConverter::float32to16(p4_.load()->M());
0019   if (unpackAfterwards) {
0020     delete p4_.exchange(nullptr);
0021     delete p4c_.exchange(nullptr);
0022     unpack();  // force the values to match with the packed ones
0023   }
0024 }
0025
0026 void pat::PackedCandidate::packVtx(bool unpackAfterwards) {
0027   reco::VertexRef pvRef = vertexRef();
0028   Point pv = pvRef.isNonnull() ? pvRef->position() : Point();
0029   float dxPV = vertex_.load()->X() - pv.X(),
0030         dyPV = vertex_.load()->Y() - pv.Y();  //, rPV = std::hypot(dxPV, dyPV);
0031   float s = std::sin(float(p4_.load()->Phi()) + dphi_),
0032         c = std::cos(float(p4_.load()->Phi() + dphi_));  // not the fastest option, but we're in reduced
0033                                                          // precision already, so let's avoid more roundoffs
0034   dxy_ = -dxPV * s + dyPV * c;
0035   // if we want to go back to the full x,y,z we need to store also
0036   // float dl = dxPV * c + dyPV * s;
0037   // float xRec = - dxy_ * s + dl * c, yRec = dxy_ * c + dl * s;
0038   float pzpt = p4_.load()->Pz() / p4_.load()->Pt();
0039   dz_ = vertex_.load()->Z() - pv.Z() - (dxPV * c + dyPV * s) * pzpt;
0040   packedDxy_ = MiniFloatConverter::float32to16(dxy_ * 100);
0041   packedDz_ = pvRef.isNonnull() ? MiniFloatConverter::float32to16(dz_ * 100)
0042                                 : int16_t(std::round(dz_ / 40.f * std::numeric_limits<int16_t>::max()));
0043   packedDPhi_ = int16_t(std::round(dphi_ / 3.2f * std::numeric_limits<int16_t>::max()));
0044   packedDEta_ = MiniFloatConverter::float32to16(deta_);
0045   packedDTrkPt_ = MiniFloatConverter::float32to16(dtrkpt_);
0046
0047   if (unpackAfterwards) {
0048     delete vertex_.exchange(nullptr);
0049     unpackVtx();
0050   }
0051 }
0052
0053 void pat::PackedCandidate::unpack() const {
0054   float pt = MiniFloatConverter::float16to32(packedPt_);
0055   double shift = (pt < 1. ? 0.1 * pt : 0.1 / pt);    // shift particle phi to break
0056                                                      // degeneracies in angular separations
0057   double sign = ((int(pt * 10) % 2 == 0) ? 1 : -1);  // introduce a pseudo-random sign of the shift
0058   double phi = int16_t(packedPhi_) * 3.2f / std::numeric_limits<int16_t>::max() +
0059                sign * shift * 3.2 / std::numeric_limits<int16_t>::max();
0060   auto p4 = std::make_unique<PolarLorentzVector>(pt,
0061                                                  int16_t(packedEta_) * 6.0f / std::numeric_limits<int16_t>::max(),
0062                                                  phi,
0063                                                  MiniFloatConverter::float16to32(packedM_));
0064   auto p4c = std::make_unique<LorentzVector>(*p4);
0065   PolarLorentzVector *expectp4 = nullptr;
0066   if (p4_.compare_exchange_strong(expectp4, p4.get())) {
0067     p4.release();
0068   }
0069
0070   // p4c_ works as the guard for unpacking so it
0071   // must be set last
0072   LorentzVector *expectp4c = nullptr;
0073   if (p4c_.compare_exchange_strong(expectp4c, p4c.get())) {
0074     p4c.release();
0075   }
0076 }
0077
0078 void pat::PackedCandidate::packCovariance(const reco::TrackBase::CovarianceMatrix &m, bool unpackAfterwards) {
0079   packedCovariance_.dptdpt = packCovarianceElement(m, 0, 0);
0080   packedCovariance_.detadeta = packCovarianceElement(m, 1, 1);
0081   packedCovariance_.dphidphi = packCovarianceElement(m, 2, 2);
0082   packedCovariance_.dxydxy = packCovarianceElement(m, 3, 3);
0083   packedCovariance_.dzdz = packCovarianceElement(m, 4, 4);
0084   packedCovariance_.dxydz = packCovarianceElement(m, 3, 4);
0085   packedCovariance_.dlambdadz = packCovarianceElement(m, 1, 4);
0086   packedCovariance_.dphidxy = packCovarianceElement(m, 2, 3);
0087   // unpack afterwards
0088   if (unpackAfterwards)
0089     unpackCovariance();
0090 }
0091
0092 void pat::PackedCandidate::unpackCovariance() const {
0093   const CovarianceParameterization &p = covarianceParameterization();
0094   if (p.isValid()) {
0095     auto m = std::make_unique<reco::TrackBase::CovarianceMatrix>();
0096     for (int i = 0; i < 5; i++)
0097       for (int j = 0; j < 5; j++) {
0098         (*m)(i, j) = 0;
0099       }
0100     unpackCovarianceElement(*m, packedCovariance_.dptdpt, 0, 0);
0101     unpackCovarianceElement(*m, packedCovariance_.detadeta, 1, 1);
0102     unpackCovarianceElement(*m, packedCovariance_.dphidphi, 2, 2);
0103     unpackCovarianceElement(*m, packedCovariance_.dxydxy, 3, 3);
0104     unpackCovarianceElement(*m, packedCovariance_.dzdz, 4, 4);
0105     unpackCovarianceElement(*m, packedCovariance_.dxydz, 3, 4);
0106     unpackCovarianceElement(*m, packedCovariance_.dlambdadz, 1, 4);
0107     unpackCovarianceElement(*m, packedCovariance_.dphidxy, 2, 3);
0108     reco::TrackBase::CovarianceMatrix *expected = nullptr;
0109     if (m_.compare_exchange_strong(expected, m.get())) {
0110       m.release();
0111     }
0112
0113   } else {
0114     throw edm::Exception(edm::errors::UnimplementedFeature)
0115         << "You do not have a valid track parameters file loaded. "
0116         << "Please check that the release version is compatible with your "
0117            "input data"
0118         << "or avoid accessing track parameter uncertainties. ";
0119   }
0120 }
0121
0122 void pat::PackedCandidate::unpackVtx() const {
0123   reco::VertexRef pvRef = vertexRef();
0124   dphi_ = int16_t(packedDPhi_) * 3.2f / std::numeric_limits<int16_t>::max(),
0125   deta_ = MiniFloatConverter::float16to32(packedDEta_);
0126   dtrkpt_ = MiniFloatConverter::float16to32(packedDTrkPt_);
0127   dxy_ = MiniFloatConverter::float16to32(packedDxy_) / 100.;
0128   dz_ = pvRef.isNonnull() ? MiniFloatConverter::float16to32(packedDz_) / 100.
0129                           : int16_t(packedDz_) * 40.f / std::numeric_limits<int16_t>::max();
0130   Point pv = pvRef.isNonnull() ? pvRef->position() : Point();
0131   float phi = p4_.load()->Phi() + dphi_, s = std::sin(phi), c = std::cos(phi);
0132   auto vertex = std::make_unique<Point>(pv.X() - dxy_ * s,
0133                                         pv.Y() + dxy_ * c,
0134                                         pv.Z() + dz_);  // for our choice of using the PCA to the PV, by definition the
0135                                                         // remaining term -(dx*cos(phi) + dy*sin(phi))*(pz/pt) is zero
0136
0137   Point *expected = nullptr;
0138   if (vertex_.compare_exchange_strong(expected, vertex.get())) {
0139     vertex.release();
0140   }
0141 }
0142
0143 pat::PackedCandidate::~PackedCandidate() {
0144   delete p4_.load();
0145   delete p4c_.load();
0146   delete vertex_.load();
0147   delete track_.load();
0148   delete m_.load();
0149 }
0150
0151 float pat::PackedCandidate::dxy(const Point &p) const {
0152   maybeUnpackBoth();
0153   const float phi = float(p4_.load()->Phi()) + dphi_;
0154   return -(vertex_.load()->X() - p.X()) * std::sin(phi) + (vertex_.load()->Y() - p.Y()) * std::cos(phi);
0155 }
0156 float pat::PackedCandidate::dz(const Point &p) const {
0157   maybeUnpackBoth();
0158   const float phi = float(p4_.load()->Phi()) + dphi_;
0159   const float pzpt = deta_ ? std::sinh(etaAtVtx()) : p4_.load()->Pz() / p4_.load()->Pt();
0160   return (vertex_.load()->Z() - p.Z()) -
0161          ((vertex_.load()->X() - p.X()) * std::cos(phi) + (vertex_.load()->Y() - p.Y()) * std::sin(phi)) * pzpt;
0162 }
0163
0164 void pat::PackedCandidate::unpackTrk() const {
0165   maybeUnpackBoth();
0166   math::RhoEtaPhiVector p3(ptTrk(), etaAtVtx(), phiAtVtx());
0167   maybeUnpackCovariance();
0168   int numberOfStripLayers = stripLayersWithMeasurement(), numberOfPixelLayers = pixelLayersWithMeasurement();
0169   int numberOfPixelHits = this->numberOfPixelHits();
0170   int numberOfHits = this->numberOfHits();
0171
0172   int ndof = numberOfHits + numberOfPixelHits - 5;
0173   LostInnerHits innerLost = lostInnerHits();
0174
0175   auto track = std::make_unique<reco::Track>(normalizedChi2_ * ndof,
0176                                              ndof,
0177                                              *vertex_,
0178                                              math::XYZVector(p3.x(), p3.y(), p3.z()),
0179                                              charge(),
0180                                              *(m_.load()),
0181                                              reco::TrackBase::undefAlgorithm,
0182                                              reco::TrackBase::loose);
0183   int i = 0;
0184   if (firstHit_ == 0) {  // Backward compatible
0185     if (innerLost == validHitInFirstPixelBarrelLayer) {
0186       track->appendTrackerHitPattern(PixelSubdetector::PixelBarrel, 1, 0, TrackingRecHit::valid);
0187       i = 1;
0188     }
0189   } else {
0190     track->appendHitPattern(firstHit_, TrackingRecHit::valid);
0191   }
0192
0193   if (firstHit_ != 0 && reco::HitPattern::pixelHitFilter(firstHit_))
0194     i = 1;
0195
0196   // add hits to match the number of laters and validHitInFirstPixelBarrelLayer
0197   if (innerLost == validHitInFirstPixelBarrelLayer) {
0198     // then to encode the number of layers, we add more hits on distinct layers
0199     // (B2, B3, B4, F1, ...)
0200     for (; i < numberOfPixelLayers; i++) {
0201       if (i <= 3) {
0202         track->appendTrackerHitPattern(PixelSubdetector::PixelBarrel, i + 1, 0, TrackingRecHit::valid);
0203       } else {
0204         track->appendTrackerHitPattern(PixelSubdetector::PixelEndcap, i - 3, 0, TrackingRecHit::valid);
0205       }
0206     }
0207   } else {
0208     // to encode the information on the layers, we add one valid hits per layer
0209     // but skipping PXB1
0210     int iOffset = 0;
0211     if (firstHit_ != 0 && reco::HitPattern::pixelHitFilter(firstHit_)) {
0212       iOffset = reco::HitPattern::getLayer(firstHit_);
0213       if (reco::HitPattern::getSubStructure(firstHit_) == PixelSubdetector::PixelEndcap)
0214         iOffset += 3;
0215     } else {
0216       iOffset = 1;
0217     }
0218     for (; i < numberOfPixelLayers; i++) {
0219       if (i + iOffset <= 2) {
0220         track->appendTrackerHitPattern(PixelSubdetector::PixelBarrel, i + iOffset + 1, 0, TrackingRecHit::valid);
0221       } else {
0222         track->appendTrackerHitPattern(PixelSubdetector::PixelEndcap, i + iOffset - 3 + 1, 0, TrackingRecHit::valid);
0223       }
0224     }
0225   }
0226   // add extra hits (overlaps, etc), all on the first layer with a hit - to
0227   // avoid increasing the layer count
0228   for (; i < numberOfPixelHits; i++) {
0229     if (firstHit_ != 0 && reco::HitPattern::pixelHitFilter(firstHit_)) {
0230       track->appendTrackerHitPattern(reco::HitPattern::getSubStructure(firstHit_),
0231                                      reco::HitPattern::getLayer(firstHit_),
0232                                      0,
0233                                      TrackingRecHit::valid);
0234     } else {
0235       track->appendTrackerHitPattern(PixelSubdetector::PixelBarrel,
0236                                      (innerLost == validHitInFirstPixelBarrelLayer ? 1 : 2),
0237                                      0,
0238                                      TrackingRecHit::valid);
0239     }
0240   }
0241   // now start adding strip layers, putting one hit on each layer so that the
0242   // hitPattern.stripLayersWithMeasurement works. we don't know what the layers
0243   // where, so we just start with TIB (4 layers), then TOB (6 layers), then TEC
0244   // (9) and then TID(3), so that we can get a number of valid strip layers up
0245   // to 4+6+9+3
0246   if (firstHit_ != 0 && reco::HitPattern::stripHitFilter(firstHit_))
0247     i += 1;
0248   int slOffset = 0;
0249   if (firstHit_ != 0 && reco::HitPattern::stripHitFilter(firstHit_)) {
0250     slOffset = reco::HitPattern::getLayer(firstHit_) - 1;
0251     if (reco::HitPattern::getSubStructure(firstHit_) == StripSubdetector::TID)
0252       slOffset += 4;
0253     if (reco::HitPattern::getSubStructure(firstHit_) == StripSubdetector::TOB)
0254       slOffset += 7;
0255     if (reco::HitPattern::getSubStructure(firstHit_) == StripSubdetector::TEC)
0256       slOffset += 13;
0257   }
0258   for (int sl = slOffset; sl < numberOfStripLayers + slOffset; ++sl, ++i) {
0259     if (sl < 4)
0260       track->appendTrackerHitPattern(StripSubdetector::TIB, sl + 1, 1, TrackingRecHit::valid);
0261     else if (sl < 4 + 3)
0262       track->appendTrackerHitPattern(StripSubdetector::TID, (sl - 4) + 1, 1, TrackingRecHit::valid);
0263     else if (sl < 7 + 6)
0264       track->appendTrackerHitPattern(StripSubdetector::TOB, (sl - 7) + 1, 1, TrackingRecHit::valid);
0265     else if (sl < 13 + 9)
0266       track->appendTrackerHitPattern(StripSubdetector::TEC, (sl - 13) + 1, 1, TrackingRecHit::valid);
0267     else
0268       break;  // wtf?
0269   }
0270   // finally we account for extra strip hits beyond the one-per-layer added
0271   // above. we put them all on TIB1, to avoid incrementing the number of
0272   // layersWithMeasurement.
0273   for (; i < numberOfHits; i++) {
0274     if (reco::HitPattern::stripHitFilter(firstHit_)) {
0275       track->appendTrackerHitPattern(reco::HitPattern::getSubStructure(firstHit_),
0276                                      reco::HitPattern::getLayer(firstHit_),
0277                                      1,
0278                                      TrackingRecHit::valid);
0279     } else {
0280       track->appendTrackerHitPattern(StripSubdetector::TIB, 1, 1, TrackingRecHit::valid);
0281     }
0282   }
0283
0284   switch (innerLost) {
0285     case validHitInFirstPixelBarrelLayer:
0286       break;
0287     case noLostInnerHits:
0288       break;
0289     case oneLostInnerHit:
0290       track->appendTrackerHitPattern(PixelSubdetector::PixelBarrel, 1, 0, TrackingRecHit::missing_inner);
0291       break;
0292     case moreLostInnerHits:
0293       track->appendTrackerHitPattern(PixelSubdetector::PixelBarrel, 1, 0, TrackingRecHit::missing_inner);
0294       track->appendTrackerHitPattern(PixelSubdetector::PixelBarrel, 2, 0, TrackingRecHit::missing_inner);
0295       break;
0296   };
0297
0298   if (trackHighPurity())
0299     track->setQuality(reco::TrackBase::highPurity);
0300
0301   reco::Track *expected = nullptr;
0302   if (track_.compare_exchange_strong(expected, track.get())) {
0303     track.release();
0304   }
0305 }
0306
0307 //// Everything below is just trivial implementations of reco::Candidate methods
0308
0309 const reco::CandidateBaseRef &pat::PackedCandidate::masterClone() const {
0310   throw cms::Exception("Invalid Reference") << "this Candidate has no master clone reference."
0311                                             << "Can't call masterClone() method.\n";
0312 }
0313
0314 bool pat::PackedCandidate::hasMasterClone() const { return false; }
0315
0316 bool pat::PackedCandidate::hasMasterClonePtr() const { return false; }
0317
0318 const reco::CandidatePtr &pat::PackedCandidate::masterClonePtr() const {
0319   throw cms::Exception("Invalid Reference") << "this Candidate has no master clone ptr."
0320                                             << "Can't call masterClonePtr() method.\n";
0321 }
0322
0323 size_t pat::PackedCandidate::numberOfDaughters() const { return 0; }
0324
0325 size_t pat::PackedCandidate::numberOfMothers() const { return 0; }
0326
0327 bool pat::PackedCandidate::overlap(const reco::Candidate &o) const {
0328   return p4() == o.p4() && vertex() == o.vertex() && charge() == o.charge();
0329   //  return  p4() == o.p4() && charge() == o.charge();
0330 }
0331
0332 const reco::Candidate *pat::PackedCandidate::daughter(size_type) const { return nullptr; }
0333
0334 const reco::Candidate *pat::PackedCandidate::mother(size_type) const { return nullptr; }
0335
0336 const reco::Candidate *pat::PackedCandidate::daughter(const std::string &) const {
0337   throw edm::Exception(edm::errors::UnimplementedFeature)
0338       << "This Candidate type does not implement daughter(std::string). "
0339       << "Please use CompositeCandidate or NamedCompositeCandidate.\n";
0340 }
0341
0342 reco::Candidate *pat::PackedCandidate::daughter(const std::string &) {
0343   throw edm::Exception(edm::errors::UnimplementedFeature)
0344       << "This Candidate type does not implement daughter(std::string). "
0345       << "Please use CompositeCandidate or NamedCompositeCandidate.\n";
0346 }
0347
0348 reco::Candidate *pat::PackedCandidate::daughter(size_type) { return nullptr; }
0349
0350 double pat::PackedCandidate::vertexChi2() const { return 0; }
0351
0352 double pat::PackedCandidate::vertexNdof() const { return 0; }
0353
0354 double pat::PackedCandidate::vertexNormalizedChi2() const { return 0; }
0355
0356 double pat::PackedCandidate::vertexCovariance(int i, int j) const {
0357   throw edm::Exception(edm::errors::UnimplementedFeature)
0358       << "reco::ConcreteCandidate does not implement vertex covariant "
0359          "matrix.\n";
0360 }
0361
0362 void pat::PackedCandidate::fillVertexCovariance(CovarianceMatrix &err) const {
0363   throw edm::Exception(edm::errors::UnimplementedFeature)
0364       << "reco::ConcreteCandidate does not implement vertex covariant "
0365          "matrix.\n";
0366 }
0367
0368 bool pat::PackedCandidate::longLived() const { return false; }
0369
0370 bool pat::PackedCandidate::massConstraint() const { return false; }
0371
0372 // puppiweight
0373 void pat::PackedCandidate::setPuppiWeight(float p, float p_nolep) {
0374   // Set both weights at once to avoid misconfigured weights if called in the
0375   // wrong order
0376   packedPuppiweight_ = std::numeric_limits<uint8_t>::max() * p;
0377   packedPuppiweightNoLepDiff_ = std::numeric_limits<int8_t>::max() * (p_nolep - p);
0378 }
0379
0380 float pat::PackedCandidate::puppiWeight() const {
0381   return 1.f * packedPuppiweight_ / std::numeric_limits<uint8_t>::max();
0382 }
0383
0384 float pat::PackedCandidate::puppiWeightNoLep() const {
0385   return 1.f * packedPuppiweightNoLepDiff_ / std::numeric_limits<int8_t>::max() +
0386          1.f * packedPuppiweight_ / std::numeric_limits<uint8_t>::max();
0387 }
0388
0389 void pat::PackedCandidate::setRawCaloFraction(float p) {
0390   if (100 * p > std::numeric_limits<uint8_t>::max())
0391     rawCaloFraction_ = std::numeric_limits<uint8_t>::max();  // Set to overflow value
0392   else
0393     rawCaloFraction_ = 100 * p;
0394 }
0395
0396 void pat::PackedCandidate::setRawHcalFraction(float p) { rawHcalFraction_ = 100 * p; }
0397
0398 void pat::PackedCandidate::setCaloFraction(float p) { caloFraction_ = 100 * p; }
0399
0400 void pat::PackedCandidate::setHcalFraction(float p) { hcalFraction_ = 100 * p; }
0401
0402 void pat::PackedCandidate::setIsIsolatedChargedHadron(bool p) { isIsolatedChargedHadron_ = p; }
0403
0404 void pat::PackedCandidate::setDTimeAssociatedPV(float aTime, float aTimeError) {
0405   if (aTime == 0 && aTimeError == 0) {
0406     packedTime_ = 0;
0407     packedTimeError_ = 0;
0408   } else if (aTimeError == 0) {
0409     packedTimeError_ = 0;
0410     packedTime_ = packTimeNoError(aTime);
0411   } else {
0412     packedTimeError_ = packTimeError(aTimeError);
0413     aTimeError = unpackTimeError(packedTimeError_);  // for reproducibility
0414     packedTime_ = packTimeWithError(aTime, aTimeError);
0415   }
0416 }
0417
0418 /// static to allow unit testing
0419 uint8_t pat::PackedCandidate::packTimeError(float timeError) {
0420   if (timeError <= 0)
0421     return 0;
0422   // log-scale packing.
0423   // for MIN_TIMEERROR = 0.002, EXPO_TIMEERROR = 5:
0424   //      minimum value 0.002 = 2ps (packed as 1)
0425   //      maximum value 0.5 ns      (packed as 255)
0426   //      constant *relative* precision of about 2%
0427   return std::max<uint8_t>(
0428       std::min(std::round(std::ldexp(std::log2(timeError / MIN_TIMEERROR), +EXPO_TIMEERROR)), 255.f), 1);
0429 }
0430 float pat::PackedCandidate::unpackTimeError(uint8_t timeError) {
0431   return timeError > 0 ? MIN_TIMEERROR * std::exp2(std::ldexp(float(timeError), -EXPO_TIMEERROR)) : -1.0f;
0432 }
0433 float pat::PackedCandidate::unpackTimeNoError(int16_t time) {
0434   if (time == 0)
0435     return 0.f;
0436   return (time > 0 ? MIN_TIME_NOERROR : -MIN_TIME_NOERROR) *
0437          std::exp2(std::ldexp(float(std::abs(time)), -EXPO_TIME_NOERROR));
0438 }
0439 int16_t pat::PackedCandidate::packTimeNoError(float time) {
0440   // encoding in log scale to store times in a large range with few bits.
0441   // for MIN_TIME_NOERROR = 0.0002 and EXPO_TIME_NOERROR = 6:
0442   //    smallest non-zero time = 0.2 ps (encoded as +/-1)
0443   //    one BX, +/- 12.5 ns, is fully covered with 11 bits (+/- 1023)
0444   //    12 bits cover by far any plausible value (+/-2047 corresponds to about
0445   //    +/- 0.8 ms!) constant *relative* ~1% precision
0446   if (std::abs(time) < MIN_TIME_NOERROR)
0447     return 0;  // prevent underflows
0448   float fpacked = std::ldexp(std::log2(std::abs(time / MIN_TIME_NOERROR)), +EXPO_TIME_NOERROR);
0449   return (time > 0 ? +1 : -1) * std::min(std::round(fpacked), 2047.f);
0450 }
0451 float pat::PackedCandidate::unpackTimeWithError(int16_t time, uint8_t timeError) {
0452   if (time % 2 == 0) {
0453     // no overflow: drop rightmost bit and unpack in units of timeError
0454     return std::ldexp(unpackTimeError(timeError), EXPO_TIME_WITHERROR) * float(time / 2);
0455   } else {
0456     // overflow: drop rightmost bit, unpack using the noError encoding
0457     return pat::PackedCandidate::unpackTimeNoError(time / 2);
0458   }
0459 }
0460 int16_t pat::PackedCandidate::packTimeWithError(float time, float timeError) {
0461   // Encode in units of timeError * 2^EXPO_TIME_WITHERROR (~1.6% if
0462   // EXPO_TIME_WITHERROR = -6) the largest value that can be stored in 14 bits +
0463   // sign bit + overflow bit is about 260 sigmas values larger than that will be
0464   // stored using the no-timeError packing (with less precision). overflows of
0465   // these kinds should happen only for particles that are late arriving,
0466   // out-of-time, or mis-reconstructed, as timeError is O(20ps) and the beam
0467   // spot witdth is O(200ps)
0468   float fpacked = std::round(time / std::ldexp(timeError, EXPO_TIME_WITHERROR));
0469   if (std::abs(fpacked) < 16383.f) {  // 16383 = (2^14 - 1) = largest absolute
0470                                       // value for a signed 15 bit integer
0471     return int16_t(fpacked) * 2;      // make it even, and fit in a signed 16 bit int
0472   } else {
0473     int16_t packed = packTimeNoError(time);    // encode
0474     return packed * 2 + (time > 0 ? +1 : -1);  // make it odd, to signal that there was an overlow
0475   }
0476 }