PYTHIA8/pythia8170/src/HelicityBasics.cxx

   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