Line data Source code
1 : // HelicityBasics.cc is a part of the PYTHIA event generator.
2 : // Copyright (C) 2015 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 "Pythia8/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 0 : 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 0 : 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 0 : w(0) = conj(w(0));
45 0 : w(1) = conj(w(1));
46 0 : w(2) = conj(w(2));
47 0 : w(3) = conj(w(3));
48 0 : return w;
49 :
50 : }
51 :
52 : //--------------------------------------------------------------------------
53 :
54 : // Permutation operator.
55 :
56 : Wave4 epsilon(Wave4 w1, Wave4 w2, Wave4 w3) {
57 :
58 0 : Wave4 w4;
59 0 : w4(0) = -(w1(1) * w2(2) * w3(3)) + (w1(1) * w2(3) * w3(2))
60 0 : + (w1(2) * w2(1) * w3(3)) - (w1(2) * w2(3) * w3(1))
61 0 : - (w1(3) * w2(1) * w3(2)) + (w1(3) * w2(2) * w3(1));
62 0 : w4(1) = -(w1(0) * w2(2) * w3(3)) + (w1(0) * w2(3) * w3(2))
63 0 : + (w1(2) * w2(0) * w3(3)) - (w1(2) * w2(3) * w3(0))
64 0 : - (w1(3) * w2(0) * w3(2)) + (w1(3) * w2(2) * w3(0));
65 0 : w4(2) = (w1(0) * w2(1) * w3(3)) - (w1(0) * w2(3) * w3(1))
66 0 : - (w1(1) * w2(0) * w3(3)) + (w1(1) * w2(3) * w3(0))
67 0 : + (w1(3) * w2(0) * w3(1)) - (w1(3) * w2(1) * w3(0));
68 0 : w4(3) = -(w1(0) * w2(1) * w3(2)) + (w1(0) * w2(2) * w3(1))
69 0 : + (w1(1) * w2(0) * w3(2)) - (w1(1) * w2(2) * w3(0))
70 0 : - (w1(2) * w2(0) * w3(1)) + (w1(2) * w2(1) * w3(0));
71 : return w4;
72 :
73 0 : }
74 :
75 : //--------------------------------------------------------------------------
76 :
77 : // Invariant squared mass for REAL Wave4 (to save time).
78 :
79 : double m2(Wave4 w) {
80 :
81 0 : return real(w(0)) * real(w(0)) - real(w(1)) * real(w(1))
82 0 : - real(w(2)) * real(w(2)) - real(w(3)) * real(w(3));
83 :
84 : }
85 :
86 : double m2(Wave4 w1, Wave4 w2) {
87 :
88 0 : return real(w1(0)) * real(w2(0)) - real(w1(1)) * real(w2(1))
89 0 : - 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 0 : os << left << setprecision(2);
100 0 : for (int i = 0; i < 4; i++) os << setw(20) << w.val[i];
101 0 : os << "\n";
102 0 : 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 0 : GammaMatrix::GammaMatrix(int mu) {
114 :
115 0 : COMPLEXZERO = complex( 0., 0.);
116 :
117 0 : if (mu == 0) {
118 0 : val[0] = 1.; val[1] = 1.; val[2] = 1.; val[3] = 1.;
119 0 : index[0] = 2; index[1] = 3; index[2] = 0; index[3] = 1;
120 :
121 0 : } else if (mu == 1) {
122 0 : val[0] = -1.; val[1] = -1.; val[2] = 1.; val[3] = 1.;
123 0 : index[0] = 3; index[1] = 2; index[2] = 1; index[3] = 0;
124 :
125 0 : } else if (mu == 2) {
126 0 : val[0] = complex(0.,-1.); val[1] = complex(0.,1.);
127 0 : val[2] = complex(0.,1.); val[3] = complex(0.,-1.);
128 0 : index[0] = 3; index[1] = 2; index[2] = 1; index[3] = 0;
129 :
130 0 : } else if (mu == 3) {
131 0 : val[0] = -1.; val[1] = 1.; val[2] = 1.; val[3] = -1.;
132 0 : index[0] = 2; index[1] = 3; index[2] = 0; index[3] = 1;
133 :
134 0 : } else if (mu == 4) {
135 0 : val[0] = 1.; val[1] = -1.; val[2] = -1.; val[3] = -1.;
136 0 : index[0] = 0; index[1] = 1; index[2] = 2; index[3] = 3;
137 :
138 0 : } else if (mu == 5) {
139 0 : val[0] = -1.; val[1] = -1.; val[2] = 1.; val[3] = 1.;
140 0 : index[0] = 0; index[1] = 1; index[2] = 2; index[3] = 3;
141 0 : }
142 :
143 0 : }
144 :
145 : //--------------------------------------------------------------------------
146 :
147 : // Wave4 * GammaMatrix.
148 :
149 : Wave4 operator*(Wave4 w, GammaMatrix g) {
150 :
151 0 : complex w0 = w(g.index[0]);
152 0 : complex w1 = w(g.index[1]);
153 0 : complex w2 = w(g.index[2]);
154 0 : complex w3 = w(g.index[3]);
155 0 : w(0) = w0 * g.val[0];
156 0 : w(1) = w1 * g.val[1];
157 0 : w(2) = w2 * g.val[2];
158 0 : w(3) = w3 * g.val[3];
159 0 : return w;
160 :
161 0 : }
162 :
163 : //--------------------------------------------------------------------------
164 :
165 : // Scalar * GammaMatrix.
166 :
167 : GammaMatrix operator*(complex s, GammaMatrix g) {
168 :
169 0 : g.val[0] = s * g.val[0];
170 0 : g.val[1] = s*g.val[1];
171 0 : g.val[2] = s * g.val[2];
172 0 : g.val[3] = s*g.val[3];
173 0 : return g;
174 :
175 : }
176 :
177 : //--------------------------------------------------------------------------
178 :
179 : // I * Scalar - Gamma5.
180 :
181 : GammaMatrix operator-(complex s, GammaMatrix g) {
182 :
183 0 : g.val[0] = s - g.val[0];
184 0 : g.val[1] = s - g.val[1];
185 0 : g.val[2] = s - g.val[2];
186 0 : g.val[3] = s - g.val[3];
187 0 : return g;
188 :
189 : }
190 :
191 : //--------------------------------------------------------------------------
192 :
193 : // I * Scalar + Gamma5.
194 :
195 : GammaMatrix operator+(complex s, GammaMatrix g) {
196 :
197 0 : g.val[0] = s + g.val[0];
198 0 : g.val[1] = s + g.val[1];
199 0 : g.val[2] = s + g.val[2];
200 0 : g.val[3] = s + g.val[3];
201 0 : return g;
202 :
203 : }
204 :
205 : //--------------------------------------------------------------------------
206 :
207 : // Print GammaMatrix.
208 :
209 : ostream& operator<< (ostream& os, GammaMatrix g) {
210 :
211 0 : os << left << setprecision(2);
212 0 : for (int i = 0; i < 4; i++) {
213 0 : for (int j = 0; j < 4; j++) os << setw(20) << g(i,j);
214 0 : os << "\n";
215 : }
216 0 : 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 : // The spinors become ill-defined for p -> -pz and the polarization vectors
229 : // become ill-defined when pT -> 0. For these special cases limits are used.
230 :
231 : //--------------------------------------------------------------------------
232 :
233 : // Return wave vector for given helicity.
234 :
235 : Wave4 HelicityParticle::wave(int h) {
236 :
237 : // Create wave vector to return.
238 0 : Wave4 w;
239 :
240 : // Fermion (spin 1/2) spinor.
241 0 : if (spinType() == 2) {
242 :
243 : // Calculate helicity independent normalization.
244 0 : double P = pAbs();
245 0 : double n = sqrtpos(2*P*(P+pz()));
246 0 : bool aligned = P + pz() == 0;
247 :
248 : // Calculate eigenspinor basis.
249 0 : vector< vector<complex> > xi(2, vector<complex>(2));
250 : // Helicity -1
251 0 : xi[0][0] = aligned ? -1 : complex(-px(),py())/n;
252 0 : xi[0][1] = aligned ? 0 : (P+pz())/n;
253 : // Helicity +1
254 0 : xi[1][0] = aligned ? 0 : (P+pz())/n;
255 0 : xi[1][1] = aligned ? 1 : complex(px(),py())/n;
256 :
257 : // Calculate helicity dependent normalization.
258 0 : vector<double> omega(2);
259 0 : omega[0] = sqrtpos(e()-P);
260 0 : omega[1] = sqrtpos(e()+P);
261 0 : vector<double> hsign(2,1);
262 0 : hsign[0] = -1;
263 :
264 : // Create particle spinor.
265 0 : if (this->id() > 0) {
266 0 : w(0) = omega[!h] * xi[h][0];
267 0 : w(1) = omega[!h] * xi[h][1];
268 0 : w(2) = omega[h] * xi[h][0];
269 0 : w(3) = omega[h] * xi[h][1];
270 :
271 : // Create anti-particle spinor.
272 0 : } else {
273 0 : w(0) = hsign[!h] * omega[h] * xi[!h][0];
274 0 : w(1) = hsign[!h] * omega[h] * xi[!h][1];
275 0 : w(2) = hsign[h] * omega[!h] * xi[!h][0];
276 0 : w(3) = hsign[h] * omega[!h] * xi[!h][1];
277 : }
278 :
279 : // Boson (spin 1) polarization vector.
280 0 : } else if (spinType() == 3) {
281 0 : double P = pAbs();
282 0 : double PT = pT();
283 :
284 : // Create helicity +1 or -1 polarization vector.
285 0 : if (h >= 0 && h <= 1) {
286 0 : double hsign = h ? -1 : 1;
287 0 : if (P == 0) {
288 0 : w(0) = 0;
289 0 : w(1) = hsign / sqrt(2);
290 0 : w(2) = complex(0, 1/sqrt(2));
291 0 : w(3) = 0;
292 0 : } else if (PT == 0) {
293 0 : w(0) = 0;
294 0 : w(1) = hsign / sqrt(2);
295 0 : w(2) = complex(0, (pz() > 0 ? 1 : -1) / sqrt(2));
296 0 : w(3) = complex(-hsign * PT / P, 0) / sqrt(2);
297 0 : } else {
298 : w(0) = 0;
299 0 : w(1) = complex(hsign * px() * pz() / (P * PT), -py() / PT) / sqrt(2);
300 0 : w(2) = complex(hsign * py() * pz() / (P * PT), px() / PT) / sqrt(2);
301 0 : w(3) = complex(-hsign * PT / P, 0) / sqrt(2);
302 : }
303 :
304 : // Create helicity 0 polarization vector (ensure boson massive).
305 0 : } else if (h == 2 && spinStates() == 3) {
306 0 : if (P == 0) {
307 0 : w(0) = 0;
308 0 : w(1) = 0;
309 0 : w(2) = 0;
310 0 : w(3) = 1;
311 0 : } else {
312 0 : w(0) = P / m();
313 0 : w(1) = px() * e() / (m() * P);
314 0 : w(2) = py() * e() / (m() * P);
315 0 : w(3) = pz() * e() / (m() * P);
316 : }
317 : }
318 :
319 : // Unknown wave function.
320 0 : } else {
321 0 : w(0) = 0;
322 0 : w(1) = 0;
323 0 : w(2) = 0;
324 0 : w(3) = 0;
325 : }
326 :
327 : // Done.
328 : return w;
329 :
330 0 : }
331 :
332 : //--------------------------------------------------------------------------
333 :
334 : // Bar of a wave function.
335 :
336 : Wave4 HelicityParticle::waveBar(int h) {
337 :
338 0 : if (spinType() == 2) return conj(wave(h)) * GammaMatrix(0);
339 0 : else return conj(wave(h));
340 :
341 0 : }
342 :
343 : //--------------------------------------------------------------------------
344 :
345 : // Normalize the rho or D matrices.
346 :
347 : void HelicityParticle::normalize(vector< vector<complex> >& matrix) {
348 :
349 0 : complex trace = 0;
350 0 : for (unsigned int i = 0; i < matrix.size(); i++) trace += matrix[i][i];
351 0 : for (unsigned int i = 0; i < matrix.size(); i++) {
352 0 : for (unsigned int j = 0; j < matrix.size(); j++) {
353 0 : if (trace != complex(0,0)) matrix[i][j] /= trace;
354 0 : else matrix[i][j] = 1 / static_cast<double>(matrix.size());
355 : }
356 : }
357 :
358 0 : }
359 :
360 : //--------------------------------------------------------------------------
361 :
362 : // Return the number of spin states.
363 :
364 : int HelicityParticle::spinStates() {
365 :
366 0 : int sT = spinType();
367 0 : if (sT == 0) return 1;
368 0 : else if (sT != 2 && m() == 0) return sT - 1;
369 0 : else return sT;
370 :
371 0 : }
372 :
373 : //==========================================================================
374 :
375 : } // end namespace Pythia8
|
