HepMC3 event record library
ReaderAscii.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // This file is part of HepMC
4 // Copyright (C) 2014-2019 The HepMC collaboration (see AUTHORS for details)
5 //
6 ///
7 /// @file ReaderAscii.cc
8 /// @brief Implementation of \b class ReaderAscii
9 ///
10 #include "HepMC3/ReaderAscii.h"
11 
12 #include "HepMC3/GenEvent.h"
13 #include "HepMC3/GenParticle.h"
14 #include "HepMC3/GenVertex.h"
15 #include "HepMC3/Units.h"
16 #include <cstring>
17 #include <sstream>
18 
19 namespace HepMC3 {
20 
21 
22 ReaderAscii::ReaderAscii(const string &filename)
23  : m_file(filename), m_stream(0), m_isstream(false)
24 {
25  if( !m_file.is_open() ) {
26  ERROR( "ReaderAscii: could not open input file: "<<filename )
27  }
28  set_run_info(make_shared<GenRunInfo>());
29 }
30 
31 
32 // Ctor for reading from stdin
33 ReaderAscii::ReaderAscii(std::istream & stream)
34  : m_stream(&stream), m_isstream(true)
35 {
36  if( !m_stream ) {
37  ERROR( "ReaderAscii: could not open input stream " )
38  }
39  set_run_info(make_shared<GenRunInfo>());
40 }
41 
42 
43 
44 ReaderAscii::~ReaderAscii() { if (!m_isstream) close(); }
45 
46 
48  if ( (!m_file.is_open()) && (!m_isstream) ) return false;
49 
50  char peek;
51  char buf[512*512];
52  bool parsed_event_header = false;
53  bool is_parsing_successful = true;
54  pair<int,int> vertices_and_particles(0,0);
55 
56  evt.clear();
57  evt.set_run_info(run_info());
58 
59  //
60  // Parse event, vertex and particle information
61  //
62  while(!failed()) {
63 
64  m_isstream ? m_stream->getline(buf,512*512) : m_file.getline(buf,512*512);
65 
66  if( strlen(buf) == 0 ) continue;
67 
68  // Check for ReaderAscii header/footer
69  if( strncmp(buf,"HepMC",5) == 0 ) {
70  if( strncmp(buf,"HepMC::Version",14) != 0 && strncmp(buf,"HepMC::Asciiv3",14)!=0 )
71  {
72  WARNING( "ReaderAscii: found unsupported expression in header. Will close the input." )
73  std::cout<<buf<<std::endl;
74  m_isstream ? m_stream->clear(ios::eofbit) : m_file.clear(ios::eofbit);
75  }
76  if(parsed_event_header) {
77  is_parsing_successful = true;
78  break;
79  }
80  continue;
81  }
82 
83  switch(buf[0]) {
84  case 'E':
85  vertices_and_particles = parse_event_information(evt,buf);
86  if (vertices_and_particles.second < 0) {
87  is_parsing_successful = false;
88  } else {
89  is_parsing_successful = true;
90  parsed_event_header = true;
91  }
92  break;
93  case 'V':
94  is_parsing_successful = parse_vertex_information(evt,buf);
95  break;
96  case 'P':
97  is_parsing_successful = parse_particle_information(evt,buf);
98  break;
99  case 'W':
100  if ( parsed_event_header )
101  is_parsing_successful = parse_weight_values(evt,buf);
102  else
103  is_parsing_successful = parse_weight_names(buf);
104  break;
105  case 'U':
106  is_parsing_successful = parse_units(evt,buf);
107  break;
108  case 'T':
109  is_parsing_successful = parse_tool(buf);
110  break;
111  case 'A':
112  if ( parsed_event_header )
113  is_parsing_successful = parse_attribute(evt,buf);
114  else
115  is_parsing_successful = parse_run_attribute(buf);
116  break;
117  default:
118  WARNING( "ReaderAscii: skipping unrecognised prefix: " << buf[0] )
119  is_parsing_successful = true;
120  break;
121  }
122 
123  if( !is_parsing_successful ) break;
124 
125  // Check for next event
126  m_isstream ? peek = m_stream->peek() : peek = m_file.peek();
127  if( parsed_event_header && peek=='E' ) break;
128  }
129 
130 
131  // Check if all particles and vertices were parsed
132  if ((int)evt.particles().size() > vertices_and_particles.second ) {
133  ERROR( "ReaderAscii: too many particles were parsed" )
134  printf("%zu vs %i expected\n",evt.particles().size(),vertices_and_particles.second );
135  is_parsing_successful = false;
136  }
137  if ((int)evt.particles().size() < vertices_and_particles.second ) {
138  ERROR( "ReaderAscii: too few particles were parsed" )
139  printf("%zu vs %i expected\n",evt.particles().size(),vertices_and_particles.second );
140  is_parsing_successful = false;
141  }
142 
143  if ((int)evt.vertices().size() > vertices_and_particles.first) {
144  ERROR( "ReaderAscii: too many vertices were parsed" )
145  printf("%zu vs %i expected\n",evt.vertices().size(),vertices_and_particles.first );
146  is_parsing_successful = false;
147  }
148 
149  if ((int)evt.vertices().size() < vertices_and_particles.first) {
150  ERROR( "ReaderAscii: too few vertices were parsed" )
151  printf("%zu vs %i expected\n",evt.vertices().size(),vertices_and_particles.first );
152  is_parsing_successful = false;
153  }
154  // Check if there were errors during parsing
155  if( !is_parsing_successful ) {
156  ERROR( "ReaderAscii: event parsing failed. Returning empty event" )
157  DEBUG( 1, "Parsing failed at line:" << endl << buf )
158 
159  evt.clear();
160  m_isstream ? m_stream->clear(ios::badbit) : m_file.clear(ios::badbit);
161 
162  return false;
163  }
164  return true;
165 }
166 
167 
168 pair<int,int> ReaderAscii::parse_event_information(GenEvent &evt, const char *buf) {
169  static const pair<int,int> err(-1,-1);
170  pair<int,int> ret(-1,-1);
171  const char *cursor = buf;
172  int event_no = 0;
173  FourVector position;
174 
175  // event number
176  if( !(cursor = strchr(cursor+1,' ')) ) return err;
177  event_no = atoi(cursor);
178  evt.set_event_number(event_no);
179 
180  // num_vertices
181  if( !(cursor = strchr(cursor+1,' ')) ) return err;
182  ret.first = atoi(cursor);
183 
184  // num_particles
185  if( !(cursor = strchr(cursor+1,' ')) ) return err;
186  ret.second = atoi(cursor);
187 
188  // check if there is position information
189  if( (cursor = strchr(cursor+1,'@')) ) {
190 
191  // x
192  if( !(cursor = strchr(cursor+1,' ')) ) return err;
193  position.setX(atof(cursor));
194 
195  // y
196  if( !(cursor = strchr(cursor+1,' ')) ) return err;
197  position.setY(atof(cursor));
198 
199  // z
200  if( !(cursor = strchr(cursor+1,' ')) ) return err;
201  position.setZ(atof(cursor));
202 
203  // t
204  if( !(cursor = strchr(cursor+1,' ')) ) return err;
205  position.setT(atof(cursor));
206  evt.shift_position_to( position );
207  }
208 
209  DEBUG( 10, "ReaderAscii: E: "<<event_no<<" ("<<ret.first<<"V, "<<ret.second<<"P)" )
210 
211  return ret;
212 }
213 
214 
215 bool ReaderAscii::parse_weight_values(GenEvent &evt, const char *buf) {
216 
217  std::istringstream iss(buf + 1);
218  vector<double> wts;
219  double w;
220  while ( iss >> w ) wts.push_back(w);
221  if ( run_info() && run_info()->weight_names().size()
222  && run_info()->weight_names().size() != wts.size() )
223  throw std::logic_error("ReaderAscii::parse_weight_values: "
224  "The number of weights ("+std::to_string((long long int)(wts.size()))+") does not match "
225  "the number weight names("+std::to_string((long long int)(run_info()->weight_names().size()))+") in the GenRunInfo object");
226  evt.weights() = wts;
227 
228  return true;
229 }
230 
231 
232 bool ReaderAscii::parse_units(GenEvent &evt, const char *buf) {
233  const char *cursor = buf;
234 
235  // momentum
236  if( !(cursor = strchr(cursor+1,' ')) ) return false;
237  ++cursor;
238  Units::MomentumUnit momentum_unit = Units::momentum_unit(cursor);
239 
240  // length
241  if( !(cursor = strchr(cursor+1,' ')) ) return false;
242  ++cursor;
243  Units::LengthUnit length_unit = Units::length_unit(cursor);
244 
245  evt.set_units(momentum_unit,length_unit);
246 
247  DEBUG( 10, "ReaderAscii: U: " << Units::name(evt.momentum_unit()) << " " << Units::name(evt.length_unit()) )
248 
249  return true;
250 }
251 
252 
253 bool ReaderAscii::parse_vertex_information(GenEvent &evt, const char *buf) {
254  GenVertexPtr data = make_shared<GenVertex>();
255  FourVector position;
256  const char *cursor = buf;
257  const char *cursor2 = nullptr;
258  int id = 0;
259  int highest_id = evt.particles().size();
260 
261  // id
262  if( !(cursor = strchr(cursor+1,' ')) ) return false;
263  id = atoi(cursor);
264 
265  // status
266  if( !(cursor = strchr(cursor+1,' ')) ) return false;
267  data->set_status( atoi(cursor) );
268 
269  // skip to the list of particles
270  if( !(cursor = strchr(cursor+1,'[')) ) return false;
271 
272  while(true) {
273  ++cursor; // skip the '[' or ',' character
274  cursor2 = cursor; // save cursor position
275  int particle_in = atoi(cursor);
276 
277  // add incoming particle to the vertex
278  if( particle_in > 0 && particle_in <= highest_id) {
279  data->add_particle_in( evt.particles()[particle_in-1] );
280  }
281  else {
282  shared_ptr<IntAttribute> existing_hc=evt.attribute<IntAttribute>("cycles");
283  if (existing_hc)
284  {
285  if (existing_hc->value()!=0&&particle_in > 0 )
286  {
287  WARNING( "ReaderAscii: event has cycles, this vertex might be repeated." )
288  //if( particle_in > 0 ) data->add_particle_in( evt.particles()[particle_in-1] );
289  }
290  }
291  else
292  return false;
293  }
294 
295  // check for next particle or end of particle list
296  if( !(cursor = strchr(cursor+1,',')) ) {
297  if( !(cursor = strchr(cursor2+1,']')) ) return false;
298  break;
299  }
300  }
301 
302  // check if there is position information
303  if( (cursor = strchr(cursor+1,'@')) ) {
304 
305  // x
306  if( !(cursor = strchr(cursor+1,' ')) ) return false;
307  position.setX(atof(cursor));
308 
309  // y
310  if( !(cursor = strchr(cursor+1,' ')) ) return false;
311  position.setY(atof(cursor));
312 
313  // z
314  if( !(cursor = strchr(cursor+1,' ')) ) return false;
315  position.setZ(atof(cursor));
316 
317  // t
318  if( !(cursor = strchr(cursor+1,' ')) ) return false;
319  position.setT(atof(cursor));
320  data->set_position( position );
321 
322  }
323 
324  DEBUG( 10, "ReaderAscii: V: "<<id<<" with "<<data->particles_in().size()<<" particles)" )
325 
326  evt.add_vertex(data);
327 
328  return true;
329 }
330 
331 
333  GenParticlePtr data = make_shared<GenParticle>();
334  FourVector momentum;
335  const char *cursor = buf;
336  int mother_id = 0;
337 
338  // verify id
339  if( !(cursor = strchr(cursor+1,' ')) ) return false;
340 
341  if( atoi(cursor) != (int)evt.particles().size() + 1 ) {
342  /// @todo Should be an exception
343  ERROR( "ReaderAscii: particle ID mismatch" )
344  return false;
345  }
346 
347  // mother id
348  if( !(cursor = strchr(cursor+1,' ')) ) return false;
349  mother_id = atoi(cursor);
350 
351  // add particle to corresponding vertex
352  if( mother_id > 0 && mother_id <= (int)evt.particles().size() ) {
353 
354  GenParticlePtr mother = evt.particles()[ mother_id-1 ];
355  GenVertexPtr vertex = mother->end_vertex();
356 
357  // create new vertex if needed
358  if( !vertex ) {
359  vertex = make_shared<GenVertex>();
360  vertex->add_particle_in(mother);
361  }
362 
363  vertex->add_particle_out(data);
364  evt.add_vertex(vertex);
365  }
366  else if( mother_id < 0 && -mother_id <= (int)evt.vertices().size() ) {
367  evt.vertices()[ (-mother_id)-1 ]->add_particle_out(data);
368  }
369 
370  // pdg id
371  if( !(cursor = strchr(cursor+1,' ')) ) return false;
372  data->set_pid( atoi(cursor) );
373 
374  // px
375  if( !(cursor = strchr(cursor+1,' ')) ) return false;
376  momentum.setPx(atof(cursor));
377 
378  // py
379  if( !(cursor = strchr(cursor+1,' ')) ) return false;
380  momentum.setPy(atof(cursor));
381 
382  // pz
383  if( !(cursor = strchr(cursor+1,' ')) ) return false;
384  momentum.setPz(atof(cursor));
385 
386  // pe
387  if( !(cursor = strchr(cursor+1,' ')) ) return false;
388  momentum.setE(atof(cursor));
389  data->set_momentum(momentum);
390 
391  // m
392  if( !(cursor = strchr(cursor+1,' ')) ) return false;
393  data->set_generated_mass( atof(cursor) );
394 
395  // status
396  if( !(cursor = strchr(cursor+1,' ')) ) return false;
397  data->set_status( atoi(cursor) );
398 
399  evt.add_particle(data);
400 
401  DEBUG( 10, "ReaderAscii: P: "<<data->id()<<" ( mother: "<<mother_id<<", pid: "<<data->pid()<<")" )
402 
403  return true;
404 }
405 
406 
407 bool ReaderAscii::parse_attribute(GenEvent &evt, const char *buf) {
408  const char *cursor = buf;
409  const char *cursor2 = buf;
410  char name[64];
411  int id = 0;
412 
413  if( !(cursor = strchr(cursor+1,' ')) ) return false;
414  id = atoi(cursor);
415 
416  if( !(cursor = strchr(cursor+1,' ')) ) return false;
417  ++cursor;
418 
419  if( !(cursor2 = strchr(cursor,' ')) ) return false;
420  sprintf(name,"%.*s", (int)(cursor2-cursor), cursor);
421 
422  cursor = cursor2+1;
423 
424  shared_ptr<Attribute> att =
425  make_shared<StringAttribute>( StringAttribute(unescape(cursor)) );
426 
427  evt.add_attribute(string(name), att, id);
428 
429  return true;
430 }
431 
432 bool ReaderAscii::parse_run_attribute(const char *buf) {
433  const char *cursor = buf;
434  const char *cursor2 = buf;
435  char name[64];
436 
437  if( !(cursor = strchr(cursor+1,' ')) ) return false;
438  ++cursor;
439 
440  if( !(cursor2 = strchr(cursor,' ')) ) return false;
441  sprintf(name,"%.*s", (int)(cursor2-cursor), cursor);
442 
443  cursor = cursor2+1;
444 
445  shared_ptr<StringAttribute> att =
446  make_shared<StringAttribute>( StringAttribute(unescape(cursor)) );
447 
448  run_info()->add_attribute(string(name), att);
449 
450  return true;
451 
452 }
453 
454 
455 bool ReaderAscii::parse_weight_names(const char *buf) {
456  const char *cursor = buf;
457 
458  if( !(cursor = strchr(cursor+1,' ')) ) return false;
459  ++cursor;
460 
461  istringstream iss(unescape(cursor));
462  vector<string> names;
463  string name;
464  while ( iss >> name ) names.push_back(name);
465 
466  run_info()->set_weight_names(names);
467 
468  return true;
469 
470 }
471 
472 bool ReaderAscii::parse_tool(const char *buf) {
473  const char *cursor = buf;
474 
475  if( !(cursor = strchr(cursor+1,' ')) ) return false;
476  ++cursor;
477  string line = unescape(cursor);
479  string::size_type pos = line.find("\n");
480  tool.name = line.substr(0, pos);
481  line = line.substr(pos + 1);
482  pos = line.find("\n");
483  tool.version = line.substr(0, pos);
484  tool.description = line.substr(pos + 1);
485  run_info()->tools().push_back(tool);
486 
487  return true;
488 
489 }
490 
491 
492 string ReaderAscii::unescape(const string& s) {
493  string ret;
494  ret.reserve(s.length());
495  for ( string::const_iterator it = s.begin(); it != s.end(); ++it ) {
496  if ( *it == '\\' ) {
497  ++it;
498  if ( *it == '|' )
499  ret += '\n';
500  else
501  ret += *it;
502  } else
503  ret += *it;
504  }
505 
506  return ret;
507 }
508 
509 
511  if( !m_file.is_open()) return;
512  m_file.close();
513 }
514 
515 
516 } // namespace HepMC3
void set_run_info(shared_ptr< GenRunInfo > run)
Set the GenRunInfo object by smart pointer.
Definition: GenEvent.h:129
const std::vector< ConstGenVertexPtr > & vertices() const
Get list of vertices (const)
Definition: GenEvent.cc:44
bool parse_run_attribute(const char *buf)
Parse run-level attribute.
Definition: ReaderAscii.cc:432
ReaderAscii(const std::string &filename)
Constructor.
Definition: ReaderAscii.cc:22
bool parse_units(GenEvent &evt, const char *buf)
Parse units.
Definition: ReaderAscii.cc:232
void add_vertex(GenVertexPtr v)
Add vertex.
Definition: GenEvent.cc:99
static LengthUnit length_unit(const std::string &name)
Get length unit based on its name.
Definition: Units.h:46
bool parse_tool(const char *buf)
Parse run-level tool information.
Definition: ReaderAscii.cc:472
#define ERROR(MESSAGE)
Macro for printing error messages.
Definition: Errors.h:23
Definition of class GenParticle.
std::pair< int, int > parse_event_information(GenEvent &evt, const char *buf)
Parse event.
Definition: ReaderAscii.cc:168
std::ifstream m_file
Input file.
Definition: ReaderAscii.h:149
Definition of class GenVertex.
const Units::LengthUnit & length_unit() const
Get length unit.
Definition: GenEvent.h:143
Attribute that holds a string.
Definition: Attribute.h:336
shared_ptr< T > attribute(const string &name, const int &id=0) const
Get attribute of type T.
Definition: GenEvent.h:381
void add_particle(GenParticlePtr p)
Add particle.
Definition: GenEvent.cc:50
static std::string name(MomentumUnit u)
Get name of momentum unit.
Definition: Units.h:56
bool failed()
Return status of the stream.
Definition: ReaderAscii.h:45
string name
The name of the tool.
Definition: GenRunInfo.h:40
bool parse_weight_names(const char *buf)
Parse run-level weight names.
Definition: ReaderAscii.cc:455
bool parse_attribute(GenEvent &evt, const char *buf)
Parse attribute.
Definition: ReaderAscii.cc:407
#define DEBUG(LEVEL, MESSAGE)
Macro for printing debug messages with appropriate debug level.
Definition: Errors.h:32
shared_ptr< GenRunInfo > run_info() const
Get the global GenRunInfo object.
Definition: Reader.h:37
LengthUnit
Position units.
Definition: Units.h:32
void setY(double yy)
Set y-component of position/displacement.
Definition: FourVector.h:68
MomentumUnit
Momentum units.
Definition: Units.h:29
const Units::MomentumUnit & momentum_unit() const
Get momentum unit.
Definition: GenEvent.h:141
void setZ(double zz)
Set z-component of position/displacement.
Definition: FourVector.h:73
Stores event-related information.
Definition: GenEvent.h:42
Generic 4-vector.
Definition: FourVector.h:34
string version
The version of the tool.
Definition: GenRunInfo.h:43
void setT(double tt)
Set time component of position/displacement.
Definition: FourVector.h:78
Interrnal struct for keeping track of tools.
Definition: GenRunInfo.h:37
bool read_event(GenEvent &evt)
Load event from file.
Definition: ReaderAscii.cc:47
Definition of class ReaderAscii.
Definition of class Units.
#define WARNING(MESSAGE)
Macro for printing warning messages.
Definition: Errors.h:26
static MomentumUnit momentum_unit(const std::string &name)
Get momentum unit based on its name.
Definition: Units.h:36
void set_units(Units::MomentumUnit new_momentum_unit, Units::LengthUnit new_length_unit)
Change event units Converts event from current units to new ones.
Definition: GenEvent.cc:368
void add_attribute(const string &name, const shared_ptr< Attribute > &att, const int &id=0)
Add event attribute to event.
Definition: GenEvent.h:202
void set_run_info(shared_ptr< GenRunInfo > run)
Set the global GenRunInfo object.
Definition: Reader.h:44
bool parse_vertex_information(GenEvent &evt, const char *buf)
Parse vertex.
Definition: ReaderAscii.cc:253
const std::vector< double > & weights() const
Get event weight values as a vector.
Definition: GenEvent.h:87
string description
Other information about how the tool was used in the run.
Definition: GenRunInfo.h:47
std::string unescape(const std::string &s)
Unsecape &#39;\&#39; and &#39; &#39; characters in string.
Definition: ReaderAscii.cc:492
Definition of class GenEvent.
bool parse_particle_information(GenEvent &evt, const char *buf)
Parse particle.
Definition: ReaderAscii.cc:332
void close()
Close file stream.
Definition: ReaderAscii.cc:510
void setX(double xx)
Set x-component of position/displacement.
Definition: FourVector.h:63
~ReaderAscii()
Destructor.
Definition: ReaderAscii.cc:44
void clear()
Remove contents of this event.
Definition: GenEvent.cc:411
const std::vector< ConstGenParticlePtr > & particles() const
Get list of particles (const)
Definition: GenEvent.cc:40
Attribute that holds an Integer implemented as an int.
Definition: Attribute.h:158
bool parse_weight_values(GenEvent &evt, const char *buf)
Parse weight value lines.
Definition: ReaderAscii.cc:215
void set_event_number(const int &num)
Set event number.
Definition: GenEvent.h:138
void shift_position_to(const FourVector &newpos)
Shift position of all vertices in the event to op.
Definition: GenEvent.h:188