PYTHIA8/pythia8175/include/HelicityMatrixElements.h

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