No OneTemporary
Actions

Size

29 KB

Subscribers

None

View Options

	diff --git a/bin/rivet-hepmz.cc b/bin/rivet-hepmz.cc
	--- a/bin/rivet-hepmz.cc
	+++ b/bin/rivet-hepmz.cc
	@@ -1,119 +1,120 @@
	#include "Rivet/Tools/RivetHepMC.hh"
	#include "Rivet/Tools/WriterCompressedAscii.hh"
	#include "../src/Core/zstr/zstr.hpp"
	#include "HepMC3/GenParticle.h"
	#include "HepMC3/GenVertex.h"
	#include "HepMC3/WriterAscii.h"
	using namespace std;

	int main(int argc, char** argv) {

	string ifile;
	string ofile;
	bool etaphi = false;
	bool strip = false;
	bool userivet = false;
	bool help = false;
	double pphi = 0.0;
	double peta = 0.0;
	double pe = 0.0;
	double pm = 0.0;
	double po = 0.0;

	for ( int iarg = 1; iarg < argc; ++iarg ) {
	string arg = argv[iarg];
	if ( arg == "-s" )
	userivet = strip = etaphi = true;
	else if ( arg == "-S" )
	userivet = strip = true;
	else if ( arg == "-c" )
	userivet = etaphi = true;
	else if ( arg == "-r" )
	userivet = true;
	else if ( arg == "-p" ) {
	if ( iarg + 1 >= argc ) {
	help = true;
	break;
	}
	arg = argv[++iarg];
	if ( arg.substr(0,4) == "phi=" ) pphi = stod(arg.substr(4));
	if ( arg.substr(0,4) == "eta=" ) peta = stod(arg.substr(4));
	if ( arg.substr(0,2) == "e=" ) pe = stod(arg.substr(2));
	if ( arg.substr(0,2) == "m=" ) pm = stod(arg.substr(2));
	if ( arg.substr(0,2) == "p=" ) po = stod(arg.substr(2));
	userivet = etaphi = true;
	}
	else if ( arg == "-h" )
	help = true;
	else {
	if ( ifile.empty() ) ifile = arg;
	else if ( ofile.empty() ) ofile = arg;
	else {
	cout << "Unknown argument '" << arg << "'" << endl;
	help = true;
	}
	}
	}
	if ( ofile.empty() ) help = true;
	if ( help ) {
	cout << "Usage: " << argv[0]
	<< " [options] <input-hepmcfile> <output-hepmcfile>" << endl;
	cout << " where options are one or more of" << endl
	<< " -r: rivet internal hepmc output" << endl
	<< " -c: compressed hepmc output (implies -r)" << endl
	<< " -s: strips hepmc from unobservable (some) particles (implies -c)" << endl
	<< " -p type=prec: precision in compressed hepmc output (implies -c)" << endl
	<< " allowed types are phi, eta, e and m. " << endl
	<< " -h: write this message and exit." << endl;

	return 1;
	}


	- Rivet::zstr::ifstream input(ifile);
	+
	+ std::shared_ptr<std::istream> istr;
	shared_ptr<ostream> output;
	if ( ofile.substr(ofile.length() - 3) == ".gz" \|\|
	ofile.substr(ofile.length() - 6) == ".hepmz" )
	output = make_shared<Rivet::zstr::ofstream>(ofile);
	else
	output = make_shared<ofstream>(ofile);

	shared_ptr<Rivet::HepMC_IO_type>
	- reader = Rivet::HepMCUtils::makeReader(input);
	+ reader = Rivet::HepMCUtils::makeReader(ifile, istr);

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

	shared_ptr<HepMC3::Writer> writer;
	if ( userivet ) {
	auto compressed = make_shared<Rivet::WriterCompressedAscii>(*output);
	if ( etaphi ) {
	compressed->use_integers();
	if ( pphi > 0.0 ) compressed->set_precision_phi(pphi);
	if ( peta > 0.0 ) compressed->set_precision_eta(peta);
	if ( pe > 0.0 ) compressed->set_precision_e(pe);
	if ( pm > 0.0 ) compressed->set_precision_m(pm);
	if ( po > 0.0 ) compressed->set_precision(po);
	}
	if ( strip ) {
	compressed->add_stripid(21);
	compressed->add_stripid(-1);
	compressed->add_stripid(1);
	compressed->add_stripid(-2);
	compressed->add_stripid(2);
	compressed->add_stripid(-3);
	compressed->add_stripid(3);
	}
	writer = compressed;
	} else {
	writer = make_shared<HepMC3::WriterAscii>(*output);
	}

	while(reader && Rivet::HepMCUtils::readEvent(reader, evt) ) {
	writer->write_event(*evt);
	}

	return 0;
	}


	diff --git a/bin/rivet-nopy.cc b/bin/rivet-nopy.cc
	--- a/bin/rivet-nopy.cc
	+++ b/bin/rivet-nopy.cc
	@@ -1,56 +1,41 @@
	#include <fstream>
	-#ifdef ENABLE_HEPMC_3
	-#include "HepMC3/GenEvent.h"
	-#include "HepMC3/GenParticle.h"
	-#include "HepMC3/GenVertex.h"
	-#include "HepMC3/GenVertex.h"
	-#include "HepMC3/ReaderFactory.h"
	-#endif
	#include "Rivet/Tools/RivetHepMC.hh"
	#include "Rivet/AnalysisHandler.hh"
	#include "Rivet/AnalysisLoader.hh"

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

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

	-#ifdef ENABLE_HEPMC_3
	- std::shared_ptr<HepMC3::Reader> reader = HepMC3::deduce_reader(argv[1]);
	- std::shared_ptr<HepMC3::GenEvent> evt = std::make_shared<HepMC3::GenEvent>();
	-
	+ std::shared_ptr<std::istream> istr;

	-
	-#else
	- std::ifstream istr(argv[1], std::ios::in);
	-
	- std::shared_ptr<Rivet::HepMC_IO_type> reader = Rivet::HepMCUtils::makeReader(istr);
	+ std::shared_ptr<Rivet::HepMC_IO_type> reader = Rivet::HepMCUtils::makeReader(argv[1], istr);

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

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




	ah.setCrossSection(std::make_pair(1.0, 0.0));
	ah.finalize();
	ah.writeData("Rivet.yoda");

	return 0;
	}
	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,107 +1,108 @@
	// -- C++ --
	#ifndef RIVET_RivetHepMC_HH
	#define RIVET_RivetHepMC_HH

	#include "Rivet/Config/RivetConfig.hh"

	#include <valarray>

	#ifdef RIVET_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::shared_ptr<HepMC_IO_type> makeReader(std::string filename,
	+ std::shared_ptr<std::istream> &istrp,
	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);
	pair<double,double> crossSection(const GenEvent & ge);
	std::valarray<double> weights(const GenEvent & ge);
	}
	}

	#endif
	diff --git a/src/Core/Run.cc b/src/Core/Run.cc
	--- a/src/Core/Run.cc
	+++ b/src/Core/Run.cc
	@@ -1,184 +1,125 @@
	// -- 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>

	-#ifdef ENABLE_HEPMC_3
	-#include "HepMC3/ReaderFactory.h"
	-#endif
	-
	using std::cout;
	using std::endl;

	namespace Rivet {
	-union magic_t {
	- uint8_t bytes[4];
	- uint32_t number;
	-};

	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 (_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;

	-#ifdef ENABLE_HEPMC_3
	-if (evtfile == "-") {
	- printf("Reading events from standard input assuming IO_GenEvent format\n");
	-#ifdef HAVE_LIBZ
	- _istr = make_shared<zstr::istream>(std::cin);
	- _hepmcReader = std::shared_ptr<RivetHepMC::Reader>((RivetHepMC::Reader) ( new RivetHepMC::ReaderAsciiHepMC2(_istr)));
	-#else
	- _hepmcReader = std::shared_ptr<RivetHepMC::Reader>((RivetHepMC::Reader*) ( new RivetHepMC::ReaderAsciiHepMC2(std::cin)));
	-#endif
	-}
	-else
	-{
	-_hepmcReader = RivetHepMC::deduce_reader(evtfile);
	-#ifdef HAVE_LIBZ
	-if (!_hepmcReader)
	-{
	- printf("No succes with deduction of file type. Test if the file is compressed.\n");
	- std::ifstream file_test(evtfile);
	- magic_t my_magic = {0x1f, 0x8b, 0x08, 0x08};
	- magic_t file_magic;
	- file_test.read((char *) file_magic.bytes, sizeof(file_magic));
	- if ( file_magic.number == my_magic.number )
	- {
	- printf("File is compressed\n");
	- printf("Reading events from compressed file assuming IO_GenEvent format\n");
	- _istr = make_shared<zstr::ifstream>(evtfile);
	- _hepmcReader = std::shared_ptr<RivetHepMC::Reader>((RivetHepMC::Reader) ( new RivetHepMC::ReaderAsciiHepMC2(_istr)));
	- }
	-}
	-#endif
	-}
	+ // Use Rivet's own file format deduction (which uses the one in
	+ // HepMC3 if needed).
	+ _hepmcReader = HepMCUtils::makeReader(evtfile, _istr, &errormessage);

	-#else
	- // 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);
	-
	- }
	-#endif
	-
	+ // Check that it worked.
	if (_hepmcReader == nullptr) {
	Log::getLog("Rivet.Run")
	- << Log::ERROR << "Read error on file '" << evtfile << "' "
	+ << Log::ERROR << "Read error in 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(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,201 +1,214 @@
	// -- C++ --

	//#include <regex>
	#include "Rivet/Tools/Utils.hh"
	#include "Rivet/Tools/RivetHepMC.hh"
	#include "Rivet/Tools/Logging.hh"
	+#include "../Core/zstr/zstr.hpp"

	/*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::shared_ptr<HepMC::IO_GenEvent> makeReader(std::string filename,
	+ std::shared_ptr<std::istream> & istrp,
	std::string *) {
	+#ifdef HAVE_LIBZ
	+ if ( filename == "-" )
	+ istrp = make_shared<Rivet::zstr::istream>(std::cin);
	+ else
	+ istrp = make_shared<Rivet::zstr::ifstream>(filename.c_str());
	+ std::istream & istr = *istrp;
	+#else
	+ if ( filename != "-" ) istrp = make_shared<std::ifstream>(filename.c_str());
	+ std::istream & istr = filename == "-"? std::cin: *istrp;
	+#endif
	+
	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;
	}

	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,185 +1,209 @@
	// -- C++ --

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

	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,
	+ shared_ptr<HepMC_IO_type> makeReader(std::string filename,
	+ std::shared_ptr<std::istream> & istrp,
	std::string * errm) {
	shared_ptr<HepMC_IO_type> ret;
	-
	+
	+#ifdef HAVE_LIBZ
	+ if ( filename == "-" )
	+ istrp = make_shared<zstr::istream>(std::cin);
	+ else
	+ istrp = make_shared<zstr::ifstream>(filename.c_str());
	+ std::istream & istr = *istrp;
	+#else
	+ if ( filename != "-" ) istrp = make_shared<std::ifstream>(filename.c_str());
	+ std::istream & istr = filename == "-"? std::cin: *istrp;
	+#endif
	+
	// 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
	+ else if ( filetype == 2 )
	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>();
	+ if ( ret ) {
	+ if ( ret->failed() ) {
	+ if ( errm ) *errm = "Problems reading from HepMC file. ";
	+ ret = shared_ptr<HepMC_IO_type>();
	+ }
	+ return ret;
	+ }
	+ if ( !ret && filename == "-" ) {
	+ if ( errm ) *errm += "Problems reading HepMC from stdin. No header found. ";
	+ return shared_ptr<HepMC_IO_type>();
	}

	+ // Now we try to reopen the file and see if we can read something.
	+ if ( errm ) *errm += "Could not deduce file format. Will ask HepMC3 to try. ";
	+ ret = RivetHepMC::deduce_reader(filename);
	+
	+
	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) {
	vector<string> ret;
	try {
	ret = ge.weight_names("");
	} catch (HepMC3::WeightError & w) {
	return vector<string>();
	}
	for ( string & w : ret )
	if ( w == "Weight" \|\| w == "0" \|\| w == "Default" )
	w = "";
	return ret;

	}

	std::valarray<double> weights(const GenEvent & ge) {
	return std::valarray<double>(&ge.weights()[0], ge.weights().size());

	}

	}
	}

File Metadata

Mime Type: text/x-diff
Expires: Sun, Feb 23, 1:59 PM (2 h, 23 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4486423
Default Alt Text: (29 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions