26 m_file(filename), m_stream(nullptr), m_isstream(false) {
28 HEPMC3_ERROR(
"ReaderAsciiHepMC2: could not open input file: " << filename )
35 : m_stream(&stream), m_isstream(true)
38 HEPMC3_ERROR(
"ReaderAsciiHepMC2: could not open input stream ")
45 : m_shared_stream(s_stream), m_stream(s_stream.get()), m_isstream(true)
48 HEPMC3_ERROR(
"ReaderAsciiHepMC2: could not open input stream ")
58 const size_t max_buffer_size = 512*512;
59 char buf[max_buffer_size];
65 if ( peek ==
'E' ) nn--;
66 if (nn < 0)
return true;
76 const size_t max_buffer_size = 512*512;
77 char buf[max_buffer_size];
78 bool parsed_event_header =
false;
79 bool is_parsing_successful =
true;
80 int parsing_result = 0;
81 unsigned int vertices_count = 0;
82 unsigned int current_vertex_particles_count = 0;
83 unsigned int current_vertex_particles_parsed = 0;
101 if ( strlen(buf) == 0 )
continue;
103 if ( strncmp(buf,
"HepMC", 5) == 0 ) {
104 if ( strncmp(buf,
"HepMC::Version", 14) != 0 && strncmp(buf,
"HepMC::IO_GenEvent", 18) != 0 )
106 HEPMC3_WARNING(
"ReaderAsciiHepMC2: found unsupported expression in header. Will close the input.")
107 std::cout <<buf << std::endl;
110 if (parsed_event_header) {
111 is_parsing_successful =
true;
119 if (parsing_result < 0) {
120 is_parsing_successful =
false;
121 HEPMC3_ERROR(
"ReaderAsciiHepMC2: HEPMC3_ERROR parsing event information")
124 vertices_count = parsing_result;
133 is_parsing_successful =
true;
135 parsed_event_header =
true;
142 if (current_vertex_particles_parsed < current_vertex_particles_count) {
143 is_parsing_successful =
false;
146 current_vertex_particles_parsed = 0;
150 if (parsing_result < 0) {
151 is_parsing_successful =
false;
152 HEPMC3_ERROR(
"ReaderAsciiHepMC2: HEPMC3_ERROR parsing vertex information")
155 current_vertex_particles_count = parsing_result;
156 is_parsing_successful =
true;
163 if (parsing_result < 0) {
164 is_parsing_successful =
false;
165 HEPMC3_ERROR(
"ReaderAsciiHepMC2: HEPMC3_ERROR parsing particle information")
168 ++current_vertex_particles_parsed;
169 is_parsing_successful =
true;
188 HEPMC3_WARNING(
"ReaderAsciiHepMC2: skipping unrecognised prefix: " << buf[0])
189 is_parsing_successful =
true;
193 if ( !is_parsing_successful )
break;
197 if ( parsed_event_header && peek ==
'E' )
break;
203 if (is_parsing_successful && current_vertex_particles_parsed < current_vertex_particles_count) {
204 HEPMC3_ERROR(
"ReaderAsciiHepMC2: not all particles parsed")
205 is_parsing_successful =
false;
208 else if (is_parsing_successful &&
m_vertex_cache.size() != vertices_count) {
209 HEPMC3_ERROR(
"ReaderAsciiHepMC2: not all vertices parsed")
210 is_parsing_successful =
false;
213 if ( !is_parsing_successful ) {
214 HEPMC3_ERROR(
"ReaderAsciiHepMC2: event parsing failed. Returning empty event")
215 HEPMC3_DEBUG(1,
"Parsing failed at line:" << std::endl << buf)
238 std::vector<GenParticlePtr> beams;
239 for (
auto p:
m_vertex_cache[i]->particles_out())
if (p->status() == 4 && !(p->end_vertex())) beams.push_back(p);
244 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());
246 if (beams.size() == 0) {
266 if (random_states_a) {
267 std::vector<long int> random_states_v = random_states_a->
value();
268 for (
size_t i = 0; i < random_states_v.size(); ++i )
269 evt.
add_attribute(
"random_states" + std::to_string((
long long unsigned int)i), std::make_shared<IntAttribute>(random_states_v[i]));
275 std::map< std::string, std::map<int, std::shared_ptr<Attribute> > > cached_attributes =
m_event_ghost->
attributes();
276 if (cached_attributes.find(
"flows") != cached_attributes.end()) {
277 std::map<int, std::shared_ptr<Attribute> > flows = cached_attributes.at(
"flows");
281 for (
auto f: flows)
if (f.first > 0 && f.first <= (
int)
m_particle_cache.size()) {
282 std::shared_ptr<VectorIntAttribute> casted = std::dynamic_pointer_cast<VectorIntAttribute>(f.second);
283 if (!casted)
continue;
284 std::vector<int> this_p_flow = casted->value();
285 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]));
290 if (cached_attributes.find(
"phi") != cached_attributes.end()) {
291 std::map<int, std::shared_ptr<Attribute> > phi = cached_attributes.at(
"phi");
295 if (cached_attributes.find(
"theta") != cached_attributes.end()) {
296 std::map<int, std::shared_ptr<Attribute> > theta = cached_attributes.at(
"theta");
300 if (cached_attributes.find(
"weights") != cached_attributes.end()) {
301 std::map<int, std::shared_ptr<Attribute> > weights = cached_attributes.at(
"weights");
303 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);
305 for (
auto f: weights)
if (f.first < 0 && f.first >= -(
int)
m_vertex_cache.size()) {
306 std::shared_ptr<VectorDoubleAttribute> casted = std::dynamic_pointer_cast<VectorDoubleAttribute>(f.second);
307 if (!casted)
continue;
308 std::vector<double> this_v_weight = casted->value();
309 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]));
314 std::shared_ptr<IntAttribute> signal_process_vertex_barcode = evt.
attribute<
IntAttribute>(
"signal_process_vertex");
315 if (signal_process_vertex_barcode) {
316 int signal_process_vertex_barcode_value = signal_process_vertex_barcode->
value();
321 std::shared_ptr<IntAttribute> signal_process_vertex = std::make_shared<IntAttribute>(
m_vertex_cache.at(i)->id());
322 evt.
add_attribute(
"signal_process_vertex", signal_process_vertex);
333 const char *cursor = buf;
334 int vertices_count = 0;
335 int random_states_size = 0;
336 int weights_size = 0;
337 std::vector<long> random_states(0);
338 std::vector<double> weights(0);
341 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
345 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
346 evt.
add_attribute(
"mpi", std::make_shared<IntAttribute>(atoi(cursor)));
349 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
350 evt.
add_attribute(
"event_scale", std::make_shared<DoubleAttribute>(atof(cursor)));
353 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
354 evt.
add_attribute(
"alphaQCD", std::make_shared<DoubleAttribute>(atof(cursor)));
357 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
358 evt.
add_attribute(
"alphaQED", std::make_shared<DoubleAttribute>(atof(cursor)));
361 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
362 evt.
add_attribute(
"signal_process_id", std::make_shared<IntAttribute>(atoi(cursor)));
365 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
366 evt.
add_attribute(
"signal_process_vertex", std::make_shared<IntAttribute>(atoi(cursor)));
369 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
370 vertices_count = atoi(cursor);
373 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
376 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
379 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
380 random_states_size = atoi(cursor);
381 random_states.resize(random_states_size);
383 for (
int i = 0; i < random_states_size; ++i ) {
384 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
385 random_states[i] = atoi(cursor);
388 if (random_states.size()) evt.
add_attribute(
"random_states", std::make_shared<VectorLongIntAttribute>(random_states));
391 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
392 weights_size = atoi(cursor);
393 weights.resize(weights_size);
395 for (
int i = 0; i < weights_size; ++i ) {
396 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
397 weights[i] = atof(cursor);
402 HEPMC3_DEBUG(10,
"ReaderAsciiHepMC2: E: " << evt.
event_number() <<
" (" << vertices_count <<
"V, " << weights_size <<
"W, " << random_states_size <<
"RS)")
404 return vertices_count;
408 const char *cursor = buf;
411 if ( !(cursor = strchr(cursor+1,
' ')) )
return false;
416 if ( !(cursor = strchr(cursor+1,
' ')) )
return false;
420 evt.
set_units(momentum_unit, length_unit);
428 GenVertexPtr data = std::make_shared<GenVertex>();
429 GenVertexPtr data_ghost = std::make_shared<GenVertex>();
430 const char *cursor = buf;
432 int num_particles_out = 0;
433 int weights_size = 0;
434 std::vector<double> weights(0);
436 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
437 barcode = atoi(cursor);
440 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
441 data->set_status(atoi(cursor));
444 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
445 double X(atof(cursor));
448 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
449 double Y(atof(cursor));
452 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
453 double Z(atof(cursor));
456 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
457 double T(atof(cursor));
461 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
464 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
465 num_particles_out = atoi(cursor);
468 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
469 weights_size = atoi(cursor);
470 weights.resize(weights_size);
472 for (
int i = 0; i < weights_size; ++i ) {
473 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
474 weights[i] = atof(cursor);
485 if (weights.size()) data_ghost->add_attribute(
"weights", std::make_shared<VectorDoubleAttribute>(weights));
489 HEPMC3_DEBUG(10,
"ReaderAsciiHepMC2: V: " << -(
int)
m_vertex_cache.size() <<
" (old barcode " << barcode <<
") " << num_particles_out <<
" particles)")
491 return num_particles_out;
495 GenParticlePtr data = std::make_shared<GenParticle>();
496 GenParticlePtr data_ghost = std::make_shared<GenParticle>();
498 const char *cursor = buf;
502 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
505 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
506 data->set_pid(atoi(cursor));
509 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
510 double Px(atof(cursor));
513 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
514 double Py(atof(cursor));
517 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
518 double Pz(atof(cursor));
521 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
522 double E(atof(cursor));
526 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
527 data->set_generated_mass(atof(cursor));
530 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
531 data->set_status(atoi(cursor));
534 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
535 double theta_v = atof(cursor);
536 if (theta_v != 0.0) data_ghost->add_attribute(
"theta", std::make_shared<DoubleAttribute>(theta_v));
539 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
540 double phi_v = atof(cursor);
541 if (phi_v != 0.0) data_ghost->add_attribute(
"phi", std::make_shared<DoubleAttribute>(phi_v));
544 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
545 end_vtx = atoi(cursor);
548 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
549 int flowsize = atoi(cursor);
551 std::map<int, int> flows;
552 for (
int i = 0; i < flowsize; i++)
554 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
555 int flowindex = atoi(cursor);
556 if ( !(cursor = strchr(cursor+1,
' ')) )
return -1;
557 int flowvalue = atoi(cursor);
558 flows[flowindex] = flowvalue;
562 std::vector<int> vectorflows;
563 for (
auto f: flows) vectorflows.push_back(f.second);
564 data_ghost->add_attribute(
"flows", std::make_shared<VectorIntAttribute>(vectorflows));
585 const char *cursor = buf;
586 std::shared_ptr<GenCrossSection> xs = std::make_shared<GenCrossSection>();
588 if ( !(cursor = strchr(cursor+1,
' ')) )
return false;
589 double xs_val = atof(cursor);
591 if ( !(cursor = strchr(cursor+1,
' ')) )
return false;
592 double xs_err = atof(cursor);
594 xs->set_cross_section(xs_val, xs_err);
601 const char *cursor = buf;
602 const char *cursor2 = buf;
604 std::vector<std::string> w_names;
609 if ( !(cursor = strchr(cursor+1,
' ')) )
return false;
610 w_count = atoi(cursor);
612 if ( w_count <= 0 )
return false;
614 w_names.resize(w_count);
616 for (
int i=0; i < w_count; ++i ) {
618 if ( !(cursor = strchr(cursor+1,
'"')) )
return false;
619 if ( !(cursor2 = strchr(cursor+1,
'"')) )
return false;
624 w_names[i].assign(cursor, cursor2-cursor);
629 run_info()->set_weight_names(w_names);
635 std::shared_ptr<GenHeavyIon> hi = std::make_shared<GenHeavyIon>();
636 const char *cursor = buf;
638 if ( !(cursor = strchr(cursor+1,
' ')) )
return false;
639 hi->Ncoll_hard = atoi(cursor);
641 if ( !(cursor = strchr(cursor+1,
' ')) )
return false;
642 hi->Npart_proj = atoi(cursor);
644 if ( !(cursor = strchr(cursor+1,
' ')) )
return false;
645 hi->Npart_targ = atoi(cursor);
647 if ( !(cursor = strchr(cursor+1,
' ')) )
return false;
648 hi->Ncoll = atoi(cursor);
650 if ( !(cursor = strchr(cursor+1,
' ')) )
return false;
651 hi->spectator_neutrons = atoi(cursor);
653 if ( !(cursor = strchr(cursor+1,
' ')) )
return false;
654 hi->spectator_protons = atoi(cursor);
656 if ( !(cursor = strchr(cursor+1,
' ')) )
return false;
657 hi->N_Nwounded_collisions = atoi(cursor);
659 if ( !(cursor = strchr(cursor+1,
' ')) )
return false;
660 hi->Nwounded_N_collisions = atoi(cursor);
662 if ( !(cursor = strchr(cursor+1,
' ')) )
return false;
663 hi->Nwounded_Nwounded_collisions = atoi(cursor);
665 if ( !(cursor = strchr(cursor+1,
' ')) )
return false;
666 hi->impact_parameter = atof(cursor);
668 if ( !(cursor = strchr(cursor+1,
' ')) )
return false;
669 hi->event_plane_angle = atof(cursor);
671 if ( !(cursor = strchr(cursor+1,
' ')) )
return false;
672 hi->eccentricity = atof(cursor);
674 if ( !(cursor = strchr(cursor+1,
' ')) )
return false;
675 hi->sigma_inel_NN = atof(cursor);
678 hi->centrality = 0.0;
686 std::shared_ptr<GenPdfInfo> pi = std::make_shared<GenPdfInfo>();
687 const char *cursor = buf;
689 if ( !(cursor = strchr(cursor+1,
' ')) )
return false;
690 pi->parton_id[0] = atoi(cursor);
692 if ( !(cursor = strchr(cursor+1,
' ')) )
return false;
693 pi->parton_id[1] = atoi(cursor);
695 if ( !(cursor = strchr(cursor+1,
' ')) )
return false;
696 pi->x[0] = atof(cursor);
698 if ( !(cursor = strchr(cursor+1,
' ')) )
return false;
699 pi->x[1] = atof(cursor);
701 if ( !(cursor = strchr(cursor+1,
' ')) )
return false;
702 pi->scale = atof(cursor);
704 if ( !(cursor = strchr(cursor+1,
' ')) )
return false;
705 pi->xf[0] = atof(cursor);
707 if ( !(cursor = strchr(cursor+1,
' ')) )
return false;
708 pi->xf[1] = atof(cursor);
712 if ( !(cursor = strchr(cursor+1,
' ')) ) pdfids =
false;
713 if (pdfids) pi->pdf_id[0] = atoi(cursor);
714 else pi->pdf_id[0] = 0;
716 if (pdfids)
if ( !(cursor = strchr(cursor+1,
' ')) ) pdfids =
false;
717 if (pdfids) pi->pdf_id[1] = atoi(cursor);
718 else pi->pdf_id[1] = 0;
728 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.
Stores event-related information.
std::shared_ptr< T > attribute(const std::string &name, const int &id=0) const
Get attribute of type T.
void add_tree(const std::vector< GenParticlePtr > &particles)
Add whole tree in topological order.
void add_vertex(GenVertexPtr v)
Add vertex.
int event_number() const
Get event number.
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.
std::map< std::string, std::map< int, std::shared_ptr< Attribute > > > attributes() const
Get a copy of the list of attributes.
void add_attribute(const std::string &name, const std::shared_ptr< Attribute > &att, const int &id=0)
void remove_attribute(const std::string &name, const int &id=0)
Remove attribute.
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.
int value() const
get the value associated to this Attribute.
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 stream.
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.
std::vector< long int > value() const
get the value associated to this Attribute.