PYTHIA8/pythia8130/src/StandardModel.cxx

   1 // StandardModel.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 the AlphaStrong class.
   7
   8 #include "StandardModel.h"
   9
  10 namespace Pythia8 {
  11
  12 //**************************************************************************
  13
  14 // The AlphaStrong class.
  15
  16 //*********
  17
  18 // Constants: could be changed here if desired, but normally should not.
  19 // These are of technical nature, as described for each.
  20
  21 // Number of iterations to determine Lambda from given alpha_s.
  22 const int AlphaStrong::NITER           = 10;
  23
  24 // Always evaluate running alpha_s above Lambda3 to avoid disaster.
  25 // Safety margin picked to freeze roughly for alpha_s = 10.
  26 const double AlphaStrong::SAFETYMARGIN1 = 1.07;
  27 const double AlphaStrong::SAFETYMARGIN2 = 1.33;
  28
  29 //*********
  30
  31 // Initialize alpha_strong calculation by finding Lambda values etc.
  32
  33 void AlphaStrong::init( double valueIn, int orderIn) {
  34
  35   // Order of alpha_s evaluation. Charm, bottom and Z0 masses.
  36   // Pick defaults if ParticleDataTable not properly initialized.
  37   valueRef = valueIn;
  38   order    = max( 0, min( 2, orderIn ) );
  39   mc                           = ParticleDataTable::m0(4);
  40   if (mc < 1. || mc > 2.) mc   = 1.5;
  41   mb                           = ParticleDataTable::m0(5);
  42   if (mb < 3. || mb > 6.) mb   = 4.8;
  43   mZ                           = ParticleDataTable::m0(23);
  44   if (mZ < 90. || mZ > 92.) mZ = 91.188;
  45
  46   // Fix alpha_s.
  47   if (order == 0) {
  48     Lambda3Save = Lambda4Save = Lambda5Save = scale2Min = 0.;
  49
  50   // First order alpha_s: match at flavour thresholds.
  51   } else if (order == 1) {
  52     Lambda5Save = mZ * exp( -6. * M_PI / (23. * valueRef) );
  53     Lambda4Save = Lambda5Save * pow(mb/Lambda5Save, 2./25.);
  54     Lambda3Save = Lambda4Save * pow(mc/Lambda4Save, 2./27.);
  55     scale2Min   = pow2(SAFETYMARGIN1 * Lambda3Save);
  56
  57   // Second order alpha_s: iterative match at flavour thresholds.
  58   } else {
  59     double b15 = 348. / 529.;
  60     double b14 = 462. / 625.;
  61     double b13 = 64. / 81.;
  62     double b25 = 224687. / 242208.;
  63     double b24 = 548575. / 426888.;
  64     double b23 = 938709. / 663552.;
  65     double logScale, loglogScale, correction, valueIter;
  66
  67     // Find Lambda_5 at m_Z.
  68     Lambda5Save = mZ * exp( -6. * M_PI / (23. * valueRef) );
  69     for (int iter = 0; iter < NITER; ++iter) {
  70       logScale    = 2. * log(mZ/Lambda5Save);
  71       loglogScale = log(logScale);
  72       correction  = 1. - b15 * loglogScale / logScale
  73         + pow2(b15 / logScale) * (pow2(loglogScale - 0.5) + b25 - 1.25);
  74       valueIter   = valueRef / correction;
  75       Lambda5Save = mZ * exp( -6. * M_PI / (23. * valueIter) );
  76     }
  77
  78     // Find Lambda_4 at m_b.
  79     double logScaleB    = 2. * log(mb/Lambda5Save);
  80     double loglogScaleB = log(logScaleB);
  81     double valueB       = 12. * M_PI / (23. * logScaleB)
  82       * (1. - b15 * loglogScaleB / logScaleB
  83         + pow2(b15 / logScaleB) * (pow2(loglogScaleB - 0.5) + b25- 1.25) );
  84     Lambda4Save         = Lambda5Save;
  85     for (int iter = 0; iter < NITER; ++iter) {
  86       logScale    = 2. * log(mb/Lambda4Save);
  87       loglogScale = log(logScale);
  88       correction  = 1. - b14 * loglogScale / logScale
  89         + pow2(b14 / logScale) * (pow2(loglogScale - 0.5) + b24 - 1.25);
  90       valueIter   = valueB / correction;
  91       Lambda4Save = mb * exp( -6. * M_PI / (25. * valueIter) );
  92     }
  93
  94     // Find Lambda_3 at m_c.
  95     double logScaleC    = 2. * log(mc/Lambda4Save);
  96     double loglogScaleC = log(logScaleC);
  97     double valueC       = 12. * M_PI / (25. * logScaleC)
  98       * (1. - b14 * loglogScaleC / logScaleC
  99         + pow2(b14 / logScaleC) * (pow2(loglogScaleC - 0.5) + b24 - 1.25) );
 100     Lambda3Save = Lambda4Save;
 101     for (int iter = 0; iter < NITER; ++iter) {
 102       logScale    = 2. * log(mc/Lambda3Save);
 103       loglogScale = log(logScale);
 104       correction  = 1. - b13 * loglogScale / logScale
 105         + pow2(b13 / logScale) * (pow2(loglogScale - 0.5) + b23 - 1.25);
 106       valueIter   = valueC / correction;
 107       Lambda3Save = mc * exp( -6. * M_PI / (27. * valueIter) );
 108     }
 109     scale2Min     = pow2(SAFETYMARGIN2 * Lambda3Save);
 110   }
 111
 112   // Save squares of mass and Lambda values as well.
 113   mc2          = pow2(mc);
 114   mb2          = pow2(mb);
 115   Lambda3Save2 = pow2(Lambda3Save);
 116   Lambda4Save2 = pow2(Lambda4Save);
 117   Lambda5Save2 = pow2(Lambda5Save);
 118   valueNow     = valueIn;
 119   scale2Now    = mZ * mZ;
 120   isInit = true;
 121
 122 }
 123
 124 //*********
 125
 126 // Calculate alpha_s value
 127
 128 double AlphaStrong::alphaS( double scale2) {
 129
 130   // Check for initialization and ensure minimal scale2 value.
 131   if (!isInit) return 0.;
 132   if (scale2 < scale2Min) scale2 = scale2Min;
 133
 134   // If equal to old scale then same answer.
 135   if (scale2 == scale2Now && (order < 2 || lastCallToFull)) return valueNow;
 136   scale2Now      = scale2;
 137   lastCallToFull = true;
 138
 139   // Fix alpha_s.
 140   if (order == 0) {
 141     valueNow = valueRef;
 142
 143   // First order alpha_s: differs by mass region.
 144   } else if (order == 1) {
 145     if (scale2 > mb2)
 146          valueNow = 12. * M_PI / (23. * log(scale2/Lambda5Save2));
 147     else if (scale2 > mc2)
 148          valueNow = 12. * M_PI / (25. * log(scale2/Lambda4Save2));
 149     else valueNow = 12. * M_PI / (27. * log(scale2/Lambda3Save2));
 150
 151   // Second order alpha_s: differs by mass region.
 152   } else {
 153     double Lambda2, b0, b1, b2;
 154     if (scale2 > mb2) {
 155       Lambda2 = Lambda5Save2;
 156       b0      = 23.;
 157       b1      = 348. / 529.;
 158       b2      = 224687. / 242208.;
 159     } else if (scale2 > mc2) {
 160       Lambda2 = Lambda4Save2;
 161       b0      = 25.;
 162       b1      = 462. / 625.;
 163       b2      = 548575. / 426888.;
 164     } else {
 165       Lambda2 = Lambda3Save2;
 166       b0      = 27.;
 167       b1      = 64. / 81.;
 168       b2      = 938709. / 663552.;
 169     }
 170     double logScale    = log(scale2/Lambda2);
 171     double loglogScale = log(logScale);
 172     valueNow = 12. * M_PI / (b0 * logScale)
 173       * ( 1. - b1 * loglogScale / logScale
 174         + pow2(b1 / logScale) * (pow2(loglogScale - 0.5) + b2 - 1.25) );
 175   }
 176
 177   // Done.
 178   return valueNow;
 179
 180 }
 181
 182 //*********
 183
 184 // Calculate alpha_s value, but only use up to first-order piece.
 185 // (To be combined with alphaS2OrdCorr.)
 186
 187 double  AlphaStrong::alphaS1Ord( double scale2) {
 188
 189   // Check for initialization and ensure minimal scale2 value.
 190   if (!isInit) return 0.;
 191   if (scale2 < scale2Min) scale2 = scale2Min;
 192
 193   // If equal to old scale then same answer.
 194   if (scale2 == scale2Now && (order < 2 || !lastCallToFull)) return valueNow;
 195   scale2Now      = scale2;
 196   lastCallToFull = false;
 197
 198   // Fix alpha_S.
 199   if (order == 0) {
 200     valueNow = valueRef;
 201
 202   // First/second order alpha_s: differs by mass region.
 203   } else {
 204     if (scale2 > mb2)
 205          valueNow = 12. * M_PI / (23. * log(scale2/Lambda5Save2));
 206     else if (scale2 > mc2)
 207          valueNow = 12. * M_PI / (25. * log(scale2/Lambda4Save2));
 208     else valueNow = 12. * M_PI / (27. * log(scale2/Lambda3Save2));
 209   }
 210
 211   // Done.
 212   return valueNow;
 213 }
 214
 215 //*********
 216
 217 // Calculates the second-order extra factor in alpha_s.
 218 // (To be combined with alphaS1Ord.)
 219
 220 double AlphaStrong::alphaS2OrdCorr( double scale2) {
 221
 222   // Check for initialization and ensure minimal scale2 value.
 223   if (!isInit) return 1.;
 224   if (scale2 < scale2Min) scale2 = scale2Min;
 225
 226   // Only meaningful for second order calculations.
 227   if (order < 2) return 1.;
 228
 229   // Second order correction term: differs by mass region.
 230   double Lambda2, b0, b1, b2;
 231   if (scale2 > mb2) {
 232     Lambda2 = Lambda5Save2;
 233     b0      = 23.;
 234     b1      = 348. / 529.;
 235     b2      = 224687. / 242208.;
 236   } else if (scale2 > mc2) {
 237     Lambda2 = Lambda4Save2;
 238     b0      = 25.;
 239     b1      = 462. / 625.;
 240     b2      = 548575. / 426888.;
 241   } else {
 242     Lambda2 = Lambda3Save2;
 243     b0      = 27.;
 244     b1      = 64. / 81.;
 245     b2      = 938709. / 663552.;
 246   }
 247   double logScale    = log(scale2/Lambda2);
 248   double loglogScale = log(logScale);
 249   return ( 1. - b1 * loglogScale / logScale
 250     + pow2(b1 / logScale) * (pow2(loglogScale - 0.5) + b2 - 1.25) );
 251
 252 }
 253
 254 //**************************************************************************
 255
 256 // The AlphaEM class.
 257
 258 //*********
 259
 260 // Definitions of static variables.
 261 double AlphaEM::alpEM0  = 0.00729735;
 262 double AlphaEM::alpEMmZ = 0.00781751;
 263 double AlphaEM::mZ2     = 8315.;
 264 // Effective thresholds for electron, muon, light quarks, tau+c, b.
 265 double AlphaEM::Q2step[5]    = {0.26e-6, 0.011, 0.25, 3.5, 90.};
 266 // Running coefficients are sum charge2 / 3 pi in pure QED, here slightly
 267 // enhanced for quarks to approximately account for QCD corrections.
 268 double AlphaEM::bRun[5]      = {0.1061, 0.2122, 0.460, 0.700, 0.725};
 269 double AlphaEM::alpEMstep[5] = {};
 270
 271 //*********
 272
 273 // Initialize alpha_EM calculation.
 274
 275 void AlphaEM::initStatic() {
 276
 277   // Read in alpha_EM value at 0 and m_Z, and mass of Z.
 278   alpEM0    = Settings::parm("StandardModel:alphaEM0");
 279   alpEMmZ   = Settings::parm("StandardModel:alphaEMmZ");
 280   double mZ = ParticleDataTable::m0(23);
 281   mZ2       = mZ * mZ;
 282
 283   // AlphaEM values at matching scales and matching b value.
 284
 285   // Step down from mZ to tau/charm threshold.
 286   alpEMstep[4] = alpEMmZ / ( 1. + alpEMmZ * bRun[4]
 287     * log(mZ2 / Q2step[4]) );
 288   alpEMstep[3] = alpEMstep[4] / ( 1. - alpEMstep[4] * bRun[3]
 289     * log(Q2step[3] / Q2step[4]) );
 290
 291   // Step up from me to light-quark threshold.
 292   alpEMstep[0] = alpEM0;
 293   alpEMstep[1] = alpEMstep[0] / ( 1. - alpEMstep[0] * bRun[0]
 294     * log(Q2step[1] / Q2step[0]) );
 295   alpEMstep[2] = alpEMstep[1] / ( 1. - alpEMstep[1] * bRun[1]
 296     * log(Q2step[2] / Q2step[1]) );
 297
 298   // Fit b in range between light-quark and tau/charm to join smoothly.
 299   bRun[2] = (1./alpEMstep[3] - 1./alpEMstep[2])
 300     / log(Q2step[2] / Q2step[3]);
 301
 302 }
 303
 304 //*********
 305
 306 // Calculate alpha_EM value
 307
 308 double AlphaEM::alphaEM( double scale2) {
 309
 310   // Fix alphaEM; for order = -1 fixed at m_Z.
 311   if (order == 0)  return alpEM0;
 312   if (order <  0)  return alpEMmZ;
 313
 314   // Running alphaEM.
 315   for (int i = 4; i >= 0; --i) if (scale2 > Q2step[i])
 316     return alpEMstep[i] / (1. - bRun[i] * alpEMstep[i]
 317       * log(scale2 / Q2step[i]) );
 318   return alpEM0;
 319
 320 }
 321
 322 //**************************************************************************
 323
 324 // The CoupEW class.
 325
 326 //*********
 327
 328 // Definitions of static variables.
 329
 330 double CoupEW::s2tW = 0.232;
 331 double CoupEW::c2tW = 0.768;
 332 double CoupEW::s2tWbar = 0.232;
 333 double CoupEW::efSave[20] = { 0., -1./3., 2./3., -1./3., 2./3., -1./3.,
 334   2./3., -1./3., 2./3., 0., 0., -1., 0., -1., 0., -1., 0., -1., 0., 0.};
 335 double CoupEW::vfSave[20] = { };
 336 double CoupEW::afSave[20] = { 0., -1., 1., -1., 1., -1., 1., -1., 1.,
 337   0., 0., -1., 1., -1., 1., -1., 1., -1., 1., 0.};
 338 double CoupEW::lfSave[20] = { };
 339 double CoupEW::rfSave[20] = { };
 340 double CoupEW::ef2Save[20] = { };
 341 double CoupEW::vf2Save[20] = { };
 342 double CoupEW::af2Save[20] = { };
 343 double CoupEW::efvfSave[20] = { };
 344 double CoupEW::vf2af2Save[20] = { };
 345
 346 //*********
 347
 348 // Initialize electroweak mixing angle and couplings.
 349
 350 void CoupEW::initStatic() {
 351
 352   // Read in electroweak mixing angle.
 353   s2tW = Settings::parm("StandardModel:sin2thetaW");
 354   c2tW = 1. - s2tW;
 355   s2tWbar = Settings::parm("StandardModel:sin2thetaWbar");
 356
 357   // Initialize electroweak couplings.
 358   for (int i = 0; i < 20; ++i) {
 359     vfSave[i]  = afSave[i] - 4. * s2tWbar * efSave[i];
 360     lfSave[i]  = afSave[i] - 2. * s2tWbar * efSave[i];
 361     rfSave[i]  =           - 2. * s2tWbar * efSave[i];
 362     ef2Save[i] = pow2(efSave[i]);
 363     vf2Save[i] = pow2(vfSave[i]);
 364     af2Save[i] = pow2(afSave[i]);
 365     efvfSave[i] = efSave[i] * vfSave[i];
 366     vf2af2Save[i] = vf2Save[i] + af2Save[i];
 367   }
 368
 369 }
 370
 371 //**************************************************************************
 372
 373 // The VCKM class.
 374
 375 //*********
 376
 377 // Definitions of static variables. Initialize to all elements zero.
 378
 379 double VCKM::Vsave[5][5]  = { };
 380 double VCKM::V2save[5][5] = { };
 381 double VCKM::V2out[20]    = { };
 382
 383 //*********
 384
 385 // Prepare to handle CKM matrix elements.
 386
 387 void VCKM::initStatic() {
 388
 389   // Read in matrix element values and store them.
 390   Vsave[1][1] = Settings::parm("StandardModel:Vud");
 391   Vsave[1][2] = Settings::parm("StandardModel:Vus");
 392   Vsave[1][3] = Settings::parm("StandardModel:Vub");
 393   Vsave[2][1] = Settings::parm("StandardModel:Vcd");
 394   Vsave[2][2] = Settings::parm("StandardModel:Vcs");
 395   Vsave[2][3] = Settings::parm("StandardModel:Vcb");
 396   Vsave[3][1] = Settings::parm("StandardModel:Vtd");
 397   Vsave[3][2] = Settings::parm("StandardModel:Vts");
 398   Vsave[3][3] = Settings::parm("StandardModel:Vtb");
 399
 400   // Also allow for the potential existence of a fourth generation.
 401   Vsave[1][4] = Settings::parm("FourthGeneration:VubPrime");
 402   Vsave[2][4] = Settings::parm("FourthGeneration:VcbPrime");
 403   Vsave[3][4] = Settings::parm("FourthGeneration:VtbPrime");
 404   Vsave[4][1] = Settings::parm("FourthGeneration:VtPrimed");
 405   Vsave[4][2] = Settings::parm("FourthGeneration:VtPrimes");
 406   Vsave[4][3] = Settings::parm("FourthGeneration:VtPrimeb");
 407   Vsave[4][4] = Settings::parm("FourthGeneration:VtPrimebPrime");
 408
 409   // Calculate squares of matrix elements.
 410   for(int i = 1; i < 5; ++i) for(int j = 1; j < 5; ++j)
 411     V2save[i][j] = pow2(Vsave[i][j]);
 412
 413   // Sum VCKM^2_out sum for given incoming flavour, excluding top as partner.
 414   V2out[1] = V2save[1][1] + V2save[2][1];
 415   V2out[2] = V2save[1][1] + V2save[1][2] + V2save[1][3];
 416   V2out[3] = V2save[1][2] + V2save[2][2];
 417   V2out[4] = V2save[2][1] + V2save[2][2] + V2save[2][3];
 418   V2out[5] = V2save[1][3] + V2save[2][3];
 419   V2out[6] = V2save[3][1] + V2save[3][2] + V2save[3][3];
 420   V2out[7] = V2save[1][4] + V2save[2][4];
 421   V2out[8] = V2save[4][1] + V2save[4][2] + V2save[4][3];
 422   for (int i = 11; i <= 18; ++i) V2out[i] = 1.;
 423
 424 }
 425
 426 //*********
 427
 428 // Return CKM value for incoming flavours (sign irrelevant).
 429
 430 double VCKM::Vid(int id1, int id2) {
 431
 432   // Use absolute sign (want to cover both f -> f' W and f fbar' -> W).
 433   int id1Abs = abs(id1);
 434   int id2Abs = abs(id2);
 435   if (id1Abs == 0 || id2Abs == 0 || (id1Abs + id2Abs)%2 != 1) return 0.;
 436
 437   // Ensure proper order before reading out from Vsave or lepton match.
 438   if (id1Abs%2 == 1) swap(id1Abs, id2Abs);
 439   if (id1Abs <= 8 && id2Abs <= 8) return Vsave[id1Abs/2][(id2Abs + 1)/2];
 440   if ( (id1Abs == 12 || id1Abs == 14 || id1Abs == 16 || id1Abs == 18)
 441     && id2Abs == id1Abs - 1 ) return 1.;
 442
 443   // No more valid cases.
 444   return 0.;
 445
 446 }
 447
 448 //*********
 449
 450 // Return squared CKM value for incoming flavours (sign irrelevant).
 451
 452 double VCKM::V2id(int id1, int id2) {
 453
 454   // Use absolute sign (want to cover both f -> f' W and f fbar' -> W).
 455   int id1Abs = abs(id1);
 456   int id2Abs = abs(id2);
 457   if (id1Abs == 0 || id2Abs == 0 || (id1Abs + id2Abs)%2 != 1) return 0.;
 458
 459   // Ensure proper order before reading out from V2save or lepton match.
 460   if (id1Abs%2 == 1) swap(id1Abs, id2Abs);
 461   if (id1Abs <= 8 && id2Abs <= 8) return V2save[id1Abs/2][(id2Abs + 1)/2];
 462   if ( (id1Abs == 12 || id1Abs == 14 || id1Abs == 16 || id1Abs == 18)
 463     && id2Abs == id1Abs - 1 ) return 1.;
 464
 465   // No more valid cases.
 466   return 0.;
 467
 468 }
 469
 470 //*********
 471
 472 // Pick an outgoing flavour for given incoming one, given CKM mixing.
 473
 474 int VCKM::V2pick(int id) {
 475
 476   // Initial values.
 477   int idIn = abs(id);
 478   int idOut = 0;
 479
 480   // Quarks: need to make random choice.
 481   if (idIn >= 1 && idIn <= 8) {
 482     double V2rndm = Rndm::flat() * V2out[idIn];
 483     if (idIn == 1) idOut = (V2rndm < V2save[1][1]) ? 2 : 4;
 484     else if (idIn == 2) idOut = (V2rndm < V2save[1][1]) ? 1
 485       : ( (V2rndm < V2save[1][1] + V2save[1][2]) ? 3 : 5 );
 486     else if (idIn == 3) idOut = (V2rndm < V2save[1][2]) ? 2 : 4;
 487     else if (idIn == 4) idOut = (V2rndm < V2save[2][1]) ? 1
 488       : ( (V2rndm < V2save[2][1] + V2save[2][2]) ? 3 : 5 );
 489     else if (idIn == 5) idOut = (V2rndm < V2save[1][3]) ? 2 : 4;
 490     else if (idIn == 6) idOut = (V2rndm < V2save[3][1]) ? 1
 491       : ( (V2rndm < V2save[3][1] + V2save[3][2]) ? 3 : 5 );
 492     else if (idIn == 7) idOut = (V2rndm < V2save[1][4]) ? 2 : 4;
 493     else if (idIn == 8) idOut = (V2rndm < V2save[4][1]) ? 1
 494       : ( (V2rndm < V2save[4][1] + V2save[4][2]) ? 3 : 5 );
 495
 496   // Leptons: unambiguous.
 497   } else if (idIn >= 11 && idIn <= 18) {
 498     if (idIn%2 == 1) idOut = idIn + 1;
 499     else idOut             = idIn - 1;
 500   }
 501
 502   // Done. Return with sign.
 503   return ( (id > 0) ? idOut : -idOut );
 504
 505 }
 506
 507 //**************************************************************************
 508
 509 } // end namespace Pythia8