StRoot  1
HelicityBasics.h
1 // HelicityBasics.h 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 // Header file for a number of helper classes used in tau decays.
7
8 #ifndef Pythia8_HelicityBasics_H
9 #define Pythia8_HelicityBasics_H
10
11 #include "Basics.h"
12 #include "Event.h"
13 #include "PythiaComplex.h"
14 #include "PythiaStdlib.h"
15
16 namespace Pythia8 {
17
18 //==========================================================================
19
20 // The Wave4 class provides a class for complex four-vector wave functions.
21 // The Wave4 class can be multiplied with the GammaMatrix class to allow
22 // for the writing of helicity matrix elements.
23
24 class Wave4 {
25
26 public:
27
28  // Constructors and destructor.
29  Wave4() {};
30  Wave4(complex v0, complex v1, complex v2, complex v3) {val[0] = v0;
31  val[1] = v1; val[2] = v2; val[3] = v3;}
32  Wave4(Vec4 v) {val[0] = v.e(); val[1] = v.px(); val[2] = v.py();
33  val[3] = v.pz();}
34  ~Wave4() {};
35
36  // Access an element of the wave vector.
37  complex& operator() (int i) {return val[i];}
38
39  // Wave4 + Wave4.
40  Wave4 operator+(Wave4 w) {return Wave4( val[0] + w.val[0],
41  val[1] + w.val[1], val[2] + w.val[2], val[3] + w.val[3]);}
42
43  // Wave4 - Wave4.
44  Wave4 operator-(Wave4 w) {return Wave4( val[0] - w.val[0],
45  val[1] - w.val[1], val[2] - w.val[2], val[3] - w.val[3]);}
46
47  // - Wave4.
48  Wave4 operator-() {return Wave4(-val[0], -val[1], -val[2], -val[3]);}
49
50  // Wave4 * Wave4.
51  complex operator*(Wave4 w) {return val[0] * w.val[0]
52  + val[1] * w.val[1] + val[2] * w.val[2] + val[3] * w.val[3];}
53
54  // Wave4 * complex.
55  Wave4 operator*(complex s) {return Wave4(val[0] * s, val[1] * s,
56  val[2] * s, val[3] * s);}
57
58  // complex * Wave4.
59  friend Wave4 operator*(complex s, const Wave4& w);
60
61  // Wave4 / complex.
62  Wave4 operator/(complex s) {return Wave4(val[0] / s, val[1] / s,
63  val[2] / s, val[3] / s);}
64
65  // Wave4 / double.
66  Wave4 operator/(double s) {return Wave4(val[0] / s, val[1] / s,
67  val[2]/s, val[3]/s);}
68
69  // Complex conjugate.
70  friend Wave4 conj(Wave4 w) ;
71
72  // Invariant squared mass for REAL Wave4 (to save time).
73  friend double m2(Wave4 w1, Wave4 w2);
74
75  // Wave4 * GammaMatrix multiplication is defined in the GammaMatrix class.
76
77  // Print a Wave4 vector.
78  friend ostream& operator<<(ostream& output, Wave4 w);
79
80 protected:
81
82  complex val[4];
83
84 };
85
86 //--------------------------------------------------------------------------
87
88 // Namespace function declarations; friends of Wave4 class.
89 Wave4 operator*(complex s, const Wave4& w);
90 Wave4 conj(Wave4 w);
91 double m2(Wave4 w1, Wave4 w2);
92 ostream& operator<< (ostream& os, Wave4 w);
93
94 //==========================================================================
95
96 // The GammaMatrix class is a special sparse matrix class used to write
97 // helicity matrix elements in conjuction with the Wave4 class. Note that
98 // only left to right multplication of Wave4 vectors with the GammaMatrix
99 // class is allowed. Additionally, subtracting a scalar from a GammaMatrix
100 // (or subtracting a GammaMatrix from a scalar) subtracts the scalar from
101 //each non-zero element of the GammaMatrix. This is designed specifically
102 // with the (1 - gamma^5) structure of matrix elements in mind.
103
104 class GammaMatrix {
105
106 public:
107
108  // Constructors and destructor.
109  GammaMatrix() {};
110  GammaMatrix(int mu);
111  ~GammaMatrix() {};
112
113  // Access an element of the matrix.
114  complex& operator() (int I, int J) {if (index[J] == I) return val[J];
115  else return COMPLEXZERO; }
116
117  // Wave4 * GammaMatrix.
118  friend Wave4 operator*(Wave4 w, GammaMatrix g);
119
120  // GammaMatrix * Scalar.
121  GammaMatrix operator*(complex s) {val[0] = s*val[0]; val[1] = s*val[1];
122  val[2] = s*val[2]; val[3] = s*val[3]; return *this;}
123
124  // Scalar * GammaMatrix.
125  friend GammaMatrix operator*(complex s, GammaMatrix g);
126
127  // Gamma5 - I * Scalar.
128  GammaMatrix operator-(complex s) {val[0] = val[0] - s; val[1] = val[1] - s;
129  val[2] = val[2] - s; val[3] = val[3] - s; return *this;}
130
131  // I * Scalar - Gamma5.
132  friend GammaMatrix operator-(complex s, GammaMatrix g);
133
134  // Gamma5 + I * Scalar
135  GammaMatrix operator+(complex s) {val[0] = val[0] + s; val[1] = val[1] + s;
136  val[2] = val[2] + s; val[3] = val[3] + s; return *this;}
137
138  // I * Scalar + Gamma5
139  friend GammaMatrix operator+(complex s, GammaMatrix g);
140
141  // << GammaMatrix.
142  friend ostream& operator<< (ostream& os, GammaMatrix g);
143
144 protected:
145
146  complex val[4];
147  int index[4];
148
149  // Need to define complex 0 as a variable for operator() to work.
150  complex COMPLEXZERO;
151
152 };
153
154 //--------------------------------------------------------------------------
155
156 // Namespace function declarations; friends of GammaMatrix class.
157 Wave4 operator*(Wave4 w, GammaMatrix g);
158 GammaMatrix operator*(complex s, GammaMatrix g);
159 GammaMatrix operator-(complex s, GammaMatrix g);
160 GammaMatrix operator+(complex s, GammaMatrix g);
161 ostream& operator<< (ostream& os, GammaMatrix g);
162
163 //==========================================================================
164
165 // Helicity particle class containing helicity information, derived from
166 // particle base class.
167
168 class HelicityParticle : public Particle {
169
170 public:
171
172  // Constructors.
173  HelicityParticle() : Particle() { direction = 1;}
174  HelicityParticle(int idIn, int statusIn = 0, int mother1In = 0,
175  int mother2In = 0, int daughter1In = 0, int daughter2In = 0,
176  int colIn = 0, int acolIn = 0, double pxIn = 0.,
177  double pyIn = 0., double pzIn = 0., double eIn = 0.,
178  double mIn = 0., double scaleIn = 0., ParticleData* ptr = 0)
179  : Particle(idIn, statusIn, mother1In, mother2In, daughter1In, daughter2In,
180  colIn, acolIn, pxIn, pyIn, pzIn, eIn, mIn, scaleIn) {
181  if (ptr) { setPDTPtr(ptr); setPDEPtr(); }
182  rho = vector< vector<complex> >(spinType(),
183  vector<complex>(spinType(), 0));
184  D = vector< vector<complex> >(spinType(),
185  vector<complex>(spinType(), 0));
186  for (int i = 0; i < spinType(); i++) { rho[i][i] = 0.5; D[i][i] = 1.;}
187  direction = 1; }
188  HelicityParticle(int idIn, int statusIn, int mother1In, int mother2In,
189  int daughter1In, int daughter2In, int colIn, int acolIn, Vec4 pIn,
190  double mIn = 0., double scaleIn = 0., ParticleData* ptr = 0)
191  : Particle(idIn, statusIn, mother1In, mother2In, daughter1In, daughter2In,
192  colIn, acolIn, pIn, mIn, scaleIn) {
193  if (ptr) { setPDTPtr(ptr); setPDEPtr();}
194  rho = vector< vector<complex> >(spinType(),
195  vector<complex>(spinType(), 0));
196  D = vector< vector<complex> >(spinType(),
197  vector<complex>(spinType(), 0));
198  for (int i = 0; i < spinType(); i++) { rho[i][i] = 0.5; D[i][i] = 1;}
199  direction = 1; }
200  HelicityParticle(const Particle& ptIn, ParticleData* ptr = 0)
201  : Particle(ptIn) {
202  if (ptr) { setPDTPtr(ptr); setPDEPtr();}
203  rho = vector< vector<complex> >(spinType(),
204  vector<complex>(spinType(), 0));
205  D = vector< vector<complex> >(spinType(),
206  vector<complex>(spinType(), 0));
207  for (int i = 0; i < spinType(); i++) { rho[i][i] = 0.5; D[i][i] = 1;}
208  direction = 1; }
209
210  // Methods.
211  Wave4 wave(int h);
212  Wave4 waveBar(int h);
213  void normalize(vector< vector<complex> >& m);
214
215  // Event record position.
216  int idx;
217
218  // Flag for whether particle is incoming (-1) or outgoing (1).
219  int direction;
220
221  // Helicity density matrix.
222  vector< vector<complex> > rho;
223
224  // Decay matrix.
225  vector< vector<complex> > D;
226
227 private:
228
229  // Constants: could only be changed in the code itself.
230  static const double TOLERANCE;
231
232 };
233
234 //==========================================================================
235
236 } // end namespace Pythia8
237
238 #endif // end Pythia8_HelicityBasics_H