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