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-2020 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  m_file << "HepMC::Version " << version() << std::endl;
43  m_file << "HepMC::IO_GenEvent-START_EVENT_LISTING" << std::endl;
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  (*m_stream) << "HepMC::Version " << version() << std::endl;
61  (*m_stream) << "HepMC::IO_GenEvent-START_EVENT_LISTING" << std::endl;
62  m_float_printf_specifier = " %." + std::to_string(m_precision) + "e";
63 }
64 
65 
67 {
68  close();
69  if ( m_buffer ) delete[] m_buffer;
70 }
71 
72 
74 {
76  if ( !m_buffer ) return;
77  auto float_printf_specifier_option = m_options.find("float_printf_specifier");
78  std::string letter=(float_printf_specifier_option != m_options.end())?float_printf_specifier_option->second.substr(0,2):"e";
79  if (letter != "e" && letter != "E" && letter != "G" && letter != "g" && letter != "f" && letter != "F" ) letter = "e";
80  m_float_printf_specifier = " %." + std::to_string(m_precision) + letter;
81  // Make sure nothing was left from previous event
82  flush();
83 
84  if ( !run_info() ) set_run_info(evt.run_info());
85  if ( evt.run_info() && run_info() != evt.run_info() ) set_run_info(evt.run_info());
86 
87 
88  std::shared_ptr<DoubleAttribute> A_event_scale = evt.attribute<DoubleAttribute>("event_scale");
89  std::shared_ptr<DoubleAttribute> A_alphaQED = evt.attribute<DoubleAttribute>("alphaQED");
90  std::shared_ptr<DoubleAttribute> A_alphaQCD = evt.attribute<DoubleAttribute>("alphaQCD");
91  std::shared_ptr<IntAttribute> A_signal_process_id = evt.attribute<IntAttribute>("signal_process_id");
92  std::shared_ptr<IntAttribute> A_mpi = evt.attribute<IntAttribute>("mpi");
93  std::shared_ptr<IntAttribute> A_signal_process_vertex = evt.attribute<IntAttribute>("signal_process_vertex");
94 
95  double event_scale = A_event_scale?(A_event_scale->value()):0.0;
96  double alphaQED = A_alphaQED?(A_alphaQED->value()):0.0;
97  double alphaQCD = A_alphaQCD?(A_alphaQCD->value()):0.0;
98  int signal_process_id = A_signal_process_id?(A_signal_process_id->value()):0;
99  int mpi = A_mpi?(A_mpi->value()):0;
100  int signal_process_vertex = A_signal_process_vertex?(A_signal_process_vertex->value()):0;
101 
102  std::vector<long> m_random_states;
103  std::shared_ptr<VectorLongIntAttribute> random_states_a = evt.attribute<VectorLongIntAttribute>("random_states");
104  if (random_states_a) {
105  m_random_states = random_states_a->value();
106  } else {
107  m_random_states.reserve(100);
108  for (int i = 0; i < 100; i++)
109  {
110  std::shared_ptr<LongAttribute> rs = evt.attribute<LongAttribute>("random_states"+std::to_string((long long unsigned int)i));
111  if (!rs) break;
112  m_random_states.push_back(rs->value());
113  }
114  }
115  // Write event info
116  //Find beam particles
117  std::vector<int> beams;
118  beams.reserve(2);
119  int idbeam = 0;
120  for (ConstGenVertexPtr v: evt.vertices())
121  {
122  for (ConstGenParticlePtr p: v->particles_in())
123  {
124  if (!p->production_vertex()) { if (p->status() == 4) beams.push_back(idbeam); idbeam++; }
125  else if (p->production_vertex()->id() == 0) { if (p->status() == 4) beams.push_back(idbeam); idbeam++; }
126  }
127  for (ConstGenParticlePtr p: v->particles_out()) { if (p->status() == 4) beams.push_back(idbeam); idbeam++; }
128  }
129  //
130  int idbeam1 = 10000;
131  int idbeam2 = 10000;
132  if (beams.size() > 0) idbeam1 += beams[0] + 1;
133  if (beams.size() > 1) idbeam2 += beams[1] + 1;
134  m_cursor += sprintf(m_cursor, "E %d %d %e %e %e %d %d %zu %i %i",
135  evt.event_number(),
136  mpi,
137  event_scale,
138  alphaQCD,
139  alphaQED,
140  signal_process_id,
141  signal_process_vertex,
142  evt.vertices().size(),
143  idbeam1, idbeam2);
144 
145  // This should be the largest single add to the buffer. Its size 11+4*11+3*22+2*11+10=153
146  flush();
147  m_cursor += sprintf(m_cursor, " %zu", m_random_states.size());
148  for (size_t q = 0; q < m_random_states.size(); q++)
149  {
150  m_cursor += sprintf(m_cursor, " %i", (int)q);
151  flush();
152  }
153  flush();
154  m_cursor += sprintf(m_cursor, " %zu", evt.weights().size());
155  if ( evt.weights().size() )
156  {
157  for (double w: evt.weights()) {
158  m_cursor += sprintf(m_cursor, m_float_printf_specifier.c_str(), w);
159  flush();
160  }
161  m_cursor += sprintf(m_cursor, "\n");
162  flush();
163  m_cursor += sprintf(m_cursor, "N %zu", evt.weights().size());
164  const std::vector<std::string> names = run_info()->weight_names();
165  for (size_t q = 0; q < evt.weights().size(); q++)
166  {
167  if (q < names.size())
168  write_string(" \""+names[q]+"\"");
169  else
170  write_string(" \""+std::to_string(q)+"\"");
171  flush();
172  }
173  }
174  m_cursor += sprintf(m_cursor, "\n");
175  flush();
176  // Write units
177  m_cursor += sprintf(m_cursor, "U %s %s\n", Units::name(evt.momentum_unit()).c_str(), Units::name(evt.length_unit()).c_str());
178  flush();
179  std::shared_ptr<GenCrossSection> cs = evt.attribute<GenCrossSection>("GenCrossSection");
180  if (cs) {
181  m_cursor += sprintf(m_cursor, ("C" + m_float_printf_specifier + m_float_printf_specifier + "\n").c_str(), cs->xsec(), cs->xsec_err() );
182  flush();
183  }
184 
185  std::shared_ptr<GenHeavyIon> hi = evt.attribute<GenHeavyIon>("GenHeavyIon");
186  if (hi) {
187  m_cursor += sprintf(m_cursor, "H %i %i %i %i %i %i %i %i %i %e %e %e %e\n",
188  hi->Ncoll_hard,
189  hi->Npart_proj,
190  hi->Npart_targ,
191  hi->Ncoll,
192  hi->spectator_neutrons,
193  hi->spectator_protons,
194  hi->N_Nwounded_collisions,
195  hi->Nwounded_N_collisions,
196  hi->Nwounded_Nwounded_collisions,
197  hi->impact_parameter,
198  hi->event_plane_angle,
199  hi->eccentricity,
200  hi->sigma_inel_NN);
201  flush();
202  }
203 
204  std::shared_ptr<GenPdfInfo> pi = evt.attribute<GenPdfInfo>("GenPdfInfo");
205  if (pi) {
206  std::string st;
207  // We use it here because the HepMC3 GenPdfInfo has the same format as in HepMC2 IO_GenEvent and get error handeling for free.
208  bool status = pi->to_string(st);
209  if ( !status )
210  {
211  HEPMC3_WARNING("WriterAsciiHepMC2::write_event: problem serializing GenPdfInfo attribute")
212  } else {
213  m_cursor += sprintf(m_cursor, "F ");
214  flush();
215  write_string(escape(st));
216  m_cursor += sprintf(m_cursor, "\n");
217  flush();
218  }
219  }
220 
221 
222  m_particle_counter = 0;
223  for (ConstGenVertexPtr v: evt.vertices() )
224  {
225  int production_vertex = 0;
226  production_vertex = v->id();
227  write_vertex(v);
228  for (ConstGenParticlePtr p: v->particles_in())
229  {
230  if (!p->production_vertex()) write_particle( p, production_vertex );
231  else
232  {
233  if (p->production_vertex()->id() == 0) write_particle( p, production_vertex );
234  }
235  }
236  for (ConstGenParticlePtr p: v->particles_out())
237  write_particle(p, production_vertex);
238  }
239 
240  // Flush rest of the buffer to file
241  forced_flush();
242 }
243 
244 
246 {
247  if ( m_buffer ) return;
248  while ( m_buffer == nullptr && m_buffer_size >= 512 ) {
249  try {
250  m_buffer = new char[ m_buffer_size ]();
251  } catch (const std::bad_alloc& e) {
252  delete[] m_buffer;
253  m_buffer_size /= 2;
254  HEPMC3_WARNING("WriterAsciiHepMC2::allocate_buffer:" << e.what() << " buffer size too large. Dividing by 2. New size: " << m_buffer_size)
255  }
256  }
257 
258  if ( !m_buffer )
259  {
260  HEPMC3_ERROR("WriterAsciiHepMC2::allocate_buffer: could not allocate buffer!")
261  return;
262  }
263 
264  m_cursor = m_buffer;
265 }
266 
267 
268 std::string WriterAsciiHepMC2::escape(const std::string& s) const
269 {
270  std::string ret;
271  ret.reserve(s.length()*2);
272  for ( std::string::const_iterator it = s.begin(); it != s.end(); ++it )
273  {
274  switch ( *it )
275  {
276  case '\\':
277  ret += "\\\\";
278  break;
279  case '\n':
280  ret += "\\|";
281  break;
282  default:
283  ret += *it;
284  }
285  }
286  return ret;
287 }
288 
289 void WriterAsciiHepMC2::write_vertex(ConstGenVertexPtr v)
290 {
291  std::vector<double> weights;
292  std::shared_ptr<VectorDoubleAttribute> weights_a = v->attribute<VectorDoubleAttribute>("weights");
293  if (weights_a) {
294  weights = weights_a->value();
295  } else {
296  weights.reserve(100);
297  for (int i = 0; i < 100; i++)
298  {
299  std::shared_ptr<DoubleAttribute> rs = v->attribute<DoubleAttribute>("weight"+std::to_string((long long unsigned int)i));
300  if (!rs) break;
301  weights.push_back(rs->value());
302  }
303  }
304  flush();
305  m_cursor += sprintf(m_cursor, "V %i %i", v->id(), v->status());
306  int orph = 0;
307  for (ConstGenParticlePtr p: v->particles_in())
308  {
309  if (!p->production_vertex()) orph++;
310  else
311  {
312  if (p->production_vertex()->id() == 0) orph++;
313  }
314  }
315  const FourVector &pos = v->position();
316  if (pos.is_zero())
317  {
318  m_cursor += sprintf(m_cursor, " 0 0 0 0");
319  }
320  else
321  {
322  m_cursor += sprintf(m_cursor, m_float_printf_specifier.c_str(), pos.x());
323  m_cursor += sprintf(m_cursor, m_float_printf_specifier.c_str(), pos.y());
324  m_cursor += sprintf(m_cursor, m_float_printf_specifier.c_str(), pos.z());
325  m_cursor += sprintf(m_cursor, m_float_printf_specifier.c_str(), pos.t());
326  }
327  m_cursor += sprintf(m_cursor, " %i %zu %zu", orph, v->particles_out().size(), weights.size());
328  flush();
329  for (size_t i = 0; i < weights.size(); i++) { m_cursor += sprintf(m_cursor, m_float_printf_specifier.c_str(), weights[i]); flush(); }
330  m_cursor += sprintf(m_cursor, "\n");
331  flush();
332 }
333 
334 
336 {
337  // The maximum size of single add to the buffer should not be larger than 256. This is a safe value as
338  // we will not allow precision larger than 24 anyway
339  if ( m_buffer + m_buffer_size < m_cursor + 512 )
340  {
341  std::ptrdiff_t length = m_cursor - m_buffer;
342  m_stream->write(m_buffer, length);
343  m_cursor = m_buffer;
344  }
345 }
346 
347 
349 {
350  std::ptrdiff_t length = m_cursor - m_buffer;
351  m_stream->write(m_buffer, length);
352  m_cursor = m_buffer;
353 }
354 
355 
357 
358 void WriterAsciiHepMC2::write_particle(ConstGenParticlePtr p, int /*second_field*/)
359 {
360  flush();
361  m_cursor += sprintf(m_cursor, "P %i", int(10001+m_particle_counter));
363  m_cursor += sprintf(m_cursor, " %i", p->pid() );
364  m_cursor += sprintf(m_cursor, m_float_printf_specifier.c_str(), p->momentum().px() );
365  m_cursor += sprintf(m_cursor, m_float_printf_specifier.c_str(), p->momentum().py());
366  m_cursor += sprintf(m_cursor, m_float_printf_specifier.c_str(), p->momentum().pz() );
367  m_cursor += sprintf(m_cursor, m_float_printf_specifier.c_str(), p->momentum().e() );
368  m_cursor += sprintf(m_cursor, m_float_printf_specifier.c_str(), p->generated_mass() );
369  m_cursor += sprintf(m_cursor, " %i", p->status() );
370  flush();
371  int ev = 0;
372  if (p->end_vertex())
373  if (p->end_vertex()->id() != 0)
374  ev = p->end_vertex()->id();
375 
376  std::shared_ptr<DoubleAttribute> A_theta = p->attribute<DoubleAttribute>("theta");
377  std:: shared_ptr<DoubleAttribute> A_phi = p->attribute<DoubleAttribute>("phi");
378  if (A_theta) m_cursor += sprintf(m_cursor, m_float_printf_specifier.c_str(), A_theta->value());
379  else m_cursor += sprintf(m_cursor, " 0");
380  if (A_phi) m_cursor += sprintf(m_cursor, m_float_printf_specifier.c_str(), A_phi->value());
381  else m_cursor += sprintf(m_cursor, " 0");
382  m_cursor += sprintf(m_cursor, " %i", ev);
383  flush();
384  std::shared_ptr<VectorIntAttribute> A_flows = p->attribute<VectorIntAttribute>("flows");
385  if (A_flows)
386  {
387  std::vector<int> flowsv = A_flows->value();
388  std::string flowss = " " + std::to_string(flowsv.size());
389  for (size_t k = 0; k < flowsv.size(); k++) { flowss += ( " " + std::to_string(k+1) + " " + std::to_string(flowsv.at(k))); }
390  flowss += "\n";
391  write_string(flowss);
392  } else {
393  std::shared_ptr<IntAttribute> A_flow1 = p->attribute<IntAttribute>("flow1");
394  std::shared_ptr<IntAttribute> A_flow2 = p->attribute<IntAttribute>("flow2");
395  std::shared_ptr<IntAttribute> A_flow3 = p->attribute<IntAttribute>("flow3");
396  int flowsize = 0;
397  if (A_flow1) flowsize++;
398  if (A_flow2) flowsize++;
399  if (A_flow3) flowsize++;
400  std::string flowss = " " + std::to_string(flowsize);
401  if (A_flow1) flowss += ( " 1 " + std::to_string(A_flow1->value()));
402  if (A_flow2) flowss += ( " 2 " + std::to_string(A_flow2->value()));
403  if (A_flow3) flowss += ( " 3 " + std::to_string(A_flow3->value()));
404  flowss += "\n";
405  write_string(flowss);
406  }
407  flush();
408 }
409 
410 
411 inline void WriterAsciiHepMC2::write_string(const std::string &str)
412 {
413  // First let's check if string will fit into the buffer
414  if ( m_buffer + m_buffer_size > m_cursor + str.length() )
415  {
416  strncpy(m_cursor, str.data(), str.length());
417  m_cursor += str.length();
418  flush();
419  }
420  // If not, flush the buffer and write the string directly
421  else
422  {
423  forced_flush();
424  m_stream->write(str.data(), str.length());
425  }
426 }
427 
428 
430 {
431  std::ofstream* ofs = dynamic_cast<std::ofstream*>(m_stream);
432  if (ofs && !ofs->is_open()) return;
433  forced_flush();
434  (*m_stream) << "HepMC::IO_GenEvent-END_EVENT_LISTING" << std::endl << std::endl;
435  if (ofs) ofs->close();
436 }
437 bool WriterAsciiHepMC2::failed() { return (bool)m_file.rdstate(); }
438 
439 void WriterAsciiHepMC2::set_precision(const int& prec ) {
440  if (prec < 2 || prec > 24) return;
441  m_precision = prec;
442 }
443 
445  return m_precision;
446 }
447 
448 void WriterAsciiHepMC2::set_buffer_size(const size_t& size ) {
449  if (m_buffer) return;
450  if (size < 1024) return;
451  m_buffer_size = size;
452 }
453 } // 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:191
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:143
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:136
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:141
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:390
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:86
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:125
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.