HepMC3 event record library
testPolarization.cc
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"
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 
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 }
GenEvent I/O serialization for structured text files.
Definition of class GenParticle.
#define M_PI
Definition of PI. Needed on some platforms.
Definition: FourVector.h:15
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
Stores event-related information.
Definition: GenEvent.h:41
Generic 4-vector.
Definition: FourVector.h:35
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
static void content(std::ostream &os, const GenEvent &event)
Print content of all GenEvent containers.
Definition: Print.cc:17