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