HepMC3 event record library
HEPEVT_Wrapper.cc
Go to the documentation of this file.
1// -*- C++ -*-
2//
3// This file is part of HepMC
4// Copyright (C) 2014-2021 The HepMC collaboration (see AUTHORS for details)
5//
6/**
7 * @file HEPEVT_Wrapper.cc
8 * @brief Implementation of helper functions used to manipulate with HEPEVT block
9 */
10#include <algorithm>
11#include <set>
12#include <vector>
13
18namespace HepMC3
19{
20
21
22/** @brief comparison of two particles */
23bool GenParticlePtr_greater::operator()(ConstGenParticlePtr lx, ConstGenParticlePtr rx) const
24{
25 /* Cannot use id as it could be different.*/
26 if (lx->pid() != rx->pid()) return (lx->pid() < rx->pid());
27 if (lx->status() != rx->status()) return (lx->status() < rx->status());
28 /*Hopefully it will reach this point not too often.*/
29 return (lx->momentum().e() < rx->momentum().e());
30}
31
32/** @brief Order vertices with equal paths. If the paths are equal, order in other quantities.
33 * We cannot use id, as it can be assigned in different way*/
34bool pair_GenVertexPtr_int_greater::operator()(const std::pair<ConstGenVertexPtr, int>& lx, const std::pair<ConstGenVertexPtr, int>& rx) const
35{
36 if (lx.second != rx.second) return (lx.second < rx.second);
37 if (lx.first->particles_in().size() != rx.first->particles_in().size()) return (lx.first->particles_in().size() < rx.first->particles_in().size());
38 if (lx.first->particles_out().size() != rx.first->particles_out().size()) return (lx.first->particles_out().size() < rx.first->particles_out().size());
39 /* The code below is usefull mainly for debug. Assures strong ordering.*/
40 std::vector<int> lx_id_in;
41 std::vector<int> rx_id_in;
42 for (ConstGenParticlePtr pp: lx.first->particles_in()) lx_id_in.push_back(pp->pid());
43 for (ConstGenParticlePtr pp: rx.first->particles_in()) rx_id_in.push_back(pp->pid());
44 std::sort(lx_id_in.begin(), lx_id_in.end());
45 std::sort(rx_id_in.begin(), rx_id_in.end());
46 for (unsigned int i = 0; i < lx_id_in.size(); i++) if (lx_id_in[i] != rx_id_in[i]) return (lx_id_in[i] < rx_id_in[i]);
47
48 std::vector<int> lx_id_out;
49 std::vector<int> rx_id_out;
50 for (ConstGenParticlePtr pp: lx.first->particles_in()) lx_id_out.push_back(pp->pid());
51 for (ConstGenParticlePtr pp: rx.first->particles_in()) rx_id_out.push_back(pp->pid());
52 std::sort(lx_id_out.begin(), lx_id_out.end());
53 std::sort(rx_id_out.begin(), rx_id_out.end());
54 for (unsigned int i = 0; i < lx_id_out.size(); i++) if (lx_id_out[i] != rx_id_out[i]) return (lx_id_out[i] < rx_id_out[i]);
55
56 std::vector<double> lx_mom_in;
57 std::vector<double> rx_mom_in;
58 for (ConstGenParticlePtr pp: lx.first->particles_in()) lx_mom_in.push_back(pp->momentum().e());
59 for (ConstGenParticlePtr pp: rx.first->particles_in()) rx_mom_in.push_back(pp->momentum().e());
60 std::sort(lx_mom_in.begin(), lx_mom_in.end());
61 std::sort(rx_mom_in.begin(), rx_mom_in.end());
62 for (unsigned int i = 0; i < lx_mom_in.size(); i++) if (lx_mom_in[i] != rx_mom_in[i]) return (lx_mom_in[i] < rx_mom_in[i]);
63
64 std::vector<double> lx_mom_out;
65 std::vector<double> rx_mom_out;
66 for (ConstGenParticlePtr pp: lx.first->particles_in()) lx_mom_out.push_back(pp->momentum().e());
67 for (ConstGenParticlePtr pp: rx.first->particles_in()) rx_mom_out.push_back(pp->momentum().e());
68 std::sort(lx_mom_out.begin(), lx_mom_out.end());
69 std::sort(rx_mom_out.begin(), rx_mom_out.end());
70 for (unsigned int i = 0; i < lx_mom_out.size(); i++) if (lx_mom_out[i] != rx_mom_out[i]) return (lx_mom_out[i] < rx_mom_out[i]);
71 /* The code above is usefull mainly for debug*/
72
73 return (lx.first < lx.first); /*This is random. This should never happen*/
74}
75/** @brief Calculates the path to the top (beam) particles */
76void calculate_longest_path_to_top(ConstGenVertexPtr v, std::map<ConstGenVertexPtr, int>& pathl)
77{
78 int p = 0;
79 for (ConstGenParticlePtr pp: v->particles_in()) {
80 ConstGenVertexPtr v2 = pp->production_vertex();
81 if (v2 == v) continue; //LOOP! THIS SHOULD NEVER HAPPEN FOR A PROPER EVENT!
82 if (!v2) p = std::max(p, 1);
83 else
84 {if (pathl.find(v2) == pathl.end()) calculate_longest_path_to_top(v2, pathl); p = std::max(p, pathl[v2]+1);}
85 }
86 pathl[v] = p;
87 return;
88}
89
90HEPMC3_EXPORT_API struct HEPEVT* hepevtptr = nullptr;
91HEPMC3_EXPORT_API std::shared_ptr<struct HEPEVT_Pointers<double> > HEPEVT_Wrapper_Runtime_Static::m_hepevtptr = nullptr;
93
94
96 m_hepevtptr = std::make_shared<struct HEPEVT_Pointers<double> >();
97 char* x = c;
98 m_hepevtptr->nevhep = (int*)x;
99 x += sizeof(int);
100 m_hepevtptr->nhep = (int*)(x);
101 x += sizeof(int);
102 m_hepevtptr->isthep = (int*)(x);
103 x += sizeof(int)*m_max_particles;
104 m_hepevtptr->idhep = (int*)(x);
105 x += sizeof(int)*m_max_particles;
106 m_hepevtptr->jmohep = (int*)(x);
107 x += sizeof(int)*m_max_particles*2;
108 m_hepevtptr->jdahep = (int*)(x);
109 x += sizeof(int)*m_max_particles*2;
110 m_hepevtptr->phep = (double*)(x);
111 x += sizeof(double)*m_max_particles*5;
112 m_hepevtptr->vhep = (double*)(x);
113}
114
115
116void HEPEVT_Wrapper_Runtime::print_hepevt( std::ostream& ostr ) const
117{
118 ostr << " Event No.: " << *(m_hepevtptr->nevhep) << std::endl;
119 ostr << " Nr Type Parent(s) Daughter(s) Px Py Pz E Inv. M." << std::endl;
120 for ( int i = 1; i <= *(m_hepevtptr->nhep); ++i )
121 {
122 print_hepevt_particle( i, ostr );
123 }
124}
125
126
127void HEPEVT_Wrapper_Runtime::print_hepevt_particle( int index, std::ostream& ostr ) const
128{
129 char buf[255];//Note: the format is fixed, so no reason for complicated treatment
130
131 sprintf(buf, "%5i %6i", index, m_hepevtptr->idhep[index-1]);
132 ostr << buf;
133 sprintf(buf, "%4i - %4i ", m_hepevtptr->jmohep[2*(index-1)], m_hepevtptr->jmohep[2*(index-1)+1]);
134 ostr << buf;
135 sprintf(buf, "%4i - %4i ", m_hepevtptr->jdahep[2*(index-1)], m_hepevtptr->jdahep[2*(index-1)+1]);
136 ostr << buf;
137 sprintf(buf, "%8.2f %8.2f %8.2f %8.2f %8.2f", m_hepevtptr->phep[5*(index-1)], m_hepevtptr->phep[5*(index-1)+1], m_hepevtptr->phep[5*(index-1)+2], m_hepevtptr->phep[5*(index-1)+3], m_hepevtptr->phep[5*(index-1)+4]);
138 ostr << buf << std::endl;
139}
140
141
143{
144 *(m_hepevtptr->nevhep) = 0;
145 *(m_hepevtptr->nhep) = 0;
146 memset(m_hepevtptr->isthep, 0, sizeof(int)*m_max_particles);
147 memset(m_hepevtptr->idhep, 0, sizeof(int)*m_max_particles);
148 memset(m_hepevtptr->jmohep, 0, sizeof(int)*m_max_particles*2);
149 memset(m_hepevtptr->jdahep, 0, sizeof(int)*m_max_particles*2);
150 memset(m_hepevtptr->phep, 0, sizeof(double)*m_max_particles*5);
151 memset(m_hepevtptr->vhep, 0, sizeof(double)*m_max_particles*4);
152}
153
154
155int HEPEVT_Wrapper_Runtime::number_parents( const int index ) const
156{
157 return (m_hepevtptr->jmohep[2*(index-1)]) ? (m_hepevtptr->jmohep[2*(index-1)+1]) ? m_hepevtptr->jmohep[2*(index-1)+1]
158 -m_hepevtptr->jmohep[2*(index-1)] : 1 : 0;
159}
160
161
162int HEPEVT_Wrapper_Runtime::number_children( const int index ) const
163{
164 return (m_hepevtptr->jdahep[2*(index-1)]) ? (m_hepevtptr->jdahep[2*(index-1)+1]) ? m_hepevtptr->jdahep[2*(index-1)+1]-m_hepevtptr->jdahep[2*(index-1)] : 1 : 0;
165}
166
168{
169 int nc = 0;
170 for ( int i = 1; i <= *(m_hepevtptr->nhep); ++i )
171 if (((m_hepevtptr->jmohep[2*(i-1)] <= index && m_hepevtptr->jmohep[2*(i-1)+1] >= index)) || (m_hepevtptr->jmohep[2*(i-1)] == index) ||
172 (m_hepevtptr->jmohep[2*(index-1)+1]==index)) nc++;
173 return nc;
174}
175
176void HEPEVT_Wrapper_Runtime::set_parents( const int index, const int firstparent, const int lastparent )
177{
178 m_hepevtptr->jmohep[2*(index-1)] = firstparent;
179 m_hepevtptr->jmohep[2*(index-1)+1] = lastparent;
180}
181
182void HEPEVT_Wrapper_Runtime::set_children( const int index, const int firstchild, const int lastchild )
183{
184 m_hepevtptr->jdahep[2*(index-1)] = firstchild;
185 m_hepevtptr->jdahep[2*(index-1)+1] = lastchild;
186}
187
188
189void HEPEVT_Wrapper_Runtime::set_momentum( const int index, const double px, const double py, const double pz, const double e )
190{
191 m_hepevtptr->phep[5*(index-1)] = px;
192 m_hepevtptr->phep[5*(index-1)+1] = py;
193 m_hepevtptr->phep[5*(index-1)+2] = pz;
194 m_hepevtptr->phep[5*(index-1)+3] = e;
195}
196
197
198void HEPEVT_Wrapper_Runtime::set_mass( const int index, double mass )
199{
200 m_hepevtptr->phep[5*(index-1)+4] = mass;
201}
202
203
204void HEPEVT_Wrapper_Runtime::set_position( const int index, const double x, const double y, const double z, const double t )
205{
206 m_hepevtptr->vhep[4*(index-1)] = x;
207 m_hepevtptr->vhep[4*(index-1)+1] = y;
208 m_hepevtptr->vhep[4*(index-1)+2] = z;
209 m_hepevtptr->vhep[4*(index-1)+3] = t;
210}
211
212
214{
215 /*AV The function should be called for a record that has correct particle ordering and mother ids.
216 As a result it produces a record with ranges where the daughters can be found.
217 Not every particle in the range will be a daughter. It is true only for proper events.
218 The return tells if the record was fixed succesfully.
219 */
220 for ( int i = 1; i <= number_entries(); i++ )
221 for ( int k=1; k <= number_entries(); k++ ) if (i != k)
222 if ((first_parent(k) <= i) && (i <= last_parent(k)))
223 set_children(i, (first_child(i) == 0 ? k : std::min(first_child(i), k)), (last_child(i) == 0 ? k : std::max(last_child(i), k)));
224 bool is_fixed = true;
225 for ( int i = 1; i <= number_entries(); i++ )
226 is_fixed = (is_fixed && (number_children_exact(i) == number_children(i)));
227 return is_fixed;
228}
229
230
232{
233 m_internal_storage.reserve(2*sizeof(int)+m_max_particles*(6*sizeof(int)+9*sizeof(double)));
235}
236
238{
239 if ( N < 1 || N > m_max_particles) return;
240 char* dest = m_internal_storage.data();
241 char* src = c;
242 memcpy(dest, src, 2*sizeof(int));
243 src += 2*sizeof(int);
244 dest += 2*sizeof(int);
245 memcpy(dest, src, N*sizeof(int));
246 src += N*sizeof(int);
247 dest += m_max_particles*sizeof(int);
248 memcpy(dest, src, N*sizeof(int));
249 src += N*sizeof(int);
250 dest += m_max_particles*sizeof(int);
251 memcpy(dest, src, 2*N*sizeof(int));
252 src += 2*N*sizeof(int);
253 dest += 2*m_max_particles*sizeof(int);
254 memcpy(dest, src, 2*N*sizeof(int));
255 src += 2*N*sizeof(int);
256 dest += 2*m_max_particles*sizeof(int);
257 memcpy(dest, src, 5*N*sizeof(double));
258 src += 5*N*sizeof(double);
259 dest += 5*m_max_particles*sizeof(double);
260 memcpy(dest, src, 4*N*sizeof(double));
261}
262
263}
Helper functions used to manipulate with HEPEVT block.
Definition of class HEPEVT_Wrapper.
HEPMC3_EXPORT_API struct HEPEVT * hepevtptr
Pointer to HEPEVT common block.
Definition of class HEPEVT_Wrapper_Runtime.
Definition of class HEPEVT_Wrapper_Runtime_Static.
static HEPMC3_EXPORT_API int m_max_particles
Block size.
static HEPMC3_EXPORT_API std::shared_ptr< struct HEPEVT_Pointers< double > > m_hepevtptr
Fortran common block HEPEVT.
int first_child(const int index) const
Get index of 1st daughter.
int last_parent(const int index) const
Get index of last mother.
int last_child(const int index) const
Get index of last daughter.
void allocate_internal_storage()
Allocates m_internal_storage storage in smart pointer to hold HEPEVT of fixed size.
int first_parent(const int index) const
Get index of 1st mother.
double pz(const int index) const
Get Z momentum.
void zero_everything()
Set all entries in HEPEVT to zero.
double py(const int index) const
Get Y momentum.
void set_parents(const int index, const int firstparent, const int lastparent)
Set parents.
int number_children_exact(const int index) const
Get number of children by counting.
void print_hepevt_particle(int index, std::ostream &ostr=std::cout) const
Print particle information.
double t(const int index) const
Get production time.
void set_momentum(const int index, const double px, const double py, const double pz, const double e)
Set 4-momentum.
void set_children(const int index, const int firstchild, const int lastchild)
Set children.
int number_parents(const int index) const
Get number of parents.
void set_position(const int index, const double x, const double y, const double z, const double t)
Set position in time-space.
double y(const int index) const
Get Y Production vertex.
bool fix_daughters()
Tries to fix list of daughters.
int number_entries() const
Get number of entries.
void set_hepevt_address(char *c)
Set Fortran block address.
double px(const int index) const
Get X momentum.
std::vector< char > m_internal_storage
Internalstorage storage. Optional.
int number_children(const int index) const
Get number of children from the range of daughters.
double z(const int index) const
Get Z Production vertex.
void print_hepevt(std::ostream &ostr=std::cout) const
Print information from HEPEVT common block.
void set_mass(const int index, double mass)
Set mass.
double x(const int index) const
Get X Production vertex.
void copy_to_internal_storage(char *c, int N)
Copies the content of foreight common block into the internal storage.
std::shared_ptr< struct HEPEVT_Pointers< double > > m_hepevtptr
Fortran common block HEPEVT.
double e(const int index) const
Get Energy.
HepMC3 main namespace.
void calculate_longest_path_to_top(ConstGenVertexPtr v, std::map< ConstGenVertexPtr, int > &pathl)
Calculates the path to the top (beam) particles.
Fortran common block HEPEVT.
bool operator()(ConstGenParticlePtr lx, ConstGenParticlePtr rx) const
comparison of two particles
bool operator()(const std::pair< ConstGenVertexPtr, int > &lx, const std::pair< ConstGenVertexPtr, int > &rx) const
Order vertices with equal paths. If the paths are equal, order in other quantities....