Fork me on GitHub

source: git/external/HepMC3/GenEvent.h@ 302624f

Last change on this file since 302624f was 95a917c, checked in by Pavel Demin <pavel.demin@…>, 4 years ago

add HepMC3 library

  • Property mode set to 100644
File size: 14.4 KB
Line 
1// -*- C++ -*-
2//
3// This file is part of HepMC
4// Copyright (C) 2014-2020 The HepMC collaboration (see AUTHORS for details)
5//
6///
7/// @file GenEvent.h
8/// @brief Definition of \b class GenEvent
9///
10#ifndef HEPMC3_GENEVENT_H
11#define HEPMC3_GENEVENT_H
12
13#include "HepMC3/Units.h"
14#include "HepMC3/GenParticle_fwd.h"
15#include "HepMC3/GenVertex_fwd.h"
16#include "HepMC3/GenPdfInfo_fwd.h"
17#include "HepMC3/GenHeavyIon_fwd.h"
18#include "HepMC3/GenCrossSection_fwd.h"
19
20#if !defined(__CINT__)
21#include "HepMC3/GenHeavyIon.h"
22#include "HepMC3/GenPdfInfo.h"
23#include "HepMC3/GenCrossSection.h"
24#include "HepMC3/GenRunInfo.h"
25#include <mutex>
26#endif // __CINT__
27
28#ifdef HEPMC3_ROOTIO
29class TBuffer;
30#endif
31
32
33namespace HepMC3 {
34
35struct GenEventData;
36
37/// @brief Stores event-related information
38///
39/// Manages event-related information.
40/// Contains lists of GenParticle and GenVertex objects
41class GenEvent {
42
43public:
44
45 /// @brief Event constructor without a run
46 GenEvent(Units::MomentumUnit momentum_unit = Units::GEV,
47 Units::LengthUnit length_unit = Units::MM);
48
49#if !defined(__CINT__)
50
51 /// @brief Constructor with associated run
52 GenEvent(std::shared_ptr<GenRunInfo> run,
53 Units::MomentumUnit momentum_unit = Units::GEV,
54 Units::LengthUnit length_unit = Units::MM);
55
56 /// @brief Copy constructor
57 GenEvent(const GenEvent&);
58
59 /// @brief Destructor
60 ~GenEvent();
61
62 /// @brief Assignment operator
63 GenEvent& operator=(const GenEvent&);
64
65 /// @name Particle and vertex access
66 //@{
67
68 /// @brief Get list of particles (const)
69 const std::vector<ConstGenParticlePtr>& particles() const;
70 /// @brief Get list of vertices (const)
71 const std::vector<ConstGenVertexPtr>& vertices() const;
72
73
74 /// @brief Get/set list of particles (non-const)
75 const std::vector<GenParticlePtr>& particles() { return m_particles; }
76 /// @brief Get/set list of vertices (non-const)
77 const std::vector<GenVertexPtr>& vertices() { return m_vertices; }
78
79 //@}
80
81
82 /// @name Event weights
83 //@{
84
85 /// Get event weight values as a vector
86 const std::vector<double>& weights() const { return m_weights; }
87 /// Get event weights as a vector (non-const)
88 std::vector<double>& weights() { return m_weights; }
89 /// Get event weight accessed by index (or the canonical/first one if there is no argument)
90 /// @note It's the user's responsibility to ensure that the given index exists!
91 double weight(const unsigned long& index=0) const { return weights().at(index); }
92 /// Get event weight accessed by weight name
93 /// @note Requires there to be an attached GenRunInfo, otherwise will throw an exception
94 /// @note It's the user's responsibility to ensure that the given name exists!
95 double weight(const std::string& name) const {
96 if (!run_info()) throw std::runtime_error("GenEvent::weight(str): named access to event weights requires the event to have a GenRunInfo");
97 return weight(run_info()->weight_index(name));
98 }
99 /// Get event weight accessed by weight name
100 /// @note Requires there to be an attached GenRunInfo, otherwise will throw an exception
101 /// @note It's the user's responsibility to ensure that the given name exists!
102 double& weight(const std::string& name) {
103 if (!run_info()) throw std::runtime_error("GenEvent::weight(str): named access to event weights requires the event to have a GenRunInfo");
104 int pos=run_info()->weight_index(name);
105 if (pos<0) throw std::runtime_error("GenEvent::weight(str): no weight with given name in this run");
106 return m_weights[pos];
107 }
108 /// Get event weight names, if there are some
109 /// @note Requires there to be an attached GenRunInfo with registered weight names, otherwise will throw an exception
110 const std::vector<std::string>& weight_names() const {
111 if (!run_info()) throw std::runtime_error("GenEvent::weight_names(): access to event weight names requires the event to have a GenRunInfo");
112 const std::vector<std::string>& weightnames = run_info()->weight_names();
113 if (weightnames.empty()) throw std::runtime_error("GenEvent::weight_names(): no event weight names are registered for this run");
114 return weightnames;
115 }
116
117 //@}
118
119
120 /// @name Auxiliary info and event metadata
121 //@{
122
123 /// @brief Get a pointer to the the GenRunInfo object.
124 std::shared_ptr<GenRunInfo> run_info() const {
125 return m_run_info;
126 }
127 /// @brief Set the GenRunInfo object by smart pointer.
128 void set_run_info(std::shared_ptr<GenRunInfo> run) {
129 m_run_info = run;
130 if ( run && !run->weight_names().empty() )
131 m_weights.resize(run->weight_names().size(), 1.0);
132 }
133
134 /// @brief Get event number
135 int event_number() const { return m_event_number; }
136 /// @brief Set event number
137 void set_event_number(const int& num) { m_event_number = num; }
138
139 /// @brief Get momentum unit
140 const Units::MomentumUnit& momentum_unit() const { return m_momentum_unit; }
141 /// @brief Get length unit
142 const Units::LengthUnit& length_unit() const { return m_length_unit; }
143 /// @brief Change event units
144 /// Converts event from current units to new ones
145 void set_units( Units::MomentumUnit new_momentum_unit, Units::LengthUnit new_length_unit);
146
147 /// @brief Get heavy ion generator additional information
148 GenHeavyIonPtr heavy_ion() { return attribute<GenHeavyIon>("GenHeavyIon"); }
149 /// @brief Get heavy ion generator additional information (const version)
150 ConstGenHeavyIonPtr heavy_ion() const { return attribute<GenHeavyIon>("GenHeavyIon"); }
151 /// @brief Set heavy ion generator additional information
152 void set_heavy_ion(GenHeavyIonPtr hi) { add_attribute("GenHeavyIon",hi); }
153
154 /// @brief Get PDF information
155 GenPdfInfoPtr pdf_info() { return attribute<GenPdfInfo>("GenPdfInfo"); }
156 /// @brief Get PDF information (const version)
157 ConstGenPdfInfoPtr pdf_info() const { return attribute<GenPdfInfo>("GenPdfInfo"); }
158 /// @brief Set PDF information
159 void set_pdf_info(GenPdfInfoPtr pi) { add_attribute("GenPdfInfo",pi); }
160
161 /// @brief Get cross-section information
162 GenCrossSectionPtr cross_section() { return attribute<GenCrossSection>("GenCrossSection"); }
163 /// @brief Get cross-section information (const version)
164 ConstGenCrossSectionPtr cross_section() const { return attribute<GenCrossSection>("GenCrossSection"); }
165 /// @brief Set cross-section information
166 void set_cross_section(GenCrossSectionPtr cs) { add_attribute("GenCrossSection",cs); }
167
168 //@}
169
170
171 /// @name Event position
172 //@{
173
174 /// Vertex representing the overall event position
175 const FourVector& event_pos() const;
176
177 /// @brief Vector of beam particles
178 std::vector<ConstGenParticlePtr> beams() const;
179
180 /// @brief Vector of beam particles
181 const std::vector<GenParticlePtr> & beams();
182
183 /// @brief Shift position of all vertices in the event by @a delta
184 void shift_position_by( const FourVector & delta );
185
186 /// @brief Shift position of all vertices in the event to @a op
187 void shift_position_to( const FourVector & newpos ) {
188 const FourVector delta = newpos - event_pos();
189 shift_position_by(delta);
190 }
191
192 /// @brief Boost event using x,y,z components of @a v as velocities
193 bool boost( const FourVector& v );
194 /// @brief Rotate event using x,y,z components of @a v as rotation angles
195 bool rotate( const FourVector& v );
196 /// @brief Change sign of @a axis
197 bool reflect(const int axis);
198
199 //@}
200
201
202 /// @name Additional attributes
203 //@{
204 /// @brief Add event attribute to event
205 ///
206 /// This will overwrite existing attribute if an attribute
207 /// with the same name is present
208 void add_attribute(const std::string &name, const std::shared_ptr<Attribute> &att, const int& id = 0) {
209 ///Disallow empty strings
210 if (name.length()==0) return;
211 if (!att) return;
212 std::lock_guard<std::recursive_mutex> lock(m_lock_attributes);
213 if (m_attributes.count(name)==0) m_attributes[name]=std::map<int, std::shared_ptr<Attribute> >();
214 m_attributes[name][id] = att;
215 att->m_event = this;
216 if ( id > 0 && id <= int(particles().size()) )
217 att->m_particle = particles()[id - 1];
218 if ( id < 0 && -id <= int(vertices().size()) )
219 att->m_vertex = vertices()[-id - 1];
220 }
221
222 /// @brief Remove attribute
223 void remove_attribute(const std::string &name, const int& id = 0);
224
225 /// @brief Get attribute of type T
226 template<class T>
227 std::shared_ptr<T> attribute(const std::string &name, const int& id = 0) const;
228
229 /// @brief Get attribute of any type as string
230 std::string attribute_as_string(const std::string &name, const int& id = 0) const;
231
232 /// @brief Get list of attribute names
233 std::vector<std::string> attribute_names( const int& id = 0) const;
234
235 /// @brief Get a copy of the list of attributes
236 /// @note To avoid thread issues, this is returns a copy. Better solution may be needed.
237 std::map< std::string, std::map<int, std::shared_ptr<Attribute> > > attributes() const {
238 std::lock_guard<std::recursive_mutex> lock(m_lock_attributes);
239 return m_attributes;
240 }
241
242 //@}
243
244
245 /// @name Particle and vertex modification
246 //@{
247
248 /// @brief Add particle
249 void add_particle( GenParticlePtr p );
250
251 /// @brief Add vertex
252 void add_vertex( GenVertexPtr v );
253
254 /// @brief Remove particle from the event
255 ///
256 /// This function will remove whole sub-tree starting from this particle
257 /// if it is the only incoming particle of this vertex.
258 /// It will also production vertex of this particle if this vertex
259 /// has no more outgoing particles
260 void remove_particle( GenParticlePtr v );
261
262 /// @brief Remove a set of particles
263 ///
264 /// This function follows rules of GenEvent::remove_particle to remove
265 /// a list of particles from the event.
266 void remove_particles( std::vector<GenParticlePtr> v );
267
268 /// @brief Remove vertex from the event
269 ///
270 /// This will remove all sub-trees of all outgoing particles of this vertex
271 void remove_vertex( GenVertexPtr v );
272
273 /// @brief Add whole tree in topological order
274 ///
275 /// This function will find the beam particles (particles
276 /// that have no production vertices or their production vertices
277 /// have no particles) and will add the whole decay tree starting from
278 /// these particles.
279 ///
280 /// @note Any particles on this list that do not belong to the tree
281 /// will be ignored.
282 void add_tree( const std::vector<GenParticlePtr> &particles );
283
284 /// @brief Reserve memory for particles and vertices
285 ///
286 /// Helps optimize event creation when size of the event is known beforehand
287 void reserve(const size_t& particles, const size_t& vertices = 0);
288
289 /// @brief Remove contents of this event
290 void clear();
291
292 //@}
293
294 /// @name Deprecated functionality
295 //@{
296
297 /// @brief Add particle by raw pointer
298 /// @deprecated Use GenEvent::add_particle( const GenParticlePtr& ) instead
299 void add_particle( GenParticle *p );
300
301 /// @brief Add vertex by raw pointer
302 /// @deprecated Use GenEvent::add_vertex( const GenVertexPtr& ) instead
303 void add_vertex ( GenVertex *v );
304
305
306 /// @brief Set incoming beam particles
307 /// @deprecated Backward compatibility
308 void set_beam_particles(GenParticlePtr p1, GenParticlePtr p2);
309
310
311 /// @brief Add particle to root vertex
312
313 void add_beam_particle(GenParticlePtr p1);
314
315
316 //@}
317
318#endif // __CINT__
319
320
321 /// @name Methods to fill GenEventData and to read it back
322 //@{
323
324 /// @brief Fill GenEventData object
325 void write_data(GenEventData &data) const;
326
327 /// @brief Fill GenEvent based on GenEventData
328 void read_data(const GenEventData &data);
329
330#ifdef HEPMC3_ROOTIO
331 /// @brief ROOT I/O streamer
332 void Streamer(TBuffer &b);
333 //@}
334#endif
335
336private:
337
338 /// @name Fields
339 //@{
340
341#if !defined(__CINT__)
342
343 /// List of particles
344 std::vector<GenParticlePtr> m_particles;
345 /// List of vertices
346 std::vector<GenVertexPtr> m_vertices;
347
348 /// Event number
349 int m_event_number;
350
351 /// Event weights
352 std::vector<double> m_weights;
353
354 /// Momentum unit
355 Units::MomentumUnit m_momentum_unit;
356 /// Length unit
357 Units::LengthUnit m_length_unit;
358
359 /// The root vertex is stored outside the normal vertices list to block user access to it
360 GenVertexPtr m_rootvertex;
361
362 /// Global run information.
363 std::shared_ptr<GenRunInfo> m_run_info;
364
365 /// @brief Map of event, particle and vertex attributes
366 ///
367 /// Keys are name and ID (0 = event, <0 = vertex, >0 = particle)
368 mutable std::map< std::string, std::map<int, std::shared_ptr<Attribute> > > m_attributes;
369
370 /// @brief Attribute map key type
371 typedef std::map< std::string, std::map<int, std::shared_ptr<Attribute> > >::value_type att_key_t;
372
373 /// @brief Attribute map value type
374 typedef std::map<int, std::shared_ptr<Attribute> >::value_type att_val_t;
375
376 /// @brief Mutex lock for the m_attibutes map.
377 mutable std::recursive_mutex m_lock_attributes;
378#endif // __CINT__
379
380 //@}
381
382};
383
384#if !defined(__CINT__)
385//
386// Template methods
387//
388template<class T>
389std::shared_ptr<T> GenEvent::attribute(const std::string &name, const int& id) const {
390 std::lock_guard<std::recursive_mutex> lock(m_lock_attributes);
391 std::map< std::string, std::map<int, std::shared_ptr<Attribute> > >::iterator i1 =
392 m_attributes.find(name);
393 if( i1 == m_attributes.end() ) {
394 if ( id == 0 && run_info() ) {
395 return run_info()->attribute<T>(name);
396 }
397 return std::shared_ptr<T>();
398 }
399
400 std::map<int, std::shared_ptr<Attribute> >::iterator i2 = i1->second.find(id);
401 if (i2 == i1->second.end() ) return std::shared_ptr<T>();
402
403 if (!i2->second->is_parsed() ) {
404
405 std::shared_ptr<T> att = std::make_shared<T>();
406 att->m_event = this;
407
408 if ( id > 0 && id <= int(particles().size()) )
409 att->m_particle = m_particles[id - 1];
410 if ( id < 0 && -id <= int(vertices().size()) )
411 att->m_vertex = m_vertices[-id - 1];
412 if ( att->from_string(i2->second->unparsed_string()) &&
413 att->init() ) {
414 // update map with new pointer
415 i2->second = att;
416 return att;
417 } else {
418 return std::shared_ptr<T>();
419 }
420 }
421 else return std::dynamic_pointer_cast<T>(i2->second);
422}
423#endif // __CINT__
424
425} // namespace HepMC3
426#endif
Note: See TracBrowser for help on using the repository browser.