PYTHIA8/pythia8140/src/StandardModel.cxx

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