No OneTemporary
Actions

Size

29 KB

Subscribers

None

View Options

	diff --git a/Shower/Dipole/Kernels/DipoleSplittingKernel.cc b/Shower/Dipole/Kernels/DipoleSplittingKernel.cc
	--- a/Shower/Dipole/Kernels/DipoleSplittingKernel.cc
	+++ b/Shower/Dipole/Kernels/DipoleSplittingKernel.cc
	@@ -1,408 +1,411 @@
	// -- C++ --
	//
	// DipoleSplittingKernel.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 DipoleSplittingKernel class.
	//

	#include "DipoleSplittingKernel.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Switch.h"

	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"

	#include "Herwig/Shower/ShowerHandler.h"

	using namespace Herwig;

	DipoleSplittingKernel::DipoleSplittingKernel()
	: HandlerBase(), theScreeningScale(0.0*GeV),
	thePresamplingPoints(2000), theMaxtry(100000),
	theFreezeGrid(500000),
	theDetuning(1.0),
	theStrictLargeN(false),
	theFactorizationScaleFactor(1.0),
	theRenormalizationScaleFactor(1.0),
	theRenormalizationScaleFreeze(1.*GeV),
	theFactorizationScaleFreeze(1.*GeV),
	theVirtualitySplittingScale(false),
	theCMWScheme(0),
	presampling(false) {}

	DipoleSplittingKernel::~DipoleSplittingKernel() {}


	+// initialize static variable out of line
	+double DipoleSplittingKernel::theMaxPDFRatio = 1000000.;
	+
	// If needed, insert default implementations of virtual function defined
	// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).


	void DipoleSplittingKernel::persistentOutput(PersistentOStream & os) const {
	os << theAlphaS << ounit(theScreeningScale,GeV)
	<< theSplittingKinematics << thePDFRatio
	<< thePresamplingPoints << theMaxtry << theFreezeGrid << theDetuning
	<< theFlavour << theMCCheck << theStrictLargeN
	<< theFactorizationScaleFactor
	<< theRenormalizationScaleFactor
	<< ounit(theRenormalizationScaleFreeze,GeV)
	<< ounit(theFactorizationScaleFreeze,GeV)
	<< theVirtualitySplittingScale<<theCMWScheme<<theUseThisKernel;
	}

	void DipoleSplittingKernel::persistentInput(PersistentIStream & is, int) {
	is >> theAlphaS >> iunit(theScreeningScale,GeV)
	>> theSplittingKinematics >> thePDFRatio
	>> thePresamplingPoints >> theMaxtry >> theFreezeGrid >> theDetuning
	>> theFlavour >> theMCCheck >> theStrictLargeN
	>> theFactorizationScaleFactor
	>> theRenormalizationScaleFactor
	>> iunit(theRenormalizationScaleFreeze,GeV)
	>> iunit(theFactorizationScaleFreeze,GeV)
	>> theVirtualitySplittingScale>>theCMWScheme>>theUseThisKernel;
	}

	double DipoleSplittingKernel::alphaPDF(const DipoleSplittingInfo& split,
	Energy optScale,
	double rScaleFactor,
	double fScaleFactor) const {

	Energy pt = optScale == ZERO ? split.lastPt() : optScale;

	Energy2 scale = ZERO;
	if ( !virtualitySplittingScale() ) {
	scale = sqr(pt) + sqr(theScreeningScale);
	} else {
	scale = sqr(splittingKinematics()->QFromPt(pt,split)) + sqr(theScreeningScale);
	}



	Energy2 fScale = sqr(theFactorizationScaleFactorfScaleFactor)scale;
	fScale = max( fScale , sqr(factorizationScaleFreeze()) );

	Energy2 rScale = sqr(theRenormalizationScaleFactorrScaleFactor)scale;
	rScale = max( rScale , sqr(renormalizationScaleFreeze()) );

	if(split.calcFixedExpansion()){
	fScale = max( sqr(split.fixedScale()) , sqr(factorizationScaleFreeze()) );
	rScale = max( sqr(split.fixedScale()) , sqr(renormalizationScaleFreeze()) );
	}

	double alphas = 1.0;
	double pdf = 1.0;

	// check if we are potentially reweighting and cache evaluations
	bool evaluatePDF = true;
	bool evaluateAlphaS = true;
	bool variations =
	!ShowerHandler::currentHandler()->showerVariations().empty() &&
	!presampling;
	if ( variations ) {

	map<double,double>::const_iterator pit = thePDFCache.find(fScaleFactor);
	evaluatePDF = (pit == thePDFCache.end());
	if ( !evaluatePDF ) {
	pdf = pit->second;
	}
	map<double,double>::const_iterator ait = theAlphaSCache.find(rScaleFactor);
	evaluateAlphaS = (ait == theAlphaSCache.end());
	if ( !evaluateAlphaS ) {
	alphas = ait->second;
	}
	}

	if ( evaluateAlphaS ){
	if (theCMWScheme==0\|\|split.calcFixedExpansion()) {
	alphas = alphaS()->value(rScale);
	}else if(theCMWScheme==1){

	alphas = alphaS()->value(rScale);
	alphas =1.+(3.(67./18.-1./6.*sqr(Constants::pi))
	-5./9.alphaS()->Nf(rScale))
	alphas/2./Constants::pi;

	}else if(theCMWScheme==2){
	double kg=exp(-(67.-3.sqr(Constants::pi)-10/3alphaS()->Nf(rScale))
	/(33.-2.*alphaS()->Nf(rScale)));
	Energy2 cmwscale2=max(kg*rScale, sqr(renormalizationScaleFreeze()) );
	alphas = alphaS()->value(cmwscale2);

	}else{
	throw Exception()
	<< "This CMW-Scheme is not implemented."
	<< Exception::abortnow;

	}
	}
	if ( evaluatePDF ) {
	if ( split.index().initialStateEmitter() ) {
	assert(pdfRatio());
	pdf *=
	split.lastEmitterZ() *
	(*pdfRatio())(split.index().emitterPDF(), fScale,
	split.index().emitterData(),split.emitterData(),
	split.emitterX(),split.lastEmitterZ());
	}

	if ( split.index().initialStateSpectator() ) {
	assert(pdfRatio());
	pdf *=
	split.lastSpectatorZ() *
	(*pdfRatio())(split.index().spectatorPDF(), fScale,
	split.index().spectatorData(),split.spectatorData(),
	split.spectatorX(),split.lastSpectatorZ());
	}
	}

	if ( evaluatePDF && variations ) {
	- thePDFCache[fScaleFactor] = pdf;
	+ thePDFCache[fScaleFactor] = min(pdf,theMaxPDFRatio);
	}

	if ( evaluateAlphaS && variations ) {
	theAlphaSCache[rScaleFactor] = alphas;
	}

	- double ret = pdf*
	+ double ret = min(pdf,theMaxPDFRatio)*
	(split.calcFixedExpansion()?
	1.:(alphas / (2.*Constants::pi)));

	if ( ret < 0. )
	ret = 0.;

	return ret;

	}

	void DipoleSplittingKernel::accept(const DipoleSplittingInfo& split,
	double, double,
	map<string,double>& weights) const {
	if ( ShowerHandler::currentHandler()->showerVariations().empty() )
	return;
	double reference = alphaPDF(split);
	assert(reference > 0.);
	for ( map<string,ShowerVariation>::const_iterator var =
	ShowerHandler::currentHandler()->showerVariations().begin();
	var != ShowerHandler::currentHandler()->showerVariations().end(); ++var ) {
	if ( ( ShowerHandler::currentHandler()->firstInteraction()
	&& var->second.firstInteraction ) \|\|
	( !ShowerHandler::currentHandler()->firstInteraction()
	&& var->second.secondaryInteractions ) ) {
	double varied = alphaPDF(split,ZERO,
	var->second.renormalizationScaleFactor,
	var->second.factorizationScaleFactor);
	if ( varied != reference ) {
	map<string,double>::iterator wi = weights.find(var->first);
	if ( wi != weights.end() )
	wi->second *= varied/reference;
	else
	weights[var->first] = varied/reference;
	}
	}
	}
	}

	void DipoleSplittingKernel::veto(const DipoleSplittingInfo& split,
	double p, double r,
	map<string,double>& weights) const {
	if ( ShowerHandler::currentHandler()->showerVariations().empty() )
	return;
	double reference = alphaPDF(split);
	// this is dangerous, but we have no other choice currently -- need to
	// carefully check for the effects; the assumption is that if the central
	// one ius zero, then so will be the variations.
	if ( reference == 0.0 )
	return;
	for ( map<string,ShowerVariation>::const_iterator var =
	ShowerHandler::currentHandler()->showerVariations().begin();
	var != ShowerHandler::currentHandler()->showerVariations().end(); ++var ) {
	if ( ( ShowerHandler::currentHandler()->firstInteraction()
	&& var->second.firstInteraction ) \|\|
	( !ShowerHandler::currentHandler()->firstInteraction()
	&& var->second.secondaryInteractions ) ) {
	double varied = alphaPDF(split,ZERO,
	var->second.renormalizationScaleFactor,
	var->second.factorizationScaleFactor);
	if ( varied != reference ) {
	map<string,double>::iterator wi = weights.find(var->first);
	if ( wi != weights.end() )
	wi->second = (r - variedp/reference) / (r-p);
	else
	weights[var->first] = (r - varied*p/reference) / (r-p);
	}
	}
	}
	}

	AbstractClassDescription<DipoleSplittingKernel>
	DipoleSplittingKernel::initDipoleSplittingKernel;
	// Definition of the static class description member.

	void DipoleSplittingKernel::Init() {

	static ClassDocumentation<DipoleSplittingKernel> documentation
	("DipoleSplittingKernel is the base class for all kernels "
	"used within the dipole shower.");

	static Reference<DipoleSplittingKernel,AlphaSBase> interfaceAlphaS
	("AlphaS",
	"The strong coupling to be used by this splitting kernel.",
	&DipoleSplittingKernel::theAlphaS, false, false, true, true, false);


	static Parameter<DipoleSplittingKernel,Energy> interfaceScreeningScale
	("ScreeningScale",
	"A colour screening scale",
	&DipoleSplittingKernel::theScreeningScale, GeV, 0.0GeV, 0.0GeV, 0*GeV,
	false, false, Interface::lowerlim);


	static Reference<DipoleSplittingKernel,DipoleSplittingKinematics>
	interfaceSplittingKinematics
	("SplittingKinematics",
	"The splitting kinematics to be used by this splitting kernel.",
	&DipoleSplittingKernel::theSplittingKinematics, false, false, true, false, false);


	static Reference<DipoleSplittingKernel,PDFRatio> interfacePDFRatio
	("PDFRatio",
	"Set the optional PDF ratio object to evaluate this kernel",
	&DipoleSplittingKernel::thePDFRatio, false, false, true, true, false);

	static Parameter<DipoleSplittingKernel,unsigned long> interfacePresamplingPoints
	("PresamplingPoints",
	"The number of points used to presample this kernel.",
	&DipoleSplittingKernel::thePresamplingPoints, 2000, 1, 0,
	false, false, Interface::lowerlim);

	static Parameter<DipoleSplittingKernel,unsigned long> interfaceMaxtry
	("Maxtry",
	"The maximum number of attempts to generate a splitting.",
	&DipoleSplittingKernel::theMaxtry, 10000, 1, 0,
	false, false, Interface::lowerlim);

	static Parameter<DipoleSplittingKernel,unsigned long> interfaceFreezeGrid
	("FreezeGrid",
	"",
	&DipoleSplittingKernel::theFreezeGrid, 500000, 1, 0,
	false, false, Interface::lowerlim);

	static Reference<DipoleSplittingKernel,ParticleData> interfaceFlavour
	("Flavour",
	"Set the flavour to be produced if ambiguous.",
	&DipoleSplittingKernel::theFlavour, false, false, true, true, false);

	static Reference<DipoleSplittingKernel,DipoleMCCheck> interfaceMCCheck
	("MCCheck",
	"[debug option] MCCheck",
	&DipoleSplittingKernel::theMCCheck, false, false, true, true, false);

	interfaceMCCheck.rank(-1);

	static Switch<DipoleSplittingKernel,bool> interfaceStrictLargeN
	("StrictLargeN",
	"Work in a strict large-N limit.",
	&DipoleSplittingKernel::theStrictLargeN, false, false, false);
	static SwitchOption interfaceStrictLargeNYes
	(interfaceStrictLargeN,
	"Yes",
	"Replace C_F -> C_A/2 where present",
	true);
	static SwitchOption interfaceStrictLargeNNo
	(interfaceStrictLargeN,
	"No",
	"Keep C_F=4/3",
	false);

	interfaceStrictLargeN.rank(-2);

	static Switch<DipoleSplittingKernel,unsigned int> interfaceCMWScheme
	("CMWScheme",
	"Use the CMW Scheme related Kg expression to the splitting",
	&DipoleSplittingKernel::theCMWScheme, 0, false, false);
	static SwitchOption interfaceCMWSchemeNo
	(interfaceCMWScheme,"No","No CMW-Scheme", 0);
	static SwitchOption interfaceCMWSchemeLinear
	(interfaceCMWScheme,"Linear",
	"Linear CMW multiplication: alpha_s(q) -> alpha_s(q)(1+K_g*alpha_s(q)/2pi )",1);
	static SwitchOption interfaceCMWSchemeFactor
	(interfaceCMWScheme,"Factor",
	"Use factor in alpha_s argument: alpha_s(q) -> alpha_s(k_gq) with kfac=exp(-(67-3pi^2-10/3Nf)/(33-2Nf)) ",2);

	static Parameter<DipoleSplittingKernel,double> interfaceFactorizationScaleFactor
	("FactorizationScaleFactor",
	"The factorization scale factor.",
	&DipoleSplittingKernel::theFactorizationScaleFactor, 1.0, 0.0, 0,
	false, false, Interface::lowerlim);

	interfaceFactorizationScaleFactor.rank(-2);

	static Parameter<DipoleSplittingKernel,double> interfaceRenormalizationScaleFactor
	("RenormalizationScaleFactor",
	"The renormalization scale factor.",
	&DipoleSplittingKernel::theRenormalizationScaleFactor, 1.0, 0.0, 0,
	false, false, Interface::lowerlim);

	interfaceRenormalizationScaleFactor.rank(-2);

	static Parameter<DipoleSplittingKernel,Energy> interfaceRenormalizationScaleFreeze
	("RenormalizationScaleFreeze",
	"The freezing scale for the renormalization scale.",
	&DipoleSplittingKernel::theRenormalizationScaleFreeze,
	GeV, 1.0GeV, 0.0GeV, 0*GeV,
	false, false, Interface::lowerlim);

	static Parameter<DipoleSplittingKernel,Energy> interfaceFactorizationScaleFreeze
	("FactorizationScaleFreeze",
	"The freezing scale for the factorization scale.",
	&DipoleSplittingKernel::theFactorizationScaleFreeze, GeV, 1.0GeV, 0.0GeV, 0*GeV,
	false, false, Interface::lowerlim);

	static Switch<DipoleSplittingKernel,bool> interfaceVirtualitySplittingScale
	("VirtualitySplittingScale",
	"Use the virtuality as the splitting scale.",
	&DipoleSplittingKernel::theVirtualitySplittingScale, false, false, false);
	static SwitchOption interfaceVirtualitySplittingScaleYes
	(interfaceVirtualitySplittingScale,
	"Yes",
	"Use vrituality.",
	true);
	static SwitchOption interfaceVirtualitySplittingScaleNo
	(interfaceVirtualitySplittingScale,
	"No",
	"Use transverse momentum.",
	false);

	static Parameter<DipoleSplittingKernel,double> interfaceDetuning
	("Detuning",
	"A value to detune the overestimate kernel.",
	&DipoleSplittingKernel::theDetuning, 1.0, 1.0, 0,
	false, false, Interface::lowerlim);


	static Switch<DipoleSplittingKernel,bool> interfaceUseThisKernel
	("UseKernel",
	"Turn On and of the Kernel.",
	&DipoleSplittingKernel::theUseThisKernel, true, false, false);
	static SwitchOption interfaceUseThisKernelYes
	(interfaceUseThisKernel,
	"Yes",
	"Use this Kernel.",
	true);
	static SwitchOption interfaceUseThisKernelNo
	(interfaceUseThisKernel,
	"No",
	"Dont use this Kernel.",
	false);


	}

	diff --git a/Shower/Dipole/Kernels/DipoleSplittingKernel.h b/Shower/Dipole/Kernels/DipoleSplittingKernel.h
	--- a/Shower/Dipole/Kernels/DipoleSplittingKernel.h
	+++ b/Shower/Dipole/Kernels/DipoleSplittingKernel.h
	@@ -1,535 +1,541 @@
	// -- C++ --
	//
	// DipoleSplittingKernel.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_DipoleSplittingKernel_H
	#define HERWIG_DipoleSplittingKernel_H
	//
	// This is the declaration of the DipoleSplittingKernel class.
	//

	#include "ThePEG/Handlers/HandlerBase.h"
	#include "ThePEG/StandardModel/AlphaSBase.h"
	#include "ThePEG/PDF/PDF.h"

	#include "Herwig/Shower/Dipole/Utility/PDFRatio.h"
	#include "Herwig/Shower/Dipole/Base/DipoleSplittingInfo.h"
	#include "Herwig/Shower/Dipole/Kinematics/DipoleSplittingKinematics.h"

	#include "ThePEG/EventRecord/RhoDMatrix.h"
	#include "Herwig/Decay/DecayMatrixElement.h"
	#include "Herwig/Decay/TwoBodyDecayMatrixElement.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* \ingroup DipoleShower
	* \author Simon Platzer
	*
	* \brief DipoleSplittingKernel is the base class for all kernels
	* used within the dipole shower.
	*
	* @see \ref DipoleSplittingKernelInterfaces "The interfaces"
	* defined for DipoleSplittingKernel.
	*/
	class DipoleSplittingKernel: public HandlerBase {

	public:

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* The default constructor.
	*/
	DipoleSplittingKernel();

	/**
	* The destructor.
	*/
	virtual ~DipoleSplittingKernel();
	//@}

	public:

	/**
	* Return the alpha_s to be used
	*/
	Ptr<AlphaSBase>::tptr alphaS() const { return theAlphaS; }

	/**
	* Set the alpha_s to be used
	*/
	void alphaS(Ptr<AlphaSBase>::tptr ap) { theAlphaS = ap; }

	/**
	* Return the splitting kinematics object
	*/
	Ptr<DipoleSplittingKinematics>::tptr splittingKinematics() const {
	return theSplittingKinematics;
	}

	/**
	* Return the mc check object
	*/
	Ptr<DipoleMCCheck>::ptr mcCheck() const { return theMCCheck; }

	/**
	* Set the splitting kinematics object
	*/
	void splittingKinematics(Ptr<DipoleSplittingKinematics>::tptr sp) {
	theSplittingKinematics = sp;
	}

	/**
	* Return the PDFRatio object
	*/
	Ptr<PDFRatio>::tptr pdfRatio() const { return thePDFRatio; }

	/**
	* Set the PDFRatio object
	*/
	void pdfRatio(Ptr<PDFRatio>::tptr sp) { thePDFRatio = sp; }

	/**
	* Return the number of additional parameter
	* random variables needed to evaluate this kernel
	* except the momentum fractions of incoming partons.
	* These will be accessible through the
	* lastSplittingParameters() container of the splitting
	* info object.
	*/
	virtual int nDimAdditional() const { return 0; }

	/**
	* Set the freezing value for the renormalization scale
	*/
	void renormalizationScaleFreeze(Energy s) { theRenormalizationScaleFreeze = s; }

	/**
	* Set the freezing value for the factorization scale
	*/
	void factorizationScaleFreeze(Energy s) { theFactorizationScaleFreeze = s; }

	/**
	* Get the freezing value for the renormalization scale
	*/
	Energy renormalizationScaleFreeze() const { return theRenormalizationScaleFreeze; }

	/**
	* Get the freezing value for the factorization scale
	*/
	Energy factorizationScaleFreeze() const { return theFactorizationScaleFreeze; }

	public:

	/**
	* Return true, if this splitting kernel
	* applies to the given dipole index.
	*/
	virtual bool canHandle(const DipoleIndex&) const = 0;

	/**
	* Return true, if this splitting kernel is
	* the same for the given index a, as the given
	* splitting kernel for index b.
	*/
	virtual bool canHandleEquivalent(const DipoleIndex& a,
	const DipoleSplittingKernel& sk,
	const DipoleIndex& b) const = 0;

	/**
	* Return the emitter data after splitting, given
	* a dipole index.
	*/
	virtual tcPDPtr emitter(const DipoleIndex&) const = 0;

	/**
	* Return the emission data after splitting, given
	* a dipole index.
	*/
	virtual tcPDPtr emission(const DipoleIndex&) const = 0;

	/**
	* Return the spectator data after splitting, given
	* a dipole index.
	*/
	virtual tcPDPtr spectator(const DipoleIndex&) const = 0;

	/**
	* Return the flavour produced, if this cannot
	* be determined from the dipole.
	*/
	PDPtr flavour() const { return theFlavour; }

	/**
	* Return true, if this splitting kernel is supposed to work in a
	* strict large-N limit, i.e. replacing C_F by C_A/2
	*/
	bool strictLargeN() const { return theStrictLargeN; }

	public:

	/**
	* Inform this splitting kernel, that it is being
	* presampled until a call to stopPresampling
	*/
	virtual void startPresampling(const DipoleIndex&) {
	presampling = true;
	}

	/**
	* Inform this splitting kernel, that it is not being
	* presampled until a call to startPresampling
	*/
	virtual void stopPresampling(const DipoleIndex&) {
	presampling = false;
	}

	/**
	* Return the number of points to presample this
	* splitting generator.
	*/
	unsigned long presamplingPoints() const { return thePresamplingPoints; }

	/**
	* Return the maximum number of trials
	* to generate a splitting.
	*/
	unsigned long maxtry() const { return theMaxtry; }

	/**
	* Return the number of accepted points after which the grid should
	* be frozen
	*/
	unsigned long freezeGrid() const { return theFreezeGrid; }

	/**
	* Set the number of accepted points after which the grid should
	* be frozen
	*/
	void freezeGrid(unsigned long n) { theFreezeGrid = n; }

	/**
	* Set a detuning factor to be applied to the sampling overestimate kernel
	*/
	void detuning(double d) { theDetuning = d; }

	/**
	* Return the detuning factor applied to the sampling overestimate kernel
	*/
	double detuning() const { return theDetuning; }

	/**
	* Evaluate this splitting kernel for the given
	* dipole splitting.
	*/
	virtual double evaluate(const DipoleSplittingInfo&) const = 0;

	/**
	* Evaluate rho_ii' V_ijk V*_i'jk / equivalent for initial-state splitting,
	* required for generating spin-correlated azimuthal angles.
	**/
	virtual vector< pair<int, Complex> > generatePhi( const DipoleSplittingInfo& dInfo, const RhoDMatrix& rho) const = 0;

	/**
	* Return the completely spin-unaveraged (i.e. spin-indexed) splitting kernel.
	**/
	virtual DecayMEPtr matrixElement(const DipoleSplittingInfo& dInfo) const = 0;

	/**
	* Clear the alphaPDF cache
	*/
	void clearAlphaPDFCache() const {
	theAlphaSCache.clear();
	thePDFCache.clear();
	}

	/**
	* Update the variations vector at the given splitting using the indicated
	* kernel and overestimate values.
	*/
	virtual void accept(const DipoleSplittingInfo&, double, double, map<string,double>&) const;

	/**
	* Update the variations vector at the given splitting using the indicated
	* kernel and overestimate values.
	*/
	virtual void veto(const DipoleSplittingInfo&, double, double, map<string,double>&) const;

	/**
	* Return true, if this kernel is capable of
	* delivering an overestimate to the kernel, and
	* of inverting the integral over the overestimate
	* w.r.t. the phasepsace provided by the given
	* DipoleSplittingInfo object.
	*/
	virtual bool haveOverestimate(const DipoleSplittingInfo&) const { return false; }

	/**
	* Return the overestimate to this splitting kernel
	* for the given dipole splitting.
	*/
	virtual double overestimate(const DipoleSplittingInfo&) const { return -1.; }

	/**
	* Invert the integral over the overestimate
	* w.r.t. the phasepsace provided by the given
	* DipoleSplittingInfo object to equal
	* the given value.
	*/
	virtual double invertOverestimateIntegral(const DipoleSplittingInfo&, double) const {
	return -1.;
	}

	/**
	* .
	*/
	bool useThisKernel() const {
	return theUseThisKernel;
	}

	public:

	/**
	* Get the factorization scale factor
	*/
	double factorizationScaleFactor() const { return theFactorizationScaleFactor; }

	/**
	* Set the factorization scale factor
	*/
	void factorizationScaleFactor(double f) { theFactorizationScaleFactor = f; }

	/**
	* Get the renormalization scale factor
	*/
	double renormalizationScaleFactor() const { return theRenormalizationScaleFactor; }

	/**
	* Set the renormalization scale factor
	*/
	void renormalizationScaleFactor(double f) { theRenormalizationScaleFactor = f; }

	protected:

	/**
	* Return the common factor of (alphas/2pi)*(pdf ratio)
	*/
	double alphaPDF(const DipoleSplittingInfo&,
	Energy optScale = ZERO,
	double rScaleFactor = 1.0,
	double fScaleFactor = 1.0) const;

	/**
	* Return true, if the virtuality of the splitting should be used as the
	* argument of alphas rather than the pt
	*/
	bool virtualitySplittingScale() const { return theVirtualitySplittingScale; }

	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();


	// If needed, insert declarations of virtual function defined in the
	// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).

	private:

	/**
	* The alpha_s to be used.
	*/
	Ptr<AlphaSBase>::ptr theAlphaS;

	/**
	* An optional 'colour screening' scale
	* for alternative intrinsic pt generation.
	*/
	Energy theScreeningScale;

	/**
	* The splitting kinematics to be used.
	*/
	Ptr<DipoleSplittingKinematics>::ptr theSplittingKinematics;

	/**
	* An optional PDF ratio object to be used
	* when evaluating this kernel.
	*/
	Ptr<PDFRatio>::ptr thePDFRatio;

	/**
	* The number of points to presample this
	* splitting generator.
	*/
	unsigned long thePresamplingPoints;

	/**
	* The maximum number of trials
	* to generate a splitting.
	*/
	unsigned long theMaxtry;
	-
	+
	+ /**
	+ * The maximum value for any pdf ratio.
	+ * TODO: JB:Should this be an interfaced value? Is there a reasobable case where it should be allowed to be bigger than 1000000.?
	+ */
	+ static double theMaxPDFRatio;
	+
	/**
	* Return the number of accepted points after which the grid should
	* be frozen
	*/
	unsigned long theFreezeGrid;

	/**
	* The detuning factor applied to the sampling overestimate kernel
	*/
	double theDetuning;

	/**
	* The flavour produced, if this cannot
	* be determined from the dipole.
	*/
	PDPtr theFlavour;

	/**
	* Pointer to a check histogram object
	*/
	Ptr<DipoleMCCheck>::ptr theMCCheck;

	/**
	* True, if this splitting kernel is supposed to work in a
	* strict large-N limit, i.e. replacing C_F by C_A/2
	*/
	bool theStrictLargeN;

	/**
	* The factorization scale factor.
	*/
	double theFactorizationScaleFactor;

	/**
	* The renormalization scale factor.
	*/
	double theRenormalizationScaleFactor;

	/**
	* A freezing value for the renormalization scale
	*/
	Energy theRenormalizationScaleFreeze;

	/**
	* A freezing value for the factorization scale
	*/
	Energy theFactorizationScaleFreeze;

	/**
	* True, if the virtuality of the splitting should be used as the
	* argument of alphas rather than the pt
	*/
	bool theVirtualitySplittingScale;

	/**
	* Implementing CMW in the kernels.
	**/

	unsigned int theCMWScheme=0;

	/**
	* Cache for alphas evaluations
	*/
	mutable map<double,double> theAlphaSCache;

	/**
	* Cache for PDF evaluations
	*/
	mutable map<double,double> thePDFCache;

	/**
	* True, if we are presampling
	*/
	bool presampling;

	/**
	* True, if the kernel should be used
	*/
	bool theUseThisKernel = true;


	private:

	/**
	* The static object used to initialize the description of this class.
	* Indicates that this is an abstract class with persistent data.
	*/
	static AbstractClassDescription<DipoleSplittingKernel> initDipoleSplittingKernel;

	/**
	* The assignment operator is private and must never be called.
	* In fact, it should not even be implemented.
	*/
	DipoleSplittingKernel & operator=(const DipoleSplittingKernel &) = delete;

	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

	/** This template specialization informs ThePEG about the
	* base classes of DipoleSplittingKernel. */
	template <>
	struct BaseClassTrait<Herwig::DipoleSplittingKernel,1> {
	/** Typedef of the first base class of DipoleSplittingKernel. */
	typedef HandlerBase NthBase;
	};

	/** This template specialization informs ThePEG about the name of
	* the DipoleSplittingKernel class and the shared object where it is defined. */
	template <>
	struct ClassTraits<Herwig::DipoleSplittingKernel>
	: public ClassTraitsBase<Herwig::DipoleSplittingKernel> {
	/** Return a platform-independent class name */
	static string className() { return "Herwig::DipoleSplittingKernel"; }
	/**
	* The name of a file containing the dynamic library where the class
	* DipoleSplittingKernel is implemented. It may also include several, space-separated,
	* libraries if the class DipoleSplittingKernel depends on other classes (base classes
	* excepted). In this case the listed libraries will be dynamically
	* linked in the order they are specified.
	*/
	static string library() { return "HwDipoleShower.so"; }
	};

	/** @endcond */

	}

	#endif /* HERWIG_DipoleSplittingKernel_H */

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 3:59 PM (1 d, 15 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3799443
Default Alt Text: (29 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions