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