HepMC3 event record library
ReaderAsciiHepMC2.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // This file is part of HepMC
4 // Copyright (C) 2014-2021 The HepMC collaboration (see AUTHORS for details)
5 //
6 /**
7  * @file ReaderAsciiHepMC2.cc
8  * @brief Implementation of \b class ReaderAsciiHepMC2
9  *
10  */
11 #include <cstring>
12 #include <cstdlib>
13 
15 
16 #include "HepMC3/GenEvent.h"
17 #include "HepMC3/GenVertex.h"
18 #include "HepMC3/GenParticle.h"
19 #include "HepMC3/GenHeavyIon.h"
20 #include "HepMC3/GenPdfInfo.h"
21 #include "HepMC3/Setup.h"
22 
23 namespace HepMC3 {
24 
25 ReaderAsciiHepMC2::ReaderAsciiHepMC2(const std::string& filename):
26  m_file(filename), m_stream(nullptr), m_isstream(false) {
27  if ( !m_file.is_open() ) {
28  HEPMC3_ERROR("ReaderAsciiHepMC2: could not open input file: " << filename )
29  }
30  set_run_info(std::make_shared<GenRunInfo>());
31  m_event_ghost = new GenEvent();
32 }
33 
35  : m_stream(&stream), m_isstream(true)
36 {
37  if ( !m_stream->good() ) {
38  HEPMC3_ERROR("ReaderAsciiHepMC2: could not open input stream ")
39  }
40  set_run_info(std::make_shared<GenRunInfo>());
41  m_event_ghost = new GenEvent();
42 }
43 
44 ReaderAsciiHepMC2::ReaderAsciiHepMC2(std::shared_ptr<std::istream> s_stream)
45  : m_shared_stream(s_stream), m_stream(s_stream.get()), m_isstream(true)
46 {
47  if ( !m_stream->good() ) {
48  HEPMC3_ERROR("ReaderAsciiHepMC2: could not open input stream ")
49  }
50  set_run_info(std::make_shared<GenRunInfo>());
51 }
52 
53 
55 
56 bool ReaderAsciiHepMC2::skip(const int n)
57 {
58  const size_t max_buffer_size = 512*512;
59  char buf[max_buffer_size];
60  int nn = n;
61  while (!failed()) {
62  char peek;
63  if ( (!m_file.is_open()) && (!m_isstream) ) return false;
64  m_isstream ? peek = m_stream->peek() : peek = m_file.peek();
65  if ( peek =='E' ) nn--;
66  if (nn < 0) return true;
67  m_isstream ? m_stream->getline(buf, max_buffer_size) : m_file.getline(buf, max_buffer_size);
68  }
69  return true;
70 }
71 
73  if ( (!m_file.is_open()) && (!m_isstream) ) return false;
74 
75  char peek;
76  const size_t max_buffer_size = 512*512;
77  char buf[max_buffer_size];
78  bool parsed_event_header = false;
79  bool is_parsing_successful = true;
80  int parsing_result = 0;
81  unsigned int vertices_count = 0;
82  unsigned int current_vertex_particles_count = 0;
83  unsigned int current_vertex_particles_parsed = 0;
84 
85  evt.clear();
86  evt.set_run_info(run_info());
87 
88  // Empty cache
89  m_vertex_cache.clear();
90  m_vertex_barcodes.clear();
91 
92  m_particle_cache.clear();
93  m_end_vertex_barcodes.clear();
94  m_particle_cache_ghost.clear();
95  m_vertex_cache_ghost.clear();
96  //
97  // Parse event, vertex and particle information
98  //
99  while (!failed()) {
100  m_isstream ? m_stream->getline(buf, max_buffer_size) : m_file.getline(buf, max_buffer_size);
101  if ( strlen(buf) == 0 ) continue;
102  // Check for IO_GenEvent header/footer
103  if ( strncmp(buf, "HepMC", 5) == 0 ) {
104  if ( strncmp(buf, "HepMC::Version", 14) != 0 && strncmp(buf, "HepMC::IO_GenEvent", 18) != 0 )
105  {
106  HEPMC3_WARNING("ReaderAsciiHepMC2: found unsupported expression in header. Will close the input.")
107  std::cout <<buf << std::endl;
108  m_isstream ? m_stream->clear(std::ios::eofbit) : m_file.clear(std::ios::eofbit);
109  }
110  if (parsed_event_header) {
111  is_parsing_successful = true;
112  break;
113  }
114  continue;
115  }
116  switch (buf[0]) {
117  case 'E':
118  parsing_result = parse_event_information(evt, buf);
119  if (parsing_result < 0) {
120  is_parsing_successful = false;
121  HEPMC3_ERROR("ReaderAsciiHepMC2: HEPMC3_ERROR parsing event information")
122  }
123  else {
124  vertices_count = parsing_result;
125  m_vertex_cache.reserve(vertices_count);
126  m_particle_cache.reserve(vertices_count*3);
127  m_vertex_cache_ghost.reserve(vertices_count);
128  m_particle_cache_ghost.reserve(vertices_count*3);
129  m_vertex_barcodes.reserve(vertices_count);
130  m_end_vertex_barcodes.reserve(vertices_count*3);
131  // Here we make a trick: reserve for this event the vertices_count*3 or the number of particles in the prev. event.
132  evt.reserve(vertices_count, m_particle_cache_ghost.capacity());
133  is_parsing_successful = true;
134  }
135  parsed_event_header = true;
136  break;
137  case 'V':
138  // If starting new vertex: verify if previous was fully parsed
139 
140  /** @bug HepMC2 files produced with Pythia8 are known to have wrong
141  information about number of particles in vertex. Hence '<' sign */
142  if (current_vertex_particles_parsed < current_vertex_particles_count) {
143  is_parsing_successful = false;
144  break;
145  }
146  current_vertex_particles_parsed = 0;
147 
148  parsing_result = parse_vertex_information(buf);
149 
150  if (parsing_result < 0) {
151  is_parsing_successful = false;
152  HEPMC3_ERROR("ReaderAsciiHepMC2: HEPMC3_ERROR parsing vertex information")
153  }
154  else {
155  current_vertex_particles_count = parsing_result;
156  is_parsing_successful = true;
157  }
158  break;
159  case 'P':
160 
161  parsing_result = parse_particle_information(buf);
162 
163  if (parsing_result < 0) {
164  is_parsing_successful = false;
165  HEPMC3_ERROR("ReaderAsciiHepMC2: HEPMC3_ERROR parsing particle information")
166  }
167  else {
168  ++current_vertex_particles_parsed;
169  is_parsing_successful = true;
170  }
171  break;
172  case 'U':
173  is_parsing_successful = parse_units(evt, buf);
174  break;
175  case 'F':
176  is_parsing_successful = parse_pdf_info(evt, buf);
177  break;
178  case 'H':
179  is_parsing_successful = parse_heavy_ion(evt, buf);
180  break;
181  case 'N':
182  is_parsing_successful = parse_weight_names(buf);
183  break;
184  case 'C':
185  is_parsing_successful = parse_xs_info(evt, buf);
186  break;
187  default:
188  HEPMC3_WARNING("ReaderAsciiHepMC2: skipping unrecognised prefix: " << buf[0])
189  is_parsing_successful = true;
190  break;
191  }
192 
193  if ( !is_parsing_successful ) break;
194 
195  // Check for next event
196  m_isstream ? peek = m_stream->peek() : peek = m_file.peek();
197  if ( parsed_event_header && peek == 'E' ) break;
198  }
199 
200  // Check if all particles in last vertex were parsed
201  /** @bug HepMC2 files produced with Pythia8 are known to have wrong
202  information about number of particles in vertex. Hence '<' sign */
203  if (is_parsing_successful && current_vertex_particles_parsed < current_vertex_particles_count) {
204  HEPMC3_ERROR("ReaderAsciiHepMC2: not all particles parsed")
205  is_parsing_successful = false;
206  }
207  // Check if all vertices were parsed
208  else if (is_parsing_successful && m_vertex_cache.size() != vertices_count) {
209  HEPMC3_ERROR("ReaderAsciiHepMC2: not all vertices parsed")
210  is_parsing_successful = false;
211  }
212 
213  if ( !is_parsing_successful ) {
214  HEPMC3_ERROR("ReaderAsciiHepMC2: event parsing failed. Returning empty event")
215  HEPMC3_DEBUG(1, "Parsing failed at line:" << std::endl << buf)
216  evt.clear();
217  m_isstream ? m_stream->clear(std::ios::badbit) : m_file.clear(std::ios::badbit);
218  return 0;
219  }
220 
221  // Restore production vertex pointers
222  for (unsigned int i = 0; i < m_particle_cache.size(); ++i) {
223  if ( !m_end_vertex_barcodes[i] ) continue;
224 
225  for (unsigned int j = 0; j < m_vertex_cache.size(); ++j) {
226  if ( m_vertex_barcodes[j] == m_end_vertex_barcodes[i] ) {
227  m_vertex_cache[j]->add_particle_in(m_particle_cache[i]);
228  break;
229  }
230  }
231  }
232 
233  // Remove vertices with no incoming particles or no outgoing particles
234  for (unsigned int i = 0; i < m_vertex_cache.size(); ++i) {
235  if ( m_vertex_cache[i]->particles_in().size() == 0 ) {
236  HEPMC3_DEBUG(30, "ReaderAsciiHepMC2::read_event - found a vertex without incoming particles: " << m_vertex_cache[i]->id() );
237  //Sometimes the root vertex has no incoming particles. Here we try to save the event.
238  std::vector<GenParticlePtr> beams;
239  for (auto p: m_vertex_cache[i]->particles_out()) if (p->status() == 4 && !(p->end_vertex())) beams.push_back(p);
240  for (auto p: beams)
241  {
242  m_vertex_cache[i]->add_particle_in(p);
243  m_vertex_cache[i]->remove_particle_out(p);
244  HEPMC3_DEBUG(30, "ReaderAsciiHepMC2::read_event - moved particle with status=4 from the outgoing to the incoming particles of vertex: " << m_vertex_cache[i]->id());
245  }
246  if (beams.size() == 0) {
247  HEPMC3_DEBUG(30, "ReaderAsciiHepMC2::read_event - removed vertex without incoming particles: " << m_vertex_cache[i]->id() );
248  m_vertex_cache[i] = nullptr;
249  }
250  }
251  else if ( m_vertex_cache[i]->particles_out().size() == 0 ) {
252  m_vertex_cache[i] = nullptr;
253  HEPMC3_DEBUG(30, "ReaderAsciiHepMC2::read_event - removed vertex without outgoing particles: " << m_vertex_cache[i]->id());
254  }
255  }
256 
257  // Reserve memory for the event
258  evt.reserve(m_particle_cache.size(), m_vertex_cache.size());
259 
260  // Add whole event tree in topological order
262 
263  if (m_options.find("event_random_states_are_separated") != m_options.end())
264  {
265  std::shared_ptr<VectorLongIntAttribute> random_states_a = evt.attribute<VectorLongIntAttribute>("random_states");
266  if (random_states_a) {
267  std::vector<long int> random_states_v = random_states_a->value();
268  for (size_t i = 0; i < random_states_v.size(); ++i )
269  evt.add_attribute("random_states" + std::to_string((long long unsigned int)i), std::make_shared<IntAttribute>(random_states_v[i]));
270  evt.remove_attribute("random_states");
271  }
272 
273  }
274 
275  std::map< std::string, std::map<int, std::shared_ptr<Attribute> > > cached_attributes = m_event_ghost->attributes();
276  if (cached_attributes.find("flows") != cached_attributes.end()) {
277  std::map<int, std::shared_ptr<Attribute> > flows = cached_attributes.at("flows");
278  if (m_options.find("particle_flows_are_separated") == m_options.end()) {
279  for (auto f: flows) if (f.first > 0 && f.first <= (int)m_particle_cache.size()) m_particle_cache[f.first-1]->add_attribute("flows", f.second);
280  } else {
281  for (auto f: flows) if (f.first > 0 && f.first <= (int)m_particle_cache.size()) {
282  std::shared_ptr<VectorIntAttribute> casted = std::dynamic_pointer_cast<VectorIntAttribute>(f.second);
283  if (!casted) continue;//Should not happen
284  std::vector<int> this_p_flow = casted->value();
285  for (size_t i = 0; i<this_p_flow.size(); i++) m_particle_cache[f.first-1]->add_attribute("flow" + std::to_string(i + 1), std::make_shared<IntAttribute>(this_p_flow[i]));
286  }
287  }
288  }
289 
290  if (cached_attributes.find("phi") != cached_attributes.end()) {
291  std::map<int, std::shared_ptr<Attribute> > phi = cached_attributes.at("phi");
292  for (auto f: phi) if (f.first > 0 &&f.first <= (int)m_particle_cache.size()) m_particle_cache[f.first-1]->add_attribute("phi", f.second);
293  }
294 
295  if (cached_attributes.find("theta") != cached_attributes.end()) {
296  std::map<int, std::shared_ptr<Attribute> > theta = cached_attributes.at("theta");
297  for (auto f: theta) if (f.first > 0 && f.first <= (int)m_particle_cache.size()) m_particle_cache[f.first-1]->add_attribute("theta", f.second);
298  }
299 
300  if (cached_attributes.find("weights") != cached_attributes.end()) {
301  std::map<int, std::shared_ptr<Attribute> > weights = cached_attributes.at("weights");
302  if (m_options.find("vertex_weights_are_separated") == m_options.end()) {
303  for (auto f: weights) if (f.first < 0 && f.first >= -(int)m_vertex_cache.size()) m_vertex_cache[-f.first-1]->add_attribute("weights", f.second);
304  } else {
305  for (auto f: weights) if (f.first < 0 && f.first >= -(int)m_vertex_cache.size()) {
306  std::shared_ptr<VectorDoubleAttribute> casted = std::dynamic_pointer_cast<VectorDoubleAttribute>(f.second);
307  if (!casted) continue;//Should not happen
308  std::vector<double> this_v_weight = casted->value();
309  for (size_t i = 0; i < this_v_weight.size(); i++) m_particle_cache[-f.first-1]->add_attribute("weight"+std::to_string(i), std::make_shared<DoubleAttribute>(this_v_weight[i]));
310  }
311  }
312  }
313 
314  std::shared_ptr<IntAttribute> signal_process_vertex_barcode = evt.attribute<IntAttribute>("signal_process_vertex");
315  if (signal_process_vertex_barcode) {
316  int signal_process_vertex_barcode_value = signal_process_vertex_barcode->value();
317  for (unsigned int i = 0; i < m_vertex_cache.size(); ++i)
318  {
319  if (i >= m_vertex_barcodes.size()) continue;//this should not happen!
320  if (signal_process_vertex_barcode_value != m_vertex_barcodes.at(i)) continue;
321  std::shared_ptr<IntAttribute> signal_process_vertex = std::make_shared<IntAttribute>(m_vertex_cache.at(i)->id());
322  evt.add_attribute("signal_process_vertex", signal_process_vertex);
323  break;
324  }
325  }
326  m_particle_cache_ghost.clear();
327  m_vertex_cache_ghost.clear();
328  m_event_ghost->clear();
329  return 1;
330 }
331 
333  const char *cursor = buf;
334  int vertices_count = 0;
335  int random_states_size = 0;
336  int weights_size = 0;
337  std::vector<long> random_states(0);
338  std::vector<double> weights(0);
339 
340  // event number
341  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
342  evt.set_event_number(atoi(cursor));
343 
344  //mpi
345  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
346  evt.add_attribute("mpi", std::make_shared<IntAttribute>(atoi(cursor)));
347 
348  //event scale
349  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
350  evt.add_attribute("event_scale", std::make_shared<DoubleAttribute>(atof(cursor)));
351 
352  //alpha_qcd
353  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
354  evt.add_attribute("alphaQCD", std::make_shared<DoubleAttribute>(atof(cursor)));
355 
356  //alpha_qed
357  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
358  evt.add_attribute("alphaQED", std::make_shared<DoubleAttribute>(atof(cursor)));
359 
360  //signal_process_id
361  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
362  evt.add_attribute("signal_process_id", std::make_shared<IntAttribute>(atoi(cursor)));
363 
364  //signal_process_vertex
365  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
366  evt.add_attribute("signal_process_vertex", std::make_shared<IntAttribute>(atoi(cursor)));
367 
368  // num_vertices
369  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
370  vertices_count = atoi(cursor);
371 
372  // SKIPPED: beam 1
373  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
374 
375  // SKIPPED: beam 2
376  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
377 
378  //random states
379  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
380  random_states_size = atoi(cursor);
381  random_states.resize(random_states_size);
382 
383  for ( int i = 0; i < random_states_size; ++i ) {
384  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
385  random_states[i] = atoi(cursor);
386  }
387 
388  if (random_states.size()) evt.add_attribute("random_states", std::make_shared<VectorLongIntAttribute>(random_states));
389 
390  // weights
391  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
392  weights_size = atoi(cursor);
393  weights.resize(weights_size);
394 
395  for ( int i = 0; i < weights_size; ++i ) {
396  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
397  weights[i] = atof(cursor);
398  }
399 
400  evt.weights() = weights;
401 
402  HEPMC3_DEBUG(10, "ReaderAsciiHepMC2: E: " << evt.event_number() << " (" << vertices_count << "V, " << weights_size << "W, " << random_states_size << "RS)")
403 
404  return vertices_count;
405 }
406 
407 bool ReaderAsciiHepMC2::parse_units(GenEvent &evt, const char *buf) {
408  const char *cursor = buf;
409 
410  // momentum
411  if ( !(cursor = strchr(cursor+1, ' ')) ) return false;
412  ++cursor;
413  Units::MomentumUnit momentum_unit = Units::momentum_unit(cursor);
414 
415  // length
416  if ( !(cursor = strchr(cursor+1, ' ')) ) return false;
417  ++cursor;
418  Units::LengthUnit length_unit = Units::length_unit(cursor);
419 
420  evt.set_units(momentum_unit, length_unit);
421 
422  HEPMC3_DEBUG(10, "ReaderAsciiHepMC2: U: " << Units::name(evt.momentum_unit()) << " " << Units::name(evt.length_unit()))
423 
424  return true;
425 }
426 
428  GenVertexPtr data = std::make_shared<GenVertex>();
429  GenVertexPtr data_ghost = std::make_shared<GenVertex>();
430  const char *cursor = buf;
431  int barcode = 0;
432  int num_particles_out = 0;
433  int weights_size = 0;
434  std::vector<double> weights(0);
435  // barcode
436  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
437  barcode = atoi(cursor);
438 
439  // status
440  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
441  data->set_status(atoi(cursor));
442 
443  // x
444  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
445  double X(atof(cursor));
446 
447  // y
448  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
449  double Y(atof(cursor));
450 
451  // z
452  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
453  double Z(atof(cursor));
454 
455  // t
456  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
457  double T(atof(cursor));
458  data->set_position(FourVector(X,Y,Z,T));
459 
460  // SKIPPED: num_orphans_in
461  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
462 
463  // num_particles_out
464  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
465  num_particles_out = atoi(cursor);
466 
467  // weights
468  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
469  weights_size = atoi(cursor);
470  weights.resize(weights_size);
471 
472  for ( int i = 0; i < weights_size; ++i ) {
473  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
474  weights[i] = atof(cursor);
475  }
476 
477 
478 
479  // Add original vertex barcode to the cache
480  m_vertex_cache.push_back(data);
481  m_vertex_barcodes.push_back(barcode);
482 
483  m_event_ghost->add_vertex(data_ghost);
484 
485  if (weights.size()) data_ghost->add_attribute("weights", std::make_shared<VectorDoubleAttribute>(weights));
486 
487  m_vertex_cache_ghost.push_back(data_ghost);
488 
489  HEPMC3_DEBUG(10, "ReaderAsciiHepMC2: V: " << -(int)m_vertex_cache.size() << " (old barcode " << barcode << ") " << num_particles_out << " particles)")
490 
491  return num_particles_out;
492 }
493 
495  GenParticlePtr data = std::make_shared<GenParticle>();
496  GenParticlePtr data_ghost = std::make_shared<GenParticle>();
497  m_event_ghost->add_particle(data_ghost);
498  const char *cursor = buf;
499  int end_vtx = 0;
500 
501  /// @note barcode is ignored
502  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
503 
504  // id
505  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
506  data->set_pid(atoi(cursor));
507 
508  // px
509  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
510  double Px(atof(cursor));
511 
512  // py
513  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
514  double Py(atof(cursor));
515 
516  // pz
517  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
518  double Pz(atof(cursor));
519 
520  // pe
521  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
522  double E(atof(cursor));
523  data->set_momentum(FourVector(Px,Py,Pz,E));
524 
525  // m
526  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
527  data->set_generated_mass(atof(cursor));
528 
529  // status
530  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
531  data->set_status(atoi(cursor));
532 
533  //theta
534  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
535  double theta_v = atof(cursor);
536  if (theta_v != 0.0) data_ghost->add_attribute("theta", std::make_shared<DoubleAttribute>(theta_v));
537 
538  //phi
539  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
540  double phi_v = atof(cursor);
541  if (phi_v != 0.0) data_ghost->add_attribute("phi", std::make_shared<DoubleAttribute>(phi_v));
542 
543  // end_vtx_code
544  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
545  end_vtx = atoi(cursor);
546 
547  //flow
548  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
549  int flowsize = atoi(cursor);
550 
551  std::map<int, int> flows;
552  for (int i = 0; i < flowsize; i++)
553  {
554  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
555  int flowindex = atoi(cursor);
556  if ( !(cursor = strchr(cursor+1, ' ')) ) return -1;
557  int flowvalue = atoi(cursor);
558  flows[flowindex] = flowvalue;
559  }
560  if (flowsize)
561  {
562  std::vector<int> vectorflows;
563  for (auto f: flows) vectorflows.push_back(f.second);
564  data_ghost->add_attribute("flows", std::make_shared<VectorIntAttribute>(vectorflows));
565  }
566  // Set prod_vtx link
567  if ( end_vtx == m_vertex_barcodes.back() ) {
568  m_vertex_cache.back()->add_particle_in(data);
569  end_vtx = 0;
570  }
571  else {
572  m_vertex_cache.back()->add_particle_out(data);
573  }
574 
575  m_particle_cache.push_back(data);
576  m_particle_cache_ghost.push_back(data_ghost);
577  m_end_vertex_barcodes.push_back(end_vtx);
578 
579  HEPMC3_DEBUG(10, "ReaderAsciiHepMC2: P: " << m_particle_cache.size() << " ( pid: " << data->pid() << ") end vertex: " << end_vtx)
580 
581  return 0;
582 }
583 
584 bool ReaderAsciiHepMC2::parse_xs_info(GenEvent &evt, const char *buf) {
585  const char *cursor = buf;
586  std::shared_ptr<GenCrossSection> xs = std::make_shared<GenCrossSection>();
587 
588  if ( !(cursor = strchr(cursor+1, ' ')) ) return false;
589  double xs_val = atof(cursor);
590 
591  if ( !(cursor = strchr(cursor+1, ' ')) ) return false;
592  double xs_err = atof(cursor);
593 
594  xs->set_cross_section(xs_val, xs_err);
595  evt.add_attribute("GenCrossSection", xs);
596 
597  return true;
598 }
599 
601  const char *cursor = buf;
602  const char *cursor2 = buf;
603  int w_count = 0;
604  std::vector<std::string> w_names;
605 
606  // Ignore weight names if no GenRunInfo object
607  if ( !run_info() ) return true;
608 
609  if ( !(cursor = strchr(cursor+1, ' ')) ) return false;
610  w_count = atoi(cursor);
611 
612  if ( w_count <= 0 ) return false;
613 
614  w_names.resize(w_count);
615 
616  for ( int i=0; i < w_count; ++i ) {
617  // Find pair of '"' characters
618  if ( !(cursor = strchr(cursor+1, '"')) ) return false;
619  if ( !(cursor2 = strchr(cursor+1, '"')) ) return false;
620 
621  // Strip cursor of leading '"' character
622  ++cursor;
623 
624  w_names[i].assign(cursor, cursor2-cursor);
625 
626  cursor = cursor2;
627  }
628 
629  run_info()->set_weight_names(w_names);
630 
631  return true;
632 }
633 
634 bool ReaderAsciiHepMC2::parse_heavy_ion(GenEvent &evt, const char *buf) {
635  std::shared_ptr<GenHeavyIon> hi = std::make_shared<GenHeavyIon>();
636  const char *cursor = buf;
637 
638  if ( !(cursor = strchr(cursor+1, ' ')) ) return false;
639  hi->Ncoll_hard = atoi(cursor);
640 
641  if ( !(cursor = strchr(cursor+1, ' ')) ) return false;
642  hi->Npart_proj = atoi(cursor);
643 
644  if ( !(cursor = strchr(cursor+1, ' ')) ) return false;
645  hi->Npart_targ = atoi(cursor);
646 
647  if ( !(cursor = strchr(cursor+1, ' ')) ) return false;
648  hi->Ncoll = atoi(cursor);
649 
650  if ( !(cursor = strchr(cursor+1, ' ')) ) return false;
651  hi->spectator_neutrons = atoi(cursor);
652 
653  if ( !(cursor = strchr(cursor+1, ' ')) ) return false;
654  hi->spectator_protons = atoi(cursor);
655 
656  if ( !(cursor = strchr(cursor+1, ' ')) ) return false;
657  hi->N_Nwounded_collisions = atoi(cursor);
658 
659  if ( !(cursor = strchr(cursor+1, ' ')) ) return false;
660  hi->Nwounded_N_collisions = atoi(cursor);
661 
662  if ( !(cursor = strchr(cursor+1, ' ')) ) return false;
663  hi->Nwounded_Nwounded_collisions = atoi(cursor);
664 
665  if ( !(cursor = strchr(cursor+1, ' ')) ) return false;
666  hi->impact_parameter = atof(cursor);
667 
668  if ( !(cursor = strchr(cursor+1, ' ')) ) return false;
669  hi->event_plane_angle = atof(cursor);
670 
671  if ( !(cursor = strchr(cursor+1, ' ')) ) return false;
672  hi->eccentricity = atof(cursor);
673 
674  if ( !(cursor = strchr(cursor+1, ' ')) ) return false;
675  hi->sigma_inel_NN = atof(cursor);
676 
677  // Not in HepMC2:
678  hi->centrality = 0.0;
679 
680  evt.add_attribute("GenHeavyIon", hi);
681 
682  return true;
683 }
684 
685 bool ReaderAsciiHepMC2::parse_pdf_info(GenEvent &evt, const char *buf) {
686  std::shared_ptr<GenPdfInfo> pi = std::make_shared<GenPdfInfo>();
687  const char *cursor = buf;
688 
689  if ( !(cursor = strchr(cursor+1, ' ')) ) return false;
690  pi->parton_id[0] = atoi(cursor);
691 
692  if ( !(cursor = strchr(cursor+1, ' ')) ) return false;
693  pi->parton_id[1] = atoi(cursor);
694 
695  if ( !(cursor = strchr(cursor+1, ' ')) ) return false;
696  pi->x[0] = atof(cursor);
697 
698  if ( !(cursor = strchr(cursor+1, ' ')) ) return false;
699  pi->x[1] = atof(cursor);
700 
701  if ( !(cursor = strchr(cursor+1, ' ')) ) return false;
702  pi->scale = atof(cursor);
703 
704  if ( !(cursor = strchr(cursor+1, ' ')) ) return false;
705  pi->xf[0] = atof(cursor);
706 
707  if ( !(cursor = strchr(cursor+1, ' ')) ) return false;
708  pi->xf[1] = atof(cursor);
709 
710  //For compatibility with original HepMC2
711  bool pdfids = true;
712  if ( !(cursor = strchr(cursor+1, ' ')) ) pdfids = false;
713  if (pdfids) pi->pdf_id[0] = atoi(cursor);
714  else pi->pdf_id[0] = 0;
715 
716  if (pdfids) if ( !(cursor = strchr(cursor+1, ' ')) ) pdfids = false;
717  if (pdfids) pi->pdf_id[1] = atoi(cursor);
718  else pi->pdf_id[1] = 0;
719 
720  evt.add_attribute("GenPdfInfo", pi);
721 
722  return true;
723 }
724 bool ReaderAsciiHepMC2::failed() { return m_isstream ? (bool)m_stream->rdstate() :(bool)m_file.rdstate(); }
725 
727  if (m_event_ghost) { m_event_ghost->clear(); delete m_event_ghost; m_event_ghost=nullptr;}
728  if ( !m_file.is_open() ) return;
729  m_file.close();
730 }
731 
732 } // namespace HepMC3
#define HEPMC3_WARNING(MESSAGE)
Macro for printing HEPMC3_HEPMC3_WARNING messages.
Definition: Errors.h:27
#define HEPMC3_DEBUG(LEVEL, MESSAGE)
Macro for printing debug messages with appropriate debug level.
Definition: Errors.h:33
#define HEPMC3_ERROR(MESSAGE)
Macro for printing error messages.
Definition: Errors.h:24
Definition of class GenEvent.
Definition of attribute class GenHeavyIon.
Definition of class GenParticle.
Definition of event attribute class GenPdfInfo.
Definition of class GenVertex.
Definition of class ReaderAsciiHepMC2.
Definition of class Setup.
Generic 4-vector.
Definition: FourVector.h:36
Stores event-related information.
Definition: GenEvent.h:41
std::shared_ptr< T > attribute(const std::string &name, const int &id=0) const
Get attribute of type T.
Definition: GenEvent.h:409
void add_tree(const std::vector< GenParticlePtr > &particles)
Add whole tree in topological order.
Definition: GenEvent.cc:265
void add_vertex(GenVertexPtr v)
Add vertex.
Definition: GenEvent.cc:96
int event_number() const
Get event number.
Definition: GenEvent.h:148
void add_particle(GenParticlePtr p)
Add particle.
Definition: GenEvent.cc:48
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:391
void set_event_number(const int &num)
Set event number.
Definition: GenEvent.h:150
std::map< std::string, std::map< int, std::shared_ptr< Attribute > > > attributes() const
Get a copy of the list of attributes.
Definition: GenEvent.h:257
void set_run_info(std::shared_ptr< GenRunInfo > run)
Set the GenRunInfo object by smart pointer.
Definition: GenEvent.h:141
void add_attribute(const std::string &name, const std::shared_ptr< Attribute > &att, const int &id=0)
Definition: GenEvent.cc:806
void remove_attribute(const std::string &name, const int &id=0)
Remove attribute.
Definition: GenEvent.cc:609
const Units::LengthUnit & length_unit() const
Get length unit.
Definition: GenEvent.h:155
void clear()
Remove contents of this event.
Definition: GenEvent.cc:599
const std::vector< double > & weights() const
Get event weight values as a vector.
Definition: GenEvent.h:98
const Units::MomentumUnit & momentum_unit() const
Get momentum unit.
Definition: GenEvent.h:153
void reserve(const size_t &particles, const size_t &vertices=0)
Reserve memory for particles and vertices.
Definition: GenEvent.cc:385
Attribute that holds an Integer implemented as an int.
Definition: Attribute.h:157
int value() const
get the value associated to this Attribute.
Definition: Attribute.h:179
bool m_isstream
toggles usage of m_file or m_stream
bool read_event(GenEvent &evt) override
Implementation of Reader::read_event.
int parse_event_information(GenEvent &evt, const char *buf)
Parse event.
std::vector< int > m_vertex_barcodes
Old vertex barcodes.
bool failed() override
Return status of the stream.
bool parse_pdf_info(GenEvent &evt, const char *buf)
Parse pdf information.
bool skip(const int) override
skip events
std::ifstream m_file
Input file.
bool parse_units(GenEvent &evt, const char *buf)
Parse units.
int parse_particle_information(const char *buf)
Parse particle.
std::vector< GenParticlePtr > m_particle_cache_ghost
Particle cache for attributes.
void close() override
Close file stream.
std::vector< GenVertexPtr > m_vertex_cache
Vertex cache.
ReaderAsciiHepMC2(const std::string &filename)
Default constructor.
std::vector< GenParticlePtr > m_particle_cache
Particle cache.
bool parse_weight_names(const char *buf)
Parse weight names.
bool parse_xs_info(GenEvent &evt, const char *buf)
Parse pdf information.
std::istream * m_stream
For ctor when reading from stream.
std::vector< GenVertexPtr > m_vertex_cache_ghost
Vertex cache for attributes.
int parse_vertex_information(const char *buf)
Parse vertex.
bool parse_heavy_ion(GenEvent &evt, const char *buf)
Parse heavy ion information.
GenEvent * m_event_ghost
To save particle and verstex attributes.
std::vector< int > m_end_vertex_barcodes
Old end vertex barcodes.
std::map< std::string, std::string > m_options
options
Definition: Reader.h:68
void set_run_info(std::shared_ptr< GenRunInfo > run)
Set the global GenRunInfo object.
Definition: Reader.h:64
std::shared_ptr< GenRunInfo > run_info() const
Get the global GenRunInfo object.
Definition: Reader.h:44
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
Attribute that holds a vector of integers of type int.
Definition: Attribute.h:1046
std::vector< long int > value() const
get the value associated to this Attribute.
Definition: Attribute.h:1072
HepMC3 main namespace.