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 }
bool operator()(ConstGenParticlePtr lx, ConstGenParticlePtr rx) const
comparison of two particles
int number_entries() const
Get number of entries.
void allocate_internal_storage()
Allocates m_internal_storage storage in smart pointer to hold HEPEVT of fixed size.
int number_children(const int index) const
Get number of children from the range of daughters.
std::shared_ptr< struct HEPEVT_Pointers< double > > m_hepevtptr
Fortran common block HEPEVT.
int number_parents(const int index) const
Get number of parents.
void set_children(const int index, const int firstchild, const int lastchild)
Set children.
Helper functions used to manipulate with HEPEVT block.
void print_hepevt(std::ostream &ostr=std::cout) const
Print information from HEPEVT common block.
void print_hepevt_particle(int index, std::ostream &ostr=std::cout) const
Print particle information.
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.
int last_child(const int index) const
Get index of last daughter.
void zero_everything()
Set all entries in HEPEVT to zero.
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.
double x(const int index) const
Get X Production vertex.
int number_children_exact(const int index) const
Get number of children by counting.
int first_parent(const int index) const
Get index of 1st mother.
int first_child(const int index) const
Get index of 1st daughter.
static HEPMC3_EXPORT_API int m_max_particles
Block size.
std::vector< char > m_internal_storage
Internalstorage storage. Optional.
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.
static HEPMC3_EXPORT_API std::shared_ptr< struct HEPEVT_Pointers< double > > m_hepevtptr
Fortran common block HEPEVT.
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.
int last_parent(const int index) const
Get index of last mother.
Definition of class HEPEVT_Wrapper.