11 #include "HepMC/ReaderAsciiHepMC2.h" 13 #include "HepMC/GenEvent.h" 14 #include "HepMC/GenVertex.h" 15 #include "HepMC/GenParticle.h" 16 #include "HepMC/GenHeavyIon.h" 17 #include "HepMC/GenPdfInfo.h" 18 #include "HepMC/Setup.h" 28 ERROR(
"ReaderAsciiHepMC2: could not open input file: "<<filename )
35 bool parsed_event_header =
false;
36 bool is_parsing_successful =
true;
37 int parsing_result = 0;
38 unsigned int vertices_count = 0;
39 unsigned int current_vertex_particles_count = 0;
40 unsigned int current_vertex_particles_parsed= 0;
58 if( strlen(buf) == 0 )
continue;
61 if( strncmp(buf,
"HepMC",5) == 0 ) {
62 if(parsed_event_header) {
63 is_parsing_successful =
true;
72 if(parsing_result<0) {
73 is_parsing_successful =
false;
74 ERROR(
"ReaderAsciiHepMC2: error parsing event information" )
77 vertices_count = parsing_result;
82 is_parsing_successful =
true;
84 parsed_event_header =
true;
91 if(current_vertex_particles_parsed < current_vertex_particles_count) {
92 is_parsing_successful =
false;
95 current_vertex_particles_parsed = 0;
99 if(parsing_result<0) {
100 is_parsing_successful =
false;
101 ERROR(
"ReaderAsciiHepMC2: error parsing vertex information" )
104 current_vertex_particles_count = parsing_result;
105 is_parsing_successful =
true;
112 if(parsing_result<0) {
113 is_parsing_successful =
false;
114 ERROR(
"ReaderAsciiHepMC2: error parsing particle information" )
117 ++current_vertex_particles_parsed;
118 is_parsing_successful =
true;
137 WARNING(
"ReaderAsciiHepMC2: skipping unrecognised prefix: " << buf[0] )
138 is_parsing_successful =
true;
142 if( !is_parsing_successful )
break;
146 if( parsed_event_header && buf[0]==
'E' )
break;
153 if( is_parsing_successful && current_vertex_particles_parsed < current_vertex_particles_count ) {
154 ERROR(
"ReaderAsciiHepMC2: not all particles parsed" )
155 is_parsing_successful =
false;
158 else if( is_parsing_successful &&
m_vertex_cache.size() != vertices_count ) {
159 ERROR(
"ReaderAsciiHepMC2: not all vertices parsed" )
160 is_parsing_successful =
false;
163 if( !is_parsing_successful ) {
164 ERROR(
"ReaderAsciiHepMC2: event parsing failed. Returning empty event" )
165 DEBUG( 1,
"Parsing failed at line:" << std::endl << buf )
168 m_file.clear(std::ios::badbit);
205 const char *cursor = buf;
207 int vertices_count = 0;
208 int random_states_size = 0;
209 int weights_size = 0;
210 std::vector<long> random_states(0);
211 std::vector<double> weights(0);
214 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
215 event_no = atoi(cursor);
219 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
222 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
225 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
228 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
231 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
234 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
237 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
238 vertices_count = atoi(cursor);
241 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
244 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
247 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
248 random_states_size = atoi(cursor);
249 random_states.resize(random_states_size);
251 for (
int i = 0; i < random_states_size; ++i ) {
252 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
253 random_states[i] = atoi(cursor);
257 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
258 weights_size = atoi(cursor);
259 weights.resize(weights_size);
261 for (
int i = 0; i < weights_size; ++i ) {
262 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
263 weights[i] = atof(cursor);
268 DEBUG( 10,
"ReaderAsciiHepMC2: E: "<<event_no<<
" ("<<vertices_count<<
"V, "<<weights_size<<
"W, "<<random_states_size<<
"RS)" )
270 return vertices_count;
274 const char *cursor = buf;
277 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
282 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
286 evt.
set_units(momentum_unit,length_unit);
296 const char *cursor = buf;
298 int num_particles_out = 0;
301 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
302 barcode = atoi(cursor);
305 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
308 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
309 position.setX(atof(cursor));
312 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
313 position.setY(atof(cursor));
316 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
317 position.setZ(atof(cursor));
320 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
321 position.setT(atof(cursor));
322 data->set_position( position );
325 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
328 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
329 num_particles_out = atoi(cursor);
337 DEBUG( 10,
"ReaderAsciiHepMC2: V: "<<-(
int)
m_vertex_cache.size()<<
" (old barcode"<<barcode<<
") "<<num_particles_out<<
" particles)" )
339 return num_particles_out;
345 const char *cursor = buf;
349 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
352 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
353 data->set_pid( atoi(cursor) );
356 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
357 momentum.setPx(atof(cursor));
360 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
361 momentum.setPy(atof(cursor));
364 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
365 momentum.setPz(atof(cursor));
368 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
369 momentum.setE(atof(cursor));
370 data->set_momentum(momentum);
373 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
374 data->set_generated_mass( atof(cursor) );
377 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
378 data->set_status( atoi(cursor) );
381 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
384 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
387 if( !(cursor = strchr(cursor+1,
' ')) )
return -1;
388 end_vtx = atoi(cursor);
404 DEBUG( 10,
"ReaderAsciiHepMC2: P: "<<
m_particle_cache.size()<<
" ( pid: "<<data->pid()<<
") end vertex: "<<end_vtx )
410 const char *cursor = buf;
411 shared_ptr<GenCrossSection> xs = make_shared<GenCrossSection>();
413 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
414 double xs_val = atof(cursor);
416 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
417 double xs_err = atof(cursor);
419 xs->set_cross_section( xs_val , xs_err);
426 const char *cursor = buf;
427 const char *cursor2 = buf;
429 vector<string> w_names;
434 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
435 w_count = atoi(cursor);
437 if( w_count <= 0 )
return false;
439 w_names.resize(w_count);
441 for(
int i=0; i < w_count; ++i ) {
443 if( !(cursor = strchr(cursor+1,
'"')) )
return false;
444 if( !(cursor2 = strchr(cursor+1,
'"')) )
return false;
449 w_names[i].assign(cursor, cursor2-cursor);
454 run_info()->set_weight_names(w_names);
460 shared_ptr<GenHeavyIon> hi = make_shared<GenHeavyIon>();
461 const char *cursor = buf;
463 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
464 hi->Ncoll_hard = atoi(cursor);
466 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
467 hi->Npart_proj = atoi(cursor);
469 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
470 hi->Npart_targ = atoi(cursor);
472 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
473 hi->Ncoll = atoi(cursor);
475 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
476 hi->spectator_neutrons = atoi(cursor);
478 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
479 hi->spectator_protons = atoi(cursor);
481 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
482 hi->N_Nwounded_collisions = atoi(cursor);
484 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
485 hi->Nwounded_N_collisions = atoi(cursor);
487 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
488 hi->Nwounded_Nwounded_collisions = atoi(cursor);
490 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
491 hi->impact_parameter = atof(cursor);
493 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
494 hi->event_plane_angle = atof(cursor);
496 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
497 hi->eccentricity = atof(cursor);
499 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
500 hi->sigma_inel_NN = atof(cursor);
503 hi->centrality = 0.0;
511 shared_ptr<GenPdfInfo> pi = make_shared<GenPdfInfo>();
512 const char *cursor = buf;
514 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
515 pi->parton_id[0] = atoi(cursor);
517 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
518 pi->parton_id[1] = atoi(cursor);
520 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
521 pi->x[0] = atof(cursor);
523 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
524 pi->x[1] = atof(cursor);
526 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
527 pi->scale = atof(cursor);
529 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
530 pi->xf[0] = atof(cursor);
532 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
533 pi->xf[1] = atof(cursor);
535 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
536 pi->pdf_id[0] = atoi(cursor);
538 if( !(cursor = strchr(cursor+1,
' ')) )
return false;
539 pi->pdf_id[1] = atoi(cursor);
547 if( !
m_file.is_open() )
return;
bool parse_heavy_ion(GenEvent &evt, const char *buf)
Parse heavy ion information.
const Units::MomentumUnit & momentum_unit() const
Get momentum unit.
ifstream m_file
Input file.
void set_event_number(int num)
Set event number.
static std::string name(MomentumUnit u)
Get name of momentum unit.
void set_units(Units::MomentumUnit new_momentum_unit, Units::LengthUnit new_length_unit)
Change event units Converts event from current units to new ones.
vector< GenParticlePtr > m_particle_cache
Particle cache.
void add_tree(const std::vector< GenParticlePtr > &particles)
Add whole tree in topological order.
ReaderAsciiHepMC2(const std::string &filename)
Default constructor.
bool read_event(GenEvent &evt)
Implementation of Reader::read_event.
bool parse_weight_names(const char *buf)
Parse weight names.
vector< int > m_vertex_barcodes
Old vertex barcodes.
bool parse_units(GenEvent &evt, const char *buf)
Parse units.
void close()
Close file stream.
vector< GenVertexPtr > m_vertex_cache
Vertex cache.
shared_ptr< GenRunInfo > run_info() const
Get the global GenRunInfo object.
void clear()
Remove contents of this event.
void set_run_info(shared_ptr< GenRunInfo > run)
Set the global GenRunInfo object.
Stores event-related information.
bool failed()
Return status of the stream.
static MomentumUnit momentum_unit(const std::string &name)
Get momentum unit based on its name.
bool parse_xs_info(GenEvent &evt, const char *buf)
Parse pdf information.
int parse_vertex_information(const char *buf)
Parse vertex.
void add_attribute(const string &name, const shared_ptr< Attribute > &att, int id=0)
Add event attribute to event.
MomentumUnit
Momentum units.
int parse_event_information(GenEvent &evt, const char *buf)
Parse event.
void reserve(unsigned int particles, unsigned int vertices=0)
Reserve memory for particles and vertices.
Definition of template class SmartPointer.
bool parse_pdf_info(GenEvent &evt, const char *buf)
Parse pdf information.
vector< int > m_end_vertex_barcodes
Old end vertex barcodes.
const Units::LengthUnit & length_unit() const
Get length unit.
int parse_particle_information(const char *buf)
Parse particle.
static LengthUnit length_unit(const std::string &name)
Get length unit based on its name.
void set_run_info(shared_ptr< GenRunInfo > run)
Set the GenRunInfo object by smart pointer.
const std::vector< double > & weights() const
Get event weight values as a vector.
LengthUnit
Position units.