StRoot  1
SusyWidthFunctions.cc
1 // SusyWidthFunctions.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2020 Torbjorn Sjostrand
3 // Authors: N. Desai, P. Skands
4 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
5 // Please respect the MCnet Guidelines, see GUIDELINES for details.
6
8 // SUSY Resonance three-body decay width classes.
9
10 #include "Pythia8/SusyResonanceWidths.h"
11 #include "Pythia8/SusyWidthFunctions.h"
12 #include "Pythia8/SusyCouplings.h"
13 #include "Pythia8/PythiaComplex.h"
14
15 namespace Pythia8 {
16
17 //==========================================================================
18
19 // The WidthFunction class.
20 // Functions to be integrated for calculating the 3-body decay widths.
21
22 //--------------------------------------------------------------------------
23
24 void WidthFunction::setPointers(Info* infoPtrIn) {
25
26  infoPtr = infoPtrIn;
27  particleDataPtr = infoPtrIn->particleDataPtr;
28  coupSUSYPtr = infoPtrIn->coupSUSYPtr;
29  coupSMPtr = infoPtr->coupSMPtr;
30 }
31
32 //==========================================================================
33
34 // The StauWidths class.
35 // Width functions for 3-body stau decays.
36
37 //--------------------------------------------------------------------------
38
39 double StauWidths::getWidth(int idResIn, int idIn){
40
41  setChannel(idResIn, idIn);
42  if (idResIn % 2 == 0) return 0.0;
43
44  // Set up a function of only x, capturing the reference to 'this'
45  auto integrand = [&](double x) { return this->f(x); };
46
47  // Integrate f
48  double width;
49  if (integrateGauss(width, integrand, 0.0, 1.0, 1.0e-3)) return width;
50  else return 0.0;
51
52 }
53
54 //--------------------------------------------------------------------------
55
56 void StauWidths::setChannel(int idResIn, int idIn) {
57
58  // Common masses.
59  idRes = abs(idResIn);
60  idIn = abs(idIn);
61  mRes = particleDataPtr->m0(idRes);
62  m1 = particleDataPtr->m0(1000022);
63  m2 = particleDataPtr->m0(idIn);
64
65  mInt = particleDataPtr->m0(15);
66  gammaInt = particleDataPtr->mWidth(15);
67
68  // Couplings etc.
69  f0 = 92.4; //MeV
70  gf = coupSMPtr->GF();
71  delm = mRes - m1;
72  cons = pow2(f0)*pow2(gf)*(pow2(delm) - pow2(m2))
73  * coupSMPtr->V2CKMid(1,1) / (128.0 * pow(mRes*M_PI,3));
74
75  if (idIn == 9000211) wparam = 1.16;
76  else if (idIn == 213) wparam = 0.808;
77  else wparam = 1.0;
78
79  double g = coupSMPtr->alphaEM(mRes * mRes);
80  int ksusy = 1000000;
81  int isl = (abs(idRes)/ksusy == 2) ? (abs(idRes)%10+1)/2 + 3
82  : (abs(idRes)%10+1)/2;
83
84  gL = g * coupSUSYPtr->LsllX[isl][3][1] / ( sqrt(2.0) * coupSUSYPtr->cosb);
85  gR = g * coupSUSYPtr->RsllX[isl][3][1] / ( sqrt(2.0) * coupSUSYPtr->cosb);
86
87  // Set function switch and internal propagators depending on decay product.
88  if (idIn == 211) fnSwitch = 1;
89  else if (idIn == 9000211 || idIn == 213) fnSwitch = 2;
90  else if (idIn == 14 || idIn == 12) {
91  m2 = particleDataPtr->m0(idIn-1);
92  fnSwitch = 3;
93  }
94  else {
95  stringstream mess;
96  mess << " unknown decay channel idIn = " << idIn;
97  infoPtr->errorMsg("Warning in StauWidths::setChannel:", mess.str() );
98  }
99
100  return;
101 }
102
103 //------------------------------------------------------------------------
104
105 double StauWidths::f(double x){
106
107  // Decay width functions documented in arXiv:1212.2886 Citron et. al.
108  double value = 0.0;
109  double qf2 = pow2(delm) - (pow2(delm) - pow2(m2)) * x;
110  double fac = 1.0 / pow3(mRes);
111  double term3 = (norm(gL) * qf2 + norm(gR) * mInt * mInt)
112  * (delm * delm + 2.0 * m1 * delm - qf2);
113  double term4 = -2.0 * real(gL * conj(gR)) * m2 * mInt * qf2;
114
115  // ~tau -> pi0 nu_tau ~chi0.
116  if (fnSwitch == 1 ) {
117  fac *= pow2(delm) - pow2(m2);
118  double term1 = sqrt((pow2(delm) - qf2) * (pow2(delm + 2 * m1) - qf2));
119  double term2 = pow2(qf2 - pow2(m2)) / qf2 / (pow2(qf2 - pow2(mInt))
120  + pow2(mInt*gammaInt));
121  value = fac * (term1 * term2 * (term3 + term4));
122  }
123
124  else if (fnSwitch == 2 ) {
125  double term1 = sqrt((pow2(delm) - qf2) * (pow2(delm + 2 * m1) - qf2));
126  double term2 = pow2(qf2 - pow2(m2)) * (qf2 + pow2(m2))
127  / (qf2 * qf2 * (pow2(qf2 - pow2(mInt)) + pow2(mInt*gammaInt)));
128  value = fac * (term1 * term2 * (term3 + term4));
129  }
130
131  else if (fnSwitch == 3 ) {
132  double qf4 = qf2 * qf2;
133  double m24 = pow2(m2*m2);
134  double term1 = sqrt((pow2(delm) - qf2) * (pow2(delm + 2 * m1) - qf2));
135  double term2a = 1.0 / (pow2(qf2 - pow2(mInt)) + pow2(mInt*gammaInt)) / qf4;
136  double term2b = 12 * m24 * qf4 * log(qf2/pow2(m2))
137  + (qf4 - m24) * (qf4 - 8 * m2 * m2 * qf2 + m24);
138  value = fac * (term1 * term2a * term2b * (term3 + term4));
139  }
140
141  else {
142  stringstream mess;
143  mess << " unknown decay channel fnSwitch = " << fnSwitch;
144  infoPtr->errorMsg("Warning in StauWidths::function:", mess.str() );
145  }
146
147  return value;
148
149 }
150
151 //==========================================================================
152
153 } //end namespace Pythia8