diff --git a/include/Rivet/Analysis.hh b/include/Rivet/Analysis.hh
--- a/include/Rivet/Analysis.hh
+++ b/include/Rivet/Analysis.hh
@@ -1,1027 +1,1031 @@
// -*- C++ -*-
#ifndef RIVET_Analysis_HH
#define RIVET_Analysis_HH
#include "Rivet/Config/RivetCommon.hh"
#include "Rivet/AnalysisInfo.hh"
#include "Rivet/Event.hh"
#include "Rivet/Projection.hh"
#include "Rivet/ProjectionApplier.hh"
#include "Rivet/ProjectionHandler.hh"
#include "Rivet/AnalysisLoader.hh"
#include "Rivet/Tools/RivetYODA.hh"
#include "Rivet/Tools/Logging.hh"
#include "Rivet/Tools/ParticleUtils.hh"
#include "Rivet/Tools/Cuts.hh"
/// @def vetoEvent
/// Preprocessor define for vetoing events, including the log message and return.
#define vetoEvent \
do { MSG_DEBUG("Vetoing event on line " << __LINE__ << " of " << __FILE__); return; } while(0)
namespace Rivet {
// Forward declaration
class AnalysisHandler;
/// @brief This is the base class of all analysis classes in Rivet.
///
/// There are
/// three virtual functions which should be implemented in base classes:
///
/// void init() is called by Rivet before a run is started. Here the
/// analysis class should book necessary histograms. The needed
/// projections should probably rather be constructed in the
/// constructor.
///
/// void analyze(const Event&) is called once for each event. Here the
/// analysis class should apply the necessary Projections and fill the
/// histograms.
///
/// void finalize() is called after a run is finished. Here the analysis
/// class should do whatever manipulations are necessary on the
/// histograms. Writing the histograms to a file is, however, done by
/// the Rivet class.
class Analysis : public ProjectionApplier {
/// The AnalysisHandler is a friend.
friend class AnalysisHandler;
public:
/// @name Standard constructors and destructors.
//@{
// /// The default constructor.
// Analysis();
/// Constructor
Analysis(const std::string& name);
/// The destructor.
virtual ~Analysis() {}
//@}
public:
/// @name Main analysis methods
//@{
/// Initialize this analysis object. A concrete class should here
/// book all necessary histograms. An overridden function must make
/// sure it first calls the base class function.
virtual void init() { }
/// Analyze one event. A concrete class should here apply the
/// necessary projections on the \a event and fill the relevant
/// histograms. An overridden function must make sure it first calls
/// the base class function.
virtual void analyze(const Event& event) = 0;
/// Finalize this analysis object. A concrete class should here make
/// all necessary operations on the histograms. Writing the
/// histograms to a file is, however, done by the Rivet class. An
/// overridden function must make sure it first calls the base class
/// function.
virtual void finalize() { }
//@}
public:
/// @name Metadata
/// Metadata is used for querying from the command line and also for
/// building web pages and the analysis pages in the Rivet manual.
//@{
/// Get the actual AnalysisInfo object in which all this metadata is stored.
const AnalysisInfo& info() const {
assert(_info && "No AnalysisInfo object :O");
return *_info;
}
/// @brief Get the name of the analysis.
///
/// By default this is computed by combining the results of the experiment,
/// year and Spires ID metadata methods and you should only override it if
/// there's a good reason why those won't work.
virtual std::string name() const {
return (info().name().empty()) ? _defaultname : info().name();
}
/// Get the Inspire ID code for this analysis.
virtual std::string inspireId() const {
return info().inspireId();
}
/// Get the SPIRES ID code for this analysis (~deprecated).
virtual std::string spiresId() const {
return info().spiresId();
}
/// @brief Names & emails of paper/analysis authors.
///
/// Names and email of authors in 'NAME \<EMAIL\>' format. The first
/// name in the list should be the primary contact person.
virtual std::vector<std::string> authors() const {
return info().authors();
}
/// @brief Get a short description of the analysis.
///
/// Short (one sentence) description used as an index entry.
/// Use @a description() to provide full descriptive paragraphs
/// of analysis details.
virtual std::string summary() const {
return info().summary();
}
/// @brief Get a full description of the analysis.
///
/// Full textual description of this analysis, what it is useful for,
/// what experimental techniques are applied, etc. Should be treated
/// as a chunk of restructuredText (http://docutils.sourceforge.net/rst.html),
/// with equations to be rendered as LaTeX with amsmath operators.
virtual std::string description() const {
return info().description();
}
/// @brief Information about the events needed as input for this analysis.
///
/// Event types, energies, kinematic cuts, particles to be considered
/// stable, etc. etc. Should be treated as a restructuredText bullet list
/// (http://docutils.sourceforge.net/rst.html)
virtual std::string runInfo() const {
return info().runInfo();
}
/// Experiment which performed and published this analysis.
virtual std::string experiment() const {
return info().experiment();
}
/// Collider on which the experiment ran.
virtual std::string collider() const {
return info().collider();
}
/// When the original experimental analysis was published.
virtual std::string year() const {
return info().year();
}
/// Journal, and preprint references.
virtual std::vector<std::string> references() const {
return info().references();
}
/// BibTeX citation key for this article.
virtual std::string bibKey() const {
return info().bibKey();
}
/// BibTeX citation entry for this article.
virtual std::string bibTeX() const {
return info().bibTeX();
}
/// Whether this analysis is trusted (in any way!)
virtual std::string status() const {
return (info().status().empty()) ? "UNVALIDATED" : info().status();
}
/// Any work to be done on this analysis.
virtual std::vector<std::string> todos() const {
return info().todos();
}
/// Return the allowed pairs of incoming beams required by this analysis.
virtual const std::vector<PdgIdPair>& requiredBeams() const {
return info().beams();
}
/// Declare the allowed pairs of incoming beams required by this analysis.
virtual Analysis& setRequiredBeams(const std::vector<PdgIdPair>& requiredBeams) {
info().setBeams(requiredBeams);
return *this;
}
/// Sets of valid beam energy pairs, in GeV
virtual const std::vector<std::pair<double, double> >& requiredEnergies() const {
return info().energies();
}
/// Declare the list of valid beam energy pairs, in GeV
virtual Analysis& setRequiredEnergies(const std::vector<std::pair<double, double> >& requiredEnergies) {
info().setEnergies(requiredEnergies);
return *this;
}
/// Return true if this analysis needs to know the process cross-section.
/// @todo Remove this and require HepMC >= 2.06
bool needsCrossSection() const {
return info().needsCrossSection();
}
/// Declare whether this analysis needs to know the process cross-section from the generator.
/// @todo Remove this and require HepMC >= 2.06
Analysis& setNeedsCrossSection(bool needed=true) {
info().setNeedsCrossSection(needed);
return *this;
}
//@}
/// @name Internal metadata modifying methods
//@{
/// Get the actual AnalysisInfo object in which all this metadata is stored (non-const).
AnalysisInfo& info() {
assert(_info && "No AnalysisInfo object :O");
return *_info;
}
//@}
/// @name Run conditions
//@{
/// Incoming beams for this run
const ParticlePair& beams() const;
/// Incoming beam IDs for this run
const PdgIdPair beamIds() const;
/// Centre of mass energy for this run
double sqrtS() const;
//@}
/// @name Analysis / beam compatibility testing
//@{
/// Check if analysis is compatible with the provided beam particle IDs and energies
bool isCompatible(const ParticlePair& beams) const;
/// Check if analysis is compatible with the provided beam particle IDs and energies
bool isCompatible(PdgId beam1, PdgId beam2, double e1, double e2) const;
/// Check if analysis is compatible with the provided beam particle IDs and energies
bool isCompatible(const PdgIdPair& beams, const std::pair<double,double>& energies) const;
//@}
/// Set the cross section from the generator
Analysis& setCrossSection(double xs);
/// Access the controlling AnalysisHandler object.
AnalysisHandler& handler() const { return *_analysishandler; }
protected:
/// Get a Log object based on the name() property of the calling analysis object.
Log& getLog() const;
/// Get the process cross-section in pb. Throws if this hasn't been set.
double crossSection() const;
/// Get the process cross-section per generated event in pb. Throws if this
/// hasn't been set.
vector<double> crossSectionPerEvent() const;
/// Get the number of events seen (via the analysis handler). Use in the
/// finalize phase only.
size_t numEvents() const;
/// Get the sum of event weights seen (via the analysis handler). Use in the
/// finalize phase only.
const vector<double>& sumOfWeights() const;
protected:
/// @name Histogram paths
//@{
/// Get the canonical histogram "directory" path for this analysis.
const std::string histoDir() const;
/// Get the canonical histogram path for the named histogram in this analysis.
const std::string histoPath(const std::string& hname) const;
/// Get the canonical histogram path for the numbered histogram in this analysis.
const std::string histoPath(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const;
/// Get the internal histogram name for given d, x and y (cf. HepData)
const std::string makeAxisCode(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const;
//@}
/// @name Histogram reference data
//@{
/// Get reference data for a named histo
/// @todo Move to the templated version when we have C++11 and can have a default fn template type
const YODA::Scatter2D& refData(const string& hname) const;
/// Get reference data for a numbered histo
/// @todo Move to the templated version when we have C++11 and can have a default fn template type
const YODA::Scatter2D& refData(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const;
/// Get reference data for a named histo
/// @todo Would be nice to just use these and ditch the S2D no-template version,
/// but we need C++11 for default args in function templates
// template <typename T=Scatter2D>
/// @todo SFINAE to ensure that the type inherits from YODA::AnalysisObject?
template <typename T>
const T& refData(const string& hname) const {
_cacheRefData();
MSG_TRACE("Using histo bin edges for " << name() << ":" << hname);
if (!_refdata[hname]) {
MSG_ERROR("Can't find reference histogram " << hname);
throw Exception("Reference data " + hname + " not found.");
}
return dynamic_cast<T&>(*_refdata[hname]);
}
/// Get reference data for a numbered histo
/// @todo Would be nice to just use these and ditch the S2D no-template version,
/// but we need C++11 for default args in function templates
// template <typename T=Scatter2D>
/// @todo SFINAE to ensure that the type inherits from YODA::AnalysisObject?
template <typename T>
const T& refData(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const {
const string hname = makeAxisCode(datasetId, xAxisId, yAxisId);
return refData(hname);
}
//@}
/// @name Counter booking
//@{
/// Book a counter.
CounterPtr bookCounter(const std::string& name,
const std::string& title="");
// const std::string& valtitle=""
/// Book a counter, using a path generated from the dataset and axis ID codes
///
/// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way.
CounterPtr bookCounter(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
const std::string& title="");
// const std::string& valtitle=""
//@}
/// @name 1D histogram booking
//@{
/// Book a 1D histogram with @a nbins uniformly distributed across the range @a lower - @a upper .
Histo1DPtr bookHisto1D(const std::string& name,
size_t nbins, double lower, double upper,
const std::string& title="",
const std::string& xtitle="",
const std::string& ytitle="");
/// Book a 1D histogram with non-uniform bins defined by the vector of bin edges @a binedges .
Histo1DPtr bookHisto1D(const std::string& name,
const std::vector<double>& binedges,
const std::string& title="",
const std::string& xtitle="",
const std::string& ytitle="");
/// Book a 1D histogram with binning from a reference scatter.
Histo1DPtr bookHisto1D(const std::string& name,
const Scatter2D& refscatter,
const std::string& title="",
const std::string& xtitle="",
const std::string& ytitle="");
/// Book a 1D histogram, using the binnings in the reference data histogram.
Histo1DPtr bookHisto1D(const std::string& name,
const std::string& title="",
const std::string& xtitle="",
const std::string& ytitle="");
/// Book a 1D histogram, using the binnings in the reference data histogram.
///
/// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way.
Histo1DPtr bookHisto1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
const std::string& title="",
const std::string& xtitle="",
const std::string& ytitle="");
//@}
/// @name 2D histogram booking
//@{
/// Book a 2D histogram with @a nxbins and @a nybins uniformly
/// distributed across the ranges @a xlower - @a xupper and @a
/// ylower - @a yupper respectively along the x- and y-axis.
Histo2DPtr bookHisto2D(const std::string& name,
size_t nxbins, double xlower, double xupper,
size_t nybins, double ylower, double yupper,
const std::string& title="",
const std::string& xtitle="",
const std::string& ytitle="",
const std::string& ztitle="");
/// Book a 2D histogram with non-uniform bins defined by the
/// vectorx of bin edges @a xbinedges and @a ybinedges.
Histo2DPtr bookHisto2D(const std::string& name,
const std::vector<double>& xbinedges,
const std::vector<double>& ybinedges,
const std::string& title="",
const std::string& xtitle="",
const std::string& ytitle="",
const std::string& ztitle="");
// /// Book a 2D histogram with binning from a reference scatter.
// Histo2DPtr bookHisto2D(const std::string& name,
// const Scatter3D& refscatter,
// const std::string& title="",
// const std::string& xtitle="",
// const std::string& ytitle="",
// const std::string& ztitle="");
// /// Book a 2D histogram, using the binnings in the reference data histogram.
// Histo2DPtr bookHisto2D(const std::string& name,
// const std::string& title="",
// const std::string& xtitle="",
// const std::string& ytitle="",
// const std::string& ztitle="");
// /// Book a 2D histogram, using the binnings in the reference data histogram.
// ///
// /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way.
// Histo2DPtr bookHisto2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
// const std::string& title="",
// const std::string& xtitle="",
// const std::string& ytitle="",
// const std::string& ztitle="");
//@}
/// @name 1D profile histogram booking
//@{
/// Book a 1D profile histogram with @a nbins uniformly distributed across the range @a lower - @a upper .
Profile1DPtr bookProfile1D(const std::string& name,
size_t nbins, double lower, double upper,
const std::string& title="",
const std::string& xtitle="",
const std::string& ytitle="");
/// Book a 1D profile histogram with non-uniform bins defined by the vector of bin edges @a binedges .
Profile1DPtr bookProfile1D(const std::string& name,
const std::vector<double>& binedges,
const std::string& title="",
const std::string& xtitle="",
const std::string& ytitle="");
/// Book a 1D profile histogram with binning from a reference scatter.
Profile1DPtr bookProfile1D(const std::string& name,
const Scatter2D& refscatter,
const std::string& title="",
const std::string& xtitle="",
const std::string& ytitle="");
/// Book a 1D profile histogram, using the binnings in the reference data histogram.
Profile1DPtr bookProfile1D(const std::string& name,
const std::string& title="",
const std::string& xtitle="",
const std::string& ytitle="");
/// Book a 1D profile histogram, using the binnings in the reference data histogram.
///
/// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way.
Profile1DPtr bookProfile1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
const std::string& title="",
const std::string& xtitle="",
const std::string& ytitle="");
//@}
/// @name 2D profile histogram booking
//@{
/// Book a 2D profile histogram with @a nxbins and @a nybins uniformly
/// distributed across the ranges @a xlower - @a xupper and @a ylower - @a
/// yupper respectively along the x- and y-axis.
Profile2DPtr bookProfile2D(const std::string& name,
size_t nxbins, double xlower, double xupper,
size_t nybins, double ylower, double yupper,
const std::string& title="",
const std::string& xtitle="",
const std::string& ytitle="",
const std::string& ztitle="");
/// Book a 2D profile histogram with non-uniform bins defined by the vectorx
/// of bin edges @a xbinedges and @a ybinedges.
Profile2DPtr bookProfile2D(const std::string& name,
const std::vector<double>& xbinedges,
const std::vector<double>& ybinedges,
const std::string& title="",
const std::string& xtitle="",
const std::string& ytitle="",
const std::string& ztitle="");
/// Book a 2D profile histogram with binning from a reference scatter.
// Profile2DPtr bookProfile2D(const std::string& name,
// const Scatter3D& refscatter,
// const std::string& title="",
// const std::string& xtitle="",
// const std::string& ytitle="",
// const std::string& ztitle="");
// /// Book a 2D profile histogram, using the binnings in the reference data histogram.
// Profile2DPtr bookProfile2D(const std::string& name,
// const std::string& title="",
// const std::string& xtitle="",
// const std::string& ytitle="",
// const std::string& ztitle="");
// /// Book a 2D profile histogram, using the binnings in the reference data histogram.
// ///
// /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way.
// Profile2DPtr bookProfile2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
// const std::string& title="",
// const std::string& xtitle="",
// const std::string& ytitle="",
// const std::string& ztitle="");
//@}
/// @name 2D scatter booking
//@{
/// @brief Book a 2-dimensional data point set with the given name.
///
/// @note Unlike histogram booking, scatter booking by default makes no
/// attempt to use reference data to pre-fill the data object. If you want
/// this, which is sometimes useful e.g. when the x-position is not really
/// meaningful and can't be extracted from the data, then set the @a
/// copy_pts parameter to true. This creates points to match the reference
/// data's x values and errors, but with the y values and errors zeroed...
/// assuming that there is a reference histo with the same name: if there
/// isn't, an exception will be thrown.
+
+ /*
Scatter2DPtr bookScatter2D(const std::string& name,
bool copy_pts=false,
const std::string& title="",
const std::string& xtitle="",
const std::string& ytitle="");
/// @brief Book a 2-dimensional data point set, using the binnings in the reference data histogram.
///
/// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way.
///
/// @note Unlike histogram booking, scatter booking by default makes no
/// attempt to use reference data to pre-fill the data object. If you want
/// this, which is sometimes useful e.g. when the x-position is not really
/// meaningful and can't be extracted from the data, then set the @a
/// copy_pts parameter to true. This creates points to match the reference
/// data's x values and errors, but with the y values and errors zeroed.
Scatter2DPtr bookScatter2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
bool copy_pts=false,
const std::string& title="",
const std::string& xtitle="",
const std::string& ytitle="");
/// @brief Book a 2-dimensional data point set with equally spaced x-points in a range.
///
/// The y values and errors will be set to 0.
Scatter2DPtr bookScatter2D(const std::string& name,
size_t npts, double lower, double upper,
const std::string& title="",
const std::string& xtitle="",
const std::string& ytitle="");
/// @brief Book a 2-dimensional data point set based on provided contiguous "bin edges".
///
/// The y values and errors will be set to 0.
Scatter2DPtr bookScatter2D(const std::string& hname,
const std::vector<double>& binedges,
const std::string& title,
const std::string& xtitle,
const std::string& ytitle);
//@}
+ */
public:
/// @name Analysis object manipulation
/// @todo Should really be protected: only public to keep BinnedHistogram happy for now...
//@{
/// Multiplicatively scale the given counter, @a cnt, by factor @s factor.
void scale(CounterPtr cnt, double factor);
/// Multiplicatively scale the given counters, @a cnts, by factor @s factor.
/// @note Constness intentional, if weird, to allow passing rvalue refs of smart ptrs (argh)
/// @todo Use SFINAE for a generic iterable of CounterPtrs
void scale(const std::vector<CounterPtr>& cnts, double factor) {
for (auto& c : cnts) scale(c, factor);
}
/// @todo YUCK!
template <std::size_t array_size>
void scale(const CounterPtr (&cnts)[array_size], double factor) {
// for (size_t i = 0; i < std::extent<decltype(cnts)>::value; ++i) scale(cnts[i], factor);
for (auto& c : cnts) scale(c, factor);
}
/// Normalize the given histogram, @a histo, to area = @a norm.
void normalize(Histo1DPtr histo, double norm=1.0, bool includeoverflows=true);
/// Normalize the given histograms, @a histos, to area = @a norm.
/// @note Constness intentional, if weird, to allow passing rvalue refs of smart ptrs (argh)
/// @todo Use SFINAE for a generic iterable of Histo1DPtrs
void normalize(const std::vector<Histo1DPtr>& histos, double norm=1.0, bool includeoverflows=true) {
for (auto& h : histos) normalize(h, norm, includeoverflows);
}
/// @todo YUCK!
template <std::size_t array_size>
void normalize(const Histo1DPtr (&histos)[array_size], double norm=1.0, bool includeoverflows=true) {
for (auto& h : histos) normalize(h, norm, includeoverflows);
}
/// Multiplicatively scale the given histogram, @a histo, by factor @s factor.
void scale(Histo1DPtr histo, double factor);
/// Multiplicatively scale the given histograms, @a histos, by factor @s factor.
/// @note Constness intentional, if weird, to allow passing rvalue refs of smart ptrs (argh)
/// @todo Use SFINAE for a generic iterable of Histo1DPtrs
void scale(const std::vector<Histo1DPtr>& histos, double factor) {
for (auto& h : histos) scale(h, factor);
}
/// @todo YUCK!
template <std::size_t array_size>
void scale(const Histo1DPtr (&histos)[array_size], double factor) {
for (auto& h : histos) scale(h, factor);
}
/// Normalize the given histogram, @a histo, to area = @a norm.
void normalize(Histo2DPtr histo, double norm=1.0, bool includeoverflows=true);
/// Normalize the given histograms, @a histos, to area = @a norm.
/// @note Constness intentional, if weird, to allow passing rvalue refs of smart ptrs (argh)
/// @todo Use SFINAE for a generic iterable of Histo2DPtrs
void normalize(const std::vector<Histo2DPtr>& histos, double norm=1.0, bool includeoverflows=true) {
for (auto& h : histos) normalize(h, norm, includeoverflows);
}
/// @todo YUCK!
template <std::size_t array_size>
void normalize(const Histo2DPtr (&histos)[array_size], double norm=1.0, bool includeoverflows=true) {
for (auto& h : histos) normalize(h, norm, includeoverflows);
}
/// Multiplicatively scale the given histogram, @a histo, by factor @s factor.
void scale(Histo2DPtr histo, double factor);
/// Multiplicatively scale the given histograms, @a histos, by factor @s factor.
/// @note Constness intentional, if weird, to allow passing rvalue refs of smart ptrs (argh)
/// @todo Use SFINAE for a generic iterable of Histo2DPtrs
void scale(const std::vector<Histo2DPtr>& histos, double factor) {
for (auto& h : histos) scale(h, factor);
}
/// @todo YUCK!
template <std::size_t array_size>
void scale(const Histo2DPtr (&histos)[array_size], double factor) {
for (auto& h : histos) scale(h, factor);
}
/// Helper for counter division.
///
/// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
void divide(CounterPtr c1, CounterPtr c2, Scatter1DPtr s) const;
/// Helper for histogram division with raw YODA objects.
///
/// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
void divide(const YODA::Counter& c1, const YODA::Counter& c2, Scatter1DPtr s) const;
/// Helper for histogram division.
///
/// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
void divide(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const;
/// Helper for histogram division with raw YODA objects.
///
/// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
void divide(const YODA::Histo1D& h1, const YODA::Histo1D& h2, Scatter2DPtr s) const;
/// Helper for profile histogram division.
///
/// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
void divide(Profile1DPtr p1, Profile1DPtr p2, Scatter2DPtr s) const;
/// Helper for profile histogram division with raw YODA objects.
///
/// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
void divide(const YODA::Profile1D& p1, const YODA::Profile1D& p2, Scatter2DPtr s) const;
/// Helper for 2D histogram division.
///
/// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
void divide(Histo2DPtr h1, Histo2DPtr h2, Scatter3DPtr s) const;
/// Helper for 2D histogram division with raw YODA objects.
///
/// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
void divide(const YODA::Histo2D& h1, const YODA::Histo2D& h2, Scatter3DPtr s) const;
/// Helper for 2D profile histogram division.
///
/// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
void divide(Profile2DPtr p1, Profile2DPtr p2, Scatter3DPtr s) const;
/// Helper for 2D profile histogram division with raw YODA objects
///
/// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
void divide(const YODA::Profile2D& p1, const YODA::Profile2D& p2, Scatter3DPtr s) const;
/// Helper for histogram efficiency calculation.
///
/// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
void efficiency(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const;
/// Helper for histogram efficiency calculation.
///
/// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
void efficiency(const YODA::Histo1D& h1, const YODA::Histo1D& h2, Scatter2DPtr s) const;
/// Helper for histogram asymmetry calculation.
///
/// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
void asymm(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const;
/// Helper for histogram asymmetry calculation.
///
/// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
void asymm(const YODA::Histo1D& h1, const YODA::Histo1D& h2, Scatter2DPtr s) const;
/// Helper for converting a differential histo to an integral one.
///
/// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
void integrate(Histo1DPtr h, Scatter2DPtr s) const;
/// Helper for converting a differential histo to an integral one.
///
/// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
void integrate(const Histo1D& h, Scatter2DPtr s) const;
//@}
public:
/// List of registered analysis data objects
- const vector<AnalysisObjectPtr>& analysisObjects() const {
+ const vector<shared_ptr<MultiweightAOPtr> >& analysisObjects() const {
return _analysisobjects;
}
protected:
/// @name Data object registration, retrieval, and removal
//@{
/// Register a data object in the histogram system
- void addAnalysisObject(AnalysisObjectPtr ao);
+ void addAnalysisObject(shared_ptr<MultiweightAOPtr> ao);
/// Get a data object from the histogram system
/// @todo Use this default function template arg in C++11
// template <typename AO=AnalysisObjectPtr>
template <typename AOPtr>
const AOPtr& getAnalysisObject(const std::string& name) const {
- foreach (const AnalysisObjectPtr& ao, analysisObjects()) {
- if (ao->path() == histoPath(name)) return dynamic_cast<AOPtr&>(ao);
+ for (const shared_ptr<MultiweightAOPtr>& ao : analysisObjects()) {
+ if ((*ao)->path() == histoPath(name)) return dynamic_cast<const AOPtr&>(*ao);
}
throw Exception("Data object " + histoPath(name) + " not found");
}
/// Get a data object from the histogram system (non-const)
/// @todo Use this default function template arg in C++11
// template <typename AO=AnalysisObjectPtr>
template <typename AOPtr>
AOPtr& getAnalysisObject(const std::string& name) {
- foreach (const AnalysisObjectPtr& ao, analysisObjects()) {
- if (ao->path() == histoPath(name)) return dynamic_cast<AOPtr&>(ao);
+ for (shared_ptr<MultiweightAOPtr>& ao : _analysisobjects) {
+ if ((*ao)->path() == histoPath(name)) return dynamic_cast<AOPtr&>(*ao);
}
+
throw Exception("Data object " + histoPath(name) + " not found");
}
/// Unregister a data object from the histogram system (by name)
void removeAnalysisObject(const std::string& path);
/// Unregister a data object from the histogram system (by pointer)
- void removeAnalysisObject(AnalysisObjectPtr ao);
+ void removeAnalysisObject(const MultiweightAOPtr& ao);
/// Get a named Histo1D object from the histogram system
const Histo1DPtr getHisto1D(const std::string& name) const {
return getAnalysisObject<Histo1DPtr>(name);
}
/// Get a named Histo1D object from the histogram system (non-const)
Histo1DPtr getHisto1D(const std::string& name) {
return getAnalysisObject<Histo1DPtr>(name);
}
/// Get a Histo1D object from the histogram system by axis ID codes (non-const)
const Histo1DPtr getHisto1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const {
return getAnalysisObject<Histo1DPtr>(makeAxisCode(datasetId, xAxisId, yAxisId));
}
/// Get a Histo1D object from the histogram system by axis ID codes (non-const)
Histo1DPtr getHisto1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) {
return getAnalysisObject<Histo1DPtr>(makeAxisCode(datasetId, xAxisId, yAxisId));
}
// /// Get a named Histo2D object from the histogram system
// const Histo2DPtr getHisto2D(const std::string& name) const {
// return getAnalysisObject<Histo2D>(name);
// }
// /// Get a named Histo2D object from the histogram system (non-const)
// Histo2DPtr getHisto2D(const std::string& name) {
// return getAnalysisObject<Histo2D>(name);
// }
// /// Get a Histo2D object from the histogram system by axis ID codes (non-const)
// const Histo2DPtr getHisto2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const {
// return getAnalysisObject<Histo2D>(makeAxisCode(datasetId, xAxisId, yAxisId));
// }
// /// Get a Histo2D object from the histogram system by axis ID codes (non-const)
// Histo2DPtr getHisto2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) {
// return getAnalysisObject<Histo2D>(makeAxisCode(datasetId, xAxisId, yAxisId));
// }
/// Get a named Profile1D object from the histogram system
const Profile1DPtr getProfile1D(const std::string& name) const {
return getAnalysisObject<Profile1DPtr>(name);
}
/// Get a named Profile1D object from the histogram system (non-const)
Profile1DPtr getProfile1D(const std::string& name) {
return getAnalysisObject<Profile1DPtr>(name);
}
/// Get a Profile1D object from the histogram system by axis ID codes (non-const)
const Profile1DPtr getProfile1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const {
return getAnalysisObject<Profile1DPtr>(makeAxisCode(datasetId, xAxisId, yAxisId));
}
/// Get a Profile1D object from the histogram system by axis ID codes (non-const)
Profile1DPtr getProfile1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) {
return getAnalysisObject<Profile1DPtr>(makeAxisCode(datasetId, xAxisId, yAxisId));
}
// /// Get a named Profile2D object from the histogram system
// const Profile2DPtr getProfile2D(const std::string& name) const {
// return getAnalysisObject<Profile2D>(name);
// }
// /// Get a named Profile2D object from the histogram system (non-const)
// Profile2DPtr getProfile2D(const std::string& name) {
// return getAnalysisObject<Profile2D>(name);
// }
// /// Get a Profile2D object from the histogram system by axis ID codes (non-const)
// const Profile2DPtr getProfile2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const {
// return getAnalysisObject<Profile2D>(makeAxisCode(datasetId, xAxisId, yAxisId));
// }
// /// Get a Profile2D object from the histogram system by axis ID codes (non-const)
// Profile2DPtr getProfile2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) {
// return getAnalysisObject<Profile2D>(makeAxisCode(datasetId, xAxisId, yAxisId));
// }
/// Get a named Scatter2D object from the histogram system
const Scatter2DPtr getScatter2D(const std::string& name) const {
return getAnalysisObject<Scatter2DPtr>(name);
}
/// Get a named Scatter2D object from the histogram system (non-const)
Scatter2DPtr getScatter2D(const std::string& name) {
return getAnalysisObject<Scatter2DPtr>(name);
}
/// Get a Scatter2D object from the histogram system by axis ID codes (non-const)
const Scatter2DPtr getScatter2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const {
return getAnalysisObject<Scatter2DPtr>(makeAxisCode(datasetId, xAxisId, yAxisId));
}
/// Get a Scatter2D object from the histogram system by axis ID codes (non-const)
Scatter2DPtr getScatter2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) {
return getAnalysisObject<Scatter2DPtr>(makeAxisCode(datasetId, xAxisId, yAxisId));
}
//@}
private:
/// Name passed to constructor (used to find .info analysis data file, and as a fallback)
string _defaultname;
/// Pointer to analysis metadata object
unique_ptr<AnalysisInfo> _info;
/// Storage of all plot objects
/// @todo Make this a map for fast lookup by path?
- vector<AnalysisObjectPtr> _analysisobjects;
+ vector<shared_ptr<MultiweightAOPtr> > _analysisobjects;
/// @name Cross-section variables
//@{
double _crossSection;
bool _gotCrossSection;
//@}
/// The controlling AnalysisHandler object.
AnalysisHandler* _analysishandler;
/// Collection of cached refdata to speed up many autobookings: the
/// reference data file should only be read once.
mutable std::map<std::string, YODA::AnalysisObjectPtr> _refdata;
private:
/// @name Utility functions
//@{
/// Get the reference data for this paper and cache it.
void _cacheRefData() const;
//@}
/// The assignment operator is private and must never be called.
/// In fact, it should not even be implemented.
Analysis& operator=(const Analysis&);
};
}
// Include definition of analysis plugin system so that analyses automatically see it when including Analysis.hh
#include "Rivet/AnalysisBuilder.hh"
/// @def DECLARE_RIVET_PLUGIN
/// Preprocessor define to prettify the global-object plugin hook mechanism.
#define DECLARE_RIVET_PLUGIN(clsname) Rivet::AnalysisBuilder<clsname> plugin_ ## clsname
/// @def DECLARE_ALIASED_RIVET_PLUGIN
/// Preprocessor define to prettify the global-object plugin hook mechanism, with an extra alias name for this analysis.
// #define DECLARE_ALIASED_RIVET_PLUGIN(clsname, alias) Rivet::AnalysisBuilder<clsname> plugin_ ## clsname ## ( ## #alias ## )
#define DECLARE_ALIASED_RIVET_PLUGIN(clsname, alias) DECLARE_RIVET_PLUGIN(clsname)( #alias )
/// @def DEFAULT_RIVET_ANA_CONSTRUCTOR
/// Preprocessor define to prettify the manky constructor with name string argument
#define DEFAULT_RIVET_ANALYSIS_CTOR(clsname) clsname() : Analysis(# clsname) {}
// DEPRECATED ALIAS
#define DEFAULT_RIVET_ANA_CONSTRUCTOR(clsname) DEFAULT_RIVET_ANALYSIS_CTOR(clsname)
#endif
diff --git a/include/Rivet/AnalysisHandler.hh b/include/Rivet/AnalysisHandler.hh
--- a/include/Rivet/AnalysisHandler.hh
+++ b/include/Rivet/AnalysisHandler.hh
@@ -1,248 +1,259 @@
// -*- 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.
const vector<double>& sumOfWeights() const {
return _sumOfWeights;
}
int numWeights() const {
return _sumOfWeights.size();
}
/// Set the active weight.
void setActiveWeight(unsigned int iWeight);
/// 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 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 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;
+ /// Get all analyses' plots as a vector of analysis objects.
+ void setWeightNames(const GenEvent& ge);
+
+ /// Do we have named weights?
+ bool haveNamedWeights();
+
+
/// 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
//@{
+ /// Weight names
+ std::vector<std::string> _weightNames;
+ std::vector<std::vector<double> > _subEventWeights;
+ size_t _numWeightTypes; // always == WeightVector.size()
+
/// 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
vector<double> _sumOfWeights;
/// 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;
/// Current event number
int _eventNumber;
//@}
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,174 +1,174 @@
// -*- C++ -*-
#ifndef RIVET_Event_HH
#define RIVET_Event_HH
#include "Rivet/Config/RivetCommon.hh"
#include "Rivet/Particle.hh"
#include "Rivet/Projection.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.
///
/// 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 Constructors and destructors.
//@{
/// Constructor from a HepMC GenEvent pointer
Event(const GenEvent* ge)
: _genevent_original(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)
: _genevent_original(&ge), _genevent(ge)
{ _init(ge); }
/// Copy constructor
Event(const Event& e)
: _genevent_original(e._genevent_original), _genevent(e._genevent)
{ }
//@}
/// @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; }
/// All the raw GenEvent particles, wrapped in Rivet::Particle objects
const Particles& allParticles() const;
/// @brief All the raw GenEvent particles, wrapped in Rivet::Particle objects, but with a Cut applied
///
/// @note Due to the cut, this returns by value, i.e. involves an expensive copy
inline Particles allParticles(const Cut& c) const {
return filter_select(allParticles(), c);
}
/// @brief All the raw GenEvent particles, wrapped in Rivet::Particle objects, but with a selection function applied
///
/// @note Due to the cut, this returns by value, i.e. involves an expensive copy
template <typename FN>
inline Particles allParticles(const FN& f) const {
return filter_select(allParticles(), f);
}
/// @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;
+ std::vector<double> weights() const;
//@}
/// @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
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 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* _genevent_original;
/// @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 neater this way.
/// @todo Change needed for HepMC3?
mutable GenEvent _genevent;
/// All the GenEvent particles, wrapped as Rivet::Particles
/// @note To be populated lazily, hence mutability
mutable Particles _particles;
/// The set of Projection objects applied so far
mutable std::set<ConstProjectionPtr> _projections;
};
}
#endif
diff --git a/include/Rivet/Tools/RivetYODA.hh b/include/Rivet/Tools/RivetYODA.hh
--- a/include/Rivet/Tools/RivetYODA.hh
+++ b/include/Rivet/Tools/RivetYODA.hh
@@ -1,303 +1,307 @@
#ifndef RIVET_RIVETYODA_HH
#define RIVET_RIVETYODA_HH
/// @author Andy Buckley
/// @date 2009-01-30
/// @author David Grellscheid
/// @date 2011-07-18
/// @author David Grellscheid
/// @date 2016-09-27
#include "Rivet/Config/RivetCommon.hh"
#include "YODA/AnalysisObject.h"
#include "YODA/Counter.h"
#include "YODA/Histo1D.h"
#include "YODA/Histo2D.h"
#include "YODA/Profile1D.h"
#include "YODA/Profile2D.h"
#include "YODA/Scatter1D.h"
#include "YODA/Scatter2D.h"
#include "YODA/Scatter3D.h"
#include <map>
namespace YODA {
typedef std::shared_ptr<YODA::AnalysisObject> AnalysisObjectPtr;
typedef std::shared_ptr<YODA::Histo1D> Histo1DPtr;
typedef std::shared_ptr<YODA::Histo2D> Histo2DPtr;
typedef std::shared_ptr<YODA::Profile1D> Profile1DPtr;
typedef std::shared_ptr<YODA::Profile2D> Profile2DPtr;
typedef std::shared_ptr<YODA::Scatter1D> Scatter1DPtr;
typedef std::shared_ptr<YODA::Scatter2D> Scatter2DPtr;
typedef std::shared_ptr<YODA::Scatter3D> Scatter3DPtr;
typedef std::shared_ptr<YODA::Counter> CounterPtr;
}
namespace Rivet {
class AnalysisObjectPtr {
public:
virtual ~AnalysisObjectPtr() {}
virtual YODA::AnalysisObject* operator->() = 0;
virtual YODA::AnalysisObject* operator->() const = 0;
virtual const YODA::AnalysisObject & operator*() const = 0;
+ bool operator ==(const AnalysisObjectPtr& p) { return (this == &p); }
+
protected:
/// @todo do we need this?
// virtual void reset() = 0;
};
/// @todo
/// implement scatter1dptr and scatter2dptr here
class Scatter1DPtr : public AnalysisObjectPtr {
public:
Scatter1DPtr() :
_scatter(YODA::Scatter1DPtr()) { }
- Scatter1DPtr(YODA::Scatter1DPtr& p) :
+ Scatter1DPtr(const YODA::Scatter1DPtr& p) :
_scatter(p) { }
bool operator!() const { return !_scatter; }
operator bool() const { return bool(_scatter); }
YODA::Scatter1D* operator->() { return _scatter.get(); }
YODA::Scatter1D* operator->() const { return _scatter.get(); }
YODA::Scatter1D & operator*() { return *_scatter; }
const YODA::Scatter1D & operator*() const { return *_scatter; }
protected:
YODA::Scatter1DPtr _scatter;
};
class Scatter2DPtr : public AnalysisObjectPtr {
public:
- Scatter2DPtr(YODA::Scatter2DPtr& p) :
+ Scatter2DPtr(const YODA::Scatter2DPtr& p) :
_scatter(p) { }
Scatter2DPtr() :
_scatter(YODA::Scatter2DPtr()) { }
bool operator!() { return !_scatter; }
operator bool() { return bool(_scatter); }
YODA::Scatter2D* operator->() { return _scatter.get(); }
YODA::Scatter2D* operator->() const { return _scatter.get(); }
YODA::Scatter2D & operator*() { return *_scatter; }
const YODA::Scatter2D & operator*() const { return *_scatter; }
protected:
YODA::Scatter2DPtr _scatter;
};
class Scatter3DPtr : public AnalysisObjectPtr {
public:
- Scatter3DPtr(YODA::Scatter3DPtr& p) :
+ Scatter3DPtr(const YODA::Scatter3DPtr& p) :
_scatter(p) { }
Scatter3DPtr() :
_scatter(YODA::Scatter3DPtr()) { }
bool operator!() { return !_scatter; }
operator bool() { return bool(_scatter); }
YODA::Scatter3D* operator->() { return _scatter.get(); }
YODA::Scatter3D* operator->() const { return _scatter.get(); }
YODA::Scatter3D & operator*() { return *_scatter; }
const YODA::Scatter3D & operator*() const { return *_scatter; }
protected:
YODA::Scatter3DPtr _scatter;
};
class MultiweightAOPtr : public AnalysisObjectPtr {
- protected:
+ public:
virtual void newSubEvent() = 0;
/// @todo
/// rename to setActive(Idx)?
virtual void setActiveWeightIdx(unsigned int iWeight) = 0;
virtual void pushToPersistent(const vector<vector<double> >& weight) = 0;
};
#define RIVETAOPTR_COMMON(YODATYPE) \
class YODATYPE##Ptr : public MultiweightAOPtr { \
public: \
\
/* @todo \
* some things are not really well-defined here \
* for instance: fill() in the finalize() method and integral() in \
* the analyze() method. \
*/ \
\
YODATYPE##Ptr(size_t len_of_weightvec, const YODA::YODATYPE& p) { \
for (size_t i = 0; i < len_of_weightvec; i++) \
_persistent.push_back(make_shared<YODA::YODATYPE>(p)); \
\
return; \
} \
\
YODA::YODATYPE##Ptr active() const { return _active; } \
\
/* @todo this probably need to loop over all? */ \
bool operator!() const {return !active();} \
operator bool() const {return bool(active());} \
\
YODA::YODATYPE* operator->() { \
return active().get(); \
} \
\
YODA::YODATYPE* operator->() const { \
return active().get(); \
} \
\
YODA::YODATYPE & operator*() { \
return *active(); \
} \
\
const YODA::YODATYPE & operator*() const { \
return *active(); \
} \
\
/* @todo \
* these need to be re-thought out. \
\
void reset() { active()->reset(); } \
*/ \
\
/* @todo \
* these probably need to loop over all? \
* do we even want to provide equality? \
*/ \
/* @todo \
* how about no. \
friend bool operator==(YODATYPE##Ptr a, YODATYPE##Ptr b){ \
if (a._persistent.size() != b._persistent.size()) \
return false; \
\
for (size_t i = 0; i < a._persistent.size(); i++) { \
if (a._persistent.at(i) != b._persistent.at(i)) { \
return false; \
} \
} \
\
return true; \
} \
\
friend bool operator!=(YODATYPE##Ptr a, YODATYPE##Ptr b){ \
return !(a == b); \
} \
\
\
friend bool operator<(YODATYPE##Ptr a, YODATYPE##Ptr b){ \
if (a._persistent.size() >= b._persistent.size()) \
return false; \
\
for (size_t i = 0; i < a._persistent.size(); i++) { \
if (*(a._persistent.at(i)) >= *(b._persistent.at(i))) { \
return false; \
} \
} \
\
return true; \
} \
*/ \
\
+ \
private: \
- \
void setActiveWeightIdx(unsigned int iWeight) { \
_active = _persistent.at(iWeight); \
return; \
} \
\
/* this is for dev only---we shouldn't need this in real runs. */ \
void unsetActiveWeight() { \
_active.reset(); \
return; \
} \
\
void newSubEvent() { \
YODA::YODATYPE##Ptr tmp = \
make_shared<YODA::YODATYPE>(_persistent[0]->clone()); \
tmp->reset(); \
_evgroup.push_back( tmp ); \
_active = _evgroup.back(); \
return; \
} \
\
const vector<YODA::YODATYPE##Ptr> & persistent() const { \
return _persistent; \
} \
\
/* to be implemented for each type */ \
void pushToPersistent(const vector<vector<double> >& weight); \
\
/* M of these, one for each weight */ \
vector<YODA::YODATYPE##Ptr> _persistent; \
\
/* N of these, one for each event in evgroup */ \
vector<YODA::YODATYPE##Ptr> _evgroup; \
\
YODA::YODATYPE##Ptr _active; \
+ \
+ friend class AnalysisHandler; \
};
RIVETAOPTR_COMMON(Histo1D)
RIVETAOPTR_COMMON(Histo2D)
RIVETAOPTR_COMMON(Profile1D)
RIVETAOPTR_COMMON(Profile2D)
RIVETAOPTR_COMMON(Counter)
using YODA::Counter;
using YODA::Histo1D;
using YODA::HistoBin1D;
using YODA::Histo2D;
using YODA::HistoBin2D;
using YODA::Profile1D;
using YODA::ProfileBin1D;
using YODA::Profile2D;
using YODA::ProfileBin2D;
using YODA::Scatter1D;
using YODA::Point1D;
using YODA::Scatter2D;
using YODA::Point2D;
using YODA::Scatter3D;
using YODA::Point3D;
/// Function to get a map of all the refdata in a paper with the
/// given @a papername.
map<string, YODA::AnalysisObjectPtr> getRefData(const string& papername);
/// @todo Also provide a Scatter3D getRefData() version?
/// Get the file system path to the reference file for this paper.
string getDatafilePath(const string& papername);
/// Return the integral over the histogram bins
/// @deprecated Prefer to directly use the histo's integral() method.
DEPRECATED("Prefer to directly use the histo's integral() method.")
inline double integral(Histo1DPtr histo) {
return histo->integral();
}
}
#endif
diff --git a/src/Core/Analysis.cc b/src/Core/Analysis.cc
--- a/src/Core/Analysis.cc
+++ b/src/Core/Analysis.cc
@@ -1,822 +1,780 @@
// -*- C++ -*-
#include "Rivet/Config/RivetCommon.hh"
#include "Rivet/Analysis.hh"
#include "Rivet/AnalysisHandler.hh"
#include "Rivet/AnalysisInfo.hh"
#include "Rivet/Tools/BeamConstraint.hh"
namespace Rivet {
Analysis::Analysis(const string& name)
: _crossSection(-1.0),
_gotCrossSection(false),
_analysishandler(NULL)
{
ProjectionApplier::_allowProjReg = false;
_defaultname = name;
unique_ptr<AnalysisInfo> ai = AnalysisInfo::make(name);
assert(ai);
_info = move(ai);
assert(_info);
}
double Analysis::sqrtS() const {
return handler().sqrtS();
}
const ParticlePair& Analysis::beams() const {
return handler().beams();
}
const PdgIdPair Analysis::beamIds() const {
return handler().beamIds();
}
const string Analysis::histoDir() const {
/// @todo Cache in a member variable
string _histoDir;
if (_histoDir.empty()) {
_histoDir = "/" + name();
if (handler().runName().length() > 0) {
_histoDir = "/" + handler().runName() + _histoDir;
}
replace_all(_histoDir, "//", "/"); //< iterates until none
}
return _histoDir;
}
const string Analysis::histoPath(const string& hname) const {
const string path = histoDir() + "/" + hname;
return path;
}
const string Analysis::histoPath(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const {
return histoDir() + "/" + makeAxisCode(datasetId, xAxisId, yAxisId);
}
const string Analysis::makeAxisCode(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const {
stringstream axisCode;
axisCode << "d";
if (datasetId < 10) axisCode << 0;
axisCode << datasetId;
axisCode << "-x";
if (xAxisId < 10) axisCode << 0;
axisCode << xAxisId;
axisCode << "-y";
if (yAxisId < 10) axisCode << 0;
axisCode << yAxisId;
return axisCode.str();
}
Log& Analysis::getLog() const {
string logname = "Rivet.Analysis." + name();
return Log::getLog(logname);
}
size_t Analysis::numEvents() const {
return handler().numEvents();
}
const vector<double>& Analysis::sumOfWeights() const {
return handler().sumOfWeights();
}
///////////////////////////////////////////
bool Analysis::isCompatible(const ParticlePair& beams) const {
return isCompatible(beams.first.pid(), beams.second.pid(),
beams.first.energy(), beams.second.energy());
}
bool Analysis::isCompatible(PdgId beam1, PdgId beam2, double e1, double e2) const {
PdgIdPair beams(beam1, beam2);
pair<double,double> energies(e1, e2);
return isCompatible(beams, energies);
}
bool Analysis::isCompatible(const PdgIdPair& beams, const pair<double,double>& energies) const {
// First check the beam IDs
bool beamIdsOk = false;
foreach (const PdgIdPair& bp, requiredBeams()) {
if (compatible(beams, bp)) {
beamIdsOk = true;
break;
}
}
if (!beamIdsOk) return false;
// Next check that the energies are compatible (within 1% or 1 GeV, whichever is larger, for a bit of UI forgiveness)
/// @todo Use some sort of standard ordering to improve comparisons, esp. when the two beams are different particles
bool beamEnergiesOk = requiredEnergies().size() > 0 ? false : true;
typedef pair<double,double> DoublePair;
foreach (const DoublePair& ep, requiredEnergies()) {
if ((fuzzyEquals(ep.first, energies.first, 0.01) && fuzzyEquals(ep.second, energies.second, 0.01)) ||
(fuzzyEquals(ep.first, energies.second, 0.01) && fuzzyEquals(ep.second, energies.first, 0.01)) ||
(abs(ep.first - energies.first) < 1*GeV && abs(ep.second - energies.second) < 1*GeV) ||
(abs(ep.first - energies.second) < 1*GeV && abs(ep.second - energies.first) < 1*GeV)) {
beamEnergiesOk = true;
break;
}
}
return beamEnergiesOk;
/// @todo Need to also check internal consistency of the analysis'
/// beam requirements with those of the projections it uses.
}
///////////////////////////////////////////
Analysis& Analysis::setCrossSection(double xs) {
_crossSection = xs;
_gotCrossSection = true;
return *this;
}
double Analysis::crossSection() const {
if (!_gotCrossSection || std::isnan(_crossSection)) {
string errMsg = "You did not set the cross section for the analysis " + name();
throw Error(errMsg);
}
return _crossSection;
}
vector<double> Analysis::crossSectionPerEvent() const {
const vector<double>& sumW = sumOfWeights();
assert(sumW.size() != 0.0);
// @todo
// is this correct?
vector<double> v(sumW);
foreach (double& x, v) {
x = _crossSection/x;
}
return v;
}
////////////////////////////////////////////////////////////
// Histogramming
void Analysis::_cacheRefData() const {
if (_refdata.empty()) {
MSG_TRACE("Getting refdata cache for paper " << name());
_refdata = getRefData(name());
}
}
const Scatter2D& Analysis::refData(const string& hname) const {
_cacheRefData();
MSG_TRACE("Using histo bin edges for " << name() << ":" << hname);
if (!_refdata[hname]) {
MSG_ERROR("Can't find reference histogram " << hname);
throw Exception("Reference data " + hname + " not found.");
}
return dynamic_cast<Scatter2D&>(*_refdata[hname]);
}
const Scatter2D& Analysis::refData(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const {
const string hname = makeAxisCode(datasetId, xAxisId, yAxisId);
return refData(hname);
}
CounterPtr Analysis::bookCounter(const string& cname,
const string& title) {
const string path = histoPath(cname);
- CounterPtr ctr(handler().numWeights(), Counter(path, title));
+ shared_ptr<CounterPtr> ctr =
+ make_shared<CounterPtr>(handler().numWeights(), Counter(path, title));
addAnalysisObject(ctr);
MSG_TRACE("Made counter " << cname << " for " << name());
- return ctr;
+ return *ctr;
}
CounterPtr Analysis::bookCounter(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
const string& title) {
// const string& xtitle,
// const string& ytitle) {
const string axisCode = makeAxisCode(datasetId, xAxisId, yAxisId);
return bookCounter(axisCode, title);
}
Histo1DPtr Analysis::bookHisto1D(const string& hname,
size_t nbins, double lower, double upper,
const string& title,
const string& xtitle,
const string& ytitle) {
const string path = histoPath(hname);
- Histo1DPtr hist(handler().numWeights(), Histo1D(nbins, lower, upper, path, title));
- addAnalysisObject(hist);
+ shared_ptr<Histo1DPtr> hp =
+ make_shared<Histo1DPtr>(handler().numWeights(), Histo1D(nbins, lower, upper, path, title));
+ addAnalysisObject(hp);
+ Histo1DPtr hist = *hp;
MSG_TRACE("Made histogram " << hname << " for " << name());
hist->setAnnotation("XLabel", xtitle);
hist->setAnnotation("YLabel", ytitle);
return hist;
}
Histo1DPtr Analysis::bookHisto1D(const string& hname,
const vector<double>& binedges,
const string& title,
const string& xtitle,
const string& ytitle) {
const string path = histoPath(hname);
- Histo1DPtr hist(handler().numWeights(), Histo1D(binedges, path, title));
- addAnalysisObject(hist);
+ shared_ptr<Histo1DPtr> hp =
+ make_shared<Histo1DPtr>(handler().numWeights(), Histo1D(binedges, path, title));
+ addAnalysisObject(hp);
+
+ Histo1DPtr hist = *hp;
MSG_TRACE("Made histogram " << hname << " for " << name());
hist->setAnnotation("XLabel", xtitle);
hist->setAnnotation("YLabel", ytitle);
return hist;
}
Histo1DPtr Analysis::bookHisto1D(const string& hname,
const Scatter2D& refscatter,
const string& title,
const string& xtitle,
const string& ytitle) {
const string path = histoPath(hname);
- Histo1DPtr hist(handler().numWeights(), Histo1D(refscatter, path));
- addAnalysisObject(hist);
+ shared_ptr<Histo1DPtr> hp =
+ make_shared<Histo1DPtr>(handler().numWeights(), Histo1D(refscatter, path));
+ addAnalysisObject(hp);
+
+ Histo1DPtr hist = *hp;
MSG_TRACE("Made histogram " << hname << " for " << name());
hist->setTitle(title);
hist->setAnnotation("XLabel", xtitle);
hist->setAnnotation("YLabel", ytitle);
return hist;
}
Histo1DPtr Analysis::bookHisto1D(const string& hname,
const string& title,
const string& xtitle,
const string& ytitle) {
const Scatter2D& refdata = refData(hname);
return bookHisto1D(hname, refdata, title, xtitle, ytitle);
}
Histo1DPtr Analysis::bookHisto1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
const string& title,
const string& xtitle,
const string& ytitle) {
const string axisCode = makeAxisCode(datasetId, xAxisId, yAxisId);
return bookHisto1D(axisCode, title, xtitle, ytitle);
}
/// @todo Add booking methods which take a path, titles and *a reference Scatter from which to book*
/////////////////
Histo2DPtr Analysis::bookHisto2D(const string& hname,
size_t nxbins, double xlower, double xupper,
size_t nybins, double ylower, double yupper,
const string& title,
const string& xtitle,
const string& ytitle,
const string& ztitle)
{
const string path = histoPath(hname);
- Histo2DPtr hist(handler().numWeights(), Histo2D(nxbins, xlower, xupper, nybins, ylower, yupper, path, title));
- addAnalysisObject(hist);
+ shared_ptr<Histo2DPtr> hp =
+ make_shared<Histo2DPtr>(handler().numWeights(), Histo2D(nxbins, xlower, xupper, nybins, ylower, yupper, path, title));
+ addAnalysisObject(hp);
+
+ Histo2DPtr hist = *hp;
MSG_TRACE("Made 2D histogram " << hname << " for " << name());
hist->setAnnotation("XLabel", xtitle);
hist->setAnnotation("YLabel", ytitle);
hist->setAnnotation("ZLabel", ztitle);
return hist;
}
Histo2DPtr Analysis::bookHisto2D(const string& hname,
const vector<double>& xbinedges,
const vector<double>& ybinedges,
const string& title,
const string& xtitle,
const string& ytitle,
const string& ztitle)
{
const string path = histoPath(hname);
- Histo2DPtr hist(handler().numWeights(), Histo2D(xbinedges, ybinedges, path, title));
- addAnalysisObject(hist);
+ shared_ptr<Histo2DPtr> hp =
+ make_shared<Histo2DPtr>(handler().numWeights(), Histo2D(xbinedges, ybinedges, path, title));
+ addAnalysisObject(hp);
+
+ Histo2DPtr hist = *hp;
MSG_TRACE("Made 2D histogram " << hname << " for " << name());
hist->setAnnotation("XLabel", xtitle);
hist->setAnnotation("YLabel", ytitle);
hist->setAnnotation("ZLabel", ztitle);
return hist;
}
- // Histo2DPtr Analysis::bookHisto2D(const string& hname,
- // const Scatter3D& refscatter,
- // const string& title="",
- // const string& xtitle="",
- // const string& ytitle="",
- // const string& ztitle="") {
- // const string path = histoPath(hname);
- // Histo2DPtr hist( new Histo2D(refscatter, path) );
- // addAnalysisObject(hist);
- // MSG_TRACE("Made 2D histogram " << hname << " for " << name());
- // hist->setTitle(title);
- // hist->setAnnotation("XLabel", xtitle);
- // hist->setAnnotation("YLabel", ytitle);
- // hist->setAnnotation("ZLabel", ztitle);
- // return hist;
- // }
-
-
- // Histo2DPtr Analysis::bookHisto2D(const string& hname,
- // const string& title,
- // const string& xtitle,
- // const string& ytitle,
- // const string& ztitle) {
- // const Scatter3D& refdata = refData(hname);
- // return bookHisto2D(hname, refdata, title, xtitle, ytitle, ztitle);
- // }
-
-
- // Histo2DPtr Analysis::bookHisto2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
- // const string& title,
- // const string& xtitle,
- // const string& ytitle,
- // const string& ztitle) {
- // const string axisCode = makeAxisCode(datasetId, xAxisId, yAxisId);
- // return bookHisto2D(axisCode, title, xtitle, ytitle, ztitle);
- // }
-
-
- /////////////////
-
-
Profile1DPtr Analysis::bookProfile1D(const string& hname,
size_t nbins, double lower, double upper,
const string& title,
const string& xtitle,
const string& ytitle) {
const string path = histoPath(hname);
- Profile1DPtr prof(handler().numWeights(), Profile1D(nbins, lower, upper, path, title));
- addAnalysisObject(prof);
+ shared_ptr<Profile1DPtr> pp =
+ make_shared<Profile1DPtr>(handler().numWeights(), Profile1D(nbins, lower, upper, path, title));
+ addAnalysisObject(pp);
+
+ Profile1DPtr prof = *pp;
MSG_TRACE("Made profile histogram " << hname << " for " << name());
prof->setAnnotation("XLabel", xtitle);
prof->setAnnotation("YLabel", ytitle);
return prof;
}
Profile1DPtr Analysis::bookProfile1D(const string& hname,
const vector<double>& binedges,
const string& title,
const string& xtitle,
const string& ytitle) {
const string path = histoPath(hname);
- Profile1DPtr prof(handler().numWeights(), Profile1D(binedges, path, title));
- addAnalysisObject(prof);
+ shared_ptr<Profile1DPtr> pp =
+ make_shared<Profile1DPtr>(handler().numWeights(), Profile1D(binedges, path, title));
+ addAnalysisObject(pp);
+
+ Profile1DPtr prof = *pp;
MSG_TRACE("Made profile histogram " << hname << " for " << name());
prof->setAnnotation("XLabel", xtitle);
prof->setAnnotation("YLabel", ytitle);
return prof;
}
Profile1DPtr Analysis::bookProfile1D(const string& hname,
const Scatter2D& refscatter,
const string& title,
const string& xtitle,
const string& ytitle) {
const string path = histoPath(hname);
- Profile1DPtr prof(handler().numWeights(), Profile1D(refscatter, path));
- addAnalysisObject(prof);
+ shared_ptr<Profile1DPtr> pp =
+ make_shared<Profile1DPtr>(handler().numWeights(), Profile1D(refscatter, path));
+ addAnalysisObject(pp);
+
+ Profile1DPtr prof = *pp;
MSG_TRACE("Made profile histogram " << hname << " for " << name());
prof->setTitle(title);
prof->setAnnotation("XLabel", xtitle);
prof->setAnnotation("YLabel", ytitle);
return prof;
}
Profile1DPtr Analysis::bookProfile1D(const string& hname,
const string& title,
const string& xtitle,
const string& ytitle) {
const Scatter2D& refdata = refData(hname);
return bookProfile1D(hname, refdata, title, xtitle, ytitle);
}
Profile1DPtr Analysis::bookProfile1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
const string& title,
const string& xtitle,
const string& ytitle) {
const string axisCode = makeAxisCode(datasetId, xAxisId, yAxisId);
return bookProfile1D(axisCode, title, xtitle, ytitle);
}
///////////////////
+ /*
Profile2DPtr Analysis::bookProfile2D(const string& hname,
size_t nxbins, double xlower, double xupper,
size_t nybins, double ylower, double yupper,
const string& title,
const string& xtitle,
const string& ytitle,
const string& ztitle)
{
const string path = histoPath(hname);
- Profile2DPtr prof(handler().numWeights(), Profile2D(nxbins, xlower, xupper, nybins, ylower, yupper, path, title));
- addAnalysisObject(prof);
+ shared_ptr<Profile2DPtr> pp =
+ make_shared<Profile2DPtr>(handler().numWeights(), Profile2D(nxbins, xlower, xupper, nybins, ylower, yupper, path, title));
+ addAnalysisObject(pp);
+
+ Profile2DPtr prof = *pp;
MSG_TRACE("Made 2D profile histogram " << hname << " for " << name());
prof->setAnnotation("XLabel", xtitle);
prof->setAnnotation("YLabel", ytitle);
prof->setAnnotation("ZLabel", ztitle);
return prof;
}
Profile2DPtr Analysis::bookProfile2D(const string& hname,
const vector<double>& xbinedges,
const vector<double>& ybinedges,
const string& title,
const string& xtitle,
const string& ytitle,
const string& ztitle)
{
const string path = histoPath(hname);
- Profile2DPtr prof(handler().numWeights(), Profile2D(xbinedges, ybinedges, path, title));
- addAnalysisObject(prof);
+ shared_ptr<Profile2DPtr> pp =
+ make_shared<Profile2DPtr>(handler().numWeights(), Profile2D(xbinedges, ybinedges, path, title));
+ addAnalysisObject(pp);
+
+ Profile2DPtr prof = *pp;
MSG_TRACE("Made 2D profile histogram " << hname << " for " << name());
prof->setAnnotation("XLabel", xtitle);
prof->setAnnotation("YLabel", ytitle);
prof->setAnnotation("ZLabel", ztitle);
return prof;
}
- // Profile2DPtr Analysis::bookProfile2D(const string& hname,
- // const Scatter3D& refscatter,
- // const string& title="",
- // const string& xtitle="",
- // const string& ytitle="",
- // const string& ztitle="") {
- // const string path = histoPath(hname);
- // Profile2DPtr prof( new Profile2D(refscatter, path) );
- // addAnalysisObject(prof);
- // MSG_TRACE("Made 2D profile histogram " << hname << " for " << name());
- // prof->setTitle(title);
- // prof->setAnnotation("XLabel", xtitle);
- // prof->setAnnotation("YLabel", ytitle);
- // prof->setAnnotation("ZLabel", ztitle);
- // return prof;
- // }
-
-
- // Profile2DPtr Analysis::bookProfile2D(const string& hname,
- // const string& title,
- // const string& xtitle,
- // const string& ytitle,
- // const string& ztitle) {
- // const Scatter3D& refdata = refData(hname);
- // return bookProfile2D(hname, refdata, title, xtitle, ytitle, ztitle);
- // }
-
-
- // Profile2DPtr Analysis::bookProfile2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
- // const string& title,
- // const string& xtitle,
- // const string& ytitle,
- // const string& ztitle) {
- // const string axisCode = makeAxisCode(datasetId, xAxisId, yAxisId);
- // return bookProfile2D(axisCode, title, xtitle, ytitle, ztitle);
- // }
-
-
- /////////////////
-
-
Scatter2DPtr Analysis::bookScatter2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
bool copy_pts,
const string& title,
const string& xtitle,
const string& ytitle) {
const string axisCode = makeAxisCode(datasetId, xAxisId, yAxisId);
return bookScatter2D(axisCode, copy_pts, title, xtitle, ytitle);
}
Scatter2DPtr Analysis::bookScatter2D(const string& hname,
bool copy_pts,
const string& title,
const string& xtitle,
const string& ytitle) {
- Scatter2DPtr s;
+ shared_ptr<Scatter2DPtr> sp = make_shared<Scatter2DPtr>(make_shared<Scatter2D>());
const string path = histoPath(hname);
if (copy_pts) {
const Scatter2D& refdata = refData(hname);
- s = make_shared<Scatter2D>(refdata, path);
- foreach (Point2D& p, s->points()) p.setY(0, 0);
+ *sp = Scatter2DPtr(make_shared<Scatter2D>(refdata, path));
+ foreach (Point2D& p, (*sp)->points()) p.setY(0, 0);
} else {
- s = make_shared<Scatter2D>(path);
+ *sp = Scatter2DPtr(make_shared<Scatter2D>(path));
}
- addAnalysisObject(s);
+ addAnalysisObject(sp);
+
+ Scatter2DPtr scat = *sp;
MSG_TRACE("Made scatter " << hname << " for " << name());
- s->setTitle(title);
- s->setAnnotation("XLabel", xtitle);
- s->setAnnotation("YLabel", ytitle);
- return s;
+ scat->setTitle(title);
+ scat->setAnnotation("XLabel", xtitle);
+ scat->setAnnotation("YLabel", ytitle);
+ return scat;
}
Scatter2DPtr Analysis::bookScatter2D(const string& hname,
size_t npts, double lower, double upper,
const string& title,
const string& xtitle,
const string& ytitle) {
const string path = histoPath(hname);
- Scatter2DPtr s = make_shared<Scatter2D>(path);
+ shared_ptr<Scatter2DPtr> sp = make_shared<Scatter2DPtr>(make_shared<Scatter2D>(path));
const double binwidth = (upper-lower)/npts;
for (size_t pt = 0; pt < npts; ++pt) {
const double bincentre = lower + (pt + 0.5) * binwidth;
- s->addPoint(bincentre, 0, binwidth/2.0, 0);
+ (*sp)->addPoint(bincentre, 0, binwidth/2.0, 0);
}
- addAnalysisObject(s);
+ addAnalysisObject(sp);
+
+ Scatter2DPtr scat = *sp;
+
MSG_TRACE("Made scatter " << hname << " for " << name());
- s->setTitle(title);
- s->setAnnotation("XLabel", xtitle);
- s->setAnnotation("YLabel", ytitle);
- return s;
+ scat->setTitle(title);
+ scat->setAnnotation("XLabel", xtitle);
+ scat->setAnnotation("YLabel", ytitle);
+ return scat;
}
Scatter2DPtr Analysis::bookScatter2D(const string& hname,
const vector<double>& binedges,
const string& title,
const string& xtitle,
const string& ytitle) {
const string path = histoPath(hname);
- Scatter2DPtr s = make_shared<Scatter2D>(path);
+ shared_ptr<Scatter2DPtr> sp =
+ make_shared<Scatter2DPtr>(make_shared<Scatter2D>(path));
for (size_t pt = 0; pt < binedges.size()-1; ++pt) {
const double bincentre = (binedges[pt] + binedges[pt+1]) / 2.0;
const double binwidth = binedges[pt+1] - binedges[pt];
- s->addPoint(bincentre, 0, binwidth/2.0, 0);
+ (*sp)->addPoint(bincentre, 0, binwidth/2.0, 0);
}
- addAnalysisObject(s);
+ addAnalysisObject(sp);
+
+ Scatter2DPtr scat = *sp;
MSG_TRACE("Made scatter " << hname << " for " << name());
- s->setTitle(title);
- s->setAnnotation("XLabel", xtitle);
- s->setAnnotation("YLabel", ytitle);
- return s;
+ scat->setTitle(title);
+ scat->setAnnotation("XLabel", xtitle);
+ scat->setAnnotation("YLabel", ytitle);
+ return scat;
}
+ */
/////////////////////
void Analysis::divide(CounterPtr c1, CounterPtr c2, Scatter1DPtr s) const {
const string path = s->path();
*s = *c1 / *c2;
s->setPath(path);
}
void Analysis::divide(const Counter& c1, const Counter& c2, Scatter1DPtr s) const {
const string path = s->path();
*s = c1 / c2;
s->setPath(path);
}
void Analysis::divide(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const {
const string path = s->path();
*s = *h1 / *h2;
s->setPath(path);
}
void Analysis::divide(const Histo1D& h1, const Histo1D& h2, Scatter2DPtr s) const {
const string path = s->path();
*s = h1 / h2;
s->setPath(path);
}
void Analysis::divide(Profile1DPtr p1, Profile1DPtr p2, Scatter2DPtr s) const {
const string path = s->path();
*s = *p1 / *p2;
s->setPath(path);
}
void Analysis::divide(const Profile1D& p1, const Profile1D& p2, Scatter2DPtr s) const {
const string path = s->path();
*s = p1 / p2;
s->setPath(path);
}
void Analysis::divide(Histo2DPtr h1, Histo2DPtr h2, Scatter3DPtr s) const {
const string path = s->path();
*s = *h1 / *h2;
s->setPath(path);
}
void Analysis::divide(const Histo2D& h1, const Histo2D& h2, Scatter3DPtr s) const {
const string path = s->path();
*s = h1 / h2;
s->setPath(path);
}
void Analysis::divide(Profile2DPtr p1, Profile2DPtr p2, Scatter3DPtr s) const {
const string path = s->path();
*s = *p1 / *p2;
s->setPath(path);
}
void Analysis::divide(const Profile2D& p1, const Profile2D& p2, Scatter3DPtr s) const {
const string path = s->path();
*s = p1 / p2;
s->setPath(path);
}
/// @todo Counter and Histo2D efficiencies and asymms
void Analysis::efficiency(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const {
const string path = s->path();
*s = YODA::efficiency(*h1, *h2);
s->setPath(path);
}
void Analysis::efficiency(const Histo1D& h1, const Histo1D& h2, Scatter2DPtr s) const {
const string path = s->path();
*s = YODA::efficiency(h1, h2);
s->setPath(path);
}
void Analysis::asymm(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const {
const string path = s->path();
*s = YODA::asymm(*h1, *h2);
s->setPath(path);
}
void Analysis::asymm(const Histo1D& h1, const Histo1D& h2, Scatter2DPtr s) const {
const string path = s->path();
*s = YODA::asymm(h1, h2);
s->setPath(path);
}
void Analysis::scale(CounterPtr cnt, double factor) {
if (!cnt) {
MSG_WARNING("Failed to scale counter=NULL in analysis " << name() << " (scale=" << factor << ")");
return;
}
if (std::isnan(factor) || std::isinf(factor)) {
MSG_WARNING("Failed to scale counter=" << cnt->path() << " in analysis: " << name() << " (invalid scale factor = " << factor << ")");
factor = 0;
}
MSG_TRACE("Scaling counter " << cnt->path() << " by factor " << factor);
try {
cnt->scaleW(factor);
} catch (YODA::Exception& we) {
MSG_WARNING("Could not scale counter " << cnt->path());
return;
}
}
void Analysis::normalize(Histo1DPtr histo, double norm, bool includeoverflows) {
if (!histo) {
MSG_WARNING("Failed to normalize histo=NULL in analysis " << name() << " (norm=" << norm << ")");
return;
}
MSG_TRACE("Normalizing histo " << histo->path() << " to " << norm);
try {
histo->normalize(norm, includeoverflows);
} catch (YODA::Exception& we) {
MSG_WARNING("Could not normalize histo " << histo->path());
return;
}
}
void Analysis::scale(Histo1DPtr histo, double factor) {
if (!histo) {
MSG_WARNING("Failed to scale histo=NULL in analysis " << name() << " (scale=" << factor << ")");
return;
}
if (std::isnan(factor) || std::isinf(factor)) {
MSG_WARNING("Failed to scale histo=" << histo->path() << " in analysis: " << name() << " (invalid scale factor = " << factor << ")");
factor = 0;
}
MSG_TRACE("Scaling histo " << histo->path() << " by factor " << factor);
try {
histo->scaleW(factor);
} catch (YODA::Exception& we) {
MSG_WARNING("Could not scale histo " << histo->path());
return;
}
}
void Analysis::normalize(Histo2DPtr histo, double norm, bool includeoverflows) {
if (!histo) {
MSG_ERROR("Failed to normalize histo=NULL in analysis " << name() << " (norm=" << norm << ")");
return;
}
MSG_TRACE("Normalizing histo " << histo->path() << " to " << norm);
try {
histo->normalize(norm, includeoverflows);
} catch (YODA::Exception& we) {
MSG_WARNING("Could not normalize histo " << histo->path());
return;
}
}
void Analysis::scale(Histo2DPtr histo, double factor) {
if (!histo) {
MSG_ERROR("Failed to scale histo=NULL in analysis " << name() << " (scale=" << factor << ")");
return;
}
if (std::isnan(factor) || std::isinf(factor)) {
MSG_ERROR("Failed to scale histo=" << histo->path() << " in analysis: " << name() << " (invalid scale factor = " << factor << ")");
factor = 0;
}
MSG_TRACE("Scaling histo " << histo->path() << " by factor " << factor);
try {
histo->scaleW(factor);
} catch (YODA::Exception& we) {
MSG_WARNING("Could not scale histo " << histo->path());
return;
}
}
void Analysis::integrate(Histo1DPtr h, Scatter2DPtr s) const {
// preserve the path info
const string path = s->path();
*s = toIntegralHisto(*h);
s->setPath(path);
}
void Analysis::integrate(const Histo1D& h, Scatter2DPtr s) const {
// preserve the path info
const string path = s->path();
*s = toIntegralHisto(h);
s->setPath(path);
}
/// @todo 2D versions of integrate... defined how, exactly?!?
//////////////////////////////////
- void Analysis::addAnalysisObject(AnalysisObjectPtr ao) {
+ void Analysis::addAnalysisObject(shared_ptr<MultiweightAOPtr> ao) {
_analysisobjects.push_back(ao);
}
void Analysis::removeAnalysisObject(const string& path) {
- for (vector<AnalysisObjectPtr>::iterator it = _analysisobjects.begin(); it != _analysisobjects.end(); ++it) {
- if ((*it)->path() == path) {
+ for (vector<shared_ptr<MultiweightAOPtr> >::iterator it = _analysisobjects.begin(); it != _analysisobjects.end(); ++it) {
+ if ((**it)->path() == path) {
_analysisobjects.erase(it);
break;
}
}
}
- void Analysis::removeAnalysisObject(AnalysisObjectPtr ao) {
- for (vector<AnalysisObjectPtr>::iterator it = _analysisobjects.begin(); it != _analysisobjects.end(); ++it) {
- if (*it == ao) {
+ void Analysis::removeAnalysisObject(const MultiweightAOPtr& ao) {
+ for (vector<shared_ptr<MultiweightAOPtr> >::iterator it = _analysisobjects.begin(); it != _analysisobjects.end(); ++it) {
+ if (**it == ao) {
_analysisobjects.erase(it);
break;
}
}
}
}
diff --git a/src/Core/AnalysisHandler.cc b/src/Core/AnalysisHandler.cc
--- a/src/Core/AnalysisHandler.cc
+++ b/src/Core/AnalysisHandler.cc
@@ -1,402 +1,422 @@
// -*- C++ -*-
#include "Rivet/Config/RivetCommon.hh"
#include "Rivet/AnalysisHandler.hh"
#include "Rivet/Analysis.hh"
#include "Rivet/Tools/ParticleName.hh"
#include "Rivet/Tools/BeamConstraint.hh"
#include "Rivet/Tools/Logging.hh"
#include "Rivet/Projections/Beam.hh"
#include "YODA/WriterYODA.h"
-#include "Rivet/Histo1D.hh"
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;
_eventNumber = ge.event_number();
setWeightNames(ge);
if (haveNamedWeights()) {
MSG_INFO("Using named weights");
}
else {
_weightNames.clear();
_weightNames.push_back("0");
MSG_INFO("NOT using named weights. Using first weight as nominal weight");
}
_numWeightTypes = _weightNames.size();
_sumOfWeights = vector<double>(_numWeightTypes, 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;
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());
}
}
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
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
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::setWeightNames(const GenEvent& ge) {
+ return;
+
+ /// @todo
+ /// this relies on boost and is removed for now
+
+ /*
/// reroute the print output to a stringstream and process
/// The iteration is done over a map in hepmc2 so this is safe
ostringstream stream;
ge.weights().print(stream); // Super lame, I know
string str = stream.str();
vector<string> temp;
- boost::split(temp, str, boost::is_any_of("="));
+ split(temp, str, is_any_of("="));
vector<string> tokens;
- boost::split(tokens, temp[0], boost::is_any_of(" "));
+ split(tokens, temp[0], is_any_of(" "));
/// Weights come in tuples: (StringName, doubleWeight) (AnotherName, anotherWeight) etc
string wname;
for (unsigned int i=0; i<tokens.size()-1;++i) {
temp.clear();
- boost::split(temp, tokens[i], boost::is_any_of(","));
+ split(temp, tokens[i], is_any_of(","));
if (temp.size()>0) {
wname = temp[0];
- trim_left_if(wname,is_any_of("("));
+ trim_left_if(wname, is_any_of("("));
_weightNames.push_back(wname);
/// Turn into debug message
MSG_INFO("Name of weight #" << i << ": " << wname);
}
}
+ */
}
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);
// won't happen for first event because _eventNumber is set in
// init()
if (_eventNumber != ge.event_number()) {
+ /// @todo
+ /// can we get away with not passing a matrix?
+
// if this is indeed a new event, push the temporary
// histograms and reset
- foreach (Rivet::Histo1DPtr& histo, _histograms) {
- histo.rivetHisto1D()->pushToPersistent(_subEventWeights);
+ for (const AnaHandle& a : _analyses) {
+ for (const shared_ptr<MultiweightAOPtr>& aoptr : a->analysisObjects()) {
+ aoptr->pushToPersistent(_subEventWeights);
+ }
}
_eventNumber = ge.event_number();
// @todo
// is this correct?
// only incremented when event number changes.
_numEvents++;
MSG_DEBUG("Event #" << _numEvents << " nominal weight sum = " << _sumOfWeights[0]);
MSG_DEBUG("Event has " << _subEventWeights.size() << " sub events.");
_subEventWeights.clear();
}
- foreach (Rivet::Histo1DPtr& histo, _histograms) {
- histo.rivetHisto1D()->newEvent();
+ for (const AnaHandle& a : _analyses) {
+ for (const shared_ptr<MultiweightAOPtr>& aoptr : a->analysisObjects()) {
+ aoptr->newSubEvent();
+ }
}
_subEventWeights.push_back(event.weights());
MSG_DEBUG("Analyzing subevent #" << _subEventWeights.size() - 1 << ".");
// @todo
// is this the correct way to sum weights over sub events?
for (unsigned int iWeight = 0; iWeight < _sumOfWeights.size(); iWeight++) {
_sumOfWeights[iWeight] += _subEventWeights.back()[iWeight];
}
// 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
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");
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.
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());
// for (const AnaHandle& a : _analyses) MSG_WARNING(a->name());
return *this;
}
AnalysisHandler& AnalysisHandler::removeAnalysis(const string& analysisname) {
std::shared_ptr<Analysis> toremove;
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;
+ vector<shared_ptr<Rivet::AnalysisObjectPtr> > AnalysisHandler::getRivetAOs() const {
+ vector<shared_ptr<AnalysisObjectPtr> > rtn;
+ for (size_t i = 0; i < _numWeights; i++) {
+ // HERE WE ARE
+ }
+ }
+
+ vector<YODA::AnalysisObjectPtr> AnalysisHandler::getData() const {
+ vector<YODA::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") );
for (const AnaHandle a : analyses()) {
vector<AnalysisObjectPtr> aos = a->analysisObjects();
// MSG_WARNING(a->name() << " " << aos.size());
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(),
[](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
+ } catch ( YODA::WriteError ) {
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;
for (AnaHandle a : _analyses) {
rtn.push_back(a->name());
}
return rtn;
}
const AnaHandle AnalysisHandler::analysis(const std::string& analysisname) const {
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) {
for (const string& aname : analysisnames) {
//MSG_DEBUG("Adding analysis '" << aname << "'");
addAnalysis(aname);
}
return *this;
}
AnalysisHandler& AnalysisHandler::removeAnalyses(const std::vector<std::string>& analysisnames) {
for (const string& aname : analysisnames) {
removeAnalysis(aname);
}
return *this;
}
bool AnalysisHandler::needCrossSection() const {
bool rtn = false;
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,114 +1,119 @@
// -*- C++ -*-
#include "Rivet/Event.hh"
#include "Rivet/Tools/BeamConstraint.hh"
#include "Rivet/Projections/Beam.hh"
#include "HepMC/GenEvent.h"
namespace Rivet {
ParticlePair Event::beams() const { return Rivet::beams(*this); }
// PdgIdPair Event::beamIds() const { return pids(beams()); }
double Event::sqrtS() const { return Rivet::sqrtS(beams()); }
double Event::asqrtS() const { return Rivet::asqrtS(beams()); }
// 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);
// }
// }
const Particles& Event::allParticles() const {
if (_particles.empty()) { //< assume that empty means no attempt yet made
for (const GenParticle* gp : particles(genEvent())) {
_particles += Particle(gp);
}
}
return _particles;
}
- double Event::weight() const {
- return (!_genevent.weights().empty()) ? _genevent.weights()[0] : 1.0;
+ vector<double> Event::weights() const {
+ vector<double> v;
+
+ for (unsigned int iw = 0; iw < _genevent.weights().size(); iw++)
+ v.push_back(_genevent.weights()[iw]);
+
+ return v;
}
}
diff --git a/src/Tools/RivetYODA.cc b/src/Tools/RivetYODA.cc
--- a/src/Tools/RivetYODA.cc
+++ b/src/Tools/RivetYODA.cc
@@ -1,115 +1,115 @@
#include "Rivet/Config/RivetCommon.hh"
#include "Rivet/Tools/RivetYODA.hh"
#include "Rivet/Tools/RivetPaths.hh"
#include "YODA/ReaderYODA.h"
#include "YODA/ReaderAIDA.h"
using namespace std;
namespace Rivet {
string getDatafilePath(const string& papername) {
/// Try to find YODA otherwise fall back to try AIDA
const string path1 = findAnalysisRefFile(papername + ".yoda");
if (!path1.empty()) return path1;
const string path2 = findAnalysisRefFile(papername + ".aida");
if (!path2.empty()) return path2;
throw Rivet::Error("Couldn't find ref data file '" + papername + ".yoda/aida" +
" in $RIVET_REF_PATH, '" + getRivetDataPath() + "', or '.'");
}
- map<string, AnalysisObjectPtr> getRefData(const string& papername) {
+ map<string, YODA::AnalysisObjectPtr> getRefData(const string& papername) {
const string datafile = getDatafilePath(papername);
// Make an appropriate data file reader and read the data objects
/// @todo Remove AIDA support some day...
YODA::Reader& reader = (datafile.find(".yoda") != string::npos) ? \
YODA::ReaderYODA::create() : YODA::ReaderAIDA::create();
vector<YODA::AnalysisObject *> aovec;
reader.read(datafile, aovec);
// Return value, to be populated
- map<string, AnalysisObjectPtr> rtn;
+ map<string, YODA::AnalysisObjectPtr> rtn;
foreach ( YODA::AnalysisObject* ao, aovec ) {
- AnalysisObjectPtr refdata(ao);
+ YODA::AnalysisObjectPtr refdata(ao);
if (!refdata) continue;
const string plotpath = refdata->path();
// Split path at "/" and only return the last field, i.e. the histogram ID
const size_t slashpos = plotpath.rfind("/");
const string plotname = (slashpos+1 < plotpath.size()) ? plotpath.substr(slashpos+1) : "";
rtn[plotname] = refdata;
}
return rtn;
}
void Histo1DPtr::pushToPersistent(const vector<vector<double> >& weight) {
// loop over event weights
for (size_t m = 0; m < _persistent.size(); ++m) {
// this is the initial event---it always exists
- YODA::Histo1DPtr sum = boost::make_shared<YODA::Histo1D>(_evgroup[0]->clone());
+ YODA::Histo1DPtr sum = make_shared<YODA::Histo1D>(_evgroup[0]->clone());
sum->scaleW( weight[0][m] );
// loop over additional subevents (maybe there aren't any)
// note loop starts at 1
for (size_t n = 1; n < _evgroup.size(); ++n) {
- YODA::Histo1DPtr tmp = boost::make_shared<YODA::Histo1D>(_evgroup[n]->clone());
+ YODA::Histo1DPtr tmp = make_shared<YODA::Histo1D>(_evgroup[n]->clone());
tmp->scaleW( weight[n][m] );
*sum += *tmp;
}
// sum typically only has one bin filled
// do we really need to fill bin.size() times?
bool filled_something = false;
/// @todo
/// this is not correct!
/// we need to trick the histogram into thinking it was
/// filled exactly onen time, even if subevents fall into
/// different bins of the histogram
/// how
foreach(const YODA::HistoBin1D& b, sum->bins()) {
if ( b.effNumEntries() != 0 && b.sumW() != 0 ) {
// @todo number of fills will be wrong
// needs to be xMid. xMean is not valid if bin contains neg.weights
_persistent[m]->fill(b.xMid(), b.sumW());
filled_something = true;
}
}
// @todo
// overflow
}
_evgroup.clear();
_active.reset();
}
void Histo2DPtr::pushToPersistent(const vector<vector<double> >& weight) {
/// @todo
return;
}
void Profile1DPtr::pushToPersistent(const vector<vector<double> >& weight) {
/// @todo
return;
}
void Profile2DPtr::pushToPersistent(const vector<vector<double> >& weight) {
/// @todo
return;
}
void Counter::pushToPersistent(const vector<vector<double> >& weight) {
/// @todo
return;
}
}