HepMC event record
ReaderAsciiHepMC2.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // This file is part of HepMC
4 // Copyright (C) 2014-2015 The HepMC collaboration (see AUTHORS for details)
5 //
6 /**
7  * @file ReaderAsciiHepMC2.cc
8  * @brief Implementation of \b class ReaderAsciiHepMC2
9  *
10  */
11 #include "HepMC/ReaderAsciiHepMC2.h"
12 
13 #include "HepMC/GenEvent.h"
14 #include "HepMC/GenVertex.h"
15 #include "HepMC/GenParticle.h"
16 #include "HepMC/GenHeavyIon.h"
17 #include "HepMC/GenPdfInfo.h"
18 #include "HepMC/Setup.h"
19 
20 #include <cstring>
21 #include <cstdlib>
22 
23 namespace HepMC {
24 
25 ReaderAsciiHepMC2::ReaderAsciiHepMC2(const std::string& filename):
26 m_file(filename) {
27  if( !m_file.is_open() ) {
28  ERROR( "ReaderAsciiHepMC2: could not open input file: "<<filename )
29  }
30  set_run_info(make_shared<GenRunInfo>());
31 }
32 
34  char buf[512];
35  bool parsed_event_header = false;
36  bool is_parsing_successful = true;
37  int parsing_result = 0;
38  unsigned int vertices_count = 0;
39  unsigned int current_vertex_particles_count = 0;
40  unsigned int current_vertex_particles_parsed= 0;
41 
42  evt.clear();
43  evt.set_run_info(run_info());
44 
45  // Empty cache
46  m_vertex_cache.clear();
47  m_vertex_barcodes.clear();
48 
49  m_particle_cache.clear();
50  m_end_vertex_barcodes.clear();
51 
52  //
53  // Parse event, vertex and particle information
54  //
55  while(!failed()) {
56  m_file.getline(buf,512);
57 
58  if( strlen(buf) == 0 ) continue;
59 
60  // Check for IO_GenEvent header/footer
61  if( strncmp(buf,"HepMC",5) == 0 ) {
62  if(parsed_event_header) {
63  is_parsing_successful = true;
64  break;
65  }
66  continue;
67  }
68 
69  switch(buf[0]) {
70  case 'E':
71  parsing_result = parse_event_information(evt,buf);
72  if(parsing_result<0) {
73  is_parsing_successful = false;
74  ERROR( "ReaderAsciiHepMC2: error parsing event information" )
75  }
76  else {
77  vertices_count = parsing_result;
78  m_vertex_cache.reserve(vertices_count);
79  m_particle_cache.reserve(vertices_count*3);
80  m_vertex_barcodes.reserve(vertices_count);
81  m_end_vertex_barcodes.reserve(vertices_count*3);
82  is_parsing_successful = true;
83  }
84  parsed_event_header = true;
85  break;
86  case 'V':
87  // If starting new vertex: verify if previous was fully parsed
88 
89  /** @bug HepMC2 files produced with Pythia8 are known to have wrong
90  information about number of particles in vertex. Hence '<' sign */
91  if(current_vertex_particles_parsed < current_vertex_particles_count) {
92  is_parsing_successful = false;
93  break;
94  }
95  current_vertex_particles_parsed = 0;
96 
97  parsing_result = parse_vertex_information(buf);
98 
99  if(parsing_result<0) {
100  is_parsing_successful = false;
101  ERROR( "ReaderAsciiHepMC2: error parsing vertex information" )
102  }
103  else {
104  current_vertex_particles_count = parsing_result;
105  is_parsing_successful = true;
106  }
107  break;
108  case 'P':
109 
110  parsing_result = parse_particle_information(buf);
111 
112  if(parsing_result<0) {
113  is_parsing_successful = false;
114  ERROR( "ReaderAsciiHepMC2: error parsing particle information" )
115  }
116  else {
117  ++current_vertex_particles_parsed;
118  is_parsing_successful = true;
119  }
120  break;
121  case 'U':
122  is_parsing_successful = parse_units(evt,buf);
123  break;
124  case 'F':
125  is_parsing_successful = parse_pdf_info(evt,buf);
126  break;
127  case 'H':
128  is_parsing_successful = parse_heavy_ion(evt,buf);
129  break;
130  case 'N':
131  is_parsing_successful = parse_weight_names(buf);
132  break;
133  case 'C':
134  is_parsing_successful = parse_xs_info(evt,buf);
135  break;
136  default:
137  WARNING( "ReaderAsciiHepMC2: skipping unrecognised prefix: " << buf[0] )
138  is_parsing_successful = true;
139  break;
140  }
141 
142  if( !is_parsing_successful ) break;
143 
144  // Check for next event
145  buf[0] = m_file.peek();
146  if( parsed_event_header && buf[0]=='E' ) break;
147  }
148 
149  // Check if all particles in last vertex were parsed
150 
151  /** @bug HepMC2 files produced with Pythia8 are known to have wrong
152  information about number of particles in vertex. Hence '<' sign */
153  if( is_parsing_successful && current_vertex_particles_parsed < current_vertex_particles_count ) {
154  ERROR( "ReaderAsciiHepMC2: not all particles parsed" )
155  is_parsing_successful = false;
156  }
157  // Check if all vertices were parsed
158  else if( is_parsing_successful && m_vertex_cache.size() != vertices_count ) {
159  ERROR( "ReaderAsciiHepMC2: not all vertices parsed" )
160  is_parsing_successful = false;
161  }
162 
163  if( !is_parsing_successful ) {
164  ERROR( "ReaderAsciiHepMC2: event parsing failed. Returning empty event" )
165  DEBUG( 1, "Parsing failed at line:" << std::endl << buf )
166 
167  evt.clear();
168  m_file.clear(std::ios::badbit);
169 
170  return 0;
171  }
172 
173  // Restore production vertex pointers
174  for(unsigned int i=0; i<m_particle_cache.size(); ++i) {
175  if( !m_end_vertex_barcodes[i] ) continue;
176 
177  for(unsigned int j=0; j<m_vertex_cache.size(); ++j) {
179  m_vertex_cache[j]->add_particle_in(m_particle_cache[i]);
180  break;
181  }
182  }
183  }
184 
185  // Remove vertices with no incoming particles or no outgoing particles
186  for(unsigned int i=0; i<m_vertex_cache.size(); ++i) {
187  if( m_vertex_cache[i]->particles_in().size() == 0 ) {
188  m_vertex_cache[i] = NULL;
189  }
190  else if( m_vertex_cache[i]->particles_out().size() == 0 ) {
191  m_vertex_cache[i] = NULL;
192  }
193  }
194 
195  // Reserve memory for the event
196  evt.reserve( m_particle_cache.size(), m_vertex_cache.size() );
197 
198  // Add whole event tree in topological order
199  evt.add_tree( m_particle_cache );
200 
201  return 1;
202 }
203 
205  const char *cursor = buf;
206  int event_no = 0;
207  int vertices_count = 0;
208  int random_states_size = 0;
209  int weights_size = 0;
210  std::vector<long> random_states(0);
211  std::vector<double> weights(0);
212 
213  // event number
214  if( !(cursor = strchr(cursor+1,' ')) ) return -1;
215  event_no = atoi(cursor);
216  evt.set_event_number(event_no);
217 
218  // SKIPPED: mpi
219  if( !(cursor = strchr(cursor+1,' ')) ) return -1;
220 
221  // SKIPPED: event scale
222  if( !(cursor = strchr(cursor+1,' ')) ) return -1;
223 
224  // SKIPPED: alpha_qcd
225  if( !(cursor = strchr(cursor+1,' ')) ) return -1;
226 
227  // SKIPPED: alpha_qed
228  if( !(cursor = strchr(cursor+1,' ')) ) return -1;
229 
230  // SKIPPED: signal_process_id
231  if( !(cursor = strchr(cursor+1,' ')) ) return -1;
232 
233  // SKIPPED: signal_process_vertex
234  if( !(cursor = strchr(cursor+1,' ')) ) return -1;
235 
236  // num_vertices
237  if( !(cursor = strchr(cursor+1,' ')) ) return -1;
238  vertices_count = atoi(cursor);
239 
240  // SKIPPED: beam 1
241  if( !(cursor = strchr(cursor+1,' ')) ) return -1;
242 
243  // SKIPPED: beam 2
244  if( !(cursor = strchr(cursor+1,' ')) ) return -1;
245 
246  // SKIPPED: random states
247  if( !(cursor = strchr(cursor+1,' ')) ) return -1;
248  random_states_size = atoi(cursor);
249  random_states.resize(random_states_size);
250 
251  for ( int i = 0; i < random_states_size; ++i ) {
252  if( !(cursor = strchr(cursor+1,' ')) ) return -1;
253  random_states[i] = atoi(cursor);
254  }
255 
256  // weights
257  if( !(cursor = strchr(cursor+1,' ')) ) return -1;
258  weights_size = atoi(cursor);
259  weights.resize(weights_size);
260 
261  for ( int i = 0; i < weights_size; ++i ) {
262  if( !(cursor = strchr(cursor+1,' ')) ) return -1;
263  weights[i] = atof(cursor);
264  }
265 
266  evt.weights() = weights;
267 
268  DEBUG( 10, "ReaderAsciiHepMC2: E: "<<event_no<<" ("<<vertices_count<<"V, "<<weights_size<<"W, "<<random_states_size<<"RS)" )
269 
270  return vertices_count;
271 }
272 
273 bool ReaderAsciiHepMC2::parse_units(GenEvent &evt, const char *buf) {
274  const char *cursor = buf;
275 
276  // momentum
277  if( !(cursor = strchr(cursor+1,' ')) ) return false;
278  ++cursor;
279  Units::MomentumUnit momentum_unit = Units::momentum_unit(cursor);
280 
281  // length
282  if( !(cursor = strchr(cursor+1,' ')) ) return false;
283  ++cursor;
284  Units::LengthUnit length_unit = Units::length_unit(cursor);
285 
286  evt.set_units(momentum_unit,length_unit);
287 
288  DEBUG( 10, "ReaderAsciiHepMC2: U: " << Units::name(evt.momentum_unit()) << " " << Units::name(evt.length_unit()) )
289 
290  return true;
291 }
292 
294  GenVertexPtr data = make_shared<GenVertex>();
295  FourVector position;
296  const char *cursor = buf;
297  int barcode = 0;
298  int num_particles_out = 0;
299 
300  // barcode
301  if( !(cursor = strchr(cursor+1,' ')) ) return -1;
302  barcode = atoi(cursor);
303 
304  // SKIPPED: id
305  if( !(cursor = strchr(cursor+1,' ')) ) return -1;
306 
307  // x
308  if( !(cursor = strchr(cursor+1,' ')) ) return -1;
309  position.setX(atof(cursor));
310 
311  // y
312  if( !(cursor = strchr(cursor+1,' ')) ) return -1;
313  position.setY(atof(cursor));
314 
315  // z
316  if( !(cursor = strchr(cursor+1,' ')) ) return -1;
317  position.setZ(atof(cursor));
318 
319  // t
320  if( !(cursor = strchr(cursor+1,' ')) ) return -1;
321  position.setT(atof(cursor));
322  data->set_position( position );
323 
324  // SKIPPED: num_orphans_in
325  if( !(cursor = strchr(cursor+1,' ')) ) return -1;
326 
327  // num_particles_out
328  if( !(cursor = strchr(cursor+1,' ')) ) return -1;
329  num_particles_out = atoi(cursor);
330 
331  // SKIPPING: weights_size, weights
332 
333  // Add original vertex barcode to the cache
334  m_vertex_cache.push_back( data );
335  m_vertex_barcodes.push_back( barcode );
336 
337  DEBUG( 10, "ReaderAsciiHepMC2: V: "<<-(int)m_vertex_cache.size()<<" (old barcode"<<barcode<<") "<<num_particles_out<<" particles)" )
338 
339  return num_particles_out;
340 }
341 
343  GenParticlePtr data = make_shared<GenParticle>();
344  FourVector momentum;
345  const char *cursor = buf;
346  int end_vtx = 0;
347 
348  /// @todo barcode ignored but maybe should be put as an attribute?...
349  if( !(cursor = strchr(cursor+1,' ')) ) return -1;
350 
351  // id
352  if( !(cursor = strchr(cursor+1,' ')) ) return -1;
353  data->set_pid( atoi(cursor) );
354 
355  // px
356  if( !(cursor = strchr(cursor+1,' ')) ) return -1;
357  momentum.setPx(atof(cursor));
358 
359  // py
360  if( !(cursor = strchr(cursor+1,' ')) ) return -1;
361  momentum.setPy(atof(cursor));
362 
363  // pz
364  if( !(cursor = strchr(cursor+1,' ')) ) return -1;
365  momentum.setPz(atof(cursor));
366 
367  // pe
368  if( !(cursor = strchr(cursor+1,' ')) ) return -1;
369  momentum.setE(atof(cursor));
370  data->set_momentum(momentum);
371 
372  // m
373  if( !(cursor = strchr(cursor+1,' ')) ) return -1;
374  data->set_generated_mass( atof(cursor) );
375 
376  // status
377  if( !(cursor = strchr(cursor+1,' ')) ) return -1;
378  data->set_status( atoi(cursor) );
379 
380  // SKIPPED: theta
381  if( !(cursor = strchr(cursor+1,' ')) ) return -1;
382 
383  // SKIPPED: phi
384  if( !(cursor = strchr(cursor+1,' ')) ) return -1;
385 
386  // end_vtx_code
387  if( !(cursor = strchr(cursor+1,' ')) ) return -1;
388  end_vtx = atoi(cursor);
389 
390  // SKIPPING: flow_size, flow patterns
391 
392  // Set prod_vtx link
393  if( end_vtx == m_vertex_barcodes.back() ) {
394  m_vertex_cache.back()->add_particle_in(data);
395  end_vtx = 0;
396  }
397  else {
398  m_vertex_cache.back()->add_particle_out(data);
399  }
400 
401  m_particle_cache.push_back( data );
402  m_end_vertex_barcodes.push_back( end_vtx );
403 
404  DEBUG( 10, "ReaderAsciiHepMC2: P: "<<m_particle_cache.size()<<" ( pid: "<<data->pid()<<") end vertex: "<<end_vtx )
405 
406  return 0;
407 }
408 
409 bool ReaderAsciiHepMC2::parse_xs_info(GenEvent &evt, const char *buf) {
410  const char *cursor = buf;
411  shared_ptr<GenCrossSection> xs = make_shared<GenCrossSection>();
412 
413  if( !(cursor = strchr(cursor+1,' ')) ) return false;
414  double xs_val = atof(cursor);
415 
416  if( !(cursor = strchr(cursor+1,' ')) ) return false;
417  double xs_err = atof(cursor);
418 
419  xs->set_cross_section( xs_val , xs_err);
420  evt.add_attribute("GenCrossSection",xs);
421 
422  return true;
423 }
424 
426  const char *cursor = buf;
427  const char *cursor2 = buf;
428  int w_count = 0;
429  vector<string> w_names;
430 
431  // Ignore weight names if no GenRunInfo object
432  if( !run_info() ) return true;
433 
434  if( !(cursor = strchr(cursor+1,' ')) ) return false;
435  w_count = atoi(cursor);
436 
437  if( w_count <= 0 ) return false;
438 
439  w_names.resize(w_count);
440 
441  for( int i=0; i < w_count; ++i ) {
442  // Find pair of '"' characters
443  if( !(cursor = strchr(cursor+1,'"')) ) return false;
444  if( !(cursor2 = strchr(cursor+1,'"')) ) return false;
445 
446  // Strip cursor of leading '"' character
447  ++cursor;
448 
449  w_names[i].assign(cursor, cursor2-cursor);
450 
451  cursor = cursor2;
452  }
453 
454  run_info()->set_weight_names(w_names);
455 
456  return true;
457 }
458 
459 bool ReaderAsciiHepMC2::parse_heavy_ion(GenEvent &evt, const char *buf) {
460  shared_ptr<GenHeavyIon> hi = make_shared<GenHeavyIon>();
461  const char *cursor = buf;
462 
463  if( !(cursor = strchr(cursor+1,' ')) ) return false;
464  hi->Ncoll_hard = atoi(cursor);
465 
466  if( !(cursor = strchr(cursor+1,' ')) ) return false;
467  hi->Npart_proj = atoi(cursor);
468 
469  if( !(cursor = strchr(cursor+1,' ')) ) return false;
470  hi->Npart_targ = atoi(cursor);
471 
472  if( !(cursor = strchr(cursor+1,' ')) ) return false;
473  hi->Ncoll = atoi(cursor);
474 
475  if( !(cursor = strchr(cursor+1,' ')) ) return false;
476  hi->spectator_neutrons = atoi(cursor);
477 
478  if( !(cursor = strchr(cursor+1,' ')) ) return false;
479  hi->spectator_protons = atoi(cursor);
480 
481  if( !(cursor = strchr(cursor+1,' ')) ) return false;
482  hi->N_Nwounded_collisions = atoi(cursor);
483 
484  if( !(cursor = strchr(cursor+1,' ')) ) return false;
485  hi->Nwounded_N_collisions = atoi(cursor);
486 
487  if( !(cursor = strchr(cursor+1,' ')) ) return false;
488  hi->Nwounded_Nwounded_collisions = atoi(cursor);
489 
490  if( !(cursor = strchr(cursor+1,' ')) ) return false;
491  hi->impact_parameter = atof(cursor);
492 
493  if( !(cursor = strchr(cursor+1,' ')) ) return false;
494  hi->event_plane_angle = atof(cursor);
495 
496  if( !(cursor = strchr(cursor+1,' ')) ) return false;
497  hi->eccentricity = atof(cursor);
498 
499  if( !(cursor = strchr(cursor+1,' ')) ) return false;
500  hi->sigma_inel_NN = atof(cursor);
501 
502  // Not in HepMC2:
503  hi->centrality = 0.0;
504 
505  evt.add_attribute("GenHeavyIon",hi);
506 
507  return true;
508 }
509 
510 bool ReaderAsciiHepMC2::parse_pdf_info(GenEvent &evt, const char *buf) {
511  shared_ptr<GenPdfInfo> pi = make_shared<GenPdfInfo>();
512  const char *cursor = buf;
513 
514  if( !(cursor = strchr(cursor+1,' ')) ) return false;
515  pi->parton_id[0] = atoi(cursor);
516 
517  if( !(cursor = strchr(cursor+1,' ')) ) return false;
518  pi->parton_id[1] = atoi(cursor);
519 
520  if( !(cursor = strchr(cursor+1,' ')) ) return false;
521  pi->x[0] = atof(cursor);
522 
523  if( !(cursor = strchr(cursor+1,' ')) ) return false;
524  pi->x[1] = atof(cursor);
525 
526  if( !(cursor = strchr(cursor+1,' ')) ) return false;
527  pi->scale = atof(cursor);
528 
529  if( !(cursor = strchr(cursor+1,' ')) ) return false;
530  pi->xf[0] = atof(cursor);
531 
532  if( !(cursor = strchr(cursor+1,' ')) ) return false;
533  pi->xf[1] = atof(cursor);
534 
535  if( !(cursor = strchr(cursor+1,' ')) ) return false;
536  pi->pdf_id[0] = atoi(cursor);
537 
538  if( !(cursor = strchr(cursor+1,' ')) ) return false;
539  pi->pdf_id[1] = atoi(cursor);
540 
541  evt.add_attribute("GenPdfInfo",pi);
542 
543  return true;
544 }
545 
547  if( !m_file.is_open() ) return;
548  m_file.close();
549 }
550 
551 } // namespace HepMC
bool parse_heavy_ion(GenEvent &evt, const char *buf)
Parse heavy ion information.
const Units::MomentumUnit & momentum_unit() const
Get momentum unit.
void set_event_number(int num)
Set event number.
static std::string name(MomentumUnit u)
Get name of momentum unit.
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:335
vector< GenParticlePtr > m_particle_cache
Particle cache.
void add_tree(const std::vector< GenParticlePtr > &particles)
Add whole tree in topological order.
Definition: GenEvent.cc:262
ReaderAsciiHepMC2(const std::string &filename)
Default constructor.
bool read_event(GenEvent &evt)
Implementation of Reader::read_event.
bool parse_weight_names(const char *buf)
Parse weight names.
vector< int > m_vertex_barcodes
Old vertex barcodes.
bool parse_units(GenEvent &evt, const char *buf)
Parse units.
void close()
Close file stream.
vector< GenVertexPtr > m_vertex_cache
Vertex cache.
shared_ptr< GenRunInfo > run_info() const
Get the global GenRunInfo object.
void clear()
Remove contents of this event.
Definition: GenEvent.cc:374
void set_run_info(shared_ptr< GenRunInfo > run)
Set the global GenRunInfo object.
Stores event-related information.
static MomentumUnit momentum_unit(const std::string &name)
Get momentum unit based on its name.
bool parse_xs_info(GenEvent &evt, const char *buf)
Parse pdf information.
int parse_vertex_information(const char *buf)
Parse vertex.
void add_attribute(const string &name, const shared_ptr< Attribute > &att, int id=0)
Add event attribute to event.
int parse_event_information(GenEvent &evt, const char *buf)
Parse event.
void reserve(unsigned int particles, unsigned int vertices=0)
Reserve memory for particles and vertices.
Definition: GenEvent.cc:329
Definition of template class SmartPointer.
bool parse_pdf_info(GenEvent &evt, const char *buf)
Parse pdf information.
vector< int > m_end_vertex_barcodes
Old end vertex barcodes.
const Units::LengthUnit & length_unit() const
Get length unit.
int parse_particle_information(const char *buf)
Parse particle.
static LengthUnit length_unit(const std::string &name)
Get length unit based on its name.
void set_run_info(shared_ptr< GenRunInfo > run)
Set the GenRunInfo object by smart pointer.
const std::vector< double > & weights() const
Get event weight values as a vector.