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
18using namespace HepMC3;
19
20
21/** Main program */
22int 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}
Definition of class GenEvent.
Definition of class GenParticle.
Definition of class GenVertex.
Definition of static class Print.
Defines helper classes to extract relatives of an input GenParticle or GenVertex.
definition of /b Selector class
Generic 4-vector.
Definition: FourVector.h:35
Stores additional information about cross-section.
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
HepMC3 main namespace.
Definition: ReaderGZ.h:28
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
std::function< bool(ConstGenParticlePtr)> Filter
type of Filter
Definition: Filter.h:17
Feature< Feature_type > abs(const Feature< Feature_type > &input)
Obtain the absolute value of a Feature. This works as you'd expect. If foo is a valid Feature,...
Definition: Feature.h:316
int main(int argc, char **argv)