HepMC3 event record library
HepMCCompatibility.h
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // This file is part of HepMC
4 // Copyright (C) 2014-2020 The HepMC collaboration (see AUTHORS for details)
5 ///
6 /// @file HepMCCompatibility.h
7 /// @brief Implementation of compatibility layer (in-memory conversion functions) between HePMC2 and HepMC3
8 ///
9 #ifndef HEPMCCOMPATIBILITY_H
10 #define HEPMCCOMPATIBILITY_H
11 
12 #include "HepMC3/GenVertex.h"
13 #include "HepMC3/GenParticle.h"
14 #include "HepMC3/GenEvent.h"
15 
16 ///Please note the HEPMC_HAS_CENTRALITY should be defined externaly
17 #include "HepMC/GenVertex.h"
18 #include "HepMC/GenParticle.h"
19 #include "HepMC/GenEvent.h"
20 
21 /** Converts HepMC3::Genevent to HepMC::Genevent */
22 HepMC::GenEvent* ConvertHepMCGenEvent_3to2(const HepMC3::GenEvent& evt)
23 {
24  HepMC::GenEvent* n=new HepMC::GenEvent();
25  HepMC::Units::LengthUnit lu;
26  HepMC::Units::MomentumUnit mu;
27  if (evt.length_unit()==HepMC3::Units::CM) lu=HepMC::Units::CM;
28  if (evt.length_unit()==HepMC3::Units::MM) lu=HepMC::Units::MM;
29  if (evt.momentum_unit()==HepMC3::Units::MEV) mu=HepMC::Units::MEV;
30  if (evt.momentum_unit()==HepMC3::Units::GEV) mu=HepMC::Units::GEV;
31  n->use_units(mu,lu);
32 
33 
34  std::shared_ptr<HepMC3::IntAttribute> A_signal_process_vertex=evt.attribute<HepMC3::IntAttribute>("signal_process_vertex");
35  int signal_process_vertex_id=0;
36  if (A_signal_process_vertex) signal_process_vertex_id=A_signal_process_vertex->value();
37 
38  std::map<HepMC3::ConstGenVertexPtr,HepMC::GenVertex*> vertexmap3to2;
39  for (auto v3: evt.vertices())
40  {
41  HepMC3::FourVector pos3=v3->position();
42  HepMC::FourVector pos2=HepMC::FourVector(pos3.x(),pos3.y(),pos3.pz(),pos3.t());
43  HepMC::GenVertex* v2= new HepMC::GenVertex(pos2);
44  std::vector<double> vweights;
45  if(v3->attribute_names().size())
46  {
47  std::shared_ptr<HepMC3::VectorDoubleAttribute> rsvec=v3->attribute<HepMC3::VectorDoubleAttribute>("weights");
48  if (rsvec) { vweights =rsvec->value(); }
49  else {
50  for (size_t ii=0; ii<100; ii++)
51  {
52  std::shared_ptr<HepMC3::DoubleAttribute> rs=v3->attribute<HepMC3::DoubleAttribute>("weight"+std::to_string((long long unsigned int)ii));
53  if (!rs) break;
54  vweights.push_back(rs->value());
55  }
56  }
57  }
58  if (vweights.empty())vweights.push_back(1.0);
59  v2->weights()=vweights;
60  v2->suggest_barcode(v3->id());
61  n->add_vertex(v2);
62  if (signal_process_vertex_id==v3->id()) n->set_signal_process_vertex(v2);
63  vertexmap3to2[v3]=v2;
64  }
65  std::map<HepMC3::ConstGenParticlePtr,HepMC::GenParticle*> particlemap3to2;
66  for (auto p3: evt.particles())
67  {
68  HepMC3::FourVector mom3=p3->momentum();
69  HepMC::FourVector mom2=HepMC::FourVector(mom3.px(),mom3.py(),mom3.pz(),mom3.e());
70  HepMC::GenParticle* p2= new HepMC::GenParticle(mom2,p3->pid(),p3->status());
71  if (p3->is_generated_mass_set())p2->set_generated_mass(p3->generated_mass());
72  /** Converts HepMC3::Genevent to HepMC::Genevent */
73  p2->suggest_barcode(10000+p3->id());
74  particlemap3to2[p3]=p2;
75 
76  auto v3production=p3->production_vertex();
77  HepMC::GenVertex* v2production=nullptr;
78  if (v3production)
79  {
80  auto v2=vertexmap3to2.find(v3production);
81  if (v2!=vertexmap3to2.end()) v2production=vertexmap3to2[v3production];
82  }
83 
84  auto v3end=p3->end_vertex();
85  HepMC::GenVertex* v2end=nullptr;
86  if (v3end)
87  {
88  auto v2=vertexmap3to2.find(v3end);
89  if (v2!=vertexmap3to2.end()) v2end=vertexmap3to2[v3end];
90  }
91 
92  if (v2production) v2production->add_particle_out(p2);
93  if (v2end) v2end->add_particle_in(p2);
94 
95  std::shared_ptr<HepMC3::DoubleAttribute> p3_theta=p3->attribute<HepMC3::DoubleAttribute>("theta");
96  std::shared_ptr<HepMC3::DoubleAttribute> p3_phi=p3->attribute<HepMC3::DoubleAttribute>("phi");
97  if (p3_phi && p3_theta)p2->set_polarization(HepMC::Polarization(p3_theta->value(),p3_phi->value()));
98  std::shared_ptr<HepMC3::VectorIntAttribute> flows = p3->attribute<HepMC3::VectorIntAttribute>("flows");
99  if (flows)
100  {
101  std::vector<int> flowvalues=flows->value();
102  for (size_t i=0; i<flowvalues.size(); i++) p2->set_flow(i+1,flowvalues.at(i));
103  } else {
104  std::shared_ptr<HepMC3::IntAttribute> p3_flow1=p3->attribute<HepMC3::IntAttribute>("flow1");
105  std::shared_ptr<HepMC3::IntAttribute> p3_flow2=p3->attribute<HepMC3::IntAttribute>("flow2");
106  std::shared_ptr<HepMC3::IntAttribute> p3_flow3=p3->attribute<HepMC3::IntAttribute>("flow3");
107  if (p3_flow1) p2->set_flow(1,p3_flow1->value());
108  if (p3_flow2) p2->set_flow(2,p3_flow2->value());
109  if (p3_flow3) p2->set_flow(3,p3_flow2->value());
110  }
111  }
112  std::vector<HepMC3::ConstGenParticlePtr> bms=evt.beams();
113  if (bms.size()>=2)
114  {
115  if (particlemap3to2.find(bms[0])!=particlemap3to2.end() && particlemap3to2.find(bms[1])!=particlemap3to2.end() )
116  {
117  n->set_beam_particles(particlemap3to2[bms[0]],particlemap3to2[bms[1]]);
118  }
119  }
120  std::shared_ptr<HepMC3::GenCrossSection> cs3 = evt.attribute<HepMC3::GenCrossSection>("GenCrossSection");
121  if(cs3) {
122  HepMC::GenCrossSection cs2;
123  cs2.set_cross_section(cs3->xsec(),cs3->xsec_err());
124  n->set_cross_section(cs2);
125  }
126  std::shared_ptr<HepMC3::GenPdfInfo> pdf3 = evt.attribute<HepMC3::GenPdfInfo>("GenPdfInfo");
127  if(pdf3)
128  {
129  HepMC::PdfInfo pdf2(
130  pdf3->parton_id[0],
131  pdf3->parton_id[1],
132  pdf3->x[0],
133  pdf3->x[1],
134  pdf3->scale,
135  pdf3->xf[0],
136  pdf3->xf[1],
137  pdf3->pdf_id[0],
138  pdf3->pdf_id[1]
139  );
140  n->set_pdf_info(pdf2);
141  }
142 
143  std::shared_ptr<HepMC3::GenHeavyIon> hi3 = evt.attribute<HepMC3::GenHeavyIon>("GenHeavyIon");
144  if(hi3)
145  {
146  HepMC::HeavyIon hi2;
147  hi2.set_Ncoll_hard(hi3->Ncoll_hard);
148  hi2.set_Npart_proj(hi3->Npart_proj);
149  hi2.set_Npart_targ(hi3->Npart_targ);
150  hi2.set_Ncoll(hi3->Ncoll);
151  hi2.set_spectator_neutrons(hi3->spectator_neutrons);
152  hi2.set_spectator_protons(hi3->spectator_protons);
153  hi2.set_N_Nwounded_collisions(hi3->N_Nwounded_collisions);
154  hi2.set_Nwounded_N_collisions(hi3->Nwounded_N_collisions);
155  hi2.set_Nwounded_Nwounded_collisions(hi3->Nwounded_Nwounded_collisions);
156  hi2.set_impact_parameter(hi3->impact_parameter);
157  hi2.set_event_plane_angle(hi3->event_plane_angle);
158  hi2.set_eccentricity(hi3->eccentricity);
159  hi2.set_sigma_inel_NN(hi3->sigma_inel_NN);
160 #ifdef HEPMC_HAS_CENTRALITY
161  hi2.set_centrality(hi3->centrality);
162 #endif
163  n->set_heavy_ion(hi2);
164  }
165  std::vector<double> wv=evt.weights();
166  std::vector<std::string> wn=evt.weight_names();
167  for (size_t i=0; i<wv.size(); i++) n->weights()[wn.at(i)]=wv.at(i);
168 
169  std::shared_ptr<HepMC3::DoubleAttribute> A_event_scale=evt.attribute<HepMC3::DoubleAttribute>("event_scale");
170  std::shared_ptr<HepMC3::DoubleAttribute> A_alphaQED=evt.attribute<HepMC3::DoubleAttribute>("alphaQED");
171  std::shared_ptr<HepMC3::DoubleAttribute> A_alphaQCD=evt.attribute<HepMC3::DoubleAttribute>("alphaQCD");
172  std::shared_ptr<HepMC3::IntAttribute> A_signal_process_id=evt.attribute<HepMC3::IntAttribute>("signal_process_id");
173  std::shared_ptr<HepMC3::IntAttribute> A_mpi=evt.attribute<HepMC3::IntAttribute>("mpi");
174 
175 
176  double event_scale=A_event_scale?(A_event_scale->value()):0.0;
177  double alphaQED=A_alphaQED?(A_alphaQED->value()):0.0;
178  double alphaQCD=A_alphaQCD?(A_alphaQCD->value()):0.0;
179  int signal_process_id=A_signal_process_id?(A_signal_process_id->value()):0;
180  int mpi=A_mpi?(A_mpi->value()):0;
181 
182  std::vector<long> random_states;
183  for (int i=0; i<100; i++)
184  {
185  std::shared_ptr<HepMC3::IntAttribute> rs=evt.attribute<HepMC3::IntAttribute>("random_states"+std::to_string((long long unsigned int)i));
186  if (!rs) break;
187  random_states.push_back(rs->value());
188  }
189  n->set_mpi( mpi );
190  n->set_signal_process_id( signal_process_id );
191  n->set_event_number( evt.event_number() );
192  n->set_random_states( random_states );
193  n->set_event_scale( event_scale );
194  n->set_alphaQCD( alphaQCD );
195  n->set_alphaQED( alphaQED );
196  return n;
197 }
198 /** Converts HepMC::Genevent to HepMC3::Genevent */
199 HepMC3::GenEvent* ConvertHepMCGenEvent_2to3(const HepMC::GenEvent& evt, std::shared_ptr<HepMC3::GenRunInfo> run )
200 {
204  if (evt.length_unit()==HepMC::Units::CM) lu=HepMC3::Units::CM;
205  if (evt.length_unit()==HepMC::Units::MM) lu=HepMC3::Units::MM;
206  if (evt.momentum_unit()==HepMC::Units::MEV) mu=HepMC3::Units::MEV;
207  if (evt.momentum_unit()==HepMC::Units::GEV) mu=HepMC3::Units::GEV;
208  n->set_units(mu,lu);
209 
210  n->set_run_info(run);
211  n->set_event_number(evt.event_number());
212  std::map<HepMC::GenVertex*,HepMC3::GenVertexPtr> vertexmap2to3;
213  for (auto v2=evt.vertices_begin(); v2!=evt.vertices_end(); v2++)
214  {
215  HepMC::FourVector pos2=(*v2)->position();
216  HepMC3::FourVector pos3=HepMC3::FourVector(pos2.x(),pos2.y(),pos2.pz(),pos2.t());
217  HepMC3::GenVertexPtr v3=std::make_shared<HepMC3::GenVertex>(pos3);
218  n->add_vertex(v3);
219  vertexmap2to3[*v2]=v3;
220  }
221 
222  std::map<HepMC::GenParticle*,HepMC3::GenParticlePtr> particlemap2to3;
223  for (auto p2=evt.particles_begin(); p2!=evt.particles_end(); p2++)
224  {
225 
226  HepMC::FourVector mom2=(*p2)->momentum();
227  HepMC3::FourVector mom3=HepMC3::FourVector(mom2.px(),mom2.py(),mom2.pz(),mom2.e());
228  HepMC3::GenParticlePtr p3= std::make_shared<HepMC3::GenParticle>(mom3,(*p2)->pdg_id(),(*p2)->status());
229  /// we set it always as there is no way to check if it is set
230  p3->set_generated_mass((*p2)->generated_mass());
231  particlemap2to3[*p2]=p3;
232 
233  auto v2production=(*p2)->production_vertex();
234  HepMC3::GenVertexPtr v3production;
235  if (v2production)
236  {
237  auto v3=vertexmap2to3.find(v2production);
238  if (v3!=vertexmap2to3.end()) v3production=vertexmap2to3[v2production];
239  }
240 
241  auto v2end=(*p2)->end_vertex();
242  HepMC3::GenVertexPtr v3end;
243  if (v2end)
244  {
245  auto v3=vertexmap2to3.find(v2end);
246  if (v3!=vertexmap2to3.end()) v3end=vertexmap2to3[v2end];
247  }
248 
249  if (v3production) v3production->add_particle_out(p3);
250  if (v3end) v3end->add_particle_in(p3);
251 
252  std::shared_ptr<HepMC3::DoubleAttribute> p3_theta=std::make_shared<HepMC3::DoubleAttribute>((*p2)->polarization().theta());
253 
254  std::shared_ptr<HepMC3::DoubleAttribute> p3_phi=std::make_shared<HepMC3::DoubleAttribute>((*p2)->polarization().phi());
255  if (std::abs(p3_phi->value())>std::numeric_limits<double>::epsilon() || std::abs(p3_theta->value())>std::numeric_limits<double>::epsilon())
256  {
257  p3->add_attribute("theta",p3_theta);
258  p3->add_attribute("phi",p3_theta);
259  }
260 
261  std::vector<int> flv;
262  auto fl=(*p2)->flow();
263  for (auto flel=fl.begin(); flel!=fl.end(); ++flel) flv.push_back((*flel).second);
264  std::shared_ptr<HepMC3::VectorIntAttribute> flows = std::make_shared<HepMC3::VectorIntAttribute>(flv);
265  if (flv.size())
266  p3->add_attribute("flows",flows);
267  }
268 
269 
270  std::shared_ptr<HepMC3::GenCrossSection> cs3 = std::make_shared<HepMC3::GenCrossSection>();
271  auto cs2=evt.cross_section();
272  if (cs2)cs3->set_cross_section(cs2->cross_section(),cs2->cross_section_error());
273  n->set_cross_section(cs3);
274 
275 
276  std::shared_ptr<HepMC3::GenPdfInfo> pdf3 = std::make_shared<HepMC3::GenPdfInfo>();
277  auto pdf2=evt.pdf_info();
278  if(pdf2)
279  {
280  pdf3->parton_id[0]=pdf2->id1();
281  pdf3->parton_id[1]=pdf2->id2();
282  pdf3->x[0]=pdf2->x1();
283  pdf3->x[1]=pdf2->x2();
284  pdf3->scale=pdf2->scalePDF();
285  pdf3->xf[0]=pdf2->pdf1();
286  pdf3->xf[1]=pdf2->pdf2();
287  pdf3->pdf_id[0]=pdf2->pdf_id1();
288  pdf3->pdf_id[1]=pdf2->pdf_id2();
289  n->set_pdf_info(pdf3);
290  }
291 
292  std::shared_ptr<HepMC3::GenHeavyIon> hi3 = std::make_shared<HepMC3::GenHeavyIon>();
293  auto hi2=evt.heavy_ion();
294  if(hi2)
295  {
296  hi3->Ncoll_hard=hi2->Ncoll_hard();
297  hi3->Npart_proj=hi2->Npart_proj();
298  hi3->Npart_targ=hi2->Npart_targ();
299  hi3->Ncoll=hi2->Ncoll();
300  hi3->spectator_neutrons=hi2->spectator_neutrons();
301  hi3->spectator_protons=hi2->spectator_protons();
302  hi3->N_Nwounded_collisions=hi2->N_Nwounded_collisions();
303  hi3->Nwounded_N_collisions=hi2->Nwounded_N_collisions();
304  hi3->Nwounded_Nwounded_collisions=hi2->Nwounded_Nwounded_collisions();
305  hi3->impact_parameter=hi2->impact_parameter();
306  hi3->event_plane_angle=hi2->event_plane_angle();
307  hi3->eccentricity=hi2->eccentricity();
308  hi3->sigma_inel_NN=hi2->sigma_inel_NN();
309 #ifdef HEPMC_HAS_CENTRALITY
310  hi3->centrality=hi2->centrality();
311 #endif
312  n->set_heavy_ion(hi3);
313  }
314  /** Yes, the desing is not always perfect */
315  std::stringstream ss;
316  evt.weights().print(ss);
317  std::string allweights=ss.str();
318  std::vector<std::string> wnames;
319  std::string token;
320  std::size_t pos=0;
321  for (;;)
322  {
323  std::size_t name_begin=allweights.find_first_not_of(" (\n",pos);
324  if (name_begin==std::string::npos) break;
325  std::size_t name_end=allweights.find_first_of(",",name_begin);
326  wnames.push_back(allweights.substr(name_begin,name_end-name_begin));
327  pos=allweights.find_first_of(")",name_end)+1;
328  }
329  n->run_info()->set_weight_names(wnames);
330  n->weights()=std::vector<double>(wnames.size(),1.0);
331  for (auto wn: wnames)
332  {
333  n->weight(wn)=evt.weights()[wn];
334  }
335 
336  std::shared_ptr<HepMC3::DoubleAttribute> A_event_scale=std::make_shared<HepMC3::DoubleAttribute>(evt.event_scale());
337  n->add_attribute("event_scale",A_event_scale);
338  std::shared_ptr<HepMC3::DoubleAttribute> A_alphaQED=std::make_shared<HepMC3::DoubleAttribute>(evt.alphaQED());
339  n->add_attribute("alphaQED",A_alphaQED);
340  std::shared_ptr<HepMC3::DoubleAttribute> A_alphaQCD=std::make_shared<HepMC3::DoubleAttribute>(evt.alphaQCD());
341  n->add_attribute("alphaQCD",A_alphaQCD);
342  std::shared_ptr<HepMC3::IntAttribute> A_signal_process_id=std::make_shared<HepMC3::IntAttribute>(evt.signal_process_id());
343  n->add_attribute("signal_process_id",A_signal_process_id);
344  std::shared_ptr<HepMC3::IntAttribute> A_mpi=std::make_shared<HepMC3::IntAttribute>(evt.mpi());
345  n->add_attribute("mpi",A_mpi);
346  std::shared_ptr<HepMC3::VectorLongIntAttribute> A_random_states= std::make_shared<HepMC3::VectorLongIntAttribute>( evt.random_states());
347  n->add_attribute("random_states",A_random_states);
348 
349  auto spv2=evt.signal_process_vertex();
350  if (vertexmap2to3.find(spv2)!=vertexmap2to3.end())
351  {
352  int signal_process_vertex_id3=vertexmap2to3[spv2]->id();
353  std::shared_ptr<HepMC3::IntAttribute> A_signal_process_vertex=std::make_shared<HepMC3::IntAttribute>(signal_process_vertex_id3);
354  n->add_attribute("signal_process_vertex",A_signal_process_vertex);
355  }
356  return n;
357 }
358 #endif
const std::vector< ConstGenVertexPtr > & vertices() const
Get list of vertices (const)
Definition: GenEvent.cc:43
void set_cross_section(GenCrossSectionPtr cs)
Set cross-section information.
Definition: GenEvent.h:166
void add_vertex(GenVertexPtr v)
Add vertex.
Definition: GenEvent.cc:97
double t() const
Time component of position/displacement.
Definition: FourVector.h:101
Definition of class GenParticle.
const std::vector< std::string > & weight_names() const
Definition: GenEvent.h:110
Definition of class GenVertex.
Attribute that holds a vector of integers of type int.
Definition: Attribute.h:1002
const Units::LengthUnit & length_unit() const
Get length unit.
Definition: GenEvent.h:142
double x() const
x-component of position/displacement
Definition: FourVector.h:80
int event_number() const
Get event number.
Definition: GenEvent.h:135
void add_attribute(const std::string &name, const std::shared_ptr< Attribute > &att, const int &id=0)
Add event attribute to event.
Definition: GenEvent.h:208
void set_cross_section(const double &xs, const double &xs_err, const long &n_acc=-1, const long &n_att=-1)
Set all fields.
std::vector< int > value() const
get the value associated to this Attribute.
Definition: Attribute.h:1028
LengthUnit
Position units.
Definition: Units.h:32
Stores additional information about PDFs.
Definition: GenPdfInfo.h:32
std::vector< ConstGenParticlePtr > beams() const
Vector of beam particles.
Definition: GenEvent.cc:420
MomentumUnit
Momentum units.
Definition: Units.h:29
const Units::MomentumUnit & momentum_unit() const
Get momentum unit.
Definition: GenEvent.h:140
double e() const
Energy component of momentum.
Definition: FourVector.h:130
Stores event-related information.
Definition: GenEvent.h:41
Generic 4-vector.
Definition: FourVector.h:35
Attribute that holds a real number as a double.
Definition: Attribute.h:242
void set_pdf_info(GenPdfInfoPtr pi)
Set PDF information.
Definition: GenEvent.h:159
double px() const
x-component of momentum
Definition: FourVector.h:109
HepMC3::GenEvent * ConvertHepMCGenEvent_2to3(const HepMC::GenEvent &evt, std::shared_ptr< HepMC3::GenRunInfo > run)
std::shared_ptr< T > attribute(const std::string &name, const int &id=0) const
Get attribute of type T.
Definition: GenEvent.h:389
void set_run_info(std::shared_ptr< GenRunInfo > run)
Set the GenRunInfo object by smart pointer.
Definition: GenEvent.h:128
Stores additional information about Heavy Ion generator.
Definition: GenHeavyIon.h:28
GenHeavyIon HeavyIon
Backward compatibility typedef.
Definition: GenHeavyIon.h:248
void set_units(Units::MomentumUnit new_momentum_unit, Units::LengthUnit new_length_unit)
Change event units Converts event from current units to new ones.
Definition: GenEvent.cc:395
HepMC::GenEvent * ConvertHepMCGenEvent_3to2(const HepMC3::GenEvent &evt)
Please note the HEPMC_HAS_CENTRALITY should be defined externaly.
double y() const
y-component of position/displacement
Definition: FourVector.h:87
const std::vector< double > & weights() const
Get event weight values as a vector.
Definition: GenEvent.h:86
double pz() const
z-component of momentum
Definition: FourVector.h:123
Definition of class GenEvent.
double py() const
y-component of momentum
Definition: FourVector.h:116
double weight(const unsigned long &index=0) const
Definition: GenEvent.h:91
const std::vector< ConstGenParticlePtr > & particles() const
Get list of particles (const)
Definition: GenEvent.cc:39
void set_heavy_ion(GenHeavyIonPtr hi)
Set heavy ion generator additional information.
Definition: GenEvent.h:152
int value() const
get the value associated to this Attribute.
Definition: Attribute.h:180
std::shared_ptr< GenRunInfo > run_info() const
Get a pointer to the the GenRunInfo object.
Definition: GenEvent.h:124
Attribute that holds an Integer implemented as an int.
Definition: Attribute.h:158
Feature< Feature_type > abs(const Feature< Feature_type > &input)
Obtain the absolute value of a Feature. This works as you&#39;d expect. If foo is a valid Feature...
Definition: Feature.h:316
Stores additional information about cross-section.
void set_event_number(const int &num)
Set event number.
Definition: GenEvent.h:137