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
std::ifstream m_file
Input file.
std::map< std::string, std::string > m_options
options
Definition: Reader.h:68
#define HEPMC3_WARNING(MESSAGE)
Macro for printing HEPMC3_HEPMC3_WARNING messages.
Definition: Errors.h:27
void add_vertex(GenVertexPtr v)
Add vertex.
Definition: GenEvent.cc:96
Attribute that holds a vector of FPs of type double.
Definition: Attribute.h:1091
bool m_isstream
toggles usage of m_file or m_stream
static LengthUnit length_unit(const std::string &name)
Get length unit based on its name.
Definition: Units.h:46
Definition of class GenParticle.
std::vector< GenParticlePtr > m_particle_cache_ghost
Particle cache for attributes.
Definition of attribute class GenHeavyIon.
int parse_particle_information(const char *buf)
Parse particle.
Definition of class GenVertex.
bool failed() override
Return status of the stream.
Attribute that holds a vector of integers of type int.
Definition: Attribute.h:1001
std::vector< double > value() const
get the value associated to this Attribute.
Definition: Attribute.h:1117
#define HEPMC3_DEBUG(LEVEL, MESSAGE)
Macro for printing debug messages with appropriate debug level.
Definition: Errors.h:33
const Units::LengthUnit & length_unit() const
Get length unit.
Definition: GenEvent.h:155
GenEvent * m_event_ghost
To save particle and verstex attributes.
bool parse_xs_info(GenEvent &evt, const char *buf)
Parse pdf information.
void add_particle(GenParticlePtr p)
Add particle.
Definition: GenEvent.cc:48
static std::string name(MomentumUnit u)
Get name of momentum unit.
Definition: Units.h:56
int event_number() const
Get event number.
Definition: GenEvent.h:148
void add_attribute(const std::string &name, const std::shared_ptr< Attribute > &att, const int &id=0)
Definition: GenEvent.cc:806
void close() override
Close file stream.
std::vector< int > value() const
get the value associated to this Attribute.
Definition: Attribute.h:1027
ReaderAsciiHepMC2(const std::string &filename)
Default constructor.
LengthUnit
Position units.
Definition: Units.h:32
Definition of class ReaderAsciiHepMC2.
MomentumUnit
Momentum units.
Definition: Units.h:29
const Units::MomentumUnit & momentum_unit() const
Get momentum unit.
Definition: GenEvent.h:153
Stores event-related information.
Definition: GenEvent.h:41
Generic 4-vector.
Definition: FourVector.h:36
bool parse_weight_names(const char *buf)
Parse weight names.
std::vector< int > m_end_vertex_barcodes
Old end vertex barcodes.
std::shared_ptr< GenRunInfo > run_info() const
Get the global GenRunInfo object.
Definition: Reader.h:44
Definition of class Setup.
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 remove_attribute(const std::string &name, const int &id=0)
Remove attribute.
Definition: GenEvent.cc:609
static MomentumUnit momentum_unit(const std::string &name)
Get momentum unit based on its name.
Definition: Units.h:36
void set_run_info(std::shared_ptr< GenRunInfo > run)
Set the GenRunInfo object by smart pointer.
Definition: GenEvent.h:141
void reserve(const size_t &particles, const size_t &vertices=0)
Reserve memory for particles and vertices.
Definition: GenEvent.cc:385
std::istream * m_stream
For ctor when reading from stream.
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
Definition of event attribute class GenPdfInfo.
void set_run_info(std::shared_ptr< GenRunInfo > run)
Set the global GenRunInfo object.
Definition: Reader.h:64
#define HEPMC3_ERROR(MESSAGE)
Macro for printing error messages.
Definition: Errors.h:24
std::vector< long int > value() const
get the value associated to this Attribute.
Definition: Attribute.h:1072
const std::vector< double > & weights() const
Get event weight values as a vector.
Definition: GenEvent.h:98
Attribute that holds a vector of integers of type int.
Definition: Attribute.h:1046
std::vector< GenVertexPtr > m_vertex_cache_ghost
Vertex cache for attributes.
Definition of class GenEvent.
bool read_event(GenEvent &evt) override
Implementation of Reader::read_event.
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
std::vector< int > m_vertex_barcodes
Old vertex barcodes.
std::vector< GenVertexPtr > m_vertex_cache
Vertex cache.
std::vector< GenParticlePtr > m_particle_cache
Particle cache.
int parse_event_information(GenEvent &evt, const char *buf)
Parse event.
void clear()
Remove contents of this event.
Definition: GenEvent.cc:599
bool parse_units(GenEvent &evt, const char *buf)
Parse units.
int value() const
get the value associated to this Attribute.
Definition: Attribute.h:179
bool parse_heavy_ion(GenEvent &evt, const char *buf)
Parse heavy ion information.
bool parse_pdf_info(GenEvent &evt, const char *buf)
Parse pdf information.
int parse_vertex_information(const char *buf)
Parse vertex.
Attribute that holds an Integer implemented as an int.
Definition: Attribute.h:157
bool skip(const int) override
skip events
void set_event_number(const int &num)
Set event number.
Definition: GenEvent.h:150