No OneTemporary
Actions

Size

45 KB

Subscribers

None

View Options

	diff --git a/include/Rivet/AnalysisHandler.hh b/include/Rivet/AnalysisHandler.hh
	--- a/include/Rivet/AnalysisHandler.hh
	+++ b/include/Rivet/AnalysisHandler.hh
	@@ -1,256 +1,256 @@
	// -- 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 {
	return _eventCounter->sumW();
	}

	- int numWeights() const {
	+ size_t numWeights() const {
	return _weightNames.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<YODA::AnalysisObjectPtr> getData() const;
	std::vector<reference_wrapper<MultiweightAOPtr> > getRivetAOs() const;
	std::vector<YODA::AnalysisObjectPtr> getYodaAOs() 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::valarray<double> > _subEventWeights;
	size_t _numWeightTypes; // always == WeightVector.size()

	/// Run name
	std::string _runname;

	mutable CounterPtr _eventCounter;

	/// 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/Tools/RivetYODA.hh b/include/Rivet/Tools/RivetYODA.hh
	--- a/include/Rivet/Tools/RivetYODA.hh
	+++ b/include/Rivet/Tools/RivetYODA.hh
	@@ -1,386 +1,386 @@
	#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::Scatter1D> Scatter1DPtr;
	typedef std::shared_ptr<YODA::Scatter2D> Scatter2DPtr;
	typedef std::shared_ptr<YODA::Scatter3D> Scatter3DPtr;
	}



	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
	/// these need to be multi-weighted eventually.
	class Scatter1DPtr : public AnalysisObjectPtr {
	public:
	Scatter1DPtr() :
	_scatter(YODA::Scatter1DPtr()) { }

	Scatter1DPtr(const YODA::Scatter1D& p) :
	_scatter(make_shared<YODA::Scatter1D>(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(const YODA::Scatter2D& p) :
	_scatter(make_shared<YODA::Scatter2D>(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(const YODA::Scatter3D& p) :
	_scatter(make_shared<YODA::Scatter3D>(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 {

	public:
	virtual void newSubEvent() = 0;

	/// @todo
	/// rename to setActive(Idx)?
	virtual void setActiveWeightIdx(unsigned int iWeight) = 0;

	virtual void pushToPersistent(const vector<valarray<double> >& weight) = 0;

	virtual YODA::AnalysisObjectPtr activeYODAPtr() const = 0;
	};


	using Weight = double;

	template <class T>
	using Fill = pair<typename T::FillType, Weight>;

	template <class T>
	using Fills = multiset<Fill<T>>;


	// TODO TODO TODO
	// need to override the old fill method too!
	// otherwise we can't intercept existing fill calls in analysis code
	//
	// need to
	// TODO TODO TODO

	template <class T>
	class TupleWrapper;

	template<>
	class TupleWrapper<YODA::Counter> : public YODA::Counter {
	public:
	typedef shared_ptr<TupleWrapper<YODA::Counter>> Ptr;
	TupleWrapper(const YODA::Counter & h) : YODA::Counter(h) {}
	// todo: do we need to deal with users using fractions directly?
	void fill( double weight=1.0, double fraction=1.0 ) {
	fills_.insert( {YODA::Counter::FillType(),weight} );
	}
	void reset() { fills_.clear(); }
	const Fills<YODA::Counter> & fills() const { return fills_; }
	private:
	// x / weight pairs
	Fills<YODA::Counter> fills_;
	};

	template<>
	class TupleWrapper<YODA::Histo1D> : public YODA::Histo1D {
	public:
	typedef shared_ptr<TupleWrapper<YODA::Histo1D>> Ptr;
	TupleWrapper(const YODA::Histo1D & h) : YODA::Histo1D(h) {}
	// todo: do we need to deal with users using fractions directly?
	void fill( double x, double weight=1.0, double fraction=1.0 ) {
	if ( std::isnan(x) ) throw YODA::RangeError("X is NaN");
	fills_.insert( { x , weight } );
	}
	void reset() { fills_.clear(); }
	const Fills<YODA::Histo1D> & fills() const { return fills_; }
	private:
	// x / weight pairs
	Fills<YODA::Histo1D> fills_;
	};

	template<>
	class TupleWrapper<YODA::Profile1D> : public YODA::Profile1D {
	public:
	typedef shared_ptr<TupleWrapper<YODA::Profile1D>> Ptr;
	TupleWrapper(const YODA::Profile1D & h) : YODA::Profile1D(h) {}
	// todo: do we need to deal with users using fractions directly?
	void fill( double x, double y, double weight=1.0, double fraction=1.0 ) {
	if ( std::isnan(x) ) throw YODA::RangeError("X is NaN");
	if ( std::isnan(y) ) throw YODA::RangeError("Y is NaN");
	- fills_.insert( { {x,y}, weight } );
	+ fills_.insert( { YODA::Profile1D::FillType{x,y}, weight } );
	}
	void reset() { fills_.clear(); }
	const Fills<YODA::Profile1D> & fills() const { return fills_; }
	private:
	// x / weight pairs
	Fills<YODA::Profile1D> fills_;
	};


	template<>
	class TupleWrapper<YODA::Histo2D> : public YODA::Histo2D {
	public:
	typedef shared_ptr<TupleWrapper<YODA::Histo2D>> Ptr;
	TupleWrapper(const YODA::Histo2D & h) : YODA::Histo2D(h) {}
	// todo: do we need to deal with users using fractions directly?
	void fill( double x, double y, double weight=1.0, double fraction=1.0 ) {
	if ( std::isnan(x) ) throw YODA::RangeError("X is NaN");
	if ( std::isnan(y) ) throw YODA::RangeError("Y is NaN");
	- fills_.insert( {{x,y}, weight} );
	+ fills_.insert( { YODA::Histo2D::FillType{x,y}, weight } );
	}
	void reset() { fills_.clear(); }
	const Fills<YODA::Histo2D> & fills() const { return fills_; }
	private:
	// x / weight pairs
	Fills<YODA::Histo2D> fills_;
	};

	template<>
	class TupleWrapper<YODA::Profile2D> : public YODA::Profile2D {
	public:
	typedef shared_ptr<TupleWrapper<YODA::Profile2D>> Ptr;
	TupleWrapper(const YODA::Profile2D & h) : YODA::Profile2D(h) {}
	// todo: do we need to deal with users using fractions directly?
	void fill( double x, double y, double z, double weight=1.0, double fraction=1.0 ) {
	if ( std::isnan(x) ) throw YODA::RangeError("X is NaN");
	if ( std::isnan(y) ) throw YODA::RangeError("Y is NaN");
	if ( std::isnan(z) ) throw YODA::RangeError("Z is NaN");
	- fills_.insert( {{x,y,z}, weight} );
	+ fills_.insert( { YODA::Profile2D::FillType{x,y,z}, weight } );
	}
	void reset() { fills_.clear(); }
	const Fills<YODA::Profile2D> & fills() const { return fills_; }
	private:
	// x / weight pairs
	Fills<YODA::Profile2D> fills_;
	};

	template <class T>
	class Wrapper : 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.
	*/

	Wrapper() = default;

	Wrapper(size_t len_of_weightvec, const T & p);

	typename T::Ptr active() const;

	/* @todo this probably need to loop over all? */
	bool operator!() const { return !active(); }

	operator bool() const { return static_cast<bool>(active()); }

	T * operator->() { return active().get(); }

	T * operator->() const { return active().get(); }

	T & operator() { return active(); }

	const T & 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==(Wrapper a, Wrapper 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!=(Wrapper a, Wrapper b){
	return !(a == b);
	}


	friend bool operator<(Wrapper a, Wrapper 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);
	}

	/* this is for dev only---we shouldn't need this in real runs. */
	void unsetActiveWeight() { _active.reset(); }

	void newSubEvent();

	virtual YODA::AnalysisObjectPtr activeYODAPtr() const { return _active; }

	const vector<typename T::Ptr> & persistent() const { return _persistent; }

	/* to be implemented for each type */
	void pushToPersistent(const vector<valarray<double> >& weight);

	/* M of these, one for each weight */
	vector<typename T::Ptr> _persistent;

	/* N of these, one for each event in evgroup */
	vector<typename TupleWrapper<T>::Ptr> _evgroup;

	typename T::Ptr _active;

	friend class AnalysisHandler;
	};


	// every object listed here needs a virtual fill method in YODA,
	// otherwise the Tuple fakery won't work.

	using Histo1DPtr = Wrapper<YODA::Histo1D>;
	using Histo2DPtr = Wrapper<YODA::Histo2D>;
	using Profile1DPtr = Wrapper<YODA::Profile1D>;
	using Profile2DPtr = Wrapper<YODA::Profile2D>;
	using CounterPtr = Wrapper<YODA::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);

	}

	#endif
	diff --git a/src/Core/AnalysisHandler.cc b/src/Core/AnalysisHandler.cc
	--- a/src/Core/AnalysisHandler.cc
	+++ b/src/Core/AnalysisHandler.cc
	@@ -1,457 +1,457 @@
	// -- 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 <regex>

	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 {


	AnalysisHandler::AnalysisHandler(const string& runname)
	: _runname(runname),
	_eventCounter(0, Counter()), _xs(NAN),
	_initialised(false), _ignoreBeams(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
	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() {
	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!");

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

	setWeightNames(ge);
	if (haveNamedWeights())
	MSG_INFO("Using named weights");
	else
	MSG_INFO("NOT using named weights. Using first weight as nominal weight");

	_numWeightTypes = _weightNames.size();
	_eventCounter = CounterPtr(_numWeightTypes, Counter("_EVTCOUNT"));

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

	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) {
	MSG_DEBUG("Name of weight #" << _weightNames.size() << ": " << temp[0]);
	_weightNames.push_back(temp[0]);
	}
	}
	}


	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
	for (const AnaHandle& a : _analyses) {
	for (const auto & ao : a->analysisObjects()) {
	ao.get().pushToPersistent(_subEventWeights);
	}
	}

	_eventCounter.pushToPersistent(_subEventWeights);

	_eventNumber = ge.event_number();

	MSG_DEBUG("Event #" << numEvents() << " weight sum = " << sumOfWeights());
	MSG_DEBUG("Event has " << _subEventWeights.size() << " sub events.");
	_subEventWeights.clear();
	}

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

	_eventCounter.newSubEvent();

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


	// 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<reference_wrapper<MultiweightAOPtr> > AnalysisHandler::getRivetAOs() const {
	vector<reference_wrapper<MultiweightAOPtr> > rtn;

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

	// Should event counter be included here?
	rtn.push_back(_eventCounter);

	return rtn;
	}

	vector<YODA::AnalysisObjectPtr> AnalysisHandler::getYodaAOs() const {
	vector<YODA::AnalysisObjectPtr> rtn;

	const auto & raos = getRivetAOs();

	for (const auto & raoref : raos) {
	MultiweightAOPtr & rao = raoref.get();
	if (rao->path().find("/TMP/") != string::npos)
	continue;

	for (size_t iW = 0; iW < numWeights(); iW++) {
	rao.setActiveWeightIdx(iW);

	// add the weight name in brackets unless we recognize a
	// nominal weight
	if (_weightNames[iW] != "Weight"
	&& _weightNames[iW] != "0"
	&& _weightNames[iW] != "Default")
	rao->setPath(rao->path() + "[" + _weightNames[iW] + "]");

	rtn.push_back(rao.activeYODAPtr());
	}
	}

	YODA::Scatter1D::Points pts; pts.insert(YODA::Point1D(_xs, _xserr));
	rtn.push_back( make_shared<Scatter1D>(pts, "/_XSEC") );

	sort(rtn.begin(), rtn.end(),
	[](YODA::AnalysisObjectPtr a, YODA::AnalysisObjectPtr b) {
	return a->path() < b->path();
	}
	);

	return rtn;
	}

	vector<YODA::AnalysisObjectPtr> AnalysisHandler::getData() const {

	return getYodaAOs();
	}


	void AnalysisHandler::writeData(const string& filename) const {
	const vector<YODA::AnalysisObjectPtr> aos = getData();
	try {
	YODA::WriterYODA::write(filename, aos.begin(), aos.end());
	} catch ( YODA::WriteError ) {
	throw UserError("Unexpected error in writing file to: " + filename);
	}
	}


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


	/*
	* why is this here?
	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/Tools/RivetYODA.cc b/src/Tools/RivetYODA.cc
	--- a/src/Tools/RivetYODA.cc
	+++ b/src/Tools/RivetYODA.cc
	@@ -1,335 +1,335 @@
	#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 {


	template <class T>
	Wrapper<T>::Wrapper(size_t len_of_weightvec, const T & p) {
	for (size_t m = 0; m < len_of_weightvec; ++m)
	_persistent.push_back(make_shared<T>(p));
	}


	template <class T>
	typename T::Ptr Wrapper<T>::active() const {
	assert(_active);
	return _active;
	}


	template <class T>
	void Wrapper<T>::newSubEvent() {
	typename TupleWrapper<T>::Ptr tmp
	= make_shared<TupleWrapper<T>>(_persistent[0]->clone());
	tmp->reset();
	_evgroup.push_back( tmp );
	_active = _evgroup.back();
	assert(_active);
	}



	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, 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, YODA::AnalysisObjectPtr> rtn;
	foreach ( YODA::AnalysisObject* ao, aovec ) {
	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;
	}
	}

	namespace {

	using Rivet::Fill;
	using Rivet::Fills;
	using Rivet::TupleWrapper;

	template <class T>
	double get_window_size(const typename T::Ptr & histo,
	typename T::BinType x) {
	// the bin index we fall in
	const auto binidx = histo->binIndexAt(x);
	// gaps, overflow, underflow don't contribute
	if ( binidx == -1 )
	return 0;


	const auto & b = histo->bin(binidx);

	// if we don't have a valid neighbouring bin,
	// we use infinite width
	typename T::Bin b1(-1.0/0.0, 1.0/0.0);

	// points in the top half compare to the upper neighbour
	if ( x > b.xMid() ) {
	- int nextidx = binidx + 1;
	+ size_t nextidx = binidx + 1;
	if ( nextidx < histo->bins().size() )
	b1 = histo->bin(nextidx);
	}
	else { // compare to the lower neighbour
	int nextidx = binidx - 1;
	if ( nextidx >= 0 )
	b1 = histo->bin(nextidx);
	}
	// the factor 2 is arbitrary, could poss. be smaller
	return min( b.width(), b1.width() ) / 2.0;
	}

	template <class T>
	typename T::BinType
	fillT2binT(typename T::FillType a) {
	return a;
	}

	template <>
	YODA::Profile1D::BinType
	fillT2binT<YODA::Profile1D>(YODA::Profile1D::FillType a) {
	return get<0>(a);
	}

	template <>
	YODA::Profile2D::BinType
	fillT2binT<YODA::Profile2D>(YODA::Profile2D::FillType a) {
	- return { get<0>(a), get<1>(a) };
	+ return YODA::Profile2D::BinType{ get<0>(a), get<1>(a) };
	}



	template <class T>
	void commit(vector<typename T::Ptr> & persistent,
	const vector< vector<Fill<T>> > & tuple,
	const vector<valarray<double>> & weights ) {

	// TODO check if all the xs are in the same bin anyway!
	// Then no windowing needed

	assert(persistent.size() == weights[0].size());

	for ( const auto & x : tuple ) {
	double maxwindow = 0.0;
	for ( const auto & xi : x ) {
	// TODO check for NOFILL here
	// persistent[0] has the same binning as all the other persistent objects
	double window = get_window_size<T>(persistent[0], fillT2binT<T>(xi.first));
	if ( window > maxwindow )
	maxwindow = window;
	}
	const double wsize = maxwindow;
	// all windows have same size

	set<double> edgeset;
	// bin edges need to be in here!
	for ( const auto & xi : x ) {
	edgeset.insert(fillT2binT<T>(xi.first) - wsize);
	edgeset.insert(fillT2binT<T>(xi.first) + wsize);
	}

	vector< std::tuple<double,valarray<double>,double> > hfill;
	double sumf = 0.0;
	auto edgit = edgeset.begin();
	double ehi = *edgit;
	while ( ++edgit != edgeset.end() ) {
	double elo = ehi;
	ehi = *edgit;
	valarray<double> sumw(0.0, persistent.size()); // need m copies of this
	bool gap = true; // Check for gaps between the sub-windows.
	- for ( int i = 0; i < x.size(); ++i ) {
	+ for ( size_t i = 0; i < x.size(); ++i ) {
	// check equals comparisons here!
	if ( fillT2binT<T>(x[i].first) + wsize >= ehi
	&&
	fillT2binT<T>(x[i].first) - wsize <= elo ) {
	sumw += x[i].second * weights[i];
	gap = false;
	}
	}
	if ( gap ) continue;
	- hfill.push_back( { (ehi + elo)/2.0, sumw, (ehi - elo) } );
	+ hfill.push_back( make_tuple( (ehi + elo)/2.0, sumw, (ehi - elo) ) );
	sumf += ehi - elo;
	}

	for ( auto f : hfill )
	for ( size_t m = 0; m < persistent.size(); ++m )
	persistent[m]->fill( get<0>(f), get<1>(f)[m], get<2>(f)/sumf );
	// Note the scaling to one single fill

	}
	}

	template<>
	void commit<YODA::Histo2D>(vector<YODA::Histo2D::Ptr> & persistent,
	const vector< vector<Fill<YODA::Histo2D>> > & tuple,
	const vector<valarray<double>> & weights)
	{}

	template<>
	void commit<YODA::Profile2D>(vector<YODA::Profile2D::Ptr> & persistent,
	const vector< vector<Fill<YODA::Profile2D>> > & tuple,
	const vector<valarray<double>> & weights)
	{}

	template <class T>
	double distance(T a, T b) {
	return abs(a - b);
	}

	template <>
	double distance<tuple<double,double> >(tuple<double,double> a, tuple<double,double> b) {
	return Rivet::sqr(get<0>(a) - get<0>(b)) + Rivet::sqr(get<1>(a) - get<1>(b));
	}



	}




	/// fills is a vector of sub-event with an ordered set of x-values of
	/// the fills in each sub-event. NOFILL should be an "impossible"
	/// value for this histogram. Returns a vector of sub-events with
	/// an ordered vector of fills (including NOFILLs) for each sub-event.
	template <class T>
	vector< vector<Fill<T> > >
	match_fills(const vector<typename TupleWrapper<T>::Ptr> & evgroup, const Fill<T> & NOFILL)
	{
	vector< vector<Fill<T> > > matched;
	// First just copy subevents into vectors and find the longest vector.
	unsigned int maxfill = 0; // length of biggest vector
	int imax = 0; // index position of biggest vector
	for ( const auto & it : evgroup ) {
	const auto & subev = it->fills();
	if ( subev.size() > maxfill ) {
	maxfill = subev.size();
	imax = matched.size();
	}
	matched.push_back(vector<Fill<T> >(subev.begin(), subev.end()));
	}
	// Now, go through all subevents with missing fills.
	const vector<Fill<T>> & full = matched[imax]; // the longest one
	for ( auto & subev : matched ) {
	if ( subev.size() == maxfill ) continue;

	// Add NOFILLs to the end;
	while ( subev.size() < maxfill ) subev.push_back(NOFILL);

	// Iterate from the back and shift all fill values backwards by
	// swapping with NOFILLs so that they better match the full
	// subevent.
	for ( int i = maxfill - 1; i >= 0; --i ) {
	if ( subev[i] == NOFILL ) continue;
	- int j = i;
	+ size_t j = i;
	while ( j + 1 < maxfill && subev[j + 1] == NOFILL &&
	distance(fillT2binT<T>(subev[j].first),
	fillT2binT<T>(full[j].first))
	>
	distance(fillT2binT<T>(subev[j].first),
	fillT2binT<T>(full[j + 1].first)) )
	{
	swap(subev[j], subev[j + 1]);
	++j;
	}
	}
	}
	// transpose
	vector<vector<Fill<T>>> result(maxfill,vector<Fill<T>>(matched.size()));
	for (size_t i = 0; i < matched.size(); ++i)
	for (size_t j = 0; j < maxfill; ++j)
	result.at(j).at(i) = matched.at(i).at(j);
	return result;
	}






	namespace Rivet {

	template <class T>
	void Wrapper<T>::pushToPersistent(const vector<valarray<double> >& weight) {

	// have we had subevents at all?
	const bool have_subevents = _evgroup.size() > 1;
	if ( ! have_subevents ) {



	assert( _evgroup.size() == 1 && weight.size() == 1 );
	// simple replay of all tuple entries
	// each recorded fill is inserted into all persistent weightname histos
	for ( size_t m = 0; m < _persistent.size(); ++m )
	for ( const auto & f : _evgroup[0]->fills() )
	_persistent[m]->fill( f.first, f.second * weight[0][m] );



	} else {



	assert( _evgroup.size() == weight.size() );
	// outer index is subevent, inner index is jets in the event
	vector<vector<Fill<T>>> linedUpXs
	= match_fills<T>(_evgroup, {typename T::FillType(), 0.0});
	commit<T>( _persistent, linedUpXs, weight );



	}
	_evgroup.clear();
	_active.reset();
	}

	template <>
	void Wrapper<YODA::Counter>::pushToPersistent(const vector<valarray<double> >& weight) {
	for ( size_t m = 0; m < _persistent.size(); ++m ) {
	for ( size_t n = 0; n < _evgroup.size(); ++n ) {
	for ( const auto & f : _evgroup[n]->fills() ) {
	_persistent[m]->fill( f.second * weight[n][m] );
	}
	}
	}
	}

	// explicitly instantiate all wrappers

	template class Wrapper<YODA::Histo1D>;
	template class Wrapper<YODA::Histo2D>;
	template class Wrapper<YODA::Profile1D>;
	template class Wrapper<YODA::Profile2D>;
	template class Wrapper<YODA::Counter>;

	}

File Metadata

Mime Type: text/x-diff
Expires: Sun, Feb 23, 2:58 PM (5 h, 46 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4486749
Default Alt Text: (45 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions