diff --git a/.hgignore b/.hgignore --- a/.hgignore +++ b/.hgignore @@ -1,80 +1,92 @@ ~$ ^(anas|test|run).*$ \.orig$ \.(o|Po|lo|Plo|la|a|so|dylib|pyc|tar\.bz2|tar\.gz|fifo|hepmc)$ ^config\.guess$ ^config\.status$ ^config\.sub$ ^config\.log$ ^configure$ ^depcomp$ ^compile$ ^ar-lib$ ^install-sh$ ^INSTALL$ ^libtool$ ^test-driver$ ^ltmain\.sh$ ^m4/ltoptions\.m4$ ^m4/ltsugar\.m4$ ^m4/ltversion\.m4$ ^m4/lt~obsolete\.m4$ ^missing$ ^autom4te\.cache$ ^include/Rivet/Config/stamp-h.$ ^tmp$ ^rivet.pc$ Makefile$ Makefile\.in$ \.(pdf|toc|bbl|blg|sty|bib|html|tex)$ Doxyfile$ \.libs$ \.deps$ #Rivet.yoda #.*/Rivet\.yoda ^bin/rivet-buildplugin$ ^bin/rivet-config$ ^bin/rivet-nopy$ ^aclocal\.m4$ ^pyext/build$ ^pyext/rivet/core.cpp$ ^pyext/rivet/rivetwrap\.py$ ^pyext/rivet/rivetwrap_wrap\.cc$ ^pyext/rivet/fix-out-of-source$ ^pyext/setup\.py$ ^rivetenv\.(sh|csh)$ ^test/test(Api|Cmp|NaN)$ ^include/Rivet/Config/DummyConfig\.hh\.in$ ^include/Rivet/Config/BuildOptions\.hh$ ^include/Rivet/Config/DummyConfig\.hh$ ^include/Rivet/Config/RivetConfig\.hh$ ^doc/.*\.(log|aux|out)$ ^doc/ana ^test/log$ ^test/.*\.log$ ^test/out\.yoda$ ^test/NaN.aida$ ^test/Rivet.yoda$ ^test/testBoost$ ^test/.*\.trs$ ^test/testMatVec$ ^test/testMath$ ^test/testBeams$ ^dev$ ^devval$ ^newanalyses$ ^.*plots$ ^a\.out$ ^a\.out\.dSYM$ ^Rivet-.\..\..$ ^local/.*$ ^(for|analyses)\d\d\d$ ^analyses/data$ +# added by aHg on Mon Feb 25 14:51:05 2019 +syntax: glob +bin/LHC-13-Jets-0.hepmc.gz + +# added by aHg on Mon Feb 25 16:16:10 2019 +syntax: glob +bin/*.gz + +# added by aHg on Thu Feb 28 11:44:29 2019 +syntax: glob +bin/hepmc-compress + # added by aHg on Tue Mar 5 14:30:12 2019 syntax: glob doc/rivet-manual.dvi # added by aHg on Tue Mar 5 14:30:21 2019 syntax: glob doc/rivet-manual.ps diff --git a/bin/Makefile.am b/bin/Makefile.am --- a/bin/Makefile.am +++ b/bin/Makefile.am @@ -1,40 +1,46 @@ bin_SCRIPTS = rivet-config rivet-buildplugin dist_bin_SCRIPTS = make-plots EXTRA_DIST = RIVETPROGS = \ rivet \ rivet-mkanalysis \ rivet-findid rivet-which \ rivet-cmphistos rivet-mkhtml if ENABLE_PYEXT dist_bin_SCRIPTS += $(RIVETPROGS) else EXTRA_DIST += $(RIVETPROGS) endif ## bash completion if ENABLE_PYEXT dist_pkgdata_DATA = rivet-completion bashcomp_dir = $(DESTDIR)$(prefix)/etc/bash_completion.d install-data-local: if [[ -d "$(bashcomp_dir)" && -w "$(bashcomp_dir)" ]]; then \ $(install_sh_DATA) rivet-completion $(bashcomp_dir)/; fi uninstall-local: rm -f $(bashcomp_dir)/rivet-completion else EXTRA_DIST += rivet-completion endif ## No-Python Rivet program noinst_PROGRAMS = rivet-nopy rivet_nopy_SOURCES = rivet-nopy.cc rivet_nopy_CPPFLAGS = -I$(top_srcdir)/include $(AM_CPPFLAGS) $(CPPFLAGS) rivet_nopy_LDADD = $(top_builddir)/src/libRivet.la if ENABLE_HEPMC_3 +noinst_PROGRAMS += hepmc-compress +hepmc_compress_SOURCES = hepmc-compress.cc +hepmc_compress_CPPFLAGS = -I$(top_srcdir)/include $(AM_CPPFLAGS) $(CPPFLAGS) +hepmc_compress_LDADD = $(top_builddir)/src/libRivet.la -lHepMC3 -lHepMC3search +hepmc_compress_LDFLAGS = -L$(HEPMC3LIBPATH) + rivet_nopy_LDFLAGS = -L$(HEPMC3LIBPATH) rivet_nopy_LDADD += -lHepMC3 -lHepMC3search else rivet_nopy_LDFLAGS = -L$(HEPMCLIBPATH) rivet_nopy_LDADD += -lHepMC endif diff --git a/bin/hepmc-compress.cc b/bin/hepmc-compress.cc new file mode 100644 --- /dev/null +++ b/bin/hepmc-compress.cc @@ -0,0 +1,54 @@ +#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) { + + if (argc < 3) { + cerr << "Usage: " << argv[0] + << " " << endl; + return 1; + } + + int outputmode = 0; + if ( argc >= 4 ) outputmode = atoi(argv[3]); + Rivet::zstr::ifstream input(argv[1]); + Rivet::zstr::ofstream output(argv[2]); + + std::shared_ptr + reader = Rivet::HepMCUtils::makeReader(input); + + std::shared_ptr + evt = make_shared(); + + shared_ptr writer; + if ( outputmode == 0 ) + writer = make_shared(output); + else { + auto compressed = make_shared(output); + if ( outputmode >= 2 ) compressed->use_integers(); + if ( outputmode == 3 ) { + 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; + } + + while(reader && Rivet::HepMCUtils::readEvent(reader, evt) ) { + cout << evt->event_number() << endl; + writer->write_event(*evt); + } + + return 0; +} + + diff --git a/include/Rivet/Makefile.am b/include/Rivet/Makefile.am --- a/include/Rivet/Makefile.am +++ b/include/Rivet/Makefile.am @@ -1,161 +1,163 @@ ## Internal headers - not to be installed nobase_dist_noinst_HEADERS = ## Public headers - to be installed nobase_pkginclude_HEADERS = ## Rivet interface nobase_pkginclude_HEADERS += \ Rivet.hh \ Run.hh \ Event.hh \ ParticleBase.hh \ Particle.fhh Particle.hh \ Jet.fhh Jet.hh \ Projection.fhh Projection.hh \ ProjectionApplier.hh \ ProjectionHandler.hh \ Analysis.hh \ AnalysisHandler.hh \ AnalysisInfo.hh \ AnalysisBuilder.hh \ AnalysisLoader.hh ## Build config stuff nobase_pkginclude_HEADERS += \ Config/RivetCommon.hh \ Config/RivetConfig.hh \ Config/BuildOptions.hh ## Projections nobase_pkginclude_HEADERS += \ Projections/AxesDefinition.hh \ Projections/Beam.hh \ Projections/BeamThrust.hh \ Projections/CentralEtHCM.hh \ Projections/ChargedFinalState.hh \ Projections/ChargedLeptons.hh \ Projections/ConstLossyFinalState.hh \ Projections/DISFinalState.hh \ Projections/DISKinematics.hh \ Projections/DISLepton.hh \ Projections/DressedLeptons.hh \ Projections/FastJets.hh \ Projections/FinalPartons.hh \ Projections/FinalState.hh \ Projections/FoxWolframMoments.hh \ Projections/FParameter.hh \ Projections/HadronicFinalState.hh \ Projections/HeavyHadrons.hh \ Projections/Hemispheres.hh \ Projections/IdentifiedFinalState.hh \ Projections/InitialQuarks.hh \ Projections/InvMassFinalState.hh \ Projections/JetAlg.hh \ Projections/JetShape.hh \ Projections/LeadingParticlesFinalState.hh \ Projections/LossyFinalState.hh \ Projections/MergedFinalState.hh \ Projections/MissingMomentum.hh \ Projections/NeutralFinalState.hh \ Projections/NonHadronicFinalState.hh \ Projections/NonPromptFinalState.hh \ Projections/ParisiTensor.hh \ Projections/ParticleFinder.hh \ Projections/PartonicTops.hh \ Projections/PrimaryHadrons.hh \ Projections/PromptFinalState.hh \ Projections/SmearedParticles.hh \ Projections/SmearedJets.hh \ Projections/SmearedMET.hh \ Projections/Sphericity.hh \ Projections/Spherocity.hh \ Projections/TauFinder.hh \ Projections/Thrust.hh \ Projections/TriggerCDFRun0Run1.hh \ Projections/TriggerCDFRun2.hh \ Projections/TriggerUA5.hh \ Projections/UnstableFinalState.hh \ Projections/VetoedFinalState.hh \ Projections/VisibleFinalState.hh \ Projections/WFinder.hh \ Projections/ZFinder.hh ## Meta-projection convenience headers nobase_pkginclude_HEADERS += \ Projections/FinalStates.hh \ Projections/Smearing.hh ## Analysis base class headers # TODO: Move to Rivet/AnalysisTools header dir? nobase_pkginclude_HEADERS += \ Analyses/MC_ParticleAnalysis.hh \ Analyses/MC_JetAnalysis.hh \ Analyses/MC_JetSplittings.hh ## Tools nobase_pkginclude_HEADERS += \ Tools/BeamConstraint.hh \ Tools/BinnedHistogram.hh \ Tools/CentralityBinner.hh \ Tools/Cmp.fhh \ Tools/Cmp.hh \ Tools/Cutflow.hh \ Tools/Cuts.fhh \ Tools/Cuts.hh \ Tools/Exceptions.hh \ Tools/JetUtils.hh \ Tools/Logging.hh \ Tools/Random.hh \ Tools/ParticleBaseUtils.hh \ Tools/ParticleIdUtils.hh \ Tools/ParticleUtils.hh \ Tools/ParticleName.hh \ Tools/PrettyPrint.hh \ + Tools/ReaderCompressedAscii.hh \ Tools/RivetPaths.hh \ Tools/RivetSTL.hh \ Tools/RivetFastJet.hh \ Tools/RivetHepMC.hh \ Tools/RivetYODA.hh \ Tools/RivetMT2.hh \ Tools/SmearingFunctions.hh \ Tools/MomentumSmearingFunctions.hh \ Tools/ParticleSmearingFunctions.hh \ Tools/JetSmearingFunctions.hh \ Tools/TypeTraits.hh \ - Tools/Utils.hh + Tools/Utils.hh \ + Tools/WriterCompressedAscii.hh nobase_dist_noinst_HEADERS += \ Tools/osdir.hh ## Maths nobase_pkginclude_HEADERS += \ Math/Matrices.hh \ Math/Vector3.hh \ Math/VectorN.hh \ Math/MatrixN.hh \ Math/MatrixDiag.hh \ Math/MathHeader.hh \ Math/Vectors.hh \ Math/LorentzTrans.hh \ Math/Matrix3.hh \ Math/MathUtils.hh \ Math/Vector4.hh \ Math/Math.hh \ Math/Units.hh \ Math/Constants.hh \ Math/eigen/util.h \ Math/eigen/regressioninternal.h \ Math/eigen/regression.h \ Math/eigen/vector.h \ Math/eigen/ludecompositionbase.h \ Math/eigen/ludecomposition.h \ Math/eigen/linearsolver.h \ Math/eigen/linearsolverbase.h \ Math/eigen/matrix.h \ Math/eigen/vectorbase.h \ Math/eigen/projective.h \ Math/eigen/matrixbase.h diff --git a/include/Rivet/Tools/ReaderCompressedAscii.hh b/include/Rivet/Tools/ReaderCompressedAscii.hh new file mode 100644 --- /dev/null +++ b/include/Rivet/Tools/ReaderCompressedAscii.hh @@ -0,0 +1,165 @@ +// -*- C++ -*- +// +// This file is part of HepMC +// Copyright (C) 2014-2019 The HepMC collaboration (see AUTHORS for details) +// +#ifndef HEPMC3_READERCOMPRESSEDASCII_H +#define HEPMC3_READERCOMPRESSEDASCII_H +/// +/// @file ReaderCompressedAscii.h +/// @brief Definition of class \b ReaderCompressedAscii +/// +/// @class HepMC3::ReaderCompressedAscii +/// @brief GenEvent I/O parsing for structured text files +/// +/// @ingroup IO +/// +#include "Rivet/Tools/RivetHepMC.hh" +#include "HepMC3/Reader.h" +#include "HepMC3/GenEvent.h" +#include +#include +#include + +namespace Rivet { + + +class ReaderCompressedAscii : public HepMC3::Reader { + +public: + + typedef HepMC3::GenParticlePtr GenParticlePtr; + typedef HepMC3::GenVertexPtr GenVertexPtr; + +public: + + /// @brief Constructor + /// @warning If file already exists, it will be cleared before writing + ReaderCompressedAscii(const std::string& filename); + + /// The ctor to read from stdin + ReaderCompressedAscii(std::istream &); + + /// @brief Destructor + ~ReaderCompressedAscii(); + + /// @brief Load event from file + /// + /// @param[out] evt Event to be filled + bool read_event(GenEvent& evt); + + /// @brief Return status of the stream + bool failed() { + return !(*m_stream); + } + + /// @brief Close file stream + void close(); + +private: + + /// @brief Unsecape '\' and '\n' characters in string + std::string unescape(const std::string& s); + + /// @name Read helpers + //@{ + + /// @brief Parse event + /// + /// Helper routine for parsing event information + /// @return vertices count and particles count for verification + std::pair parse_event_information(); + + /// @brief Parse weight value lines + /// + /// Helper routine for parsing weight value information + bool parse_weight_values(); + + /// @brief Parse precision + /// + /// Helper routine for parsing precision information + bool parse_precision(); + + /// @brief Parse vertex + /// + /// Helper routine for parsing single event information + bool parse_vertex_information(); + + /// @brief Parse particle + /// + /// Helper routine for parsing single particle information + bool parse_particle_information(); + + /// @brief Parse attribute + /// + /// Helper routine for parsing single attribute information + bool parse_attribute(); + + /// @brief Parse run-level attribute. + /// + /// Helper routine for parsing single attribute information + bool parse_run_attribute(); + + /// @brief Parse run-level weight names. + /// + /// Helper routine for parsing a line with information about + /// weight names. + bool parse_weight_names(); + + /// @brief Parse run-level tool information. + /// + /// Helper routine for parsing a line with information about + /// tools being used. + bool parse_tool(); + + /// @brief Read position information + /// + /// Reads position information from the current line and sets the + /// information in the given vertex. If no vertex is given the root + /// vertiex is assumed. + bool read_position(GenVertexPtr v = GenVertexPtr()); + + /// @brief Read momentum information + /// + /// Reads momentum and mass information from the current line and + /// sets the information in the given particle. + bool read_momentum(GenParticlePtr p); + + //@} + + +private: + + std::ifstream m_file; //!< Input file + std::istream* m_stream; //!< The stream being read from + + std::istringstream is; //!< A stream to read from the current line. + + GenEvent * m_evt; //!< The event being read in. + + double m_precision_phi; //!< Input precision in phi + double m_precision_eta; //!< Input precision in eta + double m_precision_e; //!< Input precision energy + double m_precision_m; //!< Input precision mass + bool m_using_integers; //!< Reading integers + + map m_masses; //!< Keep track of masses being read. + + /// Keep track of read particles + vector m_particles; + /// Keep track of read particles + vector m_ppvx; + /// Keep track of read vertices + map m_vertices; + /// Keep track of incoming particles to vertices + map > m_vpin; + + /** @brief Store attributes global to the run being written/read. */ + std::map< std::string, shared_ptr > m_global_attributes; + +}; + + +} // namespace HepMC3 + +#endif diff --git a/include/Rivet/Tools/WriterCompressedAscii.hh b/include/Rivet/Tools/WriterCompressedAscii.hh new file mode 100644 --- /dev/null +++ b/include/Rivet/Tools/WriterCompressedAscii.hh @@ -0,0 +1,220 @@ +// -*- C++ -*- +// +// This file is part of HepMC +// Copyright (C) 2014-2019 The HepMC collaboration (see AUTHORS for details) +// +#ifndef HEPMC3_WRITERCOMPRESSEDASCII_H +#define HEPMC3_WRITERCOMPRESSEDASCII_H +/// +/// @file WriterCompressedAscii.h +/// @brief Definition of class \b WriterCompressedAscii +/// +/// @class HepMC3::WriterCompressedAscii +/// @brief GenEvent I/O serialization for structured text files +/// +/// @ingroup IO +/// +#include "Rivet/Tools/RivetHepMC.hh" +#include "HepMC3/Writer.h" +#include "HepMC3/GenEvent.h" +#include "HepMC3/GenRunInfo.h" +#include +#include + +namespace Rivet { + +class WriterCompressedAscii : public HepMC3::Writer { + +public: + + typedef HepMC3::GenRunInfo GenRunInfo; + typedef HepMC3::FourVector FourVector; + typedef HepMC3::ConstGenVertexPtr ConstGenVertexPtr; + typedef HepMC3::ConstGenParticlePtr ConstGenParticlePtr; + + /// @brief Constructor + /// @warning If file already exists, it will be cleared before writing + WriterCompressedAscii(const std::string& filename, + shared_ptr run = shared_ptr()); + + /// @brief Constructor from ostream + WriterCompressedAscii(std::ostream& stream, + shared_ptr run = shared_ptr()); + + /// @brief Destructor + ~WriterCompressedAscii(); + + /// @brief Write event to file + /// + /// @param[in] evt Event to be serialized + void write_event(const GenEvent& evt); + + /// @brief Write the GenRunInfo object to file. + void write_run_info(); + + /// @brief Return status of the stream + bool failed() { return (bool)m_file.rdstate(); } + + /// @brief Close file stream + void close(); + + /// @brief Use cartesian coordinates + /// + /// Momenta and positions will be written out as doubles using + /// standard cartesian coordinates. + void use_doubles() { + m_use_integers = false; + } + + /// @brief Use cylindical coordinates + /// + /// Momenta and positions will be written out as integers using + /// eta-phi coordinates. + void use_integers() { + m_use_integers = true; + } + + /// @brief Add a particle id to be stripped + /// + /// Specify the PDG id of a (unobservable) particle which will be + /// attempted to be removed from the event befor writing. + void add_stripid(long pdgid) { + m_stripid.insert(pdgid); + } + + /// @brief Remove onobservable particles + /// + /// Go through the event and try to take away all intermediate + /// particles specified in add_stripid() + void strip(GenEvent & e); + + /// @brief Set output double precision + /// + /// General output precision for double + void set_precision(int prec) { + m_precision = prec; + } + + /// @brief Set output precision in phi + /// + /// Azimuth angles will be written out as integers corresponding to + /// this precision + void set_precision_phi(double prec) { + m_precision_phi = prec; + } + + /// @brief Set output precision in eta + /// + /// Pseudo repisities will be written out as integers corresponding to + /// this precision + void set_precision_eta(double prec) { + m_precision_eta = prec; + } + + /// @brief Set output precision in energy + /// + /// Energies will be written out as integers corresponding to + /// this precision (in GeV) + void set_precision_e(double prec) { + m_precision_e = prec; + } + + /// @brief Set output precision in mass + /// + /// Masses will be written out as integers corresponding to + /// this precision (in GeV) + void set_precision_m(double prec) { + m_precision_m = prec; + } + + /// @brief Return output precision for doubles. + int precision() const { + return m_precision; + } + + /// @brief Return output precision for azimuth angles. + double precision_phi() const { + return m_precision_phi; + } + + /// @brief Return output precision for pseudo rapisities. + double precision_eta() const { + return m_precision_eta; + } + + /// @brief Return output precision for energies. + double precision_e() const { + return m_precision_e; + } + + /// @brief Return output precision for masses. + double precision_m() const { + return m_precision_m; + } + + /// @brief Internal function to calculate the pseudo rapidity. + double psrap(const FourVector & p) const; + +private: + + /// @brief Escape '\' and '\n' characters in string + std::string escape(const std::string& s) const; + + //@} + + + /// @name Write helpers + //@{ + + /// @brief Inline function for writing positions + void write_position(FourVector pos); + + /// @brief Inline function for writing momenta + void write_momentum(FourVector p); + + /// @brief Inline function for writing mass of a particle + void write_mass(ConstGenParticlePtr p); + + /// @brief Write vertex + /// + /// Helper routine for writing single vertex to file + void write_vertex(ConstGenVertexPtr v); + + /// @brief Write particle + /// + /// Helper routine for writing single particle to file + void write_particle(ConstGenParticlePtr p); + + /// @brief Helper function to access the root vertex + ConstGenVertexPtr rootvertex() { + vector beams = m_current->beams(); + if ( beams.empty() ) return ConstGenVertexPtr(); + return beams[0]->production_vertex(); + } + //@} + +private: + + bool m_use_integers; //!< Compress by using intergers and + //! cylindrical coordinates + + std::ofstream m_file; //!< Output file + std::ostream* m_stream; //!< Output stream + + double m_precision_phi; //!< Output precision in phi + double m_precision_eta; //!< Output precision in eta + double m_precision_e; //!< Output precision energy + double m_precision_m; //!< Output precision mass + int m_precision; //!< General double output precision + set m_stripid; //!< Strip matching intermediate particles + map m_masses; //!< Keep track of masses being written. + + const GenEvent * m_current; //!< The event being written + std::ostringstream os; //!< Internal stream to write an entire event. + +}; + + +} // namespace HepMC3 + +#endif diff --git a/src/Tools/Makefile.am b/src/Tools/Makefile.am --- a/src/Tools/Makefile.am +++ b/src/Tools/Makefile.am @@ -1,26 +1,26 @@ noinst_LTLIBRARIES = libRivetTools.la libRivetTools_la_SOURCES = \ BinnedHistogram.cc \ Cuts.cc \ JetUtils.cc \ Random.cc \ Logging.cc \ ParticleUtils.cc \ ParticleName.cc \ RivetYODA.cc \ RivetMT2.cc \ RivetPaths.cc \ Utils.cc \ binreloc.c # mt2_bisect.cc if ENABLE_HEPMC_3 -libRivetTools_la_SOURCES += RivetHepMC_3.cc +libRivetTools_la_SOURCES += RivetHepMC_3.cc ReaderCompressedAscii.cc WriterCompressedAscii.cc else libRivetTools_la_SOURCES += RivetHepMC_2.cc endif dist_noinst_HEADERS = binreloc.h lester_mt2_bisect.hh mt2_bisect.hh libRivetTools_la_CPPFLAGS = $(AM_CPPFLAGS) -DENABLE_BINRELOC -DDEFAULTDATADIR=\"$(datadir)\" -DDEFAULTLIBDIR=\"$(libdir)\" diff --git a/src/Tools/ReaderCompressedAscii.cc b/src/Tools/ReaderCompressedAscii.cc new file mode 100644 --- /dev/null +++ b/src/Tools/ReaderCompressedAscii.cc @@ -0,0 +1,464 @@ +// -*- C++ -*- +// +// This file is part of HepMC +// Copyright (C) 2014-2019 The HepMC collaboration (see AUTHORS for details) +// +/// +/// @file ReaderCompressedAscii.cc +/// @brief Implementation of \b class ReaderCompressedAscii +/// +#include "Rivet/Tools/ReaderCompressedAscii.hh" + +#include "HepMC3/GenEvent.h" +#include "HepMC3/GenParticle.h" +#include "HepMC3/GenVertex.h" +#include "HepMC3/Units.h" +#include +#include + +namespace Rivet { + +using namespace HepMC3; + + +ReaderCompressedAscii::ReaderCompressedAscii(const string &filename) + : m_file(filename), m_stream(0), m_evt(0), m_precision_phi(0.001), + m_precision_eta(0.001), m_precision_e(0.001), m_precision_m(0.000001), + m_using_integers(false) { + if( !m_file.is_open() ) { + ERROR( "ReaderCompressedAscii: could not open input file: "<()); +} + + +// Ctor for reading from stdin +ReaderCompressedAscii::ReaderCompressedAscii(std::istream & stream) + : m_stream(&stream), m_evt(0), m_precision_phi(0.001), + m_precision_eta(0.001), m_precision_e(0.001), m_precision_m(0.000001), + m_using_integers(false) { + if( !m_stream ) { + ERROR( "ReaderCompressedAscii: could not open input stream " ) + } + set_run_info(make_shared()); +} + + + +ReaderCompressedAscii::~ReaderCompressedAscii() { } + + +bool ReaderCompressedAscii::read_event(GenEvent &evt) { + + if ( failed() ) return false; + + bool parsed_event_header = false; + bool is_parsing_successful = true; + pair vertices_and_particles(0,0); + std::string line; + + m_evt = &evt; + + evt.clear(); + evt.set_run_info(run_info()); + m_masses.clear(); + m_vertices.clear(); + m_particles.clear(); + m_ppvx.clear(); + m_vpin.clear(); + // + // Parse event, vertex and particle information + // + while(!failed()) { + + std::getline(*m_stream, line); + is = std::istringstream(line); + is.get(); // Remove the first character from the stream. + if ( line.empty() ) continue; + + if ( line.substr(0, 5) == "HepMC" ) { + if ( line.substr(0, 14) != "HepMC::Version" && + line.substr(0, 24) != "HepMC::CompressedAsciiv3" ) { + WARNING( "ReaderCompressedAscii: found unsupported expression " + "in header. Will close the input." ) + std::cout << line << std::endl; + } + continue; + } + + switch( line[0] ) { + case 'E': + vertices_and_particles = parse_event_information(); + if (vertices_and_particles.second < 0) { + is_parsing_successful = false; + } else { + is_parsing_successful = true; + parsed_event_header = true; + } + break; + case 'V': + is_parsing_successful = parse_vertex_information(); + break; + case 'P': + is_parsing_successful = parse_particle_information(); + break; + case 'W': + if ( parsed_event_header ) + is_parsing_successful = parse_weight_values(); + else + is_parsing_successful = parse_weight_names(); + break; + case 'C': + is_parsing_successful = parse_precision(); + break; + case 'T': + is_parsing_successful = parse_tool(); + break; + case 'A': + if ( parsed_event_header ) + is_parsing_successful = parse_attribute(); + else + is_parsing_successful = parse_run_attribute(); + break; + default: + WARNING( "ReaderCompressedAscii: skipping unrecognised prefix: " + << line[0] ) + is_parsing_successful = true; + break; + } + + if( !is_parsing_successful ) break; + + // Check for next event + if ( parsed_event_header && + ( m_stream->peek() == 'E' || m_stream->peek() == 'H') )break; + } + + // Set the production vertex for all particles. + for ( int ip = 0, Np = m_particles.size(); ip < Np; ++ip ) + if ( m_ppvx[ip] && m_vertices[m_ppvx[ip]] ) + m_vertices[m_ppvx[ip]]->add_particle_out(m_particles[ip]); + + // Add the incoming particles to all vertices + for ( auto iv : m_vertices ) + for ( auto ip : m_vpin[iv.first] ) iv.second->add_particle_in(m_particles[ip - 1]); + + // When all particles and vertices are connected we add all of them to the event. + for ( auto p : m_particles ) evt.add_particle(p); + for ( auto v : m_vertices ) evt.add_vertex(v.second); + + // Check if all particles and vertices were parsed + if ((int)m_evt->particles().size() > vertices_and_particles.second ) { + ERROR( "ReaderCompressedAscii: too many particles were parsed" ) + printf("%zu vs %i expected\n",m_evt->particles().size(),vertices_and_particles.second ); + is_parsing_successful = false; + } + if ((int)m_evt->particles().size() < vertices_and_particles.second ) { + ERROR( "ReaderCompressedAscii: too few particles were parsed" ) + printf("%zu vs %i expected\n",m_evt->particles().size(),vertices_and_particles.second ); + is_parsing_successful = false; + } + + if ((int)m_evt->vertices().size() > vertices_and_particles.first) { + ERROR( "ReaderCompressedAscii: too many vertices were parsed" ) + printf("%zu vs %i expected\n",m_evt->vertices().size(),vertices_and_particles.first ); + is_parsing_successful = false; + } + + if ((int)m_evt->vertices().size() < vertices_and_particles.first) { + ERROR( "ReaderCompressedAscii: too few vertices were parsed" ) + printf("%zu vs %i expected\n",m_evt->vertices().size(),vertices_and_particles.first ); + is_parsing_successful = false; + } + // Check if there were errors during parsing + if( !is_parsing_successful ) { + ERROR( "ReaderCompressedAscii: event parsing failed. Returning empty event" ) + DEBUG( 1, "Parsing failed at line:" << endl << line ) + + m_evt->clear(); + + return false; + } + return true; +} + + +pair ReaderCompressedAscii::parse_event_information() { + static const pair err(-1,-1); + pair ret(-1,-1); + int event_no = 0; + + // event number + if ( !(is >> event_no) ) return err; + m_evt->set_event_number(event_no); + + // num_vertices + if ( !(is>> ret.first) ) return err; + + // num_particles + if ( !(is>> ret.second) ) return err; + + if ( !read_position() ) return err; + + DEBUG( 10, "ReaderCompressedAscii: E: "< wts; + double w; + while ( is >> w ) wts.push_back(w); + if ( run_info() && run_info()->weight_names().size() && + run_info()->weight_names().size() != wts.size() ) + throw std::logic_error( + "ReaderCompressedAscii::parse_weight_values: " + "The number of weights ("+ + std::to_string((long long int)(wts.size()))+ + ") does not match the number weight names(" + + std::to_string((long long int)(run_info()->weight_names().size()))+ + ") in the GenRunInfo object"); + + m_evt->weights() = wts; + + return true; +} + + +bool ReaderCompressedAscii::parse_precision() { + if ( !(is >> m_precision_phi >> m_precision_eta + >> m_precision_e >> m_precision_m) ) return false; + m_using_integers = ( m_precision_phi > 0.0 ); + return true; +} + +bool ReaderCompressedAscii::parse_vertex_information() { + GenVertexPtr data = make_shared(); + + int id = 0; + if ( !(is >> id) ) return false; + + int status = 0; + if ( !(is >> status) ) return false; + data->set_status( status ); + + std::string incoming; + if ( !(is >> incoming) ) return false; + std::string::size_type i = std::string::npos; + while ( ( i = incoming.find_first_of("[,]") ) != std::string::npos ) + incoming[i] = ' '; + std::istringstream isin(incoming); + int pin = 0; + vector vpin; + while ( isin >> pin ) vpin.push_back(pin); + + if ( !read_position(data) ) return false; + + m_vertices[-id] = data; + m_vpin[-id] = vpin; + + return true; +} + + +bool ReaderCompressedAscii::parse_particle_information() { + GenParticlePtr data = make_shared(); + + int id = 0; + if ( !(is >> id) ) return false; + + int ivp = 0; + if ( !(is >> ivp) ) return false; + + int pdgid = 0; + if ( !(is >> pdgid) ) return false; + data->set_pid(pdgid); + + if ( !read_momentum(data) ) return false; + + int status = 0; + if ( !(is >> status) ) return false; + data->set_status(status); + + m_particles.push_back(data); + m_ppvx.push_back(-ivp); + + return true; +} + + +bool ReaderCompressedAscii::parse_attribute() { + int id = 0; + if ( !(is >> id ) ) return false; + + string name; + if ( !(is >> name) ) return false; + is.get(); + string contents; + if ( !std::getline(is, contents) ) return false; + shared_ptr att = + make_shared(StringAttribute(unescape(contents))); + + m_evt->add_attribute(name, att, id); + + return true; +} + +bool ReaderCompressedAscii::parse_run_attribute() { + + string name; + if ( !(is >> name) ) return false; + is.get(); + string contents; + if ( !std::getline(is, contents) ) return false; + shared_ptr att = + make_shared(StringAttribute(unescape(contents))); + + run_info()->add_attribute(name, att); + + return true; + +} + + +bool ReaderCompressedAscii::parse_weight_names() { + + vector names; + string name; + while ( is >> name ) names.push_back(name); + + run_info()->set_weight_names(names); + + return true; + +} + +bool ReaderCompressedAscii::parse_tool() { + + std::string line; + if ( !(is >> line) ) return false; + line = unescape(line); + + GenRunInfo::ToolInfo tool; + + std::string::size_type pos = line.find("\n"); + tool.name = line.substr(0, pos); + + line = line.substr(pos + 1); + pos = line.find("\n"); + tool.version = line.substr(0, pos); + + tool.description = line.substr(pos + 1); + run_info()->tools().push_back(tool); + + return true; + +} + +bool ReaderCompressedAscii::read_position(GenVertexPtr v) { + string at; + if ( !(is >> at) ) return true; + if ( at != "@" ) return false; + + if ( !m_using_integers ) { + double x = 0.0, y = 0.0, z = 0.0, t = 0.0; + if ( !(is >> x >> y >> z >> t) ) return false; + FourVector pos(x, y, z, t); + Units::convert(pos, Units::MM, m_evt->length_unit()); + v->set_position(pos); + return true; + } + + long ieta = 0; + long iphi = 0; + double p3mod = 0.0; + double t = 0.0; + if ( !(is >> ieta >> iphi >> p3mod >> t) ) return false; + + double eta = double(ieta)*m_precision_eta; + double phi = double(iphi)*m_precision_phi*M_PI; + double pt = p3mod/cosh(eta); + FourVector pos(pt*cos(phi), pt*sin(phi), pt*sinh(eta), t); + + Units::convert(pos, Units::MM, m_evt->length_unit()); + v->set_position(pos); + + return true; + +} + +bool ReaderCompressedAscii::read_momentum(GenParticlePtr p) { + if ( !m_using_integers ) { + double px = 0.0, py = 0.0, pz = 0.0, e = 0.0, m = 0.0; + if ( !(is >> px >> py >> pz >> e >> m) ) return false; + FourVector pp(px, py, pz, e); + + if ( m_evt->momentum_unit() != Units::GEV ) { + m *= 1000.0; + Units::convert(pp, Units::GEV, m_evt->momentum_unit()); + } + p->set_momentum(pp); + p->set_generated_mass(m); + return true; + } + + long iphi = 0; + long ieta = 0; + long ie = 0; + std::string sm; + if ( !(is >> ie >> ieta >> iphi >> sm) ) return false; + + double m = 0.0; + if ( sm == "*" ) { + m = m_masses[p->pid()]*m_precision_m; + } else { + m = (m_masses[p->pid()] = stol(sm))*m_precision_m; + } + + double e = double(ie)*m_precision_e; + double m2 = ( m >= 0.0? m*m: -m*m ); + double p3mod = sqrt(max(e*e - m2, 0.0)); + double eta = double(ieta)*m_precision_eta; + double phi = double(iphi)*m_precision_phi*M_PI; + double pt = p3mod/cosh(eta); + FourVector pp(pt*cos(phi), pt*sin(phi), pt*sinh(eta), e); + + if ( m_evt->momentum_unit() != Units::GEV ) { + m *= 1000.0; + Units::convert(pp, Units::GEV, m_evt->momentum_unit()); + } + + p->set_momentum(pp); + p->set_generated_mass(m); + + return true; +} + +string ReaderCompressedAscii::unescape(const string& s) { + string ret; + ret.reserve(s.length()); + for ( string::const_iterator it = s.begin(); it != s.end(); ++it ) { + if ( *it == '\\' ) { + ++it; + if ( *it == '|' ) + ret += '\n'; + else + ret += *it; + } else + ret += *it; + } + + return ret; +} + + +void ReaderCompressedAscii::close() { + if( !m_file.is_open()) return; + m_file.close(); +} + + +} // namespace HepMC3 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,111 +1,109 @@ // -*- C++ -*- #include "Rivet/Tools/RivetHepMC.hh" +#include "Rivet/Tools/ReaderCompressedAscii.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::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){ - io->read_event(*evt); - return !io->failed(); + return io->read_event(*evt) && !io->failed(); } shared_ptr makeReader(std::istream & istr, std::string * errm) { shared_ptr ret; // First scan forward and check if there is some hint as to what // kind of file we are looking att. - int nchar = 256; + int ntry = 10; 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(); + 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; } - // 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 ) + + if ( filetype == 3 ) 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; - } - - // The default is to assume it is a good old HepMC2 ascii file. - else { + else if ( filetype == 4 ) + ret = make_shared(istr); + else 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(); } return ret; } } } diff --git a/src/Tools/WriterCompressedAscii.cc b/src/Tools/WriterCompressedAscii.cc new file mode 100644 --- /dev/null +++ b/src/Tools/WriterCompressedAscii.cc @@ -0,0 +1,379 @@ +// -*- C++ -*- +// +// This file is part of HepMC +// Copyright (C) 2014-2019 The HepMC collaboration (see AUTHORS for details) +// +/// +/// @file WriterCompressedAscii.cc +/// @brief Implementation of \b class WriterCompressedAscii +/// +#include "Rivet/Tools/WriterCompressedAscii.hh" + +#include "HepMC3/Version.h" +#include "HepMC3/GenEvent.h" +#include "HepMC3/GenParticle.h" +#include "HepMC3/GenVertex.h" +#include "HepMC3/Units.h" +#include +#include //min max for VS2017 + +namespace Rivet { + +using namespace HepMC3; + + +WriterCompressedAscii::WriterCompressedAscii(const std::string &filename, shared_ptr run) + : m_use_integers(false), + m_file(filename), + m_stream(&m_file), + m_precision_phi(0.001), + m_precision_eta(0.001), + m_precision_e(0.001), + m_precision_m(0.000001), + m_precision(16), + m_current(0) { + set_run_info(run); + if ( !m_file.is_open() ) { + ERROR( "WriterCompressedAscii: could not open output file: "< run) + : m_use_integers(false), + m_file(), + m_stream(&stream), + m_precision_phi(0.001), + m_precision_eta(0.001), + m_precision_e(0.001), + m_precision_m(0.000001), + m_precision(16), + m_current(0) { + set_run_info(run); + (*m_stream) << "HepMC::Version " << version() << std::endl; + (*m_stream) << "HepMC::CompressedAsciiv3-START_EVENT_LISTING" << std::endl; + if ( run_info() ) write_run_info(); +} + + +WriterCompressedAscii::~WriterCompressedAscii() { + close(); +} + + +void WriterCompressedAscii::write_event(const GenEvent &evt) { + + if ( !m_stripid.empty() ) { + GenEvent e = evt; + cout << "#beams " << e.beams().size() << endl; + // for ( auto bp : evt.beams() ) { + // GenParticlePtr nb = e.particles()[bp->id() - 1]; + // cout << "Beam: " << bp->id() << " " << bp->pid() << " " + // << (bp->production_vertex()? bp->production_vertex()->id(): 999) << endl; + // cout << "NewBeam: " << nb->id() << " " << nb->pid() << " " + // << (nb->production_vertex()? nb->production_vertex()->id(): 999) << endl; + // } + + strip(e); + set saveid; + swap(m_stripid, saveid); + write_event(e); + swap(saveid, m_stripid); + return; + } + + m_current = &evt; + m_masses.clear(); + os = std::ostringstream(); + + if ( !run_info() ) { + set_run_info(evt.run_info()); + write_run_info(); + } else { + if ( evt.run_info() && run_info() != evt.run_info() ) { + WARNING( "WriterCompressedAscii::write_event: GenEvents contain " + "different GenRunInfo objects from - only the " + "first such object will be serialized." ) + } + } + + // Write event info + os << "E " << evt.event_number() + << " " << evt.vertices().size() + << " " << evt.particles().size(); + // Write event position if not zero + const FourVector &pos = evt.event_pos(); + if ( !pos.is_zero() ) write_position(pos); + os << endl; + + // Write conversions made for double -> ints + if ( m_use_integers ) + os << "C " << m_precision_phi << " " << m_precision_eta + << " " << m_precision_e << " " << m_precision_m << endl; + else + os << "C 0 0 0 0" << endl; + + // Write weight values if present + if ( evt.weights().size() ) { + os << "W"; + for (auto w: evt.weights()) + os << " " << w; + os << endl; + } + + // Write attributes + for ( auto vt1: evt.attributes() ) { + for ( auto vt2: vt1.second ) { + + string st; + bool status = vt2.second->to_string(st); + + if( !status ) { + WARNING( "WriterCompressedAscii::write_event: problem " + "serializing attribute: "< done; + + // Print particles + for(ConstGenParticlePtr p: evt.particles() ) { + + // Check to see if we need to write a vertex first + ConstGenVertexPtr v = p->production_vertex(); + + if ( !v ) cout << "WARMING particle " << p->id() + << " has no productions vertex!" << endl; + else if ( v->id() < 0 && done.insert(v).second ) + write_vertex(v); + + write_particle(p); + } + + bool orphans = false; + for(ConstGenVertexPtr v: evt.vertices() ) { + if ( done.insert(v).second ) { + orphans = true; + write_vertex(v); + } + } + if ( orphans ) cout << "WARMING found unconnected vertices" << endl; + + (*m_stream) << os.str(); + +} + +string WriterCompressedAscii::escape(const string& s) const { + string ret; + ret.reserve( s.length()*2 ); + for ( string::const_iterator it = s.begin(); it != s.end(); ++it ) { + switch ( *it ) { + case '\\': ret += "\\\\"; break; + case '\n': ret += "\\|"; break; + default: ret += *it; + } + } + return ret; +} + +void WriterCompressedAscii::write_vertex(ConstGenVertexPtr v) { + + os << "V " << v->id() << " " << v->status() << " "; + + bool first = true; + for ( auto p : v->particles_in() ) { + os << (first? '[': ',') << p->id(); + first = false; + } + os << "]"; + + const FourVector &pos = v->position(); + if ( !pos.is_zero() ) write_position(pos); + + os << endl; + +} + + +void WriterCompressedAscii::write_run_info() { + + // If no run info object set, create a dummy one. + if ( !run_info() ) set_run_info(make_shared()); + + vector names = run_info()->weight_names(); + + if ( !names.empty() ) { + string out = names[0]; + for ( int i = 1, N = names.size(); i < N; ++i ) + out += "\n" + names[i]; + os << "W " << escape(out) << endl; + } + + for ( int i = 0, N = run_info()->tools().size(); i < N; ++i ) { + string out = "T " + run_info()->tools()[i].name + "\n" + + run_info()->tools()[i].version + "\n" + + run_info()->tools()[i].description; + os << escape(out) << endl; + } + + + for ( auto att: run_info()->attributes() ) { + string st; + if ( ! att.second->to_string(st) ) { + WARNING ("WriterCompressedAscii::write_run_info: problem serializing attribute: "<< att.first ) + } + else { + os << "A " << att.first << " " << escape(st) << endl; + + } + } +} + +void WriterCompressedAscii:: +write_particle(ConstGenParticlePtr p) { + + ConstGenVertexPtr vp = p->production_vertex(); + + os << "P " << p->id() + << " " << (vp? vp->id(): 0) + << " " << p->pid(); + write_momentum(p->momentum()); + write_mass(p); + os << " " << p->status() << endl; + +} + +void WriterCompressedAscii::close() { + std::ofstream* ofs = dynamic_cast(m_stream); + if (ofs && !ofs->is_open()) return; + (*m_stream) << "HepMC::CompressedAsciiv3-END_EVENT_LISTING" << endl << endl; + if (ofs) ofs->close(); +} + +double WriterCompressedAscii::psrap(const FourVector & p) const { + static const double MAXETA = 100.0; + static const double MAXLOG = exp(-MAXETA); + double nom = p.p3mod() + abs(p.pz()); + double den = max(p.perp(), nom*MAXLOG); + return p.pz() > 0? log(nom/den): -log(nom/den); +} + +void WriterCompressedAscii::write_momentum(FourVector p) { + + Units::convert(p, m_current->momentum_unit(), Units::GEV); + + if ( m_use_integers ) { + os << " " << long(round(p.e()/precision_e())) + << " " << long(round(psrap(p)/precision_eta())) + << " " << long(round(p.phi()/(M_PI*precision_phi()))); + return; + } + + std::ostringstream osse; + osse << std::scientific << setprecision(precision()) + << " " << p.px() + << " " << p.py() + << " " << p.pz() + << " " << p.e(); + os << osse.str(); + +} + +void WriterCompressedAscii::write_mass(ConstGenParticlePtr p) { + + double m = p->generated_mass(); + if ( m_current->momentum_unit() != Units::GEV ) m /= 1000.0; + if ( m_use_integers ) { + long im = long(round(m/precision_m())); + auto pm = m_masses.find(p->pid()); + if ( pm == m_masses.end() || pm->second != im ) { + os << " " << im; + m_masses[p->pid()] = im; + } else { + os << " *"; + } + return; + } + std::ostringstream osse; + osse << std::scientific << setprecision(precision()) + << " " << m; + os << osse.str(); + +} + +void WriterCompressedAscii::write_position(FourVector pos) { + + Units::convert(pos, m_current->length_unit(), Units::MM); + std::ostringstream osse; + osse << std::scientific << setprecision(precision()); + + if ( m_use_integers ) { + osse << " @ " << long(psrap(pos)/precision_eta()) + << " " << long(pos.phi()/(M_PI*precision_phi())) + << " " << pos.p3mod() + << " " << pos.t(); + } else { + osse << " @ " << pos.x() + << " " << pos.y() + << " " << pos.z() + << " " << pos.t(); + } + os << osse.str(); +} + +void WriterCompressedAscii::strip(GenEvent & e) { + std::cout << "Stripping event " << e.event_number() << std::endl; + vector allparticles = e.particles(); + for ( auto & p : allparticles ) { + if ( !p->production_vertex() || !p->end_vertex() || + m_stripid.find(p->pid()) == m_stripid.end() || + 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; + } + vp->remove_particle_out(p); + ve->remove_particle_in(p); + if ( vp->particles_out().empty() ) { + for ( auto pi : vp->particles_in() ) { + vp->remove_particle_in(pi); + ve->add_particle_in(pi); + // std::cout << "Removing vertex " << vp->id() << std::endl; + e.remove_vertex(vp); + } + } + else if ( ve->particles_in().empty() ) { + for ( auto po : ve->particles_out() ) { + ve->remove_particle_out(po); + vp->add_particle_out(po); + // std::cout << "Removing vertex " << ve->id() << std::endl; + e.remove_vertex(ve); + } + } + e.remove_particle(p); + } +} + +}