StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
HelicityBasics.h
1 // HelicityBasics.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2020 Philip Ilten, Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL v2 or later, 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() : index() {};
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(), indexSave() { 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), indexSave() {
195  if (ptr) setPDEPtr( ptr->particleDataEntryPtr( idIn) );
196  initRhoD();
197  direction = 1; }
198  HelicityParticle(int idIn, int statusIn, int mother1In, int mother2In,
199  int daughter1In, int daughter2In, int colIn, int acolIn, Vec4 pIn,
200  double mIn = 0., double scaleIn = 0., ParticleData* ptr = 0)
201  : Particle(idIn, statusIn, mother1In, mother2In, daughter1In, daughter2In,
202  colIn, acolIn, pIn, mIn, scaleIn), indexSave() {
203  if (ptr) setPDEPtr( ptr->particleDataEntryPtr( idIn) );
204  initRhoD();
205  direction = 1; }
206  HelicityParticle(const Particle& ptIn, ParticleData* ptr = 0)
207  : Particle(ptIn) {
208  indexSave = ptIn.index();
209  if (ptr) setPDEPtr( ptr->particleDataEntryPtr( id()) );
210  initRhoD();
211  direction = 1; }
212 
213  // Methods.
214  Wave4 wave(int h);
215  Wave4 waveBar(int h);
216  void normalize(vector< vector<complex> >& m);
217  int spinStates();
218 
219  // Return and set mass (redefine from Particle).
220  double m() {return mSave;}
221  void m(double mIn) {mSave = mIn; initRhoD();}
222 
223  // Event record position (redefine from Particle).
224  int index() const {return indexSave;}
225  void index(int indexIn) {indexSave = indexIn;}
226 
227  // Flag for whether particle is incoming (-1) or outgoing (1).
228  int direction;
229 
230  // Helicity density matrix.
231  vector< vector<complex> > rho;
232 
233  // Decay matrix.
234  vector< vector<complex> > D;
235 
236 private:
237 
238  // Initialize the helicity density and decay matrix.
239  void initRhoD() {
240  rho = vector< vector<complex> >(spinStates(),
241  vector<complex>(spinStates(), 0));
242  D = vector< vector<complex> >(spinStates(),
243  vector<complex>(spinStates(), 0));
244  for (int i = 0; i < spinStates(); i++) {
245  rho[i][i] = 1.0 / spinStates(); D[i][i] = 1;}
246  }
247 
248  // Particle index in the event record.
249  int indexSave;
250 
251 };
252 
253 //==========================================================================
254 
255 } // end namespace Pythia8
256 
257 #endif // end Pythia8_HelicityBasics_H