22 #include "HepPDT/ParticleDataTable.hh"
28 static const char *
const s_protocolSep =
":";
34 const std::string&
name,
35 const IInterface*
parent ) :
37 m_ioBackend( nullptr )
45 "Name of the back-end we'll use to write out the HepMC::GenEvent."
46 "\nEx: ascii:hepmc.genevent.txt" );
51 "Input location of the McEventCollection to write out" );
53 declareInterface<IIOHepMcTool>(
this);
79 (*m_ioBackend) <<
"# epsilon<float> = "
80 << std::numeric_limits<float>::epsilon()
82 <<
"# dbl precision = "
83 << std::numeric_limits<double>::digits10 + 1
85 <<
"# flt precision = "
86 << std::numeric_limits<float>::digits10 + 1
89 return StatusCode::SUCCESS;
95 return StatusCode::SUCCESS;
104 return StatusCode::FAILURE;
107 if ( mcEvts->
empty() ) {
109 return StatusCode::FAILURE;
112 const HepMC::GenEvent *
evt = mcEvts->
front();
114 ATH_MSG_ERROR(
"Retrieved NULL pointer to HepMC::GenEvent !!");
115 return StatusCode::FAILURE;
127 std::ostringstream
out;
131 out.precision( std::numeric_limits<double>::digits10 + 1 );
133 out.setf(std::ios::dec,std::ios::basefield);
134 out.setf(std::ios::scientific,std::ios::floatfield);
138 long evt_vertices_size=
evt->vertices().size();
139 std::shared_ptr<HepMC3::DoubleAttribute> A_alphaQCD=
evt->attribute<HepMC3::DoubleAttribute>(
"alphaQCD");
140 double evt_alphaQCD=(A_alphaQCD?(A_alphaQCD->value()):0.0);
141 std::shared_ptr<HepMC3::DoubleAttribute> A_alphaQED=
evt->attribute<HepMC3::DoubleAttribute>(
"alphaQED");
142 double evt_alphaQED=(A_alphaQED?(A_alphaQED->value()):0.0);
143 std::shared_ptr<HepMC3::DoubleAttribute> A_event_scale=
evt->attribute<HepMC3::DoubleAttribute>(
"event_scale");
144 double evt_event_scale=(A_event_scale?(A_event_scale->value()):0.0);
145 std::shared_ptr<HepMC3::VectorLongIntAttribute> A_random_states=
evt->attribute<HepMC3::VectorLongIntAttribute>(
"random_states");
146 std::vector<long int> random_states=(A_random_states?(A_random_states->value()):std::vector<long int>());
147 long random_states_size=random_states.size();
149 long evt_vertices_size=
evt->vertices_size();
150 double evt_alphaQCD=
evt->alphaQCD();
151 double evt_alphaQED=
evt->alphaQED();
152 double evt_event_scale=
evt->event_scale();
153 std::vector<long int> random_states=
evt->random_states();
154 long random_states_size=random_states.size();
156 out <<
"# -- GenEvent -->\n";
157 out <<
"#" <<
evt->event_number()
158 <<
" " << evt_event_scale
159 <<
" " << evt_alphaQCD
160 <<
" " << evt_alphaQED
163 <<
" " << evt_vertices_size
164 <<
" " << random_states_size
167 std::copy( random_states.begin(), random_states.end(), std::ostream_iterator<long int>(
out,
" ") );
168 out <<
evt->weights().size() <<
"\n";
171 std::copy(
evt->weights().begin(),
evt->weights().end(), std::ostream_iterator<double>(
out,
" ") );
174 out <<
"#-- particles --\n";
175 for (
const auto&
p: *
evt) {
179 const HepMC::FourVector
mom =
p->momentum();
180 std::ostringstream buf;
181 buf.precision( std::numeric_limits<float>::digits10 + 1 );
182 buf.setf(std::ios::dec,std::ios::basefield);
183 buf.setf(std::ios::scientific,std::ios::floatfield);
185 const float px =
static_cast<float>(
mom.px());
186 const float py =
static_cast<float>(
mom.py());
187 const float pz =
static_cast<float>(
mom.pz());
188 const float m =
static_cast<float>(
mom.m());
190 buf <<
px <<
" " <<
py <<
" " <<
pz <<
" " <<
e <<
" " <<
m <<
"\n";
194 out <<
"# "<<
p->status()
195 <<
" " << pol.theta()
203 out <<
"#-- vertices -- \n";
205 for (
const auto&
v:
evt->vertices()) {
208 const HepMC::FourVector
pos =
v->position();
209 std::ostringstream buf;
210 buf.precision( std::numeric_limits<float>::digits10 + 1 );
211 buf.setf(std::ios::dec,std::ios::basefield);
212 buf.setf(std::ios::scientific,std::ios::floatfield);
214 buf <<
pos.x() <<
" " <<
pos.y() <<
" " <<
pos.z() <<
" " <<
pos.t() <<
"\n";
218 std::string svertexeights(
"1.0");
219 auto vertexeights=
v->attribute<HepMC3::VectorDoubleAttribute>(
"weights");
220 if (vertexeights) vertexeights->to_string(svertexeights);
221 out << svertexeights;
225 for ( HepMC::GenEvent::vertex_const_iterator
226 i =
evt->vertices_begin(),
227 iEnd =
evt->vertices_end();
230 const HepMC::GenVertex *
v = *
i;
232 out <<
"# " <<
v->barcode() <<
" " <<
v->id() <<
"\n";
234 const HepMC::FourVector
pos =
v->position();
235 std::ostringstream buf;
236 buf.precision( std::numeric_limits<float>::digits10 + 1 );
237 buf.setf(std::ios::dec,std::ios::basefield);
238 buf.setf(std::ios::scientific,std::ios::floatfield);
240 buf <<
pos.x() <<
" " <<
pos.y() <<
" " <<
pos.z() <<
" " <<
pos.t() <<
"\n";
244 std::copy(
v->weights().begin(),
v->weights().end(), std::ostream_iterator<double>(
out,
" ") );
249 out <<
"#<-- GenEvent --\n";
252 return StatusCode::SUCCESS;
262 std::string protocol =
"ascii";
263 std::string
fileName =
"hepmc.genevent.txt";
274 std::string::size_type protocolPos =
url.find(s_protocolSep);
276 if ( std::string::npos != protocolPos ) {
277 protocol =
url.substr( 0, protocolPos );
278 fileName =
url.substr( protocolPos+1, std::string::npos );
285 std::transform( protocol.begin(), protocol.end(), protocol.begin(), [](
unsigned char c){ return std::tolower(c); } );
286 if (
"ascii" == protocol ) {
290 ATH_MSG_WARNING(
"UNKNOWN protocol [" << protocol <<
"] !!" <<
endmsg <<
"Will use [ascii] instead...");