StRoot  1
HelicityBasics.h
1 // HelicityBasics.h 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
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 "Pythia8/Basics.h"
12 #include "Pythia8/Event.h"
13 #include "Pythia8/PythiaComplex.h"
14 #include "Pythia8/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 * double.
62  Wave4 operator*(double s) {return Wave4(val[0] * s, val[1] * s,
63  val[2] * s, val[3] * s);}
64
65  // double * Wave4.
66  friend Wave4 operator*(double s, const Wave4& w);
67
68  // Wave4 / complex.
69  Wave4 operator/(complex s) {return Wave4(val[0] / s, val[1] / s,
70  val[2] / s, val[3] / s);}
71
72  // Wave4 / double.
73  Wave4 operator/(double s) {return Wave4(val[0] / s, val[1] / s,
74  val[2]/s, val[3]/s);}
75
76  // Complex conjugate.
77  friend Wave4 conj(Wave4 w);
78
79  // Permutation operator.
80  friend Wave4 epsilon(Wave4 w1, Wave4 w2, Wave4 w3);
81
82  // Invariant squared mass for REAL Wave4 (to save time).
83  friend double m2(Wave4 w);
84  friend double m2(Wave4 w1, Wave4 w2);
85
86  // Wave4 * GammaMatrix multiplication is defined in the GammaMatrix class.
87
88  // Print a Wave4 vector.
89  friend ostream& operator<<(ostream& output, Wave4 w);
90
91 protected:
92
93  complex val[4];
94
95 };
96
97 //--------------------------------------------------------------------------
98
99 // Namespace function declarations; friends of Wave4 class.
100 Wave4 operator*(complex s, const Wave4& w);
101 Wave4 operator*(double s, const Wave4& w);
102 Wave4 conj(Wave4 w);
103 Wave4 epsilon(Wave4 w1, Wave4 w2, Wave4 w3);
104 double m2(Wave4 w);
105 double m2(Wave4 w1, Wave4 w2);
106 ostream& operator<< (ostream& os, Wave4 w);
107
108 //==========================================================================
109
110 // The GammaMatrix class is a special sparse matrix class used to write
111 // helicity matrix elements in conjuction with the Wave4 class. Note that
112 // only left to right multplication of Wave4 vectors with the GammaMatrix
113 // class is allowed. Additionally, subtracting a scalar from a GammaMatrix
114 // (or subtracting a GammaMatrix from a scalar) subtracts the scalar from
115 //each non-zero element of the GammaMatrix. This is designed specifically
116 // with the (1 - gamma^5) structure of matrix elements in mind.
117
118 class GammaMatrix {
119
120 public:
121
122  // Constructors and destructor.
123  GammaMatrix() {};
124  GammaMatrix(int mu);
125  ~GammaMatrix() {};
126
127  // Access an element of the matrix.
128  complex& operator() (int I, int J) {if (index[J] == I) return val[J];
129  else return COMPLEXZERO; }
130
131  // Wave4 * GammaMatrix.
132  friend Wave4 operator*(Wave4 w, GammaMatrix g);
133
134  // GammaMatrix * Scalar.
135  GammaMatrix operator*(complex s) {val[0] = s*val[0]; val[1] = s*val[1];
136  val[2] = s*val[2]; val[3] = s*val[3]; return *this;}
137
138  // Scalar * GammaMatrix.
139  friend GammaMatrix operator*(complex s, GammaMatrix g);
140
141  // Gamma5 - I * Scalar.
142  GammaMatrix operator-(complex s) {val[0] = val[0] - s; val[1] = val[1] - s;
143  val[2] = val[2] - s; val[3] = val[3] - s; return *this;}
144
145  // I * Scalar - Gamma5.
146  friend GammaMatrix operator-(complex s, GammaMatrix g);
147
148  // Gamma5 + I * Scalar
149  GammaMatrix operator+(complex s) {val[0] = val[0] + s; val[1] = val[1] + s;
150  val[2] = val[2] + s; val[3] = val[3] + s; return *this;}
151
152  // I * Scalar + Gamma5
153  friend GammaMatrix operator+(complex s, GammaMatrix g);
154
155  // << GammaMatrix.
156  friend ostream& operator<< (ostream& os, GammaMatrix g);
157
158 protected:
159
160  complex val[4];
161  int index[4];
162
163  // Need to define complex 0 as a variable for operator() to work.
164  complex COMPLEXZERO;
165
166 };
167
168 //--------------------------------------------------------------------------
169
170 // Namespace function declarations; friends of GammaMatrix class.
171 Wave4 operator*(Wave4 w, GammaMatrix g);
172 GammaMatrix operator*(complex s, GammaMatrix g);
173 GammaMatrix operator-(complex s, GammaMatrix g);
174 GammaMatrix operator+(complex s, GammaMatrix g);
175 ostream& operator<< (ostream& os, GammaMatrix g);
176
177 //==========================================================================
178
179 // Helicity particle class containing helicity information, derived from
180 // particle base class.
181
182 class HelicityParticle : public Particle {
183
184 public:
185
186  // Constructors.
187  HelicityParticle() : Particle() { direction = 1;}
188  HelicityParticle(int idIn, int statusIn = 0, int mother1In = 0,
189  int mother2In = 0, int daughter1In = 0, int daughter2In = 0,
190  int colIn = 0, int acolIn = 0, double pxIn = 0.,
191  double pyIn = 0., double pzIn = 0., double eIn = 0.,
192  double mIn = 0., double scaleIn = 0., ParticleData* ptr = 0)
193  : Particle(idIn, statusIn, mother1In, mother2In, daughter1In, daughter2In,
194  colIn, acolIn, pxIn, pyIn, pzIn, eIn, mIn, scaleIn) {
195  if (ptr) setPDEPtr( ptr->particleDataEntryPtr( idIn) );
196  rho = vector< vector<complex> >(spinStates(),
197  vector<complex>(spinStates(), 0));
198  D = vector< vector<complex> >(spinStates(),
199  vector<complex>(spinStates(), 0));
200  for (int i = 0; i < spinStates(); i++) { rho[i][i] = 0.5; D[i][i] = 1.;}
201  direction = 1; }
202  HelicityParticle(int idIn, int statusIn, int mother1In, int mother2In,
203  int daughter1In, int daughter2In, int colIn, int acolIn, Vec4 pIn,
204  double mIn = 0., double scaleIn = 0., ParticleData* ptr = 0)
205  : Particle(idIn, statusIn, mother1In, mother2In, daughter1In, daughter2In,
206  colIn, acolIn, pIn, mIn, scaleIn) {
207  if (ptr) setPDEPtr( ptr->particleDataEntryPtr( idIn) );
208  rho = vector< vector<complex> >(spinStates(),
209  vector<complex>(spinStates(), 0));
210  D = vector< vector<complex> >(spinStates(),
211  vector<complex>(spinStates(), 0));
212  for (int i = 0; i < spinStates(); i++) { rho[i][i] = 0.5; D[i][i] = 1;}
213  direction = 1; }
214  HelicityParticle(const Particle& ptIn, ParticleData* ptr = 0)
215  : Particle(ptIn) {
216  if (ptr) setPDEPtr( ptr->particleDataEntryPtr( id()) );
217  rho = vector< vector<complex> >(spinStates(),
218  vector<complex>(spinStates(), 0));
219  D = vector< vector<complex> >(spinStates(),
220  vector<complex>(spinStates(), 0));
221  for (int i = 0; i < spinStates(); i++) { rho[i][i] = 0.5; D[i][i] = 1;}
222  direction = 1; }
223
224  // Methods.
225  Wave4 wave(int h);
226  Wave4 waveBar(int h);
227  void normalize(vector< vector<complex> >& m);
228  int spinStates();
229
230  // Event record position.
231  int idx;
232
233  // Flag for whether particle is incoming (-1) or outgoing (1).
234  int direction;
235
236  // Helicity density matrix.
237  vector< vector<complex> > rho;
238
239  // Decay matrix.
240  vector< vector<complex> > D;
241
242 private:
243
244  // Constants: could only be changed in the code itself.
245  static const double TOLERANCE;
246
247 };
248
249 //==========================================================================
250
251 } // end namespace Pythia8
252
253 #endif // end Pythia8_HelicityBasics_H