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] " << 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 istr; shared_ptr output; if ( ofile.substr(ofile.length() - 3) == ".gz" || ofile.substr(ofile.length() - 6) == ".hepmz" ) output = make_shared(ofile); else output = make_shared(ofile); shared_ptr - reader = Rivet::HepMCUtils::makeReader(input); + reader = Rivet::HepMCUtils::makeReader(ifile, istr); shared_ptr evt = make_shared(); shared_ptr writer; if ( userivet ) { auto compressed = make_shared(*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(*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 -#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] << " [ ...]" << '\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 reader = HepMC3::deduce_reader(argv[1]); - std::shared_ptr evt = std::make_shared(); - + std::shared_ptr istr; - -#else - std::ifstream istr(argv[1], std::ios::in); - - std::shared_ptr reader = Rivet::HepMCUtils::makeReader(istr); + std::shared_ptr reader = Rivet::HepMCUtils::makeReader(argv[1], istr); std::shared_ptr evt = std::make_shared(); -#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 #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; /// @todo Use mcutils? namespace HepMCUtils{ ConstGenParticlePtr getParticlePtr(const RivetHepMC::GenParticle & gp); std::vector particles(ConstGenEventPtr ge); std::vector particles(const GenEvent *ge); std::vector vertices(ConstGenEventPtr ge); std::vector vertices(const GenEvent *ge); std::vector particles(ConstGenVertexPtr gv, const Relatives &relo); std::vector particles(ConstGenParticlePtr gp, const Relatives &relo); int uniqueId(ConstGenParticlePtr gp); int particles_size(ConstGenEventPtr ge); int particles_size(const GenEvent *ge); std::pair beams(const GenEvent *ge); - std::shared_ptr makeReader(std::istream &istr, + std::shared_ptr makeReader(std::string filename, + std::shared_ptr &istrp, std::string * errm = 0); bool readEvent(std::shared_ptr io, std::shared_ptr evt); void strip(GenEvent & ge, const set & stripid = {1, -1, 2, -2, 3,-3, 21}); vector weightNames(const GenEvent & ge); pair crossSection(const GenEvent & ge); std::valarray 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 #include -#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(std::cin); - _hepmcReader = std::shared_ptr((RivetHepMC::Reader*) ( new RivetHepMC::ReaderAsciiHepMC2(*_istr))); -#else - _hepmcReader = std::shared_ptr((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(evtfile); - _hepmcReader = std::shared_ptr((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(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(evtfile.c_str()); - #else - _istr = make_shared(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 #include "Rivet/Tools/Utils.hh" #include "Rivet/Tools/RivetHepMC.hh" #include "Rivet/Tools/Logging.hh" +#include "../Core/zstr/zstr.hpp" /*namespace { inline std::vector 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 particles(ConstGenEventPtr ge){ std::vector result; for(GenEvent::particle_const_iterator pi = ge->particles_begin(); pi != ge->particles_end(); ++pi){ result.push_back(*pi); } return result; } std::vector particles(const GenEvent *ge){ std::vector result; for(GenEvent::particle_const_iterator pi = ge->particles_begin(); pi != ge->particles_end(); ++pi){ result.push_back(*pi); } return result; } std::vector vertices(ConstGenEventPtr ge){ std::vector result; for(GenEvent::vertex_const_iterator vi = ge->vertices_begin(); vi != ge->vertices_end(); ++vi){ result.push_back(*vi); } return result; } std::vector vertices(const GenEvent *ge){ std::vector result; for(GenEvent::vertex_const_iterator vi = ge->vertices_begin(); vi != ge->vertices_end(); ++vi){ result.push_back(*vi); } return result; } std::vector particles(ConstGenVertexPtr gv, const Relatives &relo){ std::vector 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(gv); for (HepMC::GenVertex::particle_iterator pi = gv2->particles_begin(relo); pi != gv2->particles_end(relo); ++pi) result.push_back(const_cast(*pi)); #endif return result; } std::vector 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 beams(const GenEvent *ge){ return ge->beam_particles(); } - std::shared_ptr makeReader(std::istream &istr, + std::shared_ptr makeReader(std::string filename, + std::shared_ptr & istrp, std::string *) { +#ifdef HAVE_LIBZ + if ( filename == "-" ) + istrp = make_shared(std::cin); + else + istrp = make_shared(filename.c_str()); + std::istream & istr = *istrp; +#else + if ( filename != "-" ) istrp = make_shared(filename.c_str()); + std::istream & istr = filename == "-"? std::cin: *istrp; +#endif + return make_shared(istr); } bool readEvent(std::shared_ptr io, std::shared_ptr 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 &) {} vector 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 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 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 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 crossSection(const GenEvent & ge) { return make_pair(ge.cross_section()->cross_section(), ge.cross_section()->cross_section_error()); } std::valarray weights(const GenEvent & ge) { const size_t W = ge.weights().size(); std::valarray 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 +#include "../Core/zstr/zstr.hpp" namespace Rivet{ namespace HepMCUtils{ ConstGenParticlePtr getParticlePtr(const RivetHepMC::GenParticle & gp) { return gp.shared_from_this(); } std::vector particles(ConstGenEventPtr ge){ return ge->particles(); } std::vector particles(const GenEvent *ge){ assert(ge); return ge->particles(); } std::vector vertices(ConstGenEventPtr ge){ return ge->vertices(); } std::vector vertices(const GenEvent *ge){ assert(ge); return ge->vertices(); } std::vector particles(ConstGenVertexPtr gv, const Relatives &relo){ return relo(gv); } std::vector 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 beams(const GenEvent *ge) { std::vector beamlist = ge->beams(); if ( beamlist.size() < 2 ) { std::cerr << "CANNOT FIND ANY BEAMS!" << std::endl; return std::pair(); } return std::make_pair(beamlist[0], beamlist[1]); } bool readEvent(std::shared_ptr io, std::shared_ptr evt){ return io->read_event(*evt) && !io->failed(); } - shared_ptr makeReader(std::istream & istr, + shared_ptr makeReader(std::string filename, + std::shared_ptr & istrp, std::string * errm) { shared_ptr ret; - + +#ifdef HAVE_LIBZ + if ( filename == "-" ) + istrp = make_shared(std::cin); + else + istrp = make_shared(filename.c_str()); + std::istream & istr = *istrp; +#else + if ( filename != "-" ) istrp = make_shared(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(istr); else if ( filetype == 4 ) ret = make_shared(istr); - else + else if ( filetype == 2 ) ret = make_shared(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(); + if ( ret ) { + if ( ret->failed() ) { + if ( errm ) *errm = "Problems reading from HepMC file. "; + ret = shared_ptr(); + } + return ret; + } + if ( !ret && filename == "-" ) { + if ( errm ) *errm += "Problems reading HepMC from stdin. No header found. "; + return shared_ptr(); } + // 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 & stripid) { // std::cerr << "Stripping event " << ge.event_number() << std::endl; vector 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 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 weightNames(const GenEvent & ge) { vector ret; try { ret = ge.weight_names(""); } catch (HepMC3::WeightError & w) { return vector(); } for ( string & w : ret ) if ( w == "Weight" || w == "0" || w == "Default" ) w = ""; return ret; } std::valarray weights(const GenEvent & ge) { return std::valarray(&ge.weights()[0], ge.weights().size()); } } }