HepMC3 event record library
testPolarization.cc
1//////////////////////////////////////////////////////////////////////////
2// testPolarization.cc
3//
4// garren@fnal.gov, Oct. 2010
5// andrii.verbytskyi@mpp.mpg.org, Nov. 2018 translated into HepMC3
6//////////////////////////////////////////////////////////////////////////
7
8#include <iostream>
9#include <fstream>
10#include <vector>
11
12#include "HepMC3/Attribute.h"
13#include "HepMC3/GenEvent.h"
14#include "HepMC3/GenVertex.h"
15#include "HepMC3/GenParticle.h"
16#include "HepMC3/WriterAscii.h"
18#include "HepMC3/Print.h"
19#ifndef M_PI
20#define M_PI 3.14159265358979323846264338327950288
21#endif
22#include "HepMC3TestUtils.h"
23using namespace HepMC3;
24int 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
40
41 // declare several WriterAscii instances for comparison
42 WriterAscii xout1("testPolarization1.dat");
43 WriterAscii xout2("testPolarization2.dat");
44 // output in old format
45 WriterAsciiHepMC2 xout4( "testPolarization4.out" );
46 WriterAscii xout5( "testPolarization5.out" );
47
48 // build the graph, which will look like
49 // p7 #
50 // p1 / #
51 // \v1__p3 p5---v4 #
52 // \_v3_/ \ #
53 // / \ p8 #
54 // v2__p4 \ #
55 // / p6 #
56 // p2 #
57 //
58 // define a flow pattern as p1 -> p3 -> p6
59 // and p2 -> p4 -> p5
60 //
61
62 // First create the event container, with Signal Process 20, event number 1
63 //
64 GenEvent evt(Units::GEV,Units::MM);
65 evt.set_event_number(1);
66 evt.add_attribute("signal_process_id", std::make_shared<IntAttribute>(20));
67 // create vertex 1
68 GenVertexPtr v1 = std::make_shared<GenVertex>();
69 evt.add_vertex( v1 );
70 GenParticlePtr p1 = std::make_shared<GenParticle>( FourVector(0, 0, 7000, 7000), 2212, 3 );
71 evt.add_particle( p1 );
72 p1->add_attribute("flow1", std::make_shared<IntAttribute>(231));
73 p1->add_attribute("flow1", std::make_shared<IntAttribute>(231));
74 p1->add_attribute("theta", std::make_shared<DoubleAttribute>(std::rand()/double(RAND_MAX)*M_PI));
75 p1->add_attribute("phi", std::make_shared<DoubleAttribute>(std::rand()/double(RAND_MAX)*M_PI*2));
76
77 GenVertexPtr v2 = std::make_shared<GenVertex>();
78 evt.add_vertex( v2 );
79 GenParticlePtr p2 = std::make_shared<GenParticle>( FourVector(0, 0, -7000, 7000), 2212, 3 );
80 evt.add_particle( p2 );
81 p2->add_attribute("flow1", std::make_shared<IntAttribute>(243));
82 p2->add_attribute("theta", std::make_shared<DoubleAttribute>(std::rand()/double(RAND_MAX)*M_PI));
83 p2->add_attribute("phi", std::make_shared<DoubleAttribute>(std::rand()/double(RAND_MAX)*M_PI*2));
84 v2->add_particle_in( p2 );
85 //
86 // create the outgoing particles of v1 and v2
87 GenParticlePtr p3 = std::make_shared<GenParticle>( FourVector(.750, -1.569, 32.191, 32.238),1, 3 );
88 evt.add_particle( p3 );
89 p3->add_attribute("flow1", std::make_shared<IntAttribute>(231));
90 p3->add_attribute("theta", std::make_shared<DoubleAttribute>(std::rand()/double(RAND_MAX)*M_PI));
91 p3->add_attribute("phi", std::make_shared<DoubleAttribute>(std::rand()/double(RAND_MAX)*M_PI*2));
92 v1->add_particle_out( p3 );
93 GenParticlePtr p4 = std::make_shared<GenParticle>( FourVector(-3.047, -19., -54.629, 57.920),-2, 3 );
94 evt.add_particle( p4 );
95 p4->add_attribute("flow1", std::make_shared<IntAttribute>(243));
96 p4->add_attribute("theta", std::make_shared<DoubleAttribute>(std::rand()/double(RAND_MAX)*M_PI));
97 p4->add_attribute("phi", std::make_shared<DoubleAttribute>(std::rand()/double(RAND_MAX)*M_PI*2));
98 v2->add_particle_out( p4 );
99 //
100 // create v3
101 GenVertexPtr v3 = std::make_shared<GenVertex>();
102 evt.add_vertex( v3 );
103 v3->add_particle_in( p3 );
104 v3->add_particle_in( p4 );
105 GenParticlePtr p6 = std::make_shared<GenParticle>( FourVector(-3.813, 0.113, -1.833, 4.233 ), 22, 1 );
106 evt.add_particle( p6 );
107 p6->add_attribute("flow1", std::make_shared<IntAttribute>(231));
108 p6->add_attribute("theta", std::make_shared<DoubleAttribute>(std::rand()/double(RAND_MAX)*M_PI));
109 p6->add_attribute("phi", std::make_shared<DoubleAttribute>(std::rand()/double(RAND_MAX)*M_PI*2));
110 v3->add_particle_out( p6 );
111 GenParticlePtr p5 = std::make_shared<GenParticle>( FourVector(1.517, -20.68, -20.605, 85.925), -24, 3 );
112 evt.add_particle( p5 );
113 p5->add_attribute("flow1", std::make_shared<IntAttribute>(243));
114 p5->add_attribute("theta", std::make_shared<DoubleAttribute>(std::rand()/double(RAND_MAX)*M_PI));
115 p5->add_attribute("phi", std::make_shared<DoubleAttribute>(std::rand()/double(RAND_MAX)*M_PI*2));
116 v3->add_particle_out( p5 );
117 //
118 // create v4
119 GenVertexPtr v4 = std::make_shared<GenVertex>(FourVector(0.12, -0.3, 0.05, 0.004));
120 evt.add_vertex( v4 );
121 v4->add_particle_in( p5 );
122 GenParticlePtr p7(new GenParticle( FourVector(-2.445, 28.816, 6.082, 29.552), 1,1 ));
123 evt.add_particle( p7 );
124 v4->add_particle_out( p7 );
125 GenParticlePtr p8(new GenParticle( FourVector(3.962, -49.498, -26.687, 56.373), -2,1 ));
126 evt.add_particle( p8 );
127 v4->add_particle_out( p8 );
128 //
129 // tell the event which vertex is the signal process vertex
130 //evt.set_signal_process_vertex( v3 );
131 evt.add_attribute("signal_process_vertex", std::make_shared<IntAttribute>(v3->id()));
132 // the event is complete, we now print it out
133 Print::content(evt);
134 //we now print it out in old format
135 Print::listing(evt,8);
136 // print each particle so we can see the polarization
137 for ( GenParticlePtr ip: evt.particles()) {
138 Print::line(ip,true);
139 }
140
141 // write event
142 xout1.write_event(evt);
143 // write event in old format
144 xout4.write_event(evt);
145 // make a copy and write it
146 xout5.write_event(GenEvent(evt));
147
148 // try changing polarization
149 p2->add_attribute("theta", std::make_shared<DoubleAttribute>(std::rand()/double(RAND_MAX)*M_PI));
150 p2->add_attribute("phi", std::make_shared<DoubleAttribute>(std::rand()/double(RAND_MAX)*M_PI));
151 xout2.write_event(evt);
152
153 xout1.close();
154 xout2.close();
155 xout4.close();
156 xout5.close();
157
158 // now clean-up by deleteing all objects from memory
159 //
160 // deleting the event deletes all contained vertices, and all particles
161 // contained in those vertices
162 evt.clear();
163
164 bool passed = ((COMPARE_ASCII_FILES("testPolarization1.dat","testPolarization5.out") == 0)&&(COMPARE_ASCII_FILES("testPolarization1.dat","testPolarization2.dat") != 0));
165 if (!passed) return 1;
166 return 0;
167}
Definition of class Attribute, class IntAttribute and class StringAttribute.
Definition of class GenEvent.
Definition of class GenParticle.
Definition of class GenVertex.
Definition of static class Print.
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
GenEvent I/O serialization for structured text files.
GenEvent I/O serialization for structured text files.
Definition: WriterAscii.h:25
HepMC3 main namespace.