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