No OneTemporary
Actions

Size

42 KB

Subscribers

None

View Options

	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<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
	+ 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<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]),"
	" if negative it is solved for at initialisation such that alpha_S(Qmin)=AlphaMaxNP",
	&ShowerAlphaQCD::_qmin, GeV, 0.630882GeV, 0.330445GeV,
	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 Deleted<ShowerAlphaQCD> delInputOption
	+ ("InputOption", "The old default (1) is now the only choice");

	- static Parameter<ShowerAlphaQCD,Energy> interfaceLambdaQCD
	- ("LambdaQCD",
	- "Input value of Lambda_MSBar",
	- &ShowerAlphaQCD::_lambdain, MeV, 0.208364GeV, 100.0MeV, 500.0*MeV,
	- false, false, Interface::limited);
	+ static Deleted<ShowerAlphaQCD> delAlphaMZ
	+ ("AlphaMZ", "Renamed to AlphaIn.");

	static Parameter<ShowerAlphaQCD,double> 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<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);

	+
	+
	+
	static Parameter<ShowerAlphaQCD,Energy> 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<ShowerAlphaQCD,Energy> interfaceQuarkMasses
	("QuarkMasses",
	"The quark masses to be used instead of the masses set in the particle data.",
	&ShowerAlphaQCD::_quarkMasses, GeV, -1, 0.0GeV, 0.0GeV, 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<ZERO ) {
	Energy low = _lambda[0]+MeV;
	if(value(sqr(low))<_asMaxNP)
	Throw<InitException>() << "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<as) low = _qmin;
	}
	while ( abs(as-_asMaxNP) > _tolerance );
	}
	// final threshold is qmin
	_thresholds[0]=_qmin;
	// value of alphaS at threshold
	pair<short,Energy> 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<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;

	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_lowGeV + kstep_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<short,Energy> 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()factorsqrt(scale);
	double val(0.);
	// normal case
	if (q >= _qmin) {
	pair<short,Energy> 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<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);
	+ xtest += (alpha-alphaS(match,lamtest,nflav,_nloop))/derivativeAlphaS(match,lamtest,nflav,_nloop);
	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);
	+ 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/(b0sqr(lx));
	- else if(_nloop==2)
	- return -4.pi/(b0sqr(lx))(1.+2.b1/sqr(b0)/lx(1.-2.log(lx)));
	- else
	- return -4.pi/(b0sqr(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)+b2b0/(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/(b0lx);}
	- // two loop
	- else if(_nloop==2) {
	- return 4.pi/(b0lx)(1.-2.b1/sqr(b0)*log(lx)/lx);
	- }
	- // three loop
	- else
	- {return 4.pi/(b0lx)(1.-2.b1/sqr(b0)*log(lx)/lx +
	- 4.sqr(b1)/(sqr(sqr(b0))sqr(lx))*
	- (sqr(log(lx) - 0.5) + b2b0/(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<Energy>& 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<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 &);
	+ 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<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;

	/**
	* 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<Energy> _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/(b0sqr(lx));
	+ else if(nloop==2)
	+ return -4.pi/(b0sqr(lx))(1.+2.b1/sqr(b0)/lx(1.-2.log(lx)));
	+ else
	+ return -4.pi/(b0sqr(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)+b2b0/(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/(b0lx);}
	+ // two loop
	+ else if(nloop==2) {
	+ return 4.pi/(b0lx)(1.-2.b1/sqr(b0)*log(lx)/lx);
	+ }
	+ // three loop
	+ else
	+ {return 4.pi/(b0lx)(1.-2.b1/sqr(b0)*log(lx)/lx +
	+ 4.sqr(b1)/(sqr(sqr(b0))sqr(lx))*
	+ (sqr(log(lx) - 0.5) + b2b0/(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

File Metadata

Mime Type: text/x-diff
Expires: Sat, May 3, 5:54 AM (17 h, 24 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4982858
Default Alt Text: (42 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions