StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
ResonanceWidthsDM.cc
1 // ResonanceWidthsDM.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2018 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 // Function definitions (not found in the header) for DM resonance properties
7 // in the ResonanceS and ResonanceZp classes.
8 
9 #include "Pythia8/ParticleData.h"
10 #include "Pythia8/ResonanceWidthsDM.h"
11 #include "Pythia8/PythiaComplex.h"
12 
13 namespace Pythia8 {
14 
15 //==========================================================================
16 
17 // The ResonanceS class.
18 // Derived class for S properties. (DMmed(s=0), PDG id 54)
19 
20 //--------------------------------------------------------------------------
21 
22 // Initialize constants.
23 
24 void ResonanceS::initConstants() {
25 
26  // Locally stored properties and couplings.
27  double vq = settingsPtr->parm("Sdm:vf");
28  double vX = settingsPtr->parm("Sdm:vX");
29  double aq = settingsPtr->parm("Sdm:af");
30  double aX = settingsPtr->parm("Sdm:aX");
31 
32  gq = abs(aq) > 0 ? aq : vq;
33  gX = abs(aX) > 0 ? aX : vX;
34 
35  if (abs(aX) > 0) pScalar = true;
36  else pScalar = false;
37 
38 }
39 
40 //--------------------------------------------------------------------------
41 
42 // Calculate various common prefactors for the current mass.
43 
44 void ResonanceS::calcPreFac(bool) {
45 
46  // Common coupling factors.
47  preFac = 1.0 / (12.0 * M_PI * mRes);
48  alpS = couplingsPtr->alphaS(mHat * mHat);
49 
50 }
51 
52 //--------------------------------------------------------------------------
53 
54 // Calculate width for currently considered channel.
55 
56 void ResonanceS::calcWidth(bool) {
57 
58  // Check that above threshold.
59  if (ps == 0.) return;
60 
61  double mRat2 = pow2(mf1 / mRes);
62  double kinfac = (1 - 4 * mRat2) * (1. + 2 * mRat2);
63 
64  widNow = 0.;
65 
66  if (id1Abs == 21)
67  widNow = pow2(gq) * preFac * pow2(alpS / M_PI) * eta2gg();
68 
69  if(id1Abs < 7)
70  widNow = 3. * pow2(gq * mf1) * preFac * kinfac;
71 
72  if(id1Abs == 52)
73  widNow = pow2(gX * mf1) * preFac * kinfac;
74 
75 }
76 
77 //--------------------------------------------------------------------------
78 
79 double ResonanceS::eta2gg() {
80 
81  // Initial values.
82  complex eta = complex(0., 0.);
83  double mLoop, epsilon, root, rootLog;
84  complex phi, etaNow;
85 
86  // Loop over s, c, b, t quark flavours.
87  for (int idNow = 3; idNow < 7; ++idNow) {
88  mLoop = particleDataPtr->m0(idNow);
89  epsilon = pow2(2. * mLoop / mHat);
90 
91  // Value of loop integral.
92  if (epsilon <= 1.) {
93  root = sqrt(1. - epsilon);
94  rootLog = (epsilon < 1e-4) ? log(4. / epsilon - 2.)
95  : log( (1. + root) / (1. - root) );
96  phi = complex( -0.25 * (pow2(rootLog) - pow2(M_PI)),
97  0.5 * M_PI * rootLog );
98  }
99  else phi = complex( pow2( asin(1. / sqrt(epsilon)) ), 0.);
100 
101  // Factors that depend on Higgs and flavour type.
102  if (!pScalar) etaNow = -0.5 * epsilon
103  * (complex(1., 0.) + (1. - epsilon) * phi);
104  else etaNow = -0.5 * epsilon * phi;
105 
106 
107  // Sum up contribution and return square of absolute value.
108  eta += etaNow;
109  }
110  return (pow2(eta.real()) + pow2(eta.imag()));
111 
112 }
113 
114 //==========================================================================
115 
116 // The ResonanceZp class.
117 // Derived class for Z'^0 properties. (DMmed(s=1), PDG id 55)
118 
119 //--------------------------------------------------------------------------
120 
121 // Initialize constants.
122 
123 void ResonanceZp::initConstants() {
124 
125  // Locally stored properties and couplings.
126  gZp = settingsPtr->parm("Zp:gZp");
127  vu = settingsPtr->parm("Zp:vu");
128  vd = settingsPtr->parm("Zp:vd");
129  vl = settingsPtr->parm("Zp:vl");
130  vv = settingsPtr->parm("Zp:vv");
131  vX = settingsPtr->parm("Zp:vX");
132  au = settingsPtr->parm("Zp:au");
133  ad = settingsPtr->parm("Zp:ad");
134  al = settingsPtr->parm("Zp:al");
135  av = settingsPtr->parm("Zp:av");
136  aX = settingsPtr->parm("Zp:aX");
137 
138 }
139 
140 //--------------------------------------------------------------------------
141 
142 // Calculate various common prefactors for the current mass.
143 
144 void ResonanceZp::calcPreFac(bool) {
145 
146  // Common coupling factors.
147  preFac = pow2(gZp) * mRes / 12.0 / M_PI;
148 
149 }
150 
151 //--------------------------------------------------------------------------
152 
153 // Calculate width for currently considered channel.
154 
155 void ResonanceZp::calcWidth(bool) {
156 
157  // Check that above threshold.
158  if (ps == 0. || id1 * id2 > 0) return;
159 
160  double kinFacA = pow3(ps);
161  double kinFacV = ps * (1. + 2. * mr1);
162 
163  double fac = 0;
164  widNow = 0.;
165 
166  if (id1Abs < 7) {
167  if (id1Abs %2 ) fac = vu * vu * kinFacV + au * au * kinFacA;
168  else fac = vd * vd * kinFacV + ad * ad * kinFacA;
169  }
170  else if (id1Abs > 7 && id1Abs < 17){
171  if (id1Abs %2 ) fac = vl * vl * kinFacV + al * al * kinFacA;
172  else fac = vv * vv * kinFacV + av * av * kinFacA;
173  }
174 
175  if (id1Abs == 52) fac = vX * vX * kinFacV + aX * aX * kinFacA;
176  widNow = fac * preFac;
177  if (id1Abs < 7) widNow *= 3;
178 
179 }
180 
181 //==========================================================================
182 
183 } // end namespace Pythia8