PackedCandidate.h

Go to the documentation of this file.

#ifndef __DataFormats_PatCandidates_PackedCandidate_h__
#define __DataFormats_PatCandidates_PackedCandidate_h__

#include "DataFormats/Candidate/interface/Candidate.h"
#include "DataFormats/Candidate/interface/CandidateFwd.h"
#include "DataFormats/Common/interface/Association.h"
#include "DataFormats/Common/interface/RefVector.h"
#include "DataFormats/Math/interface/deltaPhi.h"
#include "DataFormats/PatCandidates/interface/CovarianceParameterization.h"
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/VertexReco/interface/VertexFwd.h"
#include "FWCore/Utilities/interface/thread_safety_macros.h"
#include <atomic>
#include <mutex>

/* #include "DataFormats/Math/interface/PtEtaPhiMass.h" */

// forward declare testing structure
class testPackedCandidate;

namespace pat {
class PackedCandidate : public reco::Candidate {
public:
  typedef reco::CandidateCollection daughters;
  typedef math::XYZTLorentzVector LorentzVector;
  typedef math::PtEtaPhiMLorentzVector PolarLorentzVector;
  typedef math::XYZPoint Point;
  typedef math::XYZVector Vector;

  typedef unsigned int index;
  PackedCandidate()
      : packedPt_(0), packedEta_(0), packedPhi_(0), packedM_(0), packedDxy_(0),
        packedDz_(0), packedDPhi_(0), packedDEta_(0), packedDTrkPt_(0),
        packedCovariance_(), packedPuppiweight_(0),
        packedPuppiweightNoLepDiff_(0), rawCaloFraction_(0),
        rawHcalFraction_(0), caloFraction_(0), hcalFraction_(0), packedTime_(0),
        packedTimeError_(0), isIsolatedChargedHadron_(false),
        p4_(new PolarLorentzVector(0, 0, 0, 0)),
        p4c_(new LorentzVector(0, 0, 0, 0)), vertex_(new Point(0, 0, 0)),
        dphi_(0), deta_(0), dtrkpt_(0), track_(nullptr), pdgId_(0),
        qualityFlags_(0), pvRefKey_(reco::VertexRef::invalidKey()), m_(nullptr),
        packedHits_(0), packedLayers_(0), normalizedChi2_(0),
        covarianceVersion_(0), covarianceSchema_(0), firstHit_(0) {}

  explicit PackedCandidate(const reco::Candidate &c,
                           const reco::VertexRefProd &pvRefProd,
                           reco::VertexRef::key_type pvRefKey)
      : packedPuppiweight_(0), packedPuppiweightNoLepDiff_(0),
        rawCaloFraction_(0), rawHcalFraction_(0), caloFraction_(0),
        hcalFraction_(0), packedTime_(0), packedTimeError_(0),
        isIsolatedChargedHadron_(false),
        p4_(new PolarLorentzVector(c.pt(), c.eta(), c.phi(), c.mass())),
        p4c_(new LorentzVector(*p4_)), vertex_(new Point(c.vertex())), dphi_(0),
        deta_(0), dtrkpt_(0), track_(nullptr), pdgId_(c.pdgId()),
        qualityFlags_(0), pvRefProd_(pvRefProd), pvRefKey_(pvRefKey),
        m_(nullptr), packedHits_(0), packedLayers_(0), normalizedChi2_(0),
        covarianceVersion_(0), covarianceSchema_(0), firstHit_(0) {
    packBoth();
  }

  explicit PackedCandidate(const PolarLorentzVector &p4, const Point &vtx,
                           float trkPt, float etaAtVtx, float phiAtVtx,
                           int pdgId, const reco::VertexRefProd &pvRefProd,
                           reco::VertexRef::key_type pvRefKey)
      : packedPuppiweight_(0), packedPuppiweightNoLepDiff_(0),
        rawCaloFraction_(0), rawHcalFraction_(0), caloFraction_(0),
        hcalFraction_(0), packedTime_(0), packedTimeError_(0),
        isIsolatedChargedHadron_(false), p4_(new PolarLorentzVector(p4)),
        p4c_(new LorentzVector(*p4_)), vertex_(new Point(vtx)),
        dphi_(reco::deltaPhi(phiAtVtx, p4_.load()->phi())),
        deta_(std::abs(etaAtVtx - p4_.load()->eta()) >= kMinDEtaToStore_
                  ? etaAtVtx - p4_.load()->eta()
                  : 0.),
        dtrkpt_(std::abs(trkPt - p4_.load()->pt()) >= kMinDTrkPtToStore_
                    ? trkPt - p4_.load()->pt()
                    : 0.),
        track_(nullptr), pdgId_(pdgId), qualityFlags_(0), pvRefProd_(pvRefProd),
        pvRefKey_(pvRefKey), m_(nullptr), packedHits_(0), packedLayers_(0),
        normalizedChi2_(0), covarianceVersion_(0), covarianceSchema_(0),
        firstHit_(0) {
    packBoth();
  }

  explicit PackedCandidate(const LorentzVector &p4, const Point &vtx,
                           float trkPt, float etaAtVtx, float phiAtVtx,
                           int pdgId, const reco::VertexRefProd &pvRefProd,
                           reco::VertexRef::key_type pvRefKey)
      : packedPuppiweight_(0), packedPuppiweightNoLepDiff_(0),
        rawCaloFraction_(0), rawHcalFraction_(0), caloFraction_(0),
        hcalFraction_(0), packedTime_(0), packedTimeError_(0),
        isIsolatedChargedHadron_(false),
        p4_(new PolarLorentzVector(p4.Pt(), p4.Eta(), p4.Phi(), p4.M())),
        p4c_(new LorentzVector(p4)), vertex_(new Point(vtx)),
        dphi_(reco::deltaPhi(phiAtVtx, p4_.load()->phi())),
        deta_(std::abs(etaAtVtx - p4_.load()->eta()) >= kMinDEtaToStore_
                  ? etaAtVtx - p4_.load()->eta()
                  : 0.),
        dtrkpt_(std::abs(trkPt - p4_.load()->pt()) >= kMinDTrkPtToStore_
                    ? trkPt - p4_.load()->pt()
                    : 0.),
        track_(nullptr), pdgId_(pdgId), qualityFlags_(0), pvRefProd_(pvRefProd),
        pvRefKey_(pvRefKey), m_(nullptr), packedHits_(0), packedLayers_(0),
        normalizedChi2_(0), covarianceVersion_(0), covarianceSchema_(0),
        firstHit_(0) {
    packBoth();
  }

  PackedCandidate(const PackedCandidate &iOther)
      : packedPt_(iOther.packedPt_), packedEta_(iOther.packedEta_),
        packedPhi_(iOther.packedPhi_), packedM_(iOther.packedM_),
        packedDxy_(iOther.packedDxy_), packedDz_(iOther.packedDz_),
        packedDPhi_(iOther.packedDPhi_), packedDEta_(iOther.packedDEta_),
        packedDTrkPt_(iOther.packedDTrkPt_),
        packedCovariance_(iOther.packedCovariance_),
        packedPuppiweight_(iOther.packedPuppiweight_),
        packedPuppiweightNoLepDiff_(iOther.packedPuppiweightNoLepDiff_),
        rawCaloFraction_(iOther.rawCaloFraction_),
        rawHcalFraction_(iOther.rawHcalFraction_),
        caloFraction_(iOther.caloFraction_),
        hcalFraction_(iOther.hcalFraction_), packedTime_(iOther.packedTime_),
        packedTimeError_(iOther.packedTimeError_),
        isIsolatedChargedHadron_(iOther.isIsolatedChargedHadron_),
        // Need to trigger unpacking in iOther
        p4_(new PolarLorentzVector(iOther.polarP4())),
        p4c_(new LorentzVector(iOther.p4())),
        vertex_(new Point(iOther.vertex())), dxy_(iOther.dxy_), dz_(iOther.dz_),
        dphi_(iOther.dphi_), deta_(iOther.deta_), dtrkpt_(iOther.dtrkpt_),
        track_(iOther.track_ ? new reco::Track(*iOther.track_) : nullptr),
        pdgId_(iOther.pdgId_), qualityFlags_(iOther.qualityFlags_),
        pvRefProd_(iOther.pvRefProd_), pvRefKey_(iOther.pvRefKey_),
        m_(iOther.m_ ? new reco::TrackBase::CovarianceMatrix(*iOther.m_)
                     : nullptr),
        packedHits_(iOther.packedHits_), packedLayers_(iOther.packedLayers_),
        normalizedChi2_(iOther.normalizedChi2_),
        covarianceVersion_(iOther.covarianceVersion_),
        covarianceSchema_(iOther.covarianceSchema_),
        firstHit_(iOther.firstHit_) {}

  PackedCandidate(PackedCandidate &&iOther)
      : packedPt_(iOther.packedPt_), packedEta_(iOther.packedEta_),
        packedPhi_(iOther.packedPhi_), packedM_(iOther.packedM_),
        packedDxy_(iOther.packedDxy_), packedDz_(iOther.packedDz_),
        packedDPhi_(iOther.packedDPhi_), packedDEta_(iOther.packedDEta_),
        packedDTrkPt_(iOther.packedDTrkPt_),
        packedCovariance_(iOther.packedCovariance_),
        packedPuppiweight_(iOther.packedPuppiweight_),
        packedPuppiweightNoLepDiff_(iOther.packedPuppiweightNoLepDiff_),
        rawCaloFraction_(iOther.rawCaloFraction_),
        rawHcalFraction_(iOther.rawHcalFraction_),
        caloFraction_(iOther.caloFraction_),
        hcalFraction_(iOther.hcalFraction_), packedTime_(iOther.packedTime_),
        packedTimeError_(iOther.packedTimeError_),
        isIsolatedChargedHadron_(iOther.isIsolatedChargedHadron_),
        p4_(iOther.p4_.exchange(nullptr)), p4c_(iOther.p4c_.exchange(nullptr)),
        vertex_(iOther.vertex_.exchange(nullptr)), dxy_(iOther.dxy_),
        dz_(iOther.dz_), dphi_(iOther.dphi_), deta_(iOther.deta_),
        dtrkpt_(iOther.dtrkpt_), track_(iOther.track_.exchange(nullptr)),
        pdgId_(iOther.pdgId_), qualityFlags_(iOther.qualityFlags_),
        pvRefProd_(std::move(iOther.pvRefProd_)), pvRefKey_(iOther.pvRefKey_),
        m_(iOther.m_.exchange(nullptr)), packedHits_(iOther.packedHits_),
        packedLayers_(iOther.packedLayers_),
        normalizedChi2_(iOther.normalizedChi2_),
        covarianceVersion_(iOther.covarianceVersion_),
        covarianceSchema_(iOther.covarianceSchema_),
        firstHit_(iOther.firstHit_) {}

  PackedCandidate &operator=(const PackedCandidate &iOther) {
    if (this == &iOther) {
      return *this;
    }
    packedPt_ = iOther.packedPt_;
    packedEta_ = iOther.packedEta_;
    packedPhi_ = iOther.packedPhi_;
    packedM_ = iOther.packedM_;
    packedDxy_ = iOther.packedDxy_;
    packedDz_ = iOther.packedDz_;
    packedDPhi_ = iOther.packedDPhi_;
    packedDEta_ = iOther.packedDEta_;
    packedDTrkPt_ = iOther.packedDTrkPt_;
    packedCovariance_ = iOther.packedCovariance_;
    packedPuppiweight_ = iOther.packedPuppiweight_;
    packedPuppiweightNoLepDiff_ = iOther.packedPuppiweightNoLepDiff_;
    rawCaloFraction_ = iOther.rawCaloFraction_;
    rawHcalFraction_ = iOther.rawHcalFraction_;
    caloFraction_ = iOther.caloFraction_;
    hcalFraction_ = iOther.hcalFraction_;
    packedTime_ = iOther.packedTime_;
    packedTimeError_ = iOther.packedTimeError_;
    isIsolatedChargedHadron_ = iOther.isIsolatedChargedHadron_;
    // Need to trigger unpacking in iOther
    if (p4_) {
      *p4_ = iOther.polarP4();
    } else {
      p4_.store(new PolarLorentzVector(iOther.polarP4()));
    }
    if (p4c_) {
      *p4c_ = iOther.p4();
    } else {
      p4c_.store(new LorentzVector(iOther.p4()));
    }
    if (vertex_) {
      *vertex_ = iOther.vertex();
    } else {
      vertex_.store(new Point(iOther.vertex()));
    }
    dxy_ = iOther.dxy_;
    dz_ = iOther.dz_;
    dphi_ = iOther.dphi_;
    deta_ = iOther.deta_;
    dtrkpt_ = iOther.dtrkpt_;

    if (!iOther.track_) {
      delete track_.exchange(nullptr);
    } else {
      if (!track_) {
        track_.store(new reco::Track(*iOther.track_));
      } else {
        *track_ = *(iOther.track_);
      }
    }

    pdgId_ = iOther.pdgId_;
    qualityFlags_ = iOther.qualityFlags_;
    pvRefProd_ = iOther.pvRefProd_;
    pvRefKey_ = iOther.pvRefKey_;
    if (!iOther.m_) {
      delete m_.exchange(nullptr);
    } else {
      if (!m_) {
        m_.store(new reco::Track::CovarianceMatrix(*iOther.m_));
      } else {
        *m_ = *(iOther.m_);
      }
    }

    packedHits_ = iOther.packedHits_;
    packedLayers_ = iOther.packedLayers_;
    normalizedChi2_ = iOther.normalizedChi2_;
    covarianceVersion_ = iOther.covarianceVersion_;
    covarianceSchema_ = iOther.covarianceSchema_;
    firstHit_ = iOther.firstHit_;
    return *this;
  }

  PackedCandidate &operator=(PackedCandidate &&iOther) {
    if (this == &iOther) {
      return *this;
    }
    packedPt_ = iOther.packedPt_;
    packedEta_ = iOther.packedEta_;
    packedPhi_ = iOther.packedPhi_;
    packedM_ = iOther.packedM_;
    packedDxy_ = iOther.packedDxy_;
    packedDz_ = iOther.packedDz_;
    packedDPhi_ = iOther.packedDPhi_;
    packedDEta_ = iOther.packedDEta_;
    packedDTrkPt_ = iOther.packedDTrkPt_;
    packedCovariance_ = iOther.packedCovariance_;
    packedPuppiweight_ = iOther.packedPuppiweight_;
    packedPuppiweightNoLepDiff_ = iOther.packedPuppiweightNoLepDiff_;
    rawCaloFraction_ = iOther.rawCaloFraction_;
    rawHcalFraction_ = iOther.rawHcalFraction_;
    caloFraction_ = iOther.caloFraction_;
    hcalFraction_ = iOther.hcalFraction_;
    packedTime_ = iOther.packedTime_;
    packedTimeError_ = iOther.packedTimeError_;
    isIsolatedChargedHadron_ = iOther.isIsolatedChargedHadron_;
    delete p4_.exchange(iOther.p4_.exchange(nullptr));
    delete p4c_.exchange(iOther.p4c_.exchange(nullptr));
    delete vertex_.exchange(iOther.vertex_.exchange(nullptr));
    dxy_ = iOther.dxy_;
    dz_ = iOther.dz_;
    dphi_ = iOther.dphi_;
    deta_ = iOther.deta_;
    dtrkpt_ = iOther.dtrkpt_;
    delete track_.exchange(iOther.track_.exchange(nullptr));
    pdgId_ = iOther.pdgId_;
    qualityFlags_ = iOther.qualityFlags_;
    pvRefProd_ = iOther.pvRefProd_;
    pvRefKey_ = iOther.pvRefKey_;
    delete m_.exchange(iOther.m_.exchange(nullptr));
    packedHits_ = iOther.packedHits_;
    packedLayers_ = iOther.packedLayers_;
    normalizedChi2_ = iOther.normalizedChi2_;
    covarianceVersion_ = iOther.covarianceVersion_;
    covarianceSchema_ = iOther.covarianceSchema_;
    firstHit_ = iOther.firstHit_;
    return *this;
  }

  ~PackedCandidate() override;
  size_t numberOfDaughters() const override;
  const reco::Candidate *daughter(size_type) const override;
  size_t numberOfMothers() const override;
  const reco::Candidate *mother(size_type) const override;
  reco::Candidate *daughter(size_type) override;
  reco::Candidate *daughter(const std::string &s) override;
  const reco::Candidate *daughter(const std::string &s) const override;
  size_t numberOfSourceCandidatePtrs() const override { return 0; }
  reco::CandidatePtr sourceCandidatePtr(size_type i) const override {
    return reco::CandidatePtr();
  }

  int charge() const override {
    switch (abs(pdgId_)) {
    case 211:
      return (pdgId_ > 0) - (pdgId_ < 0);
    case 11:
      return (-1) * (pdgId_ > 0) + (pdgId_ < 0); // e
    case 13:
      return (-1) * (pdgId_ > 0) + (pdgId_ < 0); // mu
    case 15:
      return (-1) * (pdgId_ > 0) + (pdgId_ < 0); // tau
    case 24:
      return (pdgId_ > 0) - (pdgId_ < 0); // W
    default:
      return 0; // FIXME: charge is not defined
    }
  }
  void setCharge(int charge) override {}
  int threeCharge() const override { return charge() * 3; }
  void setThreeCharge(int threecharge) override {}
  const LorentzVector &p4() const override {
    if (!p4c_)
      unpack();
    return *p4c_;
  }
  const PolarLorentzVector &polarP4() const override {
    if (!p4c_)
      unpack();
    return *p4_;
  }
  Vector momentum() const override {
    if (!p4c_)
      unpack();
    return p4c_.load()->Vect();
  }
  Vector boostToCM() const override {
    if (!p4c_)
      unpack();
    return p4c_.load()->BoostToCM();
  }
  double p() const override {
    if (!p4c_)
      unpack();
    return p4c_.load()->P();
  }
  double energy() const override {
    if (!p4c_)
      unpack();
    return p4c_.load()->E();
  }
  double et() const override { return (pt() <= 0) ? 0 : p4c_.load()->Et(); }
  double et2() const override { return (pt() <= 0) ? 0 : p4c_.load()->Et2(); }
  double mass() const override {
    if (!p4c_)
      unpack();
    return p4_.load()->M();
  }
  double massSqr() const override {
    if (!p4c_)
      unpack();
    auto m = p4_.load()->M();
    return m * m;
  }

  double mt() const override {
    if (!p4c_)
      unpack();
    return p4_.load()->Mt();
  }
  double mtSqr() const override {
    if (!p4c_)
      unpack();
    return p4_.load()->Mt2();
  }
  double px() const override {
    if (!p4c_)
      unpack();
    return p4c_.load()->Px();
  }
  double py() const override {
    if (!p4c_)
      unpack();
    return p4c_.load()->Py();
  }
  double pz() const override {
    if (!p4c_)
      unpack();
    return p4c_.load()->Pz();
  }
  double pt() const override {
    if (!p4c_)
      unpack();
    return p4_.load()->Pt();
  }
  double phi() const override {
    if (!p4c_)
      unpack();
    return p4_.load()->Phi();
  }

  virtual double ptTrk() const {
    maybeUnpackBoth();
    return p4_.load()->pt() + dtrkpt_;
  }
  virtual float phiAtVtx() const {
    maybeUnpackBoth();
    float ret = p4_.load()->Phi() + dphi_;
    while (ret > float(M_PI))
      ret -= 2 * float(M_PI);
    while (ret < -float(M_PI))
      ret += 2 * float(M_PI);
    return ret;
  }
  virtual float etaAtVtx() const {
    maybeUnpackBoth();
    return p4_.load()->eta() + deta_;
  }

  double theta() const override {
    if (!p4c_)
      unpack();
    return p4_.load()->Theta();
  }
  double eta() const override {
    if (!p4c_)
      unpack();
    return p4_.load()->Eta();
  }
  double rapidity() const override {
    if (!p4c_)
      unpack();
    return p4_.load()->Rapidity();
  }
  double y() const override {
    if (!p4c_)
      unpack();
    return p4_.load()->Rapidity();
  }
  void setP4(const LorentzVector &p4) override {
    maybeUnpackBoth(); // changing px,py,pz changes also mapping between dxy,dz
                       // and x,y,z
    dphi_ += polarP4().Phi() - p4.Phi();
    deta_ += polarP4().Eta() - p4.Eta();
    dtrkpt_ += polarP4().Pt() - p4.Pt();
    *p4_ = PolarLorentzVector(p4.Pt(), p4.Eta(), p4.Phi(), p4.M());
    packBoth();
  }
  void setP4(const PolarLorentzVector &p4) override {
    maybeUnpackBoth(); // changing px,py,pz changes also mapping between dxy,dz
                       // and x,y,z
    dphi_ += polarP4().Phi() - p4.Phi();
    deta_ += polarP4().Eta() - p4.Eta();
    dtrkpt_ += polarP4().Pt() - p4.Pt();
    *p4_ = p4;
    packBoth();
  }
  void setMass(double m) override {
    if (!p4c_)
      unpack();
    *p4_ = PolarLorentzVector(p4_.load()->Pt(), p4_.load()->Eta(),
                              p4_.load()->Phi(), m);
    pack();
  }
  void setPz(double pz) override {
    maybeUnpackBoth(); // changing px,py,pz changes also mapping between dxy,dz
                       // and x,y,z
    *p4c_ = LorentzVector(p4c_.load()->Px(), p4c_.load()->Py(), pz,
                          p4c_.load()->E());
    dphi_ += polarP4().Phi() - p4c_.load()->Phi();
    deta_ += polarP4().Eta() - p4c_.load()->Eta();
    dtrkpt_ += polarP4().Pt() - p4c_.load()->Pt();
    *p4_ = PolarLorentzVector(p4c_.load()->Pt(), p4c_.load()->Eta(),
                              p4c_.load()->Phi(), p4c_.load()->M());
    packBoth();
  }

  // Note: mask is also the maximum value
  enum trackHitShiftsAndMasks {
    trackPixelHitsMask = 7,
    trackStripHitsMask = 31,
    trackStripHitsShift = 3
  };

  // set number of tracker hits and layers
  virtual void setHits(const reco::Track &tk) {
    // first we count the number of layers with hits
    int numberOfPixelLayers_ = tk.hitPattern().pixelLayersWithMeasurement();
    if (numberOfPixelLayers_ > trackPixelHitsMask)
      numberOfPixelLayers_ = trackPixelHitsMask;
    int numberOfStripLayers_ = tk.hitPattern().stripLayersWithMeasurement();
    if (numberOfStripLayers_ > trackStripHitsMask)
      numberOfStripLayers_ = trackStripHitsMask;
    packedLayers_ = (numberOfPixelLayers_ & trackPixelHitsMask) |
                    (numberOfStripLayers_ << trackStripHitsShift);
    // now we count number of additional hits, beyond the one-per-layer implied
    // by packedLayers_
    int numberOfPixelHits_ =
        tk.hitPattern().numberOfValidPixelHits() - numberOfPixelLayers_;
    if (numberOfPixelHits_ > trackPixelHitsMask)
      numberOfPixelHits_ = trackPixelHitsMask;
    int numberOfStripHits_ = tk.hitPattern().numberOfValidHits() -
                             numberOfPixelHits_ - numberOfPixelLayers_ -
                             numberOfStripLayers_;
    if (numberOfStripHits_ > trackStripHitsMask)
      numberOfStripHits_ = trackStripHitsMask;

    packedHits_ = (numberOfPixelHits_ & trackPixelHitsMask) |
                  (numberOfStripHits_ << trackStripHitsShift);
  }

  virtual void
  setTrackProperties(const reco::Track &tk,
                     const reco::Track::CovarianceMatrix &covariance,
                     int quality, int covarianceVersion) {
    covarianceVersion_ = covarianceVersion;
    covarianceSchema_ = quality;
    normalizedChi2_ = tk.normalizedChi2();
    setHits(tk);
    packBoth();
    packCovariance(covariance, false);
  }

  // set track properties using quality and covarianceVersion to define the
  // level of details in the cov. matrix
  virtual void setTrackProperties(const reco::Track &tk, int quality,
                                  int covarianceVersion) {
    setTrackProperties(tk, tk.covariance(), quality, covarianceVersion);
  }

  int numberOfPixelHits() const {
    return (packedHits_ & trackPixelHitsMask) + pixelLayersWithMeasurement();
  }
  int numberOfHits() const {
    return (packedHits_ >> trackStripHitsShift) + stripLayersWithMeasurement() +
           numberOfPixelHits();
  }
  int pixelLayersWithMeasurement() const {
    return packedLayers_ & trackPixelHitsMask;
  }
  int stripLayersWithMeasurement() const {
    return (packedLayers_ >> trackStripHitsShift);
  }
  int trackerLayersWithMeasurement() const {
    return stripLayersWithMeasurement() + pixelLayersWithMeasurement();
  }
  virtual void setCovarianceVersion(int v) { covarianceVersion_ = v; }
  int covarianceVersion() const { return covarianceVersion_; }
  int covarianceSchema() const { return covarianceSchema_; }

  const Point &vertex() const override {
    maybeUnpackBoth();
    return *vertex_;
  } //{ if (fromPV_) return Point(0,0,0); else return Point(0,0,100); }
  double vx() const override {
    maybeUnpackBoth();
    return vertex_.load()->X();
  } //{ return 0; }
  double vy() const override {
    maybeUnpackBoth();
    return vertex_.load()->Y();
  } //{ return 0; }
  double vz() const override {
    maybeUnpackBoth();
    return vertex_.load()->Z();
  } //{ if (fromPV_) return 0; else return 100; }
  void setVertex(const Point &vertex) override {
    maybeUnpackBoth();
    *vertex_ = vertex;
    packVtx();
  }

  enum PVAssoc { NoPV = 0, PVLoose = 1, PVTight = 2, PVUsedInFit = 3 };
  const PVAssoc fromPV(size_t ipv = 0) const {
    reco::VertexRef pvRef = vertexRef();
    if (pvAssociationQuality() == UsedInFitTight and pvRef.key() == ipv)
      return PVUsedInFit;
    if (pvRef.key() == ipv or abs(pdgId()) == 13 or abs(pdgId()) == 11)
      return PVTight;
    if (pvAssociationQuality() == CompatibilityBTag and
        std::abs(dzAssociatedPV()) > std::abs(dz(ipv)))
      return PVTight; // it is not closest, but at least prevents the B
                      // assignment stealing
    if (pvAssociationQuality() < UsedInFitLoose or pvRef->ndof() < 4.0)
      return PVLoose;
    return NoPV;
  }

  enum PVAssociationQuality {
    NotReconstructedPrimary = 0,
    OtherDeltaZ = 1,
    CompatibilityBTag = 4,
    CompatibilityDz = 5,
    UsedInFitLoose = 6,
    UsedInFitTight = 7
  };
  const PVAssociationQuality pvAssociationQuality() const {
    return PVAssociationQuality((qualityFlags_ & assignmentQualityMask) >>
                                assignmentQualityShift);
  }
  void setAssociationQuality(PVAssociationQuality q) {
    qualityFlags_ = (qualityFlags_ & ~assignmentQualityMask) |
                    ((q << assignmentQualityShift) & assignmentQualityMask);
  }

  const reco::VertexRef vertexRef() const {
    return reco::VertexRef(pvRefProd_, pvRefKey_);
  }

  virtual float dxy() const {
    maybeUnpackBoth();
    return dxy_;
  }
  virtual float dz(size_t ipv = 0) const {
    maybeUnpackBoth();
    return dz_ + (*pvRefProd_)[pvRefKey_].position().z() -
           (*pvRefProd_)[ipv].position().z();
  }
  virtual float dzAssociatedPV() const {
    maybeUnpackBoth();
    return dz_;
  }
  virtual float dxy(const Point &p) const;
  virtual float dz(const Point &p) const;

  float dzError() const override {
    maybeUnpackCovariance();
    return sqrt((*m_.load())(4, 4));
  }
  float dxyError() const override {
    maybeUnpackCovariance();
    return sqrt((*m_.load())(3, 3));
  }

  virtual const reco::Track &pseudoTrack() const {
    if (!track_)
      unpackTrk();
    return *track_;
  }

  const reco::Track *bestTrack() const override {
    if (packedHits_ != 0 || packedLayers_ != 0) {
      maybeUnpackTrack();
      return track_.load();
    } else
      return nullptr;
  }
  bool hasTrackDetails() const {
    return (packedHits_ != 0 || packedLayers_ != 0);
  }

  bool trackHighPurity() const {
    return (qualityFlags_ & trackHighPurityMask) >> trackHighPurityShift;
  }
  void setTrackHighPurity(bool highPurity) {
    qualityFlags_ =
        (qualityFlags_ & ~trackHighPurityMask) |
        ((highPurity << trackHighPurityShift) & trackHighPurityMask);
  }

  enum LostInnerHits {
    validHitInFirstPixelBarrelLayer = -1,
    noLostInnerHits = 0, // it could still not have a hit in the first layer,
                         // e.g. if it crosses an inactive sensor
    oneLostInnerHit = 1,
    moreLostInnerHits = 2
  };
  LostInnerHits lostInnerHits() const {
    return LostInnerHits(
        int16_t((qualityFlags_ & lostInnerHitsMask) >> lostInnerHitsShift) - 1);
  }
  void setLostInnerHits(LostInnerHits hits) {
    int lost = hits;
    if (lost > 2)
      lost = 2; // protection against misuse
    lost++;     // shift so it's 0 .. 3 instead of (-1) .. 2
    qualityFlags_ = (qualityFlags_ & ~lostInnerHitsMask) |
                    ((lost << lostInnerHitsShift) & lostInnerHitsMask);
  }

  void setFirstHit(uint16_t pattern) { firstHit_ = pattern; }
  uint16_t firstHit() const { return firstHit_; }

  void setMuonID(bool isStandAlone, bool isGlobal) {
    int16_t muonFlags = isStandAlone | (2 * isGlobal);
    qualityFlags_ = (qualityFlags_ & ~muonFlagsMask) |
                    ((muonFlags << muonFlagsShift) & muonFlagsMask);
  }

  void setGoodEgamma(bool isGoodEgamma = true) {
    int16_t egFlags = (isGoodEgamma << egammaFlagsShift) & egammaFlagsMask;
    qualityFlags_ = (qualityFlags_ & ~egammaFlagsMask) | egFlags;
  }

  int pdgId() const override { return pdgId_; }
  // set PDG identifier
  void setPdgId(int pdgId) override { pdgId_ = pdgId; }
  int status() const override { return qualityFlags_; } /*FIXME*/
  void setStatus(int status) override {} /*FIXME*/
  static const unsigned int longLivedTag = 0; /*FIXME*/
  void setLongLived() override {} /*FIXME*/
  bool longLived() const override;
  static const unsigned int massConstraintTag = 0; /*FIXME*/
  void setMassConstraint() override {} /*FIXME*/
  bool massConstraint() const override;

  PackedCandidate *clone() const override { return new PackedCandidate(*this); }

  double vertexChi2() const override;
  double vertexNdof() const override;
  double vertexNormalizedChi2() const override;
  double vertexCovariance(int i, int j) const override;
  CovarianceMatrix vertexCovariance() const override {
    CovarianceMatrix m;
    fillVertexCovariance(m);
    return m;
  }
  void fillVertexCovariance(CovarianceMatrix &v) const override;
  bool hasMasterClone() const override;
  const reco::CandidateBaseRef &masterClone() const override;
  bool hasMasterClonePtr() const override;

  const reco::CandidatePtr &masterClonePtr() const override;

  template <typename Ref> Ref masterRef() const {
    return masterClone().template castTo<Ref>();
  }

  bool isElectron() const override { return false; }
  bool isMuon() const override { return false; }
  bool isStandAloneMuon() const override {
    return ((qualityFlags_ & muonFlagsMask) >> muonFlagsShift) & 1;
  }
  bool isGlobalMuon() const override {
    return ((qualityFlags_ & muonFlagsMask) >> muonFlagsShift) & 2;
  }
  bool isTrackerMuon() const override { return false; }
  bool isCaloMuon() const override { return false; }
  bool isPhoton() const override { return false; }
  bool isConvertedPhoton() const override { return false; }
  bool isJet() const override { return false; }
  bool isGoodEgamma() const { return (qualityFlags_ & egammaFlagsMask) != 0; }

  // puppiweights
  void setPuppiWeight(float p,
                      float p_nolep = 0.0); 
  float puppiWeight() const;                
  float puppiWeightNoLep() const; 

  // for the neutral fractions
  void
  setRawCaloFraction(float p); 
  float rawCaloFraction() const {
    return (rawCaloFraction_ / 100.);
  } 
  void setRawHcalFraction(
      float p); 
  float rawHcalFraction() const {
    return (rawHcalFraction_ / 100.);
  } 
  void setCaloFraction(
      float p); 
  float caloFraction() const {
    return (caloFraction_ / 100.);
  } 
  void setHcalFraction(float p); 
  float hcalFraction() const {
    return (hcalFraction_ / 100.);
  } 

  // isolated charged hadrons
  void setIsIsolatedChargedHadron(
      bool p); 
  bool isIsolatedChargedHadron() const {
    return isIsolatedChargedHadron_;
  } 

  struct PackedCovariance {
    PackedCovariance()
        : dxydxy(0), dxydz(0), dzdz(0), dlambdadz(0), dphidxy(0), dptdpt(0),
          detadeta(0), dphidphi(0) {}
    // 3D IP covariance
    uint16_t dxydxy;
    uint16_t dxydz;
    uint16_t dzdz;
    // other IP relevant elements
    uint16_t dlambdadz;
    uint16_t dphidxy;
    // other diag elements
    uint16_t dptdpt;
    uint16_t detadeta;
    uint16_t dphidphi;
  };

  virtual float time() const { return vertexRef()->t() + dtimeAssociatedPV(); }
  virtual float dtime(size_t ipv = 0) const {
    return dtimeAssociatedPV() + (*pvRefProd_)[pvRefKey_].t() -
           (*pvRefProd_)[ipv].t();
  }
  virtual float dtimeAssociatedPV() const {
    if (packedTime_ == 0)
      return 0.f;
    if (packedTimeError_ > 0)
      return unpackTimeWithError(packedTime_, packedTimeError_);
    else
      return unpackTimeNoError(packedTime_);
  }
  virtual float timeError() const { return unpackTimeError(packedTimeError_); }
  void setDTimeAssociatedPV(float aTime, float aTimeError = 0);
  void setTime(float aTime, float aTimeError = 0) {
    setDTimeAssociatedPV(aTime - vertexRef()->t(), aTimeError);
  }

private:
  void unpackCovarianceElement(reco::TrackBase::CovarianceMatrix &m,
                               uint16_t packed, int i, int j) const {
    m(i, j) = covarianceParameterization().unpack(
        packed, covarianceSchema_, i, j, pt(), eta(), numberOfHits(),
        numberOfPixelHits());
  }
  uint16_t packCovarianceElement(const reco::TrackBase::CovarianceMatrix &m,
                                 int i, int j) const {
    return covarianceParameterization().pack(m(i, j), covarianceSchema_, i, j,
                                             pt(), eta(), numberOfHits(),
                                             numberOfPixelHits());
  }

protected:
  friend class ::testPackedCandidate;
  static constexpr float kMinDEtaToStore_ = 0.001;
  static constexpr float kMinDTrkPtToStore_ = 0.001;

  uint16_t packedPt_, packedEta_, packedPhi_, packedM_;
  uint16_t packedDxy_, packedDz_, packedDPhi_, packedDEta_, packedDTrkPt_;
  PackedCovariance packedCovariance_;

  void pack(bool unpackAfterwards = true);
  void unpack() const;
  void packVtx(bool unpackAfterwards = true);
  void unpackVtx() const;
  void packCovariance(const reco::TrackBase::CovarianceMatrix &m,
                      bool unpackAfterwards = true);
  void unpackCovariance() const;
  void maybeUnpackBoth() const {
    if (!p4c_)
      unpack();
    if (!vertex_)
      unpackVtx();
  }
  void maybeUnpackTrack() const {
    if (!track_)
      unpackTrk();
  }
  void maybeUnpackCovariance() const {
    if (!m_)
      unpackCovariance();
  }
  void packBoth() {
    pack(false);
    packVtx(false);
    delete p4_.exchange(nullptr);
    delete p4c_.exchange(nullptr);
    delete vertex_.exchange(nullptr);
    unpack();
    unpackVtx();
  } // do it this way, so that we don't loose precision on the angles before
    // computing dxy,dz
  void unpackTrk() const;

  uint8_t packedPuppiweight_;
  int8_t packedPuppiweightNoLepDiff_; // storing the DIFFERENCE of (all - "no
                                      // lep") for compression optimization
  uint8_t rawCaloFraction_;
  int8_t rawHcalFraction_;
  uint8_t caloFraction_;
  int8_t hcalFraction_;
  int16_t packedTime_;
  uint8_t packedTimeError_;

  bool isIsolatedChargedHadron_;

  mutable std::atomic<PolarLorentzVector *> p4_;
  mutable std::atomic<LorentzVector *> p4c_;
  mutable std::atomic<Point *> vertex_;
  CMS_THREAD_GUARD(vertex_) mutable float dxy_, dz_, dphi_, deta_, dtrkpt_;
  mutable std::atomic<reco::Track *> track_;
  int pdgId_;
  uint16_t qualityFlags_;
  reco::VertexRefProd pvRefProd_;
  reco::VertexRef::key_type pvRefKey_;

  mutable std::atomic<reco::TrackBase::CovarianceMatrix *> m_;
  uint8_t packedHits_,
      packedLayers_; // packedLayers_ -> layers with valid hits; packedHits_ ->
                     // extra hits beyond the one-per-layer implied by
                     // packedLayers_

  uint8_t normalizedChi2_;
  uint16_t covarianceVersion_;
  uint16_t covarianceSchema_;
  CMS_THREAD_SAFE static CovarianceParameterization covarianceParameterization_;
  // static std::atomic<CovarianceParameterization*>
  // covarianceParameterization_;
  static std::once_flag covariance_load_flag;
  const CovarianceParameterization &covarianceParameterization() const {
    if (!hasTrackDetails())
      throw edm::Exception(edm::errors::InvalidReference,
                           "Trying to access covariance matrix for a "
                           "PackedCandidate for which it's not available. "
                           "Check hasTrackDetails() before!\n");
    std::call_once(covariance_load_flag,
                   [](int v) { covarianceParameterization_.load(v); },
                   covarianceVersion_);
    if (covarianceParameterization_.loadedVersion() != covarianceVersion_) {
      throw edm::Exception(edm::errors::UnimplementedFeature)
          << "Attempting to load multiple covariance version in same process. "
             "This is not supported.";
    }
    return covarianceParameterization_;
  }

  bool overlap(const reco::Candidate &) const override;
  template <typename, typename, typename> friend struct component;
  friend class ::OverlapChecker;
  friend class ShallowCloneCandidate;
  friend class ShallowClonePtrCandidate;

  enum qualityFlagsShiftsAndMasks {
    assignmentQualityMask = 0x7,
    assignmentQualityShift = 0,
    trackHighPurityMask = 0x8,
    trackHighPurityShift = 3,
    lostInnerHitsMask = 0x30,
    lostInnerHitsShift = 4,
    muonFlagsMask = 0x0600,
    muonFlagsShift = 9,
    egammaFlagsMask = 0x0800,
    egammaFlagsShift = 11
  };

  static uint8_t packTimeError(float timeError);
  static float unpackTimeError(uint8_t timeError);
  static float unpackTimeNoError(int16_t time);
  static int16_t packTimeNoError(float time);
  static float unpackTimeWithError(int16_t time, uint8_t timeError);
  static int16_t packTimeWithError(float time, float timeError);
  static constexpr float MIN_TIMEERROR =
      0.002f; // 2 ps, smallest storable non-zero uncertainty
  static constexpr float MIN_TIME_NOERROR =
      0.0002f; // 0.2 ps, smallest non-zero time that can be stored by
               // packTimeNoError
  static constexpr int EXPO_TIMEERROR =
      5; // power of 2 used in encoding timeError
  static constexpr int EXPO_TIME_NOERROR =
      6; // power of 2 used in encoding time without timeError
  static constexpr int EXPO_TIME_WITHERROR =
      -6; // power of 2 used in encoding time with timeError
public:
  uint16_t firstHit_;
};

typedef std::vector<pat::PackedCandidate> PackedCandidateCollection;
typedef edm::Ref<pat::PackedCandidateCollection> PackedCandidateRef;
typedef edm::RefVector<pat::PackedCandidateCollection> PackedCandidateRefVector;
} // namespace pat

#endif