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) 2020 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, ResonanceZp, ResonanceS1, ResonanceCha, ResonanceDM2,
8 // and ResonanceChaD classes.
9 
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 = coupSMPtr->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  // Caclulate current width.
62  double mRat2 = pow2(mf1 / mRes);
63  double kinfac = (1 - 4 * mRat2) * (1. + 2 * mRat2);
64  widNow = 0.;
65  if (id1Abs < 7) widNow = 3. * pow2(gq * mf1) * preFac * kinfac;
66  if (id1Abs == 21) widNow = pow2(gq) * preFac * pow2(alpS / M_PI) * eta2gg();
67  if (id1Abs == 52) widNow = pow2(gX * mf1) * preFac * kinfac;
68 
69 }
70 
71 //--------------------------------------------------------------------------
72 
73 // Loop integral for H -> gg coupling.
74 
75 double ResonanceS::eta2gg() {
76 
77  // Initial values.
78  complex eta = complex(0., 0.);
79  double mLoop, epsilon, root, rootLog;
80  complex phi, etaNow;
81 
82  // Loop over s, c, b, t quark flavours.
83  for (int idNow = 3; idNow < 7; ++idNow) {
84  mLoop = particleDataPtr->m0(idNow);
85  epsilon = pow2(2. * mLoop / mHat);
86 
87  // Value of loop integral.
88  if (epsilon <= 1.) {
89  root = sqrt(1. - epsilon);
90  rootLog = (epsilon < 1e-4) ? log(4. / epsilon - 2.)
91  : log( (1. + root) / (1. - root) );
92  phi = complex( -0.25 * (pow2(rootLog) - pow2(M_PI)),
93  0.5 * M_PI * rootLog );
94  }
95  else phi = complex( pow2( asin(1. / sqrt(epsilon)) ), 0.);
96 
97  // Factors that depend on Higgs and flavour type.
98  if (!pScalar) etaNow = -0.5 * epsilon
99  * (complex(1., 0.) + (1. - epsilon) * phi);
100  else etaNow = -0.5 * epsilon * phi;
101 
102 
103  // Sum up contribution and return square of absolute value.
104  eta += etaNow;
105  }
106  return (pow2(eta.real()) + pow2(eta.imag()));
107 
108 }
109 
110 //==========================================================================
111 
112 // The ResonanceZp class.
113 // Derived class for Z'^0 properties. (DMmed(s=1), PDG id 55.)
114 
115 //--------------------------------------------------------------------------
116 
117 // Initialize constants.
118 
119 void ResonanceZp::initConstants() {
120 
121  // Coupling type and common couplings.
122  kinMix = settingsPtr->flag("Zp:kineticMixing");
123  gZp = settingsPtr->parm("Zp:gZp");
124  eps = settingsPtr->parm("Zp:epsilon");
125  vX = settingsPtr->parm("Zp:vX");
126  aX = settingsPtr->parm("Zp:aX");
127 
128  // SM fermion couplings for kinetic mixing case.
129  if (kinMix) {
130  vu = eps * (2./3. + coupSMPtr->vf(2));
131  au = eps * coupSMPtr->af(2);
132  vd = eps * (-1./3. + coupSMPtr->vf(1));
133  ad = eps * coupSMPtr->af(1);
134  vl = eps * (-1. + coupSMPtr->vf(11));
135  al = eps * coupSMPtr->af(11);
136  vv = eps * coupSMPtr->vf(12);
137  av = eps * coupSMPtr->af(12);
138 
139  // SM fermion couplings set by user.
140  } else {
141  vu = settingsPtr->parm("Zp:vu");
142  vd = settingsPtr->parm("Zp:vd");
143  vl = settingsPtr->parm("Zp:vl");
144  vv = settingsPtr->parm("Zp:vv");
145  au = settingsPtr->parm("Zp:au");
146  ad = settingsPtr->parm("Zp:ad");
147  al = settingsPtr->parm("Zp:al");
148  av = settingsPtr->parm("Zp:av");
149  }
150 
151 }
152 
153 //--------------------------------------------------------------------------
154 
155 // Calculate various common prefactors for the current mass.
156 
157 void ResonanceZp::calcPreFac(bool) {
158 
159  // Common coupling factors.
160  preFac = mRes / 12.0 / M_PI;
161 
162 }
163 
164 //--------------------------------------------------------------------------
165 
166 // Calculate width for currently considered channel.
167 
168 void ResonanceZp::calcWidth(bool) {
169 
170  // Check that above threshold.
171  if (ps == 0. || id1 * id2 > 0) return;
172 
173  // Couplings to Zp.
174  double kinFacA = pow3(ps);
175  double kinFacV = ps * (1. + 2. * mr1);
176  double fac = 0;
177  widNow = 0.;
178  if (id1Abs < 7) {
179  if (id1Abs%2 == 0) fac = vu * vu * kinFacV + au * au * kinFacA;
180  else fac = vd * vd * kinFacV + ad * ad * kinFacA;
181  widNow *= 3;
182  } else if (id1Abs > 10 && id1Abs < 17) {
183  if (id1Abs%2 == 1) fac = vl * vl * kinFacV + al * al * kinFacA;
184  else fac = vv * vv * kinFacV + av * av * kinFacA;
185  } else if (id1Abs == 52) fac = vX * vX * kinFacV + aX * aX * kinFacA;
186 
187  // Set partial width.
188  double coup = pow2(gZp);
189  if (kinMix && id1Abs != 52)
190  coup = coupSMPtr->alphaEM(mRes * mRes) * 4.0 * M_PI;
191  widNow = coup * fac * preFac;
192 
193 }
194 
195 //==========================================================================
196 
197 // The ResonanceSl class.
198 // Derived class for Sl properties. (Using PDG id 56.)
199 
200 //--------------------------------------------------------------------------
201 
202 // Initialize constants.
203 
204 void ResonanceSl::initConstants() {
205 
206  // Locally stored properties and couplings.
207  yuk[0] = 0.0;
208  yuk[1] = settingsPtr->parm("DM:yuk1");
209  yuk[2] = settingsPtr->parm("DM:yuk2");
210  yuk[3] = settingsPtr->parm("DM:yuk3");
211 
212 }
213 
214 //--------------------------------------------------------------------------
215 
216 // Calculate various common prefactors for the current mass.
217 
218 void ResonanceSl::calcPreFac(bool) {
219 
220  // Common coupling factors.
221  preFac = 1. / (mRes * 16.0 * M_PI);
222 
223 }
224 
225 //--------------------------------------------------------------------------
226 
227 // Calculate width for currently considered channel.
228 
229 void ResonanceSl::calcWidth(bool) {
230 
231  // Check that above threshold.
232  if (ps == 0.) return;
233 
234  // Set partial width.
235  kinFac = (mRes * mRes - mf1 * mf1 - mf2 * mf2);
236  double coup = 0.0;
237  if(abs(id2) == 11) coup = yuk[1];
238  if(abs(id2) == 13) coup = yuk[2];
239  if(abs(id2) == 15) coup = yuk[3];
240  widNow = pow2(coup) * preFac * kinFac * ps;
241 
242  return;
243 
244 }
245 
246 //==========================================================================
247 
248 // The ResonanceCha class.
249 // Derived class for singly-charged partner properties. (Using PDG id 57.)
250 
251 //--------------------------------------------------------------------------
252 
253 void ResonanceCha::setMassMix(){
254 
255  // Check if using the Nplet model
256  // (do not over-ride DM mass if using s-channel mediators).
257  doDY = settingsPtr->flag("DM:qqbar2DY") &&
258  settingsPtr->mode("DM:DYtype") > 1;
259  if (!doDY) return;
260 
261  // Locally stored properties and couplings.
262  double M1 = settingsPtr->parm("DM:M1");
263  double M2 = settingsPtr->parm("DM:M2");
264  int type = settingsPtr->mode("DM:Nplet");
265  double Lambda = settingsPtr->parm("DM:Lambda");
266 
267  // Mixing parameters.
268  double vev = 174.0;
269  mixing = vev / Lambda;
270  if (type > 1) mixing *= sqrt(2) * vev;
271  if (type > 2) mixing *= pow2(vev) /pow2(Lambda) / sqrt(12);
272  double term1 = sqrt(pow2(M2 - M1) + pow2(mixing));
273  double sin2th = 0.5 * (1 - abs(M2 - M1) / term1);
274  mixN1 = (M1 > M2) ? sqrt(sin2th) : sqrt(1. - sin2th);
275  mixN2 = (M1 > M2) ? sqrt(1. - sin2th) : sqrt(sin2th);
276 
277  // Set particle masses.
278  double m1m = 0.5 * (M1 + M2 - term1);
279  double m2p = 0.5 * (M1 + M2 + term1);
280  double mplet = (M1 < M2) ? m2p : m1m;
281  // Neutral partners.
282  particleDataPtr->m0(52, m1m);
283  particleDataPtr->m0(58, m2p);
284  // Chargino mass with NLO mass splitting 0.16 GeV. Doubly charged mass.
285  particleDataPtr->m0(57, mplet + 0.16);
286  particleDataPtr->m0(59, mplet + 0.16 + 0.49);
287 
288  return;
289 
290 }
291 
292 //--------------------------------------------------------------------------
293 
294 // Calculate various common prefactors for the current mass.
295 
296 void ResonanceCha::calcPreFac(bool) {
297 
298  // Common coupling factors.
299  preFac = mRes / (16.0 * M_PI);
300 
301 }
302 
303 //--------------------------------------------------------------------------
304 
305 // Calculate width for currently considered channel.
306 
307 void ResonanceCha::calcWidth(bool) {
308 
309  // Check if using the Nplet model.
310  if (!doDY) return;
311 
312  // Small mass difference normal for compressed spectrum, so smaller
313  // MASSMARGIN required. Check that above threshold.
314  if (mHat < mf1 + mf2 + 0.01) return;
315 
316  // Mass parameters.
317  double dm = 0.0;
318  widNow = 0.0;
319  double mix = (abs(id2) == 58) ? mixN2 : mixN1;
320  double mpion = 0.1396;
321 
322  // Two-body decays.
323  if (mult == 2) {
324  dm = particleDataPtr->m0(57) - particleDataPtr->m0(abs(id2));
325  // Long-lived two-body decay via pion.
326  if (dm > mpion) widNow = 2.0 * pow2(mix) * 6.993e-13
327  * sqrt(1. - pow2(mpion/dm)) * pow3(dm);
328  // Two-body decay via W.
329  else if (dm > particleDataPtr->m0(24)) ;
330  else return;
331 
332  // Three-body decays.
333  } else { }
334 
335  return;
336 
337 }
338 
339 //==========================================================================
340 
341 // The ResonanceDM2 class.
342 // Derived class for neutral partner properties. (Using PDG id 58.)
343 // Not yet implemented.
344 
345 //--------------------------------------------------------------------------
346 
347 // Initialize constants.
348 
349 void ResonanceDM2::initConstants() {
350 
351  setMassMix();
352  mHiggs = particleDataPtr->m0(25);
353  wHiggs = particleDataPtr->mWidth(25);
354 
355 }
356 
357 //--------------------------------------------------------------------------
358 
359 // Calculate various common prefactors for the current mass.
360 
361 void ResonanceDM2::calcPreFac(bool) {
362 
363 }
364 
365 //--------------------------------------------------------------------------
366 
367 // Calculate width for currently considered channel.
368 
369 void ResonanceDM2::calcWidth(bool) {
370 
371  // Check if using the Nplet model
372  if (!doDY) return;
373 
374  // Check that above threshold.
375  if (ps == 0.) return;
376 
377  return;
378 
379 }
380 
381 //==========================================================================
382 
383 // The ResonanceChaD class.
384 // Derived class for doubly charged properties. (Using PDG id 59.)
385 
386 //--------------------------------------------------------------------------
387 
388 // Calculate various common prefactors for the current mass.
389 
390 void ResonanceChaD::calcPreFac(bool) {
391 
392  // Common coupling factors.
393  double dm = particleDataPtr->m0(59) - particleDataPtr->m0(57);
394  preFac = (dm > 0.) ? 4.0 * 6.993e-13 * sqrtpos(1. - pow2(0.1396/dm))
395  * pow3(dm) : 0.;
396 
397 }
398 
399 //--------------------------------------------------------------------------
400 
401 // Calculate width for currently considered channel.
402 
403 void ResonanceChaD::calcWidth(bool) {
404 
405  // Check if using the Nplet model.
406  if (!doDY) return;
407 
408  widNow = preFac;
409  return;
410 
411 }
412 
413 //==========================================================================
414 
415 } // end namespace Pythia8
Generated by   doxygen 1.8.5