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