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