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
#define HEPMC3_WARNING(MESSAGE)
Macro for printing HEPMC3_HEPMC3_WARNING messages.
Definition: Errors.h:26
#define HEPMC3_DEBUG(LEVEL, MESSAGE)
Macro for printing debug messages with appropriate debug level.
Definition: Errors.h:32
#define HEPMC3_ERROR(MESSAGE)
Macro for printing error messages.
Definition: Errors.h:23
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.
Attribute that holds a real number as a double.
Definition: Attribute.h:242
Generic 4-vector.
Definition: FourVector.h:35
void setE(double ee)
Definition: FourVector.h:134
void setT(double tt)
Definition: FourVector.h:105
void setPz(double pzz)
Definition: FourVector.h:127
void setY(double yy)
Definition: FourVector.h:91
void setPy(double pyy)
Definition: FourVector.h:120
void setX(double xx)
Definition: FourVector.h:84
void setPx(double pxx)
Definition: FourVector.h:113
void setZ(double zz)
Definition: FourVector.h:98
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:389
void add_tree(const std::vector< GenParticlePtr > &particles)
Add whole tree in topological order.
Definition: GenEvent.cc:267
void add_vertex(GenVertexPtr v)
Add vertex.
Definition: GenEvent.cc:97
void add_particle(GenParticlePtr p)
Add particle.
Definition: GenEvent.cc:49
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
void set_event_number(const int &num)
Set event number.
Definition: GenEvent.h:137
void set_run_info(std::shared_ptr< GenRunInfo > run)
Set the GenRunInfo object by smart pointer.
Definition: GenEvent.h:128
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
const Units::LengthUnit & length_unit() const
Get length unit.
Definition: GenEvent.h:142
void clear()
Remove contents of this event.
Definition: GenEvent.cc:608
const std::vector< double > & weights() const
Get event weight values as a vector.
Definition: GenEvent.h:86
const Units::MomentumUnit & momentum_unit() const
Get momentum unit.
Definition: GenEvent.h:140
void reserve(const size_t &particles, const size_t &vertices=0)
Reserve memory for particles and vertices.
Definition: GenEvent.cc:389
Attribute that holds an Integer implemented as an int.
Definition: Attribute.h:158
int value() const
get the value associated to this Attribute.
Definition: Attribute.h:180
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 stdin.
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:1002
HepMC3 main namespace.