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);
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();
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);
297 const char *cursor = buf;
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);
306 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
307 event_no = atoi(cursor);
311 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
312 std::shared_ptr<IntAttribute> mpi = std::make_shared<IntAttribute>(atoi(cursor));
316 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
317 std::shared_ptr<DoubleAttribute> event_scale = std::make_shared<DoubleAttribute>(atof(cursor));
321 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
322 std::shared_ptr<DoubleAttribute> alphaQCD = std::make_shared<DoubleAttribute>(atof(cursor));
326 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
327 std::shared_ptr<DoubleAttribute> alphaQED = std::make_shared<DoubleAttribute>(atof(cursor));
331 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
332 std::shared_ptr<IntAttribute> signal_process_id = std::make_shared<IntAttribute>(atoi(cursor));
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);
341 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
342 vertices_count = atoi(cursor);
345 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
348 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
351 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
352 random_states_size = atoi(cursor);
353 random_states.resize(random_states_size);
355 for (
int i = 0; i < random_states_size; ++i ) {
356 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
357 random_states[i] = atoi(cursor);
361 evt.
add_attribute(
"random_states",std::make_shared<VectorLongIntAttribute>(random_states));
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]));
369 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
370 weights_size = atoi(cursor);
371 weights.resize(weights_size);
373 for (
int i = 0; i < weights_size; ++i ) {
374 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
375 weights[i] = atof(cursor);
380 HEPMC3_DEBUG( 10,
"ReaderAsciiHepMC2: E: "<<event_no<<
" ("<<vertices_count<<
"V, "<<weights_size<<
"W, "<<random_states_size<<
"RS)" )
382 return vertices_count;
386 const char *cursor = buf;
389 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
394 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
398 evt.
set_units(momentum_unit,length_unit);
406 GenVertexPtr data = std::make_shared<GenVertex>();
407 GenVertexPtr data_ghost = std::make_shared<GenVertex>();
409 const char *cursor = buf;
411 int num_particles_out = 0;
412 int weights_size = 0;
413 std::vector<double> weights(0);
415 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
416 barcode = atoi(cursor);
419 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
420 data->set_status( atoi(cursor) );
423 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
424 position.
setX(atof(cursor));
427 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
428 position.
setY(atof(cursor));
431 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
432 position.
setZ(atof(cursor));
435 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
436 position.
setT(atof(cursor));
437 data->set_position( position );
440 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
443 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
444 num_particles_out = atoi(cursor);
448 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
449 weights_size = atoi(cursor);
450 weights.resize(weights_size);
452 for (
int i = 0; i < weights_size; ++i ) {
453 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
454 weights[i] = atof(cursor);
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]));
471 data_ghost->add_attribute(
"weights",std::make_shared<VectorDoubleAttribute>(weights));
472 data_ghost->add_attribute(
"weights",std::make_shared<VectorDoubleAttribute>(weights));
476 HEPMC3_DEBUG( 10,
"ReaderAsciiHepMC2: V: "<<-(
int)
m_vertex_cache.size()<<
" (old barcode"<<barcode<<
") "<<num_particles_out<<
" particles)" )
478 return num_particles_out;
482 GenParticlePtr data = std::make_shared<GenParticle>();
483 GenParticlePtr data_ghost = std::make_shared<GenParticle>();
486 const char *cursor = buf;
490 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
493 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
494 data->set_pid( atoi(cursor) );
497 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
498 momentum.
setPx(atof(cursor));
501 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
502 momentum.
setPy(atof(cursor));
505 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
506 momentum.
setPz(atof(cursor));
509 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
510 momentum.
setE(atof(cursor));
511 data->set_momentum(momentum);
514 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
515 data->set_generated_mass( atof(cursor) );
518 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
519 data->set_status( atoi(cursor) );
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);
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);
532 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
533 end_vtx = atoi(cursor);
536 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
537 int flowsize=atoi(cursor);
539 std::map<int,int> flows;
540 for (
int i=0; i<flowsize; i++)
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;
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));
556 for (
auto f: flows) data_ghost->add_attribute(
"flow"+std::to_string((
long long int)f.first),std::make_shared<IntAttribute>(f.second));
577 const char *cursor = buf;
578 std::shared_ptr<GenCrossSection> xs = std::make_shared<GenCrossSection>();
580 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
581 double xs_val = atof(cursor);
583 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
584 double xs_err = atof(cursor);
586 xs->set_cross_section( xs_val, xs_err);
593 const char *cursor = buf;
594 const char *cursor2 = buf;
596 std::vector<std::string> w_names;
601 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
602 w_count = atoi(cursor);
604 if( w_count <= 0 )
return false;
606 w_names.resize(w_count);
608 for(
int i=0; i < w_count; ++i ) {
610 if( !(cursor = strchr(cursor+1,
'"')) )
return false;
611 if( !(cursor2 = strchr(cursor+1,
'"')) )
return false;
616 w_names[i].assign(cursor, cursor2-cursor);
621 run_info()->set_weight_names(w_names);
627 std::shared_ptr<GenHeavyIon> hi = std::make_shared<GenHeavyIon>();
628 const char *cursor = buf;
630 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
631 hi->Ncoll_hard = atoi(cursor);
633 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
634 hi->Npart_proj = atoi(cursor);
636 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
637 hi->Npart_targ = atoi(cursor);
639 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
640 hi->Ncoll = atoi(cursor);
642 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
643 hi->spectator_neutrons = atoi(cursor);
645 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
646 hi->spectator_protons = atoi(cursor);
648 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
649 hi->N_Nwounded_collisions = atoi(cursor);
651 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
652 hi->Nwounded_N_collisions = atoi(cursor);
654 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
655 hi->Nwounded_Nwounded_collisions = atoi(cursor);
657 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
658 hi->impact_parameter = atof(cursor);
660 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
661 hi->event_plane_angle = atof(cursor);
663 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
664 hi->eccentricity = atof(cursor);
666 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
667 hi->sigma_inel_NN = atof(cursor);
670 hi->centrality = 0.0;
678 std::shared_ptr<GenPdfInfo> pi = std::make_shared<GenPdfInfo>();
679 const char *cursor = buf;
681 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
682 pi->parton_id[0] = atoi(cursor);
684 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
685 pi->parton_id[1] = atoi(cursor);
687 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
688 pi->x[0] = atof(cursor);
690 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
691 pi->x[1] = atof(cursor);
693 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
694 pi->scale = atof(cursor);
696 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
697 pi->xf[0] = atof(cursor);
699 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
700 pi->xf[1] = atof(cursor);
704 if( !(cursor = strchr(cursor+1,
' ')) ) pdfids=
false;
705 if(pdfids) pi->pdf_id[0] = atoi(cursor);
706 else pi->pdf_id[0] =0;
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;
720 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.
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.
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.
void add_attribute(const std::string &name, const std::shared_ptr< Attribute > &att, const int &id=0)
Add event attribute to event.
const Units::LengthUnit & length_unit() const
Get length unit.
void clear()
Remove contents of this event.
const std::vector< double > & weights() const
Get event weight values as a vector.
const Units::MomentumUnit & momentum_unit() const
Get momentum unit.
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 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::map< std::string, std::string > m_options
options
void set_run_info(std::shared_ptr< GenRunInfo > run)
Set the global GenRunInfo object.
std::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.
Attribute that holds a vector of integers of type int.