StRoot  1
HelicityBasics.cc
1 // HelicityBasics.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2014 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
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  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 1;
363  else if (sT != 2 && m() < TOLERANCE) return sT - 1;
364  else return sT;
365
366  }
367
368 //==========================================================================
369
370 } // end namespace Pythia8