1 | // HelicityBasics.cc 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 | // Function definitions (not found in the header) for helper classes |
---|
7 | // used in tau decays. |
---|
8 | |
---|
9 | #include "HelicityBasics.h" |
---|
10 | |
---|
11 | namespace Pythia8 { |
---|
12 | |
---|
13 | //========================================================================== |
---|
14 | |
---|
15 | // Wave4 class. |
---|
16 | // Friend methods to it. |
---|
17 | |
---|
18 | //-------------------------------------------------------------------------- |
---|
19 | |
---|
20 | // complex * Wave4. |
---|
21 | |
---|
22 | Wave4 operator*(complex s, const Wave4& w) { |
---|
23 | |
---|
24 | return Wave4( s * w.val[0], s * w.val[1], s * w.val[2], s * w.val[3]); |
---|
25 | |
---|
26 | } |
---|
27 | |
---|
28 | //-------------------------------------------------------------------------- |
---|
29 | |
---|
30 | // double * Wave4. |
---|
31 | |
---|
32 | Wave4 operator*(double s, const Wave4& w) { |
---|
33 | |
---|
34 | return Wave4( s * w.val[0], s * w.val[1], s * w.val[2], s * w.val[3]); |
---|
35 | |
---|
36 | } |
---|
37 | |
---|
38 | //-------------------------------------------------------------------------- |
---|
39 | |
---|
40 | // Complex conjugate. |
---|
41 | |
---|
42 | Wave4 conj(Wave4 w) { |
---|
43 | |
---|
44 | w(0) = conj(w(0)); |
---|
45 | w(1) = conj(w(1)); |
---|
46 | w(2) = conj(w(2)); |
---|
47 | w(3) = conj(w(3)); |
---|
48 | return w; |
---|
49 | |
---|
50 | } |
---|
51 | |
---|
52 | //-------------------------------------------------------------------------- |
---|
53 | |
---|
54 | // Permutation operator. |
---|
55 | |
---|
56 | Wave4 epsilon(Wave4 w1, Wave4 w2, Wave4 w3) { |
---|
57 | |
---|
58 | Wave4 w4; |
---|
59 | w4(0) = -(w1(1) * w2(2) * w3(3)) + (w1(1) * w2(3) * w3(2)) |
---|
60 | + (w1(2) * w2(1) * w3(3)) - (w1(2) * w2(3) * w3(1)) |
---|
61 | - (w1(3) * w2(1) * w3(2)) + (w1(3) * w2(2) * w3(1)); |
---|
62 | w4(1) = -(w1(0) * w2(2) * w3(3)) + (w1(0) * w2(3) * w3(2)) |
---|
63 | + (w1(2) * w2(0) * w3(3)) - (w1(2) * w2(3) * w3(0)) |
---|
64 | - (w1(3) * w2(0) * w3(2)) + (w1(3) * w2(2) * w3(0)); |
---|
65 | w4(2) = (w1(0) * w2(1) * w3(3)) - (w1(0) * w2(3) * w3(1)) |
---|
66 | - (w1(1) * w2(0) * w3(3)) + (w1(1) * w2(3) * w3(0)) |
---|
67 | + (w1(3) * w2(0) * w3(1)) - (w1(3) * w2(1) * w3(0)); |
---|
68 | w4(3) = -(w1(0) * w2(1) * w3(2)) + (w1(0) * w2(2) * w3(1)) |
---|
69 | + (w1(1) * w2(0) * w3(2)) - (w1(1) * w2(2) * w3(0)) |
---|
70 | - (w1(2) * w2(0) * w3(1)) + (w1(2) * w2(1) * w3(0)); |
---|
71 | return w4; |
---|
72 | |
---|
73 | } |
---|
74 | |
---|
75 | //-------------------------------------------------------------------------- |
---|
76 | |
---|
77 | // Invariant squared mass for REAL Wave4 (to save time). |
---|
78 | |
---|
79 | double m2(Wave4 w) { |
---|
80 | |
---|
81 | return real(w(0)) * real(w(0)) - real(w(1)) * real(w(1)) |
---|
82 | - real(w(2)) * real(w(2)) - real(w(3)) * real(w(3)); |
---|
83 | |
---|
84 | } |
---|
85 | |
---|
86 | double m2(Wave4 w1, Wave4 w2) { |
---|
87 | |
---|
88 | return real(w1(0)) * real(w2(0)) - real(w1(1)) * real(w2(1)) |
---|
89 | - real(w1(2)) * real(w2(2)) - real(w1(3)) * real(w2(3)); |
---|
90 | |
---|
91 | } |
---|
92 | |
---|
93 | //-------------------------------------------------------------------------- |
---|
94 | |
---|
95 | // Print a Wave4 vector. |
---|
96 | |
---|
97 | ostream& operator<< (ostream& os, Wave4 w) { |
---|
98 | |
---|
99 | os << left << setprecision(2); |
---|
100 | for (int i = 0; i < 4; i++) os << setw(20) << w.val[i]; |
---|
101 | os << "\n"; |
---|
102 | return os; |
---|
103 | |
---|
104 | } |
---|
105 | |
---|
106 | //========================================================================== |
---|
107 | |
---|
108 | // Constructor for the GammaMatrix class. Gamma(1) through Gamma(3) give the |
---|
109 | // corresponding gamma matrices using the Weyl basis as outlined in the HELAS |
---|
110 | // paper. Gamma(4) gives the +--- metric, while Gamma(5) gives the gamma^5 |
---|
111 | // matrix. |
---|
112 | |
---|
113 | GammaMatrix::GammaMatrix(int mu) { |
---|
114 | |
---|
115 | COMPLEXZERO = complex( 0., 0.); |
---|
116 | |
---|
117 | if (mu == 0) { |
---|
118 | val[0] = 1.; val[1] = 1.; val[2] = 1.; val[3] = 1.; |
---|
119 | index[0] = 2; index[1] = 3; index[2] = 0; index[3] = 1; |
---|
120 | |
---|
121 | } else if (mu == 1) { |
---|
122 | val[0] = -1.; val[1] = -1.; val[2] = 1.; val[3] = 1.; |
---|
123 | index[0] = 3; index[1] = 2; index[2] = 1; index[3] = 0; |
---|
124 | |
---|
125 | } else if (mu == 2) { |
---|
126 | val[0] = complex(0.,-1.); val[1] = complex(0.,1.); |
---|
127 | val[2] = complex(0.,1.); val[3] = complex(0.,-1.); |
---|
128 | index[0] = 3; index[1] = 2; index[2] = 1; index[3] = 0; |
---|
129 | |
---|
130 | } else if (mu == 3) { |
---|
131 | val[0] = -1.; val[1] = 1.; val[2] = 1.; val[3] = -1.; |
---|
132 | index[0] = 2; index[1] = 3; index[2] = 0; index[3] = 1; |
---|
133 | |
---|
134 | } else if (mu == 4) { |
---|
135 | val[0] = 1.; val[1] = -1.; val[2] = -1.; val[3] = -1.; |
---|
136 | index[0] = 0; index[1] = 1; index[2] = 2; index[3] = 3; |
---|
137 | |
---|
138 | } else if (mu == 5) { |
---|
139 | val[0] = -1.; val[1] = -1.; val[2] = 1.; val[3] = 1.; |
---|
140 | index[0] = 0; index[1] = 1; index[2] = 2; index[3] = 3; |
---|
141 | } |
---|
142 | |
---|
143 | } |
---|
144 | |
---|
145 | //-------------------------------------------------------------------------- |
---|
146 | |
---|
147 | // Wave4 * GammaMatrix. |
---|
148 | |
---|
149 | Wave4 operator*(Wave4 w, GammaMatrix g) { |
---|
150 | |
---|
151 | complex w0 = w(g.index[0]); |
---|
152 | complex w1 = w(g.index[1]); |
---|
153 | complex w2 = w(g.index[2]); |
---|
154 | complex w3 = w(g.index[3]); |
---|
155 | w(0) = w0 * g.val[0]; |
---|
156 | w(1) = w1 * g.val[1]; |
---|
157 | w(2) = w2 * g.val[2]; |
---|
158 | w(3) = w3 * g.val[3]; |
---|
159 | return w; |
---|
160 | |
---|
161 | } |
---|
162 | |
---|
163 | //-------------------------------------------------------------------------- |
---|
164 | |
---|
165 | // Scalar * GammaMatrix. |
---|
166 | |
---|
167 | GammaMatrix operator*(complex s, GammaMatrix g) { |
---|
168 | |
---|
169 | g.val[0] = s * g.val[0]; |
---|
170 | g.val[1] = s*g.val[1]; |
---|
171 | g.val[2] = s * g.val[2]; |
---|
172 | g.val[3] = s*g.val[3]; |
---|
173 | return g; |
---|
174 | |
---|
175 | } |
---|
176 | |
---|
177 | //-------------------------------------------------------------------------- |
---|
178 | |
---|
179 | // I * Scalar - Gamma5. |
---|
180 | |
---|
181 | GammaMatrix operator-(complex s, GammaMatrix g) { |
---|
182 | |
---|
183 | g.val[0] = s - g.val[0]; |
---|
184 | g.val[1] = s - g.val[1]; |
---|
185 | g.val[2] = s - g.val[2]; |
---|
186 | g.val[3] = s - g.val[3]; |
---|
187 | return g; |
---|
188 | |
---|
189 | } |
---|
190 | |
---|
191 | //-------------------------------------------------------------------------- |
---|
192 | |
---|
193 | // I * Scalar + Gamma5. |
---|
194 | |
---|
195 | GammaMatrix operator+(complex s, GammaMatrix g) { |
---|
196 | |
---|
197 | g.val[0] = s + g.val[0]; |
---|
198 | g.val[1] = s + g.val[1]; |
---|
199 | g.val[2] = s + g.val[2]; |
---|
200 | g.val[3] = s + g.val[3]; |
---|
201 | return g; |
---|
202 | |
---|
203 | } |
---|
204 | |
---|
205 | //-------------------------------------------------------------------------- |
---|
206 | |
---|
207 | // Print GammaMatrix. |
---|
208 | |
---|
209 | ostream& operator<< (ostream& os, GammaMatrix g) { |
---|
210 | |
---|
211 | os << left << setprecision(2); |
---|
212 | for (int i = 0; i < 4; i++) { |
---|
213 | for (int j = 0; j < 4; j++) os << setw(20) << g(i,j); |
---|
214 | os << "\n"; |
---|
215 | } |
---|
216 | return os; |
---|
217 | |
---|
218 | } |
---|
219 | |
---|
220 | //========================================================================== |
---|
221 | |
---|
222 | // Weyl helicity wave functions for spin 1/2 fermions and spin 1 boson. |
---|
223 | |
---|
224 | // This is the basis as given by the HELAS collaboration on page 122 of |
---|
225 | // "HELas: HELicity Amplitude Subroutines for Feynman Diagram Evaluations" |
---|
226 | // by H. Murayama, I. Watanabe, K. Hagiwara. |
---|
227 | |
---|
228 | //-------------------------------------------------------------------------- |
---|
229 | |
---|
230 | // Constants: could be changed here if desired, but normally should not. |
---|
231 | |
---|
232 | // The spinors become ill-defined for p -> -pz and the polarization vectors |
---|
233 | // become ill-defined when pT -> 0. For these special cases limits are used |
---|
234 | // and are triggered when the difference of the approached quantity falls |
---|
235 | // below the variable TOLERANCE. |
---|
236 | const double HelicityParticle::TOLERANCE = 0.000001; |
---|
237 | |
---|
238 | //-------------------------------------------------------------------------- |
---|
239 | |
---|
240 | // Return wave vector for given helicity. |
---|
241 | |
---|
242 | Wave4 HelicityParticle::wave(int h) { |
---|
243 | |
---|
244 | // Create wave vector to return. |
---|
245 | Wave4 w; |
---|
246 | |
---|
247 | // Fermion (spin 1/2) spinor. |
---|
248 | if (spinType() == 2) { |
---|
249 | |
---|
250 | // Calculate helicity independent normalization. |
---|
251 | double P = pAbs(); |
---|
252 | double n = sqrtpos(2*P*(P+pz())); |
---|
253 | bool aligned = (abs(P+pz()) < TOLERANCE); |
---|
254 | |
---|
255 | // Calculate eigenspinor basis. |
---|
256 | vector< vector<complex> > xi(2, vector<complex>(2)); |
---|
257 | // Helicity -1 |
---|
258 | xi[0][0] = aligned ? -1 : complex(-px(),py())/n; |
---|
259 | xi[0][1] = aligned ? 0 : (P+pz())/n; |
---|
260 | // Helicity +1 |
---|
261 | xi[1][0] = aligned ? 0 : (P+pz())/n; |
---|
262 | xi[1][1] = aligned ? 1 : complex(px(),py())/n; |
---|
263 | |
---|
264 | // Calculate helicity dependent normalization. |
---|
265 | vector<double> omega(2); |
---|
266 | omega[0] = sqrtpos(e()-P); |
---|
267 | omega[1] = sqrtpos(e()+P); |
---|
268 | vector<double> hsign(2,1); |
---|
269 | hsign[0] = -1; |
---|
270 | |
---|
271 | // Create particle spinor. |
---|
272 | if (this->id() > 0) { |
---|
273 | w(0) = omega[!h] * xi[h][0]; |
---|
274 | w(1) = omega[!h] * xi[h][1]; |
---|
275 | w(2) = omega[h] * xi[h][0]; |
---|
276 | w(3) = omega[h] * xi[h][1]; |
---|
277 | |
---|
278 | // Create anti-particle spinor. |
---|
279 | } else { |
---|
280 | w(0) = hsign[!h] * omega[h] * xi[!h][0]; |
---|
281 | w(1) = hsign[!h] * omega[h] * xi[!h][1]; |
---|
282 | w(2) = hsign[h] * omega[!h] * xi[!h][0]; |
---|
283 | w(3) = hsign[h] * omega[!h] * xi[!h][1]; |
---|
284 | } |
---|
285 | |
---|
286 | // Boson (spin 1) polarization vector. |
---|
287 | } else if (spinType() == 3) { |
---|
288 | double P = pAbs(); |
---|
289 | double PT = pT(); |
---|
290 | |
---|
291 | // Create helicity +1 or -1 polarization vector. |
---|
292 | if (h >= 0 && h <= 1) { |
---|
293 | double hsign = h ? -1 : 1; |
---|
294 | if (PT > TOLERANCE) { |
---|
295 | w(0) = 0; |
---|
296 | w(1) = complex(hsign * px() * pz() / (P * PT), -py() / PT); |
---|
297 | w(2) = complex(hsign * py() * pz() / (P * PT), px() / PT); |
---|
298 | w(3) = complex(-hsign * PT / P, 0); |
---|
299 | } else { |
---|
300 | w(0) = 0; |
---|
301 | w(1) = hsign * pz(); |
---|
302 | w(2) = 0; |
---|
303 | w(3) = 0; |
---|
304 | } |
---|
305 | |
---|
306 | // Create helicity 0 polarization vector (ensure boson massive). |
---|
307 | } else if (h == 2 && m() > TOLERANCE) { |
---|
308 | w(0) = P / m(); |
---|
309 | w(1) = px() * e() / (m() * P); |
---|
310 | w(2) = py() * e() / (m() * P); |
---|
311 | w(3) = pz() * e() / (m() * P); |
---|
312 | } |
---|
313 | |
---|
314 | // Unknown wave function. |
---|
315 | } else { |
---|
316 | w(0) = 0; |
---|
317 | w(1) = 0; |
---|
318 | w(2) = 0; |
---|
319 | w(3) = 0; |
---|
320 | } |
---|
321 | |
---|
322 | // Done. |
---|
323 | return w; |
---|
324 | |
---|
325 | } |
---|
326 | |
---|
327 | //-------------------------------------------------------------------------- |
---|
328 | |
---|
329 | // Bar of a wave function. |
---|
330 | |
---|
331 | Wave4 HelicityParticle::waveBar(int h) { |
---|
332 | |
---|
333 | if (spinType() == 2) return conj(wave(h)) * GammaMatrix(0); |
---|
334 | else return conj(wave(h)); |
---|
335 | |
---|
336 | } |
---|
337 | |
---|
338 | //-------------------------------------------------------------------------- |
---|
339 | |
---|
340 | // Normalize the rho or D matrices. |
---|
341 | |
---|
342 | void HelicityParticle::normalize(vector< vector<complex> >& matrix) { |
---|
343 | |
---|
344 | complex trace = 0; |
---|
345 | for (unsigned int i = 0; i < matrix.size(); i++) trace += matrix[i][i]; |
---|
346 | for (unsigned int i = 0; i < matrix.size(); i++) { |
---|
347 | for (unsigned int j = 0; j < matrix.size(); j++) { |
---|
348 | if (trace != complex(0,0)) matrix[i][j] /= trace; |
---|
349 | else matrix[i][j] = 1 / static_cast<double>(matrix.size()); |
---|
350 | } |
---|
351 | } |
---|
352 | |
---|
353 | } |
---|
354 | |
---|
355 | //-------------------------------------------------------------------------- |
---|
356 | |
---|
357 | // Return the number of spin states. |
---|
358 | |
---|
359 | int HelicityParticle::spinStates() { |
---|
360 | |
---|
361 | int sT = spinType(); |
---|
362 | if (sT == 0) return 0; |
---|
363 | else if (sT != 2 && m() < TOLERANCE) return sT - 1; |
---|
364 | else return sT; |
---|
365 | |
---|
366 | } |
---|
367 | |
---|
368 | //========================================================================== |
---|
369 | |
---|
370 | } // end namespace Pythia8 |
---|