C++ Interface to Tauola
Particle.h
1#ifndef __PARTICLE_DEF__
2#define __PARTICLE_DEF__
3/**
4 Single particle class with basic methods such as boost, rotation,
5 and angle calculation.
6*/
7
8#include <cstdio>
9#include <iostream>
10#include <math.h>
11
12namespace TauSpinner {
13
14const double ELECTRON_MASS = 0.0005111;
15const double MUON_MASS = 0.105659;
16const double TAU_MASS = 1.777;
17
18class Particle
19{
20public:
21/*
22 Constructors
23*/
24 Particle():_px(0),_py(0),_pz(0),_e(0),_pdgid(0) {}
25 Particle(double px, double py, double pz, double e, int id):_px(px),_py(py),_pz(pz),_e(e),_pdgid(id) {}
26
27public:
28/*
29 Accessors
30*/
31 double px() const { return _px; }
32 double py() const { return _py; }
33 double pz() const { return _pz; }
34 double e () const { return _e; }
35 int pdgid() const { return _pdgid; }
36
37 // Invariant mass. If mass is negative then -sqrt(-mass) is returned
38 double recalculated_mass() const
39 {
40 double p2 = _px*_px + _py*_py + _pz*_pz;
41 double e2 = _e*_e;
42 double m2 = e2 - p2;
43
44 if ( m2 < 0.0 ) return -sqrt( -m2 );
45
46 return sqrt( m2 );
47 }
48
49 void setPx(double px) { _px = px; }
50 void setPy(double py) { _py = py; }
51 void setPz(double pz) { _pz = pz; }
52 void setE (double e ) { _e = e; }
53 void print()
54 {
55 if(_pdgid) printf("%4d: %15.8e %15.8e %15.8e %15.8e | %15.8e\n",
56 pdgid(),px(), py(), pz(), e(), recalculated_mass());
57 else printf(" SUM: %15.8e %15.8e %15.8e %15.8e | %15.8e\n",
58 px(), py(), pz(), e(), recalculated_mass());
59 }
60public:
61/*
62 These methods work on above accessors - they don't have to be changed
63 even if above methods change.
64*/
65 double getAnglePhi();
66 double getAngleTheta();
67 void rotateXZ(double theta);
68 void rotateXY(double theta);
69 void boostAlongZ(double pz, double e);
70 void boostToRestFrame(Particle &p);
71 void boostFromRestFrame(Particle &p);
72private:
73 double _px,_py,_pz,_e;
74 int _pdgid;
75};
76
77inline double Particle::getAnglePhi()
78{
79 // conventions as in ANGFI(X,Y) of tauola.f and PHOAN1 of photos.f
80 // but now 'phi' in name define that it is rotation in px py
81
82 double buf = 0.;
83
84 if(fabs(py())<fabs(px()))
85 {
86 buf = atan( fabs(py()/px()) );
87 if(px()<0.) buf = M_PI-buf;
88 }
89 else buf = acos( px()/sqrt(px()*px()+py()*py()) );
90
91 if(py()<0.) buf = 2.*M_PI-buf;
92
93 return buf;
94}
95
96inline double Particle::getAngleTheta()
97{
98 // conventions as in ANGXY(X,Y) of tauola.f or PHOAN2 of photos.f
99 // but now 'theta' in name define that it is rotation in pz px
100 // note that first argument of PHOAN2 was usually z
101
102 double buf = 0.;
103
104 if(fabs(px())<fabs(pz()))
105 {
106 buf = atan( fabs(px()/pz()) );
107 if(pz()<0.) buf = M_PI-buf;
108 }
109 else buf = acos( pz()/sqrt(pz()*pz()+px()*px()) );
110
111 return buf;
112}
113
114inline void Particle::rotateXZ(double theta)
115{
116 // as PHORO2 of photos.f
117 double pX=px();
118 double pZ=pz();
119 setPx( cos(theta)*pX + sin(theta)*pZ);
120 setPz(-sin(theta)*pX + cos(theta)*pZ);
121}
122
123inline void Particle::rotateXY(double theta)
124{
125 // as PHORO3 of photos.f
126 double pX=px();
127 double pY=py();
128 setPx( cos(theta)*pX - sin(theta)*pY);
129 setPy( sin(theta)*pX + cos(theta)*pY);
130}
131
132inline void Particle::boostAlongZ(double p_pz, double p_e)
133{
134 // as PHOBO3 of photos.f
135 double m=sqrt(p_e*p_e-p_pz*p_pz);
136 double buf_pz=pz();
137 double buf_e=e();
138
139 setPz((p_e *buf_pz + p_pz*buf_e)/m);
140 setE ((p_pz*buf_pz + p_e *buf_e)/m);
141}
142
143inline void Particle::boostToRestFrame(Particle &p)
144{
145 double p_len = sqrt(p.px()*p.px()+p.py()*p.py()+p.pz()*p.pz());
146 double phi = p.getAnglePhi();
147 p.rotateXY(-phi );
148 double theta = p.getAngleTheta();
149 p.rotateXY( phi );
150
151 //Now rotate coordinates to get boost in Z direction.
152 rotateXY(-phi );
153 rotateXZ(-theta);
154 boostAlongZ(-1*p_len,p.e());
155 rotateXZ( theta);
156 rotateXY( phi );
157}
158
159inline void Particle::boostFromRestFrame(Particle &p)
160{
161 double p_len = sqrt(p.px()*p.px()+p.py()*p.py()+p.pz()*p.pz());
162 double phi = p.getAnglePhi();
163 p.rotateXY(-phi );
164 double theta = p.getAngleTheta();
165 p.rotateXY( phi );
166
167 //Now rotate coordinates to get boost in Z direction.
168 rotateXY(-phi );
169 rotateXZ(-theta);
170 boostAlongZ(p_len,p.e());
171 rotateXZ( theta);
172 rotateXY( phi );
173}
174
175} // namespace TauSpinner
176#endif