Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-05-23 03:13:04 

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