No OneTemporary
Actions

Size

62 KB

Subscribers

None

View Options

	diff --git a/bin/rivet-nopy.cc b/bin/rivet-nopy.cc
	--- a/bin/rivet-nopy.cc
	+++ b/bin/rivet-nopy.cc
	@@ -1,41 +1,42 @@
	#include "Rivet/AnalysisHandler.hh"
	#include "Rivet/AnalysisLoader.hh"
	//#include "HepMC/IO_GenEvent.h"
	#include "Rivet/Tools/RivetHepMC.hh"
	+#include <fstream>

	using namespace std;

	int main(int argc, char** argv) {
	if (argc < 2) {
	cerr << "Usage: " << argv[0] << " <hepmcfile> <ana1> [<ana2> ...]" << '\n';
	cout << "Available analyses:\n";
	for (const string& a : Rivet::AnalysisLoader::analysisNames())
	cout << " " << a << "\n";
	cout << endl;
	return 1;
	}

	Rivet::AnalysisHandler ah;
	for (int i = 2; i < argc; ++i) {
	ah.addAnalysis(argv[i]);
	}

	std::ifstream istr(argv[1], std::ios::in);

	std::shared_ptr<Rivet::HepMC_IO_type> reader = Rivet::HepMCUtils::makeReader(istr);

	std::shared_ptr<Rivet::RivetHepMC::GenEvent> evt = make_shared<Rivet::RivetHepMC::GenEvent>();

	while(reader && Rivet::HepMCUtils::readEvent(reader, evt)){
	ah.analyze(evt.get());
	evt.reset(new Rivet::RivetHepMC::GenEvent());
	}

	istr.close();

	- ah.setCrossSection(1.0, 0.0);
	+ ah.setCrossSection(make_pair(1.0, 0.0));
	ah.finalize();
	ah.writeData("Rivet.yoda");

	return 0;
	}
	diff --git a/include/Rivet/AnalysisHandler.hh b/include/Rivet/AnalysisHandler.hh
	--- a/include/Rivet/AnalysisHandler.hh
	+++ b/include/Rivet/AnalysisHandler.hh
	@@ -1,366 +1,370 @@
	// -- 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;


	/// 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 and weights
	//@{

	/// 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;

	/// @brief Access the sum of the event weights seen
	///
	/// This is the weighted equivalent of the number of events. It should only
	/// be used by external steering code or analyses in the finalize phase.
	double sumW() const { return _eventCounter->sumW(); }
	/// Access to the sum of squared-weights
	double sumW2() const { return _eventCounter->sumW2(); }

	/// Names of event weight categories
	const vector<string>& weightNames() const { return _weightNames; }

	/// Are any of the weights non-numeric?
	size_t numWeights() const { return _weightNames.size(); }

	/// Are any of the weights non-numeric?
	bool haveNamedWeights() const;

	/// Set the weight names from a GenEvent
	void setWeightNames(const GenEvent& ge);

	/// Get the index of the nominal weight-stream
	size_t defaultWeightIndex() const { return _defaultWeightIdx; }

	//@}


	/// @name Cross-sections
	//@{

	/// Get the cross-section known to the handler.
	Scatter1DPtr crossSection() const { return _xs; }

	/// Set the cross-section for the process being generated.
	- AnalysisHandler& setCrossSection(double xs, double xserr);
	+ void setCrossSection(pair<double, double> xsec);

	/// Get the nominal cross-section
	double nominalCrossSection() const {
	_xs.get()->setActiveWeightIdx(_defaultWeightIdx);
	const YODA::Scatter1D::Points& ps = _xs->points();
	if (ps.size() != 1) {
	string errMsg = "cross section missing when requesting nominal cross section";
	throw Error(errMsg);
	}
	double xs = ps[0].x();
	_xs.get()->unsetActiveWeight();
	return xs;
	}

	//@}


	/// @name Beams
	//@{

	/// 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);

	/// Setter for _skipWeights
	void skipMultiWeights(bool ignore=false);

	//@}


	/// @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::map<std::string, AnaHandle>& analysesMap() const {
	return _analyses;
	}

	/// Get the collection of currently registered analyses.
	std::vector<AnaHandle> analyses() const {
	std::vector<AnaHandle> rtn;
	rtn.reserve(_analyses.size());
	for (const auto& apair : _analyses) rtn.push_back(apair.second);
	return rtn;
	}

	/// Get a registered analysis by name.
	AnaHandle analysis(const std::string& analysisname) {
	if ( _analyses.find(analysisname) == _analyses.end() )
	throw LookupError("No analysis named '" + analysisname + "' registered in AnalysisHandler");
	try {
	return _analyses[analysisname];
	} catch (...) {
	throw LookupError("No analysis named '" + analysisname + "' registered in AnalysisHandler");
	}
	}

	/// 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 an analysis with a map of analysis options.
	AnalysisHandler& addAnalysis(const std::string& analysisname, std::map<string, string> pars);

	/// @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
	//@{

	+ /// After all subevents in an event group has been processed push
	+ /// all histo fills to the relevant histograms.
	+ void pushToPersistent();
	+
	/// Read analysis plots into the histo collection (via addData) from the named file.
	void readData(const std::string& filename);

	/// Get all multi-weight Rivet analysis object wrappers
	vector<MultiweightAOPtr> getRivetAOs() const;

	/// Get a pointer to a preloaded yoda object with the given path,
	/// or null if path is not found.
	const YODA::AnalysisObjectPtr getPreload(string path) const {
	auto it = _preloads.find(path);
	if ( it == _preloads.end() ) return nullptr;
	return it->second;
	}

	/// Write all analyses' plots (via getData) to the named file.
	void writeData(const std::string& filename) const;

	/// Tell Rivet to dump intermediate result to a file named @a
	/// dumpfile every @a period'th event. If @period is not positive,
	/// no dumping will be done.
	void dump(string dumpfile, int period) {
	_dumpPeriod = period;
	_dumpFile = dumpfile;
	}

	/// Take the vector of yoda files and merge them together using
	/// the cross section and weight information provided in each
	/// file. Each file in @a aofiles is assumed to have been produced
	/// by Rivet. By default the files are assumed to contain
	/// different processes (or the same processs but mutually
	/// exclusive cuts), but if @a equiv if ture, the files are
	/// assumed to contain output of completely equivalent (but
	/// statistically independent) Rivet runs. The corresponding
	/// analyses will be loaded and their analysis objects will be
	/// filled with the merged result. finalize() will be run on each
	/// relevant analysis. The resulting YODA file can then be rwitten
	/// out by writeData(). If delopts is non-empty, it is assumed to
	/// contain names different options to be merged into the same
	/// analysis objects.
	void mergeYodas(const vector<string> & aofiles,
	const vector<string> & delopts = vector<string>(),
	bool equiv = false);

	/// Helper function to strip specific options from data object paths.
	void stripOptions(YODA::AnalysisObjectPtr ao,
	const vector<string> & delopts) const;

	//@}


	/// Indicate which Rivet stage we're in.
	/// At the moment, only INIT is used to enable booking.
	enum class Stage { OTHER, INIT, FINALIZE };

	/// Which stage are we in?
	Stage stage() const { return _stage; }


	private:

	/// Current handler stage
	Stage _stage = Stage::OTHER;

	/// The collection of Analysis objects to be used.
	std::map<std::string, AnaHandle> _analyses;

	/// A vector of pre-loaded object which do not have a valid
	/// Analysis plugged in.
	map<string,YODA::AnalysisObjectPtr> _preloads;

	/// A vector containing copies of analysis objects after
	/// finalize() has been run.
	vector<YODA::AnalysisObjectPtr> _finalizedAOs;


	/// @name Run properties
	//@{

	/// Weight names
	std::vector<std::string> _weightNames;
	std::vector<std::valarray<double> > _subEventWeights;
	size_t _numWeightTypes; // always == WeightVector.size()

	/// Run name
	std::string _runname;

	/// Event counter
	mutable CounterPtr _eventCounter;

	/// Cross-section known to AH
	Scatter1DPtr _xs;

	/// 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;

	/// Flag to check if multiweights should be included
	bool _skipWeights;

	/// Current event number
	int _eventNumber;

	/// The index in the weight vector for the nominal weight stream
	size_t _defaultWeightIdx;

	/// Determines how often Rivet runs finalize() and writes the
	/// result to a YODA file.
	int _dumpPeriod;

	/// The name of a YODA file to which Rivet periodically dumps
	/// results.
	string _dumpFile;

	/// Flag to indicate periodic dumping is in progress
	bool _dumping;

	//@}


	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/Tools/RivetHepMC.hh b/include/Rivet/Tools/RivetHepMC.hh
	--- a/include/Rivet/Tools/RivetHepMC.hh
	+++ b/include/Rivet/Tools/RivetHepMC.hh
	@@ -1,105 +1,105 @@
	// -- C++ --
	#ifndef RIVET_RivetHepMC_HH
	#define RIVET_RivetHepMC_HH

	#include <valarray>

	#ifdef ENABLE_HEPMC_3
	#include "HepMC3/HepMC3.h"
	#include "HepMC3/Relatives.h"
	#include "HepMC3/Reader.h"

	namespace Rivet{
	namespace RivetHepMC = HepMC3;
	using RivetHepMC::ConstGenParticlePtr;
	using RivetHepMC::ConstGenVertexPtr;
	using RivetHepMC::Relatives;
	using RivetHepMC::ConstGenHeavyIonPtr;

	using HepMC_IO_type = RivetHepMC::Reader;

	using PdfInfo = RivetHepMC::GenPdfInfo;
	}

	#else
	#include "HepMC/GenEvent.h"
	#include "HepMC/GenParticle.h"
	#include "HepMC/HeavyIon.h"
	#include "HepMC/GenVertex.h"
	#include "HepMC/Version.h"
	#include "HepMC/GenRanges.h"
	#include "HepMC/IO_GenEvent.h"

	namespace Rivet{
	namespace RivetHepMC = HepMC;


	// HepMC 2.07 provides its own #defines
	typedef const HepMC::GenParticle* ConstGenParticlePtr;
	typedef const HepMC::GenVertex* ConstGenVertexPtr;
	typedef const HepMC::HeavyIon* ConstGenHeavyIonPtr;

	/// @brief Replicated the HepMC3 Relatives syntax using HepMC2 IteratorRanges
	/// This is necessary mainly because of capitalisation differences
	class Relatives{

	public:

	constexpr Relatives(HepMC::IteratorRange relo): _internal(relo){}

	constexpr HepMC::IteratorRange operator()() const {return _internal;}
	operator HepMC::IteratorRange() const {return _internal;}

	const static Relatives PARENTS;
	const static Relatives CHILDREN;
	const static Relatives ANCESTORS;
	const static Relatives DESCENDANTS;

	private:
	const HepMC::IteratorRange _internal;

	};

	using HepMC_IO_type = HepMC::IO_GenEvent;
	using PdfInfo = RivetHepMC::PdfInfo;

	}

	#endif

	#include "Rivet/Tools/RivetSTL.hh"
	#include "Rivet/Tools/Exceptions.hh"


	namespace Rivet {

	using RivetHepMC::GenEvent;
	using ConstGenEventPtr = std::shared_ptr<const GenEvent>;
	/// @todo Use mcutils?

	namespace HepMCUtils{

	ConstGenParticlePtr getParticlePtr(const RivetHepMC::GenParticle & gp);
	std::vector<ConstGenParticlePtr> particles(ConstGenEventPtr ge);
	std::vector<ConstGenParticlePtr> particles(const GenEvent *ge);
	std::vector<ConstGenVertexPtr> vertices(ConstGenEventPtr ge);
	std::vector<ConstGenVertexPtr> vertices(const GenEvent *ge);
	std::vector<ConstGenParticlePtr> particles(ConstGenVertexPtr gv, const Relatives &relo);
	std::vector<ConstGenParticlePtr> particles(ConstGenParticlePtr gp, const Relatives &relo);
	int uniqueId(ConstGenParticlePtr gp);
	int particles_size(ConstGenEventPtr ge);
	int particles_size(const GenEvent *ge);
	std::pair<ConstGenParticlePtr,ConstGenParticlePtr> beams(const GenEvent *ge);
	std::shared_ptr<HepMC_IO_type> makeReader(std::istream &istr,
	std::string * errm = 0);
	bool readEvent(std::shared_ptr<HepMC_IO_type> io,
	std::shared_ptr<GenEvent> evt);
	void strip(GenEvent & ge,
	const set<long> & stripid = {1, -1, 2, -2, 3,-3, 21});
	vector<string> weightNames(const GenEvent & ge);
	- double crossSection(const GenEvent & ge);
	+ pair<double,double> crossSection(const GenEvent & ge);
	std::valarray<double> weights(const GenEvent & ge);
	}
	}

	#endif
	diff --git a/src/Core/AnalysisHandler.cc b/src/Core/AnalysisHandler.cc
	--- a/src/Core/AnalysisHandler.cc
	+++ b/src/Core/AnalysisHandler.cc
	@@ -1,696 +1,684 @@
	// -- 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/IO.h"
	#include <iostream>

	using std::cout;
	using std::cerr;

	namespace Rivet {


	AnalysisHandler::AnalysisHandler(const string& runname)
	: _runname(runname),
	_initialised(false), _ignoreBeams(false), _skipWeights(false),
	_defaultWeightIdx(0), _dumpPeriod(0), _dumping(false)
	{ }


	AnalysisHandler::~AnalysisHandler()
	{ }


	Log& AnalysisHandler::getLog() const {
	return Log::getLog("Rivet.Analysis.Handler");
	}


	/// http://stackoverflow.com/questions/4654636/how-to-determine-if-a-string-is-a-number-with-c
	namespace {
	bool is_number(const std::string& s) {
	std::string::const_iterator it = s.begin();
	while (it != s.end() && std::isdigit(*it)) ++it;
	return !s.empty() && it == s.end();
	}
	}

	/// Check if any of the weightnames is not a number
	bool AnalysisHandler::haveNamedWeights() const {
	bool dec=false;
	for (unsigned int i=0;i<_weightNames.size();++i) {
	string s = _weightNames[i];
	if (!is_number(s)) {
	dec=true;
	break;
	}
	}
	return dec;
	}


	void AnalysisHandler::init(const GenEvent& ge) {
	if (_initialised)
	throw UserError("AnalysisHandler::init has already been called: cannot re-initialize!");

	/// @todo Should the Rivet analysis objects know about weight names?

	setRunBeams(Rivet::beams(ge));
	MSG_DEBUG("Initialising the analysis handler");
	_eventNumber = ge.event_number();

	setWeightNames(ge);
	if (_skipWeights)
	MSG_INFO("Only using nominal weight. Variation weights will be ignored.");
	else if (haveNamedWeights())
	MSG_INFO("Using named weights");
	else
	MSG_INFO("NOT using named weights. Using first weight as nominal weight");

	_eventCounter = CounterPtr(weightNames(), Counter("_EVTCOUNT"));

	// Set the cross section based on what is reported by this event.
	- if (ge.cross_section()) {
	- MSG_TRACE("Getting cross section.");
	- double xs = ge.cross_section()->cross_section();
	- double xserr = ge.cross_section()->cross_section_error();
	- setCrossSection(xs, xserr);
	- }
	+ if ( ge.cross_section() ) setCrossSection(HepMCUtils::crossSection(ge));

	// 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 ( a->info().preliminary() ) {
	MSG_WARNING("Analysis '" << a->name() << "' is preliminary: be careful, it may change and/or be renamed!");
	} else if ( a->info().obsolete() ) {
	MSG_WARNING("Analysis '" << a->name() << "' is obsolete: please update!");
	} else if (( a->info().unvalidated() ) ) {
	MSG_WARNING("Analysis '" << a->name() << "' is unvalidated: be careful, it may be broken!");
	}
	}

	// Initialize the remaining analyses
	_stage = Stage::INIT;
	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());
	}
	_stage = Stage::OTHER;
	_initialised = true;
	MSG_DEBUG("Analysis handler initialised");
	}

	void AnalysisHandler::setWeightNames(const GenEvent& ge) {
	if (!_skipWeights) _weightNames = HepMCUtils::weightNames(ge);
	if ( _weightNames.empty() ) _weightNames.push_back("");
	for ( int i = 0, N = _weightNames.size(); i < N; ++i )
	if ( _weightNames[i] == "" ) _defaultWeightIdx = i;
	}

	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 (if we're checking beams)
	if ( !_ignoreBeams ) {
	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
	bool strip = ( getEnvParam("RIVET_STRIP_HEPMC", string("NOOOO") ) != "NOOOO" );
	Event event(ge, strip);

	// set the cross section based on what is reported by this event.
	// if no cross section
	- MSG_TRACE("getting cross section.");
	- if (ge.cross_section()) {
	- MSG_TRACE("getting cross section from GenEvent.");
	- double xs = ge.cross_section()->cross_section();
	- double xserr = ge.cross_section()->cross_section_error();
	- setCrossSection(xs, xserr);
	- }
	+ if ( ge.cross_section() ) setCrossSection(HepMCUtils::crossSection(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?

	- MSG_TRACE("AnalysisHandler::analyze(): Pushing _eventCounter to persistent.");
	- _eventCounter.get()->pushToPersistent(_subEventWeights);
	- // if this is indeed a new event, push the temporary
	- // histograms and reset
	- for (const AnaHandle& a : analyses()) {
	- for (auto ao : a->analysisObjects()) {
	- MSG_TRACE("AnalysisHandler::analyze(): Pushing " << a->name() << "'s " << ao->name() << " to persistent.");
	- ao.get()->pushToPersistent(_subEventWeights);
	- }
	- MSG_TRACE("AnalysisHandler::analyze(): finished pushing " << a->name() << "'s objects to persistent.");
	- }
	+ pushToPersistent();

	- _eventNumber = ge.event_number();
	+ _eventNumber = ge.event_number();

	- MSG_DEBUG("nominal event # " << _eventCounter.get()->_persistent[0]->numEntries());
	- MSG_DEBUG("nominal sum of weights: " << _eventCounter.get()->_persistent[0]->sumW());
	- MSG_DEBUG("Event has " << _subEventWeights.size() << " sub events.");
	- _subEventWeights.clear();
	}


	MSG_TRACE("starting new sub event");
	_eventCounter.get()->newSubEvent();

	for (const AnaHandle& a : analyses()) {
	for (auto ao : a->analysisObjects()) {
	ao.get()->newSubEvent();
	}
	}

	_subEventWeights.push_back(event.weights());
	MSG_DEBUG("Analyzing subevent #" << _subEventWeights.size() - 1 << ".");

	_eventCounter->fill();
	// 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());
	}

	if ( _dumpPeriod > 0 && numEvents() > 0 && numEvents()%_dumpPeriod == 0 ) {
	MSG_INFO("Dumping intermediate results to " << _dumpFile << ".");
	_dumping = numEvents()/_dumpPeriod;
	finalize();
	writeData(_dumpFile);
	_dumping = 0;
	}

	}


	void AnalysisHandler::analyze(const GenEvent* ge) {
	if (ge == nullptr) {
	MSG_ERROR("AnalysisHandler received null pointer to GenEvent");
	//throw Error("AnalysisHandler received null pointer to GenEvent");
	}
	analyze(*ge);
	}

	+ void AnalysisHandler::pushToPersistent() {
	+ if ( _subEventWeights.empty() ) return;
	+ MSG_TRACE("AnalysisHandler::analyze(): Pushing _eventCounter to persistent.");
	+ _eventCounter.get()->pushToPersistent(_subEventWeights);
	+ for (const AnaHandle& a : analyses()) {
	+ for (auto ao : a->analysisObjects()) {
	+ MSG_TRACE("AnalysisHandler::analyze(): Pushing " << a->name()
	+ << "'s " << ao->name() << " to persistent.");
	+ ao.get()->pushToPersistent(_subEventWeights);
	+ }
	+ MSG_TRACE("AnalysisHandler::analyze(): finished pushing "
	+ << a->name() << "'s objects to persistent.");
	+ }
	+ _subEventWeights.clear();
	+ }

	void AnalysisHandler::finalize() {
	if (!_initialised) return;
	MSG_INFO("Finalising analyses");

	_stage = Stage::FINALIZE;

	// First push all analyses' objects to persistent and final
	MSG_TRACE("AnalysisHandler::finalize(): Pushing analysis objects to persistent.");
	- _eventCounter.get()->pushToPersistent(_subEventWeights);
	+ pushToPersistent();
	+
	+ // Copy all histos to finalize versions.
	_eventCounter.get()->pushToFinal();
	_xs.get()->pushToFinal();
	- for (const AnaHandle& a : analyses()) {
	- for (auto ao : a->analysisObjects()) {
	- ao.get()->pushToPersistent(_subEventWeights);
	+ for (const AnaHandle& a : analyses())
	+ for (auto ao : a->analysisObjects())
	ao.get()->pushToFinal();
	- }
	- }

	for (AnaHandle a : analyses()) {
	if ( _dumping && !a->info().reentrant() ) {
	if ( _dumping == 1 )
	MSG_INFO("Skipping finalize in periodic dump of " << a->name()
	<< " as it is not declared reentrant.");
	continue;
	}
	for (size_t iW = 0; iW < numWeights(); iW++) {
	_eventCounter.get()->setActiveFinalWeightIdx(iW);
	_xs.get()->setActiveFinalWeightIdx(iW);
	for (auto ao : a->analysisObjects())
	ao.get()->setActiveFinalWeightIdx(iW);
	try {
	MSG_TRACE("running " << a->name() << "::finalize() for weight " << iW << ".");
	a->finalize();
	} catch (const Error& err) {
	cerr << "Error in " << a->name() << "::finalize method: " << err.what() << '\n';
	exit(1);
	}
	}
	}

	// Print out number of events processed
	const int nevts = numEvents();
	MSG_INFO("Processed " << nevts << " event" << (nevts != 1 ? "s" : ""));

	_stage = Stage::OTHER;

	if ( _dumping ) return;

	// Print out MCnet boilerplate
	if (getLog().getLevel()<=20){
	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, std::map<string, string> pars) {
	// Make an option handle.
	std::string parHandle = "";
	for (map<string, string>::iterator par = pars.begin(); par != pars.end(); ++par) {
	parHandle +=":";
	parHandle += par->first + "=" + par->second;
	}
	return addAnalysis(analysisname + parHandle);
	}


	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.
	string ananame = analysisname;
	vector<string> anaopt = split(analysisname, ":");
	if ( anaopt.size() > 1 ) ananame = anaopt[0];
	AnaHandle analysis( AnalysisLoader::getAnalysis(ananame) );
	if (analysis.get() != 0) { // < Check for null analysis.
	MSG_DEBUG("Adding analysis '" << analysisname << "'");
	map<string,string> opts;
	for ( int i = 1, N = anaopt.size(); i < N; ++i ) {
	vector<string> opt = split(anaopt[i], "=");
	if ( opt.size() != 2 ) {
	MSG_WARNING("Error in option specification. Skipping analysis " << analysisname);
	return *this;
	}
	if ( !analysis->info().validOption(opt[0], opt[1]) ) {
	MSG_WARNING("Cannot set option '" << opt[0] << "' to '" << opt[1]
	<< "'. Skipping analysis " << analysisname);
	return *this;
	}
	opts[opt[0]] = opt[1];
	}
	for ( auto opt: opts) {
	analysis->_options[opt.first] = opt.second;
	analysis->_optstring += ":" + opt.first + "=" + opt.second;
	}
	for (const AnaHandle& a : analyses()) {
	if (a->name() == analysis->name() ) {
	MSG_WARNING("Analysis '" << analysisname << "' already registered: skipping duplicate");
	return *this;
	}
	}
	analysis->_analysishandler = this;
	_analyses[analysisname] = 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) {
	MSG_DEBUG("Removing analysis '" << analysisname << "'");
	if (_analyses.find(analysisname) != _analyses.end()) _analyses.erase(analysisname);
	// }
	return *this;
	}


	// void AnalysisHandler::addData(const std::vector<YODA::AnalysisObjectPtr>& aos) {
	// for (const YODA::AnalysisObjectPtr ao : aos) {
	// string path = ao->path();
	// if ( path.substr(0, 5) != "/RAW/" ) {
	// _orphanedPreloads.push_back(ao);
	// continue;
	// }

	// path = path.substr(4);
	// ao->setPath(path);
	// if (path.size() > 1) { // path > "/"
	// try {
	// const string ananame = ::split(path, "/")[0];
	// AnaHandle a = analysis(ananame);
	// /// @todo FIXXXXX
	// //MultiweightAOPtr mao = ????; /// @todo generate right Multiweight object from ao
	// //a->addAnalysisObject(mao); /// @todo Need to statistically merge...
	// } catch (const Error& e) {
	// MSG_TRACE("Adding analysis object " << path <<
	// " to the list of orphans.");
	// _orphanedPreloads.push_back(ao);
	// }
	// }
	// }
	// }


	void AnalysisHandler::stripOptions(YODA::AnalysisObjectPtr ao,
	const vector<string> & delopts) const {
	string path = ao->path();
	string ananame = split(path, "/")[0];
	vector<string> anaopts = split(ananame, ":");
	for ( int i = 1, N = anaopts.size(); i < N; ++i )
	for ( auto opt : delopts )
	if ( opt == "*" \|\| anaopts[i].find(opt + "=") == 0 )
	path.replace(path.find(":" + anaopts[i]), (":" + anaopts[i]).length(), "");
	ao->setPath(path);
	}


	void AnalysisHandler::mergeYodas(const vector<string> & aofiles,
	const vector<string> & delopts, bool equiv) {

	// Convenient typedef;
	typedef multimap<string, YODA::AnalysisObjectPtr> AOMap;

	// Store all found weights here.
	set<string> foundWeightNames;

	// Stor all found analyses.
	set<string> foundAnalyses;

	// Store all analysis objects here.
	vector<AOMap> allaos;

	// Go through all files and collect information.
	for ( auto file : aofiles ) {
	allaos.push_back(AOMap());
	AOMap & aomap = allaos.back();
	vector<YODA::AnalysisObject*> aos_raw;
	try {
	YODA::read(file, aos_raw);
	}
	catch (...) { //< YODA::ReadError&
	throw UserError("Unexpected error in reading file: " + file);
	}
	for (YODA::AnalysisObject* aor : aos_raw) {
	YODA::AnalysisObjectPtr ao(aor);
	AOPath path(ao->path());
	if ( !path )
	throw UserError("Invalid path name in file: " + file);
	if ( !path.isRaw() ) continue;

	foundWeightNames.insert(path.weight());
	// Now check if any options should be removed.
	for ( string delopt : delopts )
	if ( path.hasOption(delopt) ) path.removeOption(delopt);
	path.setPath();
	if ( path.analysisWithOptions() != "" )
	foundAnalyses.insert(path.analysisWithOptions());
	aomap.insert(make_pair(path.path(), ao));
	}
	}

	// Now make analysis handler aware of the weight names present.
	_weightNames.clear();
	_weightNames.push_back("");
	_defaultWeightIdx = 0;
	for ( string name : foundWeightNames ) _weightNames.push_back(name);

	// Then we create and initialize all analyses
	for ( string ananame : foundAnalyses ) addAnalysis(ananame);
	for (AnaHandle a : analyses() ) {
	MSG_TRACE("Initialising analysis: " << a->name());
	if ( !a->info().reentrant() )
	MSG_WARNING("Analysis " << a->name() << " has not been validated to have "
	<< "a reentrant finalize method. The merged result is unpredictable.");
	try {
	// Allow projection registration in the init phase onwards
	a->_allowProjReg = true;
	a->init();
	} catch (const Error& err) {
	cerr << "Error in " << a->name() << "::init method: " << err.what() << endl;
	exit(1);
	}
	MSG_TRACE("Done initialising analysis: " << a->name());
	}

	// Now get all booked analysis objects.
	vector<MultiweightAOPtr> raos;
	for (AnaHandle a : analyses()) {
	for (const auto & ao : a->analysisObjects()) {
	raos.push_back(ao);
	}
	}

	// Collect global weights and xcoss sections and fix scaling for
	// all files.
	_eventCounter = CounterPtr(weightNames(), Counter("_EVTCOUNT"));
	_xs = Scatter1DPtr(weightNames(), Scatter1D("_XSEC"));
	for (size_t iW = 0; iW < numWeights(); iW++) {
	_eventCounter.get()->setActiveWeightIdx(iW);
	_xs.get()->setActiveWeightIdx(iW);
	YODA::Counter & sumw = *_eventCounter;
	YODA::Scatter1D & xsec = *_xs;
	vector<YODA::Scatter1DPtr> xsecs;
	vector<YODA::CounterPtr> sows;
	for ( auto & aomap : allaos ) {
	auto xit = aomap.find(xsec.path());
	if ( xit == aomap.end() )
	xsecs.push_back(dynamic_pointer_cast<YODA::Scatter1D>(xit->second));
	else
	xsecs.push_back(YODA::Scatter1DPtr());
	xit = aomap.find(sumw.path());
	if ( xit == aomap.end() )
	sows.push_back(dynamic_pointer_cast<YODA::Counter>(xit->second));
	else
	sows.push_back(YODA::CounterPtr());
	}
	double xs = 0.0, xserr = 0.0;
	for ( int i = 0, N = sows.size(); i < N; ++i ) {
	if ( !sows[i] \|\| !xsecs[i] ) continue;
	double xseci = xsecs[i]->point(0).x();
	double xsecerri = sqr(xsecs[i]->point(0).xErrAvg());
	sumw += *sows[i];
	double effnent = sows[i]->effNumEntries();
	xs += (equiv? effnent: 1.0)*xseci;
	xserr += (equiv? sqr(effnent): 1.0)*xsecerri;
	}
	vector<double> scales(sows.size(), 1.0);
	if ( equiv ) {
	xs /= sumw.effNumEntries();
	xserr = sqrt(xserr)/sumw.effNumEntries();
	} else {
	xserr = sqrt(xserr);
	for ( int i = 0, N = sows.size(); i < N; ++i )
	scales[i] = (sumw.sumW()/sows[i]->sumW())*
	(xsecs[i]->point(0).x()/xs);
	}
	xsec.point(0) = Point1D(xs, xserr);

	// Go through alla analyses and add stuff to their analysis objects;
	for (AnaHandle a : analyses()) {
	for (const auto & ao : a->analysisObjects()) {
	ao.get()->setActiveWeightIdx(iW);
	YODA::AnalysisObjectPtr yao = ao.get()->activeYODAPtr();
	for ( int i = 0, N = sows.size(); i < N; ++i ) {
	if ( !sows[i] \|\| !xsecs[i] ) continue;
	auto range = allaos[i].equal_range(yao->path());
	for ( auto aoit = range.first; aoit != range.second; ++aoit )
	if ( !addaos(yao, aoit->second, scales[i]) )
	MSG_WARNING("Cannot merge objects with path " << yao->path()
	<<" of type " << yao->annotation("Type") );
	}
	ao.get()->unsetActiveWeight();
	}
	}
	_eventCounter.get()->unsetActiveWeight();
	_xs.get()->unsetActiveWeight();
	}

	// Finally we just have to finalize all analyses, leaving to the
	// controlling program to write it out some yoda-file.
	finalize();

	}

	void AnalysisHandler::readData(const string& filename) {
	try {
	/// @todo Use new YODA SFINAE to fill the smart ptr vector directly
	vector<YODA::AnalysisObject*> aos_raw;
	YODA::read(filename, aos_raw);
	for (YODA::AnalysisObject* aor : aos_raw)
	_preloads[aor->path()] = YODA::AnalysisObjectPtr(aor);
	} catch (...) { //< YODA::ReadError&
	throw UserError("Unexpected error in reading file: " + filename);
	}
	}

	vector<MultiweightAOPtr> AnalysisHandler::getRivetAOs() const {
	vector<MultiweightAOPtr> rtn;

	for (AnaHandle a : analyses()) {
	for (const auto & ao : a->analysisObjects()) {
	rtn.push_back(ao);
	}
	}

	rtn.push_back(_eventCounter);
	rtn.push_back(_xs);

	return rtn;
	}

	void AnalysisHandler::writeData(const string& filename) const {

	// This is where we store the OAs to be written.
	vector<YODA::AnalysisObjectPtr> output;

	// First get all multiwight AOs
	vector<MultiweightAOPtr> raos = getRivetAOs();
	output.reserve(raos.size()2numWeights());

	// Fix the oredering so that default weight is written out first.
	vector<size_t> order;
	if ( _defaultWeightIdx >= 0 && _defaultWeightIdx < numWeights() )
	order.push_back(_defaultWeightIdx);
	for ( size_t i = 0; i < numWeights(); ++i )
	if ( i != _defaultWeightIdx ) order.push_back(i);

	// Then we go through all finalized AOs one weight at a time
	for (size_t iW : order ) {
	for ( auto rao : raos ) {
	rao.get()->setActiveFinalWeightIdx(iW);
	if ( rao->path().find("/TMP/") != string::npos ) continue;
	output.push_back(rao.get()->activeYODAPtr());
	}
	}

	// Finally the RAW objects.
	for (size_t iW : order ) {
	for ( auto rao : raos ) {
	rao.get()->setActiveWeightIdx(iW);
	output.push_back(rao.get()->activeYODAPtr());
	}
	}

	try {
	YODA::write(filename, output.begin(), output.end());
	} catch (...) { //< YODA::WriteError&
	throw UserError("Unexpected error in writing file: " + filename);
	}
	}


	string AnalysisHandler::runName() const { return _runname; }
	size_t AnalysisHandler::numEvents() const { return _eventCounter->numEntries(); }


	std::vector<std::string> AnalysisHandler::analysisNames() const {
	std::vector<std::string> rtn;
	for (AnaHandle a : analyses()) {
	rtn.push_back(a->name());
	}
	return rtn;
	}


	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;
	}


	- AnalysisHandler& AnalysisHandler::setCrossSection(double xs, double xserr) {
	+ void AnalysisHandler::setCrossSection(pair<double,double> xsec) {
	_xs = Scatter1DPtr(weightNames(), Scatter1D("_XSEC"));
	_eventCounter.get()->setActiveWeightIdx(_defaultWeightIdx);
	double nomwgt = sumW();

	// The cross section of each weight variation is the nominal cross section
	// times the sumW(variation) / sumW(nominal).
	// This way the cross section will work correctly
	for (size_t iW = 0; iW < numWeights(); iW++) {
	_eventCounter.get()->setActiveWeightIdx(iW);
	double s = sumW() / nomwgt;
	_xs.get()->setActiveWeightIdx(iW);
	- _xs->addPoint(xss, xserrs);
	+ _xs->addPoint(xsec.firsts, xsec.seconds);
	}

	_eventCounter.get()->unsetActiveWeight();
	_xs.get()->unsetActiveWeight();
	- return *this;
	+ return;
	}

	AnalysisHandler& AnalysisHandler::addAnalysis(Analysis* analysis) {
	analysis->_analysishandler = this;
	// _analyses.insert(AnaHandle(analysis));
	_analyses[analysis->name()] = 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;
	}

	void AnalysisHandler::skipMultiWeights(bool ignore) {
	_skipWeights=ignore;
	}


	}
	diff --git a/src/Core/Run.cc b/src/Core/Run.cc
	--- a/src/Core/Run.cc
	+++ b/src/Core/Run.cc
	@@ -1,142 +1,143 @@
	// -- C++ --
	#include "Rivet/Run.hh"
	#include "Rivet/AnalysisHandler.hh"
	#include "Rivet/Math/MathUtils.hh"
	#include "Rivet/Tools/RivetPaths.hh"
	#include "zstr/zstr.hpp"
	#include <limits>
	#include <iostream>

	using std::cout;
	using std::endl;

	namespace Rivet {


	Run::Run(AnalysisHandler& ah)
	: _ah(ah), _fileweight(1.0), _xs(NAN)
	{ }


	Run::~Run() { }


	Run& Run::setCrossSection(const double xs) {
	_xs = xs;
	return *this;
	}


	Run& Run::setListAnalyses(const bool dolist) {
	_listAnalyses = dolist;
	return *this;
	}


	// Fill event and check for a bad read state
	bool Run::readEvent() {
	/// @todo Clear rather than new the GenEvent object per-event?
	_evt.reset(new GenEvent());
	if(!HepMCUtils::readEvent(_hepmcReader, _evt)){
	Log::getLog("Rivet.Run") << Log::DEBUG << "Read failed. End of file?" << endl;
	return false;
	}
	// Rescale event weights by file-level weight, if scaling is non-trivial
	if (!fuzzyEquals(_fileweight, 1.0)) {
	for (size_t i = 0; i < (size_t) _evt->weights().size(); ++i) {
	_evt->weights()[i] *= _fileweight;
	}
	}
	return true;
	}



	bool Run::openFile(const std::string& evtfile, double weight) {
	// Set current weight-scaling member
	_fileweight = weight;

	// In case makeReader fails.
	std::string errormessage;

	// Set up HepMC input reader objects
	if (evtfile == "-") {
	#ifdef HAVE_LIBZ
	_istr = make_shared<zstr::istream>(std::cin);
	_hepmcReader = HepMCUtils::makeReader(*_istr, &errormessage);
	#else
	_hepmcReader = HepMCUtils::makeReader(std::cin, &errormessage);
	#endif
	} else {
	if ( !fileexists(evtfile) )
	throw Error("Event file '" + evtfile + "' not found");
	#ifdef HAVE_LIBZ
	// NB. zstr auto-detects if file is deflated or plain-text
	_istr = make_shared<zstr::ifstream>(evtfile.c_str());
	#else
	_istr = make_shared<std::ifstream>(evtfile.c_str());
	#endif
	_hepmcReader = HepMCUtils::makeReader(*_istr, &errormessage);

	}

	if (_hepmcReader == nullptr) {
	Log::getLog("Rivet.Run")
	<< Log::ERROR << "Read error on file '" << evtfile << "' "
	<< errormessage << endl;
	return false;
	}
	return true;
	}


	bool Run::init(const std::string& evtfile, double weight) {
	if (!openFile(evtfile, weight)) return false;

	// Read first event to define run conditions
	bool ok = readEvent();
	if (!ok) return false;
	if(HepMCUtils::particles(_evt).size() == 0){
	Log::getLog("Rivet.Run") << Log::ERROR << "Empty first event." << endl;
	return false;
	}

	// Initialise AnalysisHandler with beam information from first event
	_ah.init(*_evt);

	// Set cross-section from command line
	if (!std::isnan(_xs)) {
	Log::getLog("Rivet.Run")
	<< Log::DEBUG << "Setting user cross-section = " << _xs << " pb" << endl;
	- _ah.setCrossSection(_xs, 0.0);
	+
	+ _ah.setCrossSection(make_pair(_xs, 0.0));
	}

	// List the chosen & compatible analyses if requested
	if (_listAnalyses) {
	for (const std::string& ana : _ah.analysisNames()) {
	cout << ana << endl;
	}
	}

	return true;
	}


	bool Run::processEvent() {
	// Analyze event
	_ah.analyze(*_evt);

	return true;
	}


	bool Run::finalize() {
	_evt.reset();

	_ah.finalize();

	return true;
	}


	}
	diff --git a/src/Tools/RivetHepMC_2.cc b/src/Tools/RivetHepMC_2.cc
	--- a/src/Tools/RivetHepMC_2.cc
	+++ b/src/Tools/RivetHepMC_2.cc
	@@ -1,200 +1,201 @@
	// -- C++ --

	//#include <regex>
	#include "Rivet/Tools/Utils.hh"
	#include "Rivet/Tools/RivetHepMC.hh"
	#include "Rivet/Tools/Logging.hh"

	/*namespace {

	inline std::vector<std::string> split(const std::string& input, const std::string& regex) {
	// passing -1 as the submatch index parameter performs splitting
	std::regex re(regex);
	std::sregex_token_iterator
	first{input.begin(), input.end(), re, -1},
	last;
	return {first, last};
	}

	}*/

	namespace Rivet{

	const Relatives Relatives::PARENTS = HepMC::parents;
	const Relatives Relatives::CHILDREN = HepMC::children;
	const Relatives Relatives::ANCESTORS = HepMC::ancestors;
	const Relatives Relatives::DESCENDANTS = HepMC::descendants;

	namespace HepMCUtils{

	ConstGenParticlePtr getParticlePtr(const RivetHepMC::GenParticle & gp) {
	return &gp;
	}
	std::vector<ConstGenParticlePtr> particles(ConstGenEventPtr ge){
	std::vector<ConstGenParticlePtr> result;
	for(GenEvent::particle_const_iterator pi = ge->particles_begin(); pi != ge->particles_end(); ++pi){
	result.push_back(*pi);
	}
	return result;
	}

	std::vector<ConstGenParticlePtr> particles(const GenEvent *ge){
	std::vector<ConstGenParticlePtr> result;
	for(GenEvent::particle_const_iterator pi = ge->particles_begin(); pi != ge->particles_end(); ++pi){
	result.push_back(*pi);
	}
	return result;
	}

	std::vector<ConstGenVertexPtr> vertices(ConstGenEventPtr ge){
	std::vector<ConstGenVertexPtr> result;
	for(GenEvent::vertex_const_iterator vi = ge->vertices_begin(); vi != ge->vertices_end(); ++vi){
	result.push_back(*vi);
	}
	return result;
	}

	std::vector<ConstGenVertexPtr> vertices(const GenEvent *ge){
	std::vector<ConstGenVertexPtr> result;
	for(GenEvent::vertex_const_iterator vi = ge->vertices_begin(); vi != ge->vertices_end(); ++vi){
	result.push_back(*vi);
	}
	return result;
	}

	std::vector<ConstGenParticlePtr> particles(ConstGenVertexPtr gv, const Relatives &relo){
	std::vector<ConstGenParticlePtr> result;
	/// @todo A particle_const_iterator on GenVertex would be nice...
	// Before HepMC 2.7.0 there were no GV::particles_const_iterators and constness consistency was all screwed up :-/
	#if HEPMC_VERSION_CODE >= 2007000
	for (HepMC::GenVertex::particle_iterator pi = gv->particles_begin(relo); pi != gv->particles_end(relo); ++pi)
	result.push_back(*pi);
	#else
	HepMC::GenVertex* gv2 = const_cast<HepMC::GenVertex*>(gv);
	for (HepMC::GenVertex::particle_iterator pi = gv2->particles_begin(relo); pi != gv2->particles_end(relo); ++pi)
	result.push_back(const_cast<ConstGenParticlePtr>(*pi));
	#endif
	return result;
	}

	std::vector<ConstGenParticlePtr> particles(ConstGenParticlePtr gp, const Relatives &relo){
	ConstGenVertexPtr vtx;

	switch(relo){
	case HepMC::parents:

	case HepMC::ancestors:
	vtx = gp->production_vertex();
	break;

	case HepMC::children:

	case HepMC::descendants:
	vtx = gp->end_vertex();
	break;

	default:

	throw std::runtime_error("Not implemented!");
	break;
	}

	return particles(vtx, relo);
	}



	int uniqueId(ConstGenParticlePtr gp){
	return gp->barcode();
	}

	int particles_size(ConstGenEventPtr ge){
	return ge->particles_size();
	}

	int particles_size(const GenEvent *ge){
	return ge->particles_size();
	}

	std::pair<ConstGenParticlePtr,ConstGenParticlePtr> beams(const GenEvent *ge){
	return ge->beam_particles();
	}

	std::shared_ptr<HepMC::IO_GenEvent> makeReader(std::istream &istr,
	std::string *) {
	return make_shared<HepMC::IO_GenEvent>(istr);
	}

	bool readEvent(std::shared_ptr<HepMC::IO_GenEvent> io, std::shared_ptr<GenEvent> evt){
	if(io->rdstate() != 0) return false;
	if(!io->fill_next_event(evt.get())) return false;
	return true;
	}

	// This functions could be filled with code doing the same stuff as
	// in the HepMC3 version of This file.
	void strip(GenEvent &, const set<long> &) {}

	vector<string> weightNames(const GenEvent & ge) {
	/// reroute the print output to a std::stringstream and process
	/// The iteration is done over a map in hepmc2 so this is safe
	vector<string> ret;

	/// Obtaining weight names using regex probably neater, but regex
	/// is not defined in GCC4.8, which is currently used by Lxplus.
	/// Attempt an alternative solution based on stringstreams:
	std::stringstream stream;
	ge.weights().print(stream);
	std::string pair; // placeholder for subtsring matches
	while (std::getline(stream, pair, ' ')) {
	if ( pair.size() < 2 ) continue;
	pair.erase(pair.begin()); // removes the "(" on the LHS
	pair.pop_back(); // removes the ")" on the RHS
	if (pair.empty()) continue;
	std::stringstream spair(pair);
	vector<string> temp;
	while (std::getline(spair, pair, ',')) {
	temp.push_back(std::move(pair));
	}
	if (temp.size() == 2) {
	// store the default weight based on weight names
	if (temp[0] == "Weight" \|\| temp[0] == "0" \|\| temp[0] == "Default") {
	ret.push_back("");
	}
	else ret.push_back(temp[0]);
	}
	}
	/// Possible future solution based on regex
	/*std::ostringstream stream;
	ge.weights().print(stream); // Super lame, I know
	string str = stream.str();

	std::regex re("(([^()]+))"); // Regex for stuff enclosed by parentheses ()
	for (std::sregex_iterator i = std::sregex_iterator(str.begin(), str.end(), re);
	i != std::sregex_iterator(); ++i ) {
	std::smatch m = *i;
	vector<string> temp = ::split(m.str(), "[,]");
	if (temp.size() ==2) {
	// store the default weight based on weight names
	if (temp[0] == "Weight" \|\| temp[0] == "0" \|\| temp[0] == "Default") {
	ret.push_back("");
	} else
	ret.push_back(temp[0]);
	}
	}*/
	return ret;
	}

	- double crossSection(const GenEvent & ge) {
	- return ge.cross_section()->cross_section();
	+ pair<double,double> crossSection(const GenEvent & ge) {
	+ return make_pair(ge.cross_section()->cross_section(),
	+ ge.cross_section()->cross_section_error());
	}

	std::valarray<double> weights(const GenEvent & ge) {
	const size_t W = ge.weights().size();
	std::valarray<double> wts(W);
	for (unsigned int iw = 0; iw < W; ++iw)
	wts[iw] = ge.weights()[iw];
	return wts;
	}
	}
	}
	diff --git a/src/Tools/RivetHepMC_3.cc b/src/Tools/RivetHepMC_3.cc
	--- a/src/Tools/RivetHepMC_3.cc
	+++ b/src/Tools/RivetHepMC_3.cc
	@@ -1,159 +1,179 @@
	// -- C++ --

	#include "Rivet/Tools/RivetHepMC.hh"
	#include "Rivet/Tools/ReaderCompressedAscii.hh"
	#include "HepMC3/ReaderAscii.h"
	#include "HepMC3/ReaderAsciiHepMC2.h"
	+#include "HepMC3/GenCrossSection.h"
	+#include <cassert>

	namespace Rivet{

	namespace HepMCUtils{

	ConstGenParticlePtr getParticlePtr(const RivetHepMC::GenParticle & gp) {
	return gp.shared_from_this();
	}

	std::vector<ConstGenParticlePtr> particles(ConstGenEventPtr ge){
	return ge->particles();
	}

	std::vector<ConstGenParticlePtr> particles(const GenEvent *ge){
	assert(ge);
	return ge->particles();
	}

	std::vector<ConstGenVertexPtr> vertices(ConstGenEventPtr ge){
	return ge->vertices();
	}

	std::vector<ConstGenVertexPtr> vertices(const GenEvent *ge){
	assert(ge);
	return ge->vertices();
	}

	std::vector<ConstGenParticlePtr> particles(ConstGenVertexPtr gv, const Relatives &relo){
	return relo(gv);
	}

	std::vector<ConstGenParticlePtr> particles(ConstGenParticlePtr gp, const Relatives &relo){
	return relo(gp);
	}

	int particles_size(ConstGenEventPtr ge){
	return particles(ge).size();
	}

	int particles_size(const GenEvent *ge){
	return particles(ge).size();
	}

	int uniqueId(ConstGenParticlePtr gp){
	return gp->id();
	}

	std::pair<ConstGenParticlePtr,ConstGenParticlePtr> beams(const GenEvent *ge) {
	std::vector<ConstGenParticlePtr> beamlist = ge->beams();
	if ( beamlist.size() < 2 ) {
	std::cerr << "CANNOT FIND ANY BEAMS!" << std::endl;
	return std::pair<ConstGenParticlePtr,ConstGenParticlePtr>();
	}
	return std::make_pair(beamlist[0], beamlist[1]);
	}

	bool readEvent(std::shared_ptr<HepMC_IO_type> io, std::shared_ptr<GenEvent> evt){
	return io->read_event(*evt) && !io->failed();
	}

	shared_ptr<HepMC_IO_type> makeReader(std::istream & istr,
	std::string * errm) {
	shared_ptr<HepMC_IO_type> ret;

	// First scan forward and check if there is some hint as to what
	// kind of file we are looking att.
	int ntry = 10;
	std::string header;
	int filetype = -1;
	while ( ntry ) {
	std::getline(istr, header);
	if ( header.empty() ) continue;
	if ( header.substr(0, 34) == "HepMC::Asciiv3-START_EVENT_LISTING" ) {
	filetype = 3;
	break;
	}
	if ( header.substr(0, 44) == "HepMC::CompressedAsciiv3-START_EVENT_LISTING" ) {
	filetype = 4;
	break;
	}
	if ( header.substr(0, 38) == "HepMC::IO_GenEvent-START_EVENT_LISTING" ) {
	filetype = 2;
	break;
	}
	--ntry;
	}


	if ( filetype == 3 )
	ret = make_shared<RivetHepMC::ReaderAscii>(istr);
	else if ( filetype == 4 )
	ret = make_shared<Rivet::ReaderCompressedAscii>(istr);
	else
	ret = make_shared<RivetHepMC::ReaderAsciiHepMC2>(istr);
	if ( filetype == 0 && errm )
	*errm += "Could not determine file type. Assuming HepMC2 file. ";
	// Check that everything was ok.
	if ( ret->failed() ) {
	if ( errm ) *errm = "Problems reading from HepMC file.";
	ret = shared_ptr<HepMC_IO_type>();
	}

	return ret;
	}

	void strip(GenEvent & ge, const set<long> & stripid) {
	// std::cerr << "Stripping event " << ge.event_number() << std::endl;
	vector<HepMC3::GenParticlePtr> allparticles = ge.particles();
	for ( auto & p : allparticles ) {
	if ( !p->production_vertex() \|\| !p->end_vertex() \|\|
	stripid.count(p->pid()) == 0 \|\|
	p->production_vertex()->id() == 0 ) continue;
	// std::cout << "Removing particle " << p->id()
	// << " (" << p->pid() << ")" << std::endl;
	HepMC3::GenVertexPtr vp = p->production_vertex();
	HepMC3::GenVertexPtr ve = p->end_vertex();
	if ( !vp \|\| !ve ) continue;
	if ( vp == ve ) continue;
	// Check if the vertices would leave particles with the sam
	// production as decay vertex - we don't want that.
	if ( ( vp->particles_out().size() == 1 && vp->particles_out()[0] == p ) \|\|
	( ve->particles_in().size() == 1 && ve->particles_in()[0] == p ) ) {
	bool loop = false;
	for ( auto pi : vp->particles_in() )
	for ( auto po : ve->particles_out() )
	if ( pi == po ) loop = true;
	if ( loop ) continue;
	}
	if ( vp->particles_in().size() == 1 &&
	( vp->particles_in()[0]->pid() > 21 &&
	vp->particles_in()[0]->pid() < 30 ) )
	continue;

	vp->remove_particle_out(p);
	ve->remove_particle_in(p);

	if ( ve->particles_in().empty() ) {
	auto prem = ve->particles_out();
	for ( auto po : prem ) vp->add_particle_out(po);
	ge.remove_vertex(ve);
	}
	else if ( vp->particles_out().empty() ) {
	auto prem = vp->particles_in();
	for ( auto pi : prem ) ve->add_particle_in(pi);
	ge.remove_vertex(vp);
	}
	ge.remove_particle(p);
	}
	}

	+ pair<double,double> crossSection(const GenEvent & ge) {
	+ // Work-around since access functions are not const.
	+ HepMC3::GenCrossSection xs = *ge.cross_section();
	+ return make_pair(xs.xsec(), xs.xsec_err());
	+ }
	+
	+ vector<string> weightNames(const GenEvent & ge) {
	+ try {
	+ return ge.weight_names("");
	+ } catch (HepMC3::WeightError & w) {
	+ return vector<string>();
	+ }
	+ }
	+
	+ std::valarray<double> weights(const GenEvent & ge) {
	+ return std::valarray<double>(&ge.weights()[0], ge.weights().size());
	+
	+ }

	}
	}
	diff --git a/test/testApi.cc b/test/testApi.cc
	--- a/test/testApi.cc
	+++ b/test/testApi.cc
	@@ -1,34 +1,35 @@
	#include "Rivet/AnalysisHandler.hh"
	#include "HepMC/GenEvent.h"
	#include "Rivet/Tools/RivetHepMC.hh"
	+#include <fstream>

	using namespace std;

	int main(int argc, char* argv[]) {
	assert(argc > 1);

	Rivet::AnalysisHandler ah;
	Rivet::Log::setLevel("Rivet", Rivet::Log::DEBUG);

	// Specify the analyses to be used
	ah.addAnalysis("EXAMPLE");
	ah.addAnalyses({{ "MC_JETS", "EXAMPLE_CUTS", "EXAMPLE_SMEAR" }});

	std::ifstream file("testApi.hepmc");
	shared_ptr<Rivet::HepMC_IO_type> reader = Rivet::HepMCUtils::makeReader(file);
	std::shared_ptr<Rivet::GenEvent> evt = make_shared<Rivet::GenEvent>();
	double sum_of_weights = 0.0;
	while ( Rivet::HepMCUtils::readEvent(reader, evt) ) {
	// Analyse current event
	ah.analyze(*evt);
	sum_of_weights += evt->weights()[0];
	}
	file.close();

	- ah.setCrossSection(1.0, 0.1);
	+ ah.setCrossSection(make_pair(1.0, 0.1));

	ah.finalize();
	ah.writeData("out.yoda");

	return 0;
	}
	diff --git a/test/testNaN.cc b/test/testNaN.cc
	--- a/test/testNaN.cc
	+++ b/test/testNaN.cc
	@@ -1,74 +1,75 @@
	#include "Rivet/AnalysisHandler.hh"
	#include "Rivet/Analysis.hh"
	#include "Rivet/Tools/RivetYODA.hh"
	#include <limits>
	#include <cmath>
	#include <iostream>
	+#include <fstream>

	using namespace std;


	class Test : public Rivet::Analysis {
	public:
	Test() : Analysis("Test") {}

	void init() {
	book(_h_test, "test", 50, 66.0, 116.0);
	}

	void analyze(const Rivet::Event & e) {
	cout << "Normal fill" << endl;
	_h_test->fill(90., 1.);

	cout << "Underflow fill" << endl;
	_h_test->fill(30.,1.);

	cout << "Overflow fill" << endl;
	_h_test->fill(130.,1.);

	cout << "Inf fill" << endl;
	try {
	_h_test->fill(numeric_limits<double>::infinity(), 1.);
	} catch (YODA::RangeError & e) {
	cerr << e.what() << '\n';
	if ( string(e.what()) != string("X is Inf") ) throw;
	}

	cout << "NaN fill" << endl;
	try {
	_h_test->fill(numeric_limits<double>::quiet_NaN(), 1.);
	} catch (YODA::RangeError & e) {
	cerr << e.what() << '\n';
	if ( string(e.what()) != string("X is NaN") ) throw;
	}
	}

	private:
	Rivet::Histo1DPtr _h_test;
	};

	DECLARE_RIVET_PLUGIN(Test);

	int main(int argc, char* argv[]) {
	assert(argc > 1);

	Rivet::AnalysisHandler rivet;
	rivet.addAnalysis("Test");

	std::ifstream file("testApi.hepmc");
	shared_ptr<Rivet::HepMC_IO_type> reader = Rivet::HepMCUtils::makeReader(file);
	std::shared_ptr<Rivet::GenEvent> evt = make_shared<Rivet::GenEvent>();
	double sum_of_weights = 0.0;
	while ( Rivet::HepMCUtils::readEvent(reader, evt) ) {
	// Analyse current event
	rivet.analyze(*evt);
	sum_of_weights += evt->weights()[0];
	}
	file.close();

	- rivet.setCrossSection(1.0, 0.1);
	+ rivet.setCrossSection(make_pair(1.0, 0.1));
	rivet.finalize();
	rivet.writeData("NaN.aida");

	return 0;
	}

File Metadata

Mime Type: text/x-diff
Expires: Sun, Feb 23, 2:47 PM (1 d, 7 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4486706
Default Alt Text: (62 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions