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 = std::make_shared<GenParticle>( FourVector( 0.0, 0.0, 7000.0, 7000.0 ),2212, 3 );
51 GenParticlePtr p2 = std::make_shared<GenParticle>( FourVector( 0.750, -1.569, 32.191, 32.238), 1, 3 );
52 GenParticlePtr p3 = std::make_shared<GenParticle>( FourVector( 0.0, 0.0, -7000.0, 7000.0 ),2212, 3 );
53 GenParticlePtr p4 = std::make_shared<GenParticle>( FourVector(-3.047,-19.0, -54.629, 57.920), -2, 3 );
54
55 GenVertexPtr v1 = std::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 = std::make_shared<GenVertex>();
64 v2->add_particle_in (p3);
65 v2->add_particle_out(p4);
66 evt.add_vertex(v2);
67
68 GenVertexPtr v3 = std::make_shared<GenVertex>();
69 v3->add_particle_in(p2);
70 v3->add_particle_in(p4);
71 evt.add_vertex(v3);
72
73 GenParticlePtr p5 = std::make_shared<GenParticle>( FourVector(-3.813, 0.113, -1.833, 4.233), 22, 1 );
74 GenParticlePtr p6 = std::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 =std:: make_shared<GenVertex>();
80 v4->add_particle_in (p6);
81 evt.add_vertex(v4);
82
83 GenParticlePtr p7 = std::make_shared<GenParticle>( FourVector(-2.445, 28.816, 6.082,29.552), 1, 1 );
84 GenParticlePtr p8 = std::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 std::shared_ptr<GenPdfInfo> pdf_info = std::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 std::shared_ptr<GenHeavyIon> heavy_ion = std::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 std::shared_ptr<GenCrossSection> cross_section = std::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 std::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 std::shared_ptr<Attribute> tool1 = std::make_shared<IntAttribute>(1);
176 std::shared_ptr<Attribute> tool999 = std::make_shared<IntAttribute>(999);
177 std::shared_ptr<Attribute> test_attribute = std::make_shared<StringAttribute>("test attribute");
178 std::shared_ptr<Attribute> test_attribute2 = std::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 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:41
HepMC3 main namespace.
int main(int argc, char **argv)