HepMC3 event record library
basic_tree.cc
1 // -*- C++ -*-
2 //
3 // This file is part of HepMC
4 // Copyright (C) 2014-2019 The HepMC collaboration (see AUTHORS for details)
5 //
6 /// @example basic_tree.cc
7 /// @brief Basic example of building HepMC3 tree by hand
8 ///
9 /// Based on HepMC2/examples/example_BuildEventFromScratch.cc
10 
11 #include "HepMC3/GenEvent.h"
12 #include "HepMC3/GenVertex.h"
13 #include "HepMC3/GenParticle.h"
14 #include "HepMC3/Print.h"
15 #include "HepMC3/Selector.h"
16 #include "HepMC3/Relatives.h"
17 
18 using namespace HepMC3;
19 
20 
21 /** Main program */
22 int main() {
23  //
24  // In this example we will place the following event into HepMC "by hand"
25  //
26  // name status pdg_id parent Px Py Pz Energy Mass
27  // 1 !p+! 3 2212 0,0 0.000 0.000 7000.000 7000.000 0.938
28  // 3 !p+! 3 2212 0,0 0.000 0.000-7000.000 7000.000 0.938
29  //=========================================================================
30  // 2 !d! 3 1 1,1 0.750 -1.569 32.191 32.238 0.000
31  // 4 !u~! 3 -2 2,2 -3.047 -19.000 -54.629 57.920 0.000
32  // 5 !W-! 3 -24 1,2 1.517 -20.68 -20.605 85.925 80.799
33  // 6 !gamma! 1 22 1,2 -3.813 0.113 -1.833 4.233 0.000
34  // 7 !d! 1 1 5,5 -2.445 28.816 6.082 29.552 0.010
35  // 8 !u~! 1 -2 5,5 3.962 -49.498 -26.687 56.373 0.006
36 
37  // now we build the graph, which will looks like
38  // p7 #
39  // p1 / #
40  // \v1__p3 p5---v4 #
41  // \_v3_/ \ #
42  // / \ p8 #
43  // v2__p4 \ #
44  // / p6 #
45  // p2 #
46  // #
47  GenEvent evt(Units::GEV,Units::MM);
48 
49  // px py pz e pdgid status
50  GenParticlePtr p1 = make_shared<GenParticle>( FourVector( 0.0, 0.0, 7000.0, 7000.0 ),2212, 3 );
51  GenParticlePtr p2 = make_shared<GenParticle>( FourVector( 0.750, -1.569, 32.191, 32.238), 1, 3 );
52  GenParticlePtr p3 = make_shared<GenParticle>( FourVector( 0.0, 0.0, -7000.0, 7000.0 ),2212, 3 );
53  GenParticlePtr p4 = make_shared<GenParticle>( FourVector(-3.047,-19.0, -54.629, 57.920), -2, 3 );
54 
55  GenVertexPtr v1 = make_shared<GenVertex>();
56  v1->add_particle_in (p1);
57  v1->add_particle_out(p2);
58  evt.add_vertex(v1);
59 
60  // Set vertex status if needed
61  v1->set_status(4);
62 
63  GenVertexPtr v2 = make_shared<GenVertex>();
64  v2->add_particle_in (p3);
65  v2->add_particle_out(p4);
66  evt.add_vertex(v2);
67 
68  GenVertexPtr v3 = make_shared<GenVertex>();
69  v3->add_particle_in(p2);
70  v3->add_particle_in(p4);
71  evt.add_vertex(v3);
72 
73  GenParticlePtr p5 = make_shared<GenParticle>( FourVector(-3.813, 0.113, -1.833, 4.233), 22, 1 );
74  GenParticlePtr p6 = make_shared<GenParticle>( FourVector( 1.517,-20.68, -20.605,85.925), -24, 3 );
75 
76  v3->add_particle_out(p5);
77  v3->add_particle_out(p6);
78 
79  GenVertexPtr v4 = make_shared<GenVertex>();
80  v4->add_particle_in (p6);
81  evt.add_vertex(v4);
82 
83  GenParticlePtr p7 = make_shared<GenParticle>( FourVector(-2.445, 28.816, 6.082,29.552), 1, 1 );
84  GenParticlePtr p8 = make_shared<GenParticle>( FourVector( 3.962,-49.498,-26.687,56.373), -2, 1 );
85 
86  v4->add_particle_out(p7);
87  v4->add_particle_out(p8);
88 
89 /* Unfortunatelly this code is not portable. TODO: make it portable.
90  //
91  // Example of use of the search engine
92  //
93 
94  // 1)
95  std::cout << std::endl << "Find all stable particles: " << std::endl;
96 
97  for(ConstGenParticlePtr p: applyFilter(StandardSelector::STATUS == 1, evt.particles())){
98  Print::line(p);
99  }
100 
101  // 2)
102  std::cout <<std::endl << "Find all ancestors of particle with id " << p5->id() << ": " << std::endl;
103 
104  for(ConstGenParticlePtr p: Relatives::ANCESTORS(p5)){
105  Print::line(p);
106  }
107 
108  // 3)
109  std::cout <<std::endl << "Find stable descendants of particle with id " << p4->id() << ": " << std::endl;
110  std::cout<<"We check both for STATUS == 1 (equivalent of IS_STABLE) and no end vertex, just to be safe" << std::endl;
111 
112  Filter has_end_vtx = [](ConstGenParticlePtr input)->bool{return (bool)input->end_vertex();};
113 
114  vector<GenParticlePtr> results3 = applyFilter(StandardSelector::STATUS==1 && has_end_vtx, Relatives::DESCENDANTS(p4));
115  for(ConstGenParticlePtr p: results3){
116  Print::line(p);
117  }
118 
119  // 3b)
120  std::cout << std::endl << "Narrow down results of previous search to quarks only: " << std::endl;
121 
122  // note the use of abs to obtain the absolute value of pdg_id :)
123  for(ConstGenParticlePtr p: applyFilter( *abs(StandardSelector::PDG_ID) <= 6, results3)){
124  Print::line(p);
125  }
126 */
127  //
128  // Example of adding event attributes
129  //
130  shared_ptr<GenPdfInfo> pdf_info = make_shared<GenPdfInfo>();
131  evt.add_attribute("GenPdfInfo",pdf_info);
132 
133  pdf_info->set(1,2,3.4,5.6,7.8,9.0,1.2,3,4);
134 
135  shared_ptr<GenHeavyIon> heavy_ion = make_shared<GenHeavyIon>();
136  evt.add_attribute("GenHeavyIon",heavy_ion);
137 
138  heavy_ion->set( 1,2,3,4,5,6,7,8,9,0.1,2.3,4.5,6.7);
139 
140  shared_ptr<GenCrossSection> cross_section = make_shared<GenCrossSection>();
141  evt.add_attribute("GenCrossSection",cross_section);
142 
143  cross_section->set_cross_section(1.2,3.4);
144 
145  //
146  // Example of manipulating the attributes
147  //
148 
149  std::cout << std::endl << " Manipulating attributes:" << std::endl;
150 
151  // get attribute
152  shared_ptr<GenCrossSection> cs = evt.attribute<GenCrossSection>("GenCrossSection");
153 
154  // if attribute exists - do something with it
155  if(cs) {
156  cs->set_cross_section(-1.0,0.0);
157  Print::line(cs);
158  }
159  else std::cout << "Problem accessing attribute!" <<std::endl;
160 
161  // remove attribute
162  evt.remove_attribute("GenCrossSection");
163  evt.remove_attribute("GenCrossSection"); // This call will do nothing
164 
165  // now this should be null
166  cs = evt.attribute<GenCrossSection>("GenCrossSection");
167 
168  if(!cs)std::cout << "Successfully removed attribute" <<std::endl;
169  else std::cout << "Problem removing attribute!" <<std::endl;
170 
171  //
172  // Example of adding attributes and finding particles with attributes
173  //
174 
175  shared_ptr<Attribute> tool1 = make_shared<IntAttribute>(1);
176  shared_ptr<Attribute> tool999 = make_shared<IntAttribute>(999);
177  shared_ptr<Attribute> test_attribute = make_shared<StringAttribute>("test attribute");
178  shared_ptr<Attribute> test_attribute2 = make_shared<StringAttribute>("test attribute2");
179 
180  p2->add_attribute( "tool" , tool1 );
181  p2->add_attribute( "other" , test_attribute );
182 
183  p4->add_attribute( "tool" , tool1 );
184 
185  p6->add_attribute( "tool" , tool999 );
186  p6->add_attribute( "other" , test_attribute2 );
187 
188  v3->add_attribute( "vtx_att" , test_attribute );
189  v4->add_attribute( "vtx_att" , test_attribute2 );
190 /* TODO: Make this code portable
191 
192  std::cout << std::endl << "Find all particles with attribute 'tool' "<< std::endl;
193  std::cout << "(should return particles 2,4,6):" << std::endl;
194 
195  /// @todo can we add some utility funcs to simplify creation of Features from Attributes and check they exist.
196  /// Features and Attributes are quite similar concepts anyway, can they be unified (but Features can also be
197  /// non-attribute-like e.g. pT, rapidity or any quantity it is possible to obtain from a particle)
198 
199  for(ConstGenParticlePtr p: applyFilter(Selector::ATTRIBUTE("tool"), evt.particles())){
200  Print::line(p);
201  }
202 
203  std::cout <<std::endl << "Find all particles with attribute 'tool' equal 1 "<< std::endl;
204  std::cout << "(should return particles 2,4):" <<std::endl;
205 
206  for(ConstGenParticlePtr p: applyFilter(Selector::ATTRIBUTE("tool") && Selector::ATTRIBUTE("tool") == tool1, evt.particles())){
207  Print::line(p);
208  }
209 
210  std::cout << std::endl << "Find all particles with a string attribute 'other' equal 'test attribute' "<< std::endl;
211  std::cout << "(should return particle 2):" << std::endl;
212 
213 
214  for(ConstGenParticlePtr p: applyFilter(Selector::ATTRIBUTE("other") && Selector::ATTRIBUTE("other") == "test_attribute", evt.particles())){
215  Print::line(p);
216  }
217 */
218 
219  std::cout << std::endl << "Offsetting event position by 5,5,5,5" << std::endl;
220 
221  evt.shift_position_by( FourVector(5,5,5,5) );
222 
223  Print::listing(evt);
224 
225  std::cout << std::endl << "Printing full content of the GenEvent object " << std::endl
226  << "(including particles and vertices in one-line format):" << std::endl << std::endl;
227 
228  Print::content(evt);
229 
230  std::cout <<std::endl << "Now: removing particle with id 6 and printing again:" <<std::endl <<std::endl;
231  evt.remove_particle(p6);
232 
233  Print::listing(evt);
234  Print::content(evt);
235 
236  std::cout <<std::endl << "Now: removing beam particles, leaving an empty event" <<std::endl <<std::endl;
237  evt.remove_particles( evt.beams() );
238 
239  Print::listing(evt);
240  Print::content(evt);
241  return 0;
242 }
definition of /b Selector class
Definition of class GenParticle.
Definition of class GenVertex.
static void listing(std::ostream &os, const GenEvent &event, unsigned short precision=2)
Print event in listing (HepMC2) format.
Definition: Print.cc:51
void set_cross_section(const double &xs, const double &xs_err, const long &n_acc=-1, const long &n_att=-1)
Set all fields.
Stores event-related information.
Definition: GenEvent.h:41
Generic 4-vector.
Definition: FourVector.h:35
Defines helper classes to extract relatives of an input GenParticle or GenVertex. ...
static void line(std::ostream &os, const GenEvent &event, bool attributes=false)
Print one-line info.
Definition: Print.cc:203
int main(int argc, char **argv)
Definition of class GenEvent.
Stores additional information about cross-section.
static void content(std::ostream &os, const GenEvent &event)
Print content of all GenEvent containers.
Definition: Print.cc:18