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, std::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(std::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, std::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(std::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  std::shared_ptr<DoubleAttribute> A_event_scale=evt.attribute<DoubleAttribute>("event_scale");
83  std::shared_ptr<DoubleAttribute> A_alphaQED=evt.attribute<DoubleAttribute>("alphaQED");
84  std::shared_ptr<DoubleAttribute> A_alphaQCD=evt.attribute<DoubleAttribute>("alphaQCD");
85  std::shared_ptr<IntAttribute> A_signal_process_id=evt.attribute<IntAttribute>("signal_process_id");
86  std::shared_ptr<IntAttribute> A_mpi=evt.attribute<IntAttribute>("mpi");
87  std::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  std::shared_ptr<IntAttribute> rs=evt.attribute<IntAttribute>("random_states"+std::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  m_cursor += sprintf(m_cursor, " %lu",evt.weights().size());
141  if ( evt.weights().size() )
142  {
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  m_cursor += sprintf(m_cursor, "N %lu",evt.weights().size());
148  std::vector<std::string> names = run_info()->weight_names();
149  for (size_t q=0; q<evt.weights().size(); q++)
150  {
151  if (q<names.size())
152  m_cursor += sprintf(m_cursor, " \"%s\"",names[q].c_str());
153  else
154  m_cursor += sprintf(m_cursor, " \"%i\"",(int)q);
155  flush();
156  }
157  }
158  m_cursor += sprintf(m_cursor, "\n");
159  flush();
160  // Write units
161  m_cursor += sprintf(m_cursor, "U %s %s\n", Units::name(evt.momentum_unit()).c_str(), Units::name(evt.length_unit()).c_str());
162  flush();
163  std::shared_ptr<GenCrossSection> cs = evt.attribute<GenCrossSection>("GenCrossSection");
164  if(cs) {m_cursor += sprintf(m_cursor, "C %.*e %.*e\n",m_precision, cs->xsec(),m_precision,cs->xsec_err()); flush(); }
165 
166 
167  // Write attributes
168  for ( auto vt1: evt.attributes() )
169  {
170  for ( auto vt2: vt1.second )
171  {
172 
173  std::string st;
174  bool status = vt2.second->to_string(st);
175 
176  if( !status )
177  {
178  HEPMC3_WARNING( "WriterAsciiHepMC2::write_event: problem serializing attribute: "<<vt1.first )
179  }
180  else
181  {
182  if (vt1.first=="GenPdfInfo")
183  {
184  m_cursor +=
185  sprintf(m_cursor, "F ");
186  flush();
187  write_string(escape(st));
188  m_cursor += sprintf(m_cursor, "\n");
189  flush();
190  }
191  }
192  }
193  }
195  for(ConstGenVertexPtr v: evt.vertices() )
196  {
197  int production_vertex = 0;
198  production_vertex=v->id();
199  write_vertex(v);
200  for(ConstGenParticlePtr p: v->particles_in())
201  {
202  if (!p->production_vertex()) write_particle( p, production_vertex );
203  else
204  {
205  if (p->production_vertex()->id()==0)write_particle( p, production_vertex );
206  }
207  }
208  for(ConstGenParticlePtr p: v->particles_out())
209  write_particle( p, production_vertex );
210  }
211 
212  // Flush rest of the buffer to file
213  forced_flush();
214 }
215 
216 
218 {
219  if ( m_buffer ) return;
220  while( m_buffer==nullptr && m_buffer_size >= 256 ) {
221  try {
222  m_buffer = new char[ m_buffer_size ]();
223  } catch (const std::bad_alloc& e) {
224  delete[] m_buffer;
225  m_buffer_size /= 2;
226  HEPMC3_WARNING( "WriterAsciiHepMC2::allocate_buffer:"<<e.what()<<" buffer size too large. Dividing by 2. New size: " << m_buffer_size )
227  }
228  }
229 
230  if ( !m_buffer )
231  {
232  HEPMC3_ERROR( "WriterAsciiHepMC2::allocate_buffer: could not allocate buffer!" )
233  return;
234  }
235 
236  m_cursor = m_buffer;
237 }
238 
239 
240 std::string WriterAsciiHepMC2::escape(const std::string& s) const
241 {
242  std:: string ret;
243  ret.reserve( s.length()*2 );
244  for ( std::string::const_iterator it = s.begin(); it != s.end(); ++it )
245  {
246  switch ( *it )
247  {
248  case '\\':
249  ret += "\\\\";
250  break;
251  case '\n':
252  ret += "\\|";
253  break;
254  default:
255  ret += *it;
256  }
257  }
258  return ret;
259 }
260 
261 void WriterAsciiHepMC2::write_vertex(ConstGenVertexPtr v)
262 {
263  std::vector<double> weights;
264  for (int i=0; i<100; i++)
265  {
266  std::shared_ptr<DoubleAttribute> rs=v->attribute<DoubleAttribute>("weight"+std::to_string((long long unsigned int)i));
267  if (!rs) break;
268  weights.push_back(rs->value());
269  }
270 
271  m_cursor += sprintf( m_cursor, "V %i %i",v->id(),v->status() );
272  flush();
273  int orph=0;
274  for(ConstGenParticlePtr p: v->particles_in())
275  {
276  if (!p->production_vertex()) orph++;
277  else
278  {
279  if (p->production_vertex()->id()==0)orph++;
280  }
281  }
282  const FourVector &pos = v->position();
283  if (pos.is_zero())
284  {
285  m_cursor += sprintf(m_cursor," 0 0 0 0");
286  }
287  else
288  {
289  m_cursor += sprintf(m_cursor," %.*e",m_precision,pos.x());
290  flush();
291  m_cursor += sprintf(m_cursor," %.*e", m_precision,pos.y());
292  flush();
293  m_cursor += sprintf(m_cursor," %.*e", m_precision,pos.z());
294  flush();
295  m_cursor += sprintf(m_cursor," %.*e", m_precision,pos.t());
296  flush();
297  }
298  m_cursor += sprintf(m_cursor," %i %lu %lu",orph,v->particles_out().size(),weights.size());
299  flush();
300  for (size_t i=0; i<weights.size(); i++) m_cursor += sprintf(m_cursor," %.*e", m_precision,weights[i]);
301  m_cursor += sprintf(m_cursor,"\n");
302  flush();
303 }
304 
305 
307 {
308  // The maximum size of single add to the buffer (other than by
309  // using WriterAsciiHepMC2::write) is 32 bytes. This is a safe value as
310  // we will not allow precision larger than 24 anyway
311  unsigned long length = m_cursor - m_buffer;
312  if ( m_buffer_size - length < 32 )
313  {
314  m_stream->write( m_buffer, length );
315  m_cursor = m_buffer;
316  }
317 }
318 
319 
321 {
322  // m_file.write( m_buffer, m_cursor-m_buffer );
323  m_stream->write( m_buffer, m_cursor - m_buffer );
324  m_cursor = m_buffer;
325 }
326 
327 
329 
330 void WriterAsciiHepMC2::write_particle(ConstGenParticlePtr p, int second_field)
331 {
332 
333  m_cursor += sprintf(m_cursor,"P %i",int(10001+m_particle_counter));
335  flush();
336  m_cursor += sprintf(m_cursor," %i", p->pid() );
337  flush();
338  m_cursor += sprintf(m_cursor," %.*e", m_precision,p->momentum().px() );
339  flush();
340  m_cursor += sprintf(m_cursor," %.*e", m_precision,p->momentum().py());
341  flush();
342  m_cursor += sprintf(m_cursor," %.*e", m_precision,p->momentum().pz() );
343  flush();
344  m_cursor += sprintf(m_cursor," %.*e", m_precision,p->momentum().e() );
345  flush();
346  m_cursor += sprintf(m_cursor," %.*e", m_precision,p->generated_mass() );
347  flush();
348  m_cursor += sprintf(m_cursor," %i", p->status() );
349  flush();
350  int ev=0;
351  if (p->end_vertex())
352  if (p->end_vertex()->id()!=0)
353  ev=p->end_vertex()->id();
354 
355  std::shared_ptr<DoubleAttribute> A_theta=p->attribute<DoubleAttribute>("theta");
356  std:: shared_ptr<DoubleAttribute> A_phi=p->attribute<DoubleAttribute>("phi");
357  if (A_theta) m_cursor += sprintf(m_cursor," %.*e", m_precision, A_theta->value());
358  else m_cursor += sprintf(m_cursor," 0");
359  flush();
360  if (A_phi) m_cursor += sprintf(m_cursor," %.*e", m_precision, A_phi->value());
361  else m_cursor += sprintf(m_cursor," 0");
362  flush();
363  m_cursor += sprintf(m_cursor," %i", ev );
364  flush();
365  std::shared_ptr<VectorIntAttribute> A_flows=p->attribute<VectorIntAttribute>("flows");
366  if (A_flows)
367  {
368  std::vector<int> flowsv=A_flows->value();
369  int flowsize=flowsv.size();
370  m_cursor += sprintf(m_cursor," %i", flowsize);
371  for (size_t k=0; k<flowsv.size(); k++) m_cursor += sprintf(m_cursor," %lu %i", k+1, flowsv.at(k));
372  m_cursor += sprintf(m_cursor,"\n");
373  flush();
374  } else {
375  std::shared_ptr<IntAttribute> A_flow1=p->attribute<IntAttribute>("flow1");
376  std::shared_ptr<IntAttribute> A_flow2=p->attribute<IntAttribute>("flow2");
377  std::shared_ptr<IntAttribute> A_flow3=p->attribute<IntAttribute>("flow3");
378  int flowsize=0;
379  if (A_flow1) flowsize++;
380  if (A_flow2) flowsize++;
381  if (A_flow3) flowsize++;
382  m_cursor += sprintf(m_cursor," %i", flowsize);
383  if (A_flow1) m_cursor += sprintf(m_cursor," 1 %i", A_flow1->value());
384  if (A_flow2) m_cursor += sprintf(m_cursor," 2 %i", A_flow2->value());
385  if (A_flow3) m_cursor += sprintf(m_cursor," 3 %i", A_flow3->value());
386  m_cursor += sprintf(m_cursor,"\n");
387  flush();
388  }
389 }
390 
391 
392 inline void WriterAsciiHepMC2::write_string( const std::string &str )
393 {
394 
395  // First let's check if string will fit into the buffer
396  unsigned long length = m_cursor-m_buffer;
397 
398  if ( m_buffer_size - length > str.length() )
399  {
400  strncpy(m_cursor,str.data(),str.length());
401  m_cursor += str.length();
402  flush();
403  }
404  // If not, flush the buffer and write the string directly
405  else
406  {
407  forced_flush();
408  m_stream->write( str.data(), str.length() );
409  }
410 }
411 
412 
414 {
415  std::ofstream* ofs = dynamic_cast<std::ofstream*>(m_stream);
416  if (ofs && !ofs->is_open()) return;
417  forced_flush();
418  (*m_stream) << "HepMC::IO_GenEvent-END_EVENT_LISTING" << std::endl << std::endl;
419  if (ofs) ofs->close();
420 }
421 bool WriterAsciiHepMC2::failed() { return (bool)m_file.rdstate(); }
422 
423 void WriterAsciiHepMC2::set_precision(const int& prec ) {
424  if (prec < 2 || prec > 24) return;
425  m_precision = prec;
426 }
427 
429  return m_precision;
430 }
431 
432 void WriterAsciiHepMC2::set_buffer_size(const size_t& size ) {
433  if (m_buffer) return;
434  if (size < 256) return;
435  m_buffer_size = size;
436 }
437 } // namespace HepMC3
const std::vector< ConstGenVertexPtr > & vertices() const
Get list of vertices (const)
Definition: GenEvent.cc:43
#define HEPMC3_WARNING(MESSAGE)
Macro for printing HEPMC3_HEPMC3_WARNING messages.
Definition: Errors.h:26
double t() const
Time component of position/displacement.
Definition: FourVector.h:101
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:192
Attribute that holds a vector of integers of type int.
Definition: Attribute.h:1002
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:142
double z() const
z-component of position/displacement
Definition: FourVector.h:94
double x() const
x-component of position/displacement
Definition: FourVector.h:80
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:135
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:1028
unsigned long m_buffer_size
Buffer size.
unsigned long m_particle_counter
Used to set bar codes.
Definition: pytypes.h:935
const Units::MomentumUnit & momentum_unit() const
Get momentum unit.
Definition: GenEvent.h:140
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: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::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:389
Definition of class WriterAsciiHepMC2.
void set_run_info(std::shared_ptr< GenRunInfo > run)
Set the global GenRunInfo object.
Definition: Writer.h:42
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:87
#define HEPMC3_ERROR(MESSAGE)
Macro for printing error messages.
Definition: Errors.h:23
const std::vector< double > & weights() const
Get event weight values as a vector.
Definition: GenEvent.h:86
int m_precision
Output precision.
Definition of class GenEvent.
void set_precision(const int &prec)
Set output precision.
std::map< std::string, std::map< int, std::shared_ptr< Attribute > > > attributes() const
Get a copy of the list of attributes.
Definition: GenEvent.h:237
std::shared_ptr< GenRunInfo > run_info() const
Get a pointer to the the GenRunInfo object.
Definition: GenEvent.h:124
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.