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