Line data Source code
1 : // HelicityMatrixElements.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 physics classes used in tau decays.
7 :
8 : #ifndef Pythia8_HelicityMatrixElements_H
9 : #define Pythia8_HelicityMatrixElements_H
10 :
11 : #include "Pythia8/Basics.h"
12 : #include "Pythia8/Event.h"
13 : #include "Pythia8/HelicityBasics.h"
14 : #include "Pythia8/PythiaComplex.h"
15 : #include "Pythia8/PythiaStdlib.h"
16 : #include "Pythia8/StandardModel.h"
17 :
18 : namespace Pythia8 {
19 :
20 : //==========================================================================
21 :
22 : // The helicity matrix element class.
23 :
24 : class HelicityMatrixElement {
25 :
26 : public:
27 :
28 : // Constructor and destructor.
29 0 : HelicityMatrixElement() {};
30 0 : virtual ~HelicityMatrixElement() {};
31 :
32 : // Initialize the physics matrices and pointers.
33 : virtual void initPointers(ParticleData*, Couplings*, Settings* = 0);
34 :
35 : // Initialize the channel.
36 : virtual HelicityMatrixElement* initChannel(vector<HelicityParticle>&);
37 :
38 : // Calculate the matrix element weight for a decay.
39 : virtual double decayWeight(vector<HelicityParticle>&);
40 :
41 : // Calculate the maximum matrix element decay weight.
42 : virtual double decayWeightMax(vector<HelicityParticle>&)
43 0 : {return DECAYWEIGHTMAX;}
44 :
45 : // Calculate the helicity matrix element.
46 0 : virtual complex calculateME(vector<int>){return complex(0,0);}
47 :
48 : // Calculate the decay matrix for a particle.
49 : virtual void calculateD(vector<HelicityParticle>&);
50 :
51 : // Calculate the density matrix for a particle.
52 : virtual void calculateRho(unsigned int, vector<HelicityParticle>&);
53 :
54 : // Set a fermion line.
55 : void setFermionLine(int, HelicityParticle&, HelicityParticle&);
56 :
57 : // Calculate Breit-Wigner's with running widths and fixed.
58 : virtual complex breitWigner(double s, double M, double G);
59 : virtual complex sBreitWigner(double m0, double m1, double s,
60 : double M, double G);
61 : virtual complex pBreitWigner(double m0, double m1, double s,
62 : double M, double G);
63 : virtual complex dBreitWigner(double m0, double m1, double s,
64 : double M, double G);
65 :
66 : protected:
67 :
68 : // Maximum decay weight. WARNING: hardcoded constant.
69 : double DECAYWEIGHTMAX;
70 :
71 : // Physics matrices.
72 : vector< GammaMatrix > gamma;
73 :
74 : // Particle map vector.
75 : vector< int > pMap;
76 :
77 : // Particle ID vector.
78 : vector< int > pID;
79 :
80 : // Particle mass vector.
81 : vector< double > pM;
82 :
83 : // Wave functions.
84 : vector< vector< Wave4 > > u;
85 :
86 : // Initialize the constants for the matrix element (called by initChannel).
87 0 : virtual void initConstants() {};
88 :
89 : // Initialize the wave functions (called by decayWeight and calculateRho/D).
90 0 : virtual void initWaves(vector<HelicityParticle>&) {};
91 :
92 : // Pointer to particle data.
93 : ParticleData* particleDataPtr;
94 :
95 : // Pointer to Standard Model constants.
96 : Couplings* couplingsPtr;
97 :
98 : // Pointer to Settings.
99 : Settings* settingsPtr;
100 :
101 : private:
102 :
103 : // Recursive sub-method to calculate the density matrix for a particle.
104 : void calculateRho(unsigned int, vector<HelicityParticle>&,
105 : vector<int>&, vector<int>&, unsigned int);
106 :
107 : // Recursive sub-method to calculate the decay matrix for a particle.
108 : void calculateD(vector<HelicityParticle>&, vector<int>&, vector<int>&,
109 : unsigned int);
110 :
111 : // Recursive sub-method to calculate the matrix element weight for a decay.
112 : void decayWeight(vector<HelicityParticle>&, vector<int>&, vector<int>&,
113 : complex&, unsigned int);
114 :
115 : // Calculate the product of the decay matrices for a hard process.
116 : complex calculateProductD(unsigned int, unsigned int,
117 : vector<HelicityParticle>&, vector<int>&, vector<int>&);
118 :
119 : // Calculate the product of the decay matrices for a decay process.
120 : complex calculateProductD(vector<HelicityParticle>&,
121 : vector<int>&, vector<int>&);
122 :
123 : };
124 :
125 : //==========================================================================
126 :
127 : // Helicity matrix element for the hard process of two fermions -> W/W' ->
128 : // two fermions.
129 :
130 0 : class HMETwoFermions2W2TwoFermions : public HelicityMatrixElement {
131 :
132 : public:
133 :
134 : void initConstants();
135 :
136 : void initWaves(vector<HelicityParticle>&);
137 :
138 : complex calculateME(vector<int>);
139 :
140 : private:
141 :
142 : // Vector and axial couplings.
143 : double p0CA, p2CA, p0CV, p2CV;
144 :
145 : };
146 :
147 : //==========================================================================
148 :
149 : // Helicity matrix element for the hard process of two fermions ->
150 : // photon/Z/Z' -> two fermions.
151 :
152 0 : class HMETwoFermions2GammaZ2TwoFermions : public HelicityMatrixElement {
153 :
154 : public:
155 :
156 : void initConstants();
157 :
158 : void initWaves(vector<HelicityParticle>&);
159 :
160 : complex calculateME(vector<int>);
161 :
162 : private:
163 :
164 : // Return gamma element for the helicity matrix element.
165 : complex calculateGammaME(vector<int>);
166 :
167 : // Return Z/Z' element for helicity matrix element.
168 : complex calculateZME(vector<int>, double, double, double, double,
169 : double, double);
170 :
171 : // Return the Z' vector or axial coupling for a fermion.
172 : double zpCoupling(int id, string type);
173 :
174 : // Vector and axial couplings.
175 : double p0CAZ, p2CAZ, p0CVZ, p2CVZ, p0CAZp, p2CAZp, p0CVZp, p2CVZp;
176 :
177 : // Weinberg angle, Z/Z' width (on-shell), Z/Z' mass (on-shell), and s.
178 : double cos2W, sin2W, zG, zM, zpG, zpM, s;
179 :
180 : // Fermion line charge.
181 : double p0Q, p2Q;
182 :
183 : // Bool whether the incoming fermions are oriented with the z-axis.
184 : bool zaxis, includeGamma, includeZ, includeZp;
185 :
186 : };
187 :
188 : //==========================================================================
189 :
190 : // Helicity matrix element for the hard process of X -> two fermions.
191 :
192 0 : class HMEX2TwoFermions : public HelicityMatrixElement {
193 :
194 : public:
195 :
196 : void initWaves(vector<HelicityParticle>&);
197 :
198 : };
199 :
200 : //==========================================================================
201 :
202 : // Helicity matrix element for the hard process of W/W' -> two fermions.
203 :
204 0 : class HMEW2TwoFermions : public HMEX2TwoFermions {
205 :
206 : public:
207 :
208 : void initConstants();
209 :
210 : complex calculateME(vector<int>);
211 :
212 : private:
213 :
214 : // Vector and axial couplings.
215 : double p2CA, p2CV;
216 :
217 : };
218 :
219 : //==========================================================================
220 :
221 : // Helicity matrix element for the hard process of photon -> two fermions.
222 :
223 0 : class HMEGamma2TwoFermions : public HMEX2TwoFermions {
224 :
225 : public:
226 :
227 : complex calculateME(vector<int>);
228 :
229 : };
230 :
231 : //==========================================================================
232 :
233 : // Helicity matrix element for the hard process of Z/Z' -> two fermions.
234 :
235 0 : class HMEZ2TwoFermions : public HMEX2TwoFermions {
236 :
237 : public:
238 :
239 : void initConstants();
240 :
241 : complex calculateME(vector<int>);
242 :
243 : private:
244 :
245 : // Return the Z' vector or axial coupling for a fermion.
246 : double zpCoupling(int id, string type);
247 :
248 : // Vector and axial couplings.
249 : double p2CA, p2CV;
250 :
251 : };
252 :
253 : //==========================================================================
254 :
255 : // Helicity matrix element for the decay of a Higgs -> two fermions.
256 :
257 : // Because the Higgs is spin zero the Higgs production mechanism is not
258 : // needed for calculating helicity density matrices. However, the CP mixing
259 : // is needed.
260 :
261 0 : class HMEHiggs2TwoFermions : public HelicityMatrixElement {
262 :
263 : public:
264 :
265 : void initConstants();
266 :
267 : void initWaves(vector<HelicityParticle>&);
268 :
269 : complex calculateME(vector<int>);
270 :
271 : private:
272 :
273 : // Coupling constants of the fermions with the Higgs.
274 : complex p2CA, p2CV;
275 :
276 : };
277 :
278 : //==========================================================================
279 :
280 : // Base class for all tau decay helicity matrix elements.
281 :
282 0 : class HMETauDecay : public HelicityMatrixElement {
283 :
284 : public:
285 :
286 : virtual void initWaves(vector<HelicityParticle>&);
287 :
288 : virtual complex calculateME(vector<int>);
289 :
290 : virtual double decayWeightMax(vector<HelicityParticle>&);
291 :
292 : protected:
293 :
294 0 : virtual void initHadronicCurrent(vector<HelicityParticle>&) {};
295 :
296 : virtual void calculateResonanceWeights(vector<double>&, vector<double>&,
297 : vector<complex>&);
298 :
299 : };
300 :
301 : //==========================================================================
302 :
303 : // Helicity matrix element for a tau decaying into a single scalar meson.
304 :
305 0 : class HMETau2Meson : public HMETauDecay {
306 :
307 : public:
308 :
309 : void initConstants();
310 :
311 : void initHadronicCurrent(vector<HelicityParticle>&);
312 :
313 : };
314 :
315 : //==========================================================================
316 :
317 : // Helicity matrix element for a tau decaying into two leptons.
318 :
319 0 : class HMETau2TwoLeptons : public HMETauDecay {
320 :
321 : public:
322 :
323 : void initConstants();
324 :
325 : void initWaves(vector<HelicityParticle>&);
326 :
327 : complex calculateME(vector<int>);
328 :
329 : };
330 :
331 : //==========================================================================
332 :
333 : // Helicity matrix element for a tau decaying into two mesons through a
334 : // vector meson resonance.
335 :
336 0 : class HMETau2TwoMesonsViaVector : public HMETauDecay {
337 :
338 : public:
339 :
340 : void initConstants();
341 :
342 : void initHadronicCurrent(vector<HelicityParticle>&);
343 :
344 : private:
345 :
346 : // Resonance masses, widths, and weights.
347 : vector<double> vecM, vecG, vecP, vecA;
348 : vector<complex> vecW;
349 :
350 : };
351 :
352 : //==========================================================================
353 :
354 : // Helicity matrix element for a tau decay into two mesons through a vector
355 : // or scalar meson resonance.
356 :
357 0 : class HMETau2TwoMesonsViaVectorScalar : public HMETauDecay {
358 :
359 : public:
360 :
361 : void initConstants();
362 :
363 : void initHadronicCurrent(vector<HelicityParticle>&);
364 :
365 : private:
366 :
367 : // Coupling to vector and scalar resonances.
368 : double scaC, vecC;
369 :
370 : // Resonance masses, widths, and weights.
371 : vector<double> scaM, scaG, scaP, scaA, vecM, vecG, vecP, vecA;
372 : vector<complex> scaW, vecW;
373 :
374 : };
375 :
376 : //==========================================================================
377 :
378 : // Helicity matrix element for a tau decay into three mesons (base class).
379 :
380 0 : class HMETau2ThreeMesons : public HMETauDecay {
381 :
382 : public:
383 :
384 : void initConstants();
385 :
386 : void initHadronicCurrent(vector<HelicityParticle>&);
387 :
388 : protected:
389 :
390 : // Decay mode of the tau.
391 : enum Mode{Pi0Pi0Pim, PimPimPip, Pi0PimK0b, PimPipKm, Pi0PimEta, PimKmKp,
392 : Pi0K0Km, KlPimKs, Pi0Pi0Km, KlKlPim, PimKsKs, PimK0bK0, Uknown};
393 : Mode mode;
394 :
395 : // Initialize decay mode and resonance constants (called by initConstants).
396 : virtual void initMode();
397 0 : virtual void initResonances() {;}
398 :
399 : // Initialize the momenta.
400 : virtual void initMomenta(vector<HelicityParticle>&);
401 :
402 : // Center of mass energies and momenta.
403 : double s1, s2, s3, s4;
404 : Wave4 q, q2, q3, q4;
405 :
406 : // Stored a1 Breit-Wigner (for speed).
407 : complex a1BW;
408 :
409 : // Form factors.
410 0 : virtual complex F1() {return complex(0, 0);}
411 0 : virtual complex F2() {return complex(0, 0);}
412 0 : virtual complex F3() {return complex(0, 0);}
413 0 : virtual complex F4() {return complex(0, 0);}
414 :
415 : // Phase space and Breit-Wigner for the a1.
416 : virtual double a1PhaseSpace(double);
417 : virtual complex a1BreitWigner(double);
418 :
419 : // Sum running p and fixed width Breit-Wigner resonances.
420 : complex T(double m0, double m1, double s,
421 : vector<double>& M, vector<double>& G, vector<double>& W);
422 : complex T(double s, vector<double>& M, vector<double>& G, vector<double>& W);
423 :
424 : };
425 :
426 : //==========================================================================
427 :
428 : // Helicity matrix element for a tau decay into three pions.
429 :
430 0 : class HMETau2ThreePions : public HMETau2ThreeMesons {
431 :
432 : private:
433 :
434 : void initResonances();
435 :
436 : // Resonance masses, widths, and weights.
437 : vector<double> rhoM, rhoG, rhoPp, rhoAp, rhoPd, rhoAd;
438 : double f0M, f0G, f0P, f0A, f2M, f2G, f2P, f2A;
439 : double sigM, sigG, sigP, sigA;
440 : vector<complex> rhoWp, rhoWd;
441 : complex f0W, f2W, sigW;
442 :
443 : // Form factors.
444 : complex F1();
445 : complex F2();
446 : complex F3();
447 :
448 : // Running width and Breit-Wigner for the a1.
449 : double a1PhaseSpace(double);
450 : complex a1BreitWigner(double);
451 :
452 : };
453 :
454 : //==========================================================================
455 :
456 : // Helicity matrix element for a tau decay into three mesons with kaons.
457 :
458 0 : class HMETau2ThreeMesonsWithKaons : public HMETau2ThreeMesons {
459 :
460 : private:
461 :
462 : void initResonances();
463 :
464 : // Resonance masses, widths, and weights.
465 : vector<double> rhoMa, rhoGa, rhoWa, rhoMv, rhoGv, rhoWv;
466 : vector<double> kstarMa, kstarGa, kstarWa, kstarMv, kstarGv, kstarWv;
467 : vector<double> k1Ma, k1Ga, k1Wa, k1Mb, k1Gb, k1Wb;
468 : vector<double> omegaM, omegaG, omegaW;
469 : double kM, piM, piW;
470 :
471 : // Form factors.
472 : complex F1();
473 : complex F2();
474 : complex F4();
475 :
476 : };
477 :
478 : //==========================================================================
479 :
480 : // Helicity matrix element for a tau decay into generic three mesons.
481 :
482 0 : class HMETau2ThreeMesonsGeneric : public HMETau2ThreeMesons {
483 :
484 : private:
485 :
486 : void initResonances();
487 :
488 : // Resonance masses, widths, and weights.
489 : vector<double> rhoMa, rhoGa, rhoWa, rhoMv, rhoGv, rhoWv;
490 : vector<double> kstarM, kstarG, kstarW, k1M, k1G, k1W;
491 : double kM, piM, piW;
492 :
493 : // Form factors.
494 : complex F1();
495 : complex F2();
496 : complex F4();
497 :
498 : };
499 :
500 : //==========================================================================
501 :
502 : // Helicity matrix element for a tau decay into two pions and a photon.
503 :
504 0 : class HMETau2TwoPionsGamma : public HMETauDecay {
505 :
506 : public:
507 :
508 : void initConstants();
509 :
510 : void initWaves(vector<HelicityParticle>&);
511 :
512 : complex calculateME(vector<int>);
513 :
514 : protected:
515 :
516 : // Resonance masses, widths, and weights.
517 : vector<double> rhoM, rhoG, rhoW, omegaM, omegaG, omegaW;
518 : double piM;
519 :
520 : // Form factor.
521 : complex F(double s, vector<double> M, vector<double> G, vector<double> W);
522 :
523 : };
524 :
525 : //==========================================================================
526 :
527 : // Helicity matrix element for a tau decay into four pions.
528 :
529 0 : class HMETau2FourPions : public HMETauDecay {
530 :
531 : public:
532 :
533 : void initConstants();
534 :
535 : void initHadronicCurrent(vector<HelicityParticle>& p);
536 :
537 : private:
538 :
539 : // G-function form factors (fits).
540 : double G(int i, double s);
541 :
542 : // T-vector functions.
543 : Wave4 t1(Wave4&, Wave4&, Wave4&, Wave4&, Wave4&);
544 : Wave4 t2(Wave4&, Wave4&, Wave4&, Wave4&, Wave4&);
545 : Wave4 t3(Wave4&, Wave4&, Wave4&, Wave4&, Wave4&);
546 :
547 : // Breit-Wigner denominators for the intermediate mesons.
548 : complex a1D(double s);
549 : complex rhoD(double s);
550 : complex sigD(double s);
551 : complex omeD(double s);
552 :
553 : // Form factors needed for the a1, rho, and omega.
554 : double a1FormFactor(double s);
555 : double rhoFormFactor1(double s);
556 : double rhoFormFactor2(double s);
557 : double omeFormFactor(double s);
558 :
559 : // Masses and widths of the intermediate mesons.
560 : double a1M, a1G, rhoM, rhoG, sigM, sigG, omeM, omeG;
561 :
562 : // Masses for the pions (charged and neutral).
563 : double picM, pinM;
564 :
565 : // Amplitude, phases, and weights for mixing.
566 : double sigA, sigP, omeA, omeP;
567 : complex sigW, omeW;
568 :
569 : // Cut-off for a1 form factor.
570 : double lambda2;
571 :
572 : };
573 :
574 : //==========================================================================
575 :
576 : // Helicity matrix element for a tau decaying into five pions.
577 :
578 0 : class HMETau2FivePions : public HMETauDecay {
579 :
580 : public:
581 :
582 : void initConstants();
583 :
584 : void initHadronicCurrent(vector<HelicityParticle>&);
585 :
586 : private:
587 :
588 : // Hadronic currents.
589 : Wave4 Ja(Wave4 &q, Wave4 &q1, Wave4 &q2, Wave4 &q3, Wave4 &q4, Wave4 &q5);
590 : Wave4 Jb(Wave4 &q, Wave4 &q1, Wave4 &q2, Wave4 &q3, Wave4 &q4, Wave4 &q5);
591 :
592 : // Simplified s-wave Breit-Wigner assuming massless products.
593 : complex breitWigner(double s, double M, double G);
594 :
595 : // Masses and widths of intermediates.
596 : double a1M, a1G, rhoM, rhoG, omegaM, omegaG, omegaW, sigmaM, sigmaG, sigmaW;
597 :
598 : };
599 :
600 : //==========================================================================
601 :
602 : // Helicity matrix element for a tau decay into flat phase space.
603 :
604 0 : class HMETau2PhaseSpace : public HMETauDecay {
605 :
606 : public:
607 :
608 0 : void initWaves(vector<HelicityParticle>&) {};
609 :
610 0 : complex calculateME(vector<int>) {return 1;}
611 :
612 0 : void calculateD(vector<HelicityParticle>&) {};
613 :
614 0 : void calculateRho(unsigned int, vector<HelicityParticle>&) {};
615 :
616 0 : double decayWeight(vector<HelicityParticle>&) {return 1.0;}
617 :
618 0 : double decayWeightMax(vector<HelicityParticle>&) {return 1.0;}
619 :
620 : };
621 :
622 : //==========================================================================
623 :
624 : } // end namespace Pythia8
625 :
626 : #endif // end Pythia8_HelicityMatrixElements_H
|
