diff --git a/include/Pythia8/MultipartonInteractions.h b/include/Pythia8/MultipartonInteractions.h index 07ea32ce044..f0ccb1067f6 100644 --- a/include/Pythia8/MultipartonInteractions.h +++ b/include/Pythia8/MultipartonInteractions.h @@ -214,7 +214,7 @@ class MultipartonInteractions { SIGMAMBLIMIT; // Initialization data, read from Settings. - bool allowRescatter, allowDoubleRes, canVetoMPI, doPartonVertex; + bool allowRescatter, allowDoubleRes, canVetoMPI, doPartonVertex, alphaSfixRun; int pTmaxMatch, alphaSorder, alphaEMorder, alphaSnfmax, bProfile, processLevel, bSelScale, rescatterMode, nQuarkIn, nSample, enhanceScreening, pT0paramMode; diff --git a/include/Pythia8/SimpleSpaceShower.h b/include/Pythia8/SimpleSpaceShower.h index 2b0e5ed63bc..c65b9cf1cca 100644 --- a/include/Pythia8/SimpleSpaceShower.h +++ b/include/Pythia8/SimpleSpaceShower.h @@ -163,7 +163,7 @@ class SimpleSpaceShower : public SpaceShower { bool doQCDshower, doQEDshowerByQ, doQEDshowerByL, useSamePTasMPI, doWeakShower, doMEcorrections, doMEafterFirst, doPhiPolAsym, doPhiPolAsymHard, doPhiIntAsym, doRapidityOrder, useFixedFacScale, - doSecondHard, canVetoEmission, hasUserHooks, alphaSuseCMW, + doSecondHard, canVetoEmission, hasUserHooks, alphaSuseCMW, alphaSfixRun, singleWeakEmission, vetoWeakJets, weakExternal, doRapidityOrderMPI, doMPI, doDipoleRecoil, doPartonVertex; int pTmaxMatch, pTdampMatch, alphaSorder, alphaSnfmax, alphaEMorder, diff --git a/include/Pythia8/SimpleTimeShower.h b/include/Pythia8/SimpleTimeShower.h index a61d11ee8d4..07ac2980aea 100644 --- a/include/Pythia8/SimpleTimeShower.h +++ b/include/Pythia8/SimpleTimeShower.h @@ -179,7 +179,7 @@ class SimpleTimeShower : public TimeShower { allowBeamRecoil, dampenBeamRecoil, recoilToColoured, useFixedFacScale, allowRescatter, canVetoEmission, doHVshower, brokenHVsym, globalRecoil, useLocalRecoilNow, doSecondHard, hasUserHooks, - singleWeakEmission, alphaSuseCMW, vetoWeakJets, allowMPIdipole, + singleWeakEmission, alphaSuseCMW, alphaSfixRun, vetoWeakJets, allowMPIdipole, weakExternal, recoilDeadCone, doDipoleRecoil, doPartonVertex; int pTmaxMatch, pTdampMatch, alphaSorder, alphaSnfmax, nGluonToQuark, weightGluonToQuark, alphaEMorder, nGammaToQuark, nGammaToLepton, diff --git a/include/Pythia8/StandardModel.h b/include/Pythia8/StandardModel.h index bd5471670af..109ed3fb4e4 100644 --- a/include/Pythia8/StandardModel.h +++ b/include/Pythia8/StandardModel.h @@ -28,7 +28,7 @@ class AlphaStrong { AlphaStrong() : isInit(false), order(0), nfmax(), Lambda3Save(0.), Lambda4Save(0.), Lambda5Save(0.), Lambda6Save(0.), Lambda3Save2(0.), Lambda4Save2(0.), Lambda5Save2(0.), Lambda6Save2(0.), - scale2Min(0.), mc(0.), mb(0.), mt(0.), mc2(0.), mb2(0.), mt2(0.), useCMW(), + scale2Min(0.), mc(0.), mb(0.), mt(0.), mc2(0.), mb2(0.), mt2(0.), useCMW(), fixRunning(), lastCallToFull(false), valueRef(0.), valueNow(0.), scale2Now(0.) {} // Destructor. @@ -36,7 +36,7 @@ class AlphaStrong { // Initialization for given value at M_Z and given order. virtual void init(double valueIn = 0.12, int orderIn = 1, int nfmaxIn = 6, - bool useCMWIn = false); + bool useCMWIn = false, bool fixRunningIn = false); // Set flavour threshold values: m_c, m_b, m_t. virtual void setThresholds(double mcIn, double mbIn, double mtIn) { @@ -81,6 +81,8 @@ class AlphaStrong { // CMW rescaling factors. bool useCMW; + // Fix alphaS running (backported from 8.309) + bool fixRunning; static const double FACCMW3, FACCMW4, FACCMW5, FACCMW6; // Safety margins to avoid getting too close to LambdaQCD. diff --git a/share/Pythia8/xmldoc/CouplingsAndScales.xml b/share/Pythia8/xmldoc/CouplingsAndScales.xml index d264d2cb44b..8d6437da331 100644 --- a/share/Pythia8/xmldoc/CouplingsAndScales.xml +++ b/share/Pythia8/xmldoc/CouplingsAndScales.xml @@ -43,6 +43,12 @@ not go to second order there is no strong reason to use this option, but there is also nothing wrong with it. + + + + +

QED interactions are regulated by the alpha_electromagnetic value at the Q^2 renormalization scale of an interaction. diff --git a/share/Pythia8/xmldoc/MultipartonInteractions.xml b/share/Pythia8/xmldoc/MultipartonInteractions.xml index 77236cad8f6..57aa8cf4077 100644 --- a/share/Pythia8/xmldoc/MultipartonInteractions.xml +++ b/share/Pythia8/xmldoc/MultipartonInteractions.xml @@ -99,6 +99,12 @@ not go to second order there is no strong reason to use this option, but there is also nothing wrong with it. + + + + +

QED interactions are regulated by the alpha_electromagnetic value at the pT^2 scale of an interaction. diff --git a/share/Pythia8/xmldoc/SpacelikeShowers.xml b/share/Pythia8/xmldoc/SpacelikeShowers.xml index 69fc990b880..e11a3df18e5 100644 --- a/share/Pythia8/xmldoc/SpacelikeShowers.xml +++ b/share/Pythia8/xmldoc/SpacelikeShowers.xml @@ -182,6 +182,12 @@ not go to second order there is no strong reason to use this option, but there is also nothing wrong with it. + + + + +

The CMW rescaling of Lambda_QCD (see the section on StandardModelParameters) diff --git a/share/Pythia8/xmldoc/TimelikeShowers.xml b/share/Pythia8/xmldoc/TimelikeShowers.xml index 16c43625646..fce3529270b 100644 --- a/share/Pythia8/xmldoc/TimelikeShowers.xml +++ b/share/Pythia8/xmldoc/TimelikeShowers.xml @@ -188,6 +188,12 @@ not go to second order there is no strong reason to use this option, but there is also nothing wrong with it. + + + + +

The CMW rescaling of Lambda_QCD (see the section on StandardModelParameters) diff --git a/src/MergingHooks.cc b/src/MergingHooks.cc index 13c4d1af72d..b6fbfaa2184 100644 --- a/src/MergingHooks.cc +++ b/src/MergingHooks.cc @@ -2049,15 +2049,17 @@ void MergingHooks::init(){ // Initialise AlphaS objects for reweighting double alphaSvalueFSR = settingsPtr->parm("TimeShower:alphaSvalue"); int alphaSorderFSR = settingsPtr->mode("TimeShower:alphaSorder"); + bool alphaSfixRunFSR= settingsPtr->flag("TimeShower:alphaSfixRun"); int alphaSnfmax = settingsPtr->mode("StandardModel:alphaSnfmax"); int alphaSuseCMWFSR= settingsPtr->flag("TimeShower:alphaSuseCMW"); AlphaS_FSRSave.init(alphaSvalueFSR, alphaSorderFSR, alphaSnfmax, - alphaSuseCMWFSR); + alphaSuseCMWFSR, alphaSfixRunFSR); double alphaSvalueISR = settingsPtr->parm("SpaceShower:alphaSvalue"); int alphaSorderISR = settingsPtr->mode("SpaceShower:alphaSorder"); + bool alphaSfixRunISR= settingsPtr->flag("SpaceShower:alphaSfixRun"); int alphaSuseCMWISR= settingsPtr->flag("SpaceShower:alphaSuseCMW"); AlphaS_ISRSave.init(alphaSvalueISR, alphaSorderISR, alphaSnfmax, - alphaSuseCMWISR); + alphaSuseCMWISR, alphaSfixRunISR); // Initialise AlphaEM objects for reweighting int alphaEMFSRorder = settingsPtr->mode("TimeShower:alphaEMorder"); diff --git a/src/MultipartonInteractions.cc b/src/MultipartonInteractions.cc index afbce90074d..c82117691e0 100644 --- a/src/MultipartonInteractions.cc +++ b/src/MultipartonInteractions.cc @@ -385,6 +385,7 @@ bool MultipartonInteractions::init( bool doMPIinit, int iDiffSysIn, // Parameters of alphaStrong generation. alphaSvalue = settings.parm("MultipartonInteractions:alphaSvalue"); alphaSorder = settings.mode("MultipartonInteractions:alphaSorder"); + alphaSfixRun = settings.flag("MultipartonInteractions:alphaSfixRun"); alphaSnfmax = settings.mode("StandardModel:alphaSnfmax"); // Parameters of alphaEM generation. @@ -502,7 +503,7 @@ bool MultipartonInteractions::init( bool doMPIinit, int iDiffSysIn, enhanceBavg = 1.; // Initialize alpha_strong generation. - alphaS.init( alphaSvalue, alphaSorder, alphaSnfmax, false); + alphaS.init( alphaSvalue, alphaSorder, alphaSnfmax, false, alphaSfixRun); double Lambda3 = alphaS.Lambda3(); // Initialize alphaEM generation. diff --git a/src/Settings.cc b/src/Settings.cc index fc631647fb8..cd72c7a6b4d 100644 --- a/src/Settings.cc +++ b/src/Settings.cc @@ -1839,6 +1839,7 @@ void Settings::resetTuneEE() { // FSR: strong coupling, IR cutoff. resetParm("TimeShower:alphaSvalue"); resetMode("TimeShower:alphaSorder"); + resetFlag("TimeShower:alphaSfixRun"); resetFlag("TimeShower:alphaSuseCMW"); resetParm("TimeShower:pTmin"); resetParm("TimeShower:pTminChgQ"); @@ -1872,6 +1873,7 @@ void Settings::resetTunePP() { // ISR: strong coupling, IR cutoff, coherence and spin correlations. resetParm("SpaceShower:alphaSvalue"); resetMode("SpaceShower:alphaSorder"); + resetParm("SpaceShower:alphaSfixRun"); resetParm("SpaceShower:alphaSuseCMW"); resetFlag("SpaceShower:samePTasMPI"); resetParm("SpaceShower:pT0Ref"); @@ -1886,6 +1888,8 @@ void Settings::resetTunePP() { // MPI: strong coupling, IR regularization, energy scaling. resetParm("MultipartonInteractions:alphaSvalue"); + resetMode("MultipartonInteractions:alphaSorder"); + resetParm("MultipartonInteractions:alphaSfixRun"); resetParm("MultipartonInteractions:pT0Ref"); resetParm("MultipartonInteractions:ecmRef"); resetParm("MultipartonInteractions:ecmPow"); diff --git a/src/SimpleSpaceShower.cc b/src/SimpleSpaceShower.cc index 3711669dcab..40eeb5ceea4 100644 --- a/src/SimpleSpaceShower.cc +++ b/src/SimpleSpaceShower.cc @@ -135,12 +135,13 @@ void SimpleSpaceShower::init( BeamParticle* beamAPtrIn, // Parameters of alphaStrong generation. alphaSvalue = settingsPtr->parm("SpaceShower:alphaSvalue"); alphaSorder = settingsPtr->mode("SpaceShower:alphaSorder"); + alphaSfixRun = settingsPtr->flag("SpaceShower:alphaSfixRun"); alphaSnfmax = settingsPtr->mode("StandardModel:alphaSnfmax"); alphaSuseCMW = settingsPtr->flag("SpaceShower:alphaSuseCMW"); alphaS2pi = 0.5 * alphaSvalue / M_PI; // Initialize alpha_strong generation. - alphaS.init( alphaSvalue, alphaSorder, alphaSnfmax, alphaSuseCMW); + alphaS.init( alphaSvalue, alphaSorder, alphaSnfmax, alphaSuseCMW, alphaSfixRun); // Lambda for 5, 4 and 3 flavours. Lambda5flav = alphaS.Lambda5(); diff --git a/src/SimpleTimeShower.cc b/src/SimpleTimeShower.cc index 55f5a812f70..afc682db631 100644 --- a/src/SimpleTimeShower.cc +++ b/src/SimpleTimeShower.cc @@ -1,5 +1,5 @@ // SimpleTimeShower.cc is a part of the PYTHIA event generator. -// Copyright (C) 2018 Torbjorn Sjostrand. +// Copyright (C) 2019 Torbjorn Sjostrand. // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details. // Please respect the MCnet Guidelines, see GUIDELINES for details. @@ -124,12 +124,13 @@ void SimpleTimeShower::init( BeamParticle* beamAPtrIn, // Parameters of alphaStrong generation. alphaSvalue = settingsPtr->parm("TimeShower:alphaSvalue"); alphaSorder = settingsPtr->mode("TimeShower:alphaSorder"); + alphaSfixRun = settingsPtr->flag("TimeShower:alphaSfixRun"); alphaSnfmax = settingsPtr->mode("StandardModel:alphaSnfmax"); alphaSuseCMW = settingsPtr->flag("TimeShower:alphaSuseCMW"); alphaS2pi = 0.5 * alphaSvalue / M_PI; // Initialize alphaStrong generation. - alphaS.init( alphaSvalue, alphaSorder, alphaSnfmax, alphaSuseCMW); + alphaS.init( alphaSvalue, alphaSorder, alphaSnfmax, alphaSuseCMW, alphaSfixRun); // Lambda for 5, 4 and 3 flavours. Lambda3flav = alphaS.Lambda3(); @@ -3077,12 +3078,16 @@ bool SimpleTimeShower::branch( Event& event, bool isInterleaved) { int iSysSelRec = dipSel->systemRec; // Sometimes need to patch up colType in junction systems. - int idRec = recBef.id(); int colTypeTmp = dipSel->colType; - if (idRec > 0 && idRec < 9 && colTypeTmp > 0) colTypeTmp = -colTypeTmp; - if (idRec > -9 && idRec < 0 && colTypeTmp < 0) colTypeTmp = -colTypeTmp; - if (idRad > 0 && idRad < 9 && colTypeTmp < 0) colTypeTmp = -colTypeTmp; - if (idRad > -9 && idRad < 0 && colTypeTmp > 0) colTypeTmp = -colTypeTmp; + int colTypeRec = particleDataPtr->colType( recBef.id() ); + // Negate colour type if recoiler is initial-state quark. + if (!recBef.isFinal()) colTypeRec = -colTypeRec; + int colTypeRad = particleDataPtr->colType( idRad ); + // Perform junction tests for all colour (anti)triplets. + if (colTypeRec == 1 && colTypeTmp > 0) colTypeTmp = -colTypeTmp; + if (colTypeRec == -1 && colTypeTmp < 0) colTypeTmp = -colTypeTmp; + if (colTypeRad == 1 && colTypeTmp < 0) colTypeTmp = -colTypeTmp; + if (colTypeRad == -1 && colTypeTmp > 0) colTypeTmp = -colTypeTmp; // Default OK for photon, photon_HV or gluon_HV emission. if (dipSel->flavour == 22 || dipSel->flavour == idHV) { diff --git a/src/StandardModel.cc b/src/StandardModel.cc index 068d164b464..4caa236ec86 100644 --- a/src/StandardModel.cc +++ b/src/StandardModel.cc @@ -40,16 +40,17 @@ const double AlphaStrong::FACCMW6 = 1.513; // Initialize alpha_strong calculation by finding Lambda values etc. void AlphaStrong::init( double valueIn, int orderIn, int nfmaxIn, - bool useCMWIn) { + bool useCMWIn, bool fixRunningIn) { // Set default mass thresholds if not already done if (mt <= 1.) setThresholds(1.5, 4.8, 171.0); // Order of alpha_s evaluation.Default values. valueRef = valueIn; - order = max( 0, min( 2, orderIn ) ); - nfmax = max(5,min(6,nfmaxIn)); + order = max( 0, min( 3, orderIn ) ); + nfmax = max( 5, min( 6, nfmaxIn ) ); useCMW = useCMWIn; + fixRunning = fixRunningIn; lastCallToFull = false; Lambda3Save = Lambda4Save = Lambda5Save = Lambda6Save = scale2Min = 0.; @@ -63,27 +64,28 @@ void AlphaStrong::init( double valueIn, int orderIn, int nfmaxIn, Lambda4Save = Lambda5Save * pow(mb/Lambda5Save, 2./25.); Lambda3Save = Lambda4Save * pow(mc/Lambda4Save, 2./27.); - // Second order alpha_s: iterative match at flavour thresholds. + // Second or third order alpha_s: iterative match at flavour thresholds. } else { - // The one-loop coefficients: b1 / b0^2 + // The two-loop coefficients: b1 / b0^2 double b16 = 234. / 441.; double b15 = 348. / 529.; double b14 = 462. / 625.; double b13 = 576. / 729.; - // The two-loop coefficients: b2 * b0 / b1^2 + // The three-loop coefficients: b2 * b0 / b1^2 double b26 = -36855. / 109512.; double b25 = 224687. / 242208.; double b24 = 548575. / 426888.; double b23 = 938709. / 663552.; double logScale, loglogScale, correction, valueIter; - // Find Lambda_5 at m_Z, starting from one-loop value + // Find Lambda_5 at m_Z, starting from one-loop value. Lambda5Save = MZ * exp( -6. * M_PI / (23. * valueRef) ); for (int iter = 0; iter < NITER; ++iter) { logScale = 2. * log(MZ/Lambda5Save); loglogScale = log(logScale); - correction = 1. - b15 * loglogScale / logScale - + pow2(b15 / logScale) * (pow2(loglogScale - 0.5) + b25 - 1.25); + correction = 1. - b15 * loglogScale / logScale; + if (order == 3 || !fixRunning) correction + += pow2(b15 / logScale) * (pow2(loglogScale - 0.5) + b25 - 1.25); valueIter = valueRef / correction; Lambda5Save = MZ * exp( -6. * M_PI / (23. * valueIter) ); } @@ -91,15 +93,17 @@ void AlphaStrong::init( double valueIn, int orderIn, int nfmaxIn, // Find Lambda_6 at m_t, by requiring alphaS(nF=6,m_t) = alphaS(nF=5,m_t) double logScaleT = 2. * log(mt/Lambda5Save); double loglogScaleT = log(logScaleT); - double valueT = 12. * M_PI / (23. * logScaleT) - * (1. - b15 * loglogScaleT / logScaleT - + pow2(b15 / logScaleT) * (pow2(loglogScaleT - 0.5) + b25 - 1.25) ); + correction = 1. - b15 * loglogScaleT / logScaleT; + if (order == 3 || !fixRunning) correction + += pow2(b15 / logScaleT) * (pow2(loglogScaleT - 0.5) + b25 - 1.25); + double valueT = 12. * M_PI / (23. * logScaleT) * correction; Lambda6Save = Lambda5Save; for (int iter = 0; iter < NITER; ++iter) { logScale = 2. * log(mt/Lambda6Save); loglogScale = log(logScale); - correction = 1. - b16 * loglogScale / logScale - + pow2(b16 / logScale) * (pow2(loglogScale - 0.5) + b26 - 1.25); + correction = 1. - b16 * loglogScale / logScale; + if (order == 3 || !fixRunning) correction + += pow2(b16 / logScale) * (pow2(loglogScale - 0.5) + b26 - 1.25); valueIter = valueT / correction; Lambda6Save = mt * exp( -6. * M_PI / (21. * valueIter) ); } @@ -107,30 +111,35 @@ void AlphaStrong::init( double valueIn, int orderIn, int nfmaxIn, // Find Lambda_4 at m_b, by requiring alphaS(nF=4,m_b) = alphaS(nF=5,m_b) double logScaleB = 2. * log(mb/Lambda5Save); double loglogScaleB = log(logScaleB); - double valueB = 12. * M_PI / (23. * logScaleB) - * (1. - b15 * loglogScaleB / logScaleB - + pow2(b15 / logScaleB) * (pow2(loglogScaleB - 0.5) + b25- 1.25) ); + correction = 1. - b15 * loglogScaleB / logScaleB; + if (order == 3 || !fixRunning) correction + += pow2(b15 / logScaleB) * (pow2(loglogScaleB - 0.5) + b25- 1.25); + double valueB = 12. * M_PI / (23. * logScaleB) * correction; Lambda4Save = Lambda5Save; for (int iter = 0; iter < NITER; ++iter) { logScale = 2. * log(mb/Lambda4Save); loglogScale = log(logScale); - correction = 1. - b14 * loglogScale / logScale - + pow2(b14 / logScale) * (pow2(loglogScale - 0.5) + b24 - 1.25); + correction = 1. - b14 * loglogScale / logScale; + if (order == 3 || !fixRunning) correction + += pow2(b14 / logScale) * (pow2(loglogScale - 0.5) + b24 - 1.25); valueIter = valueB / correction; Lambda4Save = mb * exp( -6. * M_PI / (25. * valueIter) ); - } + } + // Find Lambda_3 at m_c, by requiring alphaS(nF=3,m_c) = alphaS(nF=4,m_c) double logScaleC = 2. * log(mc/Lambda4Save); double loglogScaleC = log(logScaleC); - double valueC = 12. * M_PI / (25. * logScaleC) - * (1. - b14 * loglogScaleC / logScaleC - + pow2(b14 / logScaleC) * (pow2(loglogScaleC - 0.5) + b24 - 1.25) ); + correction = 1. - b14 * loglogScaleC / logScaleC; + if (order == 3 || !fixRunning) correction + += pow2(b14 / logScaleC) * (pow2(loglogScaleC - 0.5) + b24 - 1.25); + double valueC = 12. * M_PI / (25. * logScaleC) * correction; Lambda3Save = Lambda4Save; for (int iter = 0; iter < NITER; ++iter) { logScale = 2. * log(mc/Lambda3Save); loglogScale = log(logScale); - correction = 1. - b13 * loglogScale / logScale - + pow2(b13 / logScale) * (pow2(loglogScale - 0.5) + b23 - 1.25); + correction = 1. - b13 * loglogScale / logScale; + if (order == 3 || !fixRunning) correction + += pow2(b13 / logScale) * (pow2(loglogScale - 0.5) + b23 - 1.25); valueIter = valueC / correction; Lambda3Save = mc * exp( -6. * M_PI / (27. * valueIter) ); } @@ -146,7 +155,7 @@ void AlphaStrong::init( double valueIn, int orderIn, int nfmaxIn, // Impose SAFETYMARGINs to prevent getting too close to LambdaQCD. if (order == 1) scale2Min = pow2(SAFETYMARGIN1 * Lambda3Save); - else if (order == 2) scale2Min = pow2(SAFETYMARGIN2 * Lambda3Save); + else if (order > 1) scale2Min = pow2(SAFETYMARGIN2 * Lambda3Save); // Save squares of mass and Lambda values as well. Lambda3Save2 = pow2(Lambda3Save); @@ -217,9 +226,10 @@ double AlphaStrong::alphaS( double scale2) { } double logScale = log(scale2/Lambda2); double loglogScale = log(logScale); - valueNow = 12. * M_PI / (b0 * logScale) - * ( 1. - b1 * loglogScale / logScale - + pow2(b1 / logScale) * (pow2(loglogScale - 0.5) + b2 - 1.25) ); + double correction = 1. - b1 * loglogScale / logScale; + if (order == 3 || !fixRunning) correction + += pow2(b1 / logScale) * (pow2(loglogScale - 0.5) + b2 - 1.25); + valueNow = 12. * M_PI / (b0 * logScale) * correction; } // Done. @@ -247,7 +257,7 @@ double AlphaStrong::alphaS1Ord( double scale2) { if (order == 0) { valueNow = valueRef; - // First/second order alpha_s: differs by mass region. + // First/second/third order alpha_s: differs by mass region. } else { if (scale2 > mt2 && nfmax >= 6) valueNow = 12. * M_PI / (21. * log(scale2/Lambda6Save2)); @@ -264,7 +274,7 @@ double AlphaStrong::alphaS1Ord( double scale2) { //-------------------------------------------------------------------------- -// Calculates the second-order extra factor in alpha_s. +// Calculates the second- or third-order extra factor in alpha_s. // (To be combined with alphaS1Ord.) double AlphaStrong::alphaS2OrdCorr( double scale2) { @@ -297,8 +307,10 @@ double AlphaStrong::alphaS2OrdCorr( double scale2) { } double logScale = log(scale2/Lambda2); double loglogScale = log(logScale); - return ( 1. - b1 * loglogScale / logScale - + pow2(b1 / logScale) * (pow2(loglogScale - 0.5) + b2 - 1.25) ); + double correction = 1. - b1 * loglogScale / logScale; + if (order == 3 || !fixRunning) correction + += pow2(b1 / logScale) * (pow2(loglogScale - 0.5) + b2 - 1.25); + return correction; } @@ -435,8 +447,9 @@ void CoupSM::init(Settings& settings, Rndm* rndmPtrIn) { // Initialize the local AlphaStrong instance. double alphaSvalue = settings.parm("SigmaProcess:alphaSvalue"); int alphaSorder = settings.mode("SigmaProcess:alphaSorder"); + int alphaSfixRun = settings.flag("SigmaProcess:alphaSfixRun"); int alphaSnfmax = settings.mode("StandardModel:alphaSnfmax"); - alphaSlocal.init( alphaSvalue, alphaSorder, alphaSnfmax, false); + alphaSlocal.init( alphaSvalue, alphaSorder, alphaSnfmax, false, alphaSfixRun); // Initialize the local AlphaEM instance. int order = settings.mode("SigmaProcess:alphaEMorder"); diff --git a/src/UserHooks.cc b/src/UserHooks.cc index 0489e1544b4..ff0ba2c1707 100644 --- a/src/UserHooks.cc +++ b/src/UserHooks.cc @@ -246,15 +246,18 @@ double SuppressSmallPT::multiplySigmaBy( const SigmaProcess* sigmaProcessPtr, // alternatively as for hard processes. double alphaSvalue; int alphaSorder; + bool alphaSfixRun; int alphaSnfmax = settingsPtr->mode("StandardModel:alphaSnfmax"); if (useSameAlphaSasMPI) { alphaSvalue = settingsPtr->parm("MultipartonInteractions:alphaSvalue"); alphaSorder = settingsPtr->mode("MultipartonInteractions:alphaSorder"); + alphaSfixRun = settingsPtr->flag("MultipartonInteractions:alphaSfixRun"); } else { alphaSvalue = settingsPtr->parm("SigmaProcess:alphaSvalue"); alphaSorder = settingsPtr->mode("SigmaProcess:alphaSorder"); + alphaSfixRun = settingsPtr->flag("SigmaProcess:alphaSfixRun"); } - alphaS.init( alphaSvalue, alphaSorder, alphaSnfmax, false); + alphaS.init( alphaSvalue, alphaSorder, alphaSnfmax, false, alphaSfixRun); // Initialization finished. isInit = true;