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