Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:04:55

0001 #ifndef __AnalysisDataFormats_PackedGenParticle_h__
0002 #define __AnalysisDataFormats_PackedGenParticle_h__
0003 
0004 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
0005 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0006 #include "DataFormats/HepMCCandidate/interface/GenStatusFlags.h"
0007 #include "DataFormats/Candidate/interface/Candidate.h"
0008 #include "DataFormats/Candidate/interface/CandidateFwd.h"
0009 #include "DataFormats/Common/interface/RefVector.h"
0010 #include "DataFormats/Common/interface/Association.h"
0011 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0012 #include "DataFormats/VertexReco/interface/Vertex.h"
0013 #include "DataFormats/Math/interface/deltaPhi.h"
0014 /* #include "DataFormats/Math/interface/PtEtaPhiMass.h" */
0015 
0016 class testPackedGenParticle;
0017 
0018 namespace pat {
0019   class PackedGenParticle : public reco::Candidate {
0020   public:
0021     friend class ::testPackedGenParticle;
0022 
0023     /// collection of daughter candidates
0024     typedef reco::CandidateCollection daughters;
0025     /// Lorentz vector
0026     typedef math::XYZTLorentzVector LorentzVector;
0027     /// Lorentz vector
0028     typedef math::PtEtaPhiMLorentzVector PolarLorentzVector;
0029     /// point in the space
0030     typedef math::XYZPoint Point;
0031     /// point in the space
0032     typedef math::XYZVector Vector;
0033 
0034     typedef unsigned int index;
0035 
0036     /// default constructor
0037     PackedGenParticle()
0038         : packedPt_(0),
0039           packedY_(0),
0040           packedPhi_(0),
0041           packedM_(0),
0042           p4_(nullptr),
0043           p4c_(nullptr),
0044           vertex_(0, 0, 0),
0045           pdgId_(0),
0046           charge_(0) {}
0047     explicit PackedGenParticle(const reco::GenParticle& c)
0048         : p4_(new PolarLorentzVector(c.pt(), c.eta(), c.phi(), c.mass())),
0049           p4c_(new LorentzVector(*p4_)),
0050           vertex_(0, 0, 0),
0051           pdgId_(c.pdgId()),
0052           charge_(c.charge()),
0053           mother_(c.motherRef(0)),
0054           statusFlags_(c.statusFlags()) {
0055       pack();
0056     }
0057     explicit PackedGenParticle(const reco::GenParticle& c, const edm::Ref<reco::GenParticleCollection>& mother)
0058         : p4_(new PolarLorentzVector(c.pt(), c.eta(), c.phi(), c.mass())),
0059           p4c_(new LorentzVector(*p4_)),
0060           vertex_(0, 0, 0),
0061           pdgId_(c.pdgId()),
0062           charge_(c.charge()),
0063           mother_(mother),
0064           statusFlags_(c.statusFlags()) {
0065       pack();
0066     }
0067 
0068     PackedGenParticle(const PackedGenParticle& iOther)
0069         : packedPt_(iOther.packedPt_),
0070           packedY_(iOther.packedY_),
0071           packedPhi_(iOther.packedPhi_),
0072           packedM_(iOther.packedM_),
0073           p4_(nullptr),
0074           p4c_(nullptr),
0075           vertex_(iOther.vertex_),
0076           dxy_(iOther.dxy_),
0077           dz_(iOther.dz_),
0078           dphi_(iOther.dphi_),
0079           pdgId_(iOther.pdgId_),
0080           charge_(iOther.charge_),
0081           mother_(iOther.mother_),
0082           statusFlags_(iOther.statusFlags_) {
0083       if (iOther.p4c_) {
0084         p4_.store(new PolarLorentzVector(*iOther.p4_));
0085         p4c_.store(new LorentzVector(*iOther.p4c_));
0086       }
0087     }
0088 
0089     PackedGenParticle(PackedGenParticle&& iOther)
0090         : packedPt_(iOther.packedPt_),
0091           packedY_(iOther.packedY_),
0092           packedPhi_(iOther.packedPhi_),
0093           packedM_(iOther.packedM_),
0094           p4_(nullptr),
0095           p4c_(nullptr),
0096           vertex_(std::move(iOther.vertex_)),
0097           dxy_(iOther.dxy_),
0098           dz_(iOther.dz_),
0099           dphi_(iOther.dphi_),
0100           pdgId_(iOther.pdgId_),
0101           charge_(iOther.charge_),
0102           mother_(iOther.mother_),
0103           statusFlags_(iOther.statusFlags_) {
0104       if (iOther.p4c_) {
0105         p4_.store(p4_.exchange(nullptr));
0106         p4c_.store(p4c_.exchange(nullptr));
0107       }
0108     }
0109 
0110     PackedGenParticle& operator=(PackedGenParticle&& iOther) {
0111       if (this != &iOther) {
0112         packedPt_ = iOther.packedPt_;
0113         packedY_ = iOther.packedY_;
0114         packedPhi_ = iOther.packedPhi_;
0115         packedM_ = iOther.packedM_;
0116         if (p4c_) {
0117           delete p4_.exchange(iOther.p4_.exchange(nullptr));
0118           delete p4c_.exchange(iOther.p4c_.exchange(nullptr));
0119         } else {
0120           delete p4_.exchange(nullptr);
0121           delete p4c_.exchange(nullptr);
0122         }
0123         vertex_ = std::move(iOther.vertex_);
0124         dxy_ = iOther.dxy_;
0125         dz_ = iOther.dz_;
0126         dphi_ = iOther.dphi_;
0127         pdgId_ = iOther.pdgId_;
0128         charge_ = iOther.charge_;
0129         mother_ = iOther.mother_;
0130         statusFlags_ = iOther.statusFlags_;
0131       }
0132       return *this;
0133     }
0134 
0135     PackedGenParticle& operator=(PackedGenParticle const& iOther) {
0136       PackedGenParticle c(iOther);
0137       *this = std::move(c);
0138       return *this;
0139     }
0140 
0141     /// destructor
0142     ~PackedGenParticle() override;
0143     /// number of daughters
0144     size_t numberOfDaughters() const override;
0145     /// return daughter at a given position (throws an exception)
0146     const reco::Candidate* daughter(size_type) const override;
0147     /// number of mothers
0148     size_t numberOfMothers() const override;
0149     /// return mother at a given position (throws an exception)
0150     const reco::Candidate* mother(size_type) const override;
0151     /// direct access to the mother reference (may be null)
0152     reco::GenParticleRef motherRef() const {
0153       if (mother_.isNonnull() && mother_.isAvailable() &&
0154           mother_->status() == 1) {  //if pointing to the pruned version of myself
0155         if (mother_->numberOfMothers() > 0)
0156           return mother_->motherRef(0);  // return my mother's (that is actually myself) mother
0157         else
0158           return edm::Ref<reco::GenParticleCollection>();  // return null ref
0159       } else {
0160         return mother_;  //the stored ref is really my mother, or null, return that
0161       }
0162     }
0163     /// last surviving in pruned
0164     const reco::GenParticleRef& lastPrunedRef() const { return mother_; }
0165 
0166     /// return daughter at a given position (throws an exception)
0167     reco::Candidate* daughter(size_type) override;
0168     /// return daughter with a specified role name
0169     reco::Candidate* daughter(const std::string& s) override;
0170     /// return daughter with a specified role name
0171     const reco::Candidate* daughter(const std::string& s) const override;
0172     /// return the number of source Candidates
0173     /// ( the candidates used to construct this Candidate)
0174     size_t numberOfSourceCandidatePtrs() const override { return 0; }
0175     /// return a Ptr to one of the source Candidates
0176     /// ( the candidates used to construct this Candidate)
0177     reco::CandidatePtr sourceCandidatePtr(size_type i) const override { return reco::CandidatePtr(); }
0178 
0179     /// electric charge
0180     int charge() const override { return charge_; }
0181     /// set electric charge
0182     void setCharge(int charge) override { charge_ = charge; }
0183     /// electric charge
0184     int threeCharge() const override { return charge() * 3; }
0185     /// set electric charge
0186     void setThreeCharge(int threecharge) override {}
0187     /// four-momentum Lorentz vecto r
0188     const LorentzVector& p4() const override {
0189       if (!p4c_)
0190         unpack();
0191       return *p4c_;
0192     }
0193     /// four-momentum Lorentz vector
0194     const PolarLorentzVector& polarP4() const override {
0195       if (!p4c_)
0196         unpack();
0197       return *p4_;
0198     }
0199     /// spatial momentum vector
0200     Vector momentum() const override {
0201       if (!p4c_)
0202         unpack();
0203       return p4c_.load()->Vect();
0204     }
0205     /// boost vector to boost a Lorentz vector
0206     /// to the particle center of mass system
0207     Vector boostToCM() const override {
0208       if (!p4c_)
0209         unpack();
0210       return p4c_.load()->BoostToCM();
0211     }
0212     /// magnitude of momentum vector
0213     double p() const override {
0214       if (!p4c_)
0215         unpack();
0216       return p4c_.load()->P();
0217     }
0218     /// energy
0219     double energy() const override {
0220       if (!p4c_)
0221         unpack();
0222       return p4c_.load()->E();
0223     }
0224     /// transverse energy
0225     double et() const override { return (pt() <= 0) ? 0 : p4c_.load()->Et(); }
0226     /// transverse energy squared (use this for cuts)!
0227     double et2() const override { return (pt() <= 0) ? 0 : p4c_.load()->Et2(); }
0228     /// mass
0229     double mass() const override {
0230       if (!p4c_)
0231         unpack();
0232       return p4_.load()->M();
0233     }
0234     /// mass squared
0235     double massSqr() const override {
0236       if (!p4c_)
0237         unpack();
0238       return p4_.load()->M() * p4_.load()->M();
0239     }
0240 
0241     /// transverse mass
0242     double mt() const override {
0243       if (!p4c_)
0244         unpack();
0245       return p4_.load()->Mt();
0246     }
0247     /// transverse mass squared
0248     double mtSqr() const override {
0249       if (!p4c_)
0250         unpack();
0251       return p4_.load()->Mt2();
0252     }
0253     /// x coordinate of momentum vector
0254     double px() const override {
0255       if (!p4c_)
0256         unpack();
0257       return p4c_.load()->Px();
0258     }
0259     /// y coordinate of momentum vector
0260     double py() const override {
0261       if (!p4c_)
0262         unpack();
0263       return p4c_.load()->Py();
0264     }
0265     /// z coordinate of momentum vector
0266     double pz() const override {
0267       if (!p4c_)
0268         unpack();
0269       return p4c_.load()->Pz();
0270     }
0271     /// transverse momentum
0272     double pt() const override {
0273       if (!p4c_)
0274         unpack();
0275       return p4_.load()->Pt();
0276     }
0277     /// momentum azimuthal angle
0278     double phi() const override {
0279       if (!p4c_)
0280         unpack();
0281       return p4_.load()->Phi();
0282     }
0283     /// momentum polar angle
0284     double theta() const override {
0285       if (!p4c_)
0286         unpack();
0287       return p4_.load()->Theta();
0288     }
0289     /// momentum pseudorapidity
0290     double eta() const override {
0291       if (!p4c_)
0292         unpack();
0293       return p4_.load()->Eta();
0294     }
0295     /// rapidity
0296     double rapidity() const override {
0297       if (!p4c_)
0298         unpack();
0299       return p4_.load()->Rapidity();
0300     }
0301     /// rapidity
0302     double y() const override {
0303       if (!p4c_)
0304         unpack();
0305       return p4_.load()->Rapidity();
0306     }
0307     /// set 4-momentum
0308     void setP4(const LorentzVector& p4) override {
0309       unpack();  // changing px,py,pz changes also mapping between dxy,dz and x,y,z
0310       *p4_ = PolarLorentzVector(p4.Pt(), p4.Eta(), p4.Phi(), p4.M());
0311       pack();
0312     }
0313     /// set 4-momentum
0314     void setP4(const PolarLorentzVector& p4) override {
0315       unpack();  // changing px,py,pz changes also mapping between dxy,dz and x,y,z
0316       *p4_ = p4;
0317       pack();
0318     }
0319     /// set particle mass
0320     void setMass(double m) override {
0321       if (!p4c_)
0322         unpack();
0323       *p4_ = PolarLorentzVector(p4_.load()->Pt(), p4_.load()->Eta(), p4_.load()->Phi(), m);
0324       pack();
0325     }
0326     void setPz(double pz) override {
0327       unpack();  // changing px,py,pz changes also mapping between dxy,dz and x,y,z
0328       *p4c_ = LorentzVector(p4c_.load()->Px(), p4c_.load()->Py(), pz, p4c_.load()->E());
0329       *p4_ = PolarLorentzVector(p4c_.load()->Pt(), p4c_.load()->Eta(), p4c_.load()->Phi(), p4c_.load()->M());
0330       pack();
0331     }
0332     /// vertex position
0333     const Point& vertex() const override {
0334       return vertex_;
0335     }  //{ if (fromPV_) return Point(0,0,0); else return Point(0,0,100); }
0336     /// x coordinate of vertex position
0337     double vx() const override { return vertex_.X(); }  //{ return 0; }
0338     /// y coordinate of vertex position
0339     double vy() const override { return vertex_.Y(); }  //{ return 0; }
0340     /// z coordinate of vertex position
0341     double vz() const override { return vertex_.Z(); }  //{ if (fromPV_) return 0; else return 100; }
0342     /// set vertex
0343     void setVertex(const Point& vertex) override { vertex_ = vertex; }
0344 
0345     enum PVAssoc { NoPV = 0, PVLoose = 1, PVTight = 2, PVUsedInFit = 3 };
0346 
0347     /// dxy with respect to the PV ref
0348     virtual float dxy() const {
0349       unpack();
0350       return dxy_;
0351     }
0352     /// dz with respect to the PV ref
0353     virtual float dz() const {
0354       unpack();
0355       return dz_;
0356     }
0357     /// dxy with respect to another point
0358     virtual float dxy(const Point& p) const;
0359     /// dz  with respect to another point
0360     virtual float dz(const Point& p) const;
0361 
0362     /// PDG identifier
0363     int pdgId() const override { return pdgId_; }
0364     // set PDG identifier
0365     void setPdgId(int pdgId) override { pdgId_ = pdgId; }
0366     /// status word
0367     int status() const override { return 1; } /*FIXME*/
0368     /// set status word
0369     void setStatus(int status) override {} /*FIXME*/
0370     /// long lived flag
0371     static const unsigned int longLivedTag = 0; /*FIXME*/
0372     /// set long lived flag
0373     void setLongLived() override {} /*FIXME*/
0374     /// is long lived?
0375     bool longLived() const override;
0376     /// do mass constraint flag
0377     static const unsigned int massConstraintTag = 0; /*FIXME*/
0378     /// set mass constraint flag
0379     void setMassConstraint() override {} /*FIXME*/
0380     /// do mass constraint?
0381     bool massConstraint() const override;
0382 
0383     /// returns a clone of the Candidate object
0384     PackedGenParticle* clone() const override { return new PackedGenParticle(*this); }
0385 
0386     /// chi-squares
0387     double vertexChi2() const override;
0388     /** Number of degrees of freedom                                                                                   
0389      *  Meant to be Double32_t for soft-assignment fitters:                                                            
0390      *  tracks may contribute to the vertex with fractional weights.                                                   
0391      *  The ndof is then = to the sum of the track weights.                                                            
0392      *  see e.g. CMS NOTE-2006/032, CMS NOTE-2004/002                                                                  
0393      */
0394     double vertexNdof() const override;
0395     /// chi-squared divided by n.d.o.f.
0396     double vertexNormalizedChi2() const override;
0397     /// (i, j)-th element of error matrix, i, j = 0, ... 2
0398     double vertexCovariance(int i, int j) const override;
0399     /// return SMatrix
0400     CovarianceMatrix vertexCovariance() const override {
0401       CovarianceMatrix m;
0402       fillVertexCovariance(m);
0403       return m;
0404     }
0405     /// fill SMatrix
0406     void fillVertexCovariance(CovarianceMatrix& v) const override;
0407     /// returns true if this candidate has a reference to a master clone.
0408     /// This only happens if the concrete Candidate type is ShallowCloneCandidate
0409     bool hasMasterClone() const override;
0410     /// returns ptr to master clone, if existing.
0411     /// Throws an exception unless the concrete Candidate type is ShallowCloneCandidate
0412     const reco::CandidateBaseRef& masterClone() const override;
0413     /// returns true if this candidate has a ptr to a master clone.
0414     /// This only happens if the concrete Candidate type is ShallowClonePtrCandidate
0415     bool hasMasterClonePtr() const override;
0416     /// returns ptr to master clone, if existing.
0417     /// Throws an exception unless the concrete Candidate type is ShallowClonePtrCandidate
0418 
0419     const reco::CandidatePtr& masterClonePtr() const override;
0420 
0421     /// cast master clone reference to a concrete type
0422     template <typename Ref>
0423     Ref masterRef() const {
0424       return masterClone().template castTo<Ref>();
0425     }
0426     /// get a component
0427 
0428     bool isElectron() const override;
0429     bool isMuon() const override;
0430     bool isStandAloneMuon() const override;
0431     bool isGlobalMuon() const override;
0432     bool isTrackerMuon() const override;
0433     bool isCaloMuon() const override;
0434     bool isPhoton() const override;
0435     bool isConvertedPhoton() const override;
0436     bool isJet() const override;
0437 
0438     const reco::GenStatusFlags& statusFlags() const { return statusFlags_; }
0439     reco::GenStatusFlags& statusFlags() { return statusFlags_; }
0440 
0441     /////////////////////////////////////////////////////////////////////////////
0442     //basic set of gen status flags accessible directly here
0443     //the rest accessible through statusFlags()
0444     //(see GenStatusFlags.h for their meaning)
0445 
0446     /////////////////////////////////////////////////////////////////////////////
0447     //these are robust, generator-independent functions for categorizing
0448     //mainly final state particles, but also intermediate hadrons/taus
0449 
0450     //is particle prompt (not from hadron, muon, or tau decay) and final state
0451     bool isPromptFinalState() const { return status() == 1 && statusFlags_.isPrompt(); }
0452 
0453     //this particle is a direct decay product of a prompt tau and is final state
0454     //(eg an electron or muon from a leptonic decay of a prompt tau)
0455     bool isDirectPromptTauDecayProductFinalState() const {
0456       return status() == 1 && statusFlags_.isDirectPromptTauDecayProduct();
0457     }
0458 
0459     /////////////////////////////////////////////////////////////////////////////
0460     //these are generator history-dependent functions for tagging particles
0461     //associated with the hard process
0462     //Currently implemented for Pythia 6 and Pythia 8 status codes and history
0463     //and may not have 100% consistent meaning across all types of processes
0464     //Users are strongly encouraged to stick to the more robust flags above,
0465     //as well as the expanded set available in GenStatusFlags.h
0466 
0467     //this particle is the final state direct descendant of a hard process particle
0468     bool fromHardProcessFinalState() const { return status() == 1 && statusFlags_.fromHardProcess(); }
0469 
0470     //this particle is a direct decay product of a hardprocess tau and is final state
0471     //(eg an electron or muon from a leptonic decay of a tau from the hard process)
0472     bool isDirectHardProcessTauDecayProductFinalState() const {
0473       return status() == 1 && statusFlags_.isDirectHardProcessTauDecayProduct();
0474     }
0475 
0476   protected:
0477     uint16_t packedPt_, packedY_, packedPhi_, packedM_;
0478     void pack(bool unpackAfterwards = true);
0479     void unpack() const;
0480 
0481     /// the four vector
0482     mutable std::atomic<PolarLorentzVector*> p4_;
0483     mutable std::atomic<LorentzVector*> p4c_;
0484     /// vertex position
0485     Point vertex_;
0486     float dxy_, dz_, dphi_;
0487     /// PDG identifier
0488     int pdgId_;
0489     /// Charge
0490     int8_t charge_;
0491     ///Ref to first mother
0492     reco::GenParticleRef mother_;
0493     //status flags
0494     reco::GenStatusFlags statusFlags_;
0495 
0496     /// check overlap with another Candidate
0497     bool overlap(const reco::Candidate&) const override;
0498     template <typename, typename, typename>
0499     friend struct component;
0500     friend class ::OverlapChecker;
0501     friend class ShallowCloneCandidate;
0502     friend class ShallowClonePtrCandidate;
0503   };
0504 
0505   typedef std::vector<pat::PackedGenParticle> PackedGenParticleCollection;
0506   typedef edm::Ref<pat::PackedGenParticleCollection> PackedGenParticleRef;
0507   typedef edm::RefVector<pat::PackedGenParticleCollection> PackedGenParticleRefVector;
0508 }  // namespace pat
0509 
0510 #endif