HepMC3 event record library
ReaderAsciiHepMC2.cc
Go to the documentation of this file.
1// -*- C++ -*-
2//
3// This file is part of HepMC
4// Copyright (C) 2014-2019 The HepMC collaboration (see AUTHORS for details)
5//
6/**
7 * @file ReaderAsciiHepMC2.cc
8 * @brief Implementation of \b class ReaderAsciiHepMC2
9 *
10 */
12
13#include "HepMC3/GenEvent.h"
14#include "HepMC3/GenVertex.h"
15#include "HepMC3/GenParticle.h"
16#include "HepMC3/GenHeavyIon.h"
17#include "HepMC3/GenPdfInfo.h"
18#include "HepMC3/Setup.h"
19#include <cstring>
20#include <cstdlib>
21namespace HepMC3 {
22
23ReaderAsciiHepMC2::ReaderAsciiHepMC2(const std::string& filename):
24 m_file(filename), m_stream(0), m_isstream(false) {
25 if( !m_file.is_open() ) {
26 HEPMC3_ERROR( "ReaderAsciiHepMC2: could not open input file: "<<filename )
27 }
28 set_run_info(std::make_shared<GenRunInfo>());
30}
31// Ctor for reading from stdin
33 : m_stream(&stream), m_isstream(true)
34{
35 if( !m_stream->good() ) {
36 HEPMC3_ERROR( "ReaderAsciiHepMC2: could not open input stream " )
37 }
38 set_run_info(std::make_shared<GenRunInfo>());
40}
41
43
44bool ReaderAsciiHepMC2::skip(const int n)
45{
46 const size_t max_buffer_size=512*512;
47 char buf[max_buffer_size];
48 int nn=n;
49 char peek;
50 while(!failed()) {
51 if ( (!m_file.is_open()) && (!m_isstream) ) return false;
52 m_isstream ? peek = m_stream->peek() : peek = m_file.peek();
53 if( peek=='E' ) nn--;
54 if (nn<0) return true;
55 m_isstream ? m_stream->getline(buf,max_buffer_size) : m_file.getline(buf,max_buffer_size);
56 }
57 return true;
58}
59
61 if ( (!m_file.is_open()) && (!m_isstream) ) return false;
62
63 char peek;
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;
73
74 evt.clear();
76
77 // Empty cache
78 m_vertex_cache.clear();
79 m_vertex_barcodes.clear();
80
81 m_particle_cache.clear();
84 //
85 // Parse event, vertex and particle information
86 //
87 while(!failed()) {
88 m_isstream ? m_stream->getline(buf,max_buffer_size) : m_file.getline(buf,max_buffer_size);
89 if( strlen(buf) == 0 ) continue;
90 // Check for IO_GenEvent header/footer
91 if( strncmp(buf,"HepMC",5) == 0 ) {
92 if( strncmp(buf,"HepMC::Version",14) != 0 && strncmp(buf,"HepMC::IO_GenEvent",18)!=0 )
93 {
94 HEPMC3_WARNING( "ReaderAsciiHepMC2: found unsupported expression in header. Will close the input." )
95 std::cout<<buf<<std::endl;
96 m_isstream ? m_stream->clear(std::ios::eofbit) : m_file.clear(std::ios::eofbit);
97 }
98 if(parsed_event_header) {
99 is_parsing_successful = true;
100 break;
101 }
102 continue;
103 }
104 switch(buf[0]) {
105 case 'E':
106 parsing_result = parse_event_information(evt,buf);
107 if(parsing_result<0) {
108 is_parsing_successful = false;
109 HEPMC3_ERROR( "ReaderAsciiHepMC2: HEPMC3_ERROR parsing event information" )
110 }
111 else {
112 vertices_count = parsing_result;
113 m_vertex_cache.reserve(vertices_count);
114 m_particle_cache.reserve(vertices_count*3);
115 m_vertex_barcodes.reserve(vertices_count);
116 m_end_vertex_barcodes.reserve(vertices_count*3);
117 is_parsing_successful = true;
118 }
119 parsed_event_header = true;
120 break;
121 case 'V':
122 // If starting new vertex: verify if previous was fully parsed
123
124 /** @bug HepMC2 files produced with Pythia8 are known to have wrong
125 information about number of particles in vertex. Hence '<' sign */
126 if(current_vertex_particles_parsed < current_vertex_particles_count) {
127 is_parsing_successful = false;
128 break;
129 }
130 current_vertex_particles_parsed = 0;
131
132 parsing_result = parse_vertex_information(buf);
133
134 if(parsing_result<0) {
135 is_parsing_successful = false;
136 HEPMC3_ERROR( "ReaderAsciiHepMC2: HEPMC3_ERROR parsing vertex information" )
137 }
138 else {
139 current_vertex_particles_count = parsing_result;
140 is_parsing_successful = true;
141 }
142 break;
143 case 'P':
144
145 parsing_result = parse_particle_information(buf);
146
147 if(parsing_result<0) {
148 is_parsing_successful = false;
149 HEPMC3_ERROR( "ReaderAsciiHepMC2: HEPMC3_ERROR parsing particle information" )
150 }
151 else {
152 ++current_vertex_particles_parsed;
153 is_parsing_successful = true;
154 }
155 break;
156 case 'U':
157 is_parsing_successful = parse_units(evt,buf);
158 break;
159 case 'F':
160 is_parsing_successful = parse_pdf_info(evt,buf);
161 break;
162 case 'H':
163 is_parsing_successful = parse_heavy_ion(evt,buf);
164 break;
165 case 'N':
166 is_parsing_successful = parse_weight_names(buf);
167 break;
168 case 'C':
169 is_parsing_successful = parse_xs_info(evt,buf);
170 break;
171 default:
172 HEPMC3_WARNING( "ReaderAsciiHepMC2: skipping unrecognised prefix: " << buf[0] )
173 is_parsing_successful = true;
174 break;
175 }
176
177 if( !is_parsing_successful ) break;
178
179 // Check for next event
180 m_isstream ? peek = m_stream->peek() : peek = m_file.peek();
181 if( parsed_event_header && peek=='E' ) break;
182 }
183
184 // Check if all particles in last vertex were parsed
185 /** @bug HepMC2 files produced with Pythia8 are known to have wrong
186 information about number of particles in vertex. Hence '<' sign */
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;
190 }
191 // Check if all vertices were parsed
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;
195 }
196
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 )
200 evt.clear();
201 m_isstream ? m_stream->clear(std::ios::badbit) : m_file.clear(std::ios::badbit);
202 return 0;
203 }
204
205 // Restore production vertex pointers
206 for(unsigned int i=0; i<m_particle_cache.size(); ++i) {
207 if( !m_end_vertex_barcodes[i] ) continue;
208
209 for(unsigned int j=0; j<m_vertex_cache.size(); ++j) {
211 m_vertex_cache[j]->add_particle_in(m_particle_cache[i]);
212 break;
213 }
214 }
215 }
216
217 // Remove vertices with no incoming particles or no outgoing particles
218 for(unsigned int i=0; i<m_vertex_cache.size(); ++i) {
219 if( m_vertex_cache[i]->particles_in().size() == 0 ) {
220 HEPMC3_DEBUG( 30, "ReaderAsciiHepMC2::read_event - found a vertex without incoming particles: "<<m_vertex_cache[i]->id() );
221//Sometimes the root vertex has no incoming particles. Here we try to save the event.
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);
224 for (auto p: beams)
225 {
226 m_vertex_cache[i]->add_particle_in(p);
227 m_vertex_cache[i]->remove_particle_out(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() );
229 }
230 if (beams.size()==0) m_vertex_cache[i] = nullptr;
231 }
232 else if( m_vertex_cache[i]->particles_out().size() == 0 ) {
233 HEPMC3_DEBUG( 30, "ReaderAsciiHepMC2::read_event - found a vertex without outgouing particles: "<<m_vertex_cache[i]->id() );
234 m_vertex_cache[i] = nullptr;
235 }
236 }
237
238 // Reserve memory for the event
239 evt.reserve( m_particle_cache.size(), m_vertex_cache.size() );
240
241 // Add whole event tree in topological order
243
244 for(unsigned int i=0; i<m_particle_cache.size(); ++i) {
245 if(m_particle_cache_ghost[i]->attribute_names().size())
246 {
247 std::shared_ptr<DoubleAttribute> phi = m_particle_cache_ghost[i]->attribute<DoubleAttribute>("phi");
248 if (phi) m_particle_cache[i]->add_attribute("phi",phi);
249 std::shared_ptr<DoubleAttribute> theta = m_particle_cache_ghost[i]->attribute<DoubleAttribute>("theta");
250 if (theta) m_particle_cache[i]->add_attribute("theta",theta);
251 std::shared_ptr<IntAttribute> flow1 = m_particle_cache_ghost[i]->attribute<IntAttribute>("flow1");
252 if (m_options.find("particle_flows_are_separated")!=m_options.end())
253 {
254 if (flow1) m_particle_cache[i]->add_attribute("flow1",flow1);
255 std::shared_ptr<IntAttribute> flow2 = m_particle_cache_ghost[i]->attribute<IntAttribute>("flow2");
256 if (flow2) m_particle_cache[i]->add_attribute("flow2",flow2);
257 std::shared_ptr<IntAttribute> flow3 = m_particle_cache_ghost[i]->attribute<IntAttribute>("flow3");
258 if (flow3) m_particle_cache[i]->add_attribute("flow3",flow3);
259 }
260 else
261 {
262 std::shared_ptr<VectorIntAttribute> flows = m_particle_cache_ghost[i]->attribute<VectorIntAttribute>("flows");
263 if (flows) m_particle_cache[i]->add_attribute("flows",flows);
264 }
265 }
266 }
267
268 for(unsigned int i=0; i<m_vertex_cache.size(); ++i)
269 if(m_vertex_cache_ghost[i]->attribute_names().size())
270 {
271 for (size_t ii=0; ii<max_weights_size; ii++)
272 {
273 std::shared_ptr<DoubleAttribute> rs=m_vertex_cache_ghost[i]->attribute<DoubleAttribute>("weight"+std::to_string((long long unsigned int)ii));
274 if (!rs) break;
275 m_vertex_cache[i]->add_attribute("weight"+std::to_string((long long unsigned int)ii),rs);
276 }
277 }
279 m_vertex_cache_ghost.clear();
281 return 1;
282}
283
285 const char *cursor = buf;
286 int event_no = 0;
287 int vertices_count = 0;
288 int random_states_size = 0;
289 int weights_size = 0;
290 std::vector<long> random_states(0);
291 std::vector<double> weights(0);
292
293 // event number
294 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
295 event_no = atoi(cursor);
296 evt.set_event_number(event_no);
297
298 //mpi
299 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
300 std::shared_ptr<IntAttribute> mpi = std::make_shared<IntAttribute>(atoi(cursor));
301 evt.add_attribute("mpi",mpi);
302
303 //event scale
304 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
305 std::shared_ptr<DoubleAttribute> event_scale = std::make_shared<DoubleAttribute>(atof(cursor));
306 evt.add_attribute("event_scale",event_scale);
307
308 //alpha_qcd
309 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
310 std::shared_ptr<DoubleAttribute> alphaQCD = std::make_shared<DoubleAttribute>(atof(cursor));
311 evt.add_attribute("alphaQCD",alphaQCD);
312
313 //alpha_qed
314 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
315 std::shared_ptr<DoubleAttribute> alphaQED = std::make_shared<DoubleAttribute>(atof(cursor));
316 evt.add_attribute("alphaQED",alphaQED);
317
318 //signal_process_id
319 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
320 std::shared_ptr<IntAttribute> signal_process_id = std::make_shared<IntAttribute>(atoi(cursor));
321 evt.add_attribute("signal_process_id",signal_process_id);
322
323 //signal_process_vertex
324 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
325 std::shared_ptr<IntAttribute> signal_process_vertex = std::make_shared<IntAttribute>(atoi(cursor));
326 evt.add_attribute("signal_process_vertex",signal_process_vertex);
327
328 // num_vertices
329 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
330 vertices_count = atoi(cursor);
331
332 // SKIPPED: beam 1
333 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
334
335 // SKIPPED: beam 2
336 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
337
338 //random states
339 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
340 random_states_size = atoi(cursor);
341 random_states.resize(random_states_size);
342
343 for ( int i = 0; i < random_states_size; ++i ) {
344 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
345 random_states[i] = atoi(cursor);
346 }
347 if (m_options.find("event_random_states_are_separated")!=m_options.end())
348 {
349 evt.add_attribute("random_states",std::make_shared<VectorLongIntAttribute>(random_states));
350 }
351 else
352 {
353 for ( int i = 0; i < random_states_size; ++i )
354 evt.add_attribute("random_states"+std::to_string((long long unsigned int)i),std::make_shared<IntAttribute>(random_states[i]));
355 }
356 // weights
357 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
358 weights_size = atoi(cursor);
359 weights.resize(weights_size);
360
361 for ( int i = 0; i < weights_size; ++i ) {
362 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
363 weights[i] = atof(cursor);
364 }
365
366 evt.weights() = weights;
367
368 HEPMC3_DEBUG( 10, "ReaderAsciiHepMC2: E: "<<event_no<<" ("<<vertices_count<<"V, "<<weights_size<<"W, "<<random_states_size<<"RS)" )
369
370 return vertices_count;
371}
372
373bool ReaderAsciiHepMC2::parse_units(GenEvent &evt, const char *buf) {
374 const char *cursor = buf;
375
376 // momentum
377 if( !(cursor = strchr(cursor+1,' ')) ) return false;
378 ++cursor;
379 Units::MomentumUnit momentum_unit = Units::momentum_unit(cursor);
380
381 // length
382 if( !(cursor = strchr(cursor+1,' ')) ) return false;
383 ++cursor;
384 Units::LengthUnit length_unit = Units::length_unit(cursor);
385
386 evt.set_units(momentum_unit,length_unit);
387
388 HEPMC3_DEBUG( 10, "ReaderAsciiHepMC2: U: " << Units::name(evt.momentum_unit()) << " " << Units::name(evt.length_unit()) )
389
390 return true;
391}
392
394 GenVertexPtr data = std::make_shared<GenVertex>();
395 GenVertexPtr data_ghost = std::make_shared<GenVertex>();
396 FourVector position;
397 const char *cursor = buf;
398 int barcode = 0;
399 int num_particles_out = 0;
400 int weights_size = 0;
401 std::vector<double> weights(0);
402 // barcode
403 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
404 barcode = atoi(cursor);
405
406 // status
407 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
408 data->set_status( atoi(cursor) );
409
410 // x
411 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
412 position.setX(atof(cursor));
413
414 // y
415 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
416 position.setY(atof(cursor));
417
418 // z
419 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
420 position.setZ(atof(cursor));
421
422 // t
423 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
424 position.setT(atof(cursor));
425 data->set_position( position );
426
427 // SKIPPED: num_orphans_in
428 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
429
430 // num_particles_out
431 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
432 num_particles_out = atoi(cursor);
433
434 // weights
435
436 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
437 weights_size = atoi(cursor);
438 weights.resize(weights_size);
439
440 for ( int i = 0; i < weights_size; ++i ) {
441 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
442 weights[i] = atof(cursor);
443 }
444
445
446
447 // Add original vertex barcode to the cache
448 m_vertex_cache.push_back( data );
449 m_vertex_barcodes.push_back( barcode );
450
451 m_event_ghost->add_vertex(data_ghost);
452 if (m_options.find("vertex_weights_are_separated")!=m_options.end())
453 {
454 for ( int i = 0; i < weights_size; ++i )
455 data_ghost->add_attribute("weight"+std::to_string((long long unsigned int)i),std::make_shared<DoubleAttribute>(weights[i]));
456 }
457 else
458 {
459 data_ghost->add_attribute("weights",std::make_shared<VectorDoubleAttribute>(weights));
460 data_ghost->add_attribute("weights",std::make_shared<VectorDoubleAttribute>(weights));
461 }
462 m_vertex_cache_ghost.push_back( data_ghost );
463
464 HEPMC3_DEBUG( 10, "ReaderAsciiHepMC2: V: "<<-(int)m_vertex_cache.size()<<" (old barcode"<<barcode<<") "<<num_particles_out<<" particles)" )
465
466 return num_particles_out;
467}
468
470 GenParticlePtr data = std::make_shared<GenParticle>();
471 GenParticlePtr data_ghost = std::make_shared<GenParticle>();
472 m_event_ghost->add_particle(data_ghost);
473 FourVector momentum;
474 const char *cursor = buf;
475 int end_vtx = 0;
476
477 /// @note barcode is ignored
478 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
479
480 // id
481 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
482 data->set_pid( atoi(cursor) );
483
484 // px
485 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
486 momentum.setPx(atof(cursor));
487
488 // py
489 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
490 momentum.setPy(atof(cursor));
491
492 // pz
493 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
494 momentum.setPz(atof(cursor));
495
496 // pe
497 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
498 momentum.setE(atof(cursor));
499 data->set_momentum(momentum);
500
501 // m
502 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
503 data->set_generated_mass( atof(cursor) );
504
505 // status
506 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
507 data->set_status( atoi(cursor) );
508
509 //theta
510 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
511 std::shared_ptr<DoubleAttribute> theta = std::make_shared<DoubleAttribute>(atof(cursor));
512 if (theta->value()!=0.0) data_ghost->add_attribute("theta",theta);
513
514 //phi
515 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
516 std::shared_ptr<DoubleAttribute> phi = std::make_shared<DoubleAttribute>(atof(cursor));
517 if (phi->value()!=0.0) data_ghost->add_attribute("phi",phi);
518
519 // end_vtx_code
520 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
521 end_vtx = atoi(cursor);
522
523 //flow
524 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
525 int flowsize=atoi(cursor);
526
527 std::map<int,int> flows;
528 for (int i=0; i<flowsize; i++)
529 {
530 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
531 int flowindex=atoi(cursor);
532 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
533 int flowvalue=atoi(cursor);
534 flows[flowindex]=flowvalue;
535 }
536 if (m_options.find("particle_flows_are_separated")!=m_options.end())
537 {
538 std::vector<int> vectorflows;
539 for (auto f: flows) vectorflows.push_back(f.second);
540 data_ghost->add_attribute("flows",std::make_shared<VectorIntAttribute>(vectorflows));
541 }
542 else
543 {
544 for (auto f: flows) data_ghost->add_attribute("flow"+std::to_string((long long int)f.first),std::make_shared<IntAttribute>(f.second));
545 }
546// Set prod_vtx link
547 if( end_vtx == m_vertex_barcodes.back() ) {
548 m_vertex_cache.back()->add_particle_in(data);
549 end_vtx = 0;
550 }
551 else {
552 m_vertex_cache.back()->add_particle_out(data);
553 }
554
555 m_particle_cache.push_back( data );
556 m_particle_cache_ghost.push_back( data_ghost );
557 m_end_vertex_barcodes.push_back( end_vtx );
558
559 HEPMC3_DEBUG( 10, "ReaderAsciiHepMC2: P: "<<m_particle_cache.size()<<" ( pid: "<<data->pid()<<") end vertex: "<<end_vtx )
560
561 return 0;
562}
563
564bool ReaderAsciiHepMC2::parse_xs_info(GenEvent &evt, const char *buf) {
565 const char *cursor = buf;
566 std::shared_ptr<GenCrossSection> xs = std::make_shared<GenCrossSection>();
567
568 if( !(cursor = strchr(cursor+1,' ')) ) return false;
569 double xs_val = atof(cursor);
570
571 if( !(cursor = strchr(cursor+1,' ')) ) return false;
572 double xs_err = atof(cursor);
573
574 xs->set_cross_section( xs_val, xs_err);
575 evt.add_attribute("GenCrossSection",xs);
576
577 return true;
578}
579
581 const char *cursor = buf;
582 const char *cursor2 = buf;
583 int w_count = 0;
584 std::vector<std::string> w_names;
585
586 // Ignore weight names if no GenRunInfo object
587 if( !run_info() ) return true;
588
589 if( !(cursor = strchr(cursor+1,' ')) ) return false;
590 w_count = atoi(cursor);
591
592 if( w_count <= 0 ) return false;
593
594 w_names.resize(w_count);
595
596 for( int i=0; i < w_count; ++i ) {
597 // Find pair of '"' characters
598 if( !(cursor = strchr(cursor+1,'"')) ) return false;
599 if( !(cursor2 = strchr(cursor+1,'"')) ) return false;
600
601 // Strip cursor of leading '"' character
602 ++cursor;
603
604 w_names[i].assign(cursor, cursor2-cursor);
605
606 cursor = cursor2;
607 }
608
609 run_info()->set_weight_names(w_names);
610
611 return true;
612}
613
614bool ReaderAsciiHepMC2::parse_heavy_ion(GenEvent &evt, const char *buf) {
615 std::shared_ptr<GenHeavyIon> hi = std::make_shared<GenHeavyIon>();
616 const char *cursor = buf;
617
618 if( !(cursor = strchr(cursor+1,' ')) ) return false;
619 hi->Ncoll_hard = atoi(cursor);
620
621 if( !(cursor = strchr(cursor+1,' ')) ) return false;
622 hi->Npart_proj = atoi(cursor);
623
624 if( !(cursor = strchr(cursor+1,' ')) ) return false;
625 hi->Npart_targ = atoi(cursor);
626
627 if( !(cursor = strchr(cursor+1,' ')) ) return false;
628 hi->Ncoll = atoi(cursor);
629
630 if( !(cursor = strchr(cursor+1,' ')) ) return false;
631 hi->spectator_neutrons = atoi(cursor);
632
633 if( !(cursor = strchr(cursor+1,' ')) ) return false;
634 hi->spectator_protons = atoi(cursor);
635
636 if( !(cursor = strchr(cursor+1,' ')) ) return false;
637 hi->N_Nwounded_collisions = atoi(cursor);
638
639 if( !(cursor = strchr(cursor+1,' ')) ) return false;
640 hi->Nwounded_N_collisions = atoi(cursor);
641
642 if( !(cursor = strchr(cursor+1,' ')) ) return false;
643 hi->Nwounded_Nwounded_collisions = atoi(cursor);
644
645 if( !(cursor = strchr(cursor+1,' ')) ) return false;
646 hi->impact_parameter = atof(cursor);
647
648 if( !(cursor = strchr(cursor+1,' ')) ) return false;
649 hi->event_plane_angle = atof(cursor);
650
651 if( !(cursor = strchr(cursor+1,' ')) ) return false;
652 hi->eccentricity = atof(cursor);
653
654 if( !(cursor = strchr(cursor+1,' ')) ) return false;
655 hi->sigma_inel_NN = atof(cursor);
656
657 // Not in HepMC2:
658 hi->centrality = 0.0;
659
660 evt.add_attribute("GenHeavyIon",hi);
661
662 return true;
663}
664
665bool ReaderAsciiHepMC2::parse_pdf_info(GenEvent &evt, const char *buf) {
666 std::shared_ptr<GenPdfInfo> pi = std::make_shared<GenPdfInfo>();
667 const char *cursor = buf;
668
669 if( !(cursor = strchr(cursor+1,' ')) ) return false;
670 pi->parton_id[0] = atoi(cursor);
671
672 if( !(cursor = strchr(cursor+1,' ')) ) return false;
673 pi->parton_id[1] = atoi(cursor);
674
675 if( !(cursor = strchr(cursor+1,' ')) ) return false;
676 pi->x[0] = atof(cursor);
677
678 if( !(cursor = strchr(cursor+1,' ')) ) return false;
679 pi->x[1] = atof(cursor);
680
681 if( !(cursor = strchr(cursor+1,' ')) ) return false;
682 pi->scale = atof(cursor);
683
684 if( !(cursor = strchr(cursor+1,' ')) ) return false;
685 pi->xf[0] = atof(cursor);
686
687 if( !(cursor = strchr(cursor+1,' ')) ) return false;
688 pi->xf[1] = atof(cursor);
689
690 //For compatibility with original HepMC2
691 bool pdfids=true;
692 if( !(cursor = strchr(cursor+1,' ')) ) pdfids=false;
693 if(pdfids) pi->pdf_id[0] = atoi(cursor);
694 else pi->pdf_id[0] =0;
695
696 if(pdfids) if( !(cursor = strchr(cursor+1,' ')) ) pdfids=false;
697 if(pdfids) pi->pdf_id[1] = atoi(cursor);
698 else pi->pdf_id[1] =0;
699
700 evt.add_attribute("GenPdfInfo",pi);
701
702 return true;
703}
704bool ReaderAsciiHepMC2::failed() { return m_isstream ? (bool)m_stream->rdstate() :(bool)m_file.rdstate(); }
705
707 if (m_event_ghost) { m_event_ghost->clear(); delete m_event_ghost; m_event_ghost=nullptr;}
708 if( !m_file.is_open() ) return;
709 m_file.close();
710}
711
712} // namespace HepMC3
#define HEPMC3_WARNING(MESSAGE)
Macro for printing HEPMC3_HEPMC3_WARNING messages.
Definition: Errors.h:26
#define HEPMC3_DEBUG(LEVEL, MESSAGE)
Macro for printing debug messages with appropriate debug level.
Definition: Errors.h:32
#define HEPMC3_ERROR(MESSAGE)
Macro for printing error messages.
Definition: Errors.h:23
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.
Definition: Attribute.h:243
Generic 4-vector.
Definition: FourVector.h:35
void setE(double ee)
Definition: FourVector.h:134
void setT(double tt)
Definition: FourVector.h:105
void setPz(double pzz)
Definition: FourVector.h:127
void setY(double yy)
Definition: FourVector.h:91
void setPy(double pyy)
Definition: FourVector.h:120
void setX(double xx)
Definition: FourVector.h:84
void setPx(double pxx)
Definition: FourVector.h:113
void setZ(double zz)
Definition: FourVector.h:98
Stores event-related information.
Definition: GenEvent.h:41
void add_tree(const std::vector< GenParticlePtr > &particles)
Add whole tree in topological order.
Definition: GenEvent.cc:267
void add_vertex(GenVertexPtr v)
Add vertex.
Definition: GenEvent.cc:97
void add_particle(GenParticlePtr p)
Add particle.
Definition: GenEvent.cc:49
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: GenEvent.cc:395
void set_event_number(const int &num)
Set event number.
Definition: GenEvent.h:137
void set_run_info(std::shared_ptr< GenRunInfo > run)
Set the GenRunInfo object by smart pointer.
Definition: GenEvent.h:128
const Units::MomentumUnit & momentum_unit() const
Get momentum unit.
Definition: GenEvent.h:140
const Units::LengthUnit & length_unit() const
Get length unit.
Definition: GenEvent.h:142
void add_attribute(const std::string &name, const std::shared_ptr< Attribute > &att, const int &id=0)
Add event attribute to event.
Definition: GenEvent.h:208
const std::vector< double > & weights() const
Get event weight values as a vector.
Definition: GenEvent.h:86
void clear()
Remove contents of this event.
Definition: GenEvent.cc:608
void reserve(const size_t &particles, const size_t &vertices=0)
Reserve memory for particles and vertices.
Definition: GenEvent.cc:389
Attribute that holds an Integer implemented as an int.
Definition: Attribute.h:159
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.
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.
Definition: Reader.h:44
std::map< std::string, std::string > m_options
options
Definition: Reader.h:68
void set_run_info(std::shared_ptr< GenRunInfo > run)
Set the global GenRunInfo object.
Definition: Reader.h:64
static LengthUnit length_unit(const std::string &name)
Get length unit based on its name.
Definition: Units.h:46
LengthUnit
Position units.
Definition: Units.h:32
static std::string name(MomentumUnit u)
Get name of momentum unit.
Definition: Units.h:56
MomentumUnit
Momentum units.
Definition: Units.h:29
static MomentumUnit momentum_unit(const std::string &name)
Get momentum unit based on its name.
Definition: Units.h:36
Attribute that holds a vector of integers of type int.
Definition: Attribute.h:1003
HepMC3 main namespace.