1 // -*- C++ -*-
2 //
3 // This file is part of HepMC
4 // Copyright (C) 2014-2019 The HepMC collaboration (see AUTHORS for details)
5 //
6 #include <iostream>
7 #include <fstream>
8 #include <vector>
9
10 #include "HepMC3/Attribute.h"
11 #include "HepMC3/GenEvent.h"
12 #include "HepMC3/GenVertex.h"
13 #include "HepMC3/GenParticle.h"
14 #include "HepMC3/WriterAscii.h"
15 #include "HepMC3/WriterAsciiHepMC2.h"
16 #include "HepMC3/ReaderAscii.h"
17 #include "HepMC3/ReaderAsciiHepMC2.h"
18 #include "HepMC3/Print.h"
19 #ifndef M_PI
20 #define M_PI 3.14159265358979323846264338327950288
21 #endif
22 #include "HepMC3TestUtils.h"
23 using namespace HepMC3;
main()24 int main()
25 {
26 //
27 // In this example we will place the following event into HepMC "by hand"
28 //
29 // name status pdg_id parent Px Py Pz Energy Mass
30 // 1 !p+! 3 2212 0,0 0.000 0.000 7000.000 7000.000 0.938
31 // 2 !p+! 3 2212 0,0 0.000 0.000-7000.000 7000.000 0.938
32 //=========================================================================
33 // 3 !d! 3 1 1,1 0.750 -1.569 32.191 32.238 0.000
34 // 4 !u~! 3 -2 2,2 -3.047 -19.000 -54.629 57.920 0.000
35 // 5 !W-! 3 -24 1,2 1.517 -20.68 -20.605 85.925 80.799
36 // 6 !gamma! 1 22 1,2 -3.813 0.113 -1.833 4.233 0.000
37 // 7 !d! 1 1 5,5 -2.445 28.816 6.082 29.552 0.010
38 // 8 !u~! 1 -2 5,5 3.962 -49.498 -26.687 56.373 0.006
39 // 9 !gamma! 3 22 3,4 0.000 0.000 0.000 0.000 0.000
40
41
42 // declare several WriterAscii instances for comparison
43 WriterAscii xout1("testLoops1.out");
44 // output in old format
45 WriterAsciiHepMC2 xout2( "testLoops2.out" );
46
47 // build the graph, which will look like
48 // p7 #
49 // p1 / #
50 // \v1__p3 p5---v4 #
51 // \_v3_/ \ #
52 // / |\ p8 #
53 // v2__p4 | \ #
54 // / \ / p6 #
55 // p2 \p9_/ #
56 //
57 // define a flow pattern as p1 -> p3 -> p6
58 // and p2 -> p4 -> p5
59 //
60
61 // First create the event container, with Signal Process 20, event number 1
62 //
63 GenEvent evt(Units::GEV,Units::MM);
64 evt.set_event_number(1);
65 evt.add_attribute("signal_process_id", std::make_shared<IntAttribute>(20));
66 // create vertex 1
67 GenVertexPtr v1=std::make_shared<GenVertex>();
68 evt.add_vertex( v1 );
69 GenParticlePtr p1=std::make_shared<GenParticle>( FourVector(0,0,7000,7000),2212, 3 );
70 v1->add_particle_in( p1 );
71 p1->add_attribute("flow1", std::make_shared<IntAttribute>(231));
72 p1->add_attribute("flow1", std::make_shared<IntAttribute>(231));
73 p1->add_attribute("theta", std::make_shared<DoubleAttribute>(std::rand()/double(RAND_MAX)*M_PI));
74 p1->add_attribute("phi", std::make_shared<DoubleAttribute>(std::rand()/double(RAND_MAX)*M_PI*2));
75
76 GenVertexPtr v2=std::make_shared<GenVertex>();
77 evt.add_vertex( v2 );
78 GenParticlePtr p2=std::make_shared<GenParticle>( FourVector(0,0,-7000,7000),2212, 3 );
79 v2->add_particle_in( p2 );
80 p2->add_attribute("flow1", std::make_shared<IntAttribute>(243));
81 p2->add_attribute("theta", std::make_shared<DoubleAttribute>(std::rand()/double(RAND_MAX)*M_PI));
82 p2->add_attribute("phi", std::make_shared<DoubleAttribute>(std::rand()/double(RAND_MAX)*M_PI*2));
83 //
84 // create the outgoing particles of v1 and v2
85 GenParticlePtr p3=std::make_shared<GenParticle>( FourVector(.750,-1.569,32.191,32.238),1, 3 );
86 v1->add_particle_out( p3 );
87 p3->add_attribute("flow1", std::make_shared<IntAttribute>(231));
88 p3->add_attribute("theta", std::make_shared<DoubleAttribute>(std::rand()/double(RAND_MAX)*M_PI));
89 p3->add_attribute("phi", std::make_shared<DoubleAttribute>(std::rand()/double(RAND_MAX)*M_PI*2));
90
91 GenParticlePtr p4=std::make_shared<GenParticle>( FourVector(-3.047,-19.,-54.629,57.920),-2, 3 );
92 v2->add_particle_out( p4 );
93 p4->add_attribute("flow1", std::make_shared<IntAttribute>(243));
94 p4->add_attribute("theta", std::make_shared<DoubleAttribute>(std::rand()/double(RAND_MAX)*M_PI));
95 p4->add_attribute("phi", std::make_shared<DoubleAttribute>(std::rand()/double(RAND_MAX)*M_PI*2));
96 //
97 // create v3
98 GenVertexPtr v3=std::make_shared<GenVertex>();
99 evt.add_vertex( v3 );
100 v3->add_particle_in( p3 );
101 v3->add_particle_in( p4 );
102 GenParticlePtr p6=std::make_shared<GenParticle>( FourVector(-3.813,0.113,-1.833,4.233 ),22, 1 );
103 evt.add_particle( p6 );
104 p6->add_attribute("flow1", std::make_shared<IntAttribute>(231));
105 p6->add_attribute("theta", std::make_shared<DoubleAttribute>(std::rand()/double(RAND_MAX)*M_PI));
106 p6->add_attribute("phi", std::make_shared<DoubleAttribute>(std::rand()/double(RAND_MAX)*M_PI*2));
107 v3->add_particle_out( p6 );
108 GenParticlePtr p5=std::make_shared<GenParticle>( FourVector(1.517,-20.68,-20.605,85.925),-24, 3 );
109 v3->add_particle_out( p5 );
110 p5->add_attribute("flow1", std::make_shared<IntAttribute>(243));
111 p5->add_attribute("theta", std::make_shared<DoubleAttribute>(std::rand()/double(RAND_MAX)*M_PI));
112 p5->add_attribute("phi", std::make_shared<DoubleAttribute>(std::rand()/double(RAND_MAX)*M_PI*2));
113
114 //
115 // create v4
116 GenVertexPtr v4=std::make_shared<GenVertex>(FourVector(0.12,-0.3,0.05,0.004));
117 evt.add_vertex( v4 );
118 v4->add_particle_in( p5 );
119 GenParticlePtr p7(new GenParticle( FourVector(-2.445,28.816,6.082,29.552), 1,1 ));
120 v4->add_particle_out( p7 );
121 GenParticlePtr p8(new GenParticle( FourVector(3.962,-49.498,-26.687,56.373), -2,1 ));
122 v4->add_particle_out( p8 );
123
124
125 GenParticlePtr ploop=std::make_shared<GenParticle>( FourVector(0.0,0.0,0.0,0.0 ),21, 3 );
126 v3->add_particle_out( ploop );
127 v2->add_particle_in( ploop );
128
129 //
130 // tell the event which vertex is the signal process vertex
131 //evt.set_signal_process_vertex( v3 );
132 evt.add_attribute("signal_process_vertex", std::make_shared<IntAttribute>(v3->id()));
133 // the event is complete, we now print it out
134 Print::content(evt);
135 //we now print it out in old format
136 Print::listing(evt,8);
137 // print each particle so we can see the polarization
138 for ( ConstGenParticlePtr ip: evt.particles()) {
139 Print::line(ip,true);
140 }
141
142 // write event
143 xout1.write_event(evt);
144 // write event in old format
145 xout2.write_event(evt);
146
147 // now clean-up by deleteing all objects from memory
148 //
149 // deleting the event deletes all contained vertices, and all particles
150 // contained in those vertices
151 evt.clear();
152 xout1.close();
153 xout2.close();
154
155 int Nxin1=0;
156 ReaderAscii xin1("testLoops1.out");
157 if(xin1.failed()) {
158 xin1.close();
159 return 102;
160 }
161 while( !xin1.failed() )
162 {
163 xin1.read_event(evt);
164 if( xin1.failed() ) {
165 printf("End of file reached. Exit.\n");
166 break;
167 }
168 evt.clear();
169 Nxin1++;
170 }
171 xin1.close();
172
173 int Nxin2=0;
174 ReaderAsciiHepMC2 xin2("testLoops2.out");
175 if(xin2.failed()) {
176 xin2.close();
177 return 103;
178 }
179 while( !xin2.failed() )
180 {
181 xin2.read_event(evt);
182 if( xin2.failed() ) {
183 printf("End of file reached. Exit.\n");
184 break;
185 }
186 evt.clear();
187 Nxin2++;
188 }
189 xin2.close();
190 int ret = 10*std::abs(Nxin1-1)+std::abs(Nxin2-1);
191 return ret;
192 }
193