diff --git a/include/Rivet/Particle.hh b/include/Rivet/Particle.hh --- a/include/Rivet/Particle.hh +++ b/include/Rivet/Particle.hh @@ -1,423 +1,420 @@ // -*- C++ -*- #ifndef RIVET_Particle_HH #define RIVET_Particle_HH #include "Rivet/Particle.fhh" #include "Rivet/ParticleBase.hh" #include "Rivet/Config/RivetCommon.hh" #include "Rivet/Tools/Cuts.hh" #include "Rivet/Tools/Utils.hh" #include "Rivet/Math/LorentzTrans.hh" // NOTE: Rivet/Tools/ParticleUtils.hh included at the end #include "fastjet/PseudoJet.hh" namespace Rivet { /// Particle representation, either from a HepMC::GenEvent or reconstructed. class Particle : public ParticleBase { public: /// @name Constructors //@{ /// Default constructor. /// @note A particle without info is useless. This only exists to keep STL containers happy. Particle() : ParticleBase(), _original(0), _id(0) { } /// Constructor without GenParticle. Particle(PdgId pid, const FourMomentum& mom, const FourVector& pos=FourVector()) : ParticleBase(), _original(0), _id(pid), _momentum(mom), _origin(pos) { } /// Constructor from a HepMC GenParticle pointer. - Particle(const GenParticle* gp) + Particle(const GenParticlePtr gp) : ParticleBase(), _original(gp), _id(gp->pdg_id()), _momentum(gp->momentum()) { - const GenVertex* vprod = gp->production_vertex(); + const GenVertexPtr vprod = gp->production_vertex(); if (vprod != NULL) setOrigin(vprod->position().t(), vprod->position().x(), vprod->position().y(), vprod->position().z()); } - /// Constructor from a HepMC GenParticle. - Particle(const GenParticle& gp) - : ParticleBase(), - _original(&gp), _id(gp.pdg_id()), - _momentum(gp.momentum()) - { - const GenVertex* vprod = gp.production_vertex(); - if (vprod != NULL) - setOrigin(vprod->position().t(), vprod->position().x(), vprod->position().y(), vprod->position().z()); - } + // /// Constructor from a HepMC GenParticle. + // Particle(const GenParticle& gp) + // : ParticleBase(), + // _original(&gp), _id(gp.pdg_id()), + // _momentum(gp.momentum()) + // { + // const GenVertexPtr vprod = gp.production_vertex(); + // if (vprod != NULL) + // setOrigin(vprod->position().t(), vprod->position().x(), vprod->position().y(), vprod->position().z()); + // } //@} /// @name Other representations and implicit casts //@{ /// Converter to FastJet3 PseudoJet virtual fastjet::PseudoJet pseudojet() const { return fastjet::PseudoJet(mom().px(), mom().py(), mom().pz(), mom().E()); } /// Cast operator to FastJet3 PseudoJet operator PseudoJet () const { return pseudojet(); } /// Get a const pointer to the original GenParticle. - const GenParticle* genParticle() const { + const GenParticlePtr genParticle() const { return _original; } - /// Cast operator for conversion to GenParticle* - operator const GenParticle* () const { return genParticle(); } + /// Cast operator for conversion to GenParticlePtr + operator const GenParticlePtr () const { return genParticle(); } //@} /// @name Kinematic properties //@{ /// The momentum. const FourMomentum& momentum() const { return _momentum; } /// Set the momentum. Particle& setMomentum(const FourMomentum& momentum) { _momentum = momentum; return *this; } /// Set the momentum via components. Particle& setMomentum(double E, double px, double py, double pz) { _momentum = FourMomentum(E, px, py, pz); return *this; } /// Apply an active Lorentz transform to this particle Particle& transformBy(const LorentzTransform& lt); //@ /// @name Positional properties //@{ /// The origin position. const FourVector& origin() const { return _origin; } /// Set the origin position. Particle& setOrigin(const FourVector& position) { _origin = position; return *this; } /// Set the origin position via components. Particle& setOrigin(double t, double x, double y, double z) { _origin = FourMomentum(t, x, y, z); return *this; } //@} /// @name Particle ID code accessors //@{ /// This Particle's PDG ID code. PdgId pid() const { return _id; } /// Absolute value of the PDG ID code. PdgId abspid() const { return std::abs(_id); } - /// This Particle's PDG ID code (alias). - /// @deprecated Prefer the pid/abspid form - PdgId pdgId() const { return _id; } //@} /// @name Particle species properties //@{ /// The charge of this Particle. double charge() const { return PID::charge(pid()); } /// The absolute charge of this Particle. double abscharge() const { return PID::abscharge(pid()); } /// Three times the charge of this Particle (i.e. integer multiple of smallest quark charge). int charge3() const { return PID::charge3(pid()); } /// Alias for charge3 /// @deprecated Use charge3 int threeCharge() const { return PID::threeCharge(pid()); } /// Three times the absolute charge of this Particle (i.e. integer multiple of smallest quark charge). int abscharge3() const { return PID::abscharge3(pid()); } /// Is this a hadron? bool isHadron() const { return PID::isHadron(pid()); } /// Is this a meson? bool isMeson() const { return PID::isMeson(pid()); } /// Is this a baryon? bool isBaryon() const { return PID::isBaryon(pid()); } /// Is this a lepton? bool isLepton() const { return PID::isLepton(pid()); } /// Is this a neutrino? bool isNeutrino() const { return PID::isNeutrino(pid()); } /// Does this (hadron) contain a b quark? bool hasBottom() const { return PID::hasBottom(pid()); } /// Does this (hadron) contain a c quark? bool hasCharm() const { return PID::hasCharm(pid()); } // /// Does this (hadron) contain an s quark? // bool hasStrange() const { return PID::hasStrange(pid()); } /// Is this particle potentially visible in a detector? bool isVisible() const; //@} /// @name Ancestry properties //@{ /// @todo Add physicalAncestors, allAncestors? /// Get a list of the direct parents of the current particle (with optional selection Cut) /// /// @note This is valid in MC, but may not be answerable /// experimentally -- use this function with care when replicating /// experimental analyses! Particles parents(const Cut& c=Cuts::OPEN) const; /// Get a list of the direct parents of the current particle (with selector function) /// /// @note This is valid in MC, but may not be answerable /// experimentally -- use this function with care when replicating /// experimental analyses! template Particles parents(const FN& f) const { return filter_select(parents(), f); } /// Check whether a given PID is found in the particle's parent list /// /// @note This question is valid in MC, but may not be answerable /// experimentally -- use this function with care when replicating /// experimental analyses! /// /// @deprecated Prefer e.g. parents(Cut::pid == 123).size() bool hasParent(PdgId pid) const; /// Check whether a particle in the particle's parent list has the requested property /// /// @note This question is valid in MC, but may not be answerable /// experimentally -- use this function with care when replicating /// experimental analyses! /// /// @deprecated Prefer parents(Cut) or parents(FN) methods and .empty() template bool hasParentWith(const FN& f) const { return _hasRelativeWith(HepMC::parents, f); } bool hasParentWith(const Cut& c) const; /// Check whether a given PID is found in the particle's ancestor list /// /// @note This question is valid in MC, but may not be answerable /// experimentally -- use this function with care when replicating /// experimental analyses! bool hasAncestor(PdgId pid) const; /// Check whether a particle in the particle's ancestor list has the requested property /// /// @note This question is valid in MC, but may not be answerable /// experimentally -- use this function with care when replicating /// experimental analyses! template bool hasAncestorWith(const FN& f) const { return _hasRelativeWith(HepMC::ancestors, f); } bool hasAncestorWith(const Cut& c) const; /// @brief Determine whether the particle is from a b-hadron decay /// /// @note This question is valid in MC, but may not be perfectly answerable /// experimentally -- use this function with care when replicating /// experimental analyses! bool fromBottom() const; /// @brief Determine whether the particle is from a c-hadron decay /// /// @note If a hadron contains b and c quarks it is considered a bottom /// hadron and NOT a charm hadron. /// /// @note This question is valid in MC, but may not be perfectly answerable /// experimentally -- use this function with care when replicating /// experimental analyses! bool fromCharm() const; // /// @brief Determine whether the particle is from a s-hadron decay // /// // /// @note If a hadron contains b or c quarks as well as strange it is // /// considered a b or c hadron, but NOT a strange hadron. // /// // /// @note This question is valid in MC, but may not be perfectly answerable // /// experimentally -- use this function with care when replicating // /// experimental analyses! // bool fromStrange() const; /// @brief Determine whether the particle is from a hadron decay /// /// @note This question is valid in MC, but may not be perfectly answerable /// experimentally -- use this function with care when replicating /// experimental analyses! bool fromHadron() const; /// @brief Determine whether the particle is from a tau decay /// /// @note This question is valid in MC, but may not be perfectly answerable /// experimentally -- use this function with care when replicating /// experimental analyses! bool fromTau(bool prompt_taus_only=false) const; /// @brief Determine whether the particle is from a prompt tau decay /// /// @note This question is valid in MC, but may not be perfectly answerable /// experimentally -- use this function with care when replicating /// experimental analyses! bool fromPromptTau() const { return fromTau(true); } /// @brief Determine whether the particle is from a hadron or tau decay /// /// Specifically, walk up the ancestor chain until a status 2 hadron or /// tau is found, if at all. /// /// @note This question is valid in MC, but may not be perfectly answerable /// experimentally -- use this function with care when replicating /// experimental analyses! bool fromDecay() const { return fromHadron() || fromPromptTau(); } /// @brief Shorthand definition of 'promptness' based on set definition flags /// /// The boolean arguments allow a decay lepton to be considered prompt if /// its parent was a "real" prompt lepton. /// /// @note This one doesn't make any judgements about final-stateness bool isPrompt(bool allow_from_prompt_tau=false, bool allow_from_prompt_mu=false) const; //@} /// @name Decay info //@{ /// Whether this particle is stable according to the generator bool isStable() const; /// Get a list of the direct descendants from the current particle (with optional selection Cut) Particles children(const Cut& c=Cuts::OPEN) const; /// Get a list of the direct descendants from the current particle (with selector function) template Particles children(const FN& f) const { return filter_select(children(), f); } /// Get a list of all the descendants from the current particle (with optional selection Cut) Particles allDescendants(const Cut& c=Cuts::OPEN, bool remove_duplicates=true) const; /// Get a list of all the descendants from the current particle (with selector function) template Particles allDescendants(const FN& f, bool remove_duplicates=true) const { return filter_select(allDescendants(Cuts::OPEN, remove_duplicates), f); } /// Get a list of all the stable descendants from the current particle (with optional selection Cut) /// /// @todo Use recursion through replica-avoiding MCUtils functions to avoid bookkeeping duplicates /// @todo Insist that the current particle is post-hadronization, otherwise throw an exception? Particles stableDescendants(const Cut& c=Cuts::OPEN) const; /// Get a list of all the stable descendants from the current particle (with selector function) template Particles stableDescendants(const FN& f) const { return filter_select(stableDescendants(), f); } /// Flight length (divide by mm or cm to get the appropriate units) double flightLength() const; //@} protected: template bool _hasRelativeWith(HepMC::IteratorRange relation, const FN& f) const { - for (const GenParticle* ancestor : particles(genParticle(), relation)) { + for (const GenParticlePtr ancestor : particles(genParticle(), relation)) { if (f(Particle(ancestor))) return true; } return false; } /// A pointer to the original GenParticle from which this Particle is projected. - const GenParticle* _original; + const GenParticlePtr _original; /// The PDG ID code for this Particle. PdgId _id; /// The momentum of this particle. FourMomentum _momentum; /// The creation position of this particle. FourVector _origin; }; /// @name String representation and streaming support //@{ /// Represent a Particle as a string. std::string to_str(const Particle& p); /// Allow a Particle to be passed to an ostream. inline std::ostream& operator<<(std::ostream& os, const Particle& p) { os << to_str(p); return os; } /// Represent a ParticlePair as a string. std::string to_str(const ParticlePair& pair); /// Allow ParticlePair to be passed to an ostream. inline std::ostream& operator<<(std::ostream& os, const ParticlePair& pp) { os << to_str(pp); return os; } //@} } #include "Rivet/Tools/ParticleUtils.hh" #endif diff --git a/include/Rivet/Tools/RivetHepMC.hh b/include/Rivet/Tools/RivetHepMC.hh --- a/include/Rivet/Tools/RivetHepMC.hh +++ b/include/Rivet/Tools/RivetHepMC.hh @@ -1,206 +1,308 @@ // -*- C++ -*- #ifndef RIVET_RivetHepMC_HH #define RIVET_RivetHepMC_HH #include "HepMC/GenEvent.h" #include "HepMC/GenParticle.h" #include "HepMC/GenVertex.h" -#include "HepMC/GenRanges.h" +#include "HepMC/Version.h" #include "HepMC/IO_GenEvent.h" #include "Rivet/Tools/RivetSTL.hh" #include "Rivet/Tools/Exceptions.hh" +#if HEPMC_VERSION_CODE >= 3000000 +#include "HepMC/HepMC.h" +#else +#include "HepMC/GenRanges.h" +#endif + namespace Rivet { + /// @todo Use mcutils? + using HepMC::GenEvent; using HepMC::GenParticle; using HepMC::GenVertex; #if HEPMC_VERSION_CODE >= 3000000 - using HepMC::GenEventPtr; using HepMC::GenParticlePtr; using HepMC::GenVertexPtr; #elif HEPMC_VERSION_CODE >= 2007000 // HepMC 2.07 provides its own #defines #else - #define GenEventPtr GenEvent* #define GenParticlePtr GenParticle* #define GenVertexPtr GenVertex* #endif - /// @todo Use mcutils? - inline std::vector particles(const GenEvent* ge) { + + #if HEPMC_VERSION_CODE >= 3000000 + + /// @name Accessors from GenEvent + //@{ + + inline std::vector particles(const GenEvent* ge) { assert(ge); - std::vector rtn; + return ge->particles(); + // std::vector rtn; + // for (const GenParticlePtr p : ge->particles()) rtn.push_back(p); + // return rtn; + } + + inline std::vector vertices(const GenEvent* ge) { + assert(ge); + return ge->vertices(); + // std::vector rtn; + // for (const GenVertexPtr v : ge->vertices()) rtn.push_back(v); + // return rtn; + } + + //@} + + + /// @name Accessors from GenVertex + //@{ + + inline const vector& particles_in(const GenVertexPtr& gv) { return gv->particles_in(); } + // inline vector& particles_in(GenVertexPtr& gv) { return gv->particles_in(); } + + inline const vector& particles_out(const GenVertexPtr& gv) { return gv->particles_out(); } + // inline vector& particles_out(GenVertexPtr& gv) { return gv->particles_out(); } + + inline std::vector particles(const GenVertexPtr& gv, HepMC::IteratorRange range) { + return HepMC::FindParticles(gv, range).results(); + } + + // /// Get the direct parents or all-ancestors of GenParticle @a gp + // inline std::vector particles_in(GenParticlePtr gp, HepMC::IteratorRange range=HepMC::ancestors) { + // assert(gp); + // if (range != HepMC::parents && range != HepMC::ancestors) + // throw UserError("Requested particles_in(GenParticlePtr) with a non-'in' iterator range"); + // return particles(gp->production_vertex(), range); + // } + // /// Get the direct children or all-descendents of GenParticle @a gp + // inline std::vector particles_out(GenParticlePtr gp, HepMC::IteratorRange range=HepMC::ancestors) { + // assert(gp); + // if (range != HepMC::children && range != HepMC::descendants) + // throw UserError("Requested particles_out(GenParticlePtr) with a non-'out' iterator range"); + // return particles(gp->production_vertex(), range); + // } + + //@} + + + /// @name Accessors from GenParticle + //@{ + + /// Get any relatives of GenParticle @a gp + inline std::vector particles(GenParticlePtr gp, HepMC::IteratorRange range) { + return HepMC::FindParticles(gp, range).results(); + } + //@} + + + #else + + + /// @name Accessors from GenEvent + //@{ + + inline std::vector particles(const GenEvent* ge) { + assert(ge); + std::vector rtn; + #if HEPMC_VERSION_CODE >= 2007000 + for (const GenParticlePtr p : ge->particles()) rtn.push_back(p); + #else for (GenEvent::particle_const_iterator pi = ge->particles_begin(); pi != ge->particles_end(); ++pi) rtn.push_back(*pi); + #endif return rtn; } - inline std::vector particles(GenEvent* ge) { assert(ge); - std::vector rtn; + std::vector rtn; + #if HEPMC_VERSION_CODE >= 2007000 + for (GenParticlePtr p : ge->particles()) rtn.push_back(p); + #else for (GenEvent::particle_iterator pi = ge->particles_begin(); pi != ge->particles_end(); ++pi) rtn.push_back(*pi); - return rtn; - } - - - inline std::vector vertices(const GenEvent* ge) { - std::vector rtn; - for (GenEvent::vertex_const_iterator vi = ge->vertices_begin(); vi != ge->vertices_end(); ++vi) - rtn.push_back(*vi); - return rtn; - } - - inline std::vector vertices(GenEvent* ge) { - std::vector rtn; - for (GenEvent::vertex_iterator vi = ge->vertices_begin(); vi != ge->vertices_end(); ++vi) - rtn.push_back(*vi); - return rtn; - } - - - ////////////////////////// - - - inline std::vector particles(const GenVertex* gv, HepMC::IteratorRange range=HepMC::relatives) { - std::vector rtn; - /// @todo A particle_const_iterator on GenVertex would be nice... - // Before HepMC 2.7.0 there were no GV::particles_const_iterators and constness consistency was all screwed up :-/ - #if HEPMC_VERSION_CODE >= 2007000 - // for (GenVertex::particle_iterator pi = gv->particles_const_begin(range); pi != gv->particles_const_end(range); ++pi) - for (GenVertex::particle_iterator pi = gv->particles_begin(range); pi != gv->particles_end(range); ++pi) - rtn.push_back(*pi); - #else - GenVertex* gv2 = const_cast(gv); - for (GenVertex::particle_iterator pi = gv2->particles_begin(range); pi != gv2->particles_end(range); ++pi) - rtn.push_back(const_cast(*pi)); #endif return rtn; } - inline std::vector particles(GenVertex* gv, HepMC::IteratorRange range=HepMC::relatives) { - std::vector rtn; + + inline std::vector vertices(const GenEvent* ge) { + assert(ge); + std::vector rtn; + #if HEPMC_VERSION_CODE >= 2007000 + for (const GenVertexPtr v : ge->vertices()) rtn.push_back(v); + #else + for (GenEvent::vertex_const_iterator vi = ge->vertices_begin(); vi != ge->vertices_end(); ++vi) + rtn.push_back(*vi); + #endif + return rtn; + } + inline std::vector vertices(GenEvent* ge) { + assert(ge); + std::vector rtn; + #if HEPMC_VERSION_CODE >= 2007000 + for (const GenVertexPtr v : ge->vertices()) rtn.push_back(v); + #else + for (GenEvent::vertex_iterator vi = ge->vertices_begin(); vi != ge->vertices_end(); ++vi) + rtn.push_back(*vi); + #endif + return rtn; + } + + //@} + + + /// @name Accessors from GenVertex + //@{ + + inline std::vector particles(const GenVertexPtr gv, HepMC::IteratorRange range=HepMC::relatives) { + std::vector rtn; + /// @todo A particle_const_iterator on GenVertex would be nice... + // Before HepMC 2.7.0 there were no GV::particles_const_iterators and constness consistency was all screwed up :-/ + #if HEPMC_VERSION_CODE >= 2007000 + for (const GenParticlePtr p : gv->particles(range)) rtn.push_back(p); + #else + GenVertexPtr gv2 = const_cast(gv); + for (GenVertex::particle_iterator pi = gv2->particles_begin(range); pi != gv2->particles_end(range); ++pi) + rtn.push_back(const_cast(*pi)); + #endif + return rtn; + } + inline std::vector particles(GenVertexPtr gv, HepMC::IteratorRange range=HepMC::relatives) { + std::vector rtn; for (GenVertex::particle_iterator pi = gv->particles_begin(range); pi != gv->particles_end(range); ++pi) rtn.push_back(*pi); return rtn; } - // Get iterator ranges as wrapped begin/end pairs /// @note GenVertex _in and _out iterators are actually, secretly the same types *sigh* struct GenVertexIterRangeC { - typedef vector::const_iterator genvertex_particles_const_iterator; + typedef vector::const_iterator genvertex_particles_const_iterator; GenVertexIterRangeC(const genvertex_particles_const_iterator& begin, const genvertex_particles_const_iterator& end) : _begin(begin), _end(end) { } const genvertex_particles_const_iterator& begin() { return _begin; } const genvertex_particles_const_iterator& end() { return _end; } private: const genvertex_particles_const_iterator _begin, _end; }; - inline GenVertexIterRangeC particles_in(const GenVertex* gv) { + inline GenVertexIterRangeC particles_in(const GenVertexPtr gv) { return GenVertexIterRangeC(gv->particles_in_const_begin(), gv->particles_in_const_end()); } - inline GenVertexIterRangeC particles_out(const GenVertex* gv) { + inline GenVertexIterRangeC particles_out(const GenVertexPtr gv) { return GenVertexIterRangeC(gv->particles_out_const_begin(), gv->particles_out_const_end()); } - #if HEPMC_VERSION_CODE >= 2007000 // Get iterator ranges as wrapped begin/end pairs /// @note GenVertex _in and _out iterators are actually, secretly the same types *sigh* struct GenVertexIterRange { - typedef vector::iterator genvertex_particles_iterator; + typedef vector::iterator genvertex_particles_iterator; GenVertexIterRange(const genvertex_particles_iterator& begin, const genvertex_particles_iterator& end) : _begin(begin), _end(end) { } const genvertex_particles_iterator& begin() { return _begin; } const genvertex_particles_iterator& end() { return _end; } private: const genvertex_particles_iterator _begin, _end; }; - inline GenVertexIterRange particles_in(GenVertex* gv) { + inline GenVertexIterRange particles_in(GenVertexPtr gv) { return GenVertexIterRange(gv->particles_in_begin(), gv->particles_in_end()); } - inline GenVertexIterRange particles_out(GenVertex* gv) { + inline GenVertexIterRange particles_out(GenVertexPtr gv) { return GenVertexIterRange(gv->particles_out_begin(), gv->particles_out_end()); } #endif - ////////////////////////// - + /// @name Accessors from GenParticle + //@{ /// Get the direct parents or all-ancestors of GenParticle @a gp - inline std::vector particles_in(const GenParticle* gp, HepMC::IteratorRange range=HepMC::ancestors) { + inline std::vector _particles_in(const GenParticlePtr gp, HepMC::IteratorRange range=HepMC::ancestors) { + assert(gp); if (range != HepMC::parents && range != HepMC::ancestors) - throw UserError("Requested particles_in(GenParticle*) with a non-'in' iterator range"); - if (!gp->production_vertex()) return std::vector(); + throw UserError("Requested particles_in(GenParticlePtr) with a non-'in' iterator range"); + if (!gp->production_vertex()) return std::vector(); #if HEPMC_VERSION_CODE >= 2007000 return particles(gp->production_vertex(), range); #else // Before HepMC 2.7.0 the constness consistency of methods and their return types was all screwed up :-/ - std::vector rtn; - foreach (GenParticle* gp2, particles(gp->production_vertex(), range)) - rtn.push_back( const_cast(gp2) ); + std::vector rtn; + for (GenParticlePtr gp2 : particles(gp->production_vertex(), range)) + rtn.push_back( const_cast(gp2) ); return rtn; #endif } - /// Get the direct parents or all-ancestors of GenParticle @a gp - inline std::vector particles_in(GenParticle* gp, HepMC::IteratorRange range=HepMC::ancestors) { + inline std::vector _particles_in(GenParticlePtr gp, HepMC::IteratorRange range=HepMC::ancestors) { + assert(gp); if (range != HepMC::parents && range != HepMC::ancestors) - throw UserError("Requested particles_in(GenParticle*) with a non-'in' iterator range"); - return (gp->production_vertex()) ? particles(gp->production_vertex(), range) : std::vector(); + throw UserError("Requested particles_in(GenParticlePtr) with a non-'in' iterator range"); + return (gp->production_vertex()) ? particles(gp->production_vertex(), range) : std::vector(); } /// Get the direct children or all-descendents of GenParticle @a gp - inline std::vector particles_out(const GenParticle* gp, HepMC::IteratorRange range=HepMC::descendants) { + inline std::vector _particles_out(const GenParticlePtr gp, HepMC::IteratorRange range=HepMC::descendants) { + assert(gp); if (range != HepMC::children && range != HepMC::descendants) - throw UserError("Requested particles_out(GenParticle*) with a non-'out' iterator range"); - if (!gp->end_vertex()) return std::vector(); + throw UserError("Requested particles_out(GenParticlePtr) with a non-'out' iterator range"); + if (!gp->end_vertex()) return std::vector(); #if HEPMC_VERSION_CODE >= 2007000 return particles(gp->end_vertex(), range); #else // Before HepMC 2.7.0 the constness consistency of methods and their return types was all screwed up :-/ - std::vector rtn; - foreach (GenParticle* gp2, particles(gp->end_vertex(), range)) - rtn.push_back( const_cast(gp2) ); + std::vector rtn; + foreach (GenParticlePtr gp2, particles(gp->end_vertex(), range)) + rtn.push_back( const_cast(gp2) ); return rtn; #endif } - /// Get the direct children or all-descendents of GenParticle @a gp - inline std::vector particles_out(GenParticle* gp, HepMC::IteratorRange range=HepMC::descendants) { + inline std::vector _particles_out(GenParticlePtr gp, HepMC::IteratorRange range=HepMC::descendants) { + assert(gp); if (range != HepMC::children && range != HepMC::descendants) - throw UserError("Requested particles_out(GenParticle*) with a non-'out' iterator range"); - return (gp->end_vertex()) ? particles(gp->end_vertex(), range) : std::vector(); + throw UserError("Requested particles_out(GenParticlePtr) with a non-'out' iterator range"); + return (gp->end_vertex()) ? particles(gp->end_vertex(), range) : std::vector(); } /// Get any relatives of GenParticle @a gp - inline std::vector particles(const GenParticle* gp, HepMC::IteratorRange range=HepMC::ancestors) { + inline std::vector particles(const GenParticlePtr gp, HepMC::IteratorRange range) { if (range == HepMC::parents || range == HepMC::ancestors) - return particles_in(gp, range); + return _particles_in(gp, range); if (range == HepMC::children || range == HepMC::descendants) - return particles_in(gp, range); - throw UserError("Requested particles(GenParticle*) with an unsupported iterator range"); + return _particles_out(gp, range); + throw UserError("Requested particles(const GenParticlePtr) with an unsupported iterator range"); } + /// Get any relatives of GenParticle @a gp + inline std::vector particles(GenParticlePtr gp, HepMC::IteratorRange range) { + if (range == HepMC::parents || range == HepMC::ancestors) + return _particles_in(gp, range); + if (range == HepMC::children || range == HepMC::descendants) + return _particles_out(gp, range); + throw UserError("Requested particles(GenParticlePtr) with an unsupported iterator range"); + } + + #endif } #endif