diff --git a/.hgignore b/.hgignore
--- a/.hgignore
+++ b/.hgignore
@@ -1,72 +1,84 @@
~$
^(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
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]
+ << " <input-hepmcfile> <output-hepmcfile>" << 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<Rivet::HepMC_IO_type>
+ reader = Rivet::HepMCUtils::makeReader(input);
+
+ std::shared_ptr<Rivet::RivetHepMC::GenEvent>
+ evt = make_shared<Rivet::RivetHepMC::GenEvent>();
+
+ shared_ptr<HepMC3::Writer> writer;
+ if ( outputmode == 0 )
+ writer = make_shared<HepMC3::WriterAscii>(output);
+ else {
+ auto compressed = make_shared<Rivet::WriterCompressedAscii>(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,163 @@
+// -*- 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 <string>
+#include <fstream>
+#include <istream>
+
+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<int,int> 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<long,long> m_masses; //!< Keep track of masses being read.
+
+ /// Keep track of read particles
+ vector<GenParticlePtr> m_particles;
+ /// Keep track of read particles
+ vector< pair<int,int> > m_pvx;
+ /// Keep track of read vertices
+ map<int,GenVertexPtr> m_vertices;
+
+ /** @brief Store attributes global to the run being written/read. */
+ std::map< std::string, shared_ptr<HepMC3::Attribute> > 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 <string>
+#include <fstream>
+
+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<GenRunInfo> run = shared_ptr<GenRunInfo>());
+
+ /// @brief Constructor from ostream
+ WriterCompressedAscii(std::ostream& stream,
+ shared_ptr<GenRunInfo> run = shared_ptr<GenRunInfo>());
+
+ /// @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<ConstGenParticlePtr> 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<long> m_stripid; //!< Strip matching intermediate particles
+ map<long,long> 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 <cstring>
+#include <sstream>
+
+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: "<<filename )
+ }
+ m_stream = &m_file;
+ set_run_info(make_shared<GenRunInfo>());
+}
+
+
+// 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<GenRunInfo>());
+}
+
+
+
+ReaderCompressedAscii::~ReaderCompressedAscii() { }
+
+
+bool ReaderCompressedAscii::read_event(GenEvent &evt) {
+
+ if ( failed() ) return false;
+
+ bool parsed_event_header = false;
+ bool is_parsing_successful = true;
+ pair<int,int> 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_pvx.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;
+ }
+
+ m_vertices.erase(0);
+ for ( auto it : m_vertices ) {
+ if ( it.second ) m_evt->add_vertex(it.second);
+ }
+
+ for ( int ip = 0, Np = m_particles.size(); ip < Np; ++ip ) {
+ GenParticlePtr p = m_particles[ip];
+ int ivp = m_pvx[ip].first;
+ GenVertexPtr vp = m_vertices[-ivp];
+ if ( ivp == 0 )
+ m_evt->add_beam_particle(p);
+ else if ( vp )
+ vp->add_particle_out(p);
+ else
+ cout << "WARNING did not find production vertex for particle "
+ << ip + 1 << endl;
+
+ int ive = m_pvx[ip].second;
+ GenVertexPtr ve = m_vertices[-ive];
+ if ( ve ) ve->add_particle_in(p);
+ }
+
+
+
+ // 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<int,int> ReaderCompressedAscii::parse_event_information() {
+ static const pair<int,int> err(-1,-1);
+ pair<int,int> 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: "<<event_no<<" ("<<ret.first<<"V, "<<ret.second<<"P)" )
+
+ return ret;
+}
+
+
+bool ReaderCompressedAscii::parse_weight_values() {
+
+ vector<double> 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<GenVertex>();
+
+ int id = 0;
+ if ( !(is >> id) ) return false;
+
+ int status = 0;
+ if ( !(is >> status) ) return false;
+ data->set_status( status );
+
+ if ( !read_position(data) ) return false;
+
+ m_vertices[-id] = data;
+
+ return true;
+}
+
+
+bool ReaderCompressedAscii::parse_particle_information() {
+ GenParticlePtr data = make_shared<GenParticle>();
+
+ int id = 0;
+ if ( !(is >> id) ) return false;
+
+ int ivp = 0, ive = 0;
+ if ( !(is >> ivp >> ive) ) 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_pvx.push_back(make_pair(ivp, ive));
+
+ 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<Attribute> att =
+ make_shared<StringAttribute>(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<Attribute> att =
+ make_shared<StringAttribute>(StringAttribute(unescape(contents)));
+
+ run_info()->add_attribute(name, att);
+
+ return true;
+
+}
+
+
+bool ReaderCompressedAscii::parse_weight_names() {
+
+ vector<string> 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,108 +1,106 @@
// -*- 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<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 ) 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){
- io->read_event(*evt);
- return !io->failed();
+ return io->read_event(*evt) && !io->failed();
}
shared_ptr<HepMC_IO_type> makeReader(std::istream & istr,
std::string * errm) {
shared_ptr<HepMC_IO_type> 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<HepMC_IO_type>();
+ 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<RivetHepMC::ReaderAscii>(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<Rivet::ReaderCompressedAscii>(istr);
+ else
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>();
}
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,366 @@
+// -*- 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 <cstring>
+#include <algorithm>//min max for VS2017
+
+namespace Rivet {
+
+using namespace HepMC3;
+
+
+WriterCompressedAscii::WriterCompressedAscii(const std::string &filename, shared_ptr<GenRunInfo> 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: "<<filename )
+ } else {
+ m_file << "HepMC::Version " << version() << std::endl;
+ m_file << "HepMC::Asciiv3-START_EVENT_LISTING" << std::endl;
+ if ( run_info() ) write_run_info();
+ }
+}
+
+
+WriterCompressedAscii::WriterCompressedAscii(std::ostream &stream, shared_ptr<GenRunInfo> 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;
+ strip(e);
+ set<long> 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: "<<vt1.first )
+ }
+ else {
+ os << "A " << vt2.first
+ << " " << vt1.first
+ << " " << escape(st) << endl;
+ }
+ }
+ }
+
+ set<ConstGenVertexPtr> 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 ( 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();
+
+ 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<GenRunInfo>());
+
+ vector<string> 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();
+ ConstGenVertexPtr ve = p->end_vertex();
+
+ os << "P " << p->id()
+ << " " << (vp? vp->id(): 0)
+ << " " << (ve? ve->id(): 0)
+ << " " << p->pid();
+ write_momentum(p->momentum());
+ write_mass(p);
+ os << " " << p->status() << endl;
+
+}
+
+void WriterCompressedAscii::close() {
+ std::ofstream* ofs = dynamic_cast<std::ofstream*>(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<GenParticlePtr> 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);
+ }
+}
+
+}