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/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 "DataFormats/Math/interface/PtEtaPhiMass.h" */

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() :
      p4_(0,0,0,0), p4c_(0,0,0,0), vertex_(0,0,0), dphi_(0), pdgId_(0),
      qualityFlags_(0), pvRefKey_(reco::VertexRef::invalidKey()),
      unpacked_(false), unpackedVtx_(true), unpackedTrk_(false), dxydxy_(0),
      dzdz_(0),dxydz_(0),dlambdadz_(0), dphidxy_(0),dptdpt_(0),detadeta_(0),
      dphidphi_(0), packedHits_(0), normalizedChi2_(0) { }

    explicit PackedCandidate( const reco::Candidate & c,
                              const reco::VertexRefProd &pvRefProd,
                              reco::VertexRef::key_type pvRefKey) :
      p4_(c.pt(), c.eta(), c.phi(), c.mass()), p4c_(p4_), vertex_(c.vertex()),
      dphi_(0), pdgId_(c.pdgId()), qualityFlags_(0), pvRefProd_(pvRefProd),
      pvRefKey_(pvRefKey), unpacked_(true), unpackedVtx_(true),
      unpackedTrk_(false), dxydxy_(0),dzdz_(0),dxydz_(0), dlambdadz_(0),
      dphidxy_(0),dptdpt_(0), detadeta_(0),dphidphi_(0), packedHits_(0),
      normalizedChi2_(0) {
      packBoth();
    }

    explicit PackedCandidate( const PolarLorentzVector &p4, const Point &vtx,
                              float phiAtVtx, int pdgId,
                              const reco::VertexRefProd &pvRefProd,
                              reco::VertexRef::key_type pvRefKey) :
      p4_(p4), p4c_(p4_), vertex_(vtx),
      dphi_(reco::deltaPhi(phiAtVtx,p4_.phi())), pdgId_(pdgId),
      qualityFlags_(0), pvRefProd_(pvRefProd), pvRefKey_(pvRefKey),
      unpacked_(true), unpackedVtx_(true), unpackedTrk_(false),
      dxydxy_(0),dzdz_(0),dxydz_(0),dlambdadz_(0),dphidxy_(0),dptdpt_(0),
      detadeta_(0),dphidphi_(0),packedHits_(0),normalizedChi2_(0) {
      packBoth();
    }

    explicit PackedCandidate( const LorentzVector &p4, const Point &vtx,
                              float phiAtVtx, int pdgId,
                              const reco::VertexRefProd &pvRefProd,
                              reco::VertexRef::key_type pvRefKey) :
      p4_(p4.Pt(), p4.Eta(), p4.Phi(), p4.M()), p4c_(p4), vertex_(vtx),
      dphi_(reco::deltaPhi(phiAtVtx,p4_.phi())), pdgId_(pdgId), qualityFlags_(0),
      pvRefProd_(pvRefProd), pvRefKey_(pvRefKey), unpacked_(true),
      unpackedVtx_(true), unpackedTrk_(false), dxydxy_(0),dzdz_(0),dxydz_(0),
      dlambdadz_(0),dphidxy_(0),dptdpt_(0), detadeta_(0),dphidphi_(0),
      packedHits_(0),normalizedChi2_(0) {
      packBoth();
    }

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

    virtual int charge() const {
      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
      }
    }
    virtual void setCharge( int charge) {}
    virtual int threeCharge() const {return charge()*3;}
    virtual void setThreeCharge( int threecharge) {}
    virtual const LorentzVector & p4() const { if (!unpacked_) unpack(); return p4c_; }  
    virtual const PolarLorentzVector & polarP4() const { if (!unpacked_) unpack(); return p4_; }
    virtual Vector momentum() const  { if (!unpacked_) unpack(); return p4c_.Vect(); }
    virtual Vector boostToCM() const { if (!unpacked_) unpack(); return p4c_.BoostToCM(); }
    virtual double p() const { if (!unpacked_) unpack(); return p4c_.P(); }
    virtual double energy() const { if (!unpacked_) unpack(); return p4c_.E(); }
    double et() const { return (pt()<=0) ? 0 : p4c_.Et(); }  
    double et2() const { return (pt()<=0) ? 0 : p4c_.Et2(); }   
    virtual double mass() const { if (!unpacked_) unpack(); return p4_.M(); }
    virtual double massSqr() const { if (!unpacked_) unpack(); return p4_.M()*p4_.M(); }

    virtual double mt() const { if (!unpacked_) unpack(); return p4_.Mt(); }
    virtual double mtSqr() const { if (!unpacked_) unpack(); return p4_.Mt2(); }
    virtual double px() const { if (!unpacked_) unpack(); return p4c_.Px(); }
    virtual double py() const { if (!unpacked_) unpack(); return p4c_.Py(); }
    virtual double pz() const { if (!unpacked_) unpack(); return p4c_.Pz(); }
    virtual double pt() const { if (!unpacked_) unpack(); return p4_.Pt();}
    virtual double phi() const { if (!unpacked_) unpack(); return p4_.Phi(); }
    virtual float phiAtVtx() const { 
        maybeUnpackBoth(); 
        float ret = p4_.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 double theta() const { if (!unpacked_) unpack(); return p4_.Theta(); }
    virtual double eta() const { if (!unpacked_) unpack(); return p4_.Eta(); }
    virtual double rapidity() const { if (!unpacked_) unpack(); return p4_.Rapidity(); }
    virtual double y() const { if (!unpacked_) unpack(); return p4_.Rapidity(); }
    virtual void setP4( const LorentzVector & p4 ) { 
        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();
    }
    virtual void setP4( const PolarLorentzVector & p4 ) { 
        maybeUnpackBoth(); // changing px,py,pz changes also mapping between dxy,dz and x,y,z
        p4_ = p4; 
        packBoth();
    }
    virtual void setMass( double m ) {
      if (!unpacked_) unpack(); 
      p4_ = PolarLorentzVector(p4_.Pt(), p4_.Eta(), p4_.Phi(), m); 
      pack();
    }
    virtual void setPz( double pz ) {
      maybeUnpackBoth(); // changing px,py,pz changes also mapping between dxy,dz and x,y,z
      p4c_ = LorentzVector(p4c_.Px(), p4c_.Py(), pz, p4c_.E());
      p4_  = PolarLorentzVector(p4c_.Pt(), p4c_.Eta(), p4c_.Phi(), p4c_.M());
      packBoth();
    }

    virtual void setTrackProperties( const reco::Track & tk, const reco::Track::CovarianceMatrix & covariance) {
      dxydxy_ = covariance(3,3);
      dxydz_ = covariance(3,4);
      dzdz_ = covariance(4,4);
      dphidxy_ = covariance(2,3);
      dlambdadz_ = covariance(1,4);
      dptdpt_ = covariance(0,0)*pt()*pt();
      detadeta_ = covariance(1,1);
      dphidphi_ = covariance(2,2)*pt()*pt();

      normalizedChi2_ = tk.normalizedChi2();
      int numberOfPixelHits_ = tk.hitPattern().numberOfValidPixelHits();
      if (numberOfPixelHits_ > 7) numberOfPixelHits_ = 7;
      int numberOfStripHits_ = tk.hitPattern().numberOfValidHits() - numberOfPixelHits_;
      if (numberOfStripHits_ > 31) numberOfStripHits_ = 31;
      packedHits_ = (numberOfPixelHits_&0x7) | (numberOfStripHits_ << 3);
      packBoth();
    }

    virtual void setTrackProperties( const reco::Track & tk ) {
    setTrackProperties(tk,tk.covariance());
    }   
 
    int numberOfPixelHits() const { return packedHits_ & 0x7; }
    int numberOfHits() const { return (packedHits_ >> 3) + numberOfPixelHits(); }
    
    virtual const Point & vertex() const { maybeUnpackBoth(); return vertex_; }//{ if (fromPV_) return Point(0,0,0); else return Point(0,0,100); }
    virtual double vx() const  { maybeUnpackBoth(); return vertex_.X(); }//{ return 0; }
    virtual double vy() const  { maybeUnpackBoth(); return vertex_.Y(); }//{ return 0; }
    virtual double vz() const  { maybeUnpackBoth(); return vertex_.Z(); }//{ if (fromPV_) return 0; else return 100; }
    virtual void setVertex( const Point & vertex ) { 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 ;

    virtual float dzError() const { maybeUnpackBoth(); return sqrt(dzdz_); }
    virtual float dxyError() const { maybeUnpackBoth(); return sqrt(dxydxy_); }


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

    virtual const reco::Track * bestTrack() const {
      if (packedHits_!=0) {
        if (!unpackedTrk_) unpackTrk();
        return &track_;
      }
      else
        return nullptr;
    }

    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 setMuonID(bool isStandAlone, bool isGlobal) {
        int16_t muonFlags = isStandAlone | (2*isGlobal);
        qualityFlags_ = (qualityFlags_ & ~muonFlagsMask) | ((muonFlags << muonFlagsShift) & muonFlagsMask);
    }

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

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

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

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

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

    virtual bool isElectron() const { return false; }
    virtual bool isMuon() const { return false; }
    virtual bool isStandAloneMuon() const { return ((qualityFlags_ & muonFlagsMask) >> muonFlagsShift) & 1; }
    virtual bool isGlobalMuon() const { return ((qualityFlags_ & muonFlagsMask) >> muonFlagsShift) & 2; }
    virtual bool isTrackerMuon() const { return false; }
    virtual bool isCaloMuon() const { return false; }
    virtual bool isPhoton() const { return false; }
    virtual bool isConvertedPhoton() const { return false; }
    virtual bool isJet() const { return false; }

    // puppiweights
    void setPuppiWeight(float p, float p_nolep = 0.0);  
    float puppiWeight() const;                          
    float puppiWeightNoLep() const;                     
    
  protected:
    uint16_t packedPt_, packedEta_, packedPhi_, packedM_;
    uint16_t packedDxy_, packedDz_, packedDPhi_;
    uint16_t packedCovarianceDxyDxy_,packedCovarianceDxyDz_,packedCovarianceDzDz_;
    int8_t packedCovarianceDlambdaDz_,packedCovarianceDphiDxy_;
    int8_t packedCovarianceDptDpt_,packedCovarianceDetaDeta_,packedCovarianceDphiDphi_;
    void pack(bool unpackAfterwards=true) ;
    void unpack() const ;
    void packVtx(bool unpackAfterwards=true) ;
    void unpackVtx() const ;
    void maybeUnpackBoth() const { if (!unpacked_) unpack(); if (!unpackedVtx_) unpackVtx(); }
    void packBoth() { pack(false); packVtx(false); 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
    mutable PolarLorentzVector p4_;
    mutable LorentzVector p4c_;
    mutable Point vertex_;
    mutable float dxy_, dz_, dphi_;
    mutable reco::Track track_;
    int pdgId_;
    uint16_t qualityFlags_;
    reco::VertexRefProd pvRefProd_;
    reco::VertexRef::key_type pvRefKey_;
    // is the momentum p4 unpacked
    mutable bool unpacked_;
    // are the dxy, dz and vertex unpacked
    mutable bool unpackedVtx_;
    // is the track unpacked
    mutable bool unpackedTrk_;
    mutable float dxydxy_, dzdz_, dxydz_,dlambdadz_,dphidxy_,dptdpt_,detadeta_,dphidphi_;
    uint8_t packedHits_;
    uint8_t normalizedChi2_; 
//    uint8_t numberOfPixelHits_;
  //  uint8_t numberOfHits_;
    
    virtual bool overlap( const reco::Candidate & ) const;
    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
    };
  };

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

#endif