PackedCandidate.h

Go to the documentation of this file.

#ifndef __DataFormats_PatCandidates_PackedCandidate_h__
#define __DataFormats_PatCandidates_PackedCandidate_h__

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

/* #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),
      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), 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), 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), 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_),
      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_),
      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_;
      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_;
      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
        *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
        *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());
      *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);
    }

    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; }

    // 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 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 ;

    int8_t packedPuppiweight_;
    int8_t packedPuppiweightNoLepDiff_; // storing the DIFFERENCE of (all - "no lep") for compression optimization
    uint8_t rawCaloFraction_;
    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
    };
    
    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;
}

#endif