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
19namespace HepMC3
20{
21
22
23WriterAsciiHepMC2::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
46WriterAsciiHepMC2::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 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 std::vector<std::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 std::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 std::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
236}
237
238
239std::string WriterAsciiHepMC2::escape(const std::string& s) const
240{
241 std:: string ret;
242 ret.reserve( s.length()*2 );
243 for ( std::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
260void WriterAsciiHepMC2::write_vertex(ConstGenVertexPtr v)
261{
262 std::vector<double> weights;
263 for (int i=0; i<100; i++)
264 {
265 std::shared_ptr<DoubleAttribute> rs=v->attribute<DoubleAttribute>("weight"+std::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 );
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 );
324}
325
326
328
329void 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 std::shared_ptr<DoubleAttribute> A_theta=p->attribute<DoubleAttribute>("theta");
355 std:: 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 std::shared_ptr<IntAttribute> A_flow1=p->attribute<IntAttribute>("flow1");
366 std::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
378inline void WriterAsciiHepMC2::write_string( const std::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" << std::endl << std::endl;
405 if (ofs) ofs->close();
406}
407bool WriterAsciiHepMC2::failed() { return (bool)m_file.rdstate(); }
408
409void 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
418void 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
#define HEPMC3_WARNING(MESSAGE)
Macro for printing HEPMC3_HEPMC3_WARNING messages.
Definition: Errors.h:26
#define HEPMC3_ERROR(MESSAGE)
Macro for printing error messages.
Definition: Errors.h:23
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:243
Generic 4-vector.
Definition: FourVector.h:35
double t() const
Time component of position/displacement.
Definition: FourVector.h:101
bool is_zero() const
Check if the length of this vertex is zero.
Definition: FourVector.h:192
double x() const
x-component of position/displacement
Definition: FourVector.h:80
double y() const
y-component of position/displacement
Definition: FourVector.h:87
double z() const
z-component of position/displacement
Definition: FourVector.h:94
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:389
int event_number() const
Get event number.
Definition: GenEvent.h:135
std::shared_ptr< GenRunInfo > run_info() const
Get a pointer to the the GenRunInfo object.
Definition: GenEvent.h:124
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:140
const Units::LengthUnit & length_unit() const
Get length unit.
Definition: GenEvent.h:142
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
const std::vector< double > & weights() const
Get event weight values as a vector.
Definition: GenEvent.h:86
Attribute that holds an Integer implemented as an int.
Definition: Attribute.h:159
static std::string name(MomentumUnit u)
Get name of momentum unit.
Definition: Units.h:56
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.
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
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