PYTHIA8/pythia8175/src/StandardModel.cxx

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