HepMC3 event record library
WriterAsciiHepMC2.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 WriterAsciiHepMC2.cc
8/// @brief Implementation of \b class WriterAsciiHepMC2
9///
10#include <cstring>
11
13
14#include "HepMC3/Version.h"
15#include "HepMC3/GenEvent.h"
16#include "HepMC3/GenParticle.h"
17#include "HepMC3/GenVertex.h"
18#include "HepMC3/Units.h"
19
20namespace HepMC3
21{
22
23
24WriterAsciiHepMC2::WriterAsciiHepMC2(const std::string &filename, std::shared_ptr<GenRunInfo> run)
25 : m_file(filename),
26 m_stream(&m_file),
27 m_precision(16),
28 m_buffer(nullptr),
29 m_cursor(nullptr),
30 m_buffer_size(256*1024),
31 m_particle_counter(0)
32{
33 HEPMC3_WARNING("WriterAsciiHepMC2::WriterAsciiHepMC2: HepMC2 IO_GenEvent format is outdated. Please use HepMC3 Asciiv3 format instead.")
34 set_run_info(run);
35 if ( !run_info() ) set_run_info(std::make_shared<GenRunInfo>());
36 if ( !m_file.is_open() )
37 {
38 HEPMC3_ERROR("WriterAsciiHepMC2: could not open output file: " << filename )
39 }
40 else
41 {
42 const std::string header = "HepMC::Version " + version() + "\nHepMC::IO_GenEvent-START_EVENT_LISTING\n";
43 m_file.write(header.data(), header.length());
44 }
45 m_float_printf_specifier = " %." + std::to_string(m_precision) + "e";
46}
47
48WriterAsciiHepMC2::WriterAsciiHepMC2(std::ostream &stream, std::shared_ptr<GenRunInfo> run)
49 : m_file(),
50 m_stream(&stream),
51 m_precision(16),
52 m_buffer(nullptr),
53 m_cursor(nullptr),
54 m_buffer_size(256*1024),
55 m_particle_counter(0)
56{
57 HEPMC3_WARNING("WriterAsciiHepMC2::WriterAsciiHepMC2: HepMC2 IO_GenEvent format is outdated. Please use HepMC3 Asciiv3 format instead.")
58 set_run_info(run);
59 if ( !run_info() ) set_run_info(std::make_shared<GenRunInfo>());
60 const std::string header = "HepMC::Version " + version() + "\nHepMC::IO_GenEvent-START_EVENT_LISTING\n";
61 m_stream->write(header.data(), header.length());
62 m_float_printf_specifier = " %." + std::to_string(m_precision) + "e";
63}
64
65WriterAsciiHepMC2::WriterAsciiHepMC2(std::shared_ptr<std::ostream> s_stream, std::shared_ptr<GenRunInfo> run)
66 : m_file(),
67 m_shared_stream(s_stream),
68 m_stream(s_stream.get()),
69 m_precision(16),
70 m_buffer(nullptr),
71 m_cursor(nullptr),
72 m_buffer_size(256*1024),
73 m_particle_counter(0)
74{
75 HEPMC3_WARNING("WriterAsciiHepMC2::WriterAsciiHepMC2: HepMC2 IO_GenEvent format is outdated. Please use HepMC3 Asciiv3 format instead.")
76 set_run_info(run);
77 if ( !run_info() ) set_run_info(std::make_shared<GenRunInfo>());
78 const std::string header = "HepMC::Version " + version() + "\nHepMC::IO_GenEvent-START_EVENT_LISTING\n";
79 m_stream->write(header.data(), header.length());
80 m_float_printf_specifier = " %." + std::to_string(m_precision) + "e";
81}
82
83
85{
86 close();
87 if ( m_buffer ) delete[] m_buffer;
88}
89
90
92{
94 if ( !m_buffer ) return;
95 auto float_printf_specifier_option = m_options.find("float_printf_specifier");
96 std::string letter=(float_printf_specifier_option != m_options.end())?float_printf_specifier_option->second.substr(0,2):"e";
97 if (letter != "e" && letter != "E" && letter != "G" && letter != "g" && letter != "f" && letter != "F" ) letter = "e";
98 m_float_printf_specifier = " %." + std::to_string(m_precision) + letter;
99 // Make sure nothing was left from previous event
100 flush();
101
102 if ( !run_info() ) set_run_info(evt.run_info());
103 if ( evt.run_info() && run_info() != evt.run_info() ) set_run_info(evt.run_info());
104
105
106 std::shared_ptr<DoubleAttribute> A_event_scale = evt.attribute<DoubleAttribute>("event_scale");
107 std::shared_ptr<DoubleAttribute> A_alphaQED = evt.attribute<DoubleAttribute>("alphaQED");
108 std::shared_ptr<DoubleAttribute> A_alphaQCD = evt.attribute<DoubleAttribute>("alphaQCD");
109 std::shared_ptr<IntAttribute> A_signal_process_id = evt.attribute<IntAttribute>("signal_process_id");
110 std::shared_ptr<IntAttribute> A_mpi = evt.attribute<IntAttribute>("mpi");
111 std::shared_ptr<IntAttribute> A_signal_process_vertex = evt.attribute<IntAttribute>("signal_process_vertex");
112
113 double event_scale = A_event_scale?(A_event_scale->value()):0.0;
114 double alphaQED = A_alphaQED?(A_alphaQED->value()):0.0;
115 double alphaQCD = A_alphaQCD?(A_alphaQCD->value()):0.0;
116 int signal_process_id = A_signal_process_id?(A_signal_process_id->value()):0;
117 int mpi = A_mpi?(A_mpi->value()):0;
118 int signal_process_vertex = A_signal_process_vertex?(A_signal_process_vertex->value()):0;
119
120 std::vector<long> m_random_states;
121 std::shared_ptr<VectorLongIntAttribute> random_states_a = evt.attribute<VectorLongIntAttribute>("random_states");
122 if (random_states_a) {
123 m_random_states = random_states_a->value();
124 } else {
125 m_random_states.reserve(100);
126 for (int i = 0; i < 100; i++)
127 {
128 std::shared_ptr<LongAttribute> rs = evt.attribute<LongAttribute>("random_states"+std::to_string((long long unsigned int)i));
129 if (!rs) break;
130 m_random_states.push_back(rs->value());
131 }
132 }
133 // Write event info
134 //Find beam particles
135 std::vector<int> beams;
136 beams.reserve(2);
137 int idbeam = 0;
138 for (ConstGenVertexPtr v: evt.vertices())
139 {
140 for (ConstGenParticlePtr p: v->particles_in())
141 {
142 if (!p->production_vertex()) { if (p->status() == 4) beams.push_back(idbeam); idbeam++; }
143 else if (p->production_vertex()->id() == 0) { if (p->status() == 4) beams.push_back(idbeam); idbeam++; }
144 }
145 for (ConstGenParticlePtr p: v->particles_out()) { if (p->status() == 4) beams.push_back(idbeam); idbeam++; }
146 }
147 //
148 int idbeam1 = 10000;
149 int idbeam2 = 10000;
150 if (beams.size() > 0) idbeam1 += beams[0] + 1;
151 if (beams.size() > 1) idbeam2 += beams[1] + 1;
152 m_cursor += sprintf(m_cursor, "E %d %d %e %e %e %d %d %zu %i %i",
153 evt.event_number(),
154 mpi,
155 event_scale,
156 alphaQCD,
157 alphaQED,
158 signal_process_id,
159 signal_process_vertex,
160 evt.vertices().size(),
161 idbeam1, idbeam2);
162
163 // This should be the largest single add to the buffer. Its size 11+4*11+3*22+2*11+10=153
164 flush();
165 m_cursor += sprintf(m_cursor, " %zu", m_random_states.size());
166 for (size_t q = 0; q < m_random_states.size(); q++)
167 {
168 m_cursor += sprintf(m_cursor, " %i", (int)q);
169 flush();
170 }
171 flush();
172 m_cursor += sprintf(m_cursor, " %zu", evt.weights().size());
173 if ( evt.weights().size() )
174 {
175 for (double w: evt.weights()) {
176 m_cursor += sprintf(m_cursor, m_float_printf_specifier.c_str(), w);
177 flush();
178 }
179 m_cursor += sprintf(m_cursor, "\n");
180 flush();
181 m_cursor += sprintf(m_cursor, "N %zu", evt.weights().size());
182 const std::vector<std::string> names = run_info()->weight_names();
183 for (size_t q = 0; q < evt.weights().size(); q++)
184 {
185 if (q < names.size())
186 write_string(" \""+names[q]+"\"");
187 else
188 write_string(" \""+std::to_string(q)+"\"");
189 flush();
190 }
191 }
192 m_cursor += sprintf(m_cursor, "\n");
193 flush();
194 // Write units
195 m_cursor += sprintf(m_cursor, "U %s %s\n", Units::name(evt.momentum_unit()).c_str(), Units::name(evt.length_unit()).c_str());
196 flush();
197 std::shared_ptr<GenCrossSection> cs = evt.attribute<GenCrossSection>("GenCrossSection");
198 if (cs) {
199 m_cursor += sprintf(m_cursor, ("C" + m_float_printf_specifier + m_float_printf_specifier + "\n").c_str(), cs->xsec(), cs->xsec_err() );
200 flush();
201 }
202
203 std::shared_ptr<GenHeavyIon> hi = evt.attribute<GenHeavyIon>("GenHeavyIon");
204 if (hi) {
205 m_cursor += sprintf(m_cursor, "H %i %i %i %i %i %i %i %i %i %e %e %e %e\n",
206 hi->Ncoll_hard,
207 hi->Npart_proj,
208 hi->Npart_targ,
209 hi->Ncoll,
210 hi->spectator_neutrons,
211 hi->spectator_protons,
212 hi->N_Nwounded_collisions,
213 hi->Nwounded_N_collisions,
214 hi->Nwounded_Nwounded_collisions,
215 hi->impact_parameter,
216 hi->event_plane_angle,
217 hi->eccentricity,
218 hi->sigma_inel_NN);
219 flush();
220 }
221
222 std::shared_ptr<GenPdfInfo> pi = evt.attribute<GenPdfInfo>("GenPdfInfo");
223 if (pi) {
224 std::string st;
225 // We use it here because the HepMC3 GenPdfInfo has the same format as in HepMC2 IO_GenEvent and get error handeling for free.
226 bool status = pi->to_string(st);
227 if ( !status )
228 {
229 HEPMC3_WARNING("WriterAsciiHepMC2::write_event: problem serializing GenPdfInfo attribute")
230 } else {
231 m_cursor += sprintf(m_cursor, "F ");
232 flush();
233 write_string(escape(st));
234 m_cursor += sprintf(m_cursor, "\n");
235 flush();
236 }
237 }
238
239
241 for (ConstGenVertexPtr v: evt.vertices() )
242 {
243 int production_vertex = 0;
244 production_vertex = v->id();
245 write_vertex(v);
246 for (ConstGenParticlePtr p: v->particles_in())
247 {
248 if (!p->production_vertex()) write_particle( p, production_vertex );
249 else
250 {
251 if (p->production_vertex()->id() == 0) write_particle( p, production_vertex );
252 }
253 }
254 for (ConstGenParticlePtr p: v->particles_out())
255 write_particle(p, production_vertex);
256 }
257
258 // Flush rest of the buffer to file
259 forced_flush();
260}
261
262
264{
265 if ( m_buffer ) return;
266 while ( m_buffer == nullptr && m_buffer_size >= 512 ) {
267 try {
268 m_buffer = new char[ m_buffer_size ]();
269 } catch (const std::bad_alloc& e) {
270 delete[] m_buffer;
271 m_buffer_size /= 2;
272 HEPMC3_WARNING("WriterAsciiHepMC2::allocate_buffer:" << e.what() << " buffer size too large. Dividing by 2. New size: " << m_buffer_size)
273 }
274 }
275
276 if ( !m_buffer )
277 {
278 HEPMC3_ERROR("WriterAsciiHepMC2::allocate_buffer: could not allocate buffer!")
279 return;
280 }
281
283}
284
285
286std::string WriterAsciiHepMC2::escape(const std::string& s) const
287{
288 std::string ret;
289 ret.reserve(s.length()*2);
290 for ( std::string::const_iterator it = s.begin(); it != s.end(); ++it )
291 {
292 switch ( *it )
293 {
294 case '\\':
295 ret += "\\\\";
296 break;
297 case '\n':
298 ret += "\\|";
299 break;
300 default:
301 ret += *it;
302 }
303 }
304 return ret;
305}
306
307void WriterAsciiHepMC2::write_vertex(ConstGenVertexPtr v)
308{
309 std::vector<double> weights;
310 std::shared_ptr<VectorDoubleAttribute> weights_a = v->attribute<VectorDoubleAttribute>("weights");
311 if (weights_a) {
312 weights = weights_a->value();
313 } else {
314 weights.reserve(100);
315 for (int i = 0; i < 100; i++)
316 {
317 std::shared_ptr<DoubleAttribute> rs = v->attribute<DoubleAttribute>("weight"+std::to_string((long long unsigned int)i));
318 if (!rs) break;
319 weights.push_back(rs->value());
320 }
321 }
322 flush();
323 m_cursor += sprintf(m_cursor, "V %i %i", v->id(), v->status());
324 int orph = 0;
325 for (ConstGenParticlePtr p: v->particles_in())
326 {
327 if (!p->production_vertex()) orph++;
328 else
329 {
330 if (p->production_vertex()->id() == 0) orph++;
331 }
332 }
333 const FourVector &pos = v->position();
334 if (pos.is_zero())
335 {
336 m_cursor += sprintf(m_cursor, " 0 0 0 0");
337 }
338 else
339 {
340 m_cursor += sprintf(m_cursor, m_float_printf_specifier.c_str(), pos.x());
341 m_cursor += sprintf(m_cursor, m_float_printf_specifier.c_str(), pos.y());
342 m_cursor += sprintf(m_cursor, m_float_printf_specifier.c_str(), pos.z());
343 m_cursor += sprintf(m_cursor, m_float_printf_specifier.c_str(), pos.t());
344 }
345 m_cursor += sprintf(m_cursor, " %i %zu %zu", orph, v->particles_out().size(), weights.size());
346 flush();
347 for (size_t i = 0; i < weights.size(); i++) { m_cursor += sprintf(m_cursor, m_float_printf_specifier.c_str(), weights[i]); flush(); }
348 m_cursor += sprintf(m_cursor, "\n");
349 flush();
350}
351
352
354{
355 // The maximum size of single add to the buffer should not be larger than 256. This is a safe value as
356 // we will not allow precision larger than 24 anyway
357 if ( m_buffer + m_buffer_size < m_cursor + 512 )
358 {
359 std::ptrdiff_t length = m_cursor - m_buffer;
360 m_stream->write(m_buffer, length);
362 }
363}
364
365
367{
368 std::ptrdiff_t length = m_cursor - m_buffer;
369 m_stream->write(m_buffer, length);
371}
372
373
375
376void WriterAsciiHepMC2::write_particle(ConstGenParticlePtr p, int /*second_field*/)
377{
378 flush();
379 m_cursor += sprintf(m_cursor, "P %i", int(10001+m_particle_counter));
381 m_cursor += sprintf(m_cursor, " %i", p->pid() );
382 m_cursor += sprintf(m_cursor, m_float_printf_specifier.c_str(), p->momentum().px() );
383 m_cursor += sprintf(m_cursor, m_float_printf_specifier.c_str(), p->momentum().py());
384 m_cursor += sprintf(m_cursor, m_float_printf_specifier.c_str(), p->momentum().pz() );
385 m_cursor += sprintf(m_cursor, m_float_printf_specifier.c_str(), p->momentum().e() );
386 m_cursor += sprintf(m_cursor, m_float_printf_specifier.c_str(), p->generated_mass() );
387 m_cursor += sprintf(m_cursor, " %i", p->status() );
388 flush();
389 int ev = 0;
390 if (p->end_vertex())
391 if (p->end_vertex()->id() != 0)
392 ev = p->end_vertex()->id();
393
394 std::shared_ptr<DoubleAttribute> A_theta = p->attribute<DoubleAttribute>("theta");
395 std:: shared_ptr<DoubleAttribute> A_phi = p->attribute<DoubleAttribute>("phi");
396 if (A_theta) m_cursor += sprintf(m_cursor, m_float_printf_specifier.c_str(), A_theta->value());
397 else m_cursor += sprintf(m_cursor, " 0");
398 if (A_phi) m_cursor += sprintf(m_cursor, m_float_printf_specifier.c_str(), A_phi->value());
399 else m_cursor += sprintf(m_cursor, " 0");
400 m_cursor += sprintf(m_cursor, " %i", ev);
401 flush();
402 std::shared_ptr<VectorIntAttribute> A_flows = p->attribute<VectorIntAttribute>("flows");
403 if (A_flows)
404 {
405 std::vector<int> flowsv = A_flows->value();
406 std::string flowss = " " + std::to_string(flowsv.size());
407 for (size_t k = 0; k < flowsv.size(); k++) { flowss += ( " " + std::to_string(k+1) + " " + std::to_string(flowsv.at(k))); }
408 flowss += "\n";
409 write_string(flowss);
410 } else {
411 std::shared_ptr<IntAttribute> A_flow1 = p->attribute<IntAttribute>("flow1");
412 std::shared_ptr<IntAttribute> A_flow2 = p->attribute<IntAttribute>("flow2");
413 std::shared_ptr<IntAttribute> A_flow3 = p->attribute<IntAttribute>("flow3");
414 int flowsize = 0;
415 if (A_flow1) flowsize++;
416 if (A_flow2) flowsize++;
417 if (A_flow3) flowsize++;
418 std::string flowss = " " + std::to_string(flowsize);
419 if (A_flow1) flowss += ( " 1 " + std::to_string(A_flow1->value()));
420 if (A_flow2) flowss += ( " 2 " + std::to_string(A_flow2->value()));
421 if (A_flow3) flowss += ( " 3 " + std::to_string(A_flow3->value()));
422 flowss += "\n";
423 write_string(flowss);
424 }
425 flush();
426}
427
428
429inline void WriterAsciiHepMC2::write_string(const std::string &str)
430{
431 // First let's check if string will fit into the buffer
432 if ( m_buffer + m_buffer_size > m_cursor + str.length() )
433 {
434 strncpy(m_cursor, str.data(), str.length());
435 m_cursor += str.length();
436 flush();
437 }
438 // If not, flush the buffer and write the string directly
439 else
440 {
441 forced_flush();
442 m_stream->write(str.data(), str.length());
443 }
444}
445
446
448{
449 std::ofstream* ofs = dynamic_cast<std::ofstream*>(m_stream);
450 if (ofs && !ofs->is_open()) return;
451 forced_flush();
452 const std::string footer("HepMC::IO_GenEvent-END_EVENT_LISTING\n\n");
453 if (m_stream) m_stream->write(footer.data(),footer.length());
454 if (ofs) ofs->close();
455}
456bool WriterAsciiHepMC2::failed() { return (bool)m_file.rdstate(); }
457
458void WriterAsciiHepMC2::set_precision(const int& prec ) {
459 if (prec < 2 || prec > 24) return;
460 m_precision = prec;
461}
462
464 return m_precision;
465}
466
467void WriterAsciiHepMC2::set_buffer_size(const size_t& size ) {
468 if (m_buffer) return;
469 if (size < 1024) return;
470 m_buffer_size = size;
471}
472} // namespace HepMC3
#define HEPMC3_WARNING(MESSAGE)
Macro for printing HEPMC3_HEPMC3_WARNING messages.
Definition: Errors.h:27
#define HEPMC3_ERROR(MESSAGE)
Macro for printing error messages.
Definition: Errors.h:24
Definition of class GenEvent.
Definition of class GenParticle.
Definition of class GenVertex.
Definition of class Units.
Definition of class WriterAsciiHepMC2.
Attribute that holds a real number as a double.
Definition: Attribute.h:241
Generic 4-vector.
Definition: FourVector.h:36
double t() const
Time component of position/displacement.
Definition: FourVector.h:102
bool is_zero() const
Check if the length of this vertex is zero.
Definition: FourVector.h:193
double x() const
x-component of position/displacement
Definition: FourVector.h:81
double y() const
y-component of position/displacement
Definition: FourVector.h:88
double z() const
z-component of position/displacement
Definition: FourVector.h:95
Stores additional information about cross-section.
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
int event_number() const
Get event number.
Definition: GenEvent.h:148
std::shared_ptr< GenRunInfo > run_info() const
Get a pointer to the the GenRunInfo object.
Definition: GenEvent.h:137
const std::vector< ConstGenVertexPtr > & vertices() const
Get list of vertices (const)
Definition: GenEvent.cc:43
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
const std::vector< double > & weights() const
Get event weight values as a vector.
Definition: GenEvent.h:98
Stores additional information about Heavy Ion generator.
Definition: GenHeavyIon.h:28
Stores additional information about PDFs.
Definition: GenPdfInfo.h:32
bool to_string(std::string &att) const override
Implementation of Attribute::to_string.
Definition: GenPdfInfo.cc:51
Attribute that holds an Integer implemented as an int.
Definition: Attribute.h:157
Attribute that holds an Integer implemented as an int.
Definition: Attribute.h:198
static std::string name(MomentumUnit u)
Get name of momentum unit.
Definition: Units.h:56
Attribute that holds a vector of FPs of type double.
Definition: Attribute.h:1091
std::vector< double > value() const
get the value associated to this Attribute.
Definition: Attribute.h:1117
Attribute that holds a vector of integers of type int.
Definition: Attribute.h:1001
std::vector< int > value() const
get the value associated to this Attribute.
Definition: Attribute.h:1027
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
void set_buffer_size(const size_t &size)
Set buffer size (in bytes)
void set_precision(const int &prec)
Set output precision.
std::string escape(const std::string &s) const
Escape '\' and ' ' characters in string.
void allocate_buffer()
Attempts to allocate buffer of the chosen size.
char * m_cursor
Cursor inside stream buffer.
bool failed() override
Return status of the stream.
char * m_buffer
Stream buffer.
std::string m_float_printf_specifier
the specifier of printf used for floats
void close() override
Close file stream.
int precision() const
Return output precision.
void write_particle(ConstGenParticlePtr p, int second_field)
Write particle.
int m_precision
Output precision.
WriterAsciiHepMC2(const std::string &filename, std::shared_ptr< GenRunInfo > run=std::shared_ptr< GenRunInfo >())
Constructor.
std::ofstream m_file
Output file.
unsigned long m_particle_counter
Used to set bar codes.
void write_string(const std::string &str)
Inline function for writing strings.
unsigned long m_buffer_size
Buffer size.
void write_event(const GenEvent &evt) override
Write event to file.
void write_vertex(ConstGenVertexPtr v)
Write vertex.
void flush()
Inline function flushing buffer to output stream when close to buffer capacity.
void write_run_info()
Write the GenRunInfo object to file.
void forced_flush()
Inline function forcing flush to the output stream.
std::ostream * m_stream
Output stream.
std::shared_ptr< GenRunInfo > run_info() const
Get the global GenRunInfo object.
Definition: Writer.h:47
std::map< std::string, std::string > m_options
options
Definition: Writer.h:68
void set_run_info(std::shared_ptr< GenRunInfo > run)
Set the global GenRunInfo object.
Definition: Writer.h:42
HepMC3 main namespace.
std::string version()
Get the HepMC library version string.
Definition: Version.h:20