diff --git a/Shower/QTilde/Couplings/ShowerAlphaQCD.cc b/Shower/QTilde/Couplings/ShowerAlphaQCD.cc --- a/Shower/QTilde/Couplings/ShowerAlphaQCD.cc +++ b/Shower/QTilde/Couplings/ShowerAlphaQCD.cc @@ -1,453 +1,383 @@ // -*- C++ -*- // // ShowerAlphaQCD.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the ShowerAlphaQCD class. // #include "ShowerAlphaQCD.h" #include "ThePEG/PDT/ParticleData.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/ParVector.h" #include "ThePEG/Interface/Command.h" +#include "ThePEG/Interface/Deleted.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Utilities/Throw.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Config/Constants.h" +#include "Herwig/Utilities/AlphaS.h" using namespace Herwig; +using Herwig::Math::alphaS; +using Herwig::Math::derivativeAlphaS; DescribeClass describeShowerAlphaQCD("Herwig::ShowerAlphaQCD","HwShower.so"); IBPtr ShowerAlphaQCD::clone() const { return new_ptr(*this); } IBPtr ShowerAlphaQCD::fullclone() const { return new_ptr(*this); } void ShowerAlphaQCD::persistentOutput(PersistentOStream & os) const { - os << _asType << _asMaxNP << ounit(_qmin,GeV) << _nloop << _lambdaopt << _thresopt - << ounit(_lambdain,GeV) << _alphain << _inopt - << _tolerance << _maxtry << _alphamin + os << _asType << _asMaxNP << ounit(_qmin,GeV) << _nloop << _thresopt + << _alphain << _tolerance << _maxtry << _alphamin << ounit(_thresholds,GeV) << ounit(_lambda,GeV) << _val0 << ounit(_optInputScale,GeV) << ounit(_quarkMasses,GeV); } void ShowerAlphaQCD::persistentInput(PersistentIStream & is, int) { - is >> _asType >> _asMaxNP >> iunit(_qmin,GeV) >> _nloop >> _lambdaopt >> _thresopt - >> iunit(_lambdain,GeV) >> _alphain >> _inopt - >> _tolerance >> _maxtry >> _alphamin + is >> _asType >> _asMaxNP >> iunit(_qmin,GeV) >> _nloop >> _thresopt + >> _alphain >> _tolerance >> _maxtry >> _alphamin >> iunit(_thresholds,GeV) >> iunit(_lambda,GeV) >> _val0 >> iunit(_optInputScale,GeV) >> iunit(_quarkMasses,GeV); } void ShowerAlphaQCD::Init() { static ClassDocumentation documentation ("This (concrete) class describes the QCD alpha running."); static Switch intAsType ("NPAlphaS", "Behaviour of AlphaS in the NP region", &ShowerAlphaQCD::_asType, 1, false, false); static SwitchOption intAsTypeZero (intAsType, "Zero","zero below Q_min", 1); static SwitchOption intAsTypeConst (intAsType, "Const","const as(qmin) below Q_min", 2); static SwitchOption intAsTypeLin (intAsType, "Linear","growing linearly below Q_min", 3); static SwitchOption intAsTypeQuad (intAsType, "Quadratic","growing quadratically below Q_min", 4); static SwitchOption intAsTypeExx1 (intAsType, "Exx1", "quadratic from AlphaMaxNP down to as(Q_min)", 5); static SwitchOption intAsTypeExx2 (intAsType, "Exx2", "const = AlphaMaxNP below Q_min", 6); // default such that as(qmin) = 1 in the current parametrization. // min = Lambda3 static Parameter intQmin ("Qmin", "Q < Qmin is treated with NP parametrization as of (unit [GeV])," " if negative it is solved for at initialisation such that alpha_S(Qmin)=AlphaMaxNP", &ShowerAlphaQCD::_qmin, GeV, 0.630882*GeV, 0.330445*GeV, 100.0*GeV,false,false,false); static Parameter interfaceAlphaMaxNP ("AlphaMaxNP", "Max value of alpha in NP region, only relevant if NPAlphaS = 5,6", &ShowerAlphaQCD::_asMaxNP, 1.0, 0., 100.0, false, false, Interface::limited); static Parameter interfaceNumberOfLoops ("NumberOfLoops", "The number of loops to use in the alpha_S calculation", &ShowerAlphaQCD::_nloop, 3, 1, 3, false, false, Interface::limited); - static Switch interfaceLambdaOption - ("LambdaOption", - "Option for the calculation of the Lambda used in the simulation from the input" - " Lambda_MSbar", - &ShowerAlphaQCD::_lambdaopt, false, false, false); - static SwitchOption interfaceLambdaOptionfalse - (interfaceLambdaOption, - "Same", - "Use the same value", - false); - static SwitchOption interfaceLambdaOptionConvert - (interfaceLambdaOption, - "Convert", - "Use the conversion to the Herwig scheme from NPB349, 635", - true); + static Deleted delInputOption + ("InputOption", "The old default (1) is now the only choice"); - static Parameter interfaceLambdaQCD - ("LambdaQCD", - "Input value of Lambda_MSBar", - &ShowerAlphaQCD::_lambdain, MeV, 0.208364*GeV, 100.0*MeV, 500.0*MeV, - false, false, Interface::limited); + static Deleted delAlphaMZ + ("AlphaMZ", "Renamed to AlphaIn."); static Parameter interfaceAlphaMZ - ("AlphaMZ", - "The input value of the strong coupling at the Z mass ", + ("AlphaIn", + "The input value of the strong coupling at the chosen InputScale (default: MZ)", &ShowerAlphaQCD::_alphain, 0.118, 0.1, 0.2, false, false, Interface::limited); - static Switch interfaceInputOption - ("InputOption", - "Option for inputing the initial value of the coupling", - &ShowerAlphaQCD::_inopt, true, false, false); - static SwitchOption interfaceInputOptionAlphaMZ - (interfaceInputOption, - "AlphaMZ", - "Use the value of alpha at MZ to calculate the coupling", - true); - static SwitchOption interfaceInputOptionLambdaQCD - (interfaceInputOption, - "LambdaQCD", - "Use the input value of Lambda to calculate the coupling", - false); - static Parameter interfaceTolerance ("Tolerance", "The tolerance for discontinuities in alphaS at thresholds.", &ShowerAlphaQCD::_tolerance, 1e-10, 1e-20, 1e-4, false, false, Interface::limited); static Parameter interfaceMaximumIterations ("MaximumIterations", "The maximum number of iterations for the Newton-Raphson method to converge.", &ShowerAlphaQCD::_maxtry, 100, 10, 1000, false, false, Interface::limited); static Switch interfaceThresholdOption ("ThresholdOption", "Whether to use the consistuent or normal masses for the thresholds", &ShowerAlphaQCD::_thresopt, true, false, false); static SwitchOption interfaceThresholdOptionCurrent (interfaceThresholdOption, "Current", "Use the current masses", true); static SwitchOption interfaceThresholdOptionConstituent (interfaceThresholdOption, "Constituent", "Use the constitent masses.", false); static Command interfaceValue ("Value", "", &ShowerAlphaQCD::value, false); static Command interfacecheck ("check", "check", &ShowerAlphaQCD::check, false); + + + static Parameter interfaceInputScale ("InputScale", "An optional input scale. MZ will be used if not set.", - &ShowerAlphaQCD::_optInputScale, GeV, ZERO, ZERO, 0*GeV, + &ShowerAlphaQCD::_optInputScale, GeV, 91.1876_GeV, ZERO, 0*GeV, false, false, Interface::lowerlim); + + + static ParVector interfaceQuarkMasses ("QuarkMasses", "The quark masses to be used instead of the masses set in the particle data.", &ShowerAlphaQCD::_quarkMasses, GeV, -1, 0.0*GeV, 0.0*GeV, 0*GeV, false, false, Interface::lowerlim); } void ShowerAlphaQCD::doinit() { ShowerAlpha::doinit(); // calculate the value of 5-flavour lambda // evaluate the initial // value of Lambda from alphas if needed using Newton-Raphson - if(_inopt) { - _lambda[2]=computeLambda(_optInputScale == ZERO ? - getParticleData(ParticleID::Z0)->mass() : - _optInputScale, - _alphain,5); - } - // otherwise it was an input parameter - else{_lambda[2]=_lambdain;} - // convert lambda to the Monte Carlo scheme if needed - using Constants::pi; - if(_lambdaopt) _lambda[2] *=exp(0.5*(67.-3.*sqr(pi)-50./3.)/23.)/sqrt(2.); + _lambda[2]=computeLambda(_optInputScale,_alphain,5); + // compute the threshold matching // top threshold for(int ix=1;ix<4;++ix) { if ( _quarkMasses.empty() ) { if(_thresopt) _thresholds[ix]=getParticleData(ix+3)->mass(); else _thresholds[ix]=getParticleData(ix+3)->constituentMass(); } else { // starting at zero rather than one, cf the other alphas's _thresholds[ix] = _quarkMasses[ix+2]; } } // compute 6 flavour lambda by matching at top mass using Newton Raphson - _lambda[3]=computeLambda(_thresholds[3],alphaS(_thresholds[3],_lambda[2],5),6); + _lambda[3]=computeLambda(_thresholds[3],alphaS(_thresholds[3],_lambda[2],5,_nloop),6); // bottom threshold // compute 4 flavour lambda by matching at bottom mass using Newton Raphson - _lambda[1]=computeLambda(_thresholds[2],alphaS(_thresholds[2],_lambda[2],5),4); + _lambda[1]=computeLambda(_thresholds[2],alphaS(_thresholds[2],_lambda[2],5,_nloop),4); // charm threshold // compute 3 flavour lambda by matching at charm mass using Newton Raphson - _lambda[0]=computeLambda(_thresholds[1],alphaS(_thresholds[1],_lambda[1],4),3); + _lambda[0]=computeLambda(_thresholds[1],alphaS(_thresholds[1],_lambda[1],4,_nloop),3); // if qmin less than zero solve for alphaS(_qmin) = 1. if(_qmin() << "The value of Qmin is less than Lambda_3 in" << " ShowerAlphaQCD::doinit " << Exception::abortnow; Energy high = 10*GeV; while (value(sqr(high))>_asMaxNP) high *=2.; double as; do { _qmin = 0.5*(low+high); as = value(sqr(_qmin)); if(_asMaxNP>as) high = _qmin; else if(_asMaxNP _tolerance ); } // final threshold is qmin _thresholds[0]=_qmin; // value of alphaS at threshold pair nflam = getLamNfTwoLoop(_qmin); - _val0 = alphaS(_qmin, nflam.second, nflam.first); + _val0 = alphaS(_qmin, nflam.second, nflam.first, _nloop); // compute the maximum value of as if ( _asType < 5 ) _alphamin = _val0; else _alphamin = max(_asMaxNP, _val0); // check consistency lambda_3 < qmin if(_lambda[0]>_qmin) Throw() << "The value of Qmin is less than Lambda_3 in" << " ShowerAlphaQCD::doinit " << Exception::abortnow; } string ShowerAlphaQCD::check(string args) { doinit(); istringstream argin(args); double Q_low, Q_high; long n_steps; argin >> Q_low >> Q_high >> n_steps; string fname; argin >> fname; generator()->log() << "checking alpha_s in range [" << Q_low << "," << Q_high << "] GeV in " << n_steps << " steps.\nResults are written to " << fname << "\n"; double step_width = (Q_high-Q_low)/n_steps; ofstream out (fname.c_str()); for (long k = 0; k <= n_steps; ++k) { Energy Q = Q_low*GeV + k*step_width*GeV; out << (Q/GeV) << " " << value(Q*Q) << "\n"; } return "alpha_s check finished"; } double ShowerAlphaQCD::value(const Energy2 scale) const { Energy q = scaleFactor()*sqrt(scale); double val(0.); // normal case if (q >= _qmin) { pair nflam = getLamNfTwoLoop(q); - val = alphaS(q, nflam.second, nflam.first); + val = alphaS(q, nflam.second, nflam.first, _nloop); } // special handling if the scale is less than Qmin else { switch (_asType) { case 1: // flat, zero; the default type with no NP effects. val = 0.; break; case 2: // flat, non-zero alpha_s = alpha_s(q2min). val = _val0; break; case 3: // linear in q val = _val0*q/_qmin; break; case 4: // quadratic in q val = _val0*sqr(q/_qmin); break; case 5: // quadratic in q, starting off at asMaxNP, ending on as(qmin) val = (_val0 - _asMaxNP)*sqr(q/_qmin) + _asMaxNP; break; case 6: // just asMaxNP and constant val = _asMaxNP; break; } } return val; } double ShowerAlphaQCD::overestimateValue() const { return _alphamin; } double ShowerAlphaQCD::ratio(const Energy2 scale, double factor) const { Energy q = scaleFactor()*factor*sqrt(scale); double val(0.); // normal case if (q >= _qmin) { pair nflam = getLamNfTwoLoop(q); - val = alphaS(q, nflam.second, nflam.first); + val = alphaS(q, nflam.second, nflam.first, _nloop); } // special handling if the scale is less than Qmin else { switch (_asType) { case 1: // flat, zero; the default type with no NP effects. val = 0.; break; case 2: // flat, non-zero alpha_s = alpha_s(q2min). val = _val0; break; case 3: // linear in q val = _val0*q/_qmin; break; case 4: // quadratic in q val = _val0*sqr(q/_qmin); break; case 5: // quadratic in q, starting off at asMaxNP, ending on as(qmin) val = (_val0 - _asMaxNP)*sqr(q/_qmin) + _asMaxNP; break; case 6: // just asMaxNP and constant val = _asMaxNP; break; } } // denominator return val/_alphamin; } string ShowerAlphaQCD::value (string scale) { istringstream readscale(scale); double inScale; readscale >> inScale; Energy theScale = inScale * GeV; initialize(); ostringstream showvalue (""); showvalue << "alpha_s (" << theScale/GeV << " GeV) = " << value (sqr(theScale)); return showvalue.str(); } pair ShowerAlphaQCD::getLamNfTwoLoop(Energy q) const { short nf = 6; // get lambda and nf according to the thresholds if (q < _thresholds[1]) nf = 3; else if (q < _thresholds[2]) nf = 4; else if (q < _thresholds[3]) nf = 5; return pair(nf, _lambda[nf-3]); } Energy ShowerAlphaQCD::computeLambda(Energy match, double alpha, unsigned int nflav) const { Energy lamtest=200.0*MeV; double xtest; unsigned int ntry=0; do { ++ntry; xtest = log(sqr(match/lamtest)); - xtest += (alpha-alphaS(match,lamtest,nflav))/derivativealphaS(match,lamtest,nflav); + xtest += (alpha-alphaS(match,lamtest,nflav,_nloop))/derivativeAlphaS(match,lamtest,nflav,_nloop); Energy newLambda = match/exp(0.5*xtest); lamtest = newLambda _tolerance && ntry < _maxtry); + while(abs(alpha-alphaS(match,lamtest,nflav,_nloop)) > _tolerance && ntry < _maxtry); return lamtest; } - -double ShowerAlphaQCD::derivativealphaS(Energy q, Energy lam, int nf) const { - using Constants::pi; - double lx = log(sqr(q/lam)); - double b0 = 11. - 2./3.*nf; - double b1 = 51. - 19./3.*nf; - double b2 = 2857. - 5033./9.*nf + 325./27.*sqr(nf); - if(_nloop==1) - return -4.*pi/(b0*sqr(lx)); - else if(_nloop==2) - return -4.*pi/(b0*sqr(lx))*(1.+2.*b1/sqr(b0)/lx*(1.-2.*log(lx))); - else - return -4.*pi/(b0*sqr(lx))* - (1. + 2.*b1/sqr(b0)/lx*(1.-2.*log(lx)) - + 4.*sqr(b1)/(sqr(sqr(b0))*sqr(lx))*(1. - 2.*log(lx) - + 3.*(sqr(log(lx) - 0.5)+b2*b0/(8.*sqr(b1))-1.25))); -} - -double ShowerAlphaQCD::alphaS(Energy q, Energy lam, int nf) const { - using Constants::pi; - double lx(log(sqr(q/lam))); - double b0 = 11. - 2./3.*nf; - double b1 = 51. - 19./3.*nf; - double b2 = 2857. - 5033./9.*nf + 325./27.*sqr(nf); - // one loop - if(_nloop==1) - {return 4.*pi/(b0*lx);} - // two loop - else if(_nloop==2) { - return 4.*pi/(b0*lx)*(1.-2.*b1/sqr(b0)*log(lx)/lx); - } - // three loop - else - {return 4.*pi/(b0*lx)*(1.-2.*b1/sqr(b0)*log(lx)/lx + - 4.*sqr(b1)/(sqr(sqr(b0))*sqr(lx))* - (sqr(log(lx) - 0.5) + b2*b0/(8.*sqr(b1)) - 5./4.));} -} - diff --git a/Shower/QTilde/Couplings/ShowerAlphaQCD.h b/Shower/QTilde/Couplings/ShowerAlphaQCD.h --- a/Shower/QTilde/Couplings/ShowerAlphaQCD.h +++ b/Shower/QTilde/Couplings/ShowerAlphaQCD.h @@ -1,308 +1,276 @@ // -*- C++ -*- // // ShowerAlphaQCD.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_ShowerAlphaQCD_H #define HERWIG_ShowerAlphaQCD_H // // This is the declaration of the ShowerAlphaQCD class. // #include "Herwig/Shower/ShowerAlpha.h" namespace Herwig { using namespace ThePEG; /** \ingroup Shower * * This concrete class provides the definition of the * pure virtual function value() and overestimateValue() for the * strong coupling. * * A number of different options for the running of the coupling * and its initial definition are supported. * * @see \ref ShowerAlphaQCDInterfaces "The interfaces" * defined for ShowerAlphaQCD. */ class ShowerAlphaQCD: public ShowerAlpha { public: /** * The default constructor. */ ShowerAlphaQCD() : ShowerAlpha(), _qmin(0.630882*GeV), _asType(1), _asMaxNP(1.0), _thresholds(4), _lambda(4), - _nloop(3),_lambdaopt(false),_thresopt(false), - _lambdain(0.208364*GeV),_alphain(0.118),_inopt(true),_tolerance(1e-10), - _maxtry(100),_alphamin(0.),_val0(1.), _optInputScale(ZERO) {} + _nloop(3),_thresopt(false), + _alphain(0.118),_tolerance(1e-10), + _maxtry(100),_alphamin(0.),_val0(1.), _optInputScale(91.1876_GeV) {} public: /** * Methods to return the coupling */ //@{ /** * It returns the running coupling value evaluated at the input scale * multiplied by the scale factor scaleFactor(). * @param scale The scale * @return The coupling */ virtual double value(const Energy2 scale) const; /** * It returns the running coupling value evaluated at the input scale * multiplied by the scale factor scaleFactor(). */ virtual double overestimateValue() const; /** * Return the ratio of the coupling at the scale to the overestimated value */ virtual double ratio(const Energy2 scale,double factor =1.) const; /** * Initialize this coupling. */ virtual void initialize() { doinit(); } /** * A command to initialize the coupling and write * its value at the scale given by the argument (in GeV) */ string value(string); /** * Match thresholds and write alpha_s * specified file; arguments are * Q_low/GeV Q_high/GeV n_steps filename */ string check(string args); //@} /** * Get the value of \f$\Lambda_{\rm QCd}\f$ * @param nf number of flavours */ Energy lambdaQCD(unsigned int nf) { if (nf <= 3) return _lambda[0]; else if (nf==4 || nf==5) return _lambda[nf-3]; else return _lambda[3]; } /** * Return the quark masses to be used; if not empty these masses * should be considered instead of the ones set in the particle data * objects. */ const vector& quarkMasses() const { return _quarkMasses; } public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); //@} private: /** * Member functions which calculate the coupling */ //@{ /** - * The 1,2,3-loop parametrization of \f$\alpha_S\f$. - * @param q The scale - * @param lam \f$\Lambda_{\rm QCD}\f$ - * @param nf The number of flavours - */ - double alphaS(Energy q, Energy lam, int nf) const; - - /** - * The derivative of \f$\alpha_S\f$ with respect to \f$\ln(Q^2/\Lambda^2)\f$ - * @param q The scale - * @param lam \f$\Lambda_{\rm QCD}\f$ - * @param nf The number of flavours - */ - double derivativealphaS(Energy q, Energy lam, int nf) const; - - /** * Compute the value of \f$Lambda\f$ needed to get the input value of * the strong coupling at the scale given for the given number of flavours * using the Newton-Raphson method * @param match The scale for the coupling * @param alpha The input coupling * @param nflav The number of flavours */ Energy computeLambda(Energy match, double alpha, unsigned int nflav) const; /** * Return the value of \f$\Lambda\f$ and the number of flavours at the scale. * @param q The scale * @return The number of flavours at the scale and \f$\Lambda\f$. */ pair getLamNfTwoLoop(Energy q) const; //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ - ShowerAlphaQCD & operator=(const ShowerAlphaQCD &); + ShowerAlphaQCD & operator=(const ShowerAlphaQCD &) = delete; private: /** * Minimum value of the scale */ Energy _qmin; /** * Parameter controlling the behaviour of \f$\alpha_S\f$ in the * non-perturbative region. */ int _asType; /** * Another parameter, a possible (maximum) value of alpha in the * non-perturbative region. */ double _asMaxNP; /** * Thresholds for the different number of flavours */ vector _thresholds; /** * \f$\Lambda\f$ for the different number of flavours */ vector _lambda; /** * Option for the number of loops */ unsigned int _nloop; /** - * Option for the translation between \f$\Lambda_{\bar{MS}}\f$ and - * \f$\Lambda_{\rm Herwig}\f$ - */ - bool _lambdaopt; - - /** * Option for the threshold masses */ bool _thresopt; /** - * Input value of Lambda - */ - Energy _lambdain; - - /** * Input value of \f$alpha_S(M_Z)\f$ */ double _alphain; /** - * Option for the calculation of Lambda from input parameters - */ - bool _inopt; - - /** * Tolerance for discontinuities at the thresholds */ double _tolerance; /** * Maximum number of iterations for the Newton-Raphson method to converge */ unsigned int _maxtry; /** * The minimum value of the coupling */ double _alphamin; /** * Value of \f$\alpha_S\f$ at the minimum scale */ double _val0; /** * An optional input scale to be used for the input alphas; if zero MZ will * be used out of the particle data object. */ Energy _optInputScale; /** * The quark masses to be used; if not empty these masses should be * considered instead of the ones set in the particle data objects. */ vector _quarkMasses; }; } #endif /* HERWIG_ShowerAlphaQCD_H */ diff --git a/Utilities/AlphaS.h b/Utilities/AlphaS.h new file mode 100644 --- /dev/null +++ b/Utilities/AlphaS.h @@ -0,0 +1,72 @@ +// -*- C++ -*- +// +// AlphaS.h is a part of Herwig - A multi-purpose Monte Carlo event generator +// Copyright (C) 2018 The Herwig Collaboration +// +// Herwig is licenced under version 3 of the GPL, see COPYING for details. +// Please respect the MCnet academic guidelines, see GUIDELINES for details. +// +#ifndef HERWIG_UtilAlphaS_H +#define HERWIG_UtilAlphaS_H + +namespace Herwig { + +namespace Math { + + /** + * The derivative of \f$\alpha_S\f$ with respect to \f$\ln(Q^2/\Lambda^2)\f$ + * @param q The scale + * @param lam \f$\Lambda_{\rm QCD}\f$ + * @param nf The number of flavours + */ +inline double derivativeAlphaS(Energy q, Energy lam, + unsigned int nf, unsigned int nloop) { + using Constants::pi; + double lx = log(sqr(q/lam)); + double b0 = 11. - 2./3.*nf; + double b1 = 51. - 19./3.*nf; + double b2 = 2857. - 5033./9.*nf + 325./27.*sqr(nf); + if(nloop==1) + return -4.*pi/(b0*sqr(lx)); + else if(nloop==2) + return -4.*pi/(b0*sqr(lx))*(1.+2.*b1/sqr(b0)/lx*(1.-2.*log(lx))); + else + return -4.*pi/(b0*sqr(lx))* + (1. + 2.*b1/sqr(b0)/lx*(1.-2.*log(lx)) + + 4.*sqr(b1)/(sqr(sqr(b0))*sqr(lx))*(1. - 2.*log(lx) + + 3.*(sqr(log(lx) - 0.5)+b2*b0/(8.*sqr(b1))-1.25))); +} + + /** + * The 1,2,3-loop parametrization of \f$\alpha_S\f$. + * @param q The scale + * @param lam \f$\Lambda_{\rm QCD}\f$ + * @param nf The number of flavours + */ +inline double alphaS(Energy q, Energy lam, + unsigned int nf, unsigned int nloop) { + using Constants::pi; + double lx(log(sqr(q/lam))); + double b0 = 11. - 2./3.*nf; + double b1 = 51. - 19./3.*nf; + double b2 = 2857. - 5033./9.*nf + 325./27.*sqr(nf); + // one loop + if(nloop==1) + {return 4.*pi/(b0*lx);} + // two loop + else if(nloop==2) { + return 4.*pi/(b0*lx)*(1.-2.*b1/sqr(b0)*log(lx)/lx); + } + // three loop + else + {return 4.*pi/(b0*lx)*(1.-2.*b1/sqr(b0)*log(lx)/lx + + 4.*sqr(b1)/(sqr(sqr(b0))*sqr(lx))* + (sqr(log(lx) - 0.5) + b2*b0/(8.*sqr(b1)) - 5./4.));} +} + + +} + +} + +#endif \ No newline at end of file diff --git a/Utilities/Makefile.am b/Utilities/Makefile.am --- a/Utilities/Makefile.am +++ b/Utilities/Makefile.am @@ -1,46 +1,46 @@ SUBDIRS = XML Statistics noinst_LTLIBRARIES = libHwUtils.la libHwUtils_la_SOURCES = \ EnumParticles.h \ Interpolator.tcc Interpolator.h \ Kinematics.cc Kinematics.h \ Progress.h Progress.cc \ Maths.h Maths.cc \ StandardSelectors.cc StandardSelectors.h\ Histogram.cc Histogram.fh Histogram.h \ GaussianIntegrator.cc GaussianIntegrator.h \ GaussianIntegrator.tcc \ Statistic.h HerwigStrategy.cc HerwigStrategy.h \ GSLIntegrator.h GSLIntegrator.tcc \ GSLBisection.h GSLBisection.tcc GSLHelper.h \ expm-1.h \ -HiggsLoopFunctions.h +HiggsLoopFunctions.h AlphaS.h nodist_libHwUtils_la_SOURCES = hgstamp.inc BUILT_SOURCES = hgstamp.inc CLEANFILES = hgstamp.inc HGVERSION := $(shell hg -R $(top_srcdir) parents --template '"Herwig {node|short} ({branch})"' 2> /dev/null || echo \"$(PACKAGE_STRING)\" || true ) .PHONY: update_hgstamp hgstamp.inc: update_hgstamp @[ -f $@ ] || touch $@ @echo '$(HGVERSION)' | cmp -s $@ - || echo '$(HGVERSION)' > $@ libHwUtils_la_LIBADD = \ XML/libHwXML.la \ Statistics/libHwStatistics.la check_PROGRAMS = utilities_test utilities_test_SOURCES = \ tests/utilitiesTestsMain.cc \ tests/utilitiesTestsGlobalFixture.h \ tests/utilitiesTestsKinematics.h \ tests/utilitiesTestMaths.h \ tests/utilitiesTestsStatistic.h utilities_test_LDADD = $(BOOST_UNIT_TEST_FRAMEWORK_LIBS) $(THEPEGLIB) -ldl libHwUtils.la utilities_test_LDFLAGS = $(AM_LDFLAGS) -export-dynamic $(BOOST_UNIT_TEST_FRAMEWORK_LDFLAGS) $(THEPEGLDFLAGS) utilities_test_CPPFLAGS = $(AM_CPPFLAGS) $(BOOST_CPPFLAGS) TESTS = utilities_test diff --git a/src/defaults/Shower.in b/src/defaults/Shower.in --- a/src/defaults/Shower.in +++ b/src/defaults/Shower.in @@ -1,313 +1,312 @@ # -*- ThePEG-repository -*- ############################################################ # Setup of default parton shower # # Useful switches for users are marked near the top of # this file. # # Don't edit this file directly, but reset the switches # in your own input files! ############################################################ library HwMPI.so library HwShower.so library HwMatching.so mkdir /Herwig/Shower cd /Herwig/Shower create Herwig::QTildeShowerHandler ShowerHandler newdef ShowerHandler:MPIHandler /Herwig/UnderlyingEvent/MPIHandler newdef ShowerHandler:RemDecayer /Herwig/Partons/RemnantDecayer # use LO PDFs for Shower, can be changed later newdef ShowerHandler:PDFA /Herwig/Partons/ShowerLOPDF newdef ShowerHandler:PDFB /Herwig/Partons/ShowerLOPDF newdef ShowerHandler:PDFARemnant /Herwig/Partons/RemnantPDF newdef ShowerHandler:PDFBRemnant /Herwig/Partons/RemnantPDF ##################################### # initial setup, don't change these! ##################################### create Herwig::SplittingGenerator SplittingGenerator create Herwig::ShowerAlphaQCD AlphaQCD create Herwig::ShowerAlphaQED AlphaQED create Herwig::PartnerFinder PartnerFinder newdef PartnerFinder:PartnerMethod 1 newdef PartnerFinder:ScaleChoice 1 create Herwig::KinematicsReconstructor KinematicsReconstructor newdef KinematicsReconstructor:ReconstructionOption Colour3 newdef KinematicsReconstructor:InitialStateReconOption SofterFraction newdef KinematicsReconstructor:InitialInitialBoostOption LongTransBoost newdef /Herwig/Partons/RemnantDecayer:AlphaS AlphaQCD newdef /Herwig/Partons/RemnantDecayer:AlphaEM AlphaQED newdef ShowerHandler:PartnerFinder PartnerFinder newdef ShowerHandler:KinematicsReconstructor KinematicsReconstructor newdef ShowerHandler:SplittingGenerator SplittingGenerator newdef ShowerHandler:SpinCorrelations Yes newdef ShowerHandler:SoftCorrelations Singular ################################################################## # Intrinsic pT # # Recommended: # 1.9 GeV for Tevatron W/Z production. # 2.1 GeV for LHC W/Z production at 10 TeV # 2.2 GeV for LHC W/Z production at 14 TeV # # Set all parameters to 0 to disable ################################################################## newdef ShowerHandler:IntrinsicPtGaussian 1.3*GeV newdef ShowerHandler:IntrinsicPtBeta 0 newdef ShowerHandler:IntrinsicPtGamma 0*GeV newdef ShowerHandler:IntrinsicPtIptmax 0*GeV ############################################################# # Set up truncated shower handler. ############################################################# create Herwig::PowhegShowerHandler PowhegShowerHandler set PowhegShowerHandler:MPIHandler /Herwig/UnderlyingEvent/MPIHandler set PowhegShowerHandler:RemDecayer /Herwig/Partons/RemnantDecayer newdef PowhegShowerHandler:PDFA /Herwig/Partons/ShowerLOPDF newdef PowhegShowerHandler:PDFB /Herwig/Partons/ShowerLOPDF newdef PowhegShowerHandler:PDFARemnant /Herwig/Partons/RemnantPDF newdef PowhegShowerHandler:PDFBRemnant /Herwig/Partons/RemnantPDF newdef PowhegShowerHandler:MPIHandler /Herwig/UnderlyingEvent/MPIHandler newdef PowhegShowerHandler:RemDecayer /Herwig/Partons/RemnantDecayer newdef PowhegShowerHandler:PDFA /Herwig/Partons/ShowerLOPDF newdef PowhegShowerHandler:PDFB /Herwig/Partons/ShowerLOPDF newdef PowhegShowerHandler:PDFARemnant /Herwig/Partons/RemnantPDF newdef PowhegShowerHandler:PDFBRemnant /Herwig/Partons/RemnantPDF newdef PowhegShowerHandler:PartnerFinder PartnerFinder newdef PowhegShowerHandler:KinematicsReconstructor KinematicsReconstructor newdef PowhegShowerHandler:SplittingGenerator SplittingGenerator newdef PowhegShowerHandler:Interactions QCDandQED newdef PowhegShowerHandler:SpinCorrelations Yes newdef PowhegShowerHandler:SoftCorrelations Singular newdef PowhegShowerHandler:IntrinsicPtGaussian 1.3*GeV newdef PowhegShowerHandler:IntrinsicPtBeta 0 newdef PowhegShowerHandler:IntrinsicPtGamma 0*GeV newdef PowhegShowerHandler:IntrinsicPtIptmax 0*GeV newdef PowhegShowerHandler:ReconstructionOption OffShell5 ############################################################# # End of interesting user servicable section. # # Anything that follows below should only be touched if you # know what you're doing. # # Really. ############################################################# # # a few default values newdef ShowerHandler:MECorrMode 1 newdef ShowerHandler:ReconstructionOption OffShell5 newdef AlphaQCD:ScaleFactor 1.0 newdef AlphaQCD:NPAlphaS 2 newdef AlphaQCD:Qmin 0.935 newdef AlphaQCD:NumberOfLoops 2 -newdef AlphaQCD:InputOption 1 -newdef AlphaQCD:AlphaMZ 0.126234 +newdef AlphaQCD:AlphaIn 0.126234 # # # Lets set up all the splittings create Herwig::HalfHalfOneSplitFn QtoQGammaSplitFn set QtoQGammaSplitFn:InteractionType QED set QtoQGammaSplitFn:ColourStructure ChargedChargedNeutral set QtoQGammaSplitFn:AngularOrdered Yes create Herwig::HalfHalfOneSplitFn QtoQGSplitFn newdef QtoQGSplitFn:InteractionType QCD newdef QtoQGSplitFn:ColourStructure TripletTripletOctet set QtoQGSplitFn:AngularOrdered Yes create Herwig::OneOneOneSplitFn GtoGGSplitFn newdef GtoGGSplitFn:InteractionType QCD newdef GtoGGSplitFn:ColourStructure OctetOctetOctet set GtoGGSplitFn:AngularOrdered Yes create Herwig::OneOneOneMassiveSplitFn WtoWGammaSplitFn newdef WtoWGammaSplitFn:InteractionType QED newdef WtoWGammaSplitFn:ColourStructure ChargedChargedNeutral set WtoWGammaSplitFn:AngularOrdered Yes create Herwig::OneHalfHalfSplitFn GtoQQbarSplitFn newdef GtoQQbarSplitFn:InteractionType QCD newdef GtoQQbarSplitFn:ColourStructure OctetTripletTriplet set GtoQQbarSplitFn:AngularOrdered Yes create Herwig::OneHalfHalfSplitFn GammatoQQbarSplitFn newdef GammatoQQbarSplitFn:InteractionType QED newdef GammatoQQbarSplitFn:ColourStructure NeutralChargedCharged set GammatoQQbarSplitFn:AngularOrdered Yes create Herwig::HalfOneHalfSplitFn QtoGQSplitFn newdef QtoGQSplitFn:InteractionType QCD newdef QtoGQSplitFn:ColourStructure TripletOctetTriplet set QtoGQSplitFn:AngularOrdered Yes create Herwig::HalfOneHalfSplitFn QtoGammaQSplitFn newdef QtoGammaQSplitFn:InteractionType QED newdef QtoGammaQSplitFn:ColourStructure ChargedNeutralCharged set QtoGammaQSplitFn:AngularOrdered Yes # # Now the Sudakovs create Herwig::PTCutOff PTCutOff newdef PTCutOff:pTmin 1.222798*GeV create Herwig::SudakovFormFactor SudakovCommon newdef SudakovCommon:Alpha AlphaQCD newdef SudakovCommon:Cutoff PTCutOff newdef SudakovCommon:PDFmax 1.0 cp SudakovCommon QtoQGSudakov newdef QtoQGSudakov:SplittingFunction QtoQGSplitFn newdef QtoQGSudakov:PDFmax 1.9 cp SudakovCommon QtoQGammaSudakov set QtoQGammaSudakov:SplittingFunction QtoQGammaSplitFn set QtoQGammaSudakov:Alpha AlphaQED set QtoQGammaSudakov:PDFmax 1.9 cp QtoQGammaSudakov LtoLGammaSudakov cp PTCutOff LtoLGammaPTCutOff # Technical parameter to stop evolution. set LtoLGammaPTCutOff:pTmin 0.000001 set LtoLGammaSudakov:Cutoff LtoLGammaPTCutOff cp SudakovCommon GtoGGSudakov newdef GtoGGSudakov:SplittingFunction GtoGGSplitFn newdef GtoGGSudakov:PDFmax 2.0 cp SudakovCommon WtoWGammaSudakov newdef WtoWGammaSudakov:SplittingFunction WtoWGammaSplitFn set WtoWGammaSudakov:Alpha AlphaQED cp SudakovCommon GtoQQbarSudakov newdef GtoQQbarSudakov:SplittingFunction GtoQQbarSplitFn newdef GtoQQbarSudakov:PDFmax 120.0 cp SudakovCommon GammatoQQbarSudakov newdef GammatoQQbarSudakov:SplittingFunction GammatoQQbarSplitFn set GammatoQQbarSudakov:Alpha AlphaQED newdef GammatoQQbarSudakov:PDFmax 120.0 cp SudakovCommon GtobbbarSudakov newdef GtobbbarSudakov:SplittingFunction GtoQQbarSplitFn newdef GtobbbarSudakov:PDFmax 40000.0 cp SudakovCommon GtoccbarSudakov newdef GtoccbarSudakov:SplittingFunction GtoQQbarSplitFn newdef GtoccbarSudakov:PDFmax 2000.0 cp SudakovCommon QtoGQSudakov newdef QtoGQSudakov:SplittingFunction QtoGQSplitFn cp SudakovCommon QtoGammaQSudakov newdef QtoGammaQSudakov:SplittingFunction QtoGammaQSplitFn set QtoGammaQSudakov:Alpha AlphaQED cp SudakovCommon utoGuSudakov newdef utoGuSudakov:SplittingFunction QtoGQSplitFn newdef utoGuSudakov:PDFFactor OverOneMinusZ newdef utoGuSudakov:PDFmax 5.0 cp SudakovCommon dtoGdSudakov newdef dtoGdSudakov:SplittingFunction QtoGQSplitFn newdef dtoGdSudakov:PDFFactor OverOneMinusZ # # Now add the final splittings # do SplittingGenerator:AddFinalSplitting u->u,g; QtoQGSudakov do SplittingGenerator:AddFinalSplitting d->d,g; QtoQGSudakov do SplittingGenerator:AddFinalSplitting s->s,g; QtoQGSudakov do SplittingGenerator:AddFinalSplitting c->c,g; QtoQGSudakov do SplittingGenerator:AddFinalSplitting b->b,g; QtoQGSudakov do SplittingGenerator:AddFinalSplitting t->t,g; QtoQGSudakov # do SplittingGenerator:AddFinalSplitting g->g,g; GtoGGSudakov # do SplittingGenerator:AddFinalSplitting g->u,ubar; GtoQQbarSudakov do SplittingGenerator:AddFinalSplitting g->d,dbar; GtoQQbarSudakov do SplittingGenerator:AddFinalSplitting g->s,sbar; GtoQQbarSudakov do SplittingGenerator:AddFinalSplitting g->c,cbar; GtoccbarSudakov do SplittingGenerator:AddFinalSplitting g->b,bbar; GtobbbarSudakov do SplittingGenerator:AddFinalSplitting g->t,tbar; GtoQQbarSudakov # do SplittingGenerator:AddFinalSplitting gamma->u,ubar; GammatoQQbarSudakov do SplittingGenerator:AddFinalSplitting gamma->d,dbar; GammatoQQbarSudakov do SplittingGenerator:AddFinalSplitting gamma->s,sbar; GammatoQQbarSudakov do SplittingGenerator:AddFinalSplitting gamma->c,cbar; GammatoQQbarSudakov do SplittingGenerator:AddFinalSplitting gamma->b,bbar; GammatoQQbarSudakov do SplittingGenerator:AddFinalSplitting gamma->t,tbar; GammatoQQbarSudakov do SplittingGenerator:AddFinalSplitting gamma->e-,e+; GammatoQQbarSudakov do SplittingGenerator:AddFinalSplitting gamma->mu-,mu+; GammatoQQbarSudakov do SplittingGenerator:AddFinalSplitting gamma->tau-,tau+; GammatoQQbarSudakov # do SplittingGenerator:AddFinalSplitting u->u,gamma; QtoQGammaSudakov do SplittingGenerator:AddFinalSplitting d->d,gamma; QtoQGammaSudakov do SplittingGenerator:AddFinalSplitting s->s,gamma; QtoQGammaSudakov do SplittingGenerator:AddFinalSplitting c->c,gamma; QtoQGammaSudakov do SplittingGenerator:AddFinalSplitting b->b,gamma; QtoQGammaSudakov do SplittingGenerator:AddFinalSplitting t->t,gamma; QtoQGammaSudakov do SplittingGenerator:AddFinalSplitting e-->e-,gamma; LtoLGammaSudakov do SplittingGenerator:AddFinalSplitting mu-->mu-,gamma; LtoLGammaSudakov do SplittingGenerator:AddFinalSplitting tau-->tau-,gamma; LtoLGammaSudakov do SplittingGenerator:AddFinalSplitting W+->W+,gamma; WtoWGammaSudakov # # Now lets add the initial splittings. Remember the form a->b,c; means # that the current particle b is given and we backward branch to new # particle a which is initial state and new particle c which is final state # do SplittingGenerator:AddInitialSplitting u->u,g; QtoQGSudakov do SplittingGenerator:AddInitialSplitting d->d,g; QtoQGSudakov do SplittingGenerator:AddInitialSplitting s->s,g; QtoQGSudakov do SplittingGenerator:AddInitialSplitting c->c,g; QtoQGSudakov do SplittingGenerator:AddInitialSplitting b->b,g; QtoQGSudakov do SplittingGenerator:AddInitialSplitting u->u,gamma; QtoQGammaSudakov do SplittingGenerator:AddInitialSplitting d->d,gamma; QtoQGammaSudakov do SplittingGenerator:AddInitialSplitting s->s,gamma; QtoQGammaSudakov do SplittingGenerator:AddInitialSplitting c->c,gamma; QtoQGammaSudakov do SplittingGenerator:AddInitialSplitting b->b,gamma; QtoQGammaSudakov do SplittingGenerator:AddInitialSplitting t->t,gamma; QtoQGammaSudakov do SplittingGenerator:AddInitialSplitting g->g,g; GtoGGSudakov # do SplittingGenerator:AddInitialSplitting g->d,dbar; GtoQQbarSudakov do SplittingGenerator:AddInitialSplitting g->u,ubar; GtoQQbarSudakov do SplittingGenerator:AddInitialSplitting g->s,sbar; GtoQQbarSudakov do SplittingGenerator:AddInitialSplitting g->c,cbar; GtoccbarSudakov do SplittingGenerator:AddInitialSplitting g->b,bbar; GtobbbarSudakov # do SplittingGenerator:AddInitialSplitting gamma->d,dbar; GammatoQQbarSudakov do SplittingGenerator:AddInitialSplitting gamma->u,ubar; GammatoQQbarSudakov do SplittingGenerator:AddInitialSplitting gamma->s,sbar; GammatoQQbarSudakov do SplittingGenerator:AddInitialSplitting gamma->c,cbar; GammatoQQbarSudakov do SplittingGenerator:AddInitialSplitting gamma->b,bbar; GammatoQQbarSudakov # do SplittingGenerator:AddInitialSplitting d->g,d; dtoGdSudakov do SplittingGenerator:AddInitialSplitting u->g,u; utoGuSudakov do SplittingGenerator:AddInitialSplitting s->g,s; QtoGQSudakov do SplittingGenerator:AddInitialSplitting c->g,c; QtoGQSudakov do SplittingGenerator:AddInitialSplitting b->g,b; QtoGQSudakov do SplittingGenerator:AddInitialSplitting dbar->g,dbar; dtoGdSudakov do SplittingGenerator:AddInitialSplitting ubar->g,ubar; utoGuSudakov do SplittingGenerator:AddInitialSplitting sbar->g,sbar; QtoGQSudakov do SplittingGenerator:AddInitialSplitting cbar->g,cbar; QtoGQSudakov do SplittingGenerator:AddInitialSplitting bbar->g,bbar; QtoGQSudakov # do SplittingGenerator:AddInitialSplitting d->gamma,d; QtoGammaQSudakov do SplittingGenerator:AddInitialSplitting u->gamma,u; QtoGammaQSudakov do SplittingGenerator:AddInitialSplitting s->gamma,s; QtoGammaQSudakov do SplittingGenerator:AddInitialSplitting c->gamma,c; QtoGammaQSudakov do SplittingGenerator:AddInitialSplitting b->gamma,b; QtoGammaQSudakov do SplittingGenerator:AddInitialSplitting dbar->gamma,dbar; QtoGammaQSudakov do SplittingGenerator:AddInitialSplitting ubar->gamma,ubar; QtoGammaQSudakov do SplittingGenerator:AddInitialSplitting sbar->gamma,sbar; QtoGammaQSudakov do SplittingGenerator:AddInitialSplitting cbar->gamma,cbar; QtoGammaQSudakov do SplittingGenerator:AddInitialSplitting bbar->gamma,bbar; QtoGammaQSudakov