PackedGenParticle.h

Go to the documentation of this file.

#ifndef __AnalysisDataFormats_PackedGenParticle_h__
#define __AnalysisDataFormats_PackedGenParticle_h__

#include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
#include "DataFormats/HepMCCandidate/interface/GenStatusFlags.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" */

class testPackedGenParticle;

namespace pat {
  class PackedGenParticle : public reco::Candidate {
  public:
    friend class ::testPackedGenParticle;

    typedef reco::CandidateCollection daughters;
    typedef math::XYZTLorentzVector LorentzVector;
    typedef math::PtEtaPhiMLorentzVector PolarLorentzVector;
    typedef math::XYZPoint Point;
    typedef math::XYZVector Vector;

    typedef unsigned int index;

    PackedGenParticle()
    : packedPt_(0), packedY_(0), packedPhi_(0), packedM_(0), 
    p4_(nullptr), p4c_(nullptr), vertex_(0,0,0),  pdgId_(0), charge_(0) { }
    explicit PackedGenParticle( const reco::GenParticle & c)
    : p4_(new PolarLorentzVector(c.pt(), c.eta(), c.phi(), c.mass())), p4c_( new LorentzVector(*p4_)), vertex_(0,0,0), pdgId_(c.pdgId()), charge_(c.charge()), mother_(c.motherRef(0)),
      statusFlags_(c.statusFlags()) { pack(); }
    explicit PackedGenParticle( const reco::GenParticle & c, const edm::Ref<reco::GenParticleCollection> &  mother)
    : p4_(new PolarLorentzVector(c.pt(), c.eta(), c.phi(), c.mass())), p4c_(new LorentzVector(*p4_)), vertex_(0,0,0), pdgId_(c.pdgId()), charge_(c.charge()), mother_(mother),
      statusFlags_(c.statusFlags()) { pack(); }

    PackedGenParticle(const PackedGenParticle& iOther)
    : packedPt_(iOther.packedPt_), packedY_(iOther.packedY_), packedPhi_(iOther.packedPhi_), packedM_(iOther.packedM_),
      p4_(nullptr),p4c_(nullptr),
      vertex_(iOther.vertex_), dxy_(iOther.dxy_), dz_(iOther.dz_),dphi_(iOther.dphi_),
      pdgId_(iOther.pdgId_),charge_(iOther.charge_),mother_(iOther.mother_),
      statusFlags_(iOther.statusFlags_) {
      if(iOther.p4c_) {
        p4_.store( new PolarLorentzVector(*iOther.p4_) );
        p4c_.store( new LorentzVector(*iOther.p4c_) );
      }
    }

    PackedGenParticle(PackedGenParticle&& iOther)
    : packedPt_(iOther.packedPt_), packedY_(iOther.packedY_), packedPhi_(iOther.packedPhi_), packedM_(iOther.packedM_),
      p4_(nullptr),p4c_(nullptr),
      vertex_(std::move(iOther.vertex_)), dxy_(iOther.dxy_), dz_(iOther.dz_),dphi_(iOther.dphi_),
      pdgId_(iOther.pdgId_),charge_(iOther.charge_),mother_(std::move(iOther.mother_)),
      statusFlags_(iOther.statusFlags_) {
      if(iOther.p4c_) {
        p4_.store( p4_.exchange(nullptr) );
        p4c_.store( p4c_.exchange(nullptr) );
      }
    }

    PackedGenParticle& operator=(PackedGenParticle&& iOther) {
      if(this != &iOther) {
        packedPt_ = iOther.packedPt_;
        packedY_ = iOther.packedY_;
        packedPhi_ = iOther.packedPhi_;
        packedM_ = iOther.packedM_;
        if(p4c_) {
          delete p4_.exchange(iOther.p4_.exchange(nullptr));
          delete p4c_.exchange(iOther.p4c_.exchange(nullptr)) ;
        } else {
          delete p4_.exchange(nullptr);
          delete p4c_.exchange(nullptr);
        }
        vertex_=std::move(iOther.vertex_);
        dxy_ = iOther.dxy_;
        dz_ = iOther.dz_;
        dphi_ = iOther.dphi_;
        pdgId_ = iOther.pdgId_;
        charge_ = iOther.charge_;
        mother_ = std::move(iOther.mother_);
        statusFlags_ = iOther.statusFlags_;
      }
      return *this;
    }

    PackedGenParticle& operator=(PackedGenParticle const& iOther) {
      PackedGenParticle c(iOther);
      *this = std::move(c);
      return *this;
    }

    ~PackedGenParticle() 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::GenParticleRef motherRef() const { 
    if(mother_.isNonnull() && mother_.isAvailable()&& mother_->status()==1 ){ //if pointing to the pruned version of myself
      if(mother_->numberOfMothers() > 0) 
          return mother_->motherRef(0); // return my mother's (that is actually myself) mother
        else
          return edm::Ref<reco::GenParticleCollection>(); // return null ref
    } else { 
      return mother_; //the stored ref is really my mother, or null, return that
    }
    }
    const reco::GenParticleRef & lastPrunedRef() const { return mother_; }

    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 {
    return charge_; 
    }
    void setCharge( int charge) override {charge_=charge;}
    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(); return p4_.load()->M()*p4_.load()->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(); }
    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 { 
        unpack(); // changing px,py,pz changes also mapping between dxy,dz and x,y,z
        *p4_ = PolarLorentzVector(p4.Pt(), p4.Eta(), p4.Phi(), p4.M());
        pack();
    }
    void setP4( const PolarLorentzVector & p4 ) override { 
        unpack(); // changing px,py,pz changes also mapping between dxy,dz and x,y,z
        *p4_ = p4; 
        pack();
    }
    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 {
      unpack(); // 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());
      pack();
    }
    const Point & vertex() const override { return vertex_; }//{ if (fromPV_) return Point(0,0,0); else return Point(0,0,100); }
    double vx() const override  {  return vertex_.X(); }//{ return 0; }
    double vy() const override  {  return vertex_.Y(); }//{ return 0; }
    double vz() const override  {  return vertex_.Z(); }//{ if (fromPV_) return 0; else return 100; }
    void setVertex( const Point & vertex ) override { vertex_ = vertex;  }

    enum PVAssoc { NoPV=0, PVLoose=1, PVTight=2, PVUsedInFit=3 } ;


    virtual float dxy() const { unpack(); return dxy_; }
    virtual float dz()  const { unpack(); return dz_; }
    virtual float dxy(const Point &p) const ;
    virtual float dz(const Point &p)  const ;


    int pdgId() const override   { return pdgId_; }
    // set PDG identifier                                                                 
    void setPdgId( int pdgId ) override   { pdgId_ = pdgId; }
    int status() const override   { return 1; } /*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;

    PackedGenParticle * clone() const override  {
      return new PackedGenParticle( *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;
    bool isMuon() const override;
    bool isStandAloneMuon() const override;
    bool isGlobalMuon() const override;
    bool isTrackerMuon() const override;
    bool isCaloMuon() const override;
    bool isPhoton() const override;
    bool isConvertedPhoton() const override;
    bool isJet() const override;
    
    const reco::GenStatusFlags &statusFlags() const { return statusFlags_; }
    reco::GenStatusFlags &statusFlags() { return statusFlags_; }
    
    //basic set of gen status flags accessible directly here
    //the rest accessible through statusFlags()
    //(see GenStatusFlags.h for their meaning)
    
    //these are robust, generator-independent functions for categorizing
    //mainly final state particles, but also intermediate hadrons/taus
    
    //is particle prompt (not from hadron, muon, or tau decay) and final state
    bool isPromptFinalState() const { return status()==1 && statusFlags_.isPrompt(); }
        
    //this particle is a direct decay product of a prompt tau and is final state
    //(eg an electron or muon from a leptonic decay of a prompt tau)
    bool isDirectPromptTauDecayProductFinalState() const { return status()==1 && statusFlags_.isDirectPromptTauDecayProduct(); }
    
    //these are generator history-dependent functions for tagging particles
    //associated with the hard process
    //Currently implemented for Pythia 6 and Pythia 8 status codes and history   
    //and may not have 100% consistent meaning across all types of processes
    //Users are strongly encouraged to stick to the more robust flags above,
    //as well as the expanded set available in GenStatusFlags.h
    
    //this particle is the final state direct descendant of a hard process particle  
    bool fromHardProcessFinalState() const { return status()==1 && statusFlags_.fromHardProcess(); }
        
    //this particle is a direct decay product of a hardprocess tau and is final state
    //(eg an electron or muon from a leptonic decay of a tau from the hard process)
    bool isDirectHardProcessTauDecayProductFinalState() const { return status()==1 && statusFlags_.isDirectHardProcessTauDecayProduct(); }
        

  protected:
    uint16_t packedPt_, packedY_, packedPhi_, packedM_;
    void pack(bool unpackAfterwards=true) ;
    void unpack() const ;
 
    mutable std::atomic<PolarLorentzVector*> p4_;
    mutable std::atomic<LorentzVector*> p4c_;
    Point vertex_;
    float dxy_, dz_, dphi_;
    int pdgId_;
    int8_t charge_;
    reco::GenParticleRef mother_;
    //status flags
    reco::GenStatusFlags statusFlags_;

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

  };

  typedef std::vector<pat::PackedGenParticle> PackedGenParticleCollection;
  typedef edm::Ref<pat::PackedGenParticleCollection> PackedGenParticleRef;
  typedef edm::RefVector<pat::PackedGenParticleCollection> PackedGenParticleRefVector;
}

#endif