No OneTemporary
Actions

Size

46 KB

Subscribers

None

View Options

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

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

	-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),
	_initialised(false), _ignoreBeams(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 (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);
	}

	// 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) {
	_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);
	+ }
	+
	// 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.");
	}

	_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::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);
	_eventCounter.get()->pushToFinal();
	_xs.get()->pushToFinal();
	for (const AnaHandle& a : analyses()) {
	for (auto ao : a->analysisObjects()) {
	ao.get()->pushToPersistent(_subEventWeights);
	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
	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 = 0; iW < numWeights(); iW++) {
	+ 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 = 0; iW < numWeights(); iW++) {
	+ for (size_t iW : order ) {
	for ( auto rao : raos ) {
	- rao.get()->setActiveFinalWeightIdx(iW);
	+ 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) {
	_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);
	}

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

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


	}
	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,161 +1,174 @@
	// -- 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;
	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 = Rivet::split(m.str(), "[,]");
	+ 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();
	}

	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/RivetYODA.cc b/src/Tools/RivetYODA.cc
	--- a/src/Tools/RivetYODA.cc
	+++ b/src/Tools/RivetYODA.cc
	@@ -1,505 +1,507 @@
	#include "Rivet/Config/RivetCommon.hh"
	#include "Rivet/Tools/RivetYODA.hh"
	#include "Rivet/Tools/RivetPaths.hh"
	#include "YODA/ReaderYODA.h"
	#include "YODA/ReaderAIDA.h"

	// use execinfo for backtrace if available
	#include "DummyConfig.hh"
	#ifdef HAVE_EXECINFO_H
	#include <execinfo.h>
	#endif

	#include <regex>
	#include <sstream>

	using namespace std;

	namespace Rivet {


	template <class T>
	Wrapper<T>::~Wrapper() {}


	template <class T>
	Wrapper<T>::Wrapper(const vector<string>& weightNames, const T & p)
	{
	_basePath = p.path();
	for (const string& weightname : weightNames) {
	_persistent.push_back(make_shared<T>(p));
	_final.push_back(make_shared<T>(p));

	auto obj = _persistent.back();
	auto final = _final.back();
	if (weightname != "") {
	obj->setPath("/RAW" + obj->path() + "[" + weightname + "]");
	final->setPath(final->path() + "[" + weightname + "]");
	}
	}
	}


	template <class T>
	typename T::Ptr Wrapper<T>::active() const {
	if ( !_active ) {
	#ifdef HAVE_BACKTRACE
	void * buffer[4];
	backtrace(buffer, 4);
	backtrace_symbols_fd(buffer, 4 , 1);
	#endif
	assert(false && "No active pointer set. Was this object booked in init()?");
	}
	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" +
	" in data 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;
	for ( 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() ) {
	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 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 ( 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( 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));
	}

	}


	namespace Rivet {

	bool copyao(YODA::AnalysisObjectPtr src, YODA::AnalysisObjectPtr dst) {
	for (const std::string& a : src->annotations())
	dst->setAnnotation(a, src->annotation(a));
	if ( aocopy<Counter>(src,dst) ) return true;
	if ( aocopy<Histo1D>(src,dst) ) return true;
	if ( aocopy<Histo2D>(src,dst) ) return true;
	if ( aocopy<Profile1D>(src,dst) ) return true;
	if ( aocopy<Profile2D>(src,dst) ) return true;
	if ( aocopy<Scatter1D>(src,dst) ) return true;
	if ( aocopy<Scatter2D>(src,dst) ) return true;
	if ( aocopy<Scatter3D>(src,dst) ) return true;
	return false;
	}

	bool addaos(YODA::AnalysisObjectPtr dst, YODA::AnalysisObjectPtr src, double scale) {
	if ( aoadd<Counter>(dst,src,scale) ) return true;
	if ( aoadd<Histo1D>(dst,src,scale) ) return true;
	if ( aoadd<Histo2D>(dst,src,scale) ) return true;
	if ( aoadd<Profile1D>(dst,src,scale) ) return true;
	if ( aoadd<Profile2D>(dst,src,scale) ) return true;
	return false;
	}

	}


	namespace {

	/// 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;
	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) {
	assert( _evgroup.size() == weight.size() );

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



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



	// 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 <class T>
	void Wrapper<T>::pushToFinal() {
	for ( size_t m = 0; m < _persistent.size(); ++m ) {
	copyao(_persistent.at(m), _final.at(m));
	+ if ( _final[m]->path().substr(0,4) == "/RAW" )
	+ _final[m]->setPath(_final[m]->path().substr(4));
	}
	}



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

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

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

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

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

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

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

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

	// 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>;
	template class Wrapper<YODA::Scatter1D>;
	template class Wrapper<YODA::Scatter2D>;
	template class Wrapper<YODA::Scatter3D>;

	AOPath::AOPath(string fullpath) {
	// First checck if This is a global system path
	_path = fullpath;
	std::regex resys("^(/RAW)?/([^\\[/]+)(\\[(.+)\\])?$");
	smatch m;
	_valid = regex_search(fullpath, m, resys);
	if ( _valid ) {
	_raw = (m[1] == "/RAW");
	_name = m[2];
	_weight = m[4];
	return;
	}
	// If not, assume it is a normal analysis path.
	std::regex repath("^(/RAW)?(/REF)?/([^/:]+)(:[^/]+)?(/TMP)?/([^\\[]+)(\\[(.+)\\])?");
	_valid = regex_search(fullpath, m, repath);
	if ( !_valid ) return;
	_raw = (m[1] == "/RAW");
	_ref = (m[2] == "/REF");
	_analysis = m[3];
	_optionstring = m[4];
	_tmp = (m[5] == "/TMP");
	_name = m[6];
	_weight = m[8];
	std::regex reopt(":([^=]+)=([^:]+)");
	string s = _optionstring;
	while ( regex_search(s, m, reopt) ) {
	_options[m[1]] = m[2];
	s = m.suffix();
	}
	}

	void AOPath::fixOptionString() {
	ostringstream oss;
	for ( auto optval : _options )
	oss << ":" << optval.first << "=" << optval.second;
	_optionstring = oss.str();
	}

	string AOPath::mkPath() const {
	ostringstream oss;
	if ( isRaw() ) oss << "/RAW";
	else if ( isRef() ) oss << "/REF";
	if ( _analysis != "" ) oss << "/" << analysis();
	for ( auto optval : _options )
	oss << ":" << optval.first << "=" << optval.second;
	if ( isTmp() ) oss << "/TMP";
	oss << "/" << name();
	if ( weight() != "" )
	oss << "[" << weight() << "]";
	return oss.str();
	}

	void AOPath::debug() const {
	cout << "Full path: " << _path << endl;
	if ( !_valid ) {
	cout << "This is not a valid analysis object path" << endl << endl;
	return;
	}
	cout << "Check path: " << mkPath() << endl;
	cout << "Analysis: " << _analysis << endl;
	cout << "Name: " << _name << endl;
	cout << "Weight: " << _weight << endl;
	cout << "Properties: ";
	if ( _raw ) cout << "raw ";
	if ( _tmp ) cout << "tmp ";
	if ( _ref ) cout << "ref ";
	cout << endl;
	cout << "Options: ";
	for ( auto opt : _options )
	cout << opt.first << "->" << opt.second << " ";
	cout << endl << endl;
	}

	}

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 6:28 PM (1 d, 22 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3805593
Default Alt Text: (46 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions