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 
20 namespace HepMC3
21 {
22 
23 
24 WriterAsciiHepMC2::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 
48 WriterAsciiHepMC2::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 
65 WriterAsciiHepMC2::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 
240  m_particle_counter = 0;
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 
282  m_cursor = m_buffer;
283 }
284 
285 
286 std::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 
307 void 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);
361  m_cursor = m_buffer;
362  }
363 }
364 
365 
367 {
368  std::ptrdiff_t length = m_cursor - m_buffer;
369  m_stream->write(m_buffer, length);
370  m_cursor = m_buffer;
371 }
372 
373 
375 
376 void 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 
429 inline 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 }
456 bool WriterAsciiHepMC2::failed() { return (bool)m_file.rdstate(); }
457 
458 void 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 
467 void 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
const std::vector< ConstGenVertexPtr > & vertices() const
Get list of vertices (const)
Definition: GenEvent.cc:43
std::string m_float_printf_specifier
the specifier of printf used for floats
#define HEPMC3_WARNING(MESSAGE)
Macro for printing HEPMC3_HEPMC3_WARNING messages.
Definition: Errors.h:27
Attribute that holds a vector of FPs of type double.
Definition: Attribute.h:1091
bool to_string(std::string &att) const override
Implementation of Attribute::to_string.
Definition: GenPdfInfo.cc:51
double t() const
Time component of position/displacement.
Definition: FourVector.h:102
Definition of class GenParticle.
void close() override
Close file stream.
std::ofstream m_file
Output file.
WriterAsciiHepMC2(const std::string &filename, std::shared_ptr< GenRunInfo > run=std::shared_ptr< GenRunInfo >())
Constructor.
void write_event(const GenEvent &evt) override
Write event to file.
void forced_flush()
Inline function forcing flush to the output stream.
Definition of class GenVertex.
void write_vertex(ConstGenVertexPtr v)
Write vertex.
bool is_zero() const
Check if the length of this vertex is zero.
Definition: FourVector.h:193
Attribute that holds a vector of integers of type int.
Definition: Attribute.h:1001
std::string version()
Get the HepMC library version string.
Definition: Version.h:20
const Units::LengthUnit & length_unit() const
Get length unit.
Definition: GenEvent.h:155
double z() const
z-component of position/displacement
Definition: FourVector.h:95
double x() const
x-component of position/displacement
Definition: FourVector.h:81
static std::string name(MomentumUnit u)
Get name of momentum unit.
Definition: Units.h:56
int event_number() const
Get event number.
Definition: GenEvent.h:148
int precision() const
Return output precision.
char * m_buffer
Stream buffer.
char * m_cursor
Cursor inside stream buffer.
std::vector< int > value() const
get the value associated to this Attribute.
Definition: Attribute.h:1027
unsigned long m_buffer_size
Buffer size.
unsigned long m_particle_counter
Used to set bar codes.
Stores additional information about PDFs.
Definition: GenPdfInfo.h:32
const Units::MomentumUnit & momentum_unit() const
Get momentum unit.
Definition: GenEvent.h:153
Stores event-related information.
Definition: GenEvent.h:41
void allocate_buffer()
Attempts to allocate buffer of the chosen size.
Generic 4-vector.
Definition: FourVector.h:36
Attribute that holds a real number as a double.
Definition: Attribute.h:241
void write_string(const std::string &str)
Inline function for writing strings.
bool failed() override
Return status of the stream.
std::shared_ptr< GenRunInfo > run_info() const
Get the global GenRunInfo object.
Definition: Writer.h:47
void write_run_info()
Write the GenRunInfo object to file.
void set_buffer_size(const size_t &size)
Set buffer size (in bytes)
void write_particle(ConstGenParticlePtr p, int second_field)
Write particle.
Definition of class Units.
std::shared_ptr< T > attribute(const std::string &name, const int &id=0) const
Get attribute of type T.
Definition: GenEvent.h:409
Definition of class WriterAsciiHepMC2.
void set_run_info(std::shared_ptr< GenRunInfo > run)
Set the global GenRunInfo object.
Definition: Writer.h:42
Stores additional information about Heavy Ion generator.
Definition: GenHeavyIon.h:28
void flush()
Inline function flushing buffer to output stream when close to buffer capacity.
std::string escape(const std::string &s) const
Escape &#39;\&#39; and &#39; &#39; characters in string.
double y() const
y-component of position/displacement
Definition: FourVector.h:88
#define HEPMC3_ERROR(MESSAGE)
Macro for printing error messages.
Definition: Errors.h:24
std::vector< long int > value() const
get the value associated to this Attribute.
Definition: Attribute.h:1072
const std::vector< double > & weights() const
Get event weight values as a vector.
Definition: GenEvent.h:98
Attribute that holds a vector of integers of type int.
Definition: Attribute.h:1046
int m_precision
Output precision.
Attribute that holds an Integer implemented as an int.
Definition: Attribute.h:198
Definition of class GenEvent.
void set_precision(const int &prec)
Set output precision.
std::map< std::string, std::string > m_options
options
Definition: Writer.h:68
std::shared_ptr< GenRunInfo > run_info() const
Get a pointer to the the GenRunInfo object.
Definition: GenEvent.h:137
Attribute that holds an Integer implemented as an int.
Definition: Attribute.h:157
std::ostream * m_stream
Output stream.
Stores additional information about cross-section.