diff --git a/include/Rivet/AnalysisHandler.hh b/include/Rivet/AnalysisHandler.hh
--- a/include/Rivet/AnalysisHandler.hh
+++ b/include/Rivet/AnalysisHandler.hh
@@ -1,242 +1,242 @@
// -*- C++ -*-
#ifndef RIVET_RivetHandler_HH
#define RIVET_RivetHandler_HH
#include "Rivet/Config/RivetCommon.hh"
#include "Rivet/Particle.hh"
#include "Rivet/AnalysisLoader.hh"
#include "Rivet/Tools/RivetYODA.hh"
namespace Rivet {
// Forward declaration and smart pointer for Analysis
class Analysis;
typedef std::shared_ptr<Analysis> AnaHandle;
// Needed to make smart pointers compare equivalent in the STL set
struct CmpAnaHandle {
bool operator() (const AnaHandle& a, const AnaHandle& b) {
return a.get() < b.get();
}
};
/// A class which handles a number of analysis objects to be applied to
/// generated events. An {@link Analysis}' AnalysisHandler is also responsible
/// for handling the final writing-out of histograms.
class AnalysisHandler {
public:
/// @name Constructors and destructors. */
//@{
/// Preferred constructor, with optional run name.
AnalysisHandler(const string& runname="");
/// @brief Destructor
/// The destructor is not virtual, as this class should not be inherited from.
~AnalysisHandler();
//@}
private:
/// Get a logger object.
Log& getLog() const;
public:
/// @name Run properties
//@{
/// Get the name of this run.
string runName() const;
/// Get the number of events seen. Should only really be used by external
/// steering code or analyses in the finalize phase.
size_t numEvents() const;
/// Get the sum of the event weights seen - the weighted equivalent of the
/// number of events. Should only really be used by external steering code
/// or analyses in the finalize phase.
double sumOfWeights() const;
/// Set sum of weights. This is useful if Rivet is steered externally and
/// the analyses are run for a sub-contribution of the events
/// (but of course have to be normalised to the total sum of weights)
void setSumOfWeights(const double& sum);
/// Is cross-section information required by at least one child analysis?
bool needCrossSection() const;
/// Set the cross-section for the process being generated.
AnalysisHandler& setCrossSection(double xs);
/// Get the cross-section known to the handler.
double crossSection() const {
return _xs;
}
/// Whether the handler knows about a cross-section.
bool hasCrossSection() const;
- /// Set beams for this run
+ /// Set the beam particles for this run
AnalysisHandler& setRunBeams(const ParticlePair& beams) {
_beams = beams;
MSG_DEBUG("Setting run beams = " << beams << " @ " << sqrtS()/GeV << " GeV");
return *this;
}
- /// Get beam IDs for this run, usually determined from the first event.
- const ParticlePair& beams() const {
- return _beams;
- }
+ /// Get the beam particles for this run, usually determined from the first event.
+ const ParticlePair& beams() const { return _beams; }
/// Get beam IDs for this run, usually determined from the first event.
+ /// @deprecated Use standalone beamIds(ah.beams()), to clean AH interface
PdgIdPair beamIds() const;
/// Get energy for this run, usually determined from the first event.
+ /// @deprecated Use standalone sqrtS(ah.beams()), to clean AH interface
double sqrtS() const;
/// Setter for _ignoreBeams
void setIgnoreBeams(bool ignore=true);
//@}
/// @name Handle analyses
//@{
/// Get a list of the currently registered analyses' names.
std::vector<std::string> analysisNames() const;
/// Get the collection of currently registered analyses.
const std::set<AnaHandle, CmpAnaHandle>& analyses() const {
return _analyses;
}
/// Get a registered analysis by name.
const AnaHandle analysis(const std::string& analysisname) const;
/// Add an analysis to the run list by object
AnalysisHandler& addAnalysis(Analysis* analysis);
/// @brief Add an analysis to the run list using its name.
///
/// The actual Analysis to be used will be obtained via
/// AnalysisLoader::getAnalysis(string). If no matching analysis is found,
/// no analysis is added (i.e. the null pointer is checked and discarded.
AnalysisHandler& addAnalysis(const std::string& analysisname);
/// @brief Add analyses to the run list using their names.
///
/// The actual {@link Analysis}' to be used will be obtained via
/// AnalysisHandler::addAnalysis(string), which in turn uses
/// AnalysisLoader::getAnalysis(string). If no matching analysis is found
/// for a given name, no analysis is added, but also no error is thrown.
AnalysisHandler& addAnalyses(const std::vector<std::string>& analysisnames);
/// Remove an analysis from the run list using its name.
AnalysisHandler& removeAnalysis(const std::string& analysisname);
/// Remove analyses from the run list using their names.
AnalysisHandler& removeAnalyses(const std::vector<std::string>& analysisnames);
//@}
/// @name Main init/execute/finalise
//@{
/// Initialize a run, with the run beams taken from the example event.
void init(const GenEvent& event);
/// @brief Analyze the given \a event by reference.
///
/// This function will call the AnalysisBase::analyze() function of all
/// included analysis objects.
void analyze(const GenEvent& event);
/// @brief Analyze the given \a event by pointer.
///
/// This function will call the AnalysisBase::analyze() function of all
/// included analysis objects, after checking the event pointer validity.
void analyze(const GenEvent* event);
/// Finalize a run. This function calls the AnalysisBase::finalize()
/// functions of all included analysis objects.
void finalize();
//@}
/// @name Histogram / data object access
//@{
/// Get all analyses' plots as a vector of analysis objects.
std::vector<AnalysisObjectPtr> getData() const;
/// Write all analyses' plots to the named file.
void writeData(const std::string& filename) const;
//@}
private:
/// The collection of Analysis objects to be used.
set<AnaHandle, CmpAnaHandle> _analyses;
/// @name Run properties
//@{
/// Run name
std::string _runname;
/// Number of events seen.
/// @todo Replace by a counter
unsigned int _numEvents;
/// Sum of event weights seen.
/// @todo Replace by a counter
double _sumOfWeights, _sumOfWeightsSq;
/// Cross-section known to AH
double _xs, _xserr;
/// Beams used by this run.
ParticlePair _beams;
/// Flag to check if init has been called
bool _initialised;
/// Flag whether input event beams should be ignored in compatibility check
bool _ignoreBeams;
//@}
private:
/// The assignment operator is private and must never be called.
/// In fact, it should not even be implemented.
AnalysisHandler& operator=(const AnalysisHandler&);
/// The copy constructor is private and must never be called. In
/// fact, it should not even be implemented.
AnalysisHandler(const AnalysisHandler&);
};
}
#endif
diff --git a/include/Rivet/Event.hh b/include/Rivet/Event.hh
--- a/include/Rivet/Event.hh
+++ b/include/Rivet/Event.hh
@@ -1,136 +1,156 @@
// -*- C++ -*-
#ifndef RIVET_Event_HH
#define RIVET_Event_HH
#include "Rivet/Config/RivetCommon.hh"
+#include "Rivet/Particle.hh"
#include "Rivet/Projection.hh"
+#include "Rivet/Math/Vector3.hh"
+#include "Rivet/Math/LorentzTrans.hh"
namespace Rivet {
/// Rivet wrapper for HepMC event and Projection references.
///
- /// Event is a concrete class representing an generated event in
- /// Rivet. It is constructed given a HepMC::GenEvent, a pointer to
- /// which is kept by the Event object throughout its lifetime. The user
- /// must therefore make sure that the corresponding HepMC::GenEvent
- /// will persist at least as long as the Event object.
+ /// Event is a concrete class representing an generated event in Rivet. It is
+ /// constructed given a HepMC::GenEvent, a pointer to which is kept by the
+ /// Event object throughout its lifetime. The user must therefore make sure
+ /// that the corresponding HepMC::GenEvent will persist at least as long as
+ /// the Event object.
///
- /// In addition to the HepMC::GenEvent object the Event also keeps
- /// track of all Projections object which have been applied to the
- /// Event so far.
+ /// In addition to the HepMC::GenEvent object the Event also keeps track of
+ /// all Projection objects which have been applied to the Event so far.
class Event {
public:
- /// @name Standard constructors and destructors.
+ /// @name Constructors and destructors.
//@{
- /// Constructor from a HepMC GenEvent reference
- Event(const GenEvent& ge)
- : _originalGenEvent(&ge), _genEvent(ge)
- { _init(ge); }
-
/// Constructor from a HepMC GenEvent pointer
Event(const GenEvent* ge)
: _originalGenEvent(ge), _genEvent(*ge)
{ assert(ge); _init(*ge); }
+ /// Constructor from a HepMC GenEvent reference
+ /// @deprecated HepMC uses pointers, so we should talk to HepMC via pointers
+ Event(const GenEvent& ge)
+ : _originalGenEvent(&ge), _genEvent(ge)
+ { _init(ge); }
+
/// Copy constructor
Event(const Event& e)
: _originalGenEvent(e._originalGenEvent), _genEvent(e._genEvent)
{ }
//@}
- public:
+ /// @name Major event properties
+ //@{
+
+ // /// Get the beam particles
+ // ParticlePair beams() const;
+
+ // /// Get the beam centre-of-mass energy
+ // double sqrtS() const;
+
+ // /// Get the beam centre-of-mass energy per nucleon
+ // double asqrtS() const;
+
+ // /// Get the boost to the beam centre-of-mass
+ // Vector3 beamCMSBoost() const;
+
+ // /// Get the boost to the beam centre-of-mass
+ // LorentzTransform beamCMSTransform();
/// The generated event obtained from an external event generator
- const GenEvent* genEvent() const {
- return &_genEvent;
- }
+ const GenEvent* genEvent() const { return &_genEvent; }
/// @brief The generation weight associated with the event
///
/// @todo This needs to be revisited when we finally add the mechanism to
/// support NLO counter-events and weight vectors.
double weight() const {
return (!_genEvent.weights().empty()) ? _genEvent.weights()[0] : 1.0;
}
+ //@}
- public:
+
+ /// @name Projection running
+ //@{
/// @brief Add a projection @a p to this Event.
///
/// If an equivalent Projection has been applied before, the
/// Projection::project(const Event&) of @a p is not called and a reference
/// to the previous equivalent projection is returned. If no previous
/// Projection was found, the Projection::project(const Event&) of @a p is
/// called and a reference to @a p is returned.
template <typename PROJ>
const PROJ& applyProjection(PROJ& p) const {
const Projection* cpp(&p);
std::set<const Projection*>::const_iterator old = _projections.find(cpp);
if (old != _projections.end()) {
const Projection& pRef = **old;
return pcast<PROJ>(pRef);
}
// Add the projection via the Projection base class (only
// possible because Event is a friend of Projection)
Projection* pp = const_cast<Projection*>(cpp);
pp->project(*this);
_projections.insert(pp);
return p;
}
-
+ /// @brief Add a projection @a p to this Event by pointer.
template <typename PROJ>
const PROJ& applyProjection(PROJ* pp) const {
if (!pp) throw Error("Event::applyProjection(PROJ*): Projection pointer is null.");
return applyProjection(*pp);
}
+ //@}
+
private:
/// @brief Actual (shared) implementation of the constructors from GenEvents
- ///
- /// @todo When we can require C++11, constructors can call each other and
- /// this can be removed.
void _init(const GenEvent& ge);
- /// @brief Convert the GenEvent to use conventional alignment
- ///
- /// For example, FHerwig only produces DIS events in the unconventional
- /// hadron-lepton orientation and has to be corrected for DIS analysis
- /// portability.
- void _geNormAlignment();
+ // /// @brief Convert the GenEvent to use conventional alignment
+ // ///
+ // /// For example, FHerwig only produces DIS events in the unconventional
+ // /// hadron-lepton orientation and has to be corrected for DIS analysis
+ // /// portability.
+ // void _geNormAlignment();
/// @brief The generated event, as obtained from an external generator.
///
/// This is the original GenEvent. In practise the version seen by users
/// will often/always be a modified one.
///
/// @todo Provide access to this via an Event::originalGenEvent() method? If requested...
const GenEvent* _originalGenEvent;
/// @brief The GenEvent used by Rivet analysis projections etc.
///
/// This version may be rotated to a "normal" alignment, have
/// generator-specific particles stripped out, etc. If an analysis is
/// affected by these modifications, it is probably an unphysical analysis!
///
/// Stored as a non-pointer since it may get overwritten, and memory for
- /// copying and cleanup is much neater this way.
+ /// copying and cleanup is neater this way.
+ /// @todo Change needed for HepMC3?
GenEvent _genEvent;
/// The set of Projection objects applied so far.
mutable std::set<ConstProjectionPtr> _projections;
};
}
#endif
diff --git a/include/Rivet/Particle.hh b/include/Rivet/Particle.hh
--- a/include/Rivet/Particle.hh
+++ b/include/Rivet/Particle.hh
@@ -1,384 +1,384 @@
// -*- 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/Utils.hh"
// NOTE: Rivet/Tools/ParticleUtils.hh included at the end
#include "fastjet/PseudoJet.hh"
namespace Rivet {
- /// Representation of particles from a HepMC::GenEvent.
+ /// Particle representation, either from a HepMC::GenEvent or reconstructed.
class Particle : public ParticleBase {
public:
/// Default constructor.
/// @note A particle without info is useless. This only exists to keep STL containers happy.
Particle()
: ParticleBase(),
_original(0), _id(0)
{ }
/// Constructor without GenParticle.
Particle(PdgId pid, const FourMomentum& mom, const FourVector& pos=FourVector())
: ParticleBase(),
_original(0), _id(pid), _momentum(mom), _origin(pos)
{ }
/// Constructor from a HepMC GenParticle.
Particle(const GenParticle& gp)
: ParticleBase(),
_original(&gp), _id(gp.pdg_id()),
_momentum(gp.momentum())
{
const GenVertex* vprod = gp.production_vertex();
if (vprod != NULL)
setOrigin(vprod->position().t(), vprod->position().x(), vprod->position().y(), vprod->position().z());
}
/// Constructor from a HepMC GenParticle pointer.
Particle(const GenParticle* gp)
: ParticleBase(),
_original(gp), _id(gp->pdg_id()),
_momentum(gp->momentum())
{
const GenVertex* vprod = gp->production_vertex();
if (vprod != NULL)
setOrigin(vprod->position().t(), vprod->position().x(), vprod->position().y(), vprod->position().z());
}
public:
/// 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(); }
public:
/// @name Basic particle specific properties
//@{
/// 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(); }
/// 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;
}
/// 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()); }
/// Three times the charge of this Particle (i.e. integer multiple of smallest quark charge).
int charge3() const { return PID::charge3(pid()); }
/// The abs charge of this Particle.
double abscharge() const { return PID::abscharge(pid()); }
/// Three times the abs charge of this Particle (i.e. integer multiple of smallest quark charge).
int abscharge3() const { return PID::abscharge3(pid()); }
/// Alias for charge3
/// @deprecated Use charge3
int threeCharge() const { return charge3(); }
/// 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()); }
/// Is this particle potentially visible in a detector?
bool 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;
}
// /// Does this (hadron) contain an s quark?
// bool hasStrange() const { return PID::hasStrange(pid()); }
//@}
/// @name Ancestry properties
//@{
/// Check whether a given PID is found in the GenParticle's ancestor list
///
/// @note This question is valid in MC, but may not be answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
bool hasAncestor(PdgId pdg_id) const;
/// @brief Determine whether the particle is from a b-hadron decay
///
/// @note This question is valid in MC, but may not be perfectly answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
bool fromBottom() const;
/// @brief Determine whether the particle is from a c-hadron decay
///
/// @note If a hadron contains b and c quarks it is considered a bottom
/// hadron and NOT a charm hadron.
///
/// @note This question is valid in MC, but may not be perfectly answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
bool fromCharm() const;
// /// @brief Determine whether the particle is from a s-hadron decay
// ///
// /// @note If a hadron contains b or c quarks as well as strange it is
// /// considered a b or c hadron, but NOT a strange hadron.
// ///
// /// @note This question is valid in MC, but may not be perfectly answerable
// /// experimentally -- use this function with care when replicating
// /// experimental analyses!
// bool fromStrange() const;
/// @brief Determine whether the particle is from a hadron decay
///
/// @note This question is valid in MC, but may not be perfectly answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
bool fromHadron() const;
/// @brief Determine whether the particle is from a tau decay
///
/// @note This question is valid in MC, but may not be perfectly answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
bool fromTau(bool prompt_taus_only=false) const;
/// @brief Determine whether the particle is from a prompt tau decay
///
/// @note This question is valid in MC, but may not be perfectly answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
bool fromPromptTau() const { return fromTau(true); }
/// @brief Determine whether the particle is from a hadron or tau decay
///
/// Specifically, walk up the ancestor chain until a status 2 hadron or
/// tau is found, if at all.
///
/// @note This question is valid in MC, but may not be perfectly answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
bool fromDecay() const { return fromHadron() || fromPromptTau(); }
//@}
/// @name Decay info
//@{
/// Whether this particle is stable according to the generator
bool isStable() const {
return genParticle() != NULL && genParticle()->status() == 1 && genParticle()->end_vertex() == NULL;
}
/// Get a list of the direct descendants from the current particle
vector<Particle> children() const {
vector<Particle> rtn;
if (isStable()) return rtn;
/// @todo Remove this const mess crap when HepMC doesn't suck
HepMC::GenVertex* gv = const_cast<HepMC::GenVertex*>( genParticle()->end_vertex() );
/// @todo Would like to do this, but the range objects are broken
// foreach (const GenParticle* 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)
rtn += Particle(*it);
return rtn;
}
/// Get a list of all the descendants (including duplication of parents and children) from the current particle
/// @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?
vector<Particle> allDescendants() const {
vector<Particle> rtn;
if (isStable()) return rtn;
/// @todo Remove this const mess crap when HepMC doesn't suck
HepMC::GenVertex* gv = const_cast<HepMC::GenVertex*>( genParticle()->end_vertex() );
/// @todo Would like to do this, but the range objects are broken
// foreach (const GenParticle* gp, gv->particles(HepMC::descendants))
// rtn += Particle(gp);
for (GenVertex::particle_iterator it = gv->particles_begin(HepMC::descendants); it != gv->particles_end(HepMC::descendants); ++it)
rtn += Particle(*it);
return rtn;
}
/// Get a list of all the stable descendants from the current particle
/// @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?
vector<Particle> stableDescendants() const {
vector<Particle> rtn;
if (isStable()) return rtn;
/// @todo Remove this const mess crap when HepMC doesn't suck
HepMC::GenVertex* gv = const_cast<HepMC::GenVertex*>( genParticle()->end_vertex() );
/// @todo Would like to do this, but the range objects are broken
// foreach (const GenParticle* gp, gv->particles(HepMC::descendants))
// if (gp->status() == 1 && gp->end_vertex() == NULL)
// rtn += Particle(gp);
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)
rtn += Particle(*it);
return rtn;
}
/// Flight length (divide by mm or cm to get the appropriate units)
double 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()));
}
//@}
private:
/// 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 Unbound functions for filtering particles
//@{
/// Filter a jet collection in-place to the subset that passes the supplied Cut
Particles& filterBy(Particles& particles, const Cut& c);
/// Get a subset of the supplied particles that passes the supplied Cut
Particles filterBy(const Particles& particles, const Cut& c);
//@}
/// @name String representation
//@{
/// Print a ParticlePair as a string.
inline std::string to_str(const ParticlePair& pair) {
stringstream out;
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.
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/Projections/Beam.hh b/include/Rivet/Projections/Beam.hh
--- a/include/Rivet/Projections/Beam.hh
+++ b/include/Rivet/Projections/Beam.hh
@@ -1,137 +1,147 @@
// -*- C++ -*-
#ifndef RIVET_Beam_HH
#define RIVET_Beam_HH
#include "Rivet/Projection.hh"
#include "Rivet/Event.hh"
#include "Rivet/Particle.hh"
+#include "Rivet/Math/LorentzTrans.hh"
namespace Rivet {
/// @name Standalone beam kinematics functions
//@{
/// Get beam particles from an event
ParticlePair beams(const Event& e);
/// Get beam particle IDs from a pair of Particles
inline PdgIdPair beamIds(const ParticlePair& beams) {
return make_pair(beams.first.pid(), beams.second.pid());
}
/// Get beam particle IDs from an event
inline PdgIdPair beamIds(const Event& e) {
return beamIds(beams(e));
}
/// Get beam centre-of-mass energy from a pair of beam momenta
double sqrtS(const FourMomentum& pa, const FourMomentum& pb);
/// Get beam centre-of-mass energy from a pair of Particles
inline double sqrtS(const ParticlePair& beams) {
return sqrtS(beams.first.momentum(), beams.second.momentum());
}
/// Get beam centre-of-mass energy from an Event
inline double sqrtS(const Event& e) {
return sqrtS(beams(e));
}
/// Get per-nucleon beam centre-of-mass energy from a pair of beam momenta
/// @note Uses a nominal nucleon mass of 0.939 GeV to convert masses to A
double asqrtS(const FourMomentum& pa, const FourMomentum& pb);
/// Get per-nucleon beam centre-of-mass energy from a pair of Particles
/// @note Uses the sum of nuclear mass numbers A for each beam
double asqrtS(const ParticlePair& beams);
/// Get per-nucleon beam centre-of-mass energy from an Event
inline double asqrtS(const Event& e) {
return asqrtS(beams(e));
}
- /// Get an active Lorentz boost to the beam centre-of-mass from a pair of beam momenta
- Vector3 comBoost(const FourMomentum& pa, const FourMomentum& pb);
+ /// Get the Lorentz boost to the beam centre-of-mass system (CMS) from a pair of beam momenta
+ Vector3 cmsBoost(const FourMomentum& pa, const FourMomentum& pb);
- /// Get an active Lorentz boost to the beam centre-of-mass from a pair of Particles
- inline Vector3 comBoost(const ParticlePair& beams) {
- return comBoost(beams.first, beams.second);
+ /// Get the Lorentz boost to the beam centre-of-mass system (CMS) from a pair of Particles
+ inline Vector3 cmsBoost(const ParticlePair& beams) {
+ return cmsBoost(beams.first, beams.second);
}
- /// Get an active Lorentz boost to the beam centre-of-mass from an Event
- inline Vector3 comBoost(const Event& e) {
- return comBoost(beams(e));
+ /// Get the Lorentz boost to the beam centre-of-mass system (CMS) from an Event
+ inline Vector3 beamCMSBoost(const Event& e) {
+ return cmsBoost(beams(e));
+ }
+
+
+ /// Get the Lorentz transformation to the beam centre-of-mass system (CMS) from a pair of beam momenta
+ LorentzTransform cmsTransform(const FourMomentum& pa, const FourMomentum& pb);
+
+ /// Get the Lorentz transformation to the beam centre-of-mass system (CMS) from a pair of Particles
+ inline LorentzTransform cmsTransform(const ParticlePair& beams) {
+ return cmsTransform(beams.first, beams.second);
+ }
+
+ /// Get the Lorentz transformation to the beam centre-of-mass system (CMS) from an Event
+ inline LorentzTransform beamCMSTransform(const Event& e) {
+ return cmsTransform(beams(e));
}
//@}
/// @brief Project out the incoming beams
class Beam : public Projection {
public:
/// The default constructor
Beam() {
setName("Beam");
}
/// Clone on the heap
virtual unique_ptr<Projection> clone() const {
return unique_ptr<Projection>(new Beam(*this));
}
public:
/// The pair of beam particles in the current collision
- const ParticlePair& beams() const {
- return _theBeams;
- }
+ const ParticlePair& beams() const { return _theBeams; }
/// The pair of beam particle PDG codes in the current collision
- const PdgIdPair beamIds() const {
- return Rivet::beamIds(beams());
- }
+ const PdgIdPair beamIds() const { return Rivet::beamIds(beams()); }
/// Get centre of mass energy, \f$ \sqrt{s} \f$
double sqrtS() const { return Rivet::sqrtS(beams()); }
/// Get per-nucleon centre of mass energy, \f$ \sqrt{s}/(A_1 + A_2) \f$
double asqrtS() const { return Rivet::asqrtS(beams()); }
- /// Get an active Lorentz boost to the beam centre-of-mass
- Vector3 comBoost() const { return Rivet::comBoost(beams()); }
+ /// Get the Lorentz boost to the beam centre-of-mass
+ Vector3 cmsBoost() const { return Rivet::cmsBoost(beams()); }
+
+ /// Get the Lorentz transform to the beam centre-of-mass
+ LorentzTransform cmsTransform() const { return Rivet::cmsTransform(beams()); }
/// Get the beam interaction primary vertex (PV) position
FourVector pv() const;
- public:
-
/// Project on to the Event
virtual void project(const Event& e);
private:
/// Compare with other projections -- it's always the same, since there are no params
- virtual int compare(const Projection& UNUSED(p)) const {
- return EQUIVALENT;
- }
+ virtual int compare(const Projection& UNUSED(p)) const { return EQUIVALENT; }
/// The beam particles in the current collision
ParticlePair _theBeams;
};
}
#endif
diff --git a/src/Core/AnalysisHandler.cc b/src/Core/AnalysisHandler.cc
--- a/src/Core/AnalysisHandler.cc
+++ b/src/Core/AnalysisHandler.cc
@@ -1,332 +1,332 @@
// -*- C++ -*-
#include "Rivet/Config/RivetCommon.hh"
#include "Rivet/AnalysisHandler.hh"
#include "Rivet/Analysis.hh"
#include "Rivet/ParticleName.hh"
#include "Rivet/BeamConstraint.hh"
#include "Rivet/Tools/Logging.hh"
#include "Rivet/Projections/Beam.hh"
#include "YODA/WriterYODA.h"
namespace Rivet {
AnalysisHandler::AnalysisHandler(const string& runname)
: _runname(runname), _numEvents(0),
_sumOfWeights(0.0), _xs(NAN),
_initialised(false), _ignoreBeams(false)
{}
AnalysisHandler::~AnalysisHandler()
{}
Log& AnalysisHandler::getLog() const {
return Log::getLog("Rivet.Analysis.Handler");
}
void AnalysisHandler::init(const GenEvent& ge) {
if (_initialised)
throw UserError("AnalysisHandler::init has already been called: cannot re-initialize!");
setRunBeams(Rivet::beams(ge));
MSG_DEBUG("Initialising the analysis handler");
_numEvents = 0;
_sumOfWeights = 0.0;
_sumOfWeightsSq = 0.0;
// Check that analyses are beam-compatible, and remove those that aren't
const size_t num_anas_requested = analysisNames().size();
vector<string> anamestodelete;
- foreach (const AnaHandle a, _analyses) {
+ for (const AnaHandle a : _analyses) {
if (!_ignoreBeams && !a->isCompatible(beams())) {
//MSG_DEBUG(a->name() << " requires beams " << a->requiredBeams() << " @ " << a->requiredEnergies() << " GeV");
anamestodelete.push_back(a->name());
}
}
- foreach (const string& aname, anamestodelete) {
+ for (const string& aname : anamestodelete) {
MSG_WARNING("Analysis '" << aname << "' is incompatible with the provided beams: removing");
removeAnalysis(aname);
}
if (num_anas_requested > 0 && analysisNames().empty()) {
cerr << "All analyses were incompatible with the first event's beams\n"
<< "Exiting, since this probably wasn't intentional!" << endl;
exit(1);
}
// Warn if any analysis' status is not unblemished
- foreach (const AnaHandle a, analyses()) {
+ for (const AnaHandle a : analyses()) {
if (toUpper(a->status()) == "PRELIMINARY") {
MSG_WARNING("Analysis '" << a->name() << "' is preliminary: be careful, it may change and/or be renamed!");
} else if (toUpper(a->status()) == "OBSOLETE") {
MSG_WARNING("Analysis '" << a->name() << "' is obsolete: please update!");
} else if (toUpper(a->status()).find("UNVALIDATED") != string::npos) {
MSG_WARNING("Analysis '" << a->name() << "' is unvalidated: be careful, it may be broken!");
}
}
// Initialize the remaining analyses
- foreach (AnaHandle a, _analyses) {
+ for (AnaHandle a : _analyses) {
MSG_DEBUG("Initialising analysis: " << a->name());
try {
// Allow projection registration in the init phase onwards
a->_allowProjReg = true;
a->init();
//MSG_DEBUG("Checking consistency of analysis: " << a->name());
//a->checkConsistency();
} catch (const Error& err) {
cerr << "Error in " << a->name() << "::init method: " << err.what() << endl;
exit(1);
}
MSG_DEBUG("Done initialising analysis: " << a->name());
}
_initialised = true;
MSG_DEBUG("Analysis handler initialised");
}
void AnalysisHandler::analyze(const GenEvent& ge) {
// Call init with event as template if not already initialised
if (!_initialised) init(ge);
assert(_initialised);
// Ensure that beam details match those from the first event
const PdgIdPair beams = Rivet::beamIds(ge);
const double sqrts = Rivet::sqrtS(ge);
if (!compatible(beams, _beams) || !fuzzyEquals(sqrts, sqrtS())) {
cerr << "Event beams mismatch: "
<< PID::toBeamsString(beams) << " @ " << sqrts/GeV << " GeV" << " vs. first beams "
<< this->beams() << " @ " << this->sqrtS()/GeV << " GeV" << endl;
exit(1);
}
// Create the Rivet event wrapper
/// @todo Filter/normalize the event here
Event event(ge);
// Weights
/// @todo Drop this / just report first weight when we support multiweight events
_numEvents += 1;
_sumOfWeights += event.weight();
_sumOfWeightsSq += sqr(event.weight());
MSG_DEBUG("Event #" << _numEvents << " weight = " << event.weight());
// Cross-section
#ifdef HEPMC_HAS_CROSS_SECTION
if (ge.cross_section()) {
_xs = ge.cross_section()->cross_section();
_xserr = ge.cross_section()->cross_section_error();
}
#endif
// Run the analyses
- foreach (AnaHandle a, _analyses) {
+ for (AnaHandle a : _analyses) {
MSG_TRACE("About to run analysis " << a->name());
try {
a->analyze(event);
} catch (const Error& err) {
cerr << "Error in " << a->name() << "::analyze method: " << err.what() << endl;
exit(1);
}
MSG_TRACE("Finished running analysis " << a->name());
}
}
void AnalysisHandler::analyze(const GenEvent* ge) {
if (ge == NULL) {
MSG_ERROR("AnalysisHandler received null pointer to GenEvent");
//throw Error("AnalysisHandler received null pointer to GenEvent");
}
analyze(*ge);
}
void AnalysisHandler::finalize() {
if (!_initialised) return;
MSG_INFO("Finalising analyses");
- foreach (AnaHandle a, _analyses) {
+ for (AnaHandle a : _analyses) {
a->setCrossSection(_xs);
try {
a->finalize();
} catch (const Error& err) {
cerr << "Error in " << a->name() << "::finalize method: " << err.what() << endl;
exit(1);
}
}
// Print out number of events processed
MSG_INFO("Processed " << _numEvents << " event" << (_numEvents == 1 ? "" : "s"));
// // Delete analyses
// MSG_DEBUG("Deleting analyses");
// _analyses.clear();
// Print out MCnet boilerplate
cout << endl;
cout << "The MCnet usage guidelines apply to Rivet: see http://www.montecarlonet.org/GUIDELINES" << endl;
cout << "Please acknowledge plots made with Rivet analyses, and cite arXiv:1003.0694 (http://arxiv.org/abs/1003.0694)" << endl;
}
AnalysisHandler& AnalysisHandler::addAnalysis(const string& analysisname) {
// Check for a duplicate analysis
/// @todo Might we want to be able to run an analysis twice, with different params?
/// Requires avoiding histo tree clashes, i.e. storing the histos on the analysis objects.
- foreach (const AnaHandle& a, _analyses) {
+ for (const AnaHandle& a : _analyses) {
if (a->name() == analysisname) {
MSG_WARNING("Analysis '" << analysisname << "' already registered: skipping duplicate");
return *this;
}
}
AnaHandle analysis( AnalysisLoader::getAnalysis(analysisname) );
if (analysis.get() != 0) { // < Check for null analysis.
MSG_DEBUG("Adding analysis '" << analysisname << "'");
analysis->_analysishandler = this;
_analyses.insert(analysis);
} else {
MSG_WARNING("Analysis '" << analysisname << "' not found.");
}
// MSG_WARNING(_analyses.size());
- // foreach (const AnaHandle& a, _analyses) MSG_WARNING(a->name());
+ // for (const AnaHandle& a : _analyses) MSG_WARNING(a->name());
return *this;
}
AnalysisHandler& AnalysisHandler::removeAnalysis(const string& analysisname) {
std::shared_ptr<Analysis> toremove;
- foreach (const AnaHandle a, _analyses) {
+ for (const AnaHandle a : _analyses) {
if (a->name() == analysisname) {
toremove = a;
break;
}
}
if (toremove.get() != 0) {
MSG_DEBUG("Removing analysis '" << analysisname << "'");
_analyses.erase(toremove);
}
return *this;
}
vector<AnalysisObjectPtr> AnalysisHandler::getData() const {
vector<AnalysisObjectPtr> rtn;
rtn.push_back( make_shared<Counter>(YODA::Dbn0D(_numEvents, _sumOfWeights, _sumOfWeightsSq), "/_EVTCOUNT") );
YODA::Scatter1D::Points pts; pts.insert(YODA::Point1D(_xs, _xserr));
rtn.push_back( make_shared<Scatter1D>(pts, "/_XSEC") );
- foreach (const AnaHandle a, analyses()) {
+ for (const AnaHandle a : analyses()) {
vector<AnalysisObjectPtr> aos = a->analysisObjects();
// MSG_WARNING(a->name() << " " << aos.size());
- foreach (const AnalysisObjectPtr ao, aos) {
+ for (const AnalysisObjectPtr ao : aos) {
// Exclude paths starting with /TMP/ from final write-out
/// @todo This needs to be much more nuanced for re-entrant histogramming
if (ao->path().find("/TMP/") != string::npos) continue;
rtn.push_back(ao);
}
}
- sort(rtn.begin(), rtn.end(),
+ sort(rtn.begin(), rtn.end(),
[](AnalysisObjectPtr a, AnalysisObjectPtr b) {
return a->path() < b->path();
}
);
return rtn;
}
void AnalysisHandler::writeData(const string& filename) const {
const vector<AnalysisObjectPtr> aos = getData();
try {
YODA::WriterYODA::write(filename, aos.begin(), aos.end());
} catch (...) { /// @todo Move to specific YODA::WriteError type when YODA >= 1.5.0 is well-established
throw UserError("Unexpected error in writing file to: " + filename);
}
}
string AnalysisHandler::runName() const { return _runname; }
size_t AnalysisHandler::numEvents() const { return _numEvents; }
double AnalysisHandler::sumOfWeights() const { return _sumOfWeights; }
void AnalysisHandler::setSumOfWeights(const double& sum) {
_sumOfWeights=sum;
}
std::vector<std::string> AnalysisHandler::analysisNames() const {
std::vector<std::string> rtn;
- foreach (AnaHandle a, _analyses) {
+ for (AnaHandle a : _analyses) {
rtn.push_back(a->name());
}
return rtn;
}
const AnaHandle AnalysisHandler::analysis(const std::string& analysisname) const {
- foreach (const AnaHandle a, analyses())
+ for (const AnaHandle a : analyses())
if (a->name() == analysisname) return a;
throw Error("No analysis named '" + analysisname + "' registered in AnalysisHandler");
}
AnalysisHandler& AnalysisHandler::addAnalyses(const std::vector<std::string>& analysisnames) {
- foreach (const string& aname, analysisnames) {
+ for (const string& aname : analysisnames) {
//MSG_DEBUG("Adding analysis '" << aname << "'");
addAnalysis(aname);
}
return *this;
}
AnalysisHandler& AnalysisHandler::removeAnalyses(const std::vector<std::string>& analysisnames) {
- foreach (const string& aname, analysisnames) {
+ for (const string& aname : analysisnames) {
removeAnalysis(aname);
}
return *this;
}
bool AnalysisHandler::needCrossSection() const {
bool rtn = false;
- foreach (const AnaHandle a, _analyses) {
+ for (const AnaHandle a : _analyses) {
if (!rtn) rtn = a->needsCrossSection();
if (rtn) break;
}
return rtn;
}
AnalysisHandler& AnalysisHandler::setCrossSection(double xs) {
_xs = xs;
return *this;
}
bool AnalysisHandler::hasCrossSection() const {
return (!std::isnan(crossSection()));
}
AnalysisHandler& AnalysisHandler::addAnalysis(Analysis* analysis) {
analysis->_analysishandler = this;
_analyses.insert(AnaHandle(analysis));
return *this;
}
PdgIdPair AnalysisHandler::beamIds() const {
return Rivet::beamIds(beams());
}
double AnalysisHandler::sqrtS() const {
return Rivet::sqrtS(beams());
}
void AnalysisHandler::setIgnoreBeams(bool ignore) {
_ignoreBeams=ignore;
}
}
diff --git a/src/Core/Event.cc b/src/Core/Event.cc
--- a/src/Core/Event.cc
+++ b/src/Core/Event.cc
@@ -1,85 +1,97 @@
#include "Rivet/Event.hh"
#include "Rivet/Tools/Logging.hh"
#include "Rivet/Projections/Beam.hh"
#include "Rivet/BeamConstraint.hh"
#include "HepMC/GenEvent.h"
namespace Rivet {
+ // ParticlePair Event::beams() const { return Rivet::beams(*this); }
+
+ // double Event::sqrtS() const { return Rivet::sqrtS(*this); }
+
+ // double Event::asqrtS() const { return Rivet::asqrtS(*this); }
+
+ // Vector3 Event::beamCMSBoost() const { return Rivet::beamCMSBoost(*this); }
+
+ // LorentzTransform Event::beamCMSTransform() const { return Rivet::beamCMSTransform(*this); }
+
+
+
void Event::_init(const GenEvent& ge) {
// Use Rivet's preferred units if possible
#ifdef HEPMC_HAS_UNITS
_genEvent.use_units(HepMC::Units::GEV, HepMC::Units::MM);
#endif
// Use the conventional alignment
// _geNormAlignment();
/// @todo Filter the event to remove generator-specific particles: optional
/// behaviour? Maybe disableable in an inconvenient way, e.g. with an env
/// var, to communicate the appropriate distaste for this sort of truth
/// analysis ;-)
// Debug printout to check that copying/mangling has worked
/// @todo Enable this when HepMC has been fixed to allow printing to a stream like the Rivet logger.
//_genEvent.print();
}
// namespace { // unnamed namespace for hiding
//
// void _geRot180x(GenEvent& ge) {
// /// @todo Use nicer iterators over HepMC particles
// for (HepMC::GenEvent::particle_iterator ip = ge.particles_begin(); ip != ge.particles_end(); ++ip) {
// const HepMC::FourVector& mom = (*ip)->momentum();
// (*ip)->set_momentum(HepMC::FourVector(mom.px(), -mom.py(), -mom.pz(), mom.e()));
// }
// /// @todo Use nicer iterators over HepMC vertices
// for (HepMC::GenEvent::vertex_iterator iv = ge.vertices_begin(); iv != ge.vertices_end(); ++iv) {
// const HepMC::FourVector& pos = (*iv)->position();
// (*iv)->set_position(HepMC::FourVector(pos.x(), -pos.y(), -pos.z(), pos.t()));
// }
// }
//
// }
// void Event::_geNormAlignment() {
// if (!_genEvent.valid_beam_particles()) return;
// typedef pair<HepMC::GenParticle*, HepMC::GenParticle*> GPPair;
// GPPair bps = _genEvent.beam_particles();
//
// // Rotate e+- p and ppbar to put p along +z
// /// @todo Is there an e+ e- convention for longitudinal boosting, e.g. at B-factories? Different from LEP?
// // if (compatible(beamids, make_pdgid_pair(ELECTRON, PROTON)) ||
// // compatible(beamids, make_pdgid_pair(POSITRON, PROTON)) ||
// // compatible(beamids, make_pdgid_pair(ANTIPROTON, PROTON)) ) {
// // Log::getLog("Rivet.Event") << Log::TRACE << "May need to rotate event..." << endl;
// bool rot = false;
// const HepMC::GenParticle* plusgp = 0;
// if (bps.first->pdg_id() != PID::PROTON || bps.second->pdg_id() != PID::PROTON) {
// if (bps.first->pdg_id() == PID::PROTON) {
// plusgp = bps.first;
// } else if (bps.second->pdg_id() == PID::PROTON) {
// plusgp = bps.second;
// }
// if (plusgp && plusgp->momentum().pz() < 0) {
// rot = true;
// }
// }
//
// // Do the rotation
// if (rot) {
// if (Log::getLog("Rivet.Event").isActive(Log::TRACE)) {
// Log::getLog("Rivet.Event") << Log::TRACE << "Rotating event\n"
// << "Before rotation: "
// << bps.first->pdg_id() << "@pz=" << bps.first->momentum().pz()/GeV << ", "
// << bps.second->pdg_id() << "@pz=" << bps.second->momentum().pz()/GeV << endl;
// }
// _geRot180x(_genEvent);
// }
// }
}
diff --git a/src/Projections/Beam.cc b/src/Projections/Beam.cc
--- a/src/Projections/Beam.cc
+++ b/src/Projections/Beam.cc
@@ -1,74 +1,78 @@
// -*- C++ -*-
#include "Rivet/Config/RivetCommon.hh"
#include "Rivet/Tools/Logging.hh"
#include "Rivet/Projections/Beam.hh"
namespace Rivet {
ParticlePair beams(const Event& e) {
assert(e.genEvent()->particles_size() >= 2);
if (e.genEvent()->valid_beam_particles()) {
pair<HepMC::GenParticle*, HepMC::GenParticle*> beams = e.genEvent()->beam_particles();
assert(beams.first && beams.second);
return ParticlePair{beams.first, beams.second};
} else 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)};
}
return ParticlePair{Particle(PID::ANY, FourMomentum()), Particle(PID::ANY, FourMomentum())};
}
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, pb) / (pa.mass()/MNUCLEON + pb.mass()/MNUCLEON);
}
double asqrtS(const ParticlePair& beams) {
return sqrtS(beams) / (nuclA(beams.first) + nuclA(beams.second));
}
- Vector3 comBoost(const FourMomentum& pa, const FourMomentum& pb) {
- Vector3 rtn = -(pa.p3() + pb.p3()) / (pa.E() + pb.E());
+ Vector3 cmsBoost(const FourMomentum& pa, const FourMomentum& pb) {
+ Vector3 rtn = (pa.p3() + pb.p3()) / (pa.E() + pb.E());
return rtn;
}
+ LorentzTransform cmsTransform(const FourMomentum& pa, const FourMomentum& pb) {
+ return LorentzTransform(-cmsBoost(pa, pb));
+ }
+
/////////////////////////////////////////////
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;
}
}