StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
ResonanceWidths.cc
1 // ResonanceWidths.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2014 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 // Function definitions (not found in the header) for
7 // the ResonanceWidths class and classes derived from it.
8 
9 #include "Pythia8/ParticleData.h"
10 #include "Pythia8/ResonanceWidths.h"
11 #include "Pythia8/PythiaComplex.h"
12 
13 namespace Pythia8 {
14 
15 //==========================================================================
16 
17 // The ResonanceWidths class.
18 // Base class for the various resonances.
19 
20 //--------------------------------------------------------------------------
21 
22 // Constants: could be changed here if desired, but normally should not.
23 // These are of technical nature, as described for each.
24 
25 // Number of points in integration direction for numInt routines.
26 const int ResonanceWidths::NPOINT = 100;
27 
28 // The mass of a resonance must not be too small.
29 const double ResonanceWidths::MASSMIN = 0.4;
30 
31 // The sum of product masses must not be too close to the resonance mass.
32 const double ResonanceWidths::MASSMARGIN = 0.1;
33 
34 //--------------------------------------------------------------------------
35 
36 // Initialize data members.
37 // Calculate and store partial and total widths at the nominal mass.
38 
39 bool ResonanceWidths::init(Info* infoPtrIn, Settings* settingsPtrIn,
40  ParticleData* particleDataPtrIn, Couplings* couplingsPtrIn) {
41 
42  // Save pointers.
43  infoPtr = infoPtrIn;
44  settingsPtr = settingsPtrIn;
45  particleDataPtr = particleDataPtrIn;
46  couplingsPtr = couplingsPtrIn;
47 
48  // Perform any model dependent initialisations (pure dummy in base class).
49  bool isInit = initBSM();
50 
51  // Minimal decaying-resonance width. Minimal phase space for meMode = 103.
52  minWidth = settingsPtr->parm("ResonanceWidths:minWidth");
53  minThreshold = settingsPtr->parm("ResonanceWidths:minThreshold");
54 
55  // Pointer to particle species.
56  particlePtr = particleDataPtr->particleDataEntryPtr(idRes);
57  if (particlePtr == 0) {
58  infoPtr->errorMsg("Error in ResonanceWidths::init:"
59  " unknown resonance identity code");
60  return false;
61  }
62 
63  // Generic particles should not have meMode < 100, but allow
64  // some exceptions: not used Higgses and not used Technicolor.
65  if (idRes == 35 || idRes == 36 || idRes == 37
66  || idRes/1000000 == 3) isGeneric = false;
67 
68  // Resonance properties: antiparticle, mass, width
69  hasAntiRes = particlePtr->hasAnti();
70  mRes = particlePtr->m0();
71  GammaRes = particlePtr->mWidth();
72  m2Res = mRes*mRes;
73 
74  // A resonance cannot be too light, in particular not = 0.
75  if (mRes < MASSMIN) {
76  ostringstream idCode;
77  idCode << idRes;
78  infoPtr->errorMsg("Error in ResonanceWidths::init:"
79  " resonance mass too small", "for id = " + idCode.str(), true);
80  return false;
81  }
82 
83  // For very narrow resonances assign fictitious small width.
84  if (GammaRes < minWidth) GammaRes = 0.1 * minWidth;
85  GamMRat = (mRes == 0.) ? 0. : GammaRes / mRes;
86 
87  // Secondary decay chains by default all on.
88  openPos = 1.;
89  openNeg = 1.;
90 
91  // Allow option where on-shell width is forced to current value.
92  // Disable for mixes gamma*/Z0/Z'0
93  doForceWidth = particlePtr->doForceWidth();
94  if (idRes == 23 && settingsPtr->mode("WeakZ0:gmZmode") != 2)
95  doForceWidth = false;
96  if (idRes == 33 && settingsPtr->mode("Zprime:gmZmode") != 3)
97  doForceWidth = false;
98  forceFactor = 1.;
99 
100  // Check if we are supposed to do the width calculation
101  // (can be false e.g. if SLHA decay table should take precedence instead).
102  allowCalcWidth = isInit && allowCalc();
103  if ( allowCalcWidth ) {
104  // Initialize constants used for a resonance.
105  initConstants();
106 
107  // Calculate various common prefactors for the current mass.
108  mHat = mRes;
109  calcPreFac(true);
110  }
111 
112  // Reset quantities to sum. Declare variables inside loop.
113  double widTot = 0.;
114  double widPos = 0.;
115  double widNeg = 0.;
116  int idNow, idAnti;
117  double openSecPos, openSecNeg;
118 
119  // Loop over all decay channels. Basic properties of channel.
120  for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
121  iChannel = i;
122  onMode = particlePtr->channel(i).onMode();
123  meMode = particlePtr->channel(i).meMode();
124  mult = particlePtr->channel(i).multiplicity();
125  widNow = 0.;
126 
127  // Warn if not relevant meMode.
128  if ( meMode < 0 || meMode > 103 || (isGeneric && meMode < 100) ) {
129  infoPtr->errorMsg("Error in ResonanceWidths::init:"
130  " resonance meMode not acceptable");
131  }
132 
133  // Channels with meMode < 100 must be implemented in derived classes.
134  if (meMode < 100 && allowCalcWidth) {
135 
136  // Read out information on channel: primarily use first two.
137  id1 = particlePtr->channel(i).product(0);
138  id2 = particlePtr->channel(i).product(1);
139  id1Abs = abs(id1);
140  id2Abs = abs(id2);
141 
142  // Order first two in descending order of absolute values.
143  if (id2Abs > id1Abs) {swap( id1, id2); swap( id1Abs, id2Abs);}
144 
145  // Allow for third product to be treated in derived classes.
146  if (mult > 2) {
147  id3 = particlePtr->channel(i).product(2);
148  id3Abs = abs(id3);
149 
150  // Also order third into descending order of absolute values.
151  if (id3Abs > id2Abs) {swap( id2, id3); swap( id2Abs, id3Abs);}
152  if (id2Abs > id1Abs) {swap( id1, id2); swap( id1Abs, id2Abs);}
153  }
154 
155  // Read out masses. Calculate two-body phase space.
156  mf1 = particleDataPtr->m0(id1Abs);
157  mf2 = particleDataPtr->m0(id2Abs);
158  mr1 = pow2(mf1 / mHat);
159  mr2 = pow2(mf2 / mHat);
160  ps = (mHat < mf1 + mf2 + MASSMARGIN) ? 0.
161  : sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2 );
162  if (mult > 2) {
163  mf3 = particleDataPtr->m0(id3Abs);
164  mr3 = pow2(mf3 / mHat);
165  }
166 
167  // Let derived class calculate width for channel provided.
168  calcWidth(true);
169  }
170 
171  // Channels with meMode >= 100 are calculated based on stored values.
172  else widNow = GammaRes * particlePtr->channel(i).bRatio();
173 
174  // Find secondary open fractions of partial width.
175  openSecPos = 1.;
176  openSecNeg = 1.;
177  if (widNow > 0.) for (int j = 0; j < mult; ++j) {
178  idNow = particlePtr->channel(i).product(j);
179  idAnti = (particleDataPtr->hasAnti(idNow)) ? -idNow : idNow;
180  // Secondary widths not yet initialized for heavier states,
181  // so have to assume unit open fraction there.
182  if (idNow == 23 || abs(idNow) == 24 || idNow == 93 || abs(idNow) == 94
183  || particleDataPtr->m0(abs(idNow)) < mRes) {
184  openSecPos *= particleDataPtr->resOpenFrac(idNow);
185  openSecNeg *= particleDataPtr->resOpenFrac(idAnti);
186  }
187  }
188 
189  // Store partial widths and secondary open fractions.
190  particlePtr->channel(i).onShellWidth(widNow);
191  particlePtr->channel(i).openSec( idRes, openSecPos);
192  particlePtr->channel(i).openSec(-idRes, openSecNeg);
193 
194  // Update sum over all channnels and over open channels only.
195  widTot += widNow;
196  if (onMode == 1 || onMode == 2) widPos += widNow * openSecPos;
197  if (onMode == 1 || onMode == 3) widNeg += widNow * openSecNeg;
198  }
199 
200  // If no decay channels are open then set particle stable and done.
201  if (widTot < minWidth) {
202  particlePtr->setMayDecay(false, false);
203  particlePtr->setMWidth(0., false);
204  for (int i = 0; i < particlePtr->sizeChannels(); ++i)
205  particlePtr->channel(i).bRatio( 0., false);
206  return true;
207  }
208 
209  // Normalize branching ratios to unity.
210  double bRatio;
211  for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
212  bRatio = particlePtr->channel(i).onShellWidth() / widTot;
213  particlePtr->channel(i).bRatio( bRatio, false);
214  }
215 
216  // Optionally force total width by rescaling of all partial ones.
217  if (doForceWidth) {
218  forceFactor = GammaRes / widTot;
219  for (int i = 0; i < particlePtr->sizeChannels(); ++i)
220  particlePtr->channel(i).onShellWidthFactor( forceFactor);
221  }
222 
223  // Else update newly calculated partial width.
224  else {
225  particlePtr->setMWidth(widTot, false);
226  GammaRes = widTot;
227  }
228 
229  // Updated width-to-mass ratio. Secondary widths for open.
230  GamMRat = GammaRes / mRes;
231  openPos = widPos / widTot;
232  openNeg = widNeg / widTot;
233 
234  // Clip wings of Higgses.
235  bool isHiggs = (idRes == 25 || idRes == 35 ||idRes == 36 ||idRes == 37);
236  bool clipHiggsWings = settingsPtr->flag("Higgs:clipWings");
237  if (isHiggs && clipHiggsWings) {
238  double mMinNow = particlePtr->mMin();
239  double mMaxNow = particlePtr->mMax();
240  double wingsFac = settingsPtr->parm("Higgs:wingsFac");
241  double mMinWing = mRes - wingsFac * GammaRes;
242  double mMaxWing = mRes + wingsFac * GammaRes;
243  if (mMinWing > mMinNow) particlePtr->setMMinNoChange(mMinWing);
244  if (mMaxWing < mMaxNow || mMaxNow < mMinNow)
245  particlePtr->setMMaxNoChange(mMaxWing);
246  }
247 
248  // Done.
249  return true;
250 
251 }
252 
253 //--------------------------------------------------------------------------
254 
255 // Calculate the total width and store phase-space-weighted coupling sums.
256 
257 double ResonanceWidths::width(int idSgn, double mHatIn, int idInFlavIn,
258  bool openOnly, bool setBR, int idOutFlav1, int idOutFlav2) {
259 
260  // Calculate various prefactors for the current mass.
261  mHat = mHatIn;
262  idInFlav = idInFlavIn;
263  if (allowCalcWidth) calcPreFac(false);
264 
265  // Reset quantities to sum. Declare variables inside loop.
266  double widSum = 0.;
267  double mfSum, psOnShell;
268 
269  // Loop over all decay channels. Basic properties of channel.
270  for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
271  iChannel = i;
272  onMode = particlePtr->channel(i).onMode();
273  meMode = particlePtr->channel(i).meMode();
274  mult = particlePtr->channel(i).multiplicity();
275 
276  // Initially assume vanishing branching ratio.
277  widNow = 0.;
278  if (setBR) particlePtr->channel(i).currentBR(widNow);
279 
280  // Optionally only consider specific (two-body) decay channel.
281  // Currently only used for Higgs -> q qbar, g g or gamma gamma.
282  if (idOutFlav1 > 0 || idOutFlav2 > 0) {
283  if (mult > 2) continue;
284  if (particlePtr->channel(i).product(0) != idOutFlav1) continue;
285  if (particlePtr->channel(i).product(1) != idOutFlav2) continue;
286  }
287 
288  // Optionally only consider open channels.
289  if (openOnly) {
290  if (idSgn > 0 && onMode !=1 && onMode != 2) continue;
291  if (idSgn < 0 && onMode !=1 && onMode != 3) continue;
292  }
293 
294  // Channels with meMode < 100 must be implemented in derived classes.
295  if (meMode < 100) {
296 
297  // Read out information on channel: primarily use first two.
298  id1 = particlePtr->channel(i).product(0);
299  id2 = particlePtr->channel(i).product(1);
300  id1Abs = abs(id1);
301  id2Abs = abs(id2);
302 
303  // Order first two in descending order of absolute values.
304  if (id2Abs > id1Abs) {swap( id1, id2); swap( id1Abs, id2Abs);}
305 
306  // Allow for third product to be treated in derived classes.
307  if (mult > 2) {
308  id3 = particlePtr->channel(i).product(2);
309  id3Abs = abs(id3);
310 
311  // Also order third into descending order of absolute values.
312  if (id3Abs > id2Abs) {swap( id2, id3); swap( id2Abs, id3Abs);}
313  if (id2Abs > id1Abs) {swap( id1, id2); swap( id1Abs, id2Abs);}
314  }
315 
316  // Read out masses. Calculate two-body phase space.
317  mf1 = particleDataPtr->m0(id1Abs);
318  mf2 = particleDataPtr->m0(id2Abs);
319  mr1 = pow2(mf1 / mHat);
320  mr2 = pow2(mf2 / mHat);
321  ps = (mHat < mf1 + mf2 + MASSMARGIN) ? 0.
322  : sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2 );
323  if (mult > 2) {
324  mf3 = particleDataPtr->m0(id3Abs);
325  mr3 = pow2(mf3 / mHat);
326  }
327 
328  // Let derived class calculate width for channel provided.
329  calcWidth(false);
330  }
331 
332  // Now on to meMode >= 100. First case: no correction at all.
333  else if (meMode == 100)
334  widNow = GammaRes * particlePtr->channel(i).bRatio();
335 
336  // Correction by step at threshold.
337  else if (meMode == 101) {
338  mfSum = 0.;
339  for (int j = 0; j < mult; ++j) mfSum
340  += particleDataPtr->m0( particlePtr->channel(i).product(j) );
341  if (mfSum + MASSMARGIN < mHat)
342  widNow = GammaRes * particlePtr->channel(i).bRatio();
343  }
344 
345  // Correction by a phase space factor for two-body decays.
346  else if ( (meMode == 102 || meMode == 103) && mult == 2) {
347  mf1 = particleDataPtr->m0( particlePtr->channel(i).product(0) );
348  mf2 = particleDataPtr->m0( particlePtr->channel(i).product(1) );
349  mr1 = pow2(mf1 / mHat);
350  mr2 = pow2(mf2 / mHat);
351  ps = (mHat < mf1 + mf2 + MASSMARGIN) ? 0.
352  : sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2 );
353  mr1 = pow2(mf1 / mRes);
354  mr2 = pow2(mf2 / mRes);
355  psOnShell = (meMode == 102) ? 1. : max( minThreshold,
356  sqrtpos( pow2(1.- mr1 - mr2) - 4. * mr1 * mr2) );
357  widNow = GammaRes * particlePtr->channel(i).bRatio() * ps / psOnShell;
358  }
359 
360  // Correction by simple threshold factor for multibody decay.
361  else if (meMode == 102 || meMode == 103) {
362  mfSum = 0.;
363  for (int j = 0; j < mult; ++j) mfSum
364  += particleDataPtr->m0( particlePtr->channel(i).product(j) );
365  ps = sqrtpos(1. - mfSum / mHat);
366  psOnShell = (meMode == 102) ? 1. : max( minThreshold,
367  sqrtpos(1. - mfSum / mRes) );
368  widNow = GammaRes * particlePtr->channel(i).bRatio() * ps / psOnShell;
369  }
370 
371  // Optionally multiply by secondary widths.
372  if (openOnly) widNow *= particlePtr->channel(i).openSec(idSgn);
373 
374  // Optionally include factor to force to fixed width.
375  if (doForceWidth) widNow *= forceFactor;
376 
377  // Optionally multiply by current/nominal resonance mass??
378 
379  // Sum back up.
380  widSum += widNow;
381 
382  // Optionally store partial widths for later decay channel choice.
383  if (setBR) particlePtr->channel(i).currentBR(widNow);
384  }
385 
386  // Done.
387  return widSum;
388 
389 }
390 
391 //--------------------------------------------------------------------------
392 
393 // Numerical integration of matrix-element in two-body decay,
394 // where one particle is described by a Breit-Wigner mass distribution.
395 // Normalization to unit integral if matrix element is unity
396 // and there are no phase-space restrictions.
397 
398 double ResonanceWidths::numInt1BW(double mHatIn, double m1, double Gamma1,
399  double mMin1, double m2, int psMode) {
400 
401  // Check that phase space is open for integration.
402  if (mMin1 + m2 > mHatIn) return 0.;
403 
404  // Precalculate coefficients for Breit-Wigner selection.
405  double s1 = m1 * m1;
406  double mG1 = m1 * Gamma1;
407  double mMax1 = mHatIn - m2;
408  double atanMin1 = atan( (mMin1 * mMin1 - s1) / mG1 );
409  double atanMax1 = atan( (mMax1 * mMax1 - s1) / mG1 );
410  double atanDif1 = atanMax1 - atanMin1;
411  double wtDif1 = atanDif1 / (M_PI * NPOINT);
412 
413  // Step size in atan-mapped variable.
414  double xStep = 1. / NPOINT;
415 
416  // Variables used in loop over integration points.
417  double sum = 0.;
418  double mrNow2 = pow2(m2 / mHatIn);
419  double xNow1, sNow1, mNow1, mrNow1, psNow, value;
420 
421  // Loop with first-particle mass selection.
422  for (int ip1 = 0; ip1 < NPOINT; ++ip1) {
423  xNow1 = xStep * (ip1 + 0.5);
424  sNow1 = s1 + mG1 * tan(atanMin1 + xNow1 * atanDif1);
425  mNow1 = min( mMax1, max( mMin1, sqrtpos(sNow1) ) );
426  mrNow1 = pow2(mNow1 / mHatIn);
427 
428  // Evaluate value and add to sum. Different matrix elements.
429  psNow = sqrtpos( pow2(1. - mrNow1 - mrNow2)
430  - 4. * mrNow1 * mrNow2);
431  value = 1.;
432  if (psMode == 1) value = psNow;
433  else if (psMode == 2) value = psNow * psNow;
434  else if (psMode == 3) value = pow3(psNow);
435  else if (psMode == 5) value = psNow
436  * (pow2(1. - mrNow1 - mrNow2) + 8. * mrNow1 * mrNow2);
437  else if (psMode == 6) value = pow3(psNow);
438  sum += value;
439 
440  // End of loop over integration points. Overall normalization.
441  }
442  sum *= wtDif1;
443 
444  // Done.
445  return sum;
446 }
447 
448 //--------------------------------------------------------------------------
449 
450 // Numerical integration of matrix-element in two-body decay,
451 // where both particles are described by Breit-Wigner mass distributions.
452 // Normalization to unit integral if matrix element is unity
453 // and there are no phase-space restrictions.
454 
455 double ResonanceWidths::numInt2BW(double mHatIn, double m1, double Gamma1,
456  double mMin1, double m2, double Gamma2, double mMin2, int psMode) {
457 
458  // Check that phase space is open for integration.
459  if (mMin1 + mMin2 >= mHatIn) return 0.;
460 
461  // Precalculate coefficients for Breit-Wigner selection.
462  double s1 = m1 * m1;
463  double mG1 = m1 * Gamma1;
464  double mMax1 = mHatIn - mMin2;
465  double atanMin1 = atan( (mMin1 * mMin1 - s1) / mG1 );
466  double atanMax1 = atan( (mMax1 * mMax1 - s1) / mG1 );
467  double atanDif1 = atanMax1 - atanMin1;
468  double wtDif1 = atanDif1 / (M_PI * NPOINT);
469  double s2 = m2 * m2;
470  double mG2 = m2 * Gamma2;
471  double mMax2 = mHatIn - mMin1;
472  double atanMin2 = atan( (mMin2 * mMin2 - s2) / mG2 );
473  double atanMax2 = atan( (mMax2 * mMax2 - s2) / mG2 );
474  double atanDif2 = atanMax2 - atanMin2;
475  double wtDif2 = atanDif2 / (M_PI * NPOINT);
476 
477  // If on-shell decay forbidden then split integration range
478  // to ensure that low-mass region is not forgotten.
479  bool mustDiv = false;
480  double mDiv1 = 0.;
481  double atanDiv1 = 0.;
482  double atanDLo1 = 0.;
483  double atanDHi1 = 0.;
484  double wtDLo1 = 0.;
485  double wtDHi1 = 0.;
486  double mDiv2 = 0.;
487  double atanDiv2 = 0.;
488  double atanDLo2 = 0.;
489  double atanDHi2 = 0.;
490  double wtDLo2 = 0.;
491  double wtDHi2 = 0.;
492  if (m1 + m2 > mHatIn) {
493  mustDiv = true;
494  double tmpDiv = (mHatIn - m1 - m2) / (Gamma1 + Gamma2);
495  mDiv1 = m1 + Gamma1 * tmpDiv;
496  atanDiv1 = atan( (mDiv1 * mDiv1 - s1) / mG1 );
497  atanDLo1 = atanDiv1 - atanMin1;
498  atanDHi1 = atanMax1 - atanDiv1;
499  wtDLo1 = atanDLo1 / (M_PI * NPOINT);
500  wtDHi1 = atanDHi1 / (M_PI * NPOINT);
501  mDiv2 = m2 + Gamma2 * tmpDiv;
502  atanDiv2 = atan( (mDiv2 * mDiv2 - s2) / mG2 );
503  atanDLo2 = atanDiv2 - atanMin2;
504  atanDHi2 = atanMax2 - atanDiv2;
505  wtDLo2 = atanDLo2 / (M_PI * NPOINT);
506  wtDHi2 = atanDHi2 / (M_PI * NPOINT);
507  }
508 
509  // Step size in atan-mapped variable.
510  double xStep = 1. / NPOINT;
511  int nIter = (mustDiv) ? 2 * NPOINT : NPOINT;
512 
513  // Variables used in loop over integration points.
514  double sum = 0.;
515  double xNow1, sNow1, mNow1, mrNow1, xNow2, sNow2, mNow2, mrNow2, psNow,
516  value;
517  double wtNow1 = wtDif1;
518  double wtNow2 = wtDif2;
519 
520  // Outer loop with first-particle mass selection.
521  for (int ip1 = 0; ip1 < nIter; ++ip1) {
522  if (!mustDiv) {
523  xNow1 = xStep * (ip1 + 0.5);
524  sNow1 = s1 + mG1 * tan(atanMin1 + xNow1 * atanDif1);
525  } else if (ip1 < NPOINT) {
526  xNow1 = xStep * (ip1 + 0.5);
527  sNow1 = s1 + mG1 * tan(atanMin1 + xNow1 * atanDLo1);
528  wtNow1 = wtDLo1;
529  } else {
530  xNow1 = xStep * (ip1 - NPOINT + 0.5);
531  sNow1 = s1 + mG1 * tan(atanDiv1 + xNow1 * atanDHi1);
532  wtNow1 = wtDHi1;
533  }
534  mNow1 = min( mMax1, max( mMin1, sqrtpos(sNow1) ) );
535  mrNow1 = pow2(mNow1 / mHatIn);
536 
537  // Inner loop with second-particle mass selection.
538  for (int ip2 = 0; ip2 < nIter; ++ip2) {
539  if (!mustDiv) {
540  xNow2 = xStep * (ip2 + 0.5);
541  sNow2 = s2 + mG2 * tan(atanMin2 + xNow2 * atanDif2);
542  } else if (ip2 < NPOINT) {
543  xNow2 = xStep * (ip2 + 0.5);
544  sNow2 = s2 + mG2 * tan(atanMin2 + xNow2 * atanDLo2);
545  wtNow2 = wtDLo2;
546  } else {
547  xNow2 = xStep * (ip2 - NPOINT + 0.5);
548  sNow2 = s2 + mG2 * tan(atanDiv2 + xNow2 * atanDHi2);
549  wtNow2 = wtDHi2;
550  }
551  mNow2 = min( mMax2, max( mMin2, sqrtpos(sNow2) ) );
552  mrNow2 = pow2(mNow2 / mHatIn);
553 
554  // Check that point is inside phase space.
555  if (mNow1 + mNow2 > mHatIn) break;
556 
557  // Evaluate value and add to sum. Different matrix elements.
558  psNow = sqrtpos( pow2(1. - mrNow1 - mrNow2)
559  - 4. * mrNow1 * mrNow2);
560  value = 1.;
561  if (psMode == 1) value = psNow;
562  else if (psMode == 2) value = psNow * psNow;
563  else if (psMode == 3) value = pow3(psNow);
564  else if (psMode == 5) value = psNow
565  * (pow2(1. - mrNow1 - mrNow2) + 8. * mrNow1 * mrNow2);
566  else if (psMode == 6) value = pow3(psNow);
567  sum += value * wtNow1 * wtNow2;
568 
569  // End of second and first loop over integration points.
570  }
571  }
572 
573  // Done.
574  return sum;
575 }
576 
577 //==========================================================================
578 
579 // The ResonanceGmZ class.
580 // Derived class for gamma*/Z0 properties.
581 
582 //--------------------------------------------------------------------------
583 
584 // Initialize constants.
585 
586 void ResonanceGmZ::initConstants() {
587 
588  // Locally stored properties and couplings.
589  gmZmode = settingsPtr->mode("WeakZ0:gmZmode");
590  thetaWRat = 1. / (16. * couplingsPtr->sin2thetaW()
591  * couplingsPtr->cos2thetaW());
592 
593  // The Z0copy with id = 93 is a pure Z0.
594  if (idRes == 93) gmZmode = 2;
595 
596 }
597 
598 //--------------------------------------------------------------------------
599 
600 // Calculate various common prefactors for the current mass.
601 
602 void ResonanceGmZ::calcPreFac(bool calledFromInit) {
603 
604  // Common coupling factors.
605  alpEM = couplingsPtr->alphaEM(mHat * mHat);
606  alpS = couplingsPtr->alphaS(mHat * mHat);
607  colQ = 3. * (1. + alpS / M_PI);
608  preFac = alpEM * thetaWRat * mHat / 3.;
609 
610  // When call for incoming flavour need to consider gamma*/Z0 mix.
611  if (!calledFromInit) {
612 
613  // Couplings when an incoming fermion is specified; elso only pure Z0.
614  ei2 = 0.;
615  eivi = 0.;
616  vi2ai2 = 1.;
617  int idInFlavAbs = abs(idInFlav);
618  if (idInFlavAbs > 0 && idInFlavAbs < 19) {
619  ei2 = couplingsPtr->ef2(idInFlavAbs);
620  eivi = couplingsPtr->efvf(idInFlavAbs);
621  vi2ai2 = couplingsPtr->vf2af2(idInFlavAbs);
622  }
623 
624  // Calculate prefactors for gamma/interference/Z0 terms.
625  double sH = mHat * mHat;
626  gamNorm = ei2;
627  intNorm = 2. * eivi * thetaWRat * sH * (sH - m2Res)
628  / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
629  resNorm = vi2ai2 * pow2(thetaWRat * sH)
630  / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
631 
632  // Optionally only keep gamma* or Z0 term.
633  if (gmZmode == 1) {intNorm = 0.; resNorm = 0.;}
634  if (gmZmode == 2) {gamNorm = 0.; intNorm = 0.;}
635  }
636 
637 }
638 
639 //--------------------------------------------------------------------------
640 
641 // Calculate width for currently considered channel.
642 
643 void ResonanceGmZ::calcWidth(bool calledFromInit) {
644 
645  // Check that above threshold.
646  if (ps == 0.) return;
647 
648  // Only contributions from three fermion generations, except top.
649  if ( (id1Abs > 5 && id1Abs < 11) || id1Abs > 16 ) return;
650 
651  // At initialization only the pure Z0 should be considered.
652  if (calledFromInit) {
653 
654  // Combine kinematics with colour factor and couplings.
655  widNow = preFac * ps * (couplingsPtr->vf2(id1Abs) * (1. + 2. * mr1)
656  + couplingsPtr->af2(id1Abs) * ps*ps);
657  if (id1Abs < 6) widNow *= colQ;
658  }
659 
660  // When call for incoming flavour need to consider gamma*/Z0 mix.
661  else {
662 
663  // Kinematical factors and couplings.
664  double kinFacV = ps * (1. + 2. * mr1);
665  double ef2 = couplingsPtr->ef2(id1Abs) * kinFacV;
666  double efvf = couplingsPtr->efvf(id1Abs) * kinFacV;
667  double vf2af2 = couplingsPtr->vf2(id1Abs) * kinFacV
668  + couplingsPtr->af2(id1Abs) * pow3(ps);
669 
670  // Relative outwidths: combine instate, propagator and outstate.
671  widNow = gamNorm * ef2 + intNorm * efvf + resNorm * vf2af2;
672 
673  // Colour factor.
674  if (id1Abs < 6) widNow *= colQ;
675  }
676 
677 }
678 
679 //==========================================================================
680 
681 // The ResonanceW class.
682 // Derived class for W+- properties.
683 
684 //--------------------------------------------------------------------------
685 
686 // Initialize constants.
687 
688 void ResonanceW::initConstants() {
689 
690  // Locally stored properties and couplings.
691  thetaWRat = 1. / (12. * couplingsPtr->sin2thetaW());
692 
693 }
694 
695 //--------------------------------------------------------------------------
696 
697 // Calculate various common prefactors for the current mass.
698 
699 void ResonanceW::calcPreFac(bool) {
700 
701  // Common coupling factors.
702  alpEM = couplingsPtr->alphaEM(mHat * mHat);
703  alpS = couplingsPtr->alphaS(mHat * mHat);
704  colQ = 3. * (1. + alpS / M_PI);
705  preFac = alpEM * thetaWRat * mHat;
706 
707 }
708 
709 //--------------------------------------------------------------------------
710 
711 // Calculate width for currently considered channel.
712 
713 void ResonanceW::calcWidth(bool) {
714 
715  // Check that above threshold.
716  if (ps == 0.) return;
717 
718  // Only contributions from three fermion generations, except top.
719  if ( (id1Abs > 5 && id1Abs < 11) || id1Abs > 16 ) return;
720 
721 
722  // Combine kinematics with colour factor and couplings.
723  widNow = preFac * ps
724  * (1. - 0.5 * (mr1 + mr2) - 0.5 * pow2(mr1 - mr2));
725  if (id1Abs < 6) widNow *= colQ * couplingsPtr->V2CKMid(id1Abs, id2Abs);
726 
727 }
728 
729 //==========================================================================
730 
731 // The ResonanceTop class.
732 // Derived class for top/antitop properties.
733 
734 //--------------------------------------------------------------------------
735 
736 // Initialize constants.
737 
738 void ResonanceTop::initConstants() {
739 
740  // Locally stored properties and couplings.
741  thetaWRat = 1. / (16. * couplingsPtr->sin2thetaW());
742  m2W = pow2(particleDataPtr->m0(24));
743 
744  // Extra coupling factors for t -> H+ + b.
745  tanBeta = settingsPtr->parm("HiggsHchg:tanBeta");
746  tan2Beta = tanBeta * tanBeta;
747  mbRun = particleDataPtr->mRun( 5, particleDataPtr->m0(6) );
748 
749 }
750 
751 //--------------------------------------------------------------------------
752 
753 // Calculate various common prefactors for the current mass.
754 
755 void ResonanceTop::calcPreFac(bool) {
756 
757  // Common coupling factors.
758  alpEM = couplingsPtr->alphaEM(mHat * mHat);
759  alpS = couplingsPtr->alphaS(mHat * mHat);
760  colQ = 1. - 2.5 * alpS / M_PI;
761  preFac = alpEM * thetaWRat * pow3(mHat) / m2W;
762 
763 }
764 
765 //--------------------------------------------------------------------------
766 
767 // Calculate width for currently considered channel.
768 
769 void ResonanceTop::calcWidth(bool) {
770 
771 
772  // Check that above threshold.
773  if (ps == 0.) return;
774 
775  // Contributions from W + quark.
776  if (id1Abs == 24 && id2Abs < 6) {
777  widNow = preFac * ps
778  * ( pow2(1. - mr2) + (1. + mr2) * mr1 - 2. * mr1 * mr1 );
779 
780  // Combine with colour factor and CKM couplings.
781  widNow *= colQ * couplingsPtr->V2CKMid(6, id2Abs);
782 
783  // Contributions from H+ + quark (so far only b).
784  } else if (id1Abs == 37 && id2Abs == 5) {
785  widNow = preFac * ps * ( (1. + mr2 - mr1)
786  * (pow2(mbRun / mHat) * tan2Beta + 1. / tan2Beta)
787  + 4. * mbRun * mf2 / pow2(mHat) );
788  }
789 
790 }
791 
792 //==========================================================================
793 
794 // The ResonanceFour class.
795 // Derived class for fourth-generation properties.
796 
797 //--------------------------------------------------------------------------
798 
799 // Initialize constants.
800 
801 void ResonanceFour::initConstants() {
802 
803  // Locally stored properties and couplings.
804  thetaWRat = 1. / (16. * couplingsPtr->sin2thetaW());
805  m2W = pow2(particleDataPtr->m0(24));
806 
807 }
808 
809 //--------------------------------------------------------------------------
810 
811 // Calculate various common prefactors for the current mass.
812 
813 void ResonanceFour::calcPreFac(bool) {
814 
815  // Common coupling factors.
816  alpEM = couplingsPtr->alphaEM(mHat * mHat);
817  alpS = couplingsPtr->alphaS(mHat * mHat);
818  colQ = (idRes < 9) ? 1. - 2.5 * alpS / M_PI : 1.;
819  preFac = alpEM * thetaWRat * pow3(mHat) / m2W;
820 
821 }
822 
823 //--------------------------------------------------------------------------
824 
825 // Calculate width for currently considered channel.
826 
827 void ResonanceFour::calcWidth(bool) {
828 
829  // Only contributions from W + fermion.
830  if (id1Abs != 24 || id2Abs > 18) return;
831 
832  // Check that above threshold. Kinematical factor.
833  if (ps == 0.) return;
834  widNow = preFac * ps
835  * ( pow2(1. - mr2) + (1. + mr2) * mr1 - 2. * mr1 * mr1 );
836 
837  // Combine with colour factor and CKM couplings.
838  if (idRes < 9) widNow *= colQ * couplingsPtr->V2CKMid(idRes, id2Abs);
839 
840 }
841 
842 //==========================================================================
843 
844 // The ResonanceH class.
845 // Derived class for SM and BSM Higgs properties.
846 
847 //--------------------------------------------------------------------------
848 
849 // Constants: could be changed here if desired, but normally should not.
850 // These are of technical nature, as described for each.
851 
852 // Minimal mass for W, Z, top in integration over respective Breit-Wigner.
853 // Top constrainted by t -> W b decay, which is not seen in simple top BW.
854 const double ResonanceH::MASSMINWZ = 10.;
855 const double ResonanceH::MASSMINT = 100.;
856 
857 // Number of widths above threshold where B-W integration not needed.
858 const double ResonanceH::GAMMAMARGIN = 10.;
859 
860 //--------------------------------------------------------------------------
861 
862 // Initialize constants.
863 
864 void ResonanceH::initConstants() {
865 
866  // Locally stored properties and couplings.
867  useCubicWidth = settingsPtr->flag("Higgs:cubicWidth");
868  useRunLoopMass = settingsPtr->flag("Higgs:runningLoopMass");
869  sin2tW = couplingsPtr->sin2thetaW();
870  cos2tW = 1. - sin2tW;
871  mT = particleDataPtr->m0(6);
872  mZ = particleDataPtr->m0(23);
873  mW = particleDataPtr->m0(24);
874  mHchg = particleDataPtr->m0(37);
875  GammaT = particleDataPtr->mWidth(6);
876  GammaZ = particleDataPtr->mWidth(23);
877  GammaW = particleDataPtr->mWidth(24);
878 
879  // NLO corrections to SM Higgs width, rescaled to reference alpha_S value.
880  useNLOWidths = (higgsType == 0) && settingsPtr->flag("HiggsSM:NLOWidths");
881  rescAlpS = 0.12833 / couplingsPtr->alphaS(125. * 125.);
882  rescColQ = 1.;
883 
884  // Couplings to fermions, Z and W, depending on Higgs type.
885  coup2d = 1.;
886  coup2u = 1.;
887  coup2l = 1.;
888  coup2Z = 1.;
889  coup2W = 1.;
890  coup2Hchg = 0.;
891  coup2H1H1 = 0.;
892  coup2A3A3 = 0.;
893  coup2H1Z = 0.;
894  coup2A3Z = 0.;
895  coup2A3H1 = 0.;
896  coup2HchgW = 0.;
897  if (higgsType == 1) {
898  coup2d = settingsPtr->parm("HiggsH1:coup2d");
899  coup2u = settingsPtr->parm("HiggsH1:coup2u");
900  coup2l = settingsPtr->parm("HiggsH1:coup2l");
901  coup2Z = settingsPtr->parm("HiggsH1:coup2Z");
902  coup2W = settingsPtr->parm("HiggsH1:coup2W");
903  coup2Hchg = settingsPtr->parm("HiggsH1:coup2Hchg");
904  } else if (higgsType == 2) {
905  coup2d = settingsPtr->parm("HiggsH2:coup2d");
906  coup2u = settingsPtr->parm("HiggsH2:coup2u");
907  coup2l = settingsPtr->parm("HiggsH2:coup2l");
908  coup2Z = settingsPtr->parm("HiggsH2:coup2Z");
909  coup2W = settingsPtr->parm("HiggsH2:coup2W");
910  coup2Hchg = settingsPtr->parm("HiggsH2:coup2Hchg");
911  coup2H1H1 = settingsPtr->parm("HiggsH2:coup2H1H1");
912  coup2A3A3 = settingsPtr->parm("HiggsH2:coup2A3A3");
913  coup2H1Z = settingsPtr->parm("HiggsH2:coup2H1Z");
914  coup2A3Z = settingsPtr->parm("HiggsA3:coup2H2Z");
915  coup2A3H1 = settingsPtr->parm("HiggsH2:coup2A3H1");
916  coup2HchgW = settingsPtr->parm("HiggsH2:coup2HchgW");
917  } else if (higgsType == 3) {
918  coup2d = settingsPtr->parm("HiggsA3:coup2d");
919  coup2u = settingsPtr->parm("HiggsA3:coup2u");
920  coup2l = settingsPtr->parm("HiggsA3:coup2l");
921  coup2Z = settingsPtr->parm("HiggsA3:coup2Z");
922  coup2W = settingsPtr->parm("HiggsA3:coup2W");
923  coup2Hchg = settingsPtr->parm("HiggsA3:coup2Hchg");
924  coup2H1H1 = settingsPtr->parm("HiggsA3:coup2H1H1");
925  coup2H1Z = settingsPtr->parm("HiggsA3:coup2H1Z");
926  coup2HchgW = settingsPtr->parm("HiggsA3:coup2Hchg");
927  }
928 
929  // Initialization of threshold kinematical factor by stepwise
930  // numerical integration of H -> t tbar, Z0 Z0 and W+ W-.
931  int psModeT = (higgsType < 3) ? 3 : 1;
932  int psModeWZ = (higgsType < 3) ? 5 : 6;
933  mLowT = max( 2.02 * MASSMINT, 0.5 * mT);
934  mStepT = 0.01 * (3. * mT - mLowT);
935  mLowZ = max( 2.02 * MASSMINWZ, 0.5 * mZ);
936  mStepZ = 0.01 * (3. * mZ - mLowZ);
937  mLowW = max( 2.02 * MASSMINWZ, 0.5 * mW);
938  mStepW = 0.01 * (3. * mW - mLowW);
939  for (int i = 0; i <= 100; ++i) {
940  kinFacT[i] = numInt2BW( mLowT + i * mStepT,
941  mT, GammaT, MASSMINT, mT, GammaT, MASSMINT, psModeT);
942  kinFacZ[i] = numInt2BW( mLowZ + i * mStepZ,
943  mZ, GammaZ, MASSMINWZ, mZ, GammaZ, MASSMINWZ, psModeWZ);
944  kinFacW[i] = numInt2BW( mLowW + i * mStepW,
945  mW, GammaW, MASSMINWZ, mW, GammaW, MASSMINWZ, psModeWZ);
946  }
947 
948 }
949 
950 //--------------------------------------------------------------------------
951 
952 // Calculate various common prefactors for the current mass.
953 
954 void ResonanceH::calcPreFac(bool) {
955 
956  // Common coupling factors.
957  alpEM = couplingsPtr->alphaEM(mHat * mHat);
958  alpS = couplingsPtr->alphaS(mHat * mHat);
959  colQ = 3. * (1. + alpS / M_PI);
960  preFac = (alpEM / (8. * sin2tW)) * pow3(mHat) / pow2(mW);
961  if (useNLOWidths) rescColQ = 3. * (1. + rescAlpS * alpS / M_PI) / colQ;
962 
963 }
964 
965 //--------------------------------------------------------------------------
966 
967 // Calculate width for currently considered channel.
968 
969 void ResonanceH::calcWidth(bool) {
970 
971  // Widths of decays Higgs -> f + fbar.
972  if ( id2Abs == id1Abs && ( (id1Abs > 0 && id1Abs < 7)
973  || (id1Abs > 10 && id1Abs < 17) ) ) {
974  kinFac = 0.;
975 
976  // Check that above threshold (well above for top). Kinematical factor.
977  if ( (id1Abs != 6 && mHat > 2. * mf1 + MASSMARGIN)
978  || (id1Abs == 6 && mHat > 3. * mT ) ) {
979  // A0 behaves like beta, h0 and H0 like beta**3.
980  kinFac = (higgsType < 3) ? pow3(ps) : ps;
981  }
982 
983  // Top near or below threshold: interpolate in table.
984  else if (id1Abs == 6 && mHat > mLowT) {
985  double xTab = (mHat - mLowT) / mStepT;
986  int iTab = max( 0, min( 99, int(xTab) ) );
987  kinFac = kinFacT[iTab]
988  * pow( kinFacT[iTab + 1] / kinFacT[iTab], xTab - iTab);
989  }
990 
991  // Coupling from mass and from BSM deviation from SM.
992  double coupFac = pow2(particleDataPtr->mRun(id1Abs, mHat) / mHat);
993  if (id1Abs < 7 && id1Abs%2 == 1) coupFac *= coup2d * coup2d;
994  else if (id1Abs < 7) coupFac *= coup2u * coup2u;
995  else coupFac *= coup2l * coup2l;
996 
997  // Combine couplings and phase space with colour factor.
998  widNow = preFac * coupFac * kinFac;
999  if (id1Abs < 7) widNow *= colQ;
1000  }
1001 
1002  // Widths of decays Higgs -> g + g.
1003  else if (id1Abs == 21 && id2Abs == 21)
1004  widNow = preFac * pow2(alpS / M_PI) * eta2gg();
1005 
1006  // Widths of decays Higgs -> gamma + gamma.
1007  else if (id1Abs == 22 && id2Abs == 22)
1008  widNow = preFac * pow2(alpEM / M_PI) * 0.5 * eta2gaga();
1009 
1010  // Widths of decays Higgs -> Z0 + gamma0.
1011  else if (id1Abs == 23 && id2Abs == 22)
1012  widNow = preFac * pow2(alpEM / M_PI) * pow3(ps) * eta2gaZ();
1013 
1014  // Widths of decays Higgs (h0, H0) -> Z0 + Z0.
1015  else if (id1Abs == 23 && id2Abs == 23) {
1016  // If Higgs heavy use on-shell expression, else interpolation in table
1017  if (mHat > 3. * mZ) kinFac = (1. - 4. * mr1 + 12. * mr1 * mr1) * ps;
1018  else if (mHat > mLowZ) {
1019  double xTab = (mHat - mLowZ) / mStepZ;
1020  int iTab = max( 0, min( 99, int(xTab) ) );
1021  kinFac = kinFacZ[iTab]
1022  * pow( kinFacZ[iTab + 1] / kinFacZ[iTab], xTab - iTab );
1023  }
1024  else kinFac = 0.;
1025  // Prefactor, normally rescaled to mRes^2 * mHat rather than mHat^3.
1026  widNow = 0.25 * preFac * pow2(coup2Z) * kinFac;
1027  if (!useCubicWidth) widNow *= pow2(mRes / mHat);
1028  }
1029 
1030  // Widths of decays Higgs (h0, H0) -> W+ + W-.
1031  else if (id1Abs == 24 && id2Abs == 24) {
1032  // If Higgs heavy use on-shell expression, else interpolation in table.
1033  if (mHat > 3. * mW) kinFac = (1. - 4. * mr1 + 12. * mr1 * mr1) * ps;
1034  else if (mHat > mLowW) {
1035  double xTab = (mHat - mLowW) / mStepW;
1036  int iTab = max( 0, min( 99, int(xTab) ) );
1037  kinFac = kinFacW[iTab]
1038  * pow( kinFacW[iTab + 1] / kinFacW[iTab], xTab - iTab);
1039  }
1040  else kinFac = 0.;
1041  // Prefactor, normally rescaled to mRes^2 * mHat rather than mHat^3.
1042  widNow = 0.5 * preFac * pow2(coup2W) * kinFac;
1043  if (!useCubicWidth) widNow *= pow2(mRes / mHat);
1044  }
1045 
1046  // Widths of decays Higgs (H0) -> h0 + h0.
1047  else if (id1Abs == 25 && id2Abs == 25)
1048  widNow = 0.25 * preFac * pow4(mZ / mHat) * ps * pow2(coup2H1H1);
1049 
1050  // Widths of decays Higgs (H0) -> A0 + A0.
1051  else if (id1Abs == 36 && id2Abs == 36)
1052  widNow = 0.5 * preFac * pow4(mZ / mHat) * ps * pow2(coup2A3A3);
1053 
1054  // Widths of decays Higgs (A0) -> h0 + Z0.
1055  else if (id1Abs == 25 && id2Abs == 23)
1056  widNow = 0.5 * preFac * pow3(ps) * pow2(coup2H1Z);
1057 
1058  // Widths of decays Higgs (H0) -> A0 + Z0.
1059  else if (id1Abs == 36 && id2Abs == 23)
1060  widNow = 0.5 * preFac * pow3(ps) * pow2(coup2A3Z);
1061 
1062  // Widths of decays Higgs (H0) -> A0 + h0.
1063  else if (id1Abs == 36 && id2Abs == 25)
1064  widNow = 0.25 * preFac * pow4(mZ / mHat) * ps * pow2(coup2A3H1);
1065 
1066  // Widths of decays Higgs -> H+- + W-+.
1067  else if (id1Abs == 37 && id2Abs == 24)
1068  widNow = 0.5 * preFac * pow3(ps) * pow2(coup2HchgW);
1069 
1070  // NLO multiplicative factors for SM h0 (125 GeV) based on LHCXSWG
1071  // recommendations.
1072  if (useNLOWidths) {
1073  if (id1Abs == 21 && id2Abs == 21) widNow *= 1.47 * pow2(rescAlpS);
1074  else if (id1Abs == 22 && id2Abs == 22) widNow *= 0.88;
1075  else if (id1Abs == 22 && id2Abs == 23) widNow *= 0.95;
1076  else if (id1Abs == 23 && id2Abs == 23) widNow *= 1.10;
1077  else if (id1Abs == 24 && id2Abs == 24) widNow *= 1.09;
1078  else if (id1Abs == 5 && id2Abs == 5) widNow *= 1.07 * rescColQ;
1079  else if (id1Abs == 4 && id2Abs == 4) widNow *= 0.937 * rescColQ;
1080  else if (id1Abs == 13 && id2Abs == 13) widNow *= 0.974;
1081  else if (id1Abs == 15 && id2Abs == 15) widNow *= 0.992;
1082  }
1083 
1084 }
1085 
1086 //--------------------------------------------------------------------------
1087 
1088 // Sum up quark loop contributions in Higgs -> g + g.
1089 // Note: running quark masses are used, unlike Pythia6 (not negligible shift).
1090 
1091 double ResonanceH::eta2gg() {
1092 
1093  // Initial values.
1094  complex eta = complex(0., 0.);
1095  double mLoop, epsilon, root, rootLog;
1096  complex phi, etaNow;
1097 
1098  // Loop over s, c, b, t quark flavours.
1099  for (int idNow = 3; idNow < 7; ++idNow) {
1100  mLoop = (useRunLoopMass) ? particleDataPtr->mRun(idNow, mHat)
1101  : particleDataPtr->m0(idNow);
1102  epsilon = pow2(2. * mLoop / mHat);
1103 
1104  // Value of loop integral.
1105  if (epsilon <= 1.) {
1106  root = sqrt(1. - epsilon);
1107  rootLog = (epsilon < 1e-4) ? log(4. / epsilon - 2.)
1108  : log( (1. + root) / (1. - root) );
1109  phi = complex( -0.25 * (pow2(rootLog) - pow2(M_PI)),
1110  0.5 * M_PI * rootLog );
1111  }
1112  else phi = complex( pow2( asin(1. / sqrt(epsilon)) ), 0.);
1113 
1114  // Factors that depend on Higgs and flavour type.
1115  if (higgsType < 3) etaNow = -0.5 * epsilon
1116  * (complex(1., 0.) + (1. - epsilon) * phi);
1117  else etaNow = -0.5 * epsilon * phi;
1118  if (idNow%2 == 1) etaNow *= coup2d;
1119  else etaNow *= coup2u;
1120 
1121  // Sum up contribution and return square of absolute value.
1122  eta += etaNow;
1123  }
1124  return (pow2(eta.real()) + pow2(eta.imag()));
1125 
1126 }
1127 
1128 //--------------------------------------------------------------------------
1129 
1130 // Sum up quark, lepton, W+- and (for BSM) H+- loop contributions
1131 // in Higgs -> gamma + gamma.
1132 
1133 double ResonanceH::eta2gaga() {
1134 
1135  // Initial values.
1136  complex eta = complex(0., 0.);
1137  int idNow;
1138  double ef, mLoop, epsilon, root, rootLog;
1139  complex phi, etaNow;
1140 
1141  // Loop over s, c, b, t, mu, tau, W+-, H+- flavours.
1142  for (int idLoop = 0; idLoop < 8; ++idLoop) {
1143  if (idLoop < 4) idNow = idLoop + 3;
1144  else if (idLoop < 6) idNow = 2 * idLoop + 5;
1145  else if (idLoop < 7) idNow = 24;
1146  else idNow = 37;
1147  if (idNow == 37 && higgsType == 0) continue;
1148 
1149  // Charge and loop integral parameter.
1150  ef = (idNow < 20) ? couplingsPtr->ef(idNow) : 1.;
1151  mLoop = (useRunLoopMass) ? particleDataPtr->mRun(idNow, mHat)
1152  : particleDataPtr->m0(idNow);
1153  epsilon = pow2(2. * mLoop / mHat);
1154 
1155  // Value of loop integral.
1156  if (epsilon <= 1.) {
1157  root = sqrt(1. - epsilon);
1158  rootLog = (epsilon < 1e-4) ? log(4. / epsilon - 2.)
1159  : log( (1. + root) / (1. - root) );
1160  phi = complex( -0.25 * (pow2(rootLog) - pow2(M_PI)),
1161  0.5 * M_PI * rootLog );
1162  }
1163  else phi = complex( pow2( asin(1. / sqrt(epsilon)) ), 0.);
1164 
1165  // Expressions for quarks and leptons that depend on Higgs type.
1166  if (idNow < 17) {
1167  if (higgsType < 3) etaNow = -0.5 * epsilon
1168  * (complex(1., 0.) + (1. - epsilon) * phi);
1169  else etaNow = -0.5 * epsilon * phi;
1170  if (idNow < 7 && idNow%2 == 1) etaNow *= 3. * pow2(ef) * coup2d;
1171  else if (idNow < 7 ) etaNow *= 3. * pow2(ef) * coup2u;
1172  else etaNow *= pow2(ef) * coup2l;
1173  }
1174 
1175  // Expression for W+-.
1176  else if (idNow == 24) etaNow = (complex(0.5 + 0.75 * epsilon, 0.)
1177  + 0.75 * epsilon * (2. - epsilon) * phi) * coup2W;
1178 
1179  // Expression for H+-.
1180  else etaNow = (complex(epsilon, 0.) - epsilon * epsilon * phi)
1181  * pow2(mW / mHchg) * coup2Hchg;
1182 
1183  // Sum up contribution and return square of absolute value.
1184  eta += etaNow;
1185  }
1186  return (pow2(eta.real()) + pow2(eta.imag()));
1187 
1188 }
1189 
1190 //--------------------------------------------------------------------------
1191 
1192 // Sum up quark, lepton, W+- and (for BSM) H+- loop contributions
1193 // in Higgs -> gamma + Z0.
1194 
1195 double ResonanceH::eta2gaZ() {
1196 
1197  // Initial values.
1198  complex eta = complex(0., 0.);
1199  int idNow;
1200  double ef, vf, mLoop, epsilon, epsPrime, root, rootLog, asinEps;
1201  complex phi, psi, phiPrime, psiPrime, fXY, f1, etaNow;
1202 
1203  // Loop over s, c, b, t, mu , tau, W+-, H+- flavours.
1204  for (int idLoop = 0; idLoop < 7; ++idLoop) {
1205  if (idLoop < 4) idNow = idLoop + 3;
1206  else if (idLoop < 6) idNow = 2 * idLoop + 5;
1207  else if (idLoop < 7) idNow = 24;
1208  else idNow = 37;
1209 
1210  // Electroweak charges and loop integral parameters.
1211  ef = (idNow < 20) ? couplingsPtr->ef(idNow) : 1.;
1212  vf = (idNow < 20) ? couplingsPtr->vf(idNow) : 0.;
1213  mLoop = (useRunLoopMass) ? particleDataPtr->mRun(idNow, mHat)
1214  : particleDataPtr->m0(idNow);
1215  epsilon = pow2(2. * mLoop / mHat);
1216  epsPrime = pow2(2. * mLoop / mZ);
1217 
1218  // Value of loop integral for epsilon = 4 m^2 / sHat.
1219  if (epsilon <= 1.) {
1220  root = sqrt(1. - epsilon);
1221  rootLog = (epsilon < 1e-4) ? log(4. / epsilon - 2.)
1222  : log( (1. + root) / (1. - root) );
1223  phi = complex( -0.25 * (pow2(rootLog) - pow2(M_PI)),
1224  0.5 * M_PI * rootLog );
1225  psi = 0.5 * root * complex( rootLog, -M_PI);
1226  } else {
1227  asinEps = asin(1. / sqrt(epsilon));
1228  phi = complex( pow2(asinEps), 0.);
1229  psi = complex( sqrt(epsilon - 1.) * asinEps, 0.);
1230  }
1231 
1232  // Value of loop integral for epsilonPrime = 4 m^2 / m_Z^2.
1233  if (epsPrime <= 1.) {
1234  root = sqrt(1. - epsPrime);
1235  rootLog = (epsPrime < 1e-4) ? log(4. / epsPrime - 2.)
1236  : log( (1. + root) / (1. - root) );
1237  phiPrime = complex( -0.25 * (pow2(rootLog) - pow2(M_PI)),
1238  0.5 * M_PI * rootLog );
1239  psiPrime = 0.5 * root * complex( rootLog, -M_PI);
1240  } else {
1241  asinEps = asin(1. / sqrt(epsPrime));
1242  phiPrime = complex( pow2(asinEps), 0.);
1243  psiPrime = complex( sqrt(epsPrime - 1.) * asinEps, 0.);
1244  }
1245 
1246  // Combine the two loop integrals.
1247  fXY = (epsilon * epsPrime / (8. * pow2(epsilon - epsPrime)))
1248  * ( complex(epsilon - epsPrime, 0)
1249  + epsilon * epsPrime * (phi - phiPrime)
1250  + 2. * epsilon * (psi - psiPrime) );
1251  f1 = - (epsilon * epsPrime / (2. * (epsilon - epsPrime)))
1252  * (phi - phiPrime);
1253 
1254  // Expressions for quarks and leptons that depend on Higgs type.
1255  if (idNow < 17) {
1256  etaNow = (higgsType < 3) ? -fXY + 0.25 * f1 : 0.25 * f1;
1257  if (idNow < 7 && idNow%2 == 1) etaNow *= 3. * ef * vf * coup2d;
1258  else if (idNow < 7) etaNow *= 3. * ef * vf * coup2u;
1259  else etaNow *= ef * vf * coup2l;
1260 
1261  // Expression for W+-.
1262  } else if (idNow == 24) {
1263  double coef1 = 3. - sin2tW / cos2tW;
1264  double coefXY = (1. + 2. / epsilon) * sin2tW / cos2tW
1265  - (5. + 2. / epsilon);
1266  etaNow = -cos2tW * (coef1 * f1 + coefXY * fXY) * coup2W;
1267 
1268  // Expression for H+-.
1269  } else etaNow = (1. - 2. * sin2tW) * fXY * pow2(mW / mHchg)
1270  * coup2Hchg;
1271 
1272  // Sum up contribution and return square of absolute value.
1273  eta += etaNow;
1274  }
1275  return ( (pow2(eta.real()) + pow2(eta.imag())) / (sin2tW * cos2tW) );
1276 
1277 }
1278 
1279 //==========================================================================
1280 
1281 // The ResonanceHchg class.
1282 // Derived class for H+- properties.
1283 
1284 //--------------------------------------------------------------------------
1285 
1286 // Initialize constants.
1287 
1288 void ResonanceHchg::initConstants() {
1289 
1290  // Locally stored properties and couplings.
1291  useCubicWidth = settingsPtr->flag("Higgs:cubicWidth");
1292  thetaWRat = 1. / (8. * couplingsPtr->sin2thetaW());
1293  mW = particleDataPtr->m0(24);
1294  tanBeta = settingsPtr->parm("HiggsHchg:tanBeta");
1295  tan2Beta = tanBeta * tanBeta;
1296  coup2H1W = settingsPtr->parm("HiggsHchg:coup2H1W");
1297 
1298 }
1299 
1300 //--------------------------------------------------------------------------
1301 
1302 // Calculate various common prefactors for the current mass.
1303 
1304 void ResonanceHchg::calcPreFac(bool) {
1305 
1306  // Common coupling factors.
1307  alpEM = couplingsPtr->alphaEM(mHat * mHat);
1308  alpS = couplingsPtr->alphaS(mHat * mHat);
1309  colQ = 3. * (1. + alpS / M_PI);
1310  preFac = alpEM * thetaWRat * pow3(mHat) / pow2(mW);
1311 
1312 }
1313 
1314 //--------------------------------------------------------------------------
1315 
1316 // Calculate width for currently considered channel.
1317 
1318 void ResonanceHchg::calcWidth(bool) {
1319 
1320  // Check that above threshold.
1321  if (ps == 0.) return;
1322 
1323  // H+- decay to fermions involves running masses.
1324  if (id1Abs < 17 && (id1Abs < 7 || id1Abs > 10)) {
1325  double mRun1 = particleDataPtr->mRun(id1Abs, mHat);
1326  double mRun2 = particleDataPtr->mRun(id2Abs, mHat);
1327  double mrRunDn = pow2(mRun1 / mHat);
1328  double mrRunUp = pow2(mRun2 / mHat);
1329  if (id1Abs%2 == 0) swap( mrRunDn, mrRunUp);
1330 
1331  // Width to fermions: couplings, kinematics, colour factor.
1332  widNow = preFac * max( 0., (mrRunDn * tan2Beta + mrRunUp / tan2Beta)
1333  * (1. - mrRunDn - mrRunUp) - 4. *mrRunDn * mrRunUp ) * ps;
1334  if (id1Abs < 7) widNow *= colQ;
1335  }
1336 
1337  // H+- decay to h0 + W+-.
1338  else if (id1Abs == 25 && id2Abs == 24)
1339  widNow = 0.5 * preFac * pow3(ps) * pow2(coup2H1W);
1340 
1341 }
1342 
1343 //==========================================================================
1344 
1345 // The ResonanceZprime class.
1346 // Derived class for gamma*/Z0/Z'^0 properties.
1347 
1348 //--------------------------------------------------------------------------
1349 
1350 // Initialize constants.
1351 
1352 void ResonanceZprime::initConstants() {
1353 
1354  // Locally stored properties and couplings.
1355  gmZmode = settingsPtr->mode("Zprime:gmZmode");
1356  sin2tW = couplingsPtr->sin2thetaW();
1357  cos2tW = 1. - sin2tW;
1358  thetaWRat = 1. / (16. * sin2tW * cos2tW);
1359 
1360  // Properties of Z resonance.
1361  mZ = particleDataPtr->m0(23);
1362  GammaZ = particleDataPtr->mWidth(23);
1363  m2Z = mZ*mZ;
1364  GamMRatZ = GammaZ / mZ;
1365 
1366  // Ensure that arrays initially empty.
1367  for (int i = 0; i < 20; ++i) afZp[i] = 0.;
1368  for (int i = 0; i < 20; ++i) vfZp[i] = 0.;
1369 
1370  // Store first-generation axial and vector couplings.
1371  afZp[1] = settingsPtr->parm("Zprime:ad");
1372  afZp[2] = settingsPtr->parm("Zprime:au");
1373  afZp[11] = settingsPtr->parm("Zprime:ae");
1374  afZp[12] = settingsPtr->parm("Zprime:anue");
1375  vfZp[1] = settingsPtr->parm("Zprime:vd");
1376  vfZp[2] = settingsPtr->parm("Zprime:vu");
1377  vfZp[11] = settingsPtr->parm("Zprime:ve");
1378  vfZp[12] = settingsPtr->parm("Zprime:vnue");
1379 
1380  // Second and third generation could be carbon copy of this...
1381  if (settingsPtr->flag("Zprime:universality")) {
1382  for (int i = 3; i <= 6; ++i) {
1383  afZp[i] = afZp[i-2];
1384  vfZp[i] = vfZp[i-2];
1385  afZp[i+10] = afZp[i+8];
1386  vfZp[i+10] = vfZp[i+8];
1387  }
1388 
1389  // ... or could have different couplings.
1390  } else {
1391  afZp[3] = settingsPtr->parm("Zprime:as");
1392  afZp[4] = settingsPtr->parm("Zprime:ac");
1393  afZp[5] = settingsPtr->parm("Zprime:ab");
1394  afZp[6] = settingsPtr->parm("Zprime:at");
1395  afZp[13] = settingsPtr->parm("Zprime:amu");
1396  afZp[14] = settingsPtr->parm("Zprime:anumu");
1397  afZp[15] = settingsPtr->parm("Zprime:atau");
1398  afZp[16] = settingsPtr->parm("Zprime:anutau");
1399  vfZp[3] = settingsPtr->parm("Zprime:vs");
1400  vfZp[4] = settingsPtr->parm("Zprime:vc");
1401  vfZp[5] = settingsPtr->parm("Zprime:vb");
1402  vfZp[6] = settingsPtr->parm("Zprime:vt");
1403  vfZp[13] = settingsPtr->parm("Zprime:vmu");
1404  vfZp[14] = settingsPtr->parm("Zprime:vnumu");
1405  vfZp[15] = settingsPtr->parm("Zprime:vtau");
1406  vfZp[16] = settingsPtr->parm("Zprime:vnutau");
1407  }
1408 
1409  // Coupling for Z' -> W+ W-.
1410  coupZpWW = settingsPtr->parm("Zprime:coup2WW");
1411 
1412 }
1413 
1414 //--------------------------------------------------------------------------
1415 
1416 // Calculate various common prefactors for the current mass.
1417 
1418 void ResonanceZprime::calcPreFac(bool calledFromInit) {
1419 
1420  // Common coupling factors.
1421  alpEM = couplingsPtr->alphaEM(mHat * mHat);
1422  alpS = couplingsPtr->alphaS(mHat * mHat);
1423  colQ = 3. * (1. + alpS / M_PI);
1424  preFac = alpEM * thetaWRat * mHat / 3.;
1425 
1426  // When call for incoming flavour need to consider gamma*/Z0 mix.
1427  if (!calledFromInit) {
1428 
1429  // Couplings when an incoming fermion is specified; elso only pure Z'0.
1430  ei2 = 0.;
1431  eivi = 0.;
1432  vai2 = 0.;
1433  eivpi = 0.;
1434  vaivapi = 0.,
1435  vapi2 = 1.;
1436  int idInFlavAbs = abs(idInFlav);
1437  if (idInFlavAbs > 0 && idInFlavAbs < 19) {
1438  double ei = couplingsPtr->ef(idInFlavAbs);
1439  double ai = couplingsPtr->af(idInFlavAbs);
1440  double vi = couplingsPtr->vf(idInFlavAbs);
1441  double api = afZp[idInFlavAbs];
1442  double vpi = vfZp[idInFlavAbs];
1443  ei2 = ei * ei;
1444  eivi = ei * vi;
1445  vai2 = vi * vi + ai * ai;
1446  eivpi = ei * vpi;
1447  vaivapi = vi * vpi + ai * api;;
1448  vapi2 = vpi * vpi + api * api;
1449  }
1450 
1451  // Calculate prefactors for gamma/interference/Z0 terms.
1452  double sH = mHat * mHat;
1453  double propZ = sH / ( pow2(sH - m2Z) + pow2(sH * GamMRatZ) );
1454  double propZp = sH / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
1455  gamNorm = ei2;
1456  gamZNorm = 2. * eivi * thetaWRat * (sH - m2Z) * propZ;
1457  ZNorm = vai2 * pow2(thetaWRat) * sH * propZ;
1458  gamZpNorm = 2. * eivpi * thetaWRat * (sH - m2Res) * propZp;
1459  ZZpNorm = 2. * vaivapi * pow2(thetaWRat) * ((sH - m2Res) * (sH - m2Z)
1460  + sH * GamMRat * sH * GamMRatZ) * propZ * propZp;
1461  ZpNorm = vapi2 * pow2(thetaWRat) * sH * propZp;
1462 
1463  // Optionally only keep some of gamma*, Z0 and Z' terms.
1464  if (gmZmode == 1) {gamZNorm = 0; ZNorm = 0.; gamZpNorm = 0.;
1465  ZZpNorm = 0.; ZpNorm = 0.;}
1466  if (gmZmode == 2) {gamNorm = 0.; gamZNorm = 0.; gamZpNorm = 0.;
1467  ZZpNorm = 0.; ZpNorm = 0.;}
1468  if (gmZmode == 3) {gamNorm = 0.; gamZNorm = 0.; ZNorm = 0.;
1469  gamZpNorm = 0.; ZZpNorm = 0.;}
1470  if (gmZmode == 4) {gamZpNorm = 0.; ZZpNorm = 0.; ZpNorm = 0.;}
1471  if (gmZmode == 5) {gamZNorm = 0.; ZNorm = 0.; ZZpNorm = 0.;}
1472  if (gmZmode == 6) {gamNorm = 0.; gamZNorm = 0.; gamZpNorm = 0.;}
1473  }
1474 
1475 }
1476 
1477 //--------------------------------------------------------------------------
1478 
1479 // Calculate width for currently considered channel.
1480 
1481 void ResonanceZprime::calcWidth(bool calledFromInit) {
1482 
1483  // Check that above threshold.
1484  if (ps == 0.) return;
1485 
1486  // At initialization only the pure Z'0 should be considered.
1487  if (calledFromInit) {
1488 
1489  // Contributions from three fermion generations.
1490  if ( id1Abs < 7 || (id1Abs > 10 && id1Abs < 17) ) {
1491  double apf = afZp[id1Abs];
1492  double vpf = vfZp[id1Abs];
1493  widNow = preFac * ps * (vpf*vpf * (1. + 2. * mr1)
1494  + apf*apf * ps*ps);
1495  if (id1Abs < 7) widNow *= colQ;
1496 
1497  // Contribution from Z'0 -> W^+ W^-.
1498  } else if (id1Abs == 24) {
1499  widNow = preFac * pow2(coupZpWW * cos2tW) * pow3(ps)
1500  * (1. + mr1*mr1 + mr2*mr2 + 10. * (mr1 + mr2 + mr1 * mr2));
1501  }
1502  }
1503 
1504  // When call for incoming flavour need to consider full mix.
1505  else {
1506 
1507  // Contributions from three fermion generations.
1508  if ( id1Abs < 7 || (id1Abs > 10 && id1Abs < 17) ) {
1509 
1510  // Couplings of gamma^*/Z^0/Z'^0 to final flavour
1511  double ef = couplingsPtr->ef(id1Abs);
1512  double af = couplingsPtr->af(id1Abs);
1513  double vf = couplingsPtr->vf(id1Abs);
1514  double apf = afZp[id1Abs];
1515  double vpf = vfZp[id1Abs];
1516 
1517  // Combine couplings with kinematical factors.
1518  double kinFacA = pow3(ps);
1519  double kinFacV = ps * (1. + 2. * mr1);
1520  double ef2 = ef * ef * kinFacV;
1521  double efvf = ef * vf * kinFacV;
1522  double vaf2 = vf * vf * kinFacV + af * af * kinFacA;
1523  double efvpf = ef * vpf * kinFacV;
1524  double vafvapf = vf * vpf * kinFacV + af * apf * kinFacA;
1525  double vapf2 = vpf * vpf * kinFacV + apf * apf * kinFacA;
1526 
1527  // Relative outwidths: combine instate, propagator and outstate.
1528  widNow = gamNorm * ef2 + gamZNorm * efvf + ZNorm * vaf2
1529  + gamZpNorm * efvpf + ZZpNorm * vafvapf + ZpNorm * vapf2;
1530  if (id1Abs < 7) widNow *= colQ;
1531 
1532  // Contribution from Z'0 -> W^+ W^-.
1533  } else if (id1Abs == 24) {
1534  widNow = ZpNorm * pow2(coupZpWW * cos2tW) * pow3(ps)
1535  * (1. + mr1*mr1 + mr2*mr2 + 10. * (mr1 + mr2 + mr1 * mr2));
1536  }
1537  }
1538 
1539 }
1540 
1541 //==========================================================================
1542 
1543 // The ResonanceWprime class.
1544 // Derived class for W'+- properties.
1545 
1546 //--------------------------------------------------------------------------
1547 
1548 // Initialize constants.
1549 
1550 void ResonanceWprime::initConstants() {
1551 
1552  // Locally stored properties and couplings.
1553  thetaWRat = 1. / (12. * couplingsPtr->sin2thetaW());
1554  cos2tW = couplingsPtr->cos2thetaW();
1555 
1556  // Axial and vector couplings of fermions.
1557  aqWp = settingsPtr->parm("Wprime:aq");
1558  vqWp = settingsPtr->parm("Wprime:vq");
1559  alWp = settingsPtr->parm("Wprime:al");
1560  vlWp = settingsPtr->parm("Wprime:vl");
1561 
1562  // Coupling for W' -> W Z.
1563  coupWpWZ = settingsPtr->parm("Wprime:coup2WZ");
1564 
1565 }
1566 
1567 //--------------------------------------------------------------------------
1568 
1569 // Calculate various common prefactors for the current mass.
1570 
1571 void ResonanceWprime::calcPreFac(bool) {
1572 
1573  // Common coupling factors.
1574  alpEM = couplingsPtr->alphaEM(mHat * mHat);
1575  alpS = couplingsPtr->alphaS(mHat * mHat);
1576  colQ = 3. * (1. + alpS / M_PI);
1577  preFac = alpEM * thetaWRat * mHat;
1578 
1579 }
1580 
1581 //--------------------------------------------------------------------------
1582 
1583 // Calculate width for currently considered channel.
1584 
1585 void ResonanceWprime::calcWidth(bool) {
1586 
1587  // Check that above threshold.
1588  if (ps == 0.) return;
1589 
1590  // Decay to quarks involves colour factor and CKM matrix.
1591  if (id1Abs > 0 && id1Abs < 7) widNow
1592  = preFac * ps * 0.5 * ((vqWp*vqWp + aqWp * aqWp)
1593  + 6. * (vqWp*vqWp - aqWp * aqWp) * sqrt(mr1 *mr2))
1594  * (1. - 0.5 * (mr1 + mr2) - 0.5 * pow2(mr1 - mr2))
1595  * colQ * couplingsPtr->V2CKMid(id1Abs, id2Abs);
1596 
1597  // Decay to leptons simpler.
1598  else if (id1Abs > 10 && id1Abs < 17) widNow
1599  = preFac * ps * 0.5 * ((vlWp*vqWp + alWp * aqWp)
1600  + 6. * (vlWp*vqWp - alWp * aqWp) * sqrt(mr1 *mr2))
1601  * (1. - 0.5 * (mr1 + mr2) - 0.5 * pow2(mr1 - mr2));
1602 
1603  // Decay to W^+- Z^0.
1604  else if (id1Abs == 24 && id2Abs == 23) widNow
1605  = preFac * 0.25 * pow2(coupWpWZ) * cos2tW * (mr1 / mr2) * pow3(ps)
1606  * (1. + mr1*mr1 + mr2*mr2 + 10. * (mr1 + mr2 + mr1 * mr2));
1607 
1608 }
1609 
1610 //==========================================================================
1611 
1612 // The ResonanceRhorizontal class.
1613 // Derived class for R^0 (horizontal gauge boson) properties.
1614 
1615 //--------------------------------------------------------------------------
1616 
1617 // Initialize constants.
1618 
1619 void ResonanceRhorizontal::initConstants() {
1620 
1621  // Locally stored properties and couplings.
1622  thetaWRat = 1. / (12. * couplingsPtr->sin2thetaW());
1623 
1624 }
1625 
1626 //--------------------------------------------------------------------------
1627 
1628 // Calculate various common prefactors for the current mass.
1629 
1630 void ResonanceRhorizontal::calcPreFac(bool) {
1631 
1632  // Common coupling factors.
1633  alpEM = couplingsPtr->alphaEM(mHat * mHat);
1634  alpS = couplingsPtr->alphaS(mHat * mHat);
1635  colQ = 3. * (1. + alpS / M_PI);
1636  preFac = alpEM * thetaWRat * mHat;
1637 
1638 }
1639 
1640 //--------------------------------------------------------------------------
1641 
1642 // Calculate width for currently considered channel.
1643 
1644 void ResonanceRhorizontal::calcWidth(bool) {
1645 
1646  // Check that above threshold.
1647  if (ps == 0.) return;
1648 
1649  // R -> f fbar. Colour factor for quarks.
1650  widNow = preFac * ps * (2. - mr1 - mr2 - pow2(mr1 - mr2));
1651  if (id1Abs < 9) widNow *= colQ;
1652 
1653 }
1654 
1655 //==========================================================================
1656 
1657 // The ResonanceExcited class.
1658 // Derived class for excited-fermion properties.
1659 
1660 //--------------------------------------------------------------------------
1661 
1662 // Initialize constants.
1663 
1664 void ResonanceExcited::initConstants() {
1665 
1666  // Locally stored properties and couplings.
1667  Lambda = settingsPtr->parm("ExcitedFermion:Lambda");
1668  coupF = settingsPtr->parm("ExcitedFermion:coupF");
1669  coupFprime = settingsPtr->parm("ExcitedFermion:coupFprime");
1670  coupFcol = settingsPtr->parm("ExcitedFermion:coupFcol");
1671  sin2tW = couplingsPtr->sin2thetaW();
1672  cos2tW = 1. - sin2tW;
1673 
1674 }
1675 
1676 //--------------------------------------------------------------------------
1677 
1678 // Calculate various common prefactors for the current mass.
1679 
1680 void ResonanceExcited::calcPreFac(bool) {
1681 
1682  // Common coupling factors.
1683  alpEM = couplingsPtr->alphaEM(mHat * mHat);
1684  alpS = couplingsPtr->alphaS(mHat * mHat);
1685  preFac = pow3(mHat) / pow2(Lambda);
1686 
1687 }
1688 
1689 //--------------------------------------------------------------------------
1690 
1691 // Calculate width for currently considered channel.
1692 
1693 void ResonanceExcited::calcWidth(bool) {
1694 
1695  // Check that above threshold.
1696  if (ps == 0.) return;
1697 
1698  // f^* -> f g.
1699  if (id1Abs == 21) widNow = preFac * alpS * pow2(coupFcol) / 3.;
1700 
1701  // f^* -> f gamma.
1702  else if (id1Abs == 22) {
1703  double chgI3 = (id2Abs%2 == 0) ? 0.5 : -0.5;
1704  double chgY = (id2Abs < 9) ? 1. / 6. : -0.5;
1705  double chg = chgI3 * coupF + chgY * coupFprime;
1706  widNow = preFac * alpEM * pow2(chg) / 4.;
1707  }
1708 
1709  // f^* -> f Z^0.
1710  else if (id1Abs == 23) {
1711  double chgI3 = (id2Abs%2 == 0) ? 0.5 : -0.5;
1712  double chgY = (id2Abs < 9) ? 1. / 6. : -0.5;
1713  double chg = chgI3 * cos2tW * coupF - chgY * sin2tW * coupFprime;
1714  widNow = preFac * (alpEM * pow2(chg) / (8. * sin2tW * cos2tW))
1715  * ps*ps * (2. + mr1);
1716  }
1717 
1718  // f^* -> f' W^+-.
1719  else if (id1Abs == 24) widNow = preFac * (alpEM * pow2(coupF)
1720  / (16. * sin2tW)) * ps*ps * (2. + mr1);
1721 
1722 }
1723 
1724 //==========================================================================
1725 
1726 // The ResonanceGraviton class.
1727 // Derived class for excited Graviton properties.
1728 
1729 //--------------------------------------------------------------------------
1730 
1731 // Initialize constants.
1732 
1733 void ResonanceGraviton::initConstants() {
1734 
1735  // SMinBulk = off/on, use universal coupling (kappaMG)
1736  // or individual (Gxx) between graviton and SM particles.
1737  eDsmbulk = settingsPtr->flag("ExtraDimensionsG*:SMinBulk");
1738  eDvlvl = false;
1739  if (eDsmbulk) eDvlvl = settingsPtr->flag("ExtraDimensionsG*:VLVL");
1740  kappaMG = settingsPtr->parm("ExtraDimensionsG*:kappaMG");
1741  for (int i = 0; i < 27; ++i) eDcoupling[i] = 0.;
1742  double tmp_coup = settingsPtr->parm("ExtraDimensionsG*:Gqq");
1743  for (int i = 1; i <= 4; ++i) eDcoupling[i] = tmp_coup;
1744  eDcoupling[5] = settingsPtr->parm("ExtraDimensionsG*:Gbb");
1745  eDcoupling[6] = settingsPtr->parm("ExtraDimensionsG*:Gtt");
1746  tmp_coup = settingsPtr->parm("ExtraDimensionsG*:Gll");
1747  for (int i = 11; i <= 16; ++i) eDcoupling[i] = tmp_coup;
1748  eDcoupling[21] = settingsPtr->parm("ExtraDimensionsG*:Ggg");
1749  eDcoupling[22] = settingsPtr->parm("ExtraDimensionsG*:Ggmgm");
1750  eDcoupling[23] = settingsPtr->parm("ExtraDimensionsG*:GZZ");
1751  eDcoupling[24] = settingsPtr->parm("ExtraDimensionsG*:GWW");
1752  eDcoupling[25] = settingsPtr->parm("ExtraDimensionsG*:Ghh");
1753 
1754 }
1755 
1756 //--------------------------------------------------------------------------
1757 
1758 // Calculate various common prefactors for the current mass.
1759 
1760 void ResonanceGraviton::calcPreFac(bool) {
1761 
1762  // Common coupling factors.
1763  alpS = couplingsPtr->alphaS(mHat * mHat);
1764  colQ = 3. * (1. + alpS / M_PI);
1765  preFac = mHat / M_PI;
1766 
1767 }
1768 
1769 //--------------------------------------------------------------------------
1770 
1771 // Calculate width for currently considered channel.
1772 
1773 void ResonanceGraviton::calcWidth(bool) {
1774 
1775  // Check that above threshold.
1776  if (ps == 0.) return;
1777 
1778  // Widths to fermion pairs.
1779  if (id1Abs < 19) {
1780  widNow = preFac * pow3(ps) * (1. + 8. * mr1 / 3.) / 320.;
1781  if (id1Abs < 9) widNow *= colQ;
1782 
1783  // Widths to gluon and photon pair.
1784  } else if (id1Abs == 21) {
1785  widNow = preFac / 20.;
1786  } else if (id1Abs == 22) {
1787  widNow = preFac / 160.;
1788 
1789  // Widths to Z0 Z0 and W+ W- pair.
1790  } else if (id1Abs == 23 || id1Abs == 24) {
1791  // Longitudinal W/Z only.
1792  if (eDvlvl) {
1793  widNow = preFac * pow(ps,5) / 480.;
1794  // Transverse W/Z contributions as well.
1795  } else {
1796  widNow = preFac * ps * (13. / 12. + 14. * mr1 / 3. + 4. * mr1 * mr1)
1797  / 80.;
1798  }
1799  if (id1Abs == 23) widNow *= 0.5;
1800 
1801  // Widths to h h pair.
1802  } else if (id1Abs == 25) {
1803  widNow = preFac * pow(ps,5) / 960.;
1804  }
1805 
1806  // RS graviton coupling
1807  if (eDsmbulk) widNow *= 2. * pow2(eDcoupling[min( id1Abs, 26)] * mHat);
1808  else widNow *= pow2(kappaMG * mHat / mRes);
1809 
1810 }
1811 
1812 //==========================================================================
1813 
1814 // The ResonanceKKgluon class.
1815 // Derived class for excited kk-gluon properties.
1816 
1817 //--------------------------------------------------------------------------
1818 
1819 // Initialize constants.
1820 
1821 void ResonanceKKgluon::initConstants() {
1822 
1823  // KK-gluon gv/ga couplings and interference.
1824  for (int i = 0; i < 10; ++i) { eDgv[i] = 0.; eDga[i] = 0.; }
1825  double tmp_gL = settingsPtr->parm("ExtraDimensionsG*:KKgqL");
1826  double tmp_gR = settingsPtr->parm("ExtraDimensionsG*:KKgqR");
1827  for (int i = 1; i <= 4; ++i) {
1828  eDgv[i] = 0.5 * (tmp_gL + tmp_gR);
1829  eDga[i] = 0.5 * (tmp_gL - tmp_gR);
1830  }
1831  tmp_gL = settingsPtr->parm("ExtraDimensionsG*:KKgbL");
1832  tmp_gR = settingsPtr->parm("ExtraDimensionsG*:KKgbR");
1833  eDgv[5] = 0.5 * (tmp_gL + tmp_gR); eDga[5] = 0.5 * (tmp_gL - tmp_gR);
1834  tmp_gL = settingsPtr->parm("ExtraDimensionsG*:KKgtL");
1835  tmp_gR = settingsPtr->parm("ExtraDimensionsG*:KKgtR");
1836  eDgv[6] = 0.5 * (tmp_gL + tmp_gR); eDga[6] = 0.5 * (tmp_gL - tmp_gR);
1837  interfMode = settingsPtr->mode("ExtraDimensionsG*:KKintMode");
1838 
1839 }
1840 
1841 //--------------------------------------------------------------------------
1842 
1843 // Calculate various common prefactors for the current mass.
1844 
1845 void ResonanceKKgluon::calcPreFac(bool calledFromInit) {
1846 
1847  // Common coupling factors.
1848  alpS = couplingsPtr->alphaS(mHat * mHat);
1849  preFac = alpS * mHat / 6;
1850 
1851  // When call for incoming flavour need to consider g*/gKK mix.
1852  if (!calledFromInit) {
1853  // Calculate prefactors for g/interference/gKK terms.
1854  int idInFlavAbs = abs(idInFlav);
1855  double sH = mHat * mHat;
1856  normSM = 1;
1857  normInt = 2. * eDgv[min(idInFlavAbs, 9)] * sH * (sH - m2Res)
1858  / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
1859  normKK = ( pow2(eDgv[min(idInFlavAbs, 9)])
1860  + pow2(eDga[min(idInFlavAbs, 9)]) ) * sH * sH
1861  / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
1862 
1863  // Optionally only keep g* or gKK term.
1864  if (interfMode == 1) {normInt = 0.; normKK = 0.;}
1865  if (interfMode == 2) {normSM = 0.; normInt = 0.; normKK = 1.;}
1866  }
1867 
1868 }
1869 
1870 //--------------------------------------------------------------------------
1871 
1872 // Calculate width for currently considered channel.
1873 
1874 void ResonanceKKgluon::calcWidth(bool calledFromInit) {
1875 
1876  // Check that above threshold.
1877  if (ps == 0.) return;
1878 
1879  // Widths to quark pairs.
1880  if (id1Abs > 9) return;
1881 
1882  if (calledFromInit) {
1883  widNow = preFac * ps * (pow2(eDgv[min(id1Abs, 9)]) * (1. + 2.*mr1)
1884  + pow2(eDga[min(id1Abs, 9)]) * (1. - 4.*mr1) );
1885  } else {
1886  // Relative outwidths: combine instate, propagator and outstate.
1887  widNow = normSM * ps * (1. + 2. * mr1)
1888  + normInt * ps * eDgv[min(id1Abs, 9)] * (1. + 2. * mr1)
1889  + normKK * ps * (pow2(eDgv[min(id1Abs, 9)]) * (1. + 2.*mr1)
1890  + pow2(eDga[min(id1Abs, 9)]) * (1. - 4.*mr1) );
1891  widNow *= preFac;
1892  }
1893 
1894 }
1895 
1896 //==========================================================================
1897 
1898 // The ResonanceLeptoquark class.
1899 // Derived class for leptoquark properties.
1900 
1901 //--------------------------------------------------------------------------
1902 
1903 // Initialize constants.
1904 
1905 void ResonanceLeptoquark::initConstants() {
1906 
1907  // Locally stored properties and couplings.
1908  kCoup = settingsPtr->parm("LeptoQuark:kCoup");
1909 
1910  // Check that flavour info in decay channel is correctly set.
1911  int id1Now = particlePtr->channel(0).product(0);
1912  int id2Now = particlePtr->channel(0).product(1);
1913  if (id1Now < 1 || id1Now > 6) {
1914  infoPtr->errorMsg("Error in ResonanceLeptoquark::init:"
1915  " unallowed input quark flavour reset to u");
1916  id1Now = 2;
1917  particlePtr->channel(0).product(0, id1Now);
1918  }
1919  if (abs(id2Now) < 11 || abs(id2Now) > 16) {
1920  infoPtr->errorMsg("Error in ResonanceLeptoquark::init:"
1921  " unallowed input lepton flavour reset to e-");
1922  id2Now = 11;
1923  particlePtr->channel(0).product(1, id2Now);
1924  }
1925 
1926  // Set/overwrite charge and name of particle.
1927  bool changed = particlePtr->hasChanged();
1928  int chargeLQ = particleDataPtr->chargeType(id1Now)
1929  + particleDataPtr->chargeType(id2Now);
1930  particlePtr->setChargeType(chargeLQ);
1931  string nameLQ = "LQ_" + particleDataPtr->name(id1Now) + ","
1932  + particleDataPtr->name(id2Now);
1933  particlePtr->setNames(nameLQ, nameLQ + "bar");
1934  if (!changed) particlePtr->setHasChanged(false);
1935 
1936 }
1937 
1938 //--------------------------------------------------------------------------
1939 
1940 // Calculate various common prefactors for the current mass.
1941 
1942 void ResonanceLeptoquark::calcPreFac(bool) {
1943 
1944  // Common coupling factors.
1945  alpEM = couplingsPtr->alphaEM(mHat * mHat);
1946  preFac = 0.25 * alpEM * kCoup * mHat;
1947 
1948 }
1949 
1950 //--------------------------------------------------------------------------
1951 
1952 // Calculate width for currently considered channel.
1953 
1954 void ResonanceLeptoquark::calcWidth(bool) {
1955 
1956  // Check that above threshold.
1957  if (ps == 0.) return;
1958 
1959  // Width into lepton plus quark.
1960  if (id1Abs > 10 && id1Abs < 17 && id2Abs < 7) widNow = preFac * pow3(ps);
1961 
1962 }
1963 
1964 //==========================================================================
1965 
1966 // The ResonanceNuRight class.
1967 // Derived class for righthanded Majorana neutrino properties.
1968 
1969 //--------------------------------------------------------------------------
1970 
1971 // Initialize constants.
1972 
1973 void ResonanceNuRight::initConstants() {
1974 
1975  // Locally stored properties and couplings: righthanded W mass.
1976  thetaWRat = 1. / (768. * M_PI * pow2(couplingsPtr->sin2thetaW()));
1977  mWR = particleDataPtr->m0(9900024);
1978 
1979 }
1980 
1981 //--------------------------------------------------------------------------
1982 
1983 // Calculate various common prefactors for the current mass.
1984 
1985 void ResonanceNuRight::calcPreFac(bool) {
1986 
1987  // Common coupling factors.
1988  alpEM = couplingsPtr->alphaEM(mHat * mHat);
1989  alpS = couplingsPtr->alphaS(mHat * mHat);
1990  colQ = 3. * (1. + alpS / M_PI);
1991  preFac = pow2(alpEM) * thetaWRat * pow5(mHat) / pow4(max(mHat, mWR));
1992 
1993 }
1994 
1995 //--------------------------------------------------------------------------
1996 
1997 // Calculate width for currently considered channel.
1998 
1999 void ResonanceNuRight::calcWidth(bool) {
2000 
2001  // Check that above threshold.
2002  if (mHat < mf1 + mf2 + mf3 + MASSMARGIN) return;
2003 
2004  // Coupling part of widths to l- q qbar', l- l'+ nu_lR' and c.c.
2005  widNow = (id2Abs < 9 && id3Abs < 9)
2006  ? preFac * colQ * couplingsPtr->V2CKMid(id2, id3) : preFac;
2007 
2008  // Phase space corrections in decay. Must have y < 1.
2009  double x = (mf1 + mf2 + mf3) / mHat;
2010  double x2 = x * x;
2011  double fx = 1. - 8. * x2 + 8. * pow3(x2) - pow4(x2)
2012  - 24. * pow2(x2) * log(x);
2013  double y = min( 0.999, pow2(mHat / mWR) );
2014  double fy = ( 12. * (1. - y) * log(1. - y) + 12. * y - 6. * y*y
2015  - 2.* pow3(y) ) / pow4(y);
2016  widNow *= fx * fy;
2017 
2018 }
2019 
2020 //==========================================================================
2021 
2022 // The ResonanceZRight class.
2023 // Derived class for Z_R^0 properties.
2024 
2025 //--------------------------------------------------------------------------
2026 
2027 // Initialize constants.
2028 
2029 void ResonanceZRight::initConstants() {
2030 
2031  // Locally stored properties and couplings: righthanded W mass.
2032  sin2tW = couplingsPtr->sin2thetaW();
2033  thetaWRat = 1. / (48. * sin2tW * (1. - sin2tW) * (1. - 2. * sin2tW) );
2034 
2035 }
2036 
2037 //--------------------------------------------------------------------------
2038 
2039 // Calculate various common prefactors for the current mass.
2040 
2041 void ResonanceZRight::calcPreFac(bool) {
2042 
2043  // Common coupling factors.
2044  alpEM = couplingsPtr->alphaEM(mHat * mHat);
2045  alpS = couplingsPtr->alphaS(mHat * mHat);
2046  colQ = 3. * (1. + alpS / M_PI);
2047  preFac = alpEM * thetaWRat * mHat;
2048 
2049 }
2050 
2051 //--------------------------------------------------------------------------
2052 
2053 // Calculate width for currently considered channel.
2054 
2055 void ResonanceZRight::calcWidth(bool) {
2056 
2057  // Check that above threshold.
2058  if (ps == 0.) return;
2059 
2060  // Couplings to q qbar and l+ l-.
2061  double vf = 0.;
2062  double af = 0.;
2063  double symMaj = 1.;
2064  if (id1Abs < 9 && id1Abs%2 == 1) {
2065  af = -1. + 2. * sin2tW;
2066  vf = -1. + 4. * sin2tW / 3.;
2067  } else if (id1Abs < 9) {
2068  af = 1. - 2. * sin2tW;
2069  vf = 1. - 8. * sin2tW / 3.;
2070  } else if (id1Abs < 19 && id1Abs%2 == 1) {
2071  af = -1. + 2. * sin2tW;
2072  vf = -1. + 4. * sin2tW;
2073 
2074  // Couplings to nu_L nu_Lbar and nu_R nu_Rbar, both assumed Majoranas.
2075  } else if (id1Abs < 19) {
2076  af = -2. * sin2tW;
2077  vf = 0.;
2078  symMaj = 0.5;
2079  } else {
2080  af = 2. * (1. - sin2tW);
2081  vf = 0.;
2082  symMaj = 0.5;
2083  }
2084 
2085  // Width expression, including phase space and colour factor.
2086  widNow = preFac * (vf*vf * (1. + 2. * mr1) + af*af * ps*ps) * ps
2087  * symMaj;
2088  if (id1Abs < 9) widNow *= colQ;
2089 
2090 }
2091 
2092 //==========================================================================
2093 
2094 // The ResonanceWRight class.
2095 // Derived class for W_R+- properties.
2096 
2097 //--------------------------------------------------------------------------
2098 
2099 // Initialize constants.
2100 
2101 void ResonanceWRight::initConstants() {
2102 
2103  // Locally stored properties and couplings.
2104  thetaWRat = 1. / (12. * couplingsPtr->sin2thetaW());
2105 
2106 }
2107 
2108 //--------------------------------------------------------------------------
2109 
2110 // Calculate various common prefactors for the current mass.
2111 
2112 void ResonanceWRight::calcPreFac(bool) {
2113 
2114  // Common coupling factors.
2115  alpEM = couplingsPtr->alphaEM(mHat * mHat);
2116  alpS = couplingsPtr->alphaS(mHat * mHat);
2117  colQ = 3. * (1. + alpS / M_PI);
2118  preFac = alpEM * thetaWRat * mHat;
2119 
2120 }
2121 
2122 //--------------------------------------------------------------------------
2123 
2124 // Calculate width for currently considered channel.
2125 
2126 void ResonanceWRight::calcWidth(bool) {
2127 
2128  // Check that above threshold.
2129  if (ps == 0.) return;
2130 
2131  // Combine kinematics with colour factor and CKM couplings.
2132  widNow = preFac * (1. - 0.5 * (mr1 + mr2) - 0.5 * pow2(mr1 - mr2))
2133  * ps;
2134  if (id1Abs < 9) widNow *= colQ * couplingsPtr->V2CKMid(id1Abs, id2Abs);
2135 
2136 }
2137 
2138 //==========================================================================
2139 
2140 // The ResonanceHchgchgLeft class.
2141 // Derived class for H++/H-- (left) properties.
2142 
2143 //--------------------------------------------------------------------------
2144 
2145 // Initialize constants.
2146 
2147 void ResonanceHchgchgLeft::initConstants() {
2148 
2149  // Read in Yukawa matrix for couplings to a lepton pair.
2150  yukawa[1][1] = settingsPtr->parm("LeftRightSymmmetry:coupHee");
2151  yukawa[2][1] = settingsPtr->parm("LeftRightSymmmetry:coupHmue");
2152  yukawa[2][2] = settingsPtr->parm("LeftRightSymmmetry:coupHmumu");
2153  yukawa[3][1] = settingsPtr->parm("LeftRightSymmmetry:coupHtaue");
2154  yukawa[3][2] = settingsPtr->parm("LeftRightSymmmetry:coupHtaumu");
2155  yukawa[3][3] = settingsPtr->parm("LeftRightSymmmetry:coupHtautau");
2156 
2157  // Locally stored properties and couplings.
2158  gL = settingsPtr->parm("LeftRightSymmmetry:gL");
2159  vL = settingsPtr->parm("LeftRightSymmmetry:vL");
2160  mW = particleDataPtr->m0(24);
2161 
2162 }
2163 
2164 //--------------------------------------------------------------------------
2165 
2166 // Calculate various common prefactors for the current mass.
2167 
2168 void ResonanceHchgchgLeft::calcPreFac(bool) {
2169 
2170  // Common coupling factors.
2171  preFac = mHat / (8. * M_PI);
2172 
2173 }
2174 
2175 //--------------------------------------------------------------------------
2176 
2177 // Calculate width for currently considered channel.
2178 
2179 void ResonanceHchgchgLeft::calcWidth(bool) {
2180 
2181  // Check that above threshold.
2182  if (ps == 0.) return;
2183 
2184  // H++-- width to a pair of leptons. Combinatorial factor of 2.
2185  if (id1Abs < 17 && id2Abs < 17) {
2186  widNow = preFac * pow2(yukawa[(id1Abs-9)/2][(id2Abs-9)/2]) * ps;
2187  if (id2Abs != id1Abs) widNow *= 2.;
2188  }
2189 
2190  // H++-- width to a pair of lefthanded W's.
2191  else if (id1Abs == 24 && id2Abs == 24)
2192  widNow = preFac * 0.5 * pow2(gL*gL * vL / mW)
2193  * (3. * mr1 + 0.25 / mr1 - 1.) * ps;
2194 
2195 }
2196 
2197 //==========================================================================
2198 
2199 // The ResonanceHchgchgRight class.
2200 // Derived class for H++/H-- (right) properties.
2201 
2202 //--------------------------------------------------------------------------
2203 
2204 // Initialize constants.
2205 
2206 void ResonanceHchgchgRight::initConstants() {
2207 
2208  // Read in Yukawa matrix for couplings to a lepton pair.
2209  yukawa[1][1] = settingsPtr->parm("LeftRightSymmmetry:coupHee");
2210  yukawa[2][1] = settingsPtr->parm("LeftRightSymmmetry:coupHmue");
2211  yukawa[2][2] = settingsPtr->parm("LeftRightSymmmetry:coupHmumu");
2212  yukawa[3][1] = settingsPtr->parm("LeftRightSymmmetry:coupHtaue");
2213  yukawa[3][2] = settingsPtr->parm("LeftRightSymmmetry:coupHtaumu");
2214  yukawa[3][3] = settingsPtr->parm("LeftRightSymmmetry:coupHtautau");
2215 
2216  // Locally stored properties and couplings.
2217  idWR = 9000024;
2218  gR = settingsPtr->parm("LeftRightSymmmetry:gR");
2219 
2220 }
2221 
2222 //--------------------------------------------------------------------------
2223 
2224 // Calculate various common prefactors for the current mass.
2225 
2226 void ResonanceHchgchgRight::calcPreFac(bool) {
2227 
2228  // Common coupling factors.
2229  preFac = mHat / (8. * M_PI);
2230 
2231 }
2232 
2233 //--------------------------------------------------------------------------
2234 
2235 // Calculate width for currently considered channel.
2236 
2237 void ResonanceHchgchgRight::calcWidth(bool) {
2238 
2239  // Check that above threshold.
2240  if (ps == 0.) return;
2241 
2242  // H++-- width to a pair of leptons. Combinatorial factor of 2.
2243  if (id1Abs < 17 && id2Abs < 17) {
2244  widNow = preFac * pow2(yukawa[(id1Abs-9)/2][(id2Abs-9)/2]) * ps;
2245  if (id2Abs != id1Abs) widNow *= 2.;
2246  }
2247 
2248  // H++-- width to a pair of lefthanded W's.
2249  else if (id1Abs == idWR && id2Abs == idWR)
2250  widNow = preFac * pow2(yukawa[(id1Abs-9)/2][(id2Abs-9)/2]) * ps;
2251 
2252 }
2253 
2254 //==========================================================================
2255 
2256 } // end namespace Pythia8
Generated by   doxygen 1.8.5