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)
37 if( !m_stream->good() ) {
38 ERROR(
"ReaderAsciiHepMC2: could not open input stream " )
48 if ( (!
m_file.is_open()) && (!m_isstream) )
return false;
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;
75 m_isstream ? m_stream->getline(buf,max_buffer_size) :
m_file.getline(buf,max_buffer_size);
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;
83 m_isstream ? m_stream->clear(ios::eofbit) :
m_file.clear(ios::eofbit);
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;
167 m_isstream ? peek = m_stream->peek() : peek =
m_file.peek();
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 )
188 m_isstream ? m_stream->clear(ios::badbit) :
m_file.clear(ios::badbit);
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);
else pi->pdf_id[0] =0;
638 if(pdfids)
if( !(cursor = strchr(cursor+1,
' ')) ) pdfids=
false;
639 if(pdfids) pi->pdf_id[1] = atoi(cursor);
else pi->pdf_id[1] =0;
648 if( !
m_file.is_open() )
return;
void set_run_info(shared_ptr< GenRunInfo > run)
Set the GenRunInfo object by smart pointer.
std::ifstream m_file
Input file.
const Units::MomentumUnit & momentum_unit() const
Get momentum unit.
vector< GenVertexPtr > m_vertex_cache
Vertex cache.
const Units::LengthUnit & length_unit() const
Get length unit.
void add_vertex(GenVertexPtr v)
Add vertex.
static LengthUnit length_unit(const std::string &name)
Get length unit based on its name.
#define ERROR(MESSAGE)
Macro for printing error messages.
Definition of class GenParticle.
bool failed()
Return status of the stream.
shared_ptr< GenRunInfo > run_info() const
Get the global GenRunInfo object.
Definition of attribute class GenHeavyIon.
int parse_particle_information(const char *buf)
Parse particle.
~ReaderAsciiHepMC2()
Destructor.
Definition of class GenVertex.
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.
static std::string name(MomentumUnit u)
Get name of momentum unit.
void setPy(double pyy)
Set y-component of momentum.
ReaderAsciiHepMC2(const std::string &filename)
Default constructor.
#define DEBUG(LEVEL, MESSAGE)
Macro for printing debug messages with appropriate debug level.
void close()
Close file stream.
const std::vector< double > & weights() const
Get event weight values as a vector.
LengthUnit
Position units.
Definition of class ReaderAsciiHepMC2.
MomentumUnit
Momentum units.
Stores event-related information.
Attribute that holds a real number as a double.
bool parse_weight_names(const char *buf)
Parse weight names.
bool read_event(GenEvent &evt)
Implementation of Reader::read_event.
void setPz(double pzz)
Set z-component of momentum.
vector< int > m_vertex_barcodes
Old vertex barcodes.
vector< GenParticlePtr > m_particle_cache
Particle cache.
vector< GenVertexPtr > m_vertex_cache_ghost
Vertex cache for attributes.
Definition of class Setup.
void add_tree(const std::vector< GenParticlePtr > &particles)
Add whole tree in topological order.
#define WARNING(MESSAGE)
Macro for printing warning messages.
static MomentumUnit momentum_unit(const std::string &name)
Get momentum unit based on its name.
void reserve(const size_t &particles, const size_t &vertices=0)
Reserve memory for particles and vertices.
void setPx(double pxx)
Set x-component of momentum.
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_run_info(shared_ptr< GenRunInfo > run)
Set the global GenRunInfo object.
Definition of event attribute class GenPdfInfo.
vector< GenParticlePtr > m_particle_cache_ghost
Particle cache for attributes.
Definition of class GenEvent.
int parse_event_information(GenEvent &evt, const char *buf)
Parse event.
void clear()
Remove contents of this event.
vector< int > m_end_vertex_barcodes
Old end vertex barcodes.
bool parse_units(GenEvent &evt, const char *buf)
Parse units.
bool parse_heavy_ion(GenEvent &evt, const char *buf)
Parse heavy ion information.
void setE(double ee)
Set energy component of momentum.
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.
void set_event_number(const int &num)
Set event number.