diff --git a/include/Rivet/AnalysisHandler.hh b/include/Rivet/AnalysisHandler.hh
--- a/include/Rivet/AnalysisHandler.hh
+++ b/include/Rivet/AnalysisHandler.hh
@@ -1,242 +1,250 @@
// -*- C++ -*-
#ifndef RIVET_RivetHandler_HH
#define RIVET_RivetHandler_HH
#include "Rivet/Config/RivetCommon.hh"
#include "Rivet/Particle.hh"
#include "Rivet/AnalysisLoader.hh"
#include "Rivet/Tools/RivetYODA.hh"
namespace Rivet {
// Forward declaration and smart pointer for Analysis
class Analysis;
typedef std::shared_ptr<Analysis> AnaHandle;
// Needed to make smart pointers compare equivalent in the STL set
struct CmpAnaHandle {
bool operator() (const AnaHandle& a, const AnaHandle& b) {
return a.get() < b.get();
}
};
/// A class which handles a number of analysis objects to be applied to
/// generated events. An {@link Analysis}' AnalysisHandler is also responsible
/// for handling the final writing-out of histograms.
class AnalysisHandler {
public:
/// @name Constructors and destructors. */
//@{
/// Preferred constructor, with optional run name.
AnalysisHandler(const string& runname="");
/// @brief Destructor
/// The destructor is not virtual, as this class should not be inherited from.
~AnalysisHandler();
//@}
private:
/// Get a logger object.
Log& getLog() const;
public:
/// @name Run properties
//@{
/// Get the name of this run.
string runName() const;
/// Get the number of events seen. Should only really be used by external
/// steering code or analyses in the finalize phase.
size_t numEvents() const;
/// Get the sum of the event weights seen - the weighted equivalent of the
/// number of events. Should only really be used by external steering code
/// or analyses in the finalize phase.
double sumOfWeights() const;
/// Set sum of weights. This is useful if Rivet is steered externally and
/// the analyses are run for a sub-contribution of the events
/// (but of course have to be normalised to the total sum of weights)
void setSumOfWeights(const double& sum);
/// Is cross-section information required by at least one child analysis?
bool needCrossSection() const;
/// Set the cross-section for the process being generated.
AnalysisHandler& setCrossSection(double xs);
/// Get the cross-section known to the handler.
double crossSection() const {
return _xs;
}
/// Whether the handler knows about a cross-section.
bool hasCrossSection() const;
/// Set 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();
+ /// Method reads YODA files and adds the contained AnalysisObjects
+ /// to the member map "_readObjects".
+ void readData(const std::string& filename);
+
//@}
/// @name Histogram / data object access
//@{
/// Get all analyses' plots as a vector of analysis objects.
std::vector<AnalysisObjectPtr> getData() const;
/// Write all analyses' plots to the named file.
void writeData(const std::string& filename) const;
//@}
private:
+ /// Map containing read YODA objects.
+ std::map< std::string, AnalysisObjectPtr > _readObjects;
+ bool _haveReadData;
+
/// The collection of Analysis objects to be used.
set<AnaHandle, CmpAnaHandle> _analyses;
/// @name Run properties
//@{
/// Run name
std::string _runname;
/// Number of events seen.
/// @todo Replace by a counter
unsigned int _numEvents;
/// Sum of event weights seen.
/// @todo Replace by a counter
double _sumOfWeights, _sumOfWeightsSq;
/// Cross-section known to AH
double _xs, _xserr;
/// Beams used by this run.
ParticlePair _beams;
/// Flag to check if init has been called
bool _initialised;
/// Flag whether input event beams should be ignored in compatibility check
bool _ignoreBeams;
//@}
private:
/// The assignment operator is private and must never be called.
/// In fact, it should not even be implemented.
AnalysisHandler& operator=(const AnalysisHandler&);
/// The copy constructor is private and must never be called. In
/// fact, it should not even be implemented.
AnalysisHandler(const AnalysisHandler&);
};
}
#endif
diff --git a/src/Core/AnalysisHandler.cc b/src/Core/AnalysisHandler.cc
--- a/src/Core/AnalysisHandler.cc
+++ b/src/Core/AnalysisHandler.cc
@@ -1,333 +1,393 @@
// -*- 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 "YODA/ReaderYODA.h"
namespace Rivet {
- AnalysisHandler::AnalysisHandler(const string& runname)
- : _runname(runname), _numEvents(0),
+ AnalysisHandler::AnalysisHandler(const string& runname) // xyz Initialize introduced member "_haveReadData"
+ : _haveReadData(false), _runname(runname), _numEvents(0),
_sumOfWeights(0.0), _xs(NAN),
_initialised(false), _ignoreBeams(false)
{}
AnalysisHandler::~AnalysisHandler()
{}
Log& AnalysisHandler::getLog() const {
return Log::getLog("Rivet.Analysis.Handler");
}
void AnalysisHandler::init(const GenEvent& ge) {
if (_initialised)
throw UserError("AnalysisHandler::init has already been called: cannot re-initialize!");
setRunBeams(Rivet::beams(ge));
MSG_DEBUG("Initialising the analysis handler");
_numEvents = 0;
_sumOfWeights = 0.0;
_sumOfWeightsSq = 0.0;
// Check that analyses are beam-compatible, and remove those that aren't
const size_t num_anas_requested = analysisNames().size();
vector<string> anamestodelete;
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::analyze(const GenEvent& ge) {
// Call init with event as template if not already initialised
if (!_initialised) init(ge);
assert(_initialised);
// Ensure that beam details match those from the first event
const PdgIdPair beams = Rivet::beamIds(ge);
const double sqrts = Rivet::sqrtS(ge);
if (!compatible(beams, _beams) || !fuzzyEquals(sqrts, sqrtS())) {
cerr << "Event beams mismatch: "
<< PID::toBeamsString(beams) << " @ " << sqrts/GeV << " GeV" << " vs. first beams "
<< this->beams() << " @ " << this->sqrtS()/GeV << " GeV" << endl;
exit(1);
}
// Create the Rivet event wrapper
/// @todo Filter/normalize the event here
Event event(ge);
// Weights
/// @todo Drop this / just report first weight when we support multiweight events
_numEvents += 1;
_sumOfWeights += event.weight();
_sumOfWeightsSq += sqr(event.weight());
MSG_DEBUG("Event #" << _numEvents << " weight = " << event.weight());
// Cross-section
#ifdef HEPMC_HAS_CROSS_SECTION
if (ge.cross_section()) {
_xs = ge.cross_section()->cross_section();
_xserr = ge.cross_section()->cross_section_error();
}
#endif
// Run the analyses
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;
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
throw UserError("Unexpected error in writing file to: " + filename);
}
}
+ /// xyz Reads Objects from a given YODA file and stores those in a map. No duplicates possible
+ /// since they would be overwritten. @note Method has to be called before AnalysisHandler::init
+ /// since there potential read objects are given to the added analyses.
+ void AnalysisHandler::readData( const std::string& filename ) {
+
+ YODA::Reader &reader = YODA::mkReader( filename );
+
+ std::vector< YODA::AnalysisObject* > aos;
+
+ try {
+ reader.read( filename, aos );
+ } catch(const YODA::ReadError& err) {
+ MSG_ERROR( "Error in AnalysisHandler::readData: " << err.what() );
+ exit(1);
+ }
+ MSG_INFO( "Read data from " << filename );
+
+ for( const auto ao : aos ) {
+ //MSG_INFO( "type of ao is " << typeid( ao ).name() );
+ AnalysisObjectPtr append( ao );
+
+ int numEntries = 0;
+
+ try {
+
+ if( Histo1DPtr tmpHisto1D = dynamic_pointer_cast< YODA::Histo1D >( append ) ) {
+ numEntries = tmpHisto1D->numEntries();
+ } else if ( Histo2DPtr tmpHisto2D = dynamic_pointer_cast< YODA::Histo2D >( append ) ) {
+ numEntries = tmpHisto2D->numEntries();
+ } else if ( Profile1DPtr tmpProfile1D = dynamic_pointer_cast< YODA::Profile1D >( append ) ) {
+ numEntries = tmpProfile1D->numEntries();
+ } else if ( Profile2DPtr tmpProfile2D = dynamic_pointer_cast< YODA::Profile2D >( append ) ) {
+ numEntries = tmpProfile2D->numEntries();
+ } else if ( Scatter1DPtr tmpScatter1D = dynamic_pointer_cast< YODA::Scatter1D >( append ) ) {
+ numEntries = tmpScatter1D->numPoints();
+ } else if ( Scatter2DPtr tmpScatter2D = dynamic_pointer_cast< YODA::Scatter2D >( append ) ) {
+ numEntries = tmpScatter2D->numPoints();
+ } else if ( Scatter3DPtr tmpScatter3D = dynamic_pointer_cast< YODA::Scatter3D >( append ) ) {
+ numEntries = tmpScatter3D->numPoints();
+ }
+
+ if( numEntries <= 0 ) {
+ continue;
+ }
+
+ if ( !_readObjects[ao->path()] ) {
+ _readObjects[ao->path()] = append;
+ }
+ // in case an object with the given path is already read try to merge the new and the old one and save this version
+ // else
+ // merge...; or fail
+ } catch(...) {
+ MSG_ERROR( "Unable to read data!" );
+ throw UserError( "Unexpected error while reading data!" );
+ }
+ }
+ _haveReadData = true;
+ }
+
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;
}
}