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 
14 #include "HepMC3/HEPEVT_Helpers.h"
15 #include "HepMC3/HEPEVT_Wrapper.h"
18 namespace HepMC3
19 {
20 
21 
22 /** @brief comparison of two particles */
23 bool 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*/
34 bool 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 */
76 void 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 
90 HEPMC3_EXPORT_API struct HEPEVT* hepevtptr = nullptr;
91 HEPMC3_EXPORT_API std::shared_ptr<struct HEPEVT_Pointers<double> > HEPEVT_Wrapper_Runtime_Static::m_hepevtptr = nullptr;
92 HEPMC3_EXPORT_API int HEPEVT_Wrapper_Runtime_Static::m_max_particles = 0;
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 
116 void 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 
127 void 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 
155 int 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 
162 int 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 
176 void 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 
182 void 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 
189 void 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 
198 void HEPEVT_Wrapper_Runtime::set_mass( const int index, double mass )
199 {
200  m_hepevtptr->phep[5*(index-1)+4] = mass;
201 }
202 
203 
204 void 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 }
double x(const int index) const
Get X Production vertex.
HepMC3 main namespace.
int first_parent(const int index) const
Get index of 1st mother.
double py(const int index) const
Get Y momentum.
void allocate_internal_storage()
Allocates m_internal_storage storage in smart pointer to hold HEPEVT of fixed size.
int number_entries() const
Get number of entries.
double pz(const int index) const
Get Z momentum.
int number_parents(const int index) const
Get number of parents.
std::shared_ptr< struct HEPEVT_Pointers< double > > m_hepevtptr
Fortran common block HEPEVT.
double z(const int index) const
Get Z Production vertex.
bool operator()(ConstGenParticlePtr lx, ConstGenParticlePtr rx) const
comparison of two particles
void set_children(const int index, const int firstchild, const int lastchild)
Set children.
Helper functions used to manipulate with HEPEVT block.
void set_momentum(const int index, const double px, const double py, const double pz, const double e)
Set 4-momentum.
void set_mass(const int index, double mass)
Set mass.
void set_parents(const int index, const int firstparent, const int lastparent)
Set parents.
void set_position(const int index, const double x, const double y, const double z, const double t)
Set position in time-space.
bool fix_daughters()
Tries to fix list of daughters.
void print_hepevt(std::ostream &ostr=std::cout) const
Print information from HEPEVT common block.
void zero_everything()
Set all entries in HEPEVT to zero.
int last_child(const int index) const
Get index of last daughter.
static HEPMC3_EXPORT_API int m_max_particles
Block size.
double t(const int index) const
Get production time.
std::vector< char > m_internal_storage
Internalstorage storage. Optional.
int last_parent(const int index) const
Get index of last mother.
void calculate_longest_path_to_top(ConstGenVertexPtr v, std::map< ConstGenVertexPtr, int > &pathl)
Calculates the path to the top (beam) particles.
double e(const int index) const
Get Energy.
Fortran common block HEPEVT.
static HEPMC3_EXPORT_API std::shared_ptr< struct HEPEVT_Pointers< double > > m_hepevtptr
Fortran common block HEPEVT.
int number_children_exact(const int index) const
Get number of children by counting.
double y(const int index) const
Get Y Production vertex.
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. We cannot use id, as it can be assigned in different way.
void print_hepevt_particle(int index, std::ostream &ostr=std::cout) const
Print particle information.
double px(const int index) const
Get X momentum.
int first_child(const int index) const
Get index of 1st daughter.
HEPMC3_EXPORT_API struct HEPEVT * hepevtptr
Pointer to HEPEVT common block.
void set_hepevt_address(char *c)
Set Fortran block address.
Definition of class HEPEVT_Wrapper_Runtime_Static.
Definition of class HEPEVT_Wrapper_Runtime.
void copy_to_internal_storage(char *c, int N)
Copies the content of foreight common block into the internal storage.
Definition of class HEPEVT_Wrapper.
int number_children(const int index) const
Get number of children from the range of daughters.