1 | // HelicityBasics.h is a part of the PYTHIA event generator. |
---|
2 | // Copyright (C) 2012 Philip Ilten, Torbjorn Sjostrand. |
---|
3 | // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details. |
---|
4 | // Please respect the MCnet Guidelines, see GUIDELINES for details. |
---|
5 | |
---|
6 | // Header file for a number of helper classes used in tau decays. |
---|
7 | |
---|
8 | #ifndef Pythia8_HelicityBasics_H |
---|
9 | #define Pythia8_HelicityBasics_H |
---|
10 | |
---|
11 | #include "Basics.h" |
---|
12 | #include "Event.h" |
---|
13 | #include "PythiaComplex.h" |
---|
14 | #include "PythiaStdlib.h" |
---|
15 | |
---|
16 | namespace Pythia8 { |
---|
17 | |
---|
18 | //========================================================================== |
---|
19 | |
---|
20 | // The Wave4 class provides a class for complex four-vector wave functions. |
---|
21 | // The Wave4 class can be multiplied with the GammaMatrix class to allow |
---|
22 | // for the writing of helicity matrix elements. |
---|
23 | |
---|
24 | class Wave4 { |
---|
25 | |
---|
26 | public: |
---|
27 | |
---|
28 | // Constructors and destructor. |
---|
29 | Wave4() {}; |
---|
30 | Wave4(complex v0, complex v1, complex v2, complex v3) {val[0] = v0; |
---|
31 | val[1] = v1; val[2] = v2; val[3] = v3;} |
---|
32 | Wave4(Vec4 v) {val[0] = v.e(); val[1] = v.px(); val[2] = v.py(); |
---|
33 | val[3] = v.pz();} |
---|
34 | ~Wave4() {}; |
---|
35 | |
---|
36 | // Access an element of the wave vector. |
---|
37 | complex& operator() (int i) {return val[i];} |
---|
38 | |
---|
39 | // Wave4 + Wave4. |
---|
40 | Wave4 operator+(Wave4 w) {return Wave4( val[0] + w.val[0], |
---|
41 | val[1] + w.val[1], val[2] + w.val[2], val[3] + w.val[3]);} |
---|
42 | |
---|
43 | // Wave4 - Wave4. |
---|
44 | Wave4 operator-(Wave4 w) {return Wave4( val[0] - w.val[0], |
---|
45 | val[1] - w.val[1], val[2] - w.val[2], val[3] - w.val[3]);} |
---|
46 | |
---|
47 | // - Wave4. |
---|
48 | Wave4 operator-() {return Wave4(-val[0], -val[1], -val[2], -val[3]);} |
---|
49 | |
---|
50 | // Wave4 * Wave4. |
---|
51 | complex operator*(Wave4 w) {return val[0] * w.val[0] |
---|
52 | + val[1] * w.val[1] + val[2] * w.val[2] + val[3] * w.val[3];} |
---|
53 | |
---|
54 | // Wave4 * complex. |
---|
55 | Wave4 operator*(complex s) {return Wave4(val[0] * s, val[1] * s, |
---|
56 | val[2] * s, val[3] * s);} |
---|
57 | |
---|
58 | // complex * Wave4. |
---|
59 | friend Wave4 operator*(complex s, const Wave4& w); |
---|
60 | |
---|
61 | // Wave4 * double. |
---|
62 | Wave4 operator*(double s) {return Wave4(val[0] * s, val[1] * s, |
---|
63 | val[2] * s, val[3] * s);} |
---|
64 | |
---|
65 | // double * Wave4. |
---|
66 | friend Wave4 operator*(double s, const Wave4& w); |
---|
67 | |
---|
68 | // Wave4 / complex. |
---|
69 | Wave4 operator/(complex s) {return Wave4(val[0] / s, val[1] / s, |
---|
70 | val[2] / s, val[3] / s);} |
---|
71 | |
---|
72 | // Wave4 / double. |
---|
73 | Wave4 operator/(double s) {return Wave4(val[0] / s, val[1] / s, |
---|
74 | val[2]/s, val[3]/s);} |
---|
75 | |
---|
76 | // Complex conjugate. |
---|
77 | friend Wave4 conj(Wave4 w); |
---|
78 | |
---|
79 | // Permutation operator. |
---|
80 | friend Wave4 epsilon(Wave4 w1, Wave4 w2, Wave4 w3); |
---|
81 | |
---|
82 | // Invariant squared mass for REAL Wave4 (to save time). |
---|
83 | friend double m2(Wave4 w); |
---|
84 | friend double m2(Wave4 w1, Wave4 w2); |
---|
85 | |
---|
86 | // Wave4 * GammaMatrix multiplication is defined in the GammaMatrix class. |
---|
87 | |
---|
88 | // Print a Wave4 vector. |
---|
89 | friend ostream& operator<<(ostream& output, Wave4 w); |
---|
90 | |
---|
91 | protected: |
---|
92 | |
---|
93 | complex val[4]; |
---|
94 | |
---|
95 | }; |
---|
96 | |
---|
97 | //-------------------------------------------------------------------------- |
---|
98 | |
---|
99 | // Namespace function declarations; friends of Wave4 class. |
---|
100 | Wave4 operator*(complex s, const Wave4& w); |
---|
101 | Wave4 conj(Wave4 w); |
---|
102 | Wave4 epsilon(Wave4 w1, Wave4 w2, Wave4 w3); |
---|
103 | double m2(Wave4 w); |
---|
104 | double m2(Wave4 w1, Wave4 w2); |
---|
105 | ostream& operator<< (ostream& os, Wave4 w); |
---|
106 | |
---|
107 | //========================================================================== |
---|
108 | |
---|
109 | // The GammaMatrix class is a special sparse matrix class used to write |
---|
110 | // helicity matrix elements in conjuction with the Wave4 class. Note that |
---|
111 | // only left to right multplication of Wave4 vectors with the GammaMatrix |
---|
112 | // class is allowed. Additionally, subtracting a scalar from a GammaMatrix |
---|
113 | // (or subtracting a GammaMatrix from a scalar) subtracts the scalar from |
---|
114 | //each non-zero element of the GammaMatrix. This is designed specifically |
---|
115 | // with the (1 - gamma^5) structure of matrix elements in mind. |
---|
116 | |
---|
117 | class GammaMatrix { |
---|
118 | |
---|
119 | public: |
---|
120 | |
---|
121 | // Constructors and destructor. |
---|
122 | GammaMatrix() {}; |
---|
123 | GammaMatrix(int mu); |
---|
124 | ~GammaMatrix() {}; |
---|
125 | |
---|
126 | // Access an element of the matrix. |
---|
127 | complex& operator() (int I, int J) {if (index[J] == I) return val[J]; |
---|
128 | else return COMPLEXZERO; } |
---|
129 | |
---|
130 | // Wave4 * GammaMatrix. |
---|
131 | friend Wave4 operator*(Wave4 w, GammaMatrix g); |
---|
132 | |
---|
133 | // GammaMatrix * Scalar. |
---|
134 | GammaMatrix operator*(complex s) {val[0] = s*val[0]; val[1] = s*val[1]; |
---|
135 | val[2] = s*val[2]; val[3] = s*val[3]; return *this;} |
---|
136 | |
---|
137 | // Scalar * GammaMatrix. |
---|
138 | friend GammaMatrix operator*(complex s, GammaMatrix g); |
---|
139 | |
---|
140 | // Gamma5 - I * Scalar. |
---|
141 | GammaMatrix operator-(complex s) {val[0] = val[0] - s; val[1] = val[1] - s; |
---|
142 | val[2] = val[2] - s; val[3] = val[3] - s; return *this;} |
---|
143 | |
---|
144 | // I * Scalar - Gamma5. |
---|
145 | friend GammaMatrix operator-(complex s, GammaMatrix g); |
---|
146 | |
---|
147 | // Gamma5 + I * Scalar |
---|
148 | GammaMatrix operator+(complex s) {val[0] = val[0] + s; val[1] = val[1] + s; |
---|
149 | val[2] = val[2] + s; val[3] = val[3] + s; return *this;} |
---|
150 | |
---|
151 | // I * Scalar + Gamma5 |
---|
152 | friend GammaMatrix operator+(complex s, GammaMatrix g); |
---|
153 | |
---|
154 | // << GammaMatrix. |
---|
155 | friend ostream& operator<< (ostream& os, GammaMatrix g); |
---|
156 | |
---|
157 | protected: |
---|
158 | |
---|
159 | complex val[4]; |
---|
160 | int index[4]; |
---|
161 | |
---|
162 | // Need to define complex 0 as a variable for operator() to work. |
---|
163 | complex COMPLEXZERO; |
---|
164 | |
---|
165 | }; |
---|
166 | |
---|
167 | //-------------------------------------------------------------------------- |
---|
168 | |
---|
169 | // Namespace function declarations; friends of GammaMatrix class. |
---|
170 | Wave4 operator*(Wave4 w, GammaMatrix g); |
---|
171 | GammaMatrix operator*(complex s, GammaMatrix g); |
---|
172 | GammaMatrix operator-(complex s, GammaMatrix g); |
---|
173 | GammaMatrix operator+(complex s, GammaMatrix g); |
---|
174 | ostream& operator<< (ostream& os, GammaMatrix g); |
---|
175 | |
---|
176 | //========================================================================== |
---|
177 | |
---|
178 | // Helicity particle class containing helicity information, derived from |
---|
179 | // particle base class. |
---|
180 | |
---|
181 | class HelicityParticle : public Particle { |
---|
182 | |
---|
183 | public: |
---|
184 | |
---|
185 | // Constructors. |
---|
186 | HelicityParticle() : Particle() { direction = 1;} |
---|
187 | HelicityParticle(int idIn, int statusIn = 0, int mother1In = 0, |
---|
188 | int mother2In = 0, int daughter1In = 0, int daughter2In = 0, |
---|
189 | int colIn = 0, int acolIn = 0, double pxIn = 0., |
---|
190 | double pyIn = 0., double pzIn = 0., double eIn = 0., |
---|
191 | double mIn = 0., double scaleIn = 0., ParticleData* ptr = 0) |
---|
192 | : Particle(idIn, statusIn, mother1In, mother2In, daughter1In, daughter2In, |
---|
193 | colIn, acolIn, pxIn, pyIn, pzIn, eIn, mIn, scaleIn) { |
---|
194 | if (ptr) { setPDTPtr(ptr); setPDEPtr(); } |
---|
195 | rho = vector< vector<complex> >(spinStates(), |
---|
196 | vector<complex>(spinStates(), 0)); |
---|
197 | D = vector< vector<complex> >(spinStates(), |
---|
198 | vector<complex>(spinStates(), 0)); |
---|
199 | for (int i = 0; i < spinStates(); i++) { rho[i][i] = 0.5; D[i][i] = 1.;} |
---|
200 | direction = 1; } |
---|
201 | HelicityParticle(int idIn, int statusIn, int mother1In, int mother2In, |
---|
202 | int daughter1In, int daughter2In, int colIn, int acolIn, Vec4 pIn, |
---|
203 | double mIn = 0., double scaleIn = 0., ParticleData* ptr = 0) |
---|
204 | : Particle(idIn, statusIn, mother1In, mother2In, daughter1In, daughter2In, |
---|
205 | colIn, acolIn, pIn, mIn, scaleIn) { |
---|
206 | if (ptr) { setPDTPtr(ptr); setPDEPtr();} |
---|
207 | rho = vector< vector<complex> >(spinStates(), |
---|
208 | vector<complex>(spinStates(), 0)); |
---|
209 | D = vector< vector<complex> >(spinStates(), |
---|
210 | vector<complex>(spinStates(), 0)); |
---|
211 | for (int i = 0; i < spinStates(); i++) { rho[i][i] = 0.5; D[i][i] = 1;} |
---|
212 | direction = 1; } |
---|
213 | HelicityParticle(const Particle& ptIn, ParticleData* ptr = 0) |
---|
214 | : Particle(ptIn) { |
---|
215 | if (ptr) { setPDTPtr(ptr); setPDEPtr();} |
---|
216 | rho = vector< vector<complex> >(spinStates(), |
---|
217 | vector<complex>(spinStates(), 0)); |
---|
218 | D = vector< vector<complex> >(spinStates(), |
---|
219 | vector<complex>(spinStates(), 0)); |
---|
220 | for (int i = 0; i < spinStates(); i++) { rho[i][i] = 0.5; D[i][i] = 1;} |
---|
221 | direction = 1; } |
---|
222 | |
---|
223 | // Methods. |
---|
224 | Wave4 wave(int h); |
---|
225 | Wave4 waveBar(int h); |
---|
226 | void normalize(vector< vector<complex> >& m); |
---|
227 | int spinStates(); |
---|
228 | |
---|
229 | // Event record position. |
---|
230 | int idx; |
---|
231 | |
---|
232 | // Flag for whether particle is incoming (-1) or outgoing (1). |
---|
233 | int direction; |
---|
234 | |
---|
235 | // Helicity density matrix. |
---|
236 | vector< vector<complex> > rho; |
---|
237 | |
---|
238 | // Decay matrix. |
---|
239 | vector< vector<complex> > D; |
---|
240 | |
---|
241 | private: |
---|
242 | |
---|
243 | // Constants: could only be changed in the code itself. |
---|
244 | static const double TOLERANCE; |
---|
245 | |
---|
246 | }; |
---|
247 | |
---|
248 | //========================================================================== |
---|
249 | |
---|
250 | } // end namespace Pythia8 |
---|
251 | |
---|
252 | #endif // end Pythia8_HelicityBasics_H |
---|