24 m_file(filename), m_stream(0), m_isstream(false) {
26 HEPMC3_ERROR(
"ReaderAsciiHepMC2: could not open input file: "<<filename )
33 : m_stream(&stream), m_isstream(true)
36 HEPMC3_ERROR(
"ReaderAsciiHepMC2: could not open input stream " )
46 const size_t max_buffer_size=512*512;
47 char buf[max_buffer_size];
54 if (nn<0)
return true;
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;
89 if( strlen(buf) == 0 )
continue;
91 if( strncmp(buf,
"HepMC",5) == 0 ) {
92 if( strncmp(buf,
"HepMC::Version",14) != 0 && strncmp(buf,
"HepMC::IO_GenEvent",18)!=0 )
94 HEPMC3_WARNING(
"ReaderAsciiHepMC2: found unsupported expression in header. Will close the input." )
95 std::cout<<buf<<std::endl;
98 if(parsed_event_header) {
99 is_parsing_successful =
true;
107 if(parsing_result<0) {
108 is_parsing_successful =
false;
109 HEPMC3_ERROR(
"ReaderAsciiHepMC2: HEPMC3_ERROR parsing event information" )
112 vertices_count = parsing_result;
117 is_parsing_successful =
true;
119 parsed_event_header =
true;
126 if(current_vertex_particles_parsed < current_vertex_particles_count) {
127 is_parsing_successful =
false;
130 current_vertex_particles_parsed = 0;
134 if(parsing_result<0) {
135 is_parsing_successful =
false;
136 HEPMC3_ERROR(
"ReaderAsciiHepMC2: HEPMC3_ERROR parsing vertex information" )
139 current_vertex_particles_count = parsing_result;
140 is_parsing_successful =
true;
147 if(parsing_result<0) {
148 is_parsing_successful =
false;
149 HEPMC3_ERROR(
"ReaderAsciiHepMC2: HEPMC3_ERROR parsing particle information" )
152 ++current_vertex_particles_parsed;
153 is_parsing_successful =
true;
172 HEPMC3_WARNING(
"ReaderAsciiHepMC2: skipping unrecognised prefix: " << buf[0] )
173 is_parsing_successful =
true;
177 if( !is_parsing_successful )
break;
181 if( parsed_event_header && peek==
'E' )
break;
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;
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;
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 )
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);
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() );
271 for (
size_t ii=0; ii<max_weights_size; ii++)
275 m_vertex_cache[i]->add_attribute(
"weight"+std::to_string((
long long unsigned int)ii),rs);
285 const char *cursor = buf;
287 int vertices_count = 0;
288 int random_states_size = 0;
289 int weights_size = 0;
290 std::vector<long> random_states(0);
291 std::vector<double> weights(0);
294 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
295 event_no = atoi(cursor);
299 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
300 std::shared_ptr<IntAttribute> mpi = std::make_shared<IntAttribute>(atoi(cursor));
304 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
305 std::shared_ptr<DoubleAttribute> event_scale = std::make_shared<DoubleAttribute>(atof(cursor));
309 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
310 std::shared_ptr<DoubleAttribute> alphaQCD = std::make_shared<DoubleAttribute>(atof(cursor));
314 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
315 std::shared_ptr<DoubleAttribute> alphaQED = std::make_shared<DoubleAttribute>(atof(cursor));
319 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
320 std::shared_ptr<IntAttribute> signal_process_id = std::make_shared<IntAttribute>(atoi(cursor));
324 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
325 std::shared_ptr<IntAttribute> signal_process_vertex = std::make_shared<IntAttribute>(atoi(cursor));
326 evt.
add_attribute(
"signal_process_vertex",signal_process_vertex);
329 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
330 vertices_count = atoi(cursor);
333 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
336 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
339 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
340 random_states_size = atoi(cursor);
341 random_states.resize(random_states_size);
343 for (
int i = 0; i < random_states_size; ++i ) {
344 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
345 random_states[i] = atoi(cursor);
349 evt.
add_attribute(
"random_states",std::make_shared<VectorLongIntAttribute>(random_states));
353 for (
int i = 0; i < random_states_size; ++i )
354 evt.
add_attribute(
"random_states"+std::to_string((
long long unsigned int)i),std::make_shared<IntAttribute>(random_states[i]));
357 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
358 weights_size = atoi(cursor);
359 weights.resize(weights_size);
361 for (
int i = 0; i < weights_size; ++i ) {
362 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
363 weights[i] = atof(cursor);
368 HEPMC3_DEBUG( 10,
"ReaderAsciiHepMC2: E: "<<event_no<<
" ("<<vertices_count<<
"V, "<<weights_size<<
"W, "<<random_states_size<<
"RS)" )
370 return vertices_count;
374 const char *cursor = buf;
377 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
382 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
386 evt.
set_units(momentum_unit,length_unit);
394 GenVertexPtr data = std::make_shared<GenVertex>();
395 GenVertexPtr data_ghost = std::make_shared<GenVertex>();
397 const char *cursor = buf;
399 int num_particles_out = 0;
400 int weights_size = 0;
401 std::vector<double> weights(0);
403 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
404 barcode = atoi(cursor);
407 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
408 data->set_status( atoi(cursor) );
411 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
412 position.
setX(atof(cursor));
415 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
416 position.
setY(atof(cursor));
419 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
420 position.
setZ(atof(cursor));
423 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
424 position.
setT(atof(cursor));
425 data->set_position( position );
428 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
431 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
432 num_particles_out = atoi(cursor);
436 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
437 weights_size = atoi(cursor);
438 weights.resize(weights_size);
440 for (
int i = 0; i < weights_size; ++i ) {
441 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
442 weights[i] = atof(cursor);
454 for (
int i = 0; i < weights_size; ++i )
455 data_ghost->add_attribute(
"weight"+std::to_string((
long long unsigned int)i),std::make_shared<DoubleAttribute>(weights[i]));
459 data_ghost->add_attribute(
"weights",std::make_shared<VectorDoubleAttribute>(weights));
460 data_ghost->add_attribute(
"weights",std::make_shared<VectorDoubleAttribute>(weights));
464 HEPMC3_DEBUG( 10,
"ReaderAsciiHepMC2: V: "<<-(
int)
m_vertex_cache.size()<<
" (old barcode"<<barcode<<
") "<<num_particles_out<<
" particles)" )
466 return num_particles_out;
470 GenParticlePtr data = std::make_shared<GenParticle>();
471 GenParticlePtr data_ghost = std::make_shared<GenParticle>();
474 const char *cursor = buf;
478 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
481 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
482 data->set_pid( atoi(cursor) );
485 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
486 momentum.
setPx(atof(cursor));
489 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
490 momentum.
setPy(atof(cursor));
493 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
494 momentum.
setPz(atof(cursor));
497 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
498 momentum.
setE(atof(cursor));
499 data->set_momentum(momentum);
502 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
503 data->set_generated_mass( atof(cursor) );
506 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
507 data->set_status( atoi(cursor) );
510 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
511 std::shared_ptr<DoubleAttribute> theta = std::make_shared<DoubleAttribute>(atof(cursor));
512 if (theta->value()!=0.0) data_ghost->add_attribute(
"theta",theta);
515 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
516 std::shared_ptr<DoubleAttribute> phi = std::make_shared<DoubleAttribute>(atof(cursor));
517 if (phi->value()!=0.0) data_ghost->add_attribute(
"phi",phi);
520 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
521 end_vtx = atoi(cursor);
524 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
525 int flowsize=atoi(cursor);
527 std::map<int,int> flows;
528 for (
int i=0; i<flowsize; i++)
530 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
531 int flowindex=atoi(cursor);
532 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
533 int flowvalue=atoi(cursor);
534 flows[flowindex]=flowvalue;
538 std::vector<int> vectorflows;
539 for (
auto f: flows) vectorflows.push_back(f.second);
540 data_ghost->add_attribute(
"flows",std::make_shared<VectorIntAttribute>(vectorflows));
544 for (
auto f: flows) data_ghost->add_attribute(
"flow"+std::to_string((
long long int)f.first),std::make_shared<IntAttribute>(f.second));
565 const char *cursor = buf;
566 std::shared_ptr<GenCrossSection> xs = std::make_shared<GenCrossSection>();
568 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
569 double xs_val = atof(cursor);
571 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
572 double xs_err = atof(cursor);
574 xs->set_cross_section( xs_val, xs_err);
581 const char *cursor = buf;
582 const char *cursor2 = buf;
584 std::vector<std::string> w_names;
589 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
590 w_count = atoi(cursor);
592 if( w_count <= 0 )
return false;
594 w_names.resize(w_count);
596 for(
int i=0; i < w_count; ++i ) {
598 if( !(cursor = strchr(cursor+1,
'"')) )
return false;
599 if( !(cursor2 = strchr(cursor+1,
'"')) )
return false;
604 w_names[i].assign(cursor, cursor2-cursor);
609 run_info()->set_weight_names(w_names);
615 std::shared_ptr<GenHeavyIon> hi = std::make_shared<GenHeavyIon>();
616 const char *cursor = buf;
618 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
619 hi->Ncoll_hard = atoi(cursor);
621 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
622 hi->Npart_proj = atoi(cursor);
624 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
625 hi->Npart_targ = atoi(cursor);
627 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
628 hi->Ncoll = atoi(cursor);
630 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
631 hi->spectator_neutrons = atoi(cursor);
633 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
634 hi->spectator_protons = atoi(cursor);
636 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
637 hi->N_Nwounded_collisions = atoi(cursor);
639 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
640 hi->Nwounded_N_collisions = atoi(cursor);
642 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
643 hi->Nwounded_Nwounded_collisions = atoi(cursor);
645 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
646 hi->impact_parameter = atof(cursor);
648 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
649 hi->event_plane_angle = atof(cursor);
651 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
652 hi->eccentricity = atof(cursor);
654 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
655 hi->sigma_inel_NN = atof(cursor);
658 hi->centrality = 0.0;
666 std::shared_ptr<GenPdfInfo> pi = std::make_shared<GenPdfInfo>();
667 const char *cursor = buf;
669 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
670 pi->parton_id[0] = atoi(cursor);
672 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
673 pi->parton_id[1] = atoi(cursor);
675 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
676 pi->x[0] = atof(cursor);
678 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
679 pi->x[1] = atof(cursor);
681 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
682 pi->scale = atof(cursor);
684 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
685 pi->xf[0] = atof(cursor);
687 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
688 pi->xf[1] = atof(cursor);
692 if( !(cursor = strchr(cursor+1,
' ')) ) pdfids=
false;
693 if(pdfids) pi->pdf_id[0] = atoi(cursor);
694 else pi->pdf_id[0] =0;
696 if(pdfids)
if( !(cursor = strchr(cursor+1,
' ')) ) pdfids=
false;
697 if(pdfids) pi->pdf_id[1] = atoi(cursor);
698 else pi->pdf_id[1] =0;
708 if( !
m_file.is_open() )
return;
#define HEPMC3_WARNING(MESSAGE)
Macro for printing HEPMC3_HEPMC3_WARNING messages.
#define HEPMC3_DEBUG(LEVEL, MESSAGE)
Macro for printing debug messages with appropriate debug level.
#define HEPMC3_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 set_event_number(const int &num)
Set event number.
void set_run_info(std::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.
void add_attribute(const std::string &name, const std::shared_ptr< Attribute > &att, const int &id=0)
Add event attribute to event.
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
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.
~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.
std::vector< int > m_end_vertex_barcodes
Old end vertex barcodes.
std::shared_ptr< GenRunInfo > run_info() const
Get the global GenRunInfo object.
std::map< std::string, std::string > m_options
options
void set_run_info(std::shared_ptr< GenRunInfo > run)
Set 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.
Attribute that holds a vector of integers of type int.