ATLAS Offline Software
HepMcFloatWriterTool.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
4 
6 // HepMcFloatWriterTool.cxx
7 // Implementation file for class HepMcFloatWriterTool
8 // Author: S.Binet<binet@cern.ch>
10 
11 
12 // STL includes
13 #include <algorithm>
14 #include <cctype>
15 #include <fstream>
16 #include <sstream>
17 #include <limits>
18 
19 // FrameWork includes
20 
22 #include "HepPDT/ParticleDataTable.hh"
23 
24 // McParticleTools includes
25 #include "HepMcFloatWriterTool.h"
26 #include "AtlasHepMC/Flow.h"
28 static const char * const s_protocolSep = ":";
29 
34  const std::string& name,
35  const IInterface* parent ) :
36  base_class( type, name, parent ),
37  m_ioBackend( nullptr )
38 {
39  //
40  // Property declaration
41  //
42 
43  declareProperty( "Output",
44  m_ioBackendURL = "ascii:hepmc.genevent.txt",
45  "Name of the back-end we'll use to write out the HepMC::GenEvent."
46  "\nEx: ascii:hepmc.genevent.txt" );
47  m_ioBackendURL.declareUpdateHandler( &HepMcFloatWriterTool::setupBackend,this );
48 
49  declareProperty( "McEvents",
50  m_mcEventsName = "GEN_EVENT",
51  "Input location of the McEventCollection to write out" );
52 }
53 
57 {
58  ATH_MSG_DEBUG("Calling destructor");
59 
60  if ( m_ioBackend ) {
61  delete m_ioBackend;
62  m_ioBackend = nullptr;
63  }
64 }
65 
69 {
70  ATH_MSG_INFO("Initializing " << name() << "...");
71 
72  // setup backend
73  if ( nullptr == m_ioBackend ) {
75  }
76 
77  (*m_ioBackend) << "# epsilon<float> = "
78  << std::numeric_limits<float>::epsilon()
79  << "\n"
80  << "# dbl precision = "
81  << std::numeric_limits<double>::digits10 + 1
82  << "\n"
83  << "# flt precision = "
84  << std::numeric_limits<float>::digits10 + 1
85  << "\n";
86 
87  return StatusCode::SUCCESS;
88 }
89 
91 {
92  ATH_MSG_INFO("Finalizing " << name() << "...");
93  return StatusCode::SUCCESS;
94 }
95 
97 {
98  // retrieve the McEventCollection
99  const McEventCollection * mcEvts = nullptr;
100  if ( evtStore()->retrieve( mcEvts, m_mcEventsName ).isFailure() || nullptr == mcEvts ) {
101  ATH_MSG_ERROR("Could not retrieve a McEventCollection at [" << m_mcEventsName << "] !!");
102  return StatusCode::FAILURE;
103  }
104 
105  if ( mcEvts->empty() ) {
106  ATH_MSG_WARNING("McEventCollection at [" << m_mcEventsName << "] is EMPTY !!");
107  return StatusCode::FAILURE;
108  }
109 
110  const HepMC::GenEvent * evt = mcEvts->front();
111  if ( !evt ) {
112  ATH_MSG_ERROR("Retrieved NULL pointer to HepMC::GenEvent !!");
113  return StatusCode::FAILURE;
114  }
115 
116  return write(evt);
117 }
118 
122 
123 StatusCode HepMcFloatWriterTool::write( const HepMC::GenEvent* evt )
124 {
125  std::ostringstream out;
126 
127  // precision 8 (# digits following decimal point) is the minimum that
128  // will capture the full information stored in a float
129  out.precision( std::numeric_limits<double>::digits10 + 1 );
130  // we use decimal to store integers, because it is smaller than hex!
131  out.setf(std::ios::dec,std::ios::basefield);
132  out.setf(std::ios::scientific,std::ios::floatfield);
133 
134 
135 #ifdef HEPMC3
136  long evt_vertices_size=evt->vertices().size();
137  std::shared_ptr<HepMC3::DoubleAttribute> A_alphaQCD=evt->attribute<HepMC3::DoubleAttribute>("alphaQCD");
138  double evt_alphaQCD=(A_alphaQCD?(A_alphaQCD->value()):0.0);
139  std::shared_ptr<HepMC3::DoubleAttribute> A_alphaQED=evt->attribute<HepMC3::DoubleAttribute>("alphaQED");
140  double evt_alphaQED=(A_alphaQED?(A_alphaQED->value()):0.0);
141  std::shared_ptr<HepMC3::DoubleAttribute> A_event_scale=evt->attribute<HepMC3::DoubleAttribute>("event_scale");
142  double evt_event_scale=(A_event_scale?(A_event_scale->value()):0.0);
143  std::shared_ptr<HepMC3::VectorLongIntAttribute> A_random_states=evt->attribute<HepMC3::VectorLongIntAttribute>("random_states");
144  std::vector<long int> random_states=(A_random_states?(A_random_states->value()):std::vector<long int>());
145  long random_states_size=random_states.size();
146 #else
147  long evt_vertices_size=evt->vertices_size();
148  double evt_alphaQCD=evt->alphaQCD();
149  double evt_alphaQED=evt->alphaQED();
150  double evt_event_scale=evt->event_scale();
151  std::vector<long int> random_states=evt->random_states();
152  long random_states_size=random_states.size();
153 #endif
154  out << "# -- GenEvent -->\n";
155  out << "#" << evt->event_number()
156  << " " << evt_event_scale
157  << " " << evt_alphaQCD
158  << " " << evt_alphaQED
159  << " " << HepMC::signal_process_id(evt)
161  << " " << evt_vertices_size
162  << " " << random_states_size
163  << "\n";
164  out << "#";
165  std::copy( random_states.begin(), random_states.end(), std::ostream_iterator<long int>(out, " ") );
166  out << evt->weights().size() << "\n";
167 
168  out << "#";
169  std::copy( evt->weights().begin(), evt->weights().end(), std::ostream_iterator<double>(out, " ") );
170  out << '\n';
171 
172  out << "#-- particles --\n";
173  for (const auto& p: *evt) {
174  if ( p ) {
175  out << "# " << HepMC::barcode(p) << " " << p->pdg_id() << "\n";
176 
177  const HepMC::FourVector mom = p->momentum();
178  std::ostringstream buf;
179  buf.precision( std::numeric_limits<float>::digits10 + 1 );
180  buf.setf(std::ios::dec,std::ios::basefield);
181  buf.setf(std::ios::scientific,std::ios::floatfield);
182 
183  const float px = static_cast<float>(mom.px());
184  const float py = static_cast<float>(mom.py());
185  const float pz = static_cast<float>(mom.pz());
186  const float m = static_cast<float>(mom.m());
187  const float e = static_cast<float>(std::sqrt( std::pow( px, 2 ) + std::pow( py, 2 ) + std::pow( pz, 2 ) + std::pow( m, 2 ) ) );
188  buf << px << " " << py << " " << pz << " " << e << " " << m << "\n";
189 
190  out << buf.str();
191  auto pol=HepMC::polarization(p);
192  out << "# "<< p->status()
193  << " " << pol.theta()
194  << " " << pol.phi()
195  << " " << ( p->end_vertex() ? HepMC::barcode(p->end_vertex()) : 0 )
196  << " " << HepMC::flow(p)
197  << "\n";
198  }
199  }
200 
201  out << "#-- vertices -- \n";
202 #ifdef HEPMC3
203  for (const auto& v: evt->vertices()) {
204  if ( v ) {
205  out << "# " << HepMC::barcode(v) << " " << v->status() << "\n";
206  const HepMC::FourVector pos = v->position();
207  std::ostringstream buf;
208  buf.precision( std::numeric_limits<float>::digits10 + 1 );
209  buf.setf(std::ios::dec,std::ios::basefield);
210  buf.setf(std::ios::scientific,std::ios::floatfield);
211 
212  buf << pos.x() << " " << pos.y() << " " << pos.z() << " " << pos.t() << "\n";
213 
214  out << buf.str();
215  out << "#";
216  std::string svertexeights("1.0");
217  auto vertexeights=v->attribute<HepMC3::VectorDoubleAttribute>("weights");
218  if (vertexeights) vertexeights->to_string(svertexeights);
219  out << svertexeights;
220  out << '\n';
221  }
222 #else
223  for ( HepMC::GenEvent::vertex_const_iterator
224  i = evt->vertices_begin(),
225  iEnd = evt->vertices_end();
226  i != iEnd;
227  ++i ) {
228  const HepMC::GenVertex * v = *i;
229  if ( v ) {
230  out << "# " << v->barcode() << " " << v->id() << "\n";
231 
232  const HepMC::FourVector pos = v->position();
233  std::ostringstream buf;
234  buf.precision( std::numeric_limits<float>::digits10 + 1 );
235  buf.setf(std::ios::dec,std::ios::basefield);
236  buf.setf(std::ios::scientific,std::ios::floatfield);
237 
238  buf << pos.x() << " " << pos.y() << " " << pos.z() << " " << pos.t() << "\n";
239 
240  out << buf.str();
241  out << "#";
242  std::copy( v->weights().begin(), v->weights().end(), std::ostream_iterator<double>(out, " ") );
243  out << '\n';
244  }
245 #endif
246  }
247  out << "#<-- GenEvent --\n";
248 
249  (*m_ioBackend) << out.str() << std::flush;
250  return StatusCode::SUCCESS;
251 }
252 
256 
257 void HepMcFloatWriterTool::setupBackend( Gaudi::Details::PropertyBase& /*prop*/ )
258 {
259  // defaults
260  std::string protocol = "ascii";
261  std::string fileName = "hepmc.genevent.txt";
262 
263  // reset internal state
264  if ( m_ioBackend ) {
265  delete m_ioBackend;
266  m_ioBackend = nullptr;
267  }
268 
269  // caching URL
270  const std::string& url = m_ioBackendURL.value();
271 
272  std::string::size_type protocolPos = url.find(s_protocolSep);
273 
274  if ( std::string::npos != protocolPos ) {
275  protocol = url.substr( 0, protocolPos );
276  fileName = url.substr( protocolPos+1, std::string::npos );
277  } else {
278  //protocol = "ascii";
279  fileName = url;
280  }
281 
282  // get the protocol name in lower cases
283  std::transform( protocol.begin(), protocol.end(), protocol.begin(), [](unsigned char c){ return std::tolower(c); } );
284  if ( "ascii" == protocol ) {
285  m_ioBackend = new std::ofstream( fileName.c_str(), std::ios::out | std::ios::trunc );
286 
287  } else {
288  ATH_MSG_WARNING("UNKNOWN protocol [" << protocol << "] !!" << endmsg << "Will use [ascii] instead...");
289  protocol = "ascii";
290  m_ioBackend = new std::ofstream( fileName.c_str(), std::ios::out | std::ios::trunc );
291  }
292  ATH_MSG_DEBUG("Using protocol [" << protocol << "] and write to ["<< fileName << "]");
293 }
python.PyKernel.retrieve
def retrieve(aClass, aKey=None)
Definition: PyKernel.py:110
Flow.h
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
HepMcFloatWriterTool::initialize
StatusCode initialize()
Athena Algorithm's Hooks.
Definition: HepMcFloatWriterTool.cxx:68
test_pyathena.px
px
Definition: test_pyathena.py:18
HepMcFloatWriterTool::HepMcFloatWriterTool
HepMcFloatWriterTool()
Default constructor:
HepMC::polarization
Polarization polarization(const T &a)
Definition: Polarization.h:47
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
FullCPAlgorithmsTest_eljob.flush
flush
Definition: FullCPAlgorithmsTest_eljob.py:194
HepMC::signal_process_id
int signal_process_id(const GenEvent &e)
Definition: GenEvent.h:635
python.AthDsoLogger.out
out
Definition: AthDsoLogger.py:70
LArG4FSStartPointFilter.evt
evt
Definition: LArG4FSStartPointFilter.py:42
physics_parameters.url
string url
Definition: physics_parameters.py:27
python.CaloAddPedShiftConfig.type
type
Definition: CaloAddPedShiftConfig.py:42
HepMcFloatWriterTool::m_mcEventsName
StringProperty m_mcEventsName
Location of the McEventCollection to be written out If there is more than 1 HepMC::GenEvent in the Mc...
Definition: HepMcFloatWriterTool.h:83
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:209
ParticleGun_EoverP_Config.mom
mom
Definition: ParticleGun_EoverP_Config.py:63
HepMcFloatWriterTool::m_ioBackend
std::ostream * m_ioBackend
Abstract base class for the back-end.
Definition: HepMcFloatWriterTool.h:87
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
McEventCollection.h
HepMC::flow
int flow(const T &a, int i)
Definition: Flow.h:51
lumiFormat.i
int i
Definition: lumiFormat.py:85
endmsg
#define endmsg
Definition: AnalysisConfig_Ntuple.cxx:63
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
HepMC::barcode
int barcode(const T *p)
Definition: Barcode.h:16
HepMcFloatWriterTool::finalize
StatusCode finalize()
Definition: HepMcFloatWriterTool.cxx:90
Amg::transform
Amg::Vector3D transform(Amg::Vector3D &v, Amg::Transform3D &tr)
Transform a point from a Trasformation3D.
Definition: GeoPrimitivesHelpers.h:156
DataVector::front
const T * front() const
Access the first element in the collection as an rvalue.
Amg::pz
@ pz
Definition: GeoPrimitives.h:40
test_pyathena.parent
parent
Definition: test_pyathena.py:15
McEventCollection
This defines the McEventCollection, which is really just an ObjectVector of McEvent objects.
Definition: McEventCollection.h:32
HepMcFloatWriterTool::setupBackend
void setupBackend(Gaudi::Details::PropertyBase &ioBackendURL)
Method to configure the back-end to write out the HepMC::GenEvent.
Definition: HepMcFloatWriterTool.cxx:257
Amg::py
@ py
Definition: GeoPrimitives.h:39
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
HepMcFloatWriterTool::~HepMcFloatWriterTool
virtual ~HepMcFloatWriterTool()
Destructor:
Definition: HepMcFloatWriterTool.cxx:56
HepMcFloatWriterTool::write
StatusCode write(const HepMC::GenEvent *evt)
Process the HepMC::GenEvent through the I/O backend.
Definition: HepMcFloatWriterTool.cxx:123
python.LumiBlobConversion.pos
pos
Definition: LumiBlobConversion.py:16
python.PyAthena.v
v
Definition: PyAthena.py:154
HepMcFloatWriterTool::m_ioBackendURL
StringProperty m_ioBackendURL
URL of the I/O back-end (only "ASCII" for now...) glued with the name of the output file name.
Definition: HepMcFloatWriterTool.h:77
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
HepMcFloatWriterTool.h
calibdata.copy
bool copy
Definition: calibdata.py:26
jobOptions.fileName
fileName
Definition: jobOptions.SuperChic_ALP2.py:39
Polarization.h
pow
constexpr int pow(int base, int exp) noexcept
Definition: ap_fixedTest.cxx:15
python.compressB64.c
def c
Definition: compressB64.py:93
DataVector::empty
bool empty() const noexcept
Returns true if the collection is empty.
HepMcFloatWriterTool::execute
StatusCode execute()
Definition: HepMcFloatWriterTool.cxx:96
python.SystemOfUnits.m
float m
Definition: SystemOfUnits.py:106
HepMC::signal_process_vertex
GenVertex * signal_process_vertex(const GenEvent *e)
Definition: GenEvent.h:625