Page MenuHomeHEPForge
Files F7877106

No OneTemporary
Actions

View Options

diff --git a/Shower/Couplings/ShowerAlphaQCD.cc b/Shower/Couplings/ShowerAlphaQCD.cc
--- a/Shower/Couplings/ShowerAlphaQCD.cc
+++ b/Shower/Couplings/ShowerAlphaQCD.cc
@@ -1,372 +1,412 @@
// -*- C++ -*-
//
// ShowerAlphaQCD.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 The Herwig Collaboration
//
// Herwig is licenced under version 2 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/Command.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"
using namespace Herwig;
DescribeClass<ShowerAlphaQCD,ShowerAlpha>
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
<< ounit(_thresholds,GeV) << ounit(_lambda,GeV);
}
void ShowerAlphaQCD::persistentInput(PersistentIStream & is, int) {
is >> _asType >> _asMaxNP >> iunit(_qmin,GeV) >> _nloop >> _lambdaopt >> _thresopt
>> iunit(_lambdain,GeV) >> _alphain >> _inopt
>> _tolerance >> _maxtry >> _alphamin
>> iunit(_thresholds,GeV) >> iunit(_lambda,GeV);
}
void ShowerAlphaQCD::Init() {
static ClassDocumentation<ShowerAlphaQCD> documentation
("This (concrete) class describes the QCD alpha running.");
static Switch<ShowerAlphaQCD, int> 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<ShowerAlphaQCD,Energy> intQmin
("Qmin", "Q < Qmin is treated with NP parametrization as of (unit [GeV])",
&ShowerAlphaQCD::_qmin, GeV, 0.630882*GeV, 0.330445*GeV,
100.0*GeV,false,false,false);
static Parameter<ShowerAlphaQCD,double> 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<ShowerAlphaQCD,unsigned int> interfaceNumberOfLoops
("NumberOfLoops",
"The number of loops to use in the alpha_S calculation",
&ShowerAlphaQCD::_nloop, 3, 1, 3,
false, false, Interface::limited);
static Switch<ShowerAlphaQCD,bool> 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 Parameter<ShowerAlphaQCD,Energy> interfaceLambdaQCD
("LambdaQCD",
"Input value of Lambda_MSBar",
&ShowerAlphaQCD::_lambdain, MeV, 0.208364*GeV, 100.0*MeV, 500.0*MeV,
false, false, Interface::limited);
static Parameter<ShowerAlphaQCD,double> interfaceAlphaMZ
("AlphaMZ",
"The input value of the strong coupling at the Z mass ",
&ShowerAlphaQCD::_alphain, 0.118, 0.1, 0.2,
false, false, Interface::limited);
static Switch<ShowerAlphaQCD,bool> 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<ShowerAlphaQCD,double> interfaceTolerance
("Tolerance",
"The tolerance for discontinuities in alphaS at thresholds.",
&ShowerAlphaQCD::_tolerance, 1e-10, 1e-20, 1e-4,
false, false, Interface::limited);
static Parameter<ShowerAlphaQCD,unsigned int> 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<ShowerAlphaQCD,bool> 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<ShowerAlphaQCD> interfaceValue
("Value",
"",
&ShowerAlphaQCD::value, false);
+ static Command<ShowerAlphaQCD> interfacecheck
+ ("check",
+ "check",
+ &ShowerAlphaQCD::check, false);
+
+
}
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(getParticleData(ParticleID::Z0)->mass(),_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.);
// compute the threshold matching
// top threshold
for(int ix=1;ix<4;++ix) {
if(_thresopt)
_thresholds[ix]=getParticleData(ix+3)->mass();
else
_thresholds[ix]=getParticleData(ix+3)->constituentMass();
}
// compute 6 flavour lambda by matching at top mass using Newton Raphson
_lambda[3]=computeLambda(_thresholds[3],alphaS(_thresholds[3],_lambda[2],5),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);
// 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);
// final threshold is qmin
_thresholds[0]=_qmin;
// compute the maximum value of as
if ( _asType < 5 ) _alphamin = value(sqr(_qmin)+1.0e-8*sqr(MeV)); // approx as = 1
else _alphamin = max(_asMaxNP, value(sqr(_qmin)+1.0e-8*sqr(MeV)));
// check consistency lambda_3 < qmin
if(_lambda[0]>_qmin)
Throw<InitException>() << "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;
+
+ cout << "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 {
pair<short,Energy> nflam;
Energy q = sqrt(scale);
double val(0.);
// special handling if the scale is less than Qmin
if (q < _qmin) {
nflam = getLamNfTwoLoop(_qmin);
double val0 = alphaS(_qmin, nflam.second, nflam.first);
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;
}
} else {
// the 'ordinary' case
nflam = getLamNfTwoLoop(q);
val = alphaS(q, nflam.second, nflam.first);
}
return scaleFactor() * val;
}
double ShowerAlphaQCD::overestimateValue() const {
return scaleFactor() * _alphamin;
}
double ShowerAlphaQCD::ratio(const Energy2 scale) const {
pair<short,Energy> nflam;
Energy q = sqrt(scale);
double val(0.);
// special handling if the scale is less than Qmin
if (q < _qmin) {
nflam = getLamNfTwoLoop(_qmin);
double val0 = alphaS(_qmin, nflam.second, nflam.first);
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;
}
} else {
// the 'ordinary' case
nflam = getLamNfTwoLoop(q);
val = alphaS(q, nflam.second, nflam.first);
}
// 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<short, Energy> 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<short,Energy>(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);
Energy newLambda = match/exp(0.5*xtest);
lamtest = newLambda<match ? newLambda : 0.5*(lamtest+match);
}
while(abs(alpha-alphaS(match,lamtest,nflav)) > _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/Couplings/ShowerAlphaQCD.h b/Shower/Couplings/ShowerAlphaQCD.h
--- a/Shower/Couplings/ShowerAlphaQCD.h
+++ b/Shower/Couplings/ShowerAlphaQCD.h
@@ -1,278 +1,285 @@
// -*- C++ -*-
//
// ShowerAlphaQCD.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 The Herwig Collaboration
//
// Herwig is licenced under version 2 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 "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.) {}
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) 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];
}
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<short, Energy> 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 &);
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<Energy> _thresholds;
/**
* \f$\Lambda\f$ for the different number of flavours
*/
vector<Energy> _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;
};
}
#endif /* HERWIG_ShowerAlphaQCD_H */

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 2:50 PM (1 d, 14 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3804825
Default Alt Text
(20 KB)

Event Timeline