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,94 +1,96 @@ // -*- C++ -*- #ifndef RIVET_RivetHepMC_HH #define RIVET_RivetHepMC_HH #ifdef ENABLE_HEPMC_3 #include "HepMC3/HepMC3.h" #include "HepMC3/Relatives.h" #include "HepMC3/Reader.h" namespace Rivet{ namespace RivetHepMC = HepMC3; using RivetHepMC::ConstGenParticlePtr; using RivetHepMC::ConstGenVertexPtr; using RivetHepMC::Relatives; using RivetHepMC::ConstGenHeavyIonPtr; using HepMC_IO_type = RivetHepMC::Reader; using PdfInfo = RivetHepMC::GenPdfInfo; } #else #include "HepMC/GenEvent.h" #include "HepMC/GenParticle.h" #include "HepMC/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 #define ConstGenParticlePtr const HepMC::GenParticle* #define ConstGenVertexPtr const HepMC::GenVertex* #define ConstGenHeavyIonPtr const HepMC::HeavyIon* /// @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{ 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::vector beams(const GenEvent *ge); - std::shared_ptr makeReader(std::istream &istr); - bool readEvent(std::shared_ptr io, std::shared_ptr evt); + std::shared_ptr makeReader(std::istream &istr, + std::string * errm = 0); + bool readEvent(std::shared_ptr io, + std::shared_ptr evt); } } #endif diff --git a/src/Core/Run.cc b/src/Core/Run.cc --- a/src/Core/Run.cc +++ b/src/Core/Run.cc @@ -1,173 +1,176 @@ // -*- C++ -*- #include "Rivet/Run.hh" #include "Rivet/AnalysisHandler.hh" #include "Rivet/Math/MathUtils.hh" #include "Rivet/Tools/RivetPaths.hh" /// @todo reinstate zlib once HepMC3 stream reading is ok #include "zstr/zstr.hpp" #include namespace Rivet { Run::Run(AnalysisHandler& ah) : _ah(ah), _fileweight(1.0), _xs(NAN) { } Run::~Run() { } Run& Run::setCrossSection(const double xs) { _xs = xs; return *this; } double Run::crossSection() const { return _ah.crossSection(); } Run& Run::setListAnalyses(const bool dolist) { _listAnalyses = dolist; return *this; } // Fill event and check for a bad read state bool Run::readEvent() { /// @todo Clear rather than new the GenEvent object per-event? _evt.reset(new GenEvent()); if(!HepMCUtils::readEvent(_hepmcReader, _evt)){ Log::getLog("Rivet.Run") << Log::DEBUG << "Read failed. End of file?" << endl; return false; } // Rescale event weights by file-level weight, if scaling is non-trivial if (!fuzzyEquals(_fileweight, 1.0)) { for (size_t i = 0; i < (size_t) _evt->weights().size(); ++i) { _evt->weights()[i] *= _fileweight; } } return true; } - /* - // Fill event and check for a bad read state --- to skip, maybe HEPMC3 will have a better way - bool Run::skipEvent() { - if (_io->rdstate() != 0 || !_io->fill_next_event(_evt.get()) ) { - Log::getLog("Rivet.Run") << Log::DEBUG << "Read failed. End of file?" << endl; - return false; - } - return true; - } -*/ + bool Run::openFile(const std::string& evtfile, double weight) { // Set current weight-scaling member _fileweight = weight; + // In case makeReader fails. + std::string errormessage; + // Set up HepMC input reader objects if (evtfile == "-") { - _hepmcReader = HepMCUtils::makeReader(std::cin); +#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 + 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.reset(new zstr::ifstream(evtfile.c_str())); - #else - _istr.reset(new std::fstream(evtfile.c_str(), std::ios::in)); - #endif - _hepmcReader = HepMCUtils::makeReader(*_istr); - + _istr = make_shared(evtfile.c_str()); +#else + _istr = make_shared(evtfile.c_str()); +#endif + _hepmcReader = HepMCUtils::makeReader(*_istr, &errormessage); + } + if (_hepmcReader == nullptr) { - Log::getLog("Rivet.Run") << Log::ERROR << "Read error on file " << evtfile << endl; + Log::getLog("Rivet.Run") + << Log::ERROR << "Read error on file '" << evtfile << "' " + << errormessage << endl; return false; } return true; } bool Run::init(const std::string& evtfile, double weight) { if (!openFile(evtfile, weight)) return false; // Read first event to define run conditions bool ok = readEvent(); if (!ok) return false; if(HepMCUtils::particles(_evt).size() == 0){ Log::getLog("Rivet.Run") << Log::ERROR << "Empty first event." << endl; return false; } // Initialise AnalysisHandler with beam information from first event _ah.init(*_evt); // Set cross-section from command line if (!std::isnan(_xs)) { Log::getLog("Rivet.Run") << Log::DEBUG << "Setting user cross-section = " << _xs << " pb" << endl; _ah.setCrossSection(_xs); } // List the chosen & compatible analyses if requested if (_listAnalyses) { for (const std::string& ana : _ah.analysisNames()) { cout << ana << endl; } } return true; } bool Run::processEvent() { // Set cross-section if found in event and not from command line #if defined ENABLE_HEPMC_3 if (std::isnan(_xs) && _evt->cross_section()) { const double xs = _evt->cross_section()->xsec(); ///< in pb Log::getLog("Rivet.Run") << Log::DEBUG << "Setting cross-section = " << xs << " pb" << endl; _ah.setCrossSection(xs); } #elif defined HEPMC_HAS_CROSS_SECTION if (std::isnan(_xs) && _evt->cross_section()) { const double xs = _evt->cross_section()->cross_section(); ///< in pb Log::getLog("Rivet.Run") << Log::DEBUG << "Setting cross-section = " << xs << " pb" << endl; _ah.setCrossSection(xs); } #endif // Complain about absence of cross-section if required! if (_ah.needCrossSection() && !_ah.hasCrossSection()) { Log::getLog("Rivet.Run") << Log::ERROR << "Total cross-section needed for at least one of the analyses. " << "Please set it (on the command line with '-x' if using the 'rivet' program)" << endl; return false; } // Analyze event _ah.analyze(*_evt); return true; } bool Run::finalize() { _evt.reset(); /// @todo reinstate for HepMC3 //_istr.reset(); if (!std::isnan(_xs)) _ah.setCrossSection(_xs); _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,116 +1,117 @@ // -*- C++ -*- #include "Rivet/Tools/RivetHepMC.hh" 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{ 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::vector beams(const GenEvent *ge){ pair beams = ge->beam_particles(); return std::vector{beams.first, beams.second}; } - std::shared_ptr makeReader(std::istream &istr){ + std::shared_ptr makeReader(std::istream &istr, + std::string *) { 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; } } } 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,91 +1,106 @@ // -*- C++ -*- #include "Rivet/Tools/RivetHepMC.hh" #include "HepMC3/ReaderAscii.h" #include "HepMC3/ReaderAsciiHepMC2.h" namespace Rivet{ namespace HepMCUtils{ 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::vector beams(const GenEvent *ge){ return ge->beams(); } bool readEvent(std::shared_ptr io, std::shared_ptr evt){ io->read_event(*evt); return !io->failed(); } - shared_ptr makeReader(std::istream &istr){ - if(&istr == &std::cin) return make_shared(istr); + shared_ptr makeReader(std::istream & istr, + std::string * errm) { + shared_ptr ret; - istr.seekg(istr.beg); - std::string line1, line2; - - while(line1.empty()){ - std::getline(istr, line1); + // First scan forward and check if there is some hint as to what + // kind of file we are looking att. + int nchar = 256; + std::string header; + while ( nchar-- && !istr.eof() ) + header += char(istr.get()); + // If this stream was too short to contain any reasonable number + // of events, just give up. + if ( !istr ) { + if ( errm ) *errm = "Could not find HepMC header information"; + return shared_ptr(); } - - while(line2.empty()){ - std::getline(istr, line2); + + // Now reset the stream to its original state ... + for ( int i = header.length() - 1; i >= 0; --i ) + istr.putback(header[i]); + + // ... and check which kind of format it was and create the + // corresponding reader. First try the HepM3 ascii format. + if ( header.find("HepMC::Asciiv3-START_EVENT_LISTING") != + std::string::npos ) + ret = make_shared(istr); + + // Check if the file is written by WriterRoot or WriterRootTree. + else if ( header.substr(0, 4) == "root" ) { + if ( errm ) *errm = "Rivet cancurrently not read HepMC root files."; + return ret; } - - istr.seekg(istr.beg); - - // if this is absent it doesn't appear to be a HepMC file :( - if(line1.find("HepMC::Version") == std::string::npos) return nullptr; - - shared_ptr result; - - // Looks like the new HepMC 3 format! - if(line2.find("HepMC::Asciiv3") != std::string::npos){ - result = make_shared(istr); - }else{ - // assume old HepMC 2 format from here - result = make_shared(istr); + + // The default is to assume it is a good old HepMC2 ascii file. + else { + ret = make_shared(istr); } - - if(result->failed()) result.reset((RivetHepMC::Reader*)nullptr); - return result; + + // Check that everything was ok. + if ( ret->failed() ) { + if ( errm ) *errm = "Problems reading from HepMC file."; + ret = shared_ptr(); + } + + return ret; } + } }