Line data Source code
1 : // HelicityBasics.h is a part of the PYTHIA event generator.
2 : // Copyright (C) 2015 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 0 : Wave4() {};
30 0 : Wave4(complex v0, complex v1, complex v2, complex v3) {val[0] = v0;
31 0 : val[1] = v1; val[2] = v2; val[3] = v3;}
32 0 : Wave4(Vec4 v) {val[0] = v.e(); val[1] = v.px(); val[2] = v.py();
33 0 : val[3] = v.pz();}
34 0 : ~Wave4() {};
35 :
36 : // Access an element of the wave vector.
37 0 : complex& operator() (int i) {return val[i];}
38 :
39 : // Wave4 + Wave4.
40 0 : Wave4 operator+(Wave4 w) {return Wave4( val[0] + w.val[0],
41 0 : val[1] + w.val[1], val[2] + w.val[2], val[3] + w.val[3]);}
42 :
43 : // Wave4 - Wave4.
44 0 : Wave4 operator-(Wave4 w) {return Wave4( val[0] - w.val[0],
45 0 : 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 0 : complex operator*(Wave4 w) {return val[0] * w.val[0]
52 0 : + val[1] * w.val[1] + val[2] * w.val[2] + val[3] * w.val[3];}
53 :
54 : // Wave4 * complex.
55 0 : Wave4 operator*(complex s) {return Wave4(val[0] * s, val[1] * s,
56 0 : val[2] * s, val[3] * s);}
57 :
58 : // complex * Wave4.
59 : friend Wave4 operator*(complex s, const Wave4& w);
60 :
61 : // Wave4 * double.
62 0 : Wave4 operator*(double s) {return Wave4(val[0] * s, val[1] * s,
63 0 : val[2] * s, val[3] * s);}
64 :
65 : // double * Wave4.
66 : friend Wave4 operator*(double s, const Wave4& w);
67 :
68 : // Wave4 / complex.
69 0 : Wave4 operator/(complex s) {return Wave4(val[0] / s, val[1] / s,
70 0 : 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 0 : ~GammaMatrix() {};
126 :
127 : // Access an element of the matrix.
128 0 : complex& operator() (int I, int J) {if (index[J] == I) return val[J];
129 0 : 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 0 : class HelicityParticle : public Particle {
183 :
184 : public:
185 :
186 : // Constructors.
187 0 : HelicityParticle() : Particle() { direction = 1;}
188 0 : 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 0 : : Particle(idIn, statusIn, mother1In, mother2In, daughter1In, daughter2In,
194 0 : colIn, acolIn, pxIn, pyIn, pzIn, eIn, mIn, scaleIn) {
195 0 : if (ptr) setPDEPtr( ptr->particleDataEntryPtr( idIn) );
196 0 : initRhoD();
197 0 : direction = 1; }
198 0 : 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 0 : : Particle(idIn, statusIn, mother1In, mother2In, daughter1In, daughter2In,
202 0 : colIn, acolIn, pIn, mIn, scaleIn) {
203 0 : if (ptr) setPDEPtr( ptr->particleDataEntryPtr( idIn) );
204 0 : initRhoD();
205 0 : direction = 1; }
206 0 : HelicityParticle(const Particle& ptIn, ParticleData* ptr = 0)
207 0 : : Particle(ptIn) {
208 0 : indexSave = ptIn.index();
209 0 : if (ptr) setPDEPtr( ptr->particleDataEntryPtr( id()) );
210 0 : initRhoD();
211 0 : 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 0 : double m() {return mSave;}
221 0 : void m(double mIn) {mSave = mIn; initRhoD();}
222 :
223 : // Event record position (redefine from Particle).
224 0 : int index() const {return indexSave;}
225 0 : 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 0 : rho = vector< vector<complex> >(spinStates(),
241 0 : vector<complex>(spinStates(), 0));
242 0 : D = vector< vector<complex> >(spinStates(),
243 0 : vector<complex>(spinStates(), 0));
244 0 : for (int i = 0; i < spinStates(); i++) {
245 0 : rho[i][i] = 1.0 / spinStates(); D[i][i] = 1;}
246 0 : }
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
|
