HepMC3 event record library
testLoops.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 #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"
16 #include "HepMC3/ReaderAscii.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;
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 }
GenEvent I/O serialization for structured text files.
HepMC3 main namespace.
GenEvent I/O parsing for structured text files.
Definition: ReaderAscii.h:29
Definition of class GenParticle.
#define M_PI
Definition of PI. Needed on some platforms.
Definition: FourVector.h:16
Definition of class GenVertex.
Definition of class WriterAscii.
Stores particle-related information.
Definition: GenParticle.h:31
static void listing(std::ostream &os, const GenEvent &event, unsigned short precision=2)
Print event in listing (HepMC2) format.
Definition: Print.cc:50
Parser for HepMC2 I/O files.
Definition of class ReaderAsciiHepMC2.
Stores event-related information.
Definition: GenEvent.h:41
Generic 4-vector.
Definition: FourVector.h:36
Definition of class ReaderAscii.
static void line(std::ostream &os, const GenEvent &event, bool attributes=false)
Print one-line info.
Definition: Print.cc:202
Definition of class WriterAsciiHepMC2.
int main(int argc, char **argv)
Definition of class GenEvent.
Definition of class Attribute, class IntAttribute and class StringAttribute.
GenEvent I/O serialization for structured text files.
Definition: WriterAscii.h:25
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:323
static void content(std::ostream &os, const GenEvent &event)
Print content of all GenEvent containers.
Definition: Print.cc:17