26 m_file(filename), m_stream(0), m_isstream(false) {
28 ERROR(
"ReaderAsciiHepMC2: could not open input file: "<<filename )
35 : m_stream(&stream), m_isstream(true)
38 ERROR(
"ReaderAsciiHepMC2: could not open input stream " )
51 const size_t max_buffer_size=512*512;
52 const size_t max_weights_size=256;
53 char buf[max_buffer_size];
54 bool parsed_event_header =
false;
55 bool is_parsing_successful =
true;
56 int parsing_result = 0;
57 unsigned int vertices_count = 0;
58 unsigned int current_vertex_particles_count = 0;
59 unsigned int current_vertex_particles_parsed= 0;
76 if( strlen(buf) == 0 )
continue;
78 if( strncmp(buf,
"HepMC",5) == 0 ) {
79 if( strncmp(buf,
"HepMC::Version",14) != 0 && strncmp(buf,
"HepMC::IO_GenEvent",18)!=0 )
81 WARNING(
"ReaderAsciiHepMC2: found unsupported expression in header. Will close the input." )
82 std::cout<<buf<<std::endl;
85 if(parsed_event_header) {
86 is_parsing_successful =
true;
94 if(parsing_result<0) {
95 is_parsing_successful =
false;
96 ERROR(
"ReaderAsciiHepMC2: error parsing event information" )
99 vertices_count = parsing_result;
104 is_parsing_successful =
true;
106 parsed_event_header =
true;
113 if(current_vertex_particles_parsed < current_vertex_particles_count) {
114 is_parsing_successful =
false;
117 current_vertex_particles_parsed = 0;
121 if(parsing_result<0) {
122 is_parsing_successful =
false;
123 ERROR(
"ReaderAsciiHepMC2: error parsing vertex information" )
126 current_vertex_particles_count = parsing_result;
127 is_parsing_successful =
true;
134 if(parsing_result<0) {
135 is_parsing_successful =
false;
136 ERROR(
"ReaderAsciiHepMC2: error parsing particle information" )
139 ++current_vertex_particles_parsed;
140 is_parsing_successful =
true;
159 WARNING(
"ReaderAsciiHepMC2: skipping unrecognised prefix: " << buf[0] )
160 is_parsing_successful =
true;
164 if( !is_parsing_successful )
break;
168 if( parsed_event_header && peek==
'E' )
break;
174 if( is_parsing_successful && current_vertex_particles_parsed < current_vertex_particles_count ) {
175 ERROR(
"ReaderAsciiHepMC2: not all particles parsed" )
176 is_parsing_successful =
false;
179 else if( is_parsing_successful &&
m_vertex_cache.size() != vertices_count ) {
180 ERROR(
"ReaderAsciiHepMC2: not all vertices parsed" )
181 is_parsing_successful =
false;
184 if( !is_parsing_successful ) {
185 ERROR(
"ReaderAsciiHepMC2: event parsing failed. Returning empty event" )
186 DEBUG( 1,
"Parsing failed at line:" << std::endl << buf )
237 for (
size_t ii=0; ii<max_weights_size; ii++)
241 m_vertex_cache[i]->add_attribute(
"weight"+to_string((
long long unsigned int)ii),rs);
251 const char *cursor = buf;
253 int vertices_count = 0;
254 int random_states_size = 0;
255 int weights_size = 0;
256 std::vector<long> random_states(0);
257 std::vector<double> weights(0);
260 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
261 event_no = atoi(cursor);
265 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
266 shared_ptr<IntAttribute> mpi = make_shared<IntAttribute>(atoi(cursor));
270 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
271 shared_ptr<DoubleAttribute> event_scale = make_shared<DoubleAttribute>(atof(cursor));
275 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
276 shared_ptr<DoubleAttribute> alphaQCD = make_shared<DoubleAttribute>(atof(cursor));
280 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
281 shared_ptr<DoubleAttribute> alphaQED = make_shared<DoubleAttribute>(atof(cursor));
285 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
286 shared_ptr<IntAttribute> signal_process_id = make_shared<IntAttribute>(atoi(cursor));
290 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
291 shared_ptr<IntAttribute> signal_process_vertex = make_shared<IntAttribute>(atoi(cursor));
292 evt.
add_attribute(
"signal_process_vertex",signal_process_vertex);
295 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
296 vertices_count = atoi(cursor);
299 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
302 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
305 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
306 random_states_size = atoi(cursor);
307 random_states.resize(random_states_size);
309 for (
int i = 0; i < random_states_size; ++i ) {
310 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
311 random_states[i] = atoi(cursor);
314 for (
int i = 0; i < random_states_size; ++i )
315 evt.
add_attribute(
"random_states"+to_string((
long long unsigned int)i),make_shared<IntAttribute>(random_states[i]));
318 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
319 weights_size = atoi(cursor);
320 weights.resize(weights_size);
322 for (
int i = 0; i < weights_size; ++i ) {
323 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
324 weights[i] = atof(cursor);
329 DEBUG( 10,
"ReaderAsciiHepMC2: E: "<<event_no<<
" ("<<vertices_count<<
"V, "<<weights_size<<
"W, "<<random_states_size<<
"RS)" )
331 return vertices_count;
335 const char *cursor = buf;
338 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
343 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
347 evt.
set_units(momentum_unit,length_unit);
355 GenVertexPtr data = make_shared<GenVertex>();
356 GenVertexPtr data_ghost = make_shared<GenVertex>();
358 const char *cursor = buf;
360 int num_particles_out = 0;
361 int weights_size = 0;
362 std::vector<double> weights(0);
364 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
365 barcode = atoi(cursor);
368 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
369 data->set_status( atoi(cursor) );
372 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
373 position.
setX(atof(cursor));
376 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
377 position.
setY(atof(cursor));
380 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
381 position.
setZ(atof(cursor));
384 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
385 position.
setT(atof(cursor));
386 data->set_position( position );
389 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
392 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
393 num_particles_out = atoi(cursor);
397 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
398 weights_size = atoi(cursor);
399 weights.resize(weights_size);
401 for (
int i = 0; i < weights_size; ++i ) {
402 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
403 weights[i] = atof(cursor);
413 for (
int i = 0; i < weights_size; ++i )
414 data_ghost->add_attribute(
"weight"+to_string((
long long unsigned int)i),make_shared<DoubleAttribute>(weights[i]));
417 DEBUG( 10,
"ReaderAsciiHepMC2: V: "<<-(
int)
m_vertex_cache.size()<<
" (old barcode"<<barcode<<
") "<<num_particles_out<<
" particles)" )
419 return num_particles_out;
423 GenParticlePtr data = make_shared<GenParticle>();
424 GenParticlePtr data_ghost = make_shared<GenParticle>();
427 const char *cursor = buf;
431 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
434 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
435 data->set_pid( atoi(cursor) );
438 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
439 momentum.
setPx(atof(cursor));
442 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
443 momentum.
setPy(atof(cursor));
446 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
447 momentum.
setPz(atof(cursor));
450 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
451 momentum.
setE(atof(cursor));
452 data->set_momentum(momentum);
455 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
456 data->set_generated_mass( atof(cursor) );
459 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
460 data->set_status( atoi(cursor) );
463 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
464 shared_ptr<DoubleAttribute> theta = make_shared<DoubleAttribute>(atof(cursor));
465 if (theta->value()!=0.0) data_ghost->add_attribute(
"theta",theta);
468 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
469 shared_ptr<DoubleAttribute> phi = make_shared<DoubleAttribute>(atof(cursor));
470 if (phi->value()!=0.0) data_ghost->add_attribute(
"phi",phi);
473 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
474 end_vtx = atoi(cursor);
477 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
478 int flowsize=atoi(cursor);
480 for (
int i=0; i<flowsize; i++)
482 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
483 int flowindex=atoi(cursor);
484 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
485 int flowvalue=atoi(cursor);
486 data_ghost->add_attribute(
"flow"+to_string((
long long int)flowindex),make_shared<IntAttribute>(flowvalue));
502 DEBUG( 10,
"ReaderAsciiHepMC2: P: "<<
m_particle_cache.size()<<
" ( pid: "<<data->pid()<<
") end vertex: "<<end_vtx )
508 const char *cursor = buf;
509 shared_ptr<GenCrossSection> xs = make_shared<GenCrossSection>();
511 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
512 double xs_val = atof(cursor);
514 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
515 double xs_err = atof(cursor);
517 xs->set_cross_section( xs_val , xs_err);
524 const char *cursor = buf;
525 const char *cursor2 = buf;
527 vector<string> w_names;
532 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
533 w_count = atoi(cursor);
535 if( w_count <= 0 )
return false;
537 w_names.resize(w_count);
539 for(
int i=0; i < w_count; ++i ) {
541 if( !(cursor = strchr(cursor+1,
'"')) )
return false;
542 if( !(cursor2 = strchr(cursor+1,
'"')) )
return false;
547 w_names[i].assign(cursor, cursor2-cursor);
552 run_info()->set_weight_names(w_names);
558 shared_ptr<GenHeavyIon> hi = make_shared<GenHeavyIon>();
559 const char *cursor = buf;
561 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
562 hi->Ncoll_hard = atoi(cursor);
564 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
565 hi->Npart_proj = atoi(cursor);
567 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
568 hi->Npart_targ = atoi(cursor);
570 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
571 hi->Ncoll = atoi(cursor);
573 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
574 hi->spectator_neutrons = atoi(cursor);
576 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
577 hi->spectator_protons = atoi(cursor);
579 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
580 hi->N_Nwounded_collisions = atoi(cursor);
582 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
583 hi->Nwounded_N_collisions = atoi(cursor);
585 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
586 hi->Nwounded_Nwounded_collisions = atoi(cursor);
588 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
589 hi->impact_parameter = atof(cursor);
591 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
592 hi->event_plane_angle = atof(cursor);
594 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
595 hi->eccentricity = atof(cursor);
597 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
598 hi->sigma_inel_NN = atof(cursor);
601 hi->centrality = 0.0;
609 shared_ptr<GenPdfInfo> pi = make_shared<GenPdfInfo>();
610 const char *cursor = buf;
612 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
613 pi->parton_id[0] = atoi(cursor);
615 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
616 pi->parton_id[1] = atoi(cursor);
618 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
619 pi->x[0] = atof(cursor);
621 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
622 pi->x[1] = atof(cursor);
624 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
625 pi->scale = atof(cursor);
627 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
628 pi->xf[0] = atof(cursor);
630 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
631 pi->xf[1] = atof(cursor);
635 if( !(cursor = strchr(cursor+1,
' ')) ) pdfids=
false;
636 if(pdfids) pi->pdf_id[0] = atoi(cursor);
637 else pi->pdf_id[0] =0;
639 if(pdfids)
if( !(cursor = strchr(cursor+1,
' ')) ) pdfids=
false;
640 if(pdfids) pi->pdf_id[1] = atoi(cursor);
641 else pi->pdf_id[1] =0;
651 if( !
m_file.is_open() )
return;
#define WARNING(MESSAGE)
Macro for printing warning messages.
#define DEBUG(LEVEL, MESSAGE)
Macro for printing debug messages with appropriate debug level.
#define ERROR(MESSAGE)
Macro for printing error messages.
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.
Stores event-related information.
void add_tree(const std::vector< GenParticlePtr > &particles)
Add whole tree in topological order.
void add_vertex(GenVertexPtr v)
Add vertex.
void add_particle(GenParticlePtr p)
Add particle.
void set_units(Units::MomentumUnit new_momentum_unit, Units::LengthUnit new_length_unit)
Change event units Converts event from current units to new ones.
void add_attribute(const string &name, const shared_ptr< Attribute > &att, const int &id=0)
Add event attribute to event.
void set_event_number(const int &num)
Set event number.
void set_run_info(shared_ptr< GenRunInfo > run)
Set the GenRunInfo object by smart pointer.
const Units::MomentumUnit & momentum_unit() const
Get momentum unit.
const Units::LengthUnit & length_unit() const
Get length unit.
const std::vector< double > & weights() const
Get event weight values as a vector.
void clear()
Remove contents of this event.
void reserve(const size_t &particles, const size_t &vertices=0)
Reserve memory for particles and vertices.
Attribute that holds an Integer implemented as an int.
bool m_isstream
toggles usage of m_file or m_stream
vector< int > m_vertex_barcodes
Old vertex barcodes.
bool read_event(GenEvent &evt)
Implementation of Reader::read_event.
int parse_event_information(GenEvent &evt, const char *buf)
Parse event.
vector< int > m_end_vertex_barcodes
Old end vertex barcodes.
bool parse_pdf_info(GenEvent &evt, const char *buf)
Parse pdf information.
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.
void close()
Close file stream.
ReaderAsciiHepMC2(const std::string &filename)
Default constructor.
vector< GenParticlePtr > m_particle_cache
Particle cache.
bool parse_weight_names(const char *buf)
Parse weight names.
bool failed()
Return status of the stream.
bool parse_xs_info(GenEvent &evt, const char *buf)
Parse pdf information.
vector< GenVertexPtr > m_vertex_cache
Vertex cache.
std::istream * m_stream
For ctor when reading from stdin.
vector< GenParticlePtr > m_particle_cache_ghost
Particle cache for attributes.
vector< GenVertexPtr > m_vertex_cache_ghost
Vertex cache for attributes.
~ReaderAsciiHepMC2()
Destructor.
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.
void set_run_info(shared_ptr< GenRunInfo > run)
Set the global GenRunInfo object.
shared_ptr< GenRunInfo > run_info() const
Get the global GenRunInfo object.
static LengthUnit length_unit(const std::string &name)
Get length unit based on its name.
LengthUnit
Position units.
static std::string name(MomentumUnit u)
Get name of momentum unit.
MomentumUnit
Momentum units.
static MomentumUnit momentum_unit(const std::string &name)
Get momentum unit based on its name.