diff --git a/include/Rivet/Particle.hh b/include/Rivet/Particle.hh --- a/include/Rivet/Particle.hh +++ b/include/Rivet/Particle.hh @@ -1,617 +1,616 @@ // -*- 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(nullptr), _id(PID::ANY) + _original((GenParticle*) nullptr), _id(PID::ANY) { } /// Constructor without GenParticle. Particle(PdgId pid, const FourMomentum& mom, const FourVector& pos=FourVector()) : ParticleBase(), _original(nullptr), _id(pid), _momentum(mom), _origin(pos) { } /// Constructor from a HepMC GenParticle pointer. - Particle(const GenParticle* gp) + Particle(ConstGenParticlePtr 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 != nullptr) { 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(); + const GenVertexPtr vprod = gp.production_vertex(); if (vprod != nullptr) { 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 { + ConstGenParticlePtr genParticle() const { return _original; } /// Cast operator for conversion to GenParticle* - operator const GenParticle* () const { return genParticle(); } + operator ConstGenParticlePtr () 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 //@{ /// 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! Particles parents(const ParticleSelector& f) const { return filter_select(parents(), f); } /// Check whether any 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! bool hasParentWith(const ParticleSelector& f) const { return !parents(f).empty(); } /// Check whether any 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! bool hasParentWith(const Cut& c) const; /// Check whether any particle in the particle's parent list does not have 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! bool hasParentWithout(const ParticleSelector& f) const { return hasParentWith([&](const Particle& p){ return !f(p); }); } /// Check whether any particle in the particle's parent list does not have 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! bool hasParentWithout(const Cut& c) const; /// 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. hasParentWith(Cut::pid == 123) //DEPRECATED("Prefer e.g. hasParentWith(Cut::pid == 123)"); bool hasParent(PdgId pid) const; /// Get a list of the ancestors of the current particle (with optional selection Cut) /// /// @note By default only physical ancestors, with status=2, are returned. /// /// @note This is valid in MC, but may not be answerable experimentally -- /// use this function with care when replicating experimental analyses! Particles ancestors(const Cut& c=Cuts::OPEN, bool only_physical=true) const; /// Get a list of the direct parents of the current particle (with selector function) /// /// @note By default only physical ancestors, with status=2, are returned. /// /// @note This is valid in MC, but may not be answerable experimentally -- /// use this function with care when replicating experimental analyses! Particles ancestors(const ParticleSelector& f, bool only_physical=true) const { return filter_select(ancestors(Cuts::OPEN, only_physical), f); } - bool hasParentWith(const Cut& c) const; - + /// Check whether any 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! bool hasAncestorWith(const ParticleSelector& f, bool only_physical=true) const { return !ancestors(f, only_physical).empty(); } /// Check whether any 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! bool hasAncestorWith(const Cut& c, bool only_physical=true) const; /// Check whether any particle in the particle's ancestor list does not have 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! bool hasAncestorWithout(const ParticleSelector& f, bool only_physical=true) const { return hasAncestorWith([&](const Particle& p){ return !f(p); }, only_physical); } /// Check whether any particle in the particle's ancestor list does not have 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! bool hasAncestorWithout(const Cut& c, bool only_physical=true) 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! /// /// @deprecated Prefer hasAncestorWith(Cuts::pid == pid) etc. //DEPRECATED("Prefer e.g. hasAncestorWith(Cut::pid == 123)"); bool hasAncestor(PdgId pid, bool only_physical=true) 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 tau which decayed hadronically /// /// @note This question is valid in MC, but may not be perfectly answerable /// experimentally -- use this function with care when replicating /// experimental analyses! bool fromHadronicTau(bool prompt_taus_only=false) const; /// @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! /// /// @deprecated Too vague: use fromHadron or fromHadronicTau bool fromDecay() const { return fromHadron() || fromPromptTau(); } /// @brief Shorthand definition of 'promptness' based on set definition flags /// /// A "direct" particle is one directly connected to the hard process. It is a /// preferred alias for "prompt", since it has no confusing implications about /// distinguishability by timing information. /// /// The boolean arguments allow a decay lepton to be considered direct if /// its parent was a "real" direct lepton. /// /// @note This one doesn't make any judgements about final-stateness bool isDirect(bool allow_from_direct_tau=false, bool allow_from_direct_mu=false) const; /// Alias for isDirect bool isPrompt(bool allow_from_prompt_tau=false, bool allow_from_prompt_mu=false) const { return isDirect(allow_from_prompt_tau, allow_from_prompt_mu); } //@} /// @name Decay info //@{ /// Whether this particle is stable according to the generator bool isStable() const; /// @todo isDecayed? How to restrict to physical particles? /// 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) Particles children(const ParticleSelector& f) const { return filter_select(children(), f); } /// Check whether any direct child of this particle 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! bool hasChildWith(const ParticleSelector& f) const { return !children(f).empty(); } /// Check whether any direct child of this particle 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! bool hasChildWith(const Cut& c) const; /// Check whether any direct child of this particle does not have 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! bool hasChildWithout(const ParticleSelector& f) const { return hasChildWith([&](const Particle& p){ return !f(p); }); } /// Check whether any direct child of this particle does not have 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! bool hasChildWithout(const Cut& c) const; /// 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) Particles allDescendants(const ParticleSelector& f, bool remove_duplicates=true) const { return filter_select(allDescendants(Cuts::OPEN, remove_duplicates), f); } /// Check whether any descendant of this particle 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! bool hasDescendantWith(const ParticleSelector& f, bool remove_duplicates=true) const { return !allDescendants(f, remove_duplicates).empty(); } /// Check whether any descendant of this particle 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! bool hasDescendantWith(const Cut& c, bool remove_duplicates=true) const; /// Check whether any descendant of this particle does not have 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! bool hasDescendantWithout(const ParticleSelector& f, bool remove_duplicates=true) const { return hasDescendantWith([&](const Particle& p){ return !f(p); }, remove_duplicates); } /// Check whether any descendant of this particle does not have 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! bool hasDescendantWithout(const Cut& c, bool remove_duplicates=true) const; /// 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) Particles stableDescendants(const ParticleSelector& f) const { return filter_select(stableDescendants(), f); } /// Check whether any stable descendant of this particle 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! bool hasStableDescendantWith(const ParticleSelector& f) const { return !stableDescendants(f).empty(); } /// Check whether any stable descendant of this particle 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! bool hasStableDescendantWith(const Cut& c) const; /// Check whether any stable descendant of this particle does not have 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! bool hasStableDescendantWithout(const ParticleSelector& f) const { return hasStableDescendantWith([&](const Particle& p){ return !f(p); }); } /// Check whether any stable descendant of this particle does not have 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! bool hasStableDescendantWithout(const Cut& c) const; /// Flight length (divide by mm or cm to get the appropriate units) double flightLength() const; //@} /// @name Duplicate testing //@{ /// @brief Determine whether a particle is the first in a decay chain to meet the function requirement inline bool isFirstWith(const ParticleSelector& f) const { if (!f(*this)) return false; //< This doesn't even meet f, let alone being the last to do so if (any(parents(), f)) return false; //< If a direct parent has this property, this isn't the first return true; } /// @brief Determine whether a particle is the first in a decay chain not to meet the function requirement inline bool isFirstWithout(const ParticleSelector& f) const { return isFirstWith([&](const Particle& p){ return !f(p); }); } /// @brief Determine whether a particle is the last in a decay chain to meet the function requirement inline bool isLastWith(const ParticleSelector& f) const { if (!f(*this)) return false; //< This doesn't even meet f, let alone being the last to do so if (any(children(), f)) return false; //< If a child has this property, this isn't the last return true; } /// @brief Determine whether a particle is the last in a decay chain not to meet the function requirement inline bool isLastWithout(const ParticleSelector& f) const { return isLastWith([&](const Particle& p){ return !f(p); }); } //@} protected: /// A pointer to the original GenParticle from which this Particle is projected. - const GenParticle* _original; + ConstGenParticlePtr _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,308 +1,310 @@ // -*- C++ -*- #ifndef RIVET_RivetHepMC_HH #define RIVET_RivetHepMC_HH #include "HepMC/GenEvent.h" #include "HepMC/GenParticle.h" #include "HepMC/GenVertex.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::GenParticlePtr; + using HepMC::ConstGenParticlePtr; using HepMC::GenVertexPtr; #elif HEPMC_VERSION_CODE >= 2007000 // HepMC 2.07 provides its own #defines #else #define GenParticlePtr GenParticle* #define GenVertexPtr GenVertex* #endif #if HEPMC_VERSION_CODE >= 3000000 /// @name Accessors from GenEvent //@{ inline std::vector particles(const GenEvent* ge) { assert(ge); 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 //@{ + /// @todo are these really necessary? Why not call GenVertex::particles_[in, out] directly? 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; #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); #endif return 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; 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 GenVertexPtr gv) { return GenVertexIterRangeC(gv->particles_in_const_begin(), gv->particles_in_const_end()); } 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; 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(GenVertexPtr gv) { return GenVertexIterRange(gv->particles_in_begin(), gv->particles_in_end()); } 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 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"); 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; 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(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 (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 GenParticlePtr gp, HepMC::IteratorRange range=HepMC::descendants) { assert(gp); if (range != HepMC::children && range != HepMC::descendants) 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 (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(GenParticlePtr gp, HepMC::IteratorRange range=HepMC::descendants) { assert(gp); if (range != HepMC::children && range != HepMC::descendants) 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 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(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 diff --git a/src/Core/Particle.cc b/src/Core/Particle.cc --- a/src/Core/Particle.cc +++ b/src/Core/Particle.cc @@ -1,290 +1,294 @@ #include "Rivet/Particle.hh" #include "Rivet/Tools/Cuts.hh" #include "Rivet/Tools/ParticleIdUtils.hh" namespace Rivet { Particle& Particle::transformBy(const LorentzTransform& lt) { _momentum = lt.transform(_momentum); return *this; } bool Particle::isVisible() const { // Charged particles are visible if ( PID::threeCharge(pid()) != 0 ) return true; // Neutral hadrons are visible if ( PID::isHadron(pid()) ) return true; // Photons are visible if ( pid() == PID::PHOTON ) return true; // Gluons are visible (for parton level analyses) if ( pid() == PID::GLUON ) return true; // Everything else is invisible return false; } bool Particle::isStable() const { return genParticle() != NULL && genParticle()->status() == 1 && genParticle()->end_vertex() == NULL; } vector Particle::ancestors(const Cut& c, bool physical_only) const { + + /// @todo: use HepMC FindParticles when it actually works properly for const objects + //HepMC::FindParticles searcher(genParticle(), HepMC::ANCESTORS, HepMC::STATUS==1); + vector rtn; - /// @todo Remove this const mess crap when HepMC doesn't suck - GenVertexPtr gv = const_cast( genParticle()->production_vertex() ); + HepMC::ConstGenVertexPtr gv = genParticle()->production_vertex(); if (gv == NULL) return rtn; - /// @todo Would like to do this, but the range objects are broken - // foreach (const GenParticlePtr gp, gv->particles(HepMC::children)) - // rtn += Particle(gp); - for (GenVertex::particle_iterator it = gv->particles_begin(HepMC::ancestors); it != gv->particles_end(HepMC::ancestors); ++it) { - if (physical_only && (*it)->status() != 1 && (*it)->status() != 2) continue; - const Particle p(*it); - if (c != Cuts::OPEN && !c->accept(p)) continue; + + vector ancestors = genParticle()->ancestors(); + + for(const auto &a: ancestors){ + if(physical_only && a->status() != 1 && a->status() != 2) continue; + const Particle p(a); + if(c != Cuts::OPEN && !c->accept(p)) continue; rtn += p; } + return rtn; } vector Particle::parents(const Cut& c) const { vector rtn; #if HEPMC_VERSION_CODE >= 3000000 for (const GenParticlePtr gp : genParticle()->parents()) { const Particle p(gp); if (c != Cuts::OPEN && !c->accept(p)) continue; rtn += p; } #else GenVertexPtr gv = const_cast( genParticle()->production_vertex() ); if (gv == NULL) return rtn; for (GenVertex::particle_iterator it = gv->particles_begin(HepMC::parents); it != gv->particles_end(HepMC::parents); ++it) { const Particle p(*it); if (c != Cuts::OPEN && !c->accept(p)) continue; rtn += p; } #endif return rtn; } vector Particle::children(const Cut& c) const { vector rtn; if (isStable()) return rtn; #if HEPMC_VERSION_CODE >= 3000000 for (const GenParticlePtr gp : genParticle()->children()) { const Particle p(gp); if (c != Cuts::OPEN && !c->accept(p)) continue; rtn += p; } #else GenVertexPtr gv = const_cast( genParticle()->end_vertex() ); if (gv == NULL) return rtn; /// @todo Would like to do this, but the range objects are broken // foreach (const GenParticlePtr gp, gv->particles(HepMC::children)) // rtn += Particle(gp); for (GenVertex::particle_iterator it = gv->particles_begin(HepMC::children); it != gv->particles_end(HepMC::children); ++it) { const Particle p(*it); if (c != Cuts::OPEN && !c->accept(p)) continue; rtn += p; } #endif return rtn; } /// @todo Insist that the current particle is post-hadronization, otherwise throw an exception? /// @todo Use recursion through replica-avoiding functions to avoid bookkeeping duplicates vector Particle::allDescendants(const Cut& c, bool remove_duplicates) const { vector rtn; if (isStable()) return rtn; #if HEPMC_VERSION_CODE >= 3000000 for (const GenParticlePtr gp : genParticle()->descendants()) { const Particle p(gp); if (c != Cuts::OPEN && !c->accept(p)) continue; if (remove_duplicates) { bool dup = false; for (const GenParticlePtr gp2 : gp->children()) { if (gp->pid() == gp2->pid()) { dup = true; break; } } if (dup) continue; } rtn += p; } #else GenVertexPtr gv = const_cast( genParticle()->end_vertex() ); if (gv == NULL) return rtn; /// @todo Would like to do this, but the range objects are broken // foreach (const GenParticlePtr gp, gv->particles(HepMC::descendants)) for (GenVertex::particle_iterator it = gv->particles_begin(HepMC::descendants); it != gv->particles_end(HepMC::descendants); ++it) { const Particle p(*it); if (c != Cuts::OPEN && !c->accept(p)) continue; if (remove_duplicates && (*it)->end_vertex() != NULL) { // size_t n = 0; ///< @todo Only remove 1-to-1 duplicates? bool dup = false; /// @todo Yuck, HepMC for (GenVertex::particle_iterator it2 = (*it)->end_vertex()->particles_begin(HepMC::children); it2 != (*it)->end_vertex()->particles_end(HepMC::children); ++it2) { // n += 1; if (n > 1) break; if ((*it)->pdg_id() == (*it2)->pdg_id()) { dup = true; break; } } if (dup) continue; } rtn += p; } return rtn; } /// @todo Insist that the current particle is post-hadronization, otherwise throw an exception? vector Particle::stableDescendants(const Cut& c) const { vector rtn; if (isStable()) return rtn; /// @todo Remove this const mess crap when HepMC doesn't suck GenVertexPtr gv = const_cast( genParticle()->end_vertex() ); if (gv == NULL) return rtn; /// @todo Would like to do this, but the range objects are broken // foreach (const GenParticlePtr gp, gv->particles(HepMC::descendants)) for (GenVertex::particle_iterator it = gv->particles_begin(HepMC::descendants); it != gv->particles_end(HepMC::descendants); ++it) { // if ((*it)->status() != 1 || (*it)->end_vertex() != NULL) continue; const Particle p(*it); if (!p.isStable()) continue; if (c != Cuts::OPEN && !c->accept(p)) continue; rtn += p; } #endif return rtn; } double Particle::flightLength() const { if (isStable()) return -1; if (genParticle() == NULL) return 0; if (genParticle()->production_vertex() == NULL) return 0; const HepMC::FourVector v1 = genParticle()->production_vertex()->position(); const HepMC::FourVector v2 = genParticle()->end_vertex()->position(); return sqrt(sqr(v2.x()-v1.x()) + sqr(v2.y()-v1.y()) + sqr(v2.z()-v1.z())); } bool Particle::hasParent(PdgId pid) const { return hasParentWith(hasPID(pid)); } bool Particle::hasParentWith(const Cut& c) const { return hasParentWith([&](const Particle& p){return c->accept(p);}); } bool Particle::hasAncestor(PdgId pid, bool only_physical) const { return hasAncestorWith(hasPID(pid), only_physical); } bool Particle::hasAncestorWith(const Cut& c, bool only_physical) const { return hasAncestorWith([&](const Particle& p){return c->accept(p);}, only_physical); } bool Particle::hasChildWith(const Cut& c) const { return hasChildWith([&](const Particle& p){return c->accept(p);}); } bool Particle::hasDescendantWith(const Cut& c, bool remove_duplicates) const { return hasDescendantWith([&](const Particle& p){return c->accept(p);}, remove_duplicates); } bool Particle::hasStableDescendantWith(const Cut& c) const { return hasStableDescendantWith([&](const Particle& p){return c->accept(p);}); } bool Particle::fromBottom() const { return hasAncestorWith([](const Particle& p){ return p.genParticle()->status() == 2 && p.isHadron() && p.hasBottom(); }); } bool Particle::fromCharm() const { return hasAncestorWith([](const Particle& p){ return p.genParticle()->status() == 2 && p.isHadron() && p.hasCharm(); }); } bool Particle::fromHadron() const { return hasAncestorWith([](const Particle& p){ return p.genParticle()->status() == 2 && p.isHadron(); }); } bool Particle::fromTau(bool prompt_taus_only) const { if (prompt_taus_only && fromHadron()) return false; return hasAncestorWith([](const Particle& p){ return p.genParticle()->status() == 2 && isTau(p); }); } bool Particle::fromHadronicTau(bool prompt_taus_only) const { return hasAncestorWith([&](const Particle& p){ return p.genParticle()->status() == 2 && isTau(p) && (!prompt_taus_only || p.isPrompt()) && hasHadronicDecay(p); }); } bool Particle::isDirect(bool allow_from_direct_tau, bool allow_from_direct_mu) const { if (genParticle() == NULL) return false; // no HepMC connection, give up! Throw UserError exception? const GenVertexPtr prodVtx = genParticle()->production_vertex(); if (prodVtx == NULL) return false; // orphaned particle, has to be assume false const pair beams = prodVtx->parent_event()->beam_particles(); /// @todo Would be nicer to be able to write this recursively up the chain, exiting as soon as a parton or string/cluster is seen for (const GenParticlePtr ancestor : Rivet::particles(prodVtx, HepMC::ancestors)) { const PdgId pid = ancestor->pdg_id(); if (ancestor->status() != 2) continue; // no non-standard statuses or beams to be used in decision making if (ancestor == beams.first || ancestor == beams.second) continue; // PYTHIA6 uses status 2 for beams, I think... (sigh) if (PID::isParton(pid)) continue; // PYTHIA6 also uses status 2 for some partons, I think... (sigh) if (PID::isHadron(pid)) return false; // direct particles can't be from hadron decays if (abs(pid) == PID::TAU && abspid() != PID::TAU && !allow_from_direct_tau) return false; // allow or ban particles from tau decays (permitting tau copies) if (abs(pid) == PID::MUON && abspid() != PID::MUON && !allow_from_direct_mu) return false; // allow or ban particles from muon decays (permitting muon copies) } return true; } /////////////////////// string to_str(const Particle& p) { string pname; try { pname = PID::toParticleName(p.pid()); } catch (...) { pname = "PID=" + to_str(p.pid()); } stringstream out; out << pname << " @ " << p.momentum() << " GeV"; return out.str(); } string to_str(const ParticlePair& pair) { stringstream out; out << "[" << pair.first << ", " << pair.second << "]"; // out << "[" // << PID::toParticleName(pair.first.pid()) << " @ " // << pair.first.momentum().E()/GeV << " GeV, " // << PID::toParticleName(pair.second.pid()) << " @ " // << pair.second.momentum().E()/GeV << " GeV]"; return out.str(); } } diff --git a/src/Projections/Beam.cc b/src/Projections/Beam.cc --- a/src/Projections/Beam.cc +++ b/src/Projections/Beam.cc @@ -1,162 +1,166 @@ // -*- C++ -*- #include "Rivet/Config/RivetCommon.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/Beam.hh" namespace Rivet { ParticlePair beams(const Event& e) { #if HEPMC_VERSION_CODE >= 3000000 // First try the official way: ask the GenEvent for the beam pointers assert(e.genEvent()->particles().size() >= 2); const vector beams = e.genEvent()->beams(); if (beams.size() == 2 && beams[0] && beams[1]) { return ParticlePair{beams[0], beams[1]}; } // Ok, that failed: let's find the status = 4 particles by hand const vector pstat4s = e.allParticles([](const Particle& p){ return p.genParticle()->status() == 4; }); if (pstat4s.size() >= 2) { return ParticlePair{pstat4s[0], pstat4s[1]}; } + + /// There are no barcodes in HepMC3 + /// @todo implement some other fallback rubric? + /* // Hmm, this sucks. Last guess is that barcodes 1 and 2 are the beams if (e.genEvent()->barcode_to_particle(1) && e.genEvent()->barcode_to_particle(2)) { return ParticlePair{e.genEvent()->barcode_to_particle(1), e.genEvent()->barcode_to_particle(2)}; } - + */ #else // First try the official way: ask the GenEvent for the beam pointers assert(e.genEvent()->particles_size() >= 2); if (e.genEvent()->valid_beam_particles()) { pair beams = e.genEvent()->beam_particles(); assert(beams.first && beams.second); return ParticlePair{beams.first, beams.second}; } // Ok, that failed: let's find the status = 4 particles by hand const vector pstat4s = e.allParticles([](const Particle& p){ return p.genParticle()->status() == 4; }); if (pstat4s.size() >= 2) { return ParticlePair{pstat4s[0], pstat4s[1]}; } // Hmm, this sucks. Last guess is that barcodes 1 and 2 are the beams if (e.genEvent()->barcode_to_particle(1) && e.genEvent()->barcode_to_particle(2)) { return ParticlePair{e.genEvent()->barcode_to_particle(1), e.genEvent()->barcode_to_particle(2)}; } #endif // Give up: return null beams return ParticlePair{Particle(), Particle()}; } double sqrtS(const FourMomentum& pa, const FourMomentum& pb) { const double mom1 = pa.pz(); const double e1 = pa.E(); const double mom2 = pb.pz(); const double e2 = pb.E(); double sqrts = sqrt( sqr(e1+e2) - sqr(mom1+mom2) ); return sqrts; } double asqrtS(const FourMomentum& pa, const FourMomentum& pb) { const static double MNUCLEON = 939*MeV; //< nominal nucleon mass return sqrtS(pa/(pa.mass()/MNUCLEON), pb/(pb.mass()/MNUCLEON)); } double asqrtS(const ParticlePair& beams) { return sqrtS(beams.first.mom()/nuclA(beams.first), beams.second.mom()/nuclA(beams.second)); } FourMomentum acmsBoostVec(const FourMomentum& pa, const FourMomentum& pb) { const static double MNUCLEON = 939*MeV; //< nominal nucleon mass const double na = pa.mass()/MNUCLEON, nb = pb.mass()/MNUCLEON; return cmsBoostVec(pa/na, pb/nb); } FourMomentum acmsBoostVec(const ParticlePair& beams) { return cmsBoostVec(beams.first.mom()/nuclA(beams.first), beams.second.mom()/nuclA(beams.second)); } Vector3 cmsBetaVec(const FourMomentum& pa, const FourMomentum& pb) { // const Vector3 rtn = (pa.p3() + pb.p3()) / (pa.E() + pb.E()); const Vector3 rtn = (pa + pb).betaVec(); return rtn; } Vector3 acmsBetaVec(const FourMomentum& pa, const FourMomentum& pb) { const static double MNUCLEON = 939*MeV; //< nominal nucleon mass const Vector3 rtn = cmsBetaVec(pa/(pa.mass()/MNUCLEON), pb/(pb.mass()/MNUCLEON)); return rtn; } Vector3 acmsBetaVec(const ParticlePair& beams) { const Vector3 rtn = cmsBetaVec(beams.first.mom()/nuclA(beams.first), beams.second.mom()/nuclA(beams.second)); return rtn; } Vector3 cmsGammaVec(const FourMomentum& pa, const FourMomentum& pb) { // const Vector3 rtn = (pa + pb).gammaVec(); const double gamma = (pa.E() + pb.E()) / sqrt( sqr(pa.mass()) + sqr(pb.mass()) + 2*(pa.E()*pb.E() - dot(pa.p3(), pb.p3())) ); const Vector3 rtn = gamma * (pa.p3() + pb.p3()).unit(); return rtn; } Vector3 acmsGammaVec(const FourMomentum& pa, const FourMomentum& pb) { const static double MNUCLEON = 939*MeV; //< nominal nucleon mass Vector3 rtn = cmsGammaVec(pa/(pa.mass()/MNUCLEON), pb/(pb.mass()/MNUCLEON)); return rtn; } Vector3 acmsGammaVec(const ParticlePair& beams) { Vector3 rtn = cmsGammaVec(beams.first.mom()/nuclA(beams.first), beams.second.mom()/nuclA(beams.second)); return rtn; } LorentzTransform cmsTransform(const FourMomentum& pa, const FourMomentum& pb) { /// @todo Automatically choose to construct from beta or gamma according to which is more precise? return LorentzTransform::mkFrameTransformFromGamma(cmsGammaVec(pa, pb)); } LorentzTransform acmsTransform(const FourMomentum& pa, const FourMomentum& pb) { /// @todo Automatically choose to construct from beta or gamma according to which is more precise? return LorentzTransform::mkFrameTransformFromGamma(acmsGammaVec(pa, pb)); } LorentzTransform acmsTransform(const ParticlePair& beams) { return LorentzTransform::mkFrameTransformFromGamma(acmsGammaVec(beams)); } ///////////////////////////////////////////// void Beam::project(const Event& e) { _theBeams = Rivet::beams(e); MSG_DEBUG("Beam particles = " << _theBeams << " => sqrt(s) = " << sqrtS()/GeV << " GeV"); } FourVector Beam::pv() const { HepMC::FourVector v1, v2; const ParticlePair bpair = beams(); if (bpair.first.genParticle() && bpair.first.genParticle()->end_vertex()) v1 = bpair.first.genParticle()->end_vertex()->position(); if (bpair.second.genParticle() && bpair.second.genParticle()->end_vertex()) v2 = bpair.second.genParticle()->end_vertex()->position(); const FourVector rtn = (v1 == v2) ? FourVector(v1.t(), v1.x(), v1.y(), v1.z()) : FourVector(); MSG_DEBUG("Beam PV 4-position = " << rtn); return rtn; } } diff --git a/src/Projections/DISFinalState.cc b/src/Projections/DISFinalState.cc --- a/src/Projections/DISFinalState.cc +++ b/src/Projections/DISFinalState.cc @@ -1,32 +1,32 @@ // -*- C++ -*- #include "Rivet/Projections/DISFinalState.hh" namespace Rivet { void DISFinalState::project(const Event& e) { const DISKinematics& diskin = apply(e, "Kinematics"); LorentzTransform hcmboost; //< Null boost = LAB frame by default if (_boosttype == HCM) hcmboost = diskin.boostHCM(); else if (_boosttype == BREIT) hcmboost = diskin.boostBreit(); const DISLepton& dislep = diskin.apply(e, "Lepton"); const FinalState& fs = apply(e, "FS"); // Fill the particle list with all particles _other_ than the DIS scattered // lepton, with momenta boosted into the appropriate frame. _theParticles.clear(); _theParticles.reserve(fs.particles().size()-1); - const GenParticlePtr dislepGP = dislep.out().genParticle(); + ConstGenParticlePtr dislepGP = dislep.out().genParticle(); // const GenParticlePtr dislepIN = dislep.in().genParticle(); for (const Particle& p : fs.particles()) { ///< Ensure that we skip the DIS lepton Particle temp = p; if (_boosttype != LAB) temp.setMomentum(hcmboost.transform(temp.momentum())); if (p.genParticle() != dislepGP) _theParticles.push_back(temp); } } } diff --git a/src/Projections/DISLepton.cc b/src/Projections/DISLepton.cc --- a/src/Projections/DISLepton.cc +++ b/src/Projections/DISLepton.cc @@ -1,59 +1,60 @@ // -*- C++ -*- #include "Rivet/Projections/DISLepton.hh" namespace Rivet { int DISLepton::compare(const Projection& p) const { const DISLepton& other = pcast(p); return mkNamedPCmp(other, "Beam") || mkNamedPCmp(other, "FS"); } void DISLepton::project(const Event& e) { const ParticlePair& inc = applyProjection(e, "Beam").beams(); bool firstIsLepton = PID::isLepton(inc.first.pid()); bool secondIsLepton = PID::isLepton(inc.second.pid()); if (firstIsLepton && !secondIsLepton) { _incoming = inc.first; } else if (!firstIsLepton && secondIsLepton) { _incoming = inc.second; } else { //eek! throw Error("DISLepton projector could not find the correct beam."); } - const GenParticle* current_l=_incoming.genParticle(); + + ConstGenParticlePtr current_l=_incoming.genParticle(); bool found_next_vertex = true; while (found_next_vertex) { found_next_vertex = false; if (!current_l->end_vertex()) break; std::vector out_n; std::vector out_c; for (const GenParticle* pp : particles_out(current_l, HepMC::children)) { if (current_l->pdg_id() == pp->pdg_id()) out_n.push_back(pp); //+-1 should allow neutrino to electron and electron to neutrino if (std::abs(std::abs(current_l->pdg_id()) - std::abs(pp->pdg_id())) == 1) out_c.push_back(pp); } if (out_n.empty() && out_c.empty()) { MSG_WARNING("DISLepton projector: no electron/lepton in the new vertex."); break; } if (out_c.size() + out_n.size() > 1) { MSG_WARNING("DISLepton projector: more than one electron/lepton in the new vertex."); break; } if (out_c.size() == 1) current_l = out_c.front(); if (out_n.size() == 1) current_l = out_n.front(); found_next_vertex = true; } if (current_l == NULL) throw Error("DISLepton projector could not find the scattered lepton."); _outgoing = Particle(current_l); if (_outgoing.charge() == _incoming.charge()) _charged=0; else _charged = 0; // We consider only electric charge } }