1 //////////////////////////////////////////////////////////////////////////
2 // testPolarization.cc
3 //
4 // garren@fnal.gov, Oct. 2010
5 // andrii.verbytskyi@mpp.mpg.gov, Nov. 2018
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"
17 #include "HepMC3/WriterAsciiHepMC2.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
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 }
168