//==========================================================================
// AIDA Detector description implementation
//--------------------------------------------------------------------------
// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
// All rights reserved.
//
// For the licensing terms see $DD4hepINSTALL/LICENSE.
// For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
//
// Author : M.Frank
//
//==========================================================================
/** \addtogroup Geant4EventReader
*
@{
\package Geant4EventReaderHepMC
* \brief Plugin to read HepMC2 ASCII files
*
*
@}
*/
// Framework include files
#include <DDG4/IoStreams.h>
#include <DDG4/Geant4InputAction.h>
// C/C++ include files
/// Namespace for the AIDA detector description toolkit
namespace dd4hep {
/// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
namespace sim {
/// HepMC namespace declaration
namespace HepMC {
/// HepMC EventStream class used internally by the Geant4EventReaderHepMC plugin
class EventStream;
}
/// Class to populate Geant4 primaries from HepMC(2) files.
/**
* Class to populate Geant4 primary particles and vertices from a
* file in HepMC format (ASCII)
*
* For details also see:
* http://hepmc.web.cern.ch/hepmc/ReaderAsciiHepMC2_8cc_source.html
*
* \author P.Kostka (main author)
* \author M.Frank (code reshuffeling into new DDG4 scheme)
* \version 1.0
* \ingroup DD4HEP_SIMULATION
*/
class Geant4EventReaderHepMC : public Geant4EventReader {
typedef boost::iostreams::stream<dd4hep_file_source<int> > in_stream;
//typedef boost::iostreams::stream<dd4hep_file_source<TFile*> > in_stream;
typedef HepMC::EventStream EventStream;
protected:
in_stream m_input;
EventStream* m_events;
public:
/// Initializing constructor
explicit Geant4EventReaderHepMC(const std::string& nam);
/// Default destructor
virtual ~Geant4EventReaderHepMC();
/// Read an event and fill a vector of MCParticles.
virtual EventReaderStatus readParticles(int event_number,
Vertices& vertices,
std::vector<Particle*>& particles) override;
virtual EventReaderStatus moveToEvent(int event_number) override;
virtual EventReaderStatus skipEvent() override { return EVENT_READER_OK; }
};
} /* End namespace sim */
} /* End namespace dd4hep */
//====================================================================
// AIDA Detector description implementation
//--------------------------------------------------------------------
//
//====================================================================
// #include <DDG4/Geant4EventReaderHepMC"
// Framework include files
#include <DDG4/Factories.h>
#include <DD4hep/Printout.h>
#include <DDG4/Geant4Primary.h>
#include <CLHEP/Units/SystemOfUnits.h>
#include <CLHEP/Units/PhysicalConstants.h>
// C/C++ include files
#include <cerrno>
#include <climits>
#include <algorithm>
using namespace dd4hep::sim;
using PropertyMask = dd4hep::detail::ReferenceBitMask<int>;
// Factory entry
DECLARE_GEANT4_EVENT_READER(Geant4EventReaderHepMC)
/// Namespace for the AIDA detector description toolkit
namespace dd4hep {
/// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
namespace sim {
/// HepMC namespace declaration
namespace HepMC {
/// HepMC EventHeader class used internally by the Geant4EventReaderHepMC plugin
/*
* \author P.Kostka (main author)
* \author M.Frank (code reshuffeling into new DDG4 scheme)
* \version 1.0
* \ingroup DD4HEP_SIMULATION
*/
class EventHeader {
public:
int id;
int num_vertices;
int bp1, bp2;
int signal_process_id;
int signal_process_vertex;
float scale;
float alpha_qcd;
float alpha_qed;
std::vector<float> weights;
std::vector<long> random;
/// Default constructor
EventHeader() : id(0), num_vertices(0), bp1(0), bp2(0),
signal_process_id(0), signal_process_vertex(0),
scale(0.0), alpha_qcd(0.0), alpha_qed(0.0), weights(), random() {}
};
/// The known_io enum is used to track which type of input is being read
enum known_io { gen=1, ascii, extascii, ascii_pdt, extascii_pdt };
/// HepMC EventStream class used internally by the Geant4EventReaderHepMC plugin
/*
* \author P.Kostka (main author)
* \author M.Frank (code reshuffeling into new DDG4 scheme)
* \version 1.0
* \ingroup DD4HEP_SIMULATION
*/
class EventStream {
public:
typedef std::map<int,Geant4Vertex*> Vertices;
typedef std::map<int,Geant4Particle*> Particles;
std::istream& instream;
// io information
std::string key;
double mom_unit, pos_unit;
int io_type;
float xsection, xsection_err;
EventHeader header;
Vertices m_vertices;
Particles m_particles;
/// Default constructor
EventStream(std::istream& in) : instream(in), mom_unit(0.0), pos_unit(0.0),
io_type(0), xsection(0.0), xsection_err(0.0)
{ use_default_units(); }
/// Check if data stream is in proper state and has data
bool ok() const;
Geant4Vertex* vertex(int i);
Particles& particles() { return m_particles; }
Vertices& vertices() { return m_vertices; }
void set_io(int typ, const std::string& k)
{ io_type = typ; key = k; }
void use_default_units()
{ mom_unit = CLHEP::MeV; pos_unit = CLHEP::mm; }
bool read();
void clear();
};
char get_input(std::istream& is, std::istringstream& iline);
int read_until_event_end(std::istream & is);
int read_weight_names(EventStream &, std::istringstream& iline);
int read_particle(EventStream &info, std::istringstream& iline, Geant4Particle * p);
int read_vertex(EventStream &info, std::istream& is, std::istringstream & iline);
int read_event_header(EventStream &info, std::istringstream & input, EventHeader& header);
int read_cross_section(EventStream &info, std::istringstream & input);
int read_units(EventStream &info, std::istringstream & input);
int read_heavy_ion(EventStream &, std::istringstream & input);
int read_pdf(EventStream &, std::istringstream & input);
Geant4Vertex* vertex(EventStream& info, int i);
void fix_particles(EventStream &info);
}
}
}
// For debugging specify here vertex numbers
//#define DD4HEP_DEBUG_HEP_MC_VERTEX -390
//#define DD4HEP_DEBUG_HEP_MC_PARTICLE 418
/// Initializing constructor
Geant4EventReaderHepMC::Geant4EventReaderHepMC(const std::string& nam)
: Geant4EventReader(nam), m_input(), m_events(0)
{
// Now open the input file:
m_input.open(nam.c_str(), BOOST_IOS::in|BOOST_IOS::binary);
if ( not m_input.is_open() ) {
except("+++ Failed to open input stream: %s Error:%s.", nam.c_str(), ::strerror(errno));
}
m_events = new HepMC::EventStream(m_input);
}
/// Default destructor
Geant4EventReaderHepMC::~Geant4EventReaderHepMC() {
delete m_events;
m_events = 0;
m_input.close();
}
/// skipEvents if required
Geant4EventReader::EventReaderStatus
Geant4EventReaderHepMC::moveToEvent(int event_number) {
if( m_currEvent < event_number && event_number != 0 ) {
printout(INFO,"EventReaderHepMC::moveToEvent","Current event:%d Skipping the next %d events",
m_currEvent, event_number);
while ( m_currEvent < event_number ) {
if ( not m_events->read() ) return EVENT_READER_ERROR;
++m_currEvent;
}
}
printout(DEBUG,"EventReaderHepMC::moveToEvent","Current event number: %d",m_currEvent);
return EVENT_READER_OK;
}
/// Read an event and fill a vector of MCParticles.
Geant4EventReaderHepMC::EventReaderStatus
Geant4EventReaderHepMC::readParticles(int /* ev_id */,
Vertices& vertices,
Particles& output) {
//fg: for now we create exactly one event vertex here ( as before )
// this needs revisiting as HepMC allows to have more than one vertex ...
Geant4Vertex* primary_vertex = new Geant4Vertex ;
vertices.emplace_back( primary_vertex );
primary_vertex->x = 0;
primary_vertex->y = 0;
primary_vertex->z = 0;
if ( !m_events->ok() ) {
vertices.clear();
output.clear();
return EVENT_READER_EOF;
}
else if ( m_events->read() ) {
EventStream::Particles& parts = m_events->particles();
Position pos(primary_vertex->x,primary_vertex->y,primary_vertex->z);
output.reserve(parts.size());
transform(parts.begin(),parts.end(),back_inserter(output),detail::reference2nd(parts));
m_events->clear();
if (pos.mag2() > std::numeric_limits<double>::epsilon() ) {
for(Particles::iterator k=output.begin(); k != output.end(); ++k) {
Geant4ParticleHandle p(*k);
p->vsx += pos.x();
p->vsy += pos.y();
p->vsz += pos.z();
p->vex += pos.x();
p->vey += pos.y();
p->vez += pos.z();
}
}
for(Particles::const_iterator k=output.begin(); k != output.end(); ++k) {
Geant4ParticleHandle p(*k);
printout(VERBOSE,m_name,
"+++ %s ID:%3d status:%08X typ:%9d Mom:(%+.2e,%+.2e,%+.2e)[MeV] "
"time: %+.2e [ns] #Dau:%3d #Par:%1d",
"",p->id,p->status,p->pdgID,
p->psx/CLHEP::MeV,p->psy/CLHEP::MeV,p->psz/CLHEP::MeV,p->time/CLHEP::ns,
p->daughters.size(),
p->parents.size());
//output.emplace_back(p);
//add particles to the 'primary vertex'
if ( p->parents.size() == 0 ) {
PropertyMask status(p->status);
if ( status.isSet(G4PARTICLE_GEN_EMPTY) || status.isSet(G4PARTICLE_GEN_DOCUMENTATION) )
primary_vertex->in.insert(p->id); // Beam particles and primary quarks etc.
else
primary_vertex->out.insert(p->id); // Stuff, to be given to Geant4 together with daughters
}
}
++m_currEvent;
return EVENT_READER_OK;
}
vertices.clear();
output.clear();
return EVENT_READER_EOF;
}
void HepMC::fix_particles(EventStream& info) {
EventStream::Particles& parts = info.particles();
EventStream::Vertices& verts = info.vertices();
EventStream::Particles::iterator i;
std::set<int>::const_iterator id, ip;
for(i=parts.begin(); i != parts.end(); ++i) {
Geant4ParticleHandle p((*i).second);
int end_vtx_id = p->secondaries;
p->secondaries = 0;
Geant4Vertex* v = vertex(info,end_vtx_id);
#if defined(DD4HEP_DEBUG_HEP_MC_VERTEX)
if ( end_vtx_id == DD4HEP_DEBUG_HEP_MC_VERTEX ) {
std::cout << "End-vertex:" << end_vtx_id << std::endl;
}
#endif
if ( v ) {
p->vex = v->x;
p->vey = v->y;
p->vez = v->z;
v->in.insert(p->id);
for(id=v->out.begin(); id!=v->out.end();++id) {
EventStream::Particles::iterator ipp = parts.find(*id);
Geant4Particle* dau = ipp != parts.end() ? (*ipp).second : 0;
if ( !dau )
std::cout << "ERROR: Invalid daughter particle: " << *id << std::endl;
else
dau->parents.insert(p->id);
p->daughters.insert(*id);
}
}
}
for(const auto& iv : verts) {
Geant4Vertex* v = iv.second;
for (int pout : v->out) {
EventStream::Particles::iterator ipp = parts.find(pout);
if ( ipp != parts.end() ) {
Geant4Particle* p = (*ipp).second;
for (int d : v->in) {
p->parents.insert(d);
}
}
}
}
/// Particles originating from the beam (=no parents) must be
/// be stripped off their parents and the status set to G4PARTICLE_GEN_DECAYED!
std::vector<Geant4Particle*> beam;
for(const auto& ipart : parts) {
Geant4ParticleHandle p(ipart.second);
if ( p->parents.size() == 0 ) {
for(int d : p->daughters) {
Geant4Particle *pp = parts[d];
beam.emplace_back(pp);
}
}
}
for(auto* ipp : beam) {
//std::cout << "Clear parents of " << (*ipp)->id << std::endl;
ipp->parents.clear();
ipp->status = G4PARTICLE_GEN_DECAYED;
}
}
Geant4Vertex* HepMC::vertex(EventStream& info, int i) {
EventStream::Vertices::iterator it=info.vertices().find(i);
return (it==info.vertices().end()) ? 0 : (*it).second;
}
char HepMC::get_input(std::istream& is, std::istringstream& iline) {
char value = is.peek();
if ( !is ) { // make sure the stream is valid
std::cerr << "StreamHelpers: setting badbit." << std::endl;
is.clear(std::ios::badbit);
return -1;
}
std::string line, firstc;
getline(is,line);
if ( !is ) { // make sure the stream is valid
std::cerr << "StreamHelpers: setting badbit." << std::endl;
is.clear(std::ios::badbit);
return -1;
}
iline.clear();
iline.str(line);
if ( line.length()>0 && line[1] == ' ' ) iline >> firstc;
return iline ? value : -1;
}
int HepMC::read_until_event_end(std::istream & is) {
std::string line;
while ( is ) {
char val = is.peek();
if( val == 'E' ) { // next event
return 1;
}
getline(is,line);
if ( !is.good() ) return 0;
}
return 0;
}
int HepMC::read_weight_names(EventStream&, std::istringstream&) {
#if 0
int HepMC::read_weight_names(EventStream& info, std::istringstream& iline)
size_t name_size = 0;
iline >> name_size;
info.weights.names.clear();
info.weights.weights.clear();
std::string name;
WeightContainer namedWeight;
std::string::size_type i1 = line.find("\""), i2, len = line.size();
for(size_t ii = 0; ii < name_size; ++ii) {
// weight names may contain blanks
if(i1 >= len) {
std::cout << "debug: attempting to read past the end of the named weight line " << std::endl;
std::cout << "debug: We should never get here" << std::endl;
std::cout << "debug: Looking for the end of this event" << std::endl;
read_until_event_end(is);
}
i2 = line.find("\"",i1+1);
name = line.substr(i1+1,i2-i1-1);
if ( namedWeight.names.find(name) == namedWeight.names.end() )
namedWeight.names[name] = namedWeight.names.size();
namedWeight.weights[ii] = info.weights.weights[ii];
i1 = line.find("\"",i2+1);
}
info.weights.weights = namedWeight.weights;
info.weights.names = namedWeight.names;
#endif
return 1;
}
int HepMC::read_particle(EventStream &info, std::istringstream& input, Geant4Particle * p) {
float ene = 0., theta = 0., phi = 0;
int size = 0, stat=0;
PropertyMask status(p->status);
// check that the input stream is still OK after reading item
input >> p->id >> p->pdgID >> p->psx >> p->psy >> p->psz >> ene;
p->id = info.particles().size();
#if defined(DD4HEP_DEBUG_HEP_MC_PARTICLE)
if ( p->id == DD4HEP_DEBUG_HEP_MC_PARTICLE ) {
std::cout << "Particle id: " << p->id << std::endl;
}
#endif
p->charge = 0;
p->psx *= info.mom_unit;
p->psy *= info.mom_unit;
p->psz *= info.mom_unit;
ene *= info.mom_unit;
if ( !input )
return 0;
else if ( info.io_type != ascii ) {
input >> p->mass;
p->mass *= info.mom_unit;
}
else {
p->mass = std::sqrt(fabs(ene*ene - (p->psx*p->psx + p->psy*p->psy + p->psz*p->psz)));
}
// Reuse here the secondaries to store the end-vertex ID
input >> stat >> theta >> phi >> p->secondaries >> size;
if( !input ) {
return 0;
}
//
// Generator status
// Simulator status 0 until simulator acts on it
status.clear();
if ( stat == 0 ) status.set(G4PARTICLE_GEN_EMPTY);
else if ( stat == 0x1 ) status.set(G4PARTICLE_GEN_STABLE);
else if ( stat == 0x2 ) status.set(G4PARTICLE_GEN_DECAYED);
else if ( stat == 0x3 ) status.set(G4PARTICLE_GEN_DOCUMENTATION);
else if ( stat == 0x4 ) status.set(G4PARTICLE_GEN_DOCUMENTATION);
else if ( stat == 0xB ) status.set(G4PARTICLE_GEN_DOCUMENTATION);
else status.set(G4PARTICLE_GEN_OTHER);
/// If there is an end vertex, the particle already decayed
if ( p->secondaries != 0 ) {
status.set(G4PARTICLE_GEN_DECAYED);
}
/// Keep a copy of the full generator status
p->genStatus = stat&G4PARTICLE_GEN_STATUS_MASK;
// read flow patterns if any exist. Protect against tainted readings.
size = std::min(size,100);
for (int i = 0; i < size; ++i ) {
input >> p->colorFlow[0] >> p->colorFlow[1];
if(!input) return 0;
}
return 1;
}
int HepMC::read_vertex(EventStream &info, std::istream& is, std::istringstream & input) {
int id=0, dummy = 0, num_orphans_in=0, num_particles_out=0, weights_size=0;
std::vector<float> weights;
Geant4Vertex* v = new Geant4Vertex();
Geant4Particle* p;
if( !input ) {
delete v;
return 0;
}
input >> id >> dummy >> v->x >> v->y >> v->z >> v->time
>> num_orphans_in >> num_particles_out >> weights_size;
if( !input ) {
delete v;
return 0;
}
if ( weights_size < 0 || weights_size > USHRT_MAX ) {
delete v;
return 0;
}
#if defined(DD4HEP_DEBUG_HEP_MC_VERTEX)
if ( id == DD4HEP_DEBUG_HEP_MC_VERTEX ) {
printout(ALWAYS,"HepMC","++ Created Vertex ID=%d",id);
}
#endif
v->x *= info.pos_unit;
v->y *= info.pos_unit;
v->z *= info.pos_unit;
for (int i1 = 0; i1 < weights_size; ++i1) {
float value = 0e0;
input >> value;
weights.emplace_back(value);
if( !input ) {
delete v;
return 0;
}
}
info.vertices().emplace(id,v);
for(char value = is.peek(); value=='P'; value=is.peek()) {
value = get_input(is,input);
if( !input || value < 0 )
return 0;
read_particle(info, input,p = new Geant4Particle());
if( !input ) {
printout(ERROR,"HepMC","++ Vertex %d Failed to daughter read particle!",id);
delete p;
return 0;
}
info.particles().emplace(p->id,p);
p->pex = p->psx;
p->pey = p->psy;
p->pez = p->psz;
if ( --num_orphans_in >= 0 ) {
v->in.insert(p->id);
p->vex = v->x;
p->vey = v->y;
p->vez = v->z;
#if defined(DD4HEP_DEBUG_HEP_MC_VERTEX)
if ( id == DD4HEP_DEBUG_HEP_MC_VERTEX ) {
printout(ALWAYS,"HepMC","++ Vertex %d Add INGOING Particle: %d",id,p->id);
}
#endif
}
else if ( num_particles_out >= 0 ) {
v->out.insert(p->id);
p->vsx = v->x;
p->vsy = v->y;
p->vsz = v->z;
#if defined(DD4HEP_DEBUG_HEP_MC_VERTEX)
if ( id == DD4HEP_DEBUG_HEP_MC_VERTEX ) {
printout(ALWAYS,"HepMC","++ Vertex %d Add OUTGOING Particle: %d",id,p->id);
}
#endif
}
else {
delete p;
except("HepMC", "Invalid number of particles....");
}
}
return 1;
}
int HepMC::read_event_header(EventStream &info, std::istringstream & input, EventHeader& header) {
// read values into temp variables, then fill GenEvent
int size = 0;
input >> header.id;
if( info.io_type == gen || info.io_type == extascii ) {
int nmpi = -1;
input >> nmpi;
if( input.fail() ) return 0;
//MSF set_mpi( nmpi );
}
input >> header.scale;
input >> header.alpha_qcd;
input >> header.alpha_qed;
input >> header.signal_process_id;
input >> header.signal_process_vertex;
input >> header.num_vertices;
if( info.io_type == gen || info.io_type == extascii )
input >> header.bp1 >> header.bp2;
printout(DEBUG,"HepMC","++ Event header: %s",input.str().c_str());
input >> size;
input.clear();
if( input.fail() )
return 0;
if( size < 0 || size > USHRT_MAX )
return 0;
for(int i = 0; i < size; ++i ) {
long val = 0e0;
input >> val;
header.random.emplace_back(val);
if( input.fail() ) return 0;
}
input >> size;
if( input.fail() )
return 0;
if( size < 0 || size > USHRT_MAX )
return 0;
std::vector<float> wgt;
for( int ii = 0; ii < size; ++ii ) {
float val = 0e0;
input >> val;
wgt.emplace_back(val);
if( input.fail() ) return 0;
}
// weight names will be added later if they exist
if( !wgt.empty() ) header.weights = std::move(wgt);
return 1;
}
int HepMC::read_cross_section(EventStream &info, std::istringstream & input) {
input >> info.xsection >> info.xsection_err;
return input.fail() ? 0 : 1;
}
int HepMC::read_units(EventStream &info, std::istringstream & input) {
if( info.io_type == gen ) {
std::string mom, pos;
input >> mom >> pos;
if ( !input.fail() ) {
if ( mom == "KEV" ) info.mom_unit = CLHEP::keV;
else if ( mom == "MEV" ) info.mom_unit = CLHEP::MeV;
else if ( mom == "GEV" ) info.mom_unit = CLHEP::GeV;
else if ( mom == "TEV" ) info.mom_unit = CLHEP::TeV;
if ( pos == "MM" ) info.pos_unit = CLHEP::mm;
else if ( pos == "CM" ) info.pos_unit = CLHEP::cm;
else if ( pos == "M" ) info.pos_unit = CLHEP::m;
}
}
return input.fail() ? 0 : 1;
}
int HepMC::read_heavy_ion(EventStream &, std::istringstream & input) {
// read values into temp variables, then create a new HeavyIon object
int nh =0, np =0, nt =0, nc =0,
neut = 0, prot = 0, nw =0, nwn =0, nwnw =0;
float impact = 0., plane = 0., xcen = 0., inel = 0.;
input >> nh >> np >> nt >> nc >> neut >> prot >> nw >> nwn >> nwnw;
input >> impact >> plane >> xcen >> inel;
/*
std::cerr << "Reading heavy ion, but igoring data!" << std::endl;
ion->set_Ncoll_hard(nh);
ion->set_Npart_proj(np);
ion->set_Npart_targ(nt);
ion->set_Ncoll(nc);
ion->set_spectator_neutrons(neut);
ion->set_spectator_protons(prot);
ion->set_N_Nwounded_collisions(nw);
ion->set_Nwounded_N_collisions(nwn);
ion->set_Nwounded_Nwounded_collisions(nwnw);
ion->set_impact_parameter(impact);
ion->set_event_plane_angle(plane);
ion->set_eccentricity(xcen);
ion->set_sigma_inel_NN(inel);
*/
return input.fail() ? 0 : 1;
}
int HepMC::read_pdf(EventStream &, std::istringstream & input) {
// read values into temp variables, then create a new PdfInfo object
int id1 =0, id2 =0;
double x1 = 0., x2 = 0., scale = 0., pdf1 = 0., pdf2 = 0.;
input >> id1 ;
if ( input.fail() )
return 0;
// check now for empty PdfInfo line
if( id1 == 0 )
return 0;
// continue reading
input >> id2 >> x1 >> x2 >> scale >> pdf1 >> pdf2;
if ( input.fail() )
return 0;
/*
std::cerr << "Reading pdf, but igoring data!" << std::endl;
pdf->set_id1( id1 );
pdf->set_id2( id2 );
pdf->set_x1( x1 );
pdf->set_x2( x2 );
pdf->set_scalePDF( scale );
pdf->set_pdf1( pdf1 );
pdf->set_pdf2( pdf2 );
*/
// check to see if we are at the end of the line
if( !input.eof() ) {
int pdf_id1=0, pdf_id2=0;
input >> pdf_id1 >> pdf_id2;
/*
pdf->set_pdf_id1( pdf_id1 );
pdf->set_pdf_id2( pdf_id2 );
*/
}
return input.fail() ? 0 : 1;
}
/// Check if data stream is in proper state and has data
bool HepMC::EventStream::ok() const {
// make sure the stream is good
if ( instream.eof() || instream.fail() ) {
instream.clear(std::ios::badbit);
return false;
}
return true;
}
void HepMC::EventStream::clear() {
detail::releaseObjects(m_vertices);
detail::releaseObjects(m_particles);
}
bool HepMC::EventStream::read() {
EventStream& info = *this;
bool event_read = false;
detail::releaseObjects(vertices());
detail::releaseObjects(particles());
while( instream.good() ) {
char value = instream.peek();
std::istringstream input_line;
if ( value == 'E' && event_read )
break;
else if ( instream.eof() && event_read )
goto Done;
else if ( instream.eof() )
return false;
else if ( value=='#' || ::isspace(value) ) {
get_input(instream,input_line);
continue;
}
value = get_input(instream,input_line);
// On failure switch to end
if( !input_line || value < 0 )
goto Skip;
switch( value ) {
case 'H': {
int iotype = 0;
std::string key_value;
input_line >> key_value;
// search for event listing key before first event only.
key_value = key_value.substr(0,key_value.find('\r'));
if ( key_value == "H" && (this->io_type == gen || this->io_type == extascii) ) {
read_heavy_ion(info, input_line);
break;
}
else if( key_value == "HepMC::IO_GenEvent-START_EVENT_LISTING" )
this->set_io(gen,key_value);
else if( key_value == "HepMC::IO_Ascii-START_EVENT_LISTING" )
this->set_io(ascii,key_value);
else if( key_value == "HepMC::IO_ExtendedAscii-START_EVENT_LISTING" )
this->set_io(extascii,key_value);
else if( key_value == "HepMC::IO_Ascii-START_PARTICLE_DATA" )
this->set_io(ascii_pdt,key_value);
else if( key_value == "HepMC::IO_ExtendedAscii-START_PARTICLE_DATA" )
this->set_io(extascii_pdt,key_value);
else if( key_value == "HepMC::IO_GenEvent-END_EVENT_LISTING" )
iotype = gen;
else if( key_value == "HepMC::IO_Ascii-END_EVENT_LISTING" )
iotype = ascii;
else if( key_value == "HepMC::IO_ExtendedAscii-END_EVENT_LISTING" )
iotype = extascii;
else if( key_value == "HepMC::IO_Ascii-END_PARTICLE_DATA" )
iotype = ascii_pdt;
else if( key_value == "HepMC::IO_ExtendedAscii-END_PARTICLE_DATA" )
iotype = extascii_pdt;
if( iotype != 0 && this->io_type != iotype ) {
std::cerr << "GenEvent::find_end_key: iotype keys have changed. "
<< "MALFORMED INPUT" << std::endl;
instream.clear(std::ios::badbit);
return false;
}
else if ( iotype != 0 ) {
break;
}
continue;
}
case 'E': // deal with the event line
if ( !read_event_header(info, input_line, this->header) )
goto Skip;
event_read = true;
continue;
case 'N': // get weight names
if ( !read_weight_names(info, input_line) )
goto Skip;
continue;
case 'U': // get unit information if it exists
if ( !read_units(info, input_line) )
goto Skip;
continue;
case 'C': // we have a GenCrossSection line
if ( !read_cross_section(info, input_line) )
goto Skip;
continue;
case 'V': // Read vertex with particles
if ( !read_vertex(info, instream, input_line) )
goto Skip;
continue;
case 'F': // Read PDF
if ( !read_pdf(info, input_line) )
goto Skip;
continue;
case 'P': // we should not find this line
std::cerr << "streaming input: found unexpected Particle line." << std::endl;
continue;
default: // ignore everything else
continue;
}
continue;
Skip:
printout(WARNING,"HepMC::EventStream","+++ Skip event with ID: %d",this->header.id);
detail::releaseObjects(vertices());
detail::releaseObjects(particles());
read_until_event_end(instream);
event_read = false;
if ( instream.eof() ) return false;
}
if( not instream.good() ) return false;
Done:
fix_particles(info);
detail::releaseObjects(vertices());
return true;
}