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;
int event_number() const
Get event number.
std::ifstream m_file
Input file.
const Units::MomentumUnit & momentum_unit() const
Get momentum unit.
std::map< std::string, std::string > m_options
options
const Units::LengthUnit & length_unit() const
Get length unit.
#define HEPMC3_WARNING(MESSAGE)
Macro for printing HEPMC3_HEPMC3_WARNING messages.
void add_vertex(GenVertexPtr v)
Add vertex.
Attribute that holds a vector of FPs of type double.
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 of class GenParticle.
std::vector< GenParticlePtr > m_particle_cache_ghost
Particle cache for attributes.
std::vector< double > value() const
get the value associated to this Attribute.
Definition of attribute class GenHeavyIon.
int parse_particle_information(const char *buf)
Parse particle.
~ReaderAsciiHepMC2()
Destructor.
Definition of class GenVertex.
bool failed() override
Return status of the stream.
Attribute that holds a vector of integers of type int.
#define HEPMC3_DEBUG(LEVEL, MESSAGE)
Macro for printing debug messages with appropriate debug level.
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.
int value() const
get the value associated to this Attribute.
static std::string name(MomentumUnit u)
Get name of momentum unit.
void add_attribute(const std::string &name, const std::shared_ptr< Attribute > &att, const int &id=0)
void close() override
Close file stream.
ReaderAsciiHepMC2(const std::string &filename)
Default constructor.
const std::vector< double > & weights() const
Get event weight values as a vector.
std::shared_ptr< T > attribute(const std::string &name, const int &id=0) const
Get attribute of type T.
LengthUnit
Position units.
Definition of class ReaderAsciiHepMC2.
MomentumUnit
Momentum units.
Stores event-related information.
bool parse_weight_names(const char *buf)
Parse weight names.
std::vector< int > m_end_vertex_barcodes
Old end vertex barcodes.
Definition of class Setup.
void add_tree(const std::vector< GenParticlePtr > &particles)
Add whole tree in topological order.
std::shared_ptr< GenRunInfo > run_info() const
Get the global GenRunInfo object.
void remove_attribute(const std::string &name, const int &id=0)
Remove attribute.
static MomentumUnit momentum_unit(const std::string &name)
Get momentum unit based on its name.
void set_run_info(std::shared_ptr< GenRunInfo > run)
Set the GenRunInfo object by smart pointer.
void reserve(const size_t &particles, const size_t &vertices=0)
Reserve memory for particles and vertices.
std::istream * m_stream
For ctor when reading from stream.
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 of event attribute class GenPdfInfo.
void set_run_info(std::shared_ptr< GenRunInfo > run)
Set the global GenRunInfo object.
#define HEPMC3_ERROR(MESSAGE)
Macro for printing error messages.
Attribute that holds a vector of integers of type int.
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< long int > value() const
get the value associated to this Attribute.
std::map< std::string, std::map< int, std::shared_ptr< Attribute > > > attributes() const
Get a copy of the list of attributes.
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.
std::vector< int > value() const
get the value associated to this Attribute.
void clear()
Remove contents of this event.
bool parse_units(GenEvent &evt, const char *buf)
Parse units.
bool parse_heavy_ion(GenEvent &evt, const char *buf)
Parse heavy ion information.
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.
bool skip(const int) override
skip events
void set_event_number(const int &num)
Set event number.