PYTHIA8/pythia8130/src/ResonanceWidths.cxx

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