diff --git a/include/Rivet/Tools/ReaderCompressedAscii.hh b/include/Rivet/Tools/ReaderCompressedAscii.hh --- a/include/Rivet/Tools/ReaderCompressedAscii.hh +++ b/include/Rivet/Tools/ReaderCompressedAscii.hh @@ -1,163 +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< pair > m_pvx; + 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/src/Tools/ReaderCompressedAscii.cc b/src/Tools/ReaderCompressedAscii.cc --- a/src/Tools/ReaderCompressedAscii.cc +++ b/src/Tools/ReaderCompressedAscii.cc @@ -1,464 +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_pvx.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; } - m_vertices.erase(0); - for ( auto it : m_vertices ) { - if ( it.second ) m_evt->add_vertex(it.second); - } + // 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]); - 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; + // 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]); - int ive = m_pvx[ip].second; - GenVertexPtr ve = m_vertices[-ive]; - if ( ve ) ve->add_particle_in(p); - } - - + // 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, ive = 0; - if ( !(is >> ivp >> ive) ) 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_pvx.push_back(make_pair(ivp, ive)); + 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/WriterCompressedAscii.cc b/src/Tools/WriterCompressedAscii.cc --- a/src/Tools/WriterCompressedAscii.cc +++ b/src/Tools/WriterCompressedAscii.cc @@ -1,366 +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 ( done.insert(v).second ) + 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(); + 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(); - 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(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); } } }