diff --git a/include/Rivet/Jet.hh b/include/Rivet/Jet.hh
--- a/include/Rivet/Jet.hh
+++ b/include/Rivet/Jet.hh
@@ -1,271 +1,265 @@
// -*- C++ -*-
#ifndef RIVET_Jet_HH
#define RIVET_Jet_HH
#include "Rivet/Config/RivetCommon.hh"
#include "Rivet/Jet.fhh"
#include "Rivet/Particle.hh"
#include "Rivet/Tools/Cuts.hh"
#include "Rivet/Tools/Utils.hh"
#include "Rivet/Tools/RivetFastJet.hh"
#include "Rivet/Math/LorentzTrans.hh"
#include <numeric>
namespace Rivet {
/// @brief Representation of a clustered jet of particles.
class Jet : public ParticleBase {
public:
/// @name Constructors
//@{
/// Constructor from a FastJet PseudoJet, with optional full particle constituents information.
Jet(const fastjet::PseudoJet& pj, const Particles& particles=Particles(), const Particles& tags=Particles()) {
setState(pj, particles, tags);
}
/// Set the jet data, with optional full particle information.
Jet(const FourMomentum& pjet, const Particles& particles=Particles(), const Particles& tags=Particles()) {
setState(pjet, particles, tags);
}
/// Default constructor -- only for STL storability
Jet() { clear(); }
//@}
/// @name Access jet constituents
//@{
/// Number of particles in this jet.
size_t size() const { return _particles.size(); }
/// Get the particles in this jet.
Particles& particles() { return _particles; }
/// Get the particles in this jet (const version)
const Particles& particles() const { return _particles; }
/// Get the particles in this jet which pass a cut (const)
const Particles particles(const Cut& c) const { return filter_select(_particles, c); }
/// Get the particles in this jet which pass a filtering functor (const)
const Particles particles(const ParticleSelector& s) const { return filter_select(_particles, s); }
/// Get the particles in this jet (FastJet-like alias)
Particles& constituents() { return particles(); }
/// Get the particles in this jet (FastJet-like alias, const version)
const Particles& constituents() const { return particles(); }
/// Get the particles in this jet which pass a cut (FastJet-like alias, const)
const Particles constituents(const Cut& c) const { return particles(c); }
/// Get the particles in this jet which pass a filtering functor (FastJet-like alias, const)
const Particles constituents(const ParticleSelector& s) const { return particles(s); }
/// Check whether this jet contains a particular particle.
bool containsParticle(const Particle& particle) const;
/// Nicer alias for containsParticleId
bool containsPID(const Particle& particle) const { return containsParticle(particle); }
/// Check whether this jet contains a certain particle type.
bool containsParticleId(PdgId pid) const;
/// Nicer alias for containsParticleId
bool containsPID(PdgId pid) const { return containsParticleId(pid); }
/// Check whether this jet contains at least one of certain particle types.
bool containsParticleId(const vector<PdgId>& pids) const;
/// Nicer alias for containsParticleId
bool containsPID(const vector<PdgId>& pids) const { return containsParticleId(pids); }
//@}
/// @name Tagging
///
/// @note General sources of tag particles are planned. The default jet finding
/// adds b-hadron, c-hadron, and tau tags by ghost association.
//@{
/// @brief Particles which have been tag-matched to this jet
Particles& tags() { return _tags; }
/// @brief Particles which have been tag-matched to this jet (const version)
const Particles& tags() const { return _tags; }
/// @brief Particles which have been tag-matched to this jet _and_ pass a selector function
///
/// @note Note the less efficient return by value, due to the filtering.
Particles tags(const ParticleSelector& f) const { return filter_select(tags(), f); }
/// @brief Particles which have been tag-matched to this jet _and_ pass a Cut
///
/// @note Note the less efficient return by value, due to the cut-pass filtering.
Particles tags(const Cut& c) const;
/// @brief b particles which have been tag-matched to this jet (and pass an optional Cut)
///
/// The default jet finding adds b-hadron tags by ghost association.
Particles bTags(const Cut& c=Cuts::open()) const;
/// @brief b particles which have been tag-matched to this jet _and_ pass a selector function
Particles bTags(const ParticleSelector& f) const { return filter_select(bTags(), f); }
/// Does this jet have at least one b-tag (that passes an optional Cut)?
bool bTagged(const Cut& c=Cuts::open()) const { return !bTags(c).empty(); }
/// Does this jet have at least one b-tag (that passes the supplied selector function)?
bool bTagged(const ParticleSelector& f) const { return !bTags(f).empty(); }
/// @brief c (and not b) particles which have been tag-matched to this jet (and pass an optional Cut)
///
/// The default jet finding adds c-hadron tags by ghost association.
Particles cTags(const Cut& c=Cuts::open()) const;
/// @brief c (and not b) particles which have been tag-matched to this jet and pass a selector function
Particles cTags(const ParticleSelector& f) const { return filter_select(cTags(), f); }
/// Does this jet have at least one c-tag (that passes an optional Cut)?
bool cTagged(const Cut& c=Cuts::open()) const { return !cTags(c).empty(); }
/// Does this jet have at least one c-tag (that passes the supplied selector function)?
bool cTagged(const ParticleSelector& f) const { return !cTags(f).empty(); }
/// @brief Tau particles which have been tag-matched to this jet (and pass an optional Cut)
///
/// The default jet finding adds tau tags by ghost association.
Particles tauTags(const Cut& c=Cuts::open()) const;
/// @brief Tau particles which have been tag-matched to this jet and pass a selector function
Particles tauTags(const ParticleSelector& f) const { return filter_select(tauTags(), f); }
/// Does this jet have at least one tau-tag (that passes an optional Cut)?
bool tauTagged(const Cut& c=Cuts::open()) const { return !tauTags(c).empty(); }
/// Does this jet have at least one tau-tag (that passes the supplied selector function)?
bool tauTagged(const ParticleSelector& f) const { return !tauTags(f).empty(); }
/// @brief Check whether this jet contains a bottom-flavoured hadron.
///
/// @deprecated The bTags() or bTagged() function is probably what you want
/// for tagging. This one ignores the tags() list and draws conclusions
/// based directly on the jet constituents; the other gives a much better match
/// to typical experimental methods.
///
/// @note The decision is made by first trying to find a bottom-flavoured particle
/// in the particles list. Most likely this will fail unless bottom hadrons
/// are set stable. If @a include_decay_products is true (the default), a
/// fallback is attempted, using the post-hadronization ancestor history of
/// all constituents.
DEPRECATED("Prefer the bTags() or bTagged() function")
bool containsBottom(bool include_decay_products=true) const;
/// @brief Check whether this jet contains a charm-flavoured hadron.
///
/// @deprecated The cTags() or cTagged() function is probably what you want
/// for tagging. This one ignores the tags() list and draws conclusions
/// based directly on the jet constituents; the other gives a much better match
/// to typical experimental methods.
///
/// @note The decision is made by first trying to find a charm-flavoured particle
/// in the particles list. Most likely this will fail unless charmed hadrons
/// are set stable. If @a include_decay_products is true (the default), a
/// fallback is attempted, using the post-hadronization ancestor history of
/// all constituents.
DEPRECATED("Prefer the cTags() or cTagged() function")
bool containsCharm(bool include_decay_products=true) const;
//@}
/// @name Effective jet 4-vector properties
//@{
/// Get equivalent single momentum four-vector.
const FourMomentum& momentum() const { return _momentum; }
/// Apply an active Lorentz transform to this jet
/// @note The Rivet jet momentum, constituent particles, and tag particles will be modified.
/// @warning The FastJet cluster sequence and pseudojets will not be modified: don't use them after transformation!
Jet& transformBy(const LorentzTransform& lt);
/// Get the total energy of this jet.
double totalEnergy() const { return momentum().E(); }
/// Get the energy carried in this jet by neutral particles.
double neutralEnergy() const;
/// Get the energy carried in this jet by hadrons.
double hadronicEnergy() const;
//@}
/// @name Interaction with FastJet
//@{
/// Access the internal FastJet3 PseudoJet (as a const reference)
const fastjet::PseudoJet& pseudojet() const { return _pseudojet; }
/// Cast operator to FastJet3 PseudoJet (as a const reference)
operator const fastjet::PseudoJet& () const { return pseudojet(); }
//@}
/// @name Set the jet constituents and properties
//@{
/// @brief Set the jet data from a FastJet PseudoJet, with optional particle constituents and tags lists.
///
/// @note The particles() list will be extracted from PseudoJet constituents
/// by default, making use of an attached user info if one is found.
Jet& setState(const fastjet::PseudoJet& pj, const Particles& particles=Particles(), const Particles& tags=Particles());
/// Set all the jet data, with optional full particle constituent and tag information.
Jet& setState(const FourMomentum& mom, const Particles& particles, const Particles& tags=Particles());
/// @brief Set the particles collection with full particle information.
///
/// If set, this overrides particle info extracted from the PseudoJet
Jet& setParticles(const Particles& particles);
Jet& setConstituents(const Particles& particles) { return setParticles(particles); }
/// Reset this jet as empty.
Jet& clear();
//@}
private:
/// FJ3 PseudoJet member to unify PseudoJet and Jet
fastjet::PseudoJet _pseudojet;
/// Full constituent particle information. (Filled from PseudoJet if possible.)
/// @todo Make these mutable or similar? Add a flag to force a cache rebuild?
Particles _particles;
/// Particles used to tag this jet (can be anything, but c and b hadrons are the most common)
Particles _tags;
/// Effective jet 4-vector (just for caching)
mutable FourMomentum _momentum;
};
/// @name String representation and streaming support
//@{
- /// Represent a Jet as a string.
- std::string to_str(const Jet& j);
-
/// Allow a Jet to be passed to an ostream.
- inline std::ostream& operator<<(std::ostream& os, const Jet& j) {
- os << to_str(j);
- return os;
- }
+ std::ostream& operator << (std::ostream& os, const Jet& j);
//@}
}
#include "Rivet/Tools/JetUtils.hh"
#endif
diff --git a/include/Rivet/Particle.hh b/include/Rivet/Particle.hh
--- a/include/Rivet/Particle.hh
+++ b/include/Rivet/Particle.hh
@@ -1,616 +1,603 @@
// -*- 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)
{ }
/// 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)
: ParticleBase(),
_original(gp), _id(gp->pdg_id()),
_momentum(gp->momentum())
{
const GenVertex* 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();
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 {
return _original;
}
/// Cast operator for conversion to GenParticle*
operator const GenParticle* () 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);
}
/// 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;
/// 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);
+ std::ostream& operator << (std::ostream& os, const Particle& p);
/// Allow ParticlePair to be passed to an ostream.
- inline std::ostream& operator<<(std::ostream& os, const ParticlePair& pp) {
- os << to_str(pp);
- return os;
- }
+ std::ostream& operator << (std::ostream& os, const ParticlePair& pp);
//@}
}
#include "Rivet/Tools/ParticleUtils.hh"
#endif
diff --git a/src/Core/Jet.cc b/src/Core/Jet.cc
--- a/src/Core/Jet.cc
+++ b/src/Core/Jet.cc
@@ -1,192 +1,202 @@
#include "Rivet/Jet.hh"
#include "Rivet/Tools/Cuts.hh"
#include "Rivet/Tools/ParticleName.hh"
#include "Rivet/Tools/Logging.hh"
#include "Rivet/Tools/ParticleIdUtils.hh"
namespace Rivet {
Jet& Jet::clear() {
_momentum = FourMomentum();
_pseudojet.reset(0,0,0,0);
_particles.clear();
return *this;
}
Jet& Jet::setState(const FourMomentum& mom, const Particles& particles, const Particles& tags) {
clear();
_momentum = mom;
_pseudojet = fastjet::PseudoJet(mom.px(), mom.py(), mom.pz(), mom.E());
_particles = particles;
_tags = tags;
return *this;
}
Jet& Jet::setState(const fastjet::PseudoJet& pj, const Particles& particles, const Particles& tags) {
clear();
_pseudojet = pj;
_momentum = FourMomentum(pj.e(), pj.px(), pj.py(), pj.pz());
_particles = particles;
_tags = tags;
// if (_particles.empty()) {
// foreach (const fastjet::PseudoJet pjc, _pseudojet.constituents()) {
// // If there is no attached user info, we can't create a meaningful particle, so skip
// if (!pjc.has_user_info<RivetFJInfo>()) continue;
// const RivetFJInfo& fjinfo = pjc.user_info<RivetFJInfo>();
// // Don't add ghosts to the particles list
// if (fjinfo.isGhost) continue;
// // Otherwise construct a Particle from the PseudoJet, preferably from an associated GenParticle
// ?if (fjinfo.genParticle != NULL) {
// _particles.push_back(Particle(fjinfo.genParticle));
// } else {
// if (fjinfo.pid == 0) continue; // skip if there is a null PID entry in the FJ info
// const FourMomentum pjcmom(pjc.e(), pjc.px(), pjc.py(), pjc.pz());
// _particles.push_back(Particle(fjinfo.pid, pjcmom));
// }
// }
// }
return *this;
}
Jet& Jet::setParticles(const Particles& particles) {
_particles = particles;
return *this;
}
bool Jet::containsParticle(const Particle& particle) const {
const int barcode = particle.genParticle()->barcode();
foreach (const Particle& p, particles()) {
if (p.genParticle()->barcode() == barcode) return true;
}
return false;
}
bool Jet::containsParticleId(PdgId pid) const {
foreach (const Particle& p, particles()) {
if (p.pid() == pid) return true;
}
return false;
}
bool Jet::containsParticleId(const vector<PdgId>& pids) const {
foreach (const Particle& p, particles()) {
foreach (PdgId pid, pids) {
if (p.pid() == pid) return true;
}
}
return false;
}
/// @todo Jet::containsMatch(Matcher m) { ... if m(pid) return true; ... }
Jet& Jet::transformBy(const LorentzTransform& lt) {
_momentum = lt.transform(_momentum);
for (Particle& p : _particles) p.transformBy(lt);
for (Particle& t : _tags) t.transformBy(lt);
_pseudojet.reset(_momentum.px(), _momentum.py(), _momentum.pz(), _momentum.E()); //< lose ClusterSeq etc.
return *this;
}
double Jet::neutralEnergy() const {
double e_neutral = 0.0;
foreach (const Particle& p, particles()) {
const PdgId pid = p.pid();
if (PID::threeCharge(pid) == 0) {
e_neutral += p.E();
}
}
return e_neutral;
}
double Jet::hadronicEnergy() const {
double e_hadr = 0.0;
foreach (const Particle& p, particles()) {
const PdgId pid = p.pid();
if (PID::isHadron(pid)) {
e_hadr += p.E();
}
}
return e_hadr;
}
bool Jet::containsCharm(bool include_decay_products) const {
foreach (const Particle& p, particles()) {
const PdgId pid = p.pid();
if (abs(pid) == PID::CQUARK) return true;
if (PID::isHadron(pid) && PID::hasCharm(pid)) return true;
if (include_decay_products) {
const HepMC::GenVertex* gv = p.genParticle()->production_vertex();
if (gv) {
foreach (const GenParticle* pi, Rivet::particles(gv, HepMC::ancestors)) {
const PdgId pid2 = pi->pdg_id();
if (PID::isHadron(pid2) && PID::hasCharm(pid2)) return true;
}
}
}
}
return false;
}
bool Jet::containsBottom(bool include_decay_products) const {
foreach (const Particle& p, particles()) {
const PdgId pid = p.pid();
if (abs(pid) == PID::BQUARK) return true;
if (PID::isHadron(pid) && PID::hasBottom(pid)) return true;
if (include_decay_products) {
const HepMC::GenVertex* gv = p.genParticle()->production_vertex();
if (gv) {
foreach (const GenParticle* pi, Rivet::particles(gv, HepMC::ancestors)) {
const PdgId pid2 = pi->pdg_id();
if (PID::isHadron(pid2) && PID::hasBottom(pid2)) return true;
}
}
}
}
return false;
}
Particles Jet::tags(const Cut& c) const {
return filter_select(tags(), c);
}
Particles Jet::bTags(const Cut& c) const {
Particles rtn;
for (const Particle& tp : tags()) {
if (hasBottom(tp) && c->accept(tp)) rtn.push_back(tp);
}
return rtn;
}
Particles Jet::cTags(const Cut& c) const {
Particles rtn;
for (const Particle& tp : tags()) {
/// @todo Is making b and c tags exclusive the right thing to do?
if (hasCharm(tp) && !hasBottom(tp) && c->accept(tp)) rtn.push_back(tp);
}
return rtn;
}
Particles Jet::tauTags(const Cut& c) const {
Particles rtn;
for (const Particle& tp : tags()) {
if (isTau(tp) && c->accept(tp)) rtn.push_back(tp);
}
return rtn;
}
+ /// Allow a Jet to be passed to an ostream.
+ std::ostream& operator << (std::ostream& os, const Jet& j) {
+ os << "Jet<" << j.mom()/GeV << " GeV; Nparticles=" << j.size() << "; ";
+ os << "bTag=" << boolalpha << j.bTagged() << ", ";
+ os << "cTag=" << boolalpha << j.cTagged() << ", ";
+ os << "tauTag=" << boolalpha << j.tauTagged() << ">";
+ return os;
+ }
+
+
}
diff --git a/src/Core/Particle.cc b/src/Core/Particle.cc
--- a/src/Core/Particle.cc
+++ b/src/Core/Particle.cc
@@ -1,263 +1,257 @@
#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> Particle::ancestors(const Cut& c, bool physical_only) const {
vector<Particle> rtn;
/// @todo Remove this const mess crap when HepMC doesn't suck
GenVertexPtr gv = const_cast<GenVertexPtr>( 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;
rtn += p;
}
return rtn;
}
vector<Particle> Particle::parents(const Cut& c) const {
vector<Particle> rtn;
/// @todo Remove this const mess crap when HepMC doesn't suck
GenVertexPtr gv = const_cast<GenVertexPtr>( 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::parents); it != gv->particles_end(HepMC::parents); ++it) {
const Particle p(*it);
if (c != Cuts::OPEN && !c->accept(p)) continue;
rtn += p;
}
return rtn;
}
vector<Particle> Particle::children(const Cut& c) const {
vector<Particle> rtn;
if (isStable()) return rtn;
/// @todo Remove this const mess crap when HepMC doesn't suck
GenVertexPtr gv = const_cast<GenVertexPtr>( 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;
}
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> Particle::allDescendants(const Cut& c, bool remove_duplicates) const {
vector<Particle> rtn;
if (isStable()) return rtn;
/// @todo Remove this const mess crap when HepMC doesn't suck
GenVertexPtr gv = const_cast<GenVertexPtr>( 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> Particle::stableDescendants(const Cut& c) const {
vector<Particle> rtn;
if (isStable()) return rtn;
/// @todo Remove this const mess crap when HepMC doesn't suck
GenVertexPtr gv = const_cast<GenVertexPtr>( 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;
}
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<GenParticlePtr, GenParticlePtr> 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) {
+ /// Allow a Particle to be passed to an ostream.
+ std::ostream& operator << (std::ostream& os, 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();
+ os << "Particle<" << pname << " @ " << p.momentum()/GeV << " GeV>";
+ return os;
}
- 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();
+ /// Allow ParticlePair to be passed to an ostream.
+ std::ostream& operator << (std::ostream& os, const ParticlePair& pp) {
+ os << "[" << pp.first << ", " << pp.second << "]";
+ return os;
}
}