HepMC event record
basic_tree.cc
1 // -*- C++ -*-
2 //
3 // This file is part of HepMC
4 // Copyright (C) 2014 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 "HepMC/GenEvent.h"
12 #include "HepMC/GenVertex.h"
13 #include "HepMC/GenParticle.h"
14 #include "HepMC/Search/FindParticles.h"
15 #include "HepMC/Print.h"
16 
17 using namespace HepMC;
18 using namespace std;
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<HepMC::GenParticle>( FourVector( 0.0, 0.0, 7000.0, 7000.0 ),2212, 3 );
51  GenParticlePtr p2 = make_shared<HepMC::GenParticle>( FourVector( 0.750, -1.569, 32.191, 32.238), 1, 3 );
52  GenParticlePtr p3 = make_shared<HepMC::GenParticle>( FourVector( 0.0, 0.0, -7000.0, 7000.0 ),2212, 3 );
53  GenParticlePtr p4 = make_shared<HepMC::GenParticle>( FourVector(-3.047,-19.0, -54.629, 57.920), -2, 3 );
54 
55  GenVertexPtr v1 = make_shared<HepMC::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<HepMC::GenVertex>();
64  v2->add_particle_in (p3);
65  v2->add_particle_out(p4);
66  evt.add_vertex(v2);
67 
68  GenVertexPtr v3 = make_shared<HepMC::GenVertex>();
69  v3->add_particle_in(p2);
70  v3->add_particle_in(p4);
71  evt.add_vertex(v3);
72 
73  GenParticlePtr p5 = make_shared<HepMC::GenParticle>( FourVector(-3.813, 0.113, -1.833, 4.233), 22, 1 );
74  GenParticlePtr p6 = make_shared<HepMC::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<HepMC::GenVertex>();
80  v4->add_particle_in (p6);
81  evt.add_vertex(v4);
82 
83  GenParticlePtr p7 = make_shared<HepMC::GenParticle>( FourVector(-2.445, 28.816, 6.082,29.552), 1, 1 );
84  GenParticlePtr p8 = make_shared<HepMC::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  //
90  // Example of use of the search engine
91  //
92 
93  // 1)
94  cout << endl << "Find all stable particles: " << endl;
95 
96  FindParticles search(evt, FIND_ALL, IS_STABLE);
97 
98  FOREACH( const GenParticlePtr &p, search.results() ) {
99  Print::line(p);
100  }
101 
102  // 2)
103  cout << endl << "Find all ancestors of particle with id " << p5->id() << ": " << endl;
104 
105  FindParticles search2(p5, FIND_ALL_ANCESTORS);
106 
107  FOREACH( const GenParticlePtr &p, search2.results() ) {
108  Print::line(p);
109  }
110 
111  // 3)
112  cout << endl << "Find stable descendants of particle with id " << p4->id() << ": " << endl;
113  cout<<"We check both for STATUS == 1 (equivalent of IS_STABLE) and no end vertex, just to be safe" << endl;
114 
115  FindParticles search3(p4, FIND_ALL_DESCENDANTS, STATUS == 1 && !HAS_END_VERTEX);
116 
117  FOREACH( const GenParticlePtr &p, search3.results() ) {
118  Print::line(p);
119  }
120 
121  // 3b)
122  cout << endl << "Narrow down results of previous search to quarks only: " << endl;
123 
124  search3.narrow_down( PDG_ID >= -6 && PDG_ID <= 6 );
125 
126  FOREACH( const GenParticlePtr &p, search3.results() ) {
127  Print::line(p);
128  }
129 
130  //
131  // Example of adding event attributes
132  //
133  shared_ptr<GenPdfInfo> pdf_info = make_shared<GenPdfInfo>();
134  evt.add_attribute("GenPdfInfo",pdf_info);
135 
136  pdf_info->set(1,2,3.4,5.6,7.8,9.0,1.2,3,4);
137 
138  shared_ptr<GenHeavyIon> heavy_ion = make_shared<GenHeavyIon>();
139  evt.add_attribute("GenHeavyIon",heavy_ion);
140 
141  heavy_ion->set( 1,2,3,4,5,6,7,8,9,0.1,2.3,4.5,6.7);
142 
143  shared_ptr<GenCrossSection> cross_section = make_shared<GenCrossSection>();
144  evt.add_attribute("GenCrossSection",cross_section);
145 
146  cross_section->set_cross_section(1.2,3.4);
147 
148  //
149  // Example of manipulating the attributes
150  //
151 
152  cout << endl << " Manipulating attributes:" << endl;
153 
154  // get attribute
155  shared_ptr<GenCrossSection> cs = evt.attribute<GenCrossSection>("GenCrossSection");
156 
157  // if attribute exists - do something with it
158  if(cs) {
159  cs->set_cross_section(-1.0,0.0);
160  Print::line(cs);
161  }
162  else cout << "Problem accessing attribute!" << endl;
163 
164  // remove attribute
165  evt.remove_attribute("GenCrossSection");
166  evt.remove_attribute("GenCrossSection"); // This call will do nothing
167 
168  // now this should be null
169  cs = evt.attribute<GenCrossSection>("GenCrossSection");
170 
171  if(!cs) cout << "Successfully removed attribute" << endl;
172  else cout << "Problem removing attribute!" << endl;
173 
174  //
175  // Example of adding attributes and finding particles with attributes
176  //
177 
178  shared_ptr<Attribute> tool1 = make_shared<IntAttribute>(1);
179  shared_ptr<Attribute> tool999 = make_shared<IntAttribute>(999);
180  shared_ptr<Attribute> test_attribute = make_shared<StringAttribute>("test attribute");
181  shared_ptr<Attribute> test_attribute2 = make_shared<StringAttribute>("test attribute2");
182 
183  p2->add_attribute( "tool" , tool1 );
184  p2->add_attribute( "other" , test_attribute );
185 
186  p4->add_attribute( "tool" , tool1 );
187 
188  p6->add_attribute( "tool" , tool999 );
189  p6->add_attribute( "other" , test_attribute2 );
190 
191  v3->add_attribute( "vtx_att" , test_attribute );
192  v4->add_attribute( "vtx_att" , test_attribute2 );
193 
194  cout << endl << "Find all particles with attribute 'tool' "<< endl;
195  cout << "(should return particles 2,4,6):" << endl;
196 
197  FindParticles search_attributes(evt, FIND_ALL, ATTRIBUTE("tool") );
198 
199  FOREACH( const GenParticlePtr &p, search_attributes.results() ) {
200  Print::line(p);
201  }
202 
203  cout << endl << "Find all particles with attribute 'tool' equal 1 "<< endl;
204  cout << "(should return particles 2,4):" << endl;
205 
206  FindParticles search_attributes2(evt, FIND_ALL, ATTRIBUTE("tool") == tool1 );
207 
208  FOREACH( const GenParticlePtr &p, search_attributes2.results() ) {
209  Print::line(p);
210  }
211 
212  cout << endl << "Find all particles with a string attribute 'other' equal 'test attribute' "<< endl;
213  cout << "(should return particle 2):" << endl;
214 
215  FindParticles search_attributes3(evt, FIND_ALL, ATTRIBUTE("other") == "test attribute" );
216 
217  FOREACH( const GenParticlePtr &p, search_attributes3.results() ) {
218  Print::line(p);
219  }
220 
221  cout << endl << "Offsetting event position by 5,5,5,5" << endl;
222 
223  evt.shift_position_by( FourVector(5,5,5,5) );
224 
225  Print::listing(evt);
226 
227  cout << endl << "Printing full content of the GenEvent object " << endl
228  << "(including particles and vertices in one-line format):" << endl << endl;
229 
230  Print::content(evt);
231 
232  cout << endl << "Now: removing particle with id 6 and printing again:" << endl << endl;
233  evt.remove_particle(p6);
234 
235  Print::listing(evt);
236  Print::content(evt);
237 
238  cout << endl << "Now: removing beam particles, leaving an empty event" << endl << endl;
239  evt.remove_particles( evt.beams() );
240 
241  Print::listing(evt);
242  Print::content(evt);
243  return 0;
244 }
static void listing(const GenEvent &event, unsigned short precision=2)
Print event in listing (HepMC2) format.
Definition: Print.cc:57
Stores additional information about cross-section.
STL namespace.
static void line(const GenEvent &event)
Print one-line info.
static const Filter HAS_END_VERTEX
Filter for checking if HepMC::GenParticle::end_vertex() != NULL.
Stores event-related information.
static const Filter IS_STABLE
Filter for checking if particle is stable.
static const FilterBase PDG_ID
Filter base for filtering GenParticle::pid()
static void content(const GenEvent &event)
Print content of all GenEvent containers.
Definition: Print.cc:20
int main(int argc, char **argv)
Definition of template class SmartPointer.
static const FilterBase STATUS
Filter base for filtering GenParticle::status()
void set_cross_section(double xs, double xs_err, long n_acc=-1, long n_att=-1)
Set all fields.