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
19namespace HepMC3 {
20
21
22ReaderAscii::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
33ReaderAscii::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();
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
190pair<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
237bool 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
254bool 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
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
435bool 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
460bool 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
483bool 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
500bool 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
520string 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
537bool 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
#define WARNING(MESSAGE)
Macro for printing warning messages.
Definition Errors.h:26
#define DEBUG(LEVEL, MESSAGE)
Macro for printing debug messages with appropriate debug level.
Definition Errors.h:32
#define ERROR(MESSAGE)
Macro for printing error messages.
Definition Errors.h:23
Definition of class GenEvent.
Definition of class GenParticle.
Definition of class GenVertex.
Definition of class ReaderAscii.
Definition of class Units.
Generic 4-vector.
Definition FourVector.h:35
void setE(double ee)
Definition FourVector.h:116
void setT(double tt)
Definition FourVector.h:87
void setPz(double pzz)
Definition FourVector.h:109
void setY(double yy)
Definition FourVector.h:73
void setPy(double pyy)
Definition FourVector.h:102
void setX(double xx)
Definition FourVector.h:66
void setPx(double pxx)
Definition FourVector.h:95
void setZ(double zz)
Definition FourVector.h:80
Stores event-related information.
Definition GenEvent.h:42
void add_vertex(GenVertexPtr v)
Add vertex.
Definition GenEvent.cc:98
void shift_position_to(const FourVector &newpos)
Shift position of all vertices in the event to op.
Definition GenEvent.h:188
void add_particle(GenParticlePtr p)
Add particle.
Definition GenEvent.cc:50
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_event_number(const int &num)
Set event number.
Definition GenEvent.h:138
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
const Units::MomentumUnit & momentum_unit() const
Get momentum unit.
Definition GenEvent.h:141
const Units::LengthUnit & length_unit() const
Get length unit.
Definition GenEvent.h:143
const std::vector< double > & weights() const
Get event weight values as a vector.
Definition GenEvent.h:87
void clear()
Remove contents of this event.
Definition GenEvent.cc:610
const std::vector< ConstGenParticlePtr > & particles() const
Get list of particles (const)
Definition GenEvent.cc:40
bool parse_tool(const char *buf)
Parse run-level tool information.
bool m_isstream
toggles usage of m_file or m_stream
bool parse_weight_values(GenEvent &evt, const char *buf)
Parse weight value lines.
std::map< GenParticlePtr, int > m_forward_daughters
Temp storage for prod vertex ids.
bool read_event(GenEvent &evt)
Load event from file.
std::ifstream m_file
Input file.
bool parse_units(GenEvent &evt, const char *buf)
Parse units.
bool parse_particle_information(GenEvent &evt, const char *buf)
Parse particle.
void close()
Close file stream.
std::map< GenVertexPtr, std::set< int > > m_forward_mothers
Temp storage for outgoing particle ids.
bool parse_attribute(GenEvent &evt, const char *buf)
Parse attribute.
std::pair< int, int > parse_event_information(GenEvent &evt, const char *buf)
Parse event.
bool parse_vertex_information(GenEvent &evt, const char *buf)
Parse vertex.
std::string unescape(const std::string &s)
Unsecape '\' and ' ' characters in string.
~ReaderAscii()
Destructor.
bool parse_weight_names(const char *buf)
Parse run-level weight names.
bool failed()
Return status of the stream.
ReaderAscii(const std::string &filename)
Constructor.
std::istream * m_stream
For ctor when reading from stdin.
bool parse_run_attribute(const char *buf)
Parse run-level attribute.
void set_run_info(shared_ptr< GenRunInfo > run)
Set the global GenRunInfo object.
Definition Reader.h:46
shared_ptr< GenRunInfo > run_info() const
Get the global GenRunInfo object.
Definition Reader.h:39
Attribute that holds a string.
Definition Attribute.h:336
static LengthUnit length_unit(const std::string &name)
Get length unit based on its name.
Definition Units.h:46
LengthUnit
Position units.
Definition Units.h:32
static std::string name(MomentumUnit u)
Get name of momentum unit.
Definition Units.h:56
MomentumUnit
Momentum units.
Definition Units.h:29
static MomentumUnit momentum_unit(const std::string &name)
Get momentum unit based on its name.
Definition Units.h:36
HepMC3 main namespace.
Definition ReaderGZ.h:28
Interrnal struct for keeping track of tools.
Definition GenRunInfo.h:37
string description
Other information about how the tool was used in the run.
Definition GenRunInfo.h:47
string name
The name of the tool.
Definition GenRunInfo.h:40
string version
The version of the tool.
Definition GenRunInfo.h:43