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