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  //
90  // Example of use of the search engine
91  //
92 
93  // 1)
94  std::cout << std::endl << "Find all stable particles: " << std::endl;
95 
96  for(ConstGenParticlePtr p: applyFilter(Selector::STATUS == 1, evt.particles())){
97  Print::line(p);
98  }
99 
100  // 2)
101  std::cout <<std::endl << "Find all ancestors of particle with id " << p5->id() << ": " << std::endl;
102 
103  for(ConstGenParticlePtr p: Relatives::ANCESTORS(p5)){
104  Print::line(p);
105  }
106 
107  // 3)
108  std::cout <<std::endl << "Find stable descendants of particle with id " << p4->id() << ": " << std::endl;
109  std::cout<<"We check both for STATUS == 1 (equivalent of IS_STABLE) and no end vertex, just to be safe" << std::endl;
110 
111  Filter has_end_vtx = [](ConstGenParticlePtr input)->bool{return (bool)input->end_vertex();};
112 
113  vector<GenParticlePtr> results3 = applyFilter(Selector::STATUS==1 && has_end_vtx, Relatives::DESCENDANTS(p4));
114  for(ConstGenParticlePtr p: results3){
115  Print::line(p);
116  }
117 
118  // 3b)
119  std::cout << std::endl << "Narrow down results of previous search to quarks only: " << std::endl;
120 
121  // note the use of abs to obtain the absolute value of pdg_id :)
122  for(ConstGenParticlePtr p: applyFilter( *abs(Selector::PDG_ID) <= 6, results3)){
123  Print::line(p);
124  }
125 
126  //
127  // Example of adding event attributes
128  //
129  shared_ptr<GenPdfInfo> pdf_info = make_shared<GenPdfInfo>();
130  evt.add_attribute("GenPdfInfo",pdf_info);
131 
132  pdf_info->set(1,2,3.4,5.6,7.8,9.0,1.2,3,4);
133 
134  shared_ptr<GenHeavyIon> heavy_ion = make_shared<GenHeavyIon>();
135  evt.add_attribute("GenHeavyIon",heavy_ion);
136 
137  heavy_ion->set( 1,2,3,4,5,6,7,8,9,0.1,2.3,4.5,6.7);
138 
139  shared_ptr<GenCrossSection> cross_section = make_shared<GenCrossSection>();
140  evt.add_attribute("GenCrossSection",cross_section);
141 
142  cross_section->set_cross_section(1.2,3.4);
143 
144  //
145  // Example of manipulating the attributes
146  //
147 
148  std::cout << std::endl << " Manipulating attributes:" << std::endl;
149 
150  // get attribute
151  shared_ptr<GenCrossSection> cs = evt.attribute<GenCrossSection>("GenCrossSection");
152 
153  // if attribute exists - do something with it
154  if(cs) {
155  cs->set_cross_section(-1.0,0.0);
156  Print::line(cs);
157  }
158  else std::cout << "Problem accessing attribute!" <<std::endl;
159 
160  // remove attribute
161  evt.remove_attribute("GenCrossSection");
162  evt.remove_attribute("GenCrossSection"); // This call will do nothing
163 
164  // now this should be null
165  cs = evt.attribute<GenCrossSection>("GenCrossSection");
166 
167  if(!cs)std::cout << "Successfully removed attribute" <<std::endl;
168  else std::cout << "Problem removing attribute!" <<std::endl;
169 
170  //
171  // Example of adding attributes and finding particles with attributes
172  //
173 
174  shared_ptr<Attribute> tool1 = make_shared<IntAttribute>(1);
175  shared_ptr<Attribute> tool999 = make_shared<IntAttribute>(999);
176  shared_ptr<Attribute> test_attribute = make_shared<StringAttribute>("test attribute");
177  shared_ptr<Attribute> test_attribute2 = make_shared<StringAttribute>("test attribute2");
178 
179  p2->add_attribute( "tool" , tool1 );
180  p2->add_attribute( "other" , test_attribute );
181 
182  p4->add_attribute( "tool" , tool1 );
183 
184  p6->add_attribute( "tool" , tool999 );
185  p6->add_attribute( "other" , test_attribute2 );
186 
187  v3->add_attribute( "vtx_att" , test_attribute );
188  v4->add_attribute( "vtx_att" , test_attribute2 );
189 
190  std::cout << std::endl << "Find all particles with attribute 'tool' "<< std::endl;
191  std::cout << "(should return particles 2,4,6):" << std::endl;
192 
193  /// @todo can we add some utility funcs to simplify creation of Features from Attributes and check they exist.
194  /// Features and Attributes are quite similar concepts anyway, can they be unified (but Features can also be
195  /// non-attribute-like e.g. pT, rapidity or any quantity it is possible to obtain from a particle)
196 
197  for(ConstGenParticlePtr p: applyFilter(Selector::ATTRIBUTE("tool"), evt.particles())){
198  Print::line(p);
199  }
200 
201  std::cout <<std::endl << "Find all particles with attribute 'tool' equal 1 "<< std::endl;
202  std::cout << "(should return particles 2,4):" <<std::endl;
203 
204  for(ConstGenParticlePtr p: applyFilter(Selector::ATTRIBUTE("tool") && Selector::ATTRIBUTE("tool") == tool1, evt.particles())){
205  Print::line(p);
206  }
207 
208  std::cout << std::endl << "Find all particles with a string attribute 'other' equal 'test attribute' "<< std::endl;
209  std::cout << "(should return particle 2):" << std::endl;
210 
211 
212  for(ConstGenParticlePtr p: applyFilter(Selector::ATTRIBUTE("other") && Selector::ATTRIBUTE("other") == "test_attribute", evt.particles())){
213  Print::line(p);
214  }
215 
216  std::cout << std::endl << "Offsetting event position by 5,5,5,5" << std::endl;
217 
218  evt.shift_position_by( FourVector(5,5,5,5) );
219 
220  Print::listing(evt);
221 
222  std::cout << std::endl << "Printing full content of the GenEvent object " << std::endl
223  << "(including particles and vertices in one-line format):" << std::endl << std::endl;
224 
225  Print::content(evt);
226 
227  std::cout <<std::endl << "Now: removing particle with id 6 and printing again:" <<std::endl <<std::endl;
228  evt.remove_particle(p6);
229 
230  Print::listing(evt);
231  Print::content(evt);
232 
233  std::cout <<std::endl << "Now: removing beam particles, leaving an empty event" <<std::endl <<std::endl;
234  evt.remove_particles( evt.beams() );
235 
236  Print::listing(evt);
237  Print::content(evt);
238  return 0;
239 }
HepMC3 main namespace.
Definition: ReaderGZ.h:28
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:42
Generic 4-vector.
Definition: FourVector.h:35
Defines helper classes to extract relatives of an input GenParticle or GenVertex. ...
std::function< bool(ConstGenParticlePtr)> Filter
type of Filter
Definition: Filter.h:17
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.
vector< GenParticlePtr > applyFilter(const Filter &filter, const vector< GenParticlePtr > &particles)
Apply a Filter to a list of GenParticles Returns a vector of GenParticles that satisfy the Filter...
Definition: Filter.h:21
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.
static void content(std::ostream &os, const GenEvent &event)
Print content of all GenEvent containers.
Definition: Print.cc:18