HepMC event record
include/HepMC/HEPEVT_Wrapper.h
1 // -*- C++ -*-
2 //
3 // This file is part of HepMC
4 // Copyright (C) 2014-2015 The HepMC collaboration (see AUTHORS for details)
5 //
6 #ifndef HEPMC_HEPEVT_WRAPPER_H
7 #define HEPMC_HEPEVT_WRAPPER_H
8 /**
9  * @file HEPEVT_Wrapper.h
10  * @brief Definition of \b class HEPEVT_Wrapper
11  *
12  * @class HepMC::HEPEVT_Wrapper
13  * @brief An interface to HEPEVT common block
14  *
15  * @note This header file does not include HEPEVT definition, only declaration.
16  * Including this wrapper requires that HEPEVT is defined somewhere
17  * in the project (most likely as FORTRAN common block).
18  *
19  * @note Make sure that HEPEVT definition in project matches this definition
20  * (NMXHEP, double precision, etc.) Change this definition if necessary.
21  *
22  * @todo Do we just make write_event and fill_next_event instead?
23  */
24 
25 #ifndef HEPMC_HEPEVT_NMXHEP
26 /** Default number of particles in the HEPEVT structure */
27 #define HEPMC_HEPEVT_NMXHEP 10000
28 #endif
29 
30 #ifndef HEPMC_HEPEVT_PRECISION
31 /** Default precision of the 4-momentum, time-space position and mass */
32 #define HEPMC_HEPEVT_PRECISION double
33 #endif
34 
35 /* This definition of HEPEVT corresponds to FORTRAN definition:
36 
37  PARAMETER (NMXHEP=10000)
38  COMMON /HEPEVT/ NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
39  & JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
40  INTEGER NEVHEP,NHEP,ISTHEP,IDHEP,JMOHEP,JDAHEP
41  DOUBLE PRECISION PHEP,VHEP
42 */
43 
44 static const int NMXHEP = HEPMC_HEPEVT_NMXHEP; //!< Number of particles in the HEPEVT structure
45 typedef HEPMC_HEPEVT_PRECISION momentum_t; //!< Precision of the 4-momentum, time-space position and mass
46 
47 /** @struct HEPEVT
48  * @brief C structure representing Fortran common block HEPEVT
49  * T. Sjöstrand et al., "A proposed standard event record",
50  * in `Z physics at LEP 1', eds. G. Altarelli, R. Kleiss and C. Verzegnassi,
51  * Geneva, Switzerland, September 4-5, 1989, CERN 89-08 (Geneva, 1989), Vol. 3, p. 327
52  * Disk representation is given by Fortran WRITE/READ format.
53  */
54 struct HEPEVT
55 {
56  int nevhep; //!< Event number
57  int nhep; //!< Number of entries in the event
58  int isthep[NMXHEP]; //!< Status code
59  int idhep [NMXHEP]; //!< PDG ID
60  int jmohep[NMXHEP][2]; //!< Pointer to position of 1st and 2nd (or last!) mother
61  int jdahep[NMXHEP][2]; //!< Pointer to position of 1nd and 2nd (or last!) daughter
62  momentum_t phep [NMXHEP][5]; //!< Momentum: px, py, pz, e, m
63  momentum_t vhep [NMXHEP][4]; //!< Time-space position: x, y, z, t
64 }; //!< Fortran common block HEPEVT
65 
66 
67 
68 #include <iostream>
69 #include <cstdio>
70 #include <cstring> // memset
71 #ifndef HEPEVT_WRAPPER_HEADER_ONLY
72 #include "HepMC/GenEvent.h"
73 #endif
74 using std::endl;
75 
76 namespace HepMC
77 {
78 extern struct HEPEVT* hepevtptr;
79 
80 class HEPEVT_Wrapper
81 {
82 
83 //
84 // Functions
85 //
86 public:
87  /** @brief Print information from HEPEVT common block */
88  static void print_hepevt( std::ostream& ostr = std::cout );
89 
90  /** @brief Print particle information */
91  static void print_hepevt_particle( int index, std::ostream& ostr = std::cout );
92 
93  /** @brief Check for problems with HEPEVT common block */
94  static bool check_hepevt_consistency( std::ostream& ostr = std::cout );
95 
96  /** @brief Set all entries in HEPEVT to zero */
97  static void zero_everything();
98 #ifndef HEPEVT_WRAPPER_HEADER_ONLY
99  /** @brief Convert GenEvent to HEPEVT*/
100  static bool GenEvent_to_HEPEVT( const GenEvent* evt );
101  /** @brief Convert HEPEVT to GenEvent*/
102  static bool HEPEVT_to_GenEvent( GenEvent* evt );
103 #endif
104  /** @brief Tries to fix list of daughters */
105  static bool fix_daughters();
106 //
107 // Accessors
108 //
109 public:
110  static void set_hepevt_address(char *c) { hepevtptr=(struct HEPEVT*)c; } //!< Set Fortran block address
111  static int max_number_entries() { return NMXHEP; } //!< Block size
112  static int event_number() { return hepevtptr->nevhep; } //!< Get event number
113  static int number_entries() { return hepevtptr->nhep; } //!< Get number of entries
114  static int status( int index ) { return hepevtptr->isthep[index-1]; } //!< Get status code
115  static int id( int index ) { return hepevtptr->idhep[index-1]; } //!< Get PDG particle id
116  static int first_parent( int index ) { return hepevtptr->jmohep[index-1][0]; } //!< Get index of 1st mother
117  static int last_parent( int index ) { return hepevtptr->jmohep[index-1][1]; } //!< Get index of last mother
118  static int first_child( int index ) { return hepevtptr->jdahep[index-1][0]; } //!< Get index of 1st daughter
119  static int last_child( int index ) { return hepevtptr->jdahep[index-1][1]; } //!< Get index of last daughter
120  static double px( int index ) { return hepevtptr->phep[index-1][0]; } //!< Get X momentum
121  static double py( int index ) { return hepevtptr->phep[index-1][1]; } //!< Get Y momentum
122  static double pz( int index ) { return hepevtptr->phep[index-1][2]; } //!< Get Z momentum
123  static double e( int index ) { return hepevtptr->phep[index-1][3]; } //!< Get Energy
124  static double m( int index ) { return hepevtptr->phep[index-1][4]; } //!< Get generated mass
125  static double x( int index ) { return hepevtptr->vhep[index-1][0]; } //!< Get X Production vertex
126  static double y( int index ) { return hepevtptr->vhep[index-1][1]; } //!< Get Y Production vertex
127  static double z( int index ) { return hepevtptr->vhep[index-1][2]; } //!< Get Z Production vertex
128  static double t( int index ) { return hepevtptr->vhep[index-1][3]; } //!< Get production time
129  static int number_parents( int index ); //!< Get number of parents
130  static int number_children( int index ); //!< Get number of children from the range of daughters
131  static int number_children_exact( int index ); //!< Get number of children by counting
132 
133  static void set_event_number( int evtno ) { hepevtptr->nevhep = evtno; } //!< Set event number
134  static void set_number_entries( int noentries ) { hepevtptr->nhep = noentries; } //!< Set number of entries
135  static void set_status( int index, int status ) { hepevtptr->isthep[index-1] = status; } //!< Set status code
136  static void set_id( int index, int id ) { hepevtptr->idhep[index-1] = id; } //!< Set PDG particle id
137  static void set_parents( int index, int firstparent, int lastparent ); //!< Set parents
138  static void set_children( int index, int firstchild, int lastchild ); //!< Set children
139  static void set_momentum( int index, double px, double py, double pz, double e ); //!< Set 4-momentum
140  static void set_mass( int index, double mass ); //!< Set mass
141  static void set_position( int index, double x, double y, double z, double t ); //!< Set position in time-space
142 };
143 
144 //
145 // inline definitions
146 //
147 inline void HEPEVT_Wrapper::print_hepevt( std::ostream& ostr )
148 {
149  ostr << " Event No.: " << hepevtptr->nevhep << endl;
150  ostr<< " Nr Type Parent(s) Daughter(s) Px Py Pz E Inv. M." << endl;
151  for( int i=1; i<=hepevtptr->nhep; ++i )
152  {
154  }
155 }
156 
157 inline void HEPEVT_Wrapper::print_hepevt_particle( int index, std::ostream& ostr )
158 {
159  char buf[255];
160 
161  sprintf(buf,"%5i %6i",index,hepevtptr->idhep[index-1]);
162  ostr << buf;
163  sprintf(buf,"%4i - %4i ",hepevtptr->jmohep[index-1][0],hepevtptr->jmohep[index-1][1]);
164  ostr << buf;
165  sprintf(buf,"%4i - %4i ",hepevtptr->jdahep[index-1][0],hepevtptr->jdahep[index-1][1]);
166  ostr << buf;
167  // print the rest of particle info
168  sprintf(buf,"%8.2f %8.2f %8.2f %8.2f %8.2f",hepevtptr->phep[index-1][0],hepevtptr->phep[index-1][1],hepevtptr->phep[index-1][2],hepevtptr->phep[index-1][3],hepevtptr->phep[index-1][4]);
169  ostr << buf << endl;
170 }
171 
172 
173 inline bool HEPEVT_Wrapper::check_hepevt_consistency( std::ostream& /*ostr*/ )
174 {
175  //!< @todo HEPEVT_Wrapper::check_hepevt_consistency unimplemented!
176  printf("HEPEVT_Wrapper::check_hepevt_consistency unimplemented!\n");
177  return true;
178 }
179 
181 {
182  memset(hepevtptr,0,sizeof(struct HEPEVT));
183 }
184 
185 inline int HEPEVT_Wrapper::number_parents( int index )
186 {
187  return (hepevtptr->jmohep[index-1][0]) ? (hepevtptr->jmohep[index-1][1]) ? hepevtptr->jmohep[index-1][1]-hepevtptr->jmohep[index-1][0] : 1 : 0;
188 }
189 
190 inline int HEPEVT_Wrapper::number_children( int index )
191 {
192  return (hepevtptr->jdahep[index-1][0]) ? (hepevtptr->jdahep[index-1][1]) ? hepevtptr->jdahep[index-1][1]-hepevtptr->jdahep[index-1][0] : 1 : 0;
193 }
194 
195 inline int HEPEVT_Wrapper::number_children_exact( int index )
196 {
197  int nc=0;
198  for( int i=1; i<=hepevtptr->nhep; ++i )
199  if (((hepevtptr->jmohep[i-1][0]<=index&&hepevtptr->jmohep[i-1][1]>=index))||(hepevtptr->jmohep[i-1][0]==index)||(hepevtptr->jmohep[i-1][1]==index)) nc++;
200  return nc;
201 }
202 
203 
204 
205 inline void HEPEVT_Wrapper::set_parents( int index, int firstparent, int lastparent )
206 {
207  hepevtptr->jmohep[index-1][0] = firstparent;
208  hepevtptr->jmohep[index-1][1] = lastparent;
209 }
210 
211 inline void HEPEVT_Wrapper::set_children( int index, int firstchild, int lastchild )
212 {
213  hepevtptr->jdahep[index-1][0] = firstchild;
214  hepevtptr->jdahep[index-1][1] = lastchild;
215 }
216 
217 inline void HEPEVT_Wrapper::set_momentum( int index, double px, double py, double pz, double e )
218 {
219  hepevtptr->phep[index-1][0] = px;
220  hepevtptr->phep[index-1][1] = py;
221  hepevtptr->phep[index-1][2] = pz;
222  hepevtptr->phep[index-1][3] = e;
223 }
224 
225 inline void HEPEVT_Wrapper::set_mass( int index, double mass )
226 {
227  hepevtptr->phep[index-1][4] = mass;
228 }
229 
230 inline void HEPEVT_Wrapper::set_position( int index, double x, double y, double z, double t )
231 {
232  hepevtptr->vhep[index-1][0] = x;
233  hepevtptr->vhep[index-1][1] = y;
234  hepevtptr->vhep[index-1][2] = z;
235  hepevtptr->vhep[index-1][3] = t;
236 }
237 
238 
239 
240 inline bool HEPEVT_Wrapper::fix_daughters()
241 {
242  /*AV The function should be called for a record that has correct particle ordering and mother ids.
243  As a result it produces a record with ranges where the daughters can be found.
244  Not every particle in the range will be a daughter. It is true only for proper events.
245  The return tells if the record was fixed succesfully.
246  */
247 
248  for( int i=1; i<=HEPEVT_Wrapper::number_entries(); i++ )
249  for( int k=1; k<=HEPEVT_Wrapper::number_entries(); k++ ) if (i!=k)
252  bool is_fixed=true;
253  for( int i=1; i<=HEPEVT_Wrapper::number_entries(); i++ )
255  return is_fixed;
256 
257  return true;
258 }
259 
260 
261 
262 
263 
264 
265 
266 } // namespace HepMC
267 
268 #endif
static double py(int index)
Get Y momentum.
static void set_number_entries(int noentries)
Set number of entries.
static int status(int index)
Get status code.
static bool GenEvent_to_HEPEVT(const GenEvent *evt)
Convert GenEvent to HEPEVT.
static void set_status(int index, int status)
Set status code.
static double y(int index)
Get Y Production vertex.
static double e(int index)
Get Energy.
static double px(int index)
Get X momentum.
static void set_parents(int index, int firstparent, int lastparent)
Set parents.
static double pz(int index)
Get Z momentum.
static void set_hepevt_address(char *c)
Set Fortran block address.
static int first_child(int index)
Get index of 1st daughter.
momentum_t phep[NMXHEP][5]
Momentum: px, py, pz, e, m.
int nhep
Number of entries in the event.
static double z(int index)
Get Z Production vertex.
static int max_number_entries()
Block size.
momentum_t vhep[NMXHEP][4]
Time-space position: x, y, z, t.
int jmohep[NMXHEP][2]
Pointer to position of 1st and 2nd (or last!) mother.
int isthep[NMXHEP]
Status code.
static void print_hepevt(std::ostream &ostr=std::cout)
Print information from HEPEVT common block.
static int number_parents(int index)
Get number of parents.
static double t(int index)
Get production time.
static void set_id(int index, int id)
Set PDG particle id.
static int last_child(int index)
Get index of last daughter.
static bool check_hepevt_consistency(std::ostream &ostr=std::cout)
Check for problems with HEPEVT common block.
Fortran common block HEPEVT.
static int event_number()
Get event number.
static void set_event_number(int evtno)
Set event number.
static double x(int index)
Get X Production vertex.
static int first_parent(int index)
Get index of 1st mother.
static int id(int index)
Get PDG particle id.
static double m(int index)
Get generated mass.
static void print_hepevt_particle(int index, std::ostream &ostr=std::cout)
Print particle information.
static void set_children(int index, int firstchild, int lastchild)
Set children.
static void set_mass(int index, double mass)
Set mass.
static int number_entries()
Get number of entries.
static int last_parent(int index)
Get index of last mother.
static int number_children(int index)
Get number of children from the range of daughters.
static bool fix_daughters()
Tries to fix list of daughters.
Definition of template class SmartPointer.
static int number_children_exact(int index)
Get number of children by counting.
static void zero_everything()
Set all entries in HEPEVT to zero.
static void set_position(int index, double x, double y, double z, double t)
Set position in time-space.
static void set_momentum(int index, double px, double py, double pz, double e)
Set 4-momentum.
static bool HEPEVT_to_GenEvent(GenEvent *evt)
Convert HEPEVT to GenEvent.
int jdahep[NMXHEP][2]
Pointer to position of 1nd and 2nd (or last!) daughter.