Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-01-31 03:06:22 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!