No OneTemporary
Actions

Size

76 KB

Subscribers

None

View Options

	diff --git a/DipoleShower/Base/DipoleSplittingGenerator.cc b/DipoleShower/Base/DipoleSplittingGenerator.cc
	--- a/DipoleShower/Base/DipoleSplittingGenerator.cc
	+++ b/DipoleShower/Base/DipoleSplittingGenerator.cc
	@@ -1,721 +1,724 @@
	// -- C++ --
	//
	// DipoleSplittingGenerator.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2007 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 DipoleSplittingGenerator class.
	//
	#include <config.h>
	#include "DipoleSplittingGenerator.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Repository/EventGenerator.h"

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

	#include "Herwig/DipoleShower/DipoleShowerHandler.h"

	using namespace Herwig;

	DipoleSplittingGenerator::DipoleSplittingGenerator()
	: HandlerBase(),
	theExponentialGenerator(0), prepared(false), presampling(false),
	theDoCompensate(false) {
	if ( ShowerHandler::currentHandler() )
	setGenerator(ShowerHandler::currentHandler()->generator());
	}

	DipoleSplittingGenerator::~DipoleSplittingGenerator() {
	if ( theExponentialGenerator ) {
	delete theExponentialGenerator;
	theExponentialGenerator = 0;
	}
	}

	IBPtr DipoleSplittingGenerator::clone() const {
	return new_ptr(*this);
	}

	IBPtr DipoleSplittingGenerator::fullclone() const {
	return new_ptr(*this);
	}

	void DipoleSplittingGenerator::wrap(Ptr<DipoleSplittingGenerator>::ptr other) {
	assert(!prepared);
	theOtherGenerator = other;
	}

	void DipoleSplittingGenerator::prepare(const DipoleSplittingInfo& sp) {

	generatedSplitting = sp;

	generatedSplitting.splittingKinematics(splittingKernel()->splittingKinematics());
	generatedSplitting.splittingParameters().resize(splittingKernel()->nDimAdditional());

	if ( wrapping() ) {
	generatedSplitting.emitterData(theSplittingKernel->emitter(generatedSplitting.index()));
	generatedSplitting.spectatorData(theSplittingKernel->spectator(generatedSplitting.index()));
	generatedSplitting.emissionData(theSplittingKernel->emission(generatedSplitting.index()));
	parameters.resize(theOtherGenerator->nDim());
	prepared = true;
	return;
	}

	generatedSplitting.emitterData(splittingKernel()->emitter(generatedSplitting.index()));
	generatedSplitting.spectatorData(splittingKernel()->spectator(generatedSplitting.index()));
	generatedSplitting.emissionData(splittingKernel()->emission(generatedSplitting.index()));

	presampledSplitting = generatedSplitting;

	prepared = true;

	parameters.resize(nDim());

	theExponentialGenerator =
	new exsample::exponential_generator<DipoleSplittingGenerator,UseRandom>();

	theExponentialGenerator->sampling_parameters().maxtry = maxtry();
	theExponentialGenerator->sampling_parameters().presampling_points = presamplingPoints();
	theExponentialGenerator->sampling_parameters().freeze_grid = freezeGrid();

	theExponentialGenerator->docompensate(theDoCompensate);
	theExponentialGenerator->function(this);
	theExponentialGenerator->initialize();

	}

	void DipoleSplittingGenerator::fixParameters(const DipoleSplittingInfo& sp,
	Energy optHardPt) {

	assert(generator());

	assert(!presampling);
	assert(prepared);

	assert(sp.index() == generatedSplitting.index());

	generatedSplitting.scale(sp.scale());
	parameters[3] = sp.scale()/generator()->maximumCMEnergy();

	generatedSplitting.hardPt(sp.hardPt());

	parameters[0] = splittingKinematics()->ptToRandom(optHardPt == ZERO ?
	generatedSplitting.hardPt() :
	min(generatedSplitting.hardPt(),optHardPt),
	sp.scale(),
	sp.emitterX(), sp.spectatorX(),
	generatedSplitting.index(),
	*splittingKernel());

	size_t shift = 4;

	if ( generatedSplitting.index().emitterPDF().pdf() &&
	generatedSplitting.index().spectatorPDF().pdf() ) {
	generatedSplitting.emitterX(sp.emitterX());
	generatedSplitting.spectatorX(sp.spectatorX());
	parameters[4] = sp.emitterX();
	parameters[5] = sp.spectatorX();
	shift += 2;
	}

	if ( generatedSplitting.index().emitterPDF().pdf() &&
	!generatedSplitting.index().spectatorPDF().pdf() ) {
	generatedSplitting.emitterX(sp.emitterX());
	parameters[4] = sp.emitterX();
	++shift;
	}

	if ( !generatedSplitting.index().emitterPDF().pdf() &&
	generatedSplitting.index().spectatorPDF().pdf() ) {
	generatedSplitting.spectatorX(sp.spectatorX());
	parameters[4] = sp.spectatorX();
	++shift;
	}

	if ( splittingReweight() ) {
	parameters[shift] = splittingReweight()->evaluate(sp);
	++shift;
	}

	if ( splittingKernel()->nDimAdditional() )
	copy(sp.lastSplittingParameters().begin(),sp.lastSplittingParameters().end(),parameters.begin()+shift);

	if ( sp.emitter() )
	generatedSplitting.emitter(sp.emitter());

	if ( sp.spectator() )
	generatedSplitting.spectator(sp.spectator());

	}

	int DipoleSplittingGenerator::nDim() const {

	assert(!wrapping());
	assert(prepared);

	int ret = 4; // 0 pt, 1 z, 2 phi, 3 scale, 4/5 xs + parameters

	if ( generatedSplitting.index().emitterPDF().pdf() ) {
	++ret;
	}

	if ( generatedSplitting.index().spectatorPDF().pdf() ) {
	++ret;
	}

	if ( splittingReweight() )
	++ret;

	ret += splittingKernel()->nDimAdditional();

	return ret;

	}

	const vector<bool>& DipoleSplittingGenerator::sampleFlags() {

	assert(!wrapping());

	if ( !theFlags.empty() )
	return theFlags;

	theFlags.resize(nDim(),false);
	theFlags[0] = true; theFlags[1] = true; theFlags[2] = true; // 0 pt, 1 z, 2 phi
	return theFlags;
	}

	const pair<vector<double>,vector<double> >& DipoleSplittingGenerator::support() {

	assert(!wrapping());

	if ( !theSupport.first.empty() )
	return theSupport;

	vector<double> lower(nDim(),0.);
	vector<double> upper(nDim(),1.);

	pair<double,double> kSupport =
	generatedSplitting.splittingKinematics()->kappaSupport(generatedSplitting);

	pair<double,double> xSupport =
	generatedSplitting.splittingKinematics()->xiSupport(generatedSplitting);

	lower[0] = kSupport.first;
	lower[1] = xSupport.first;

	upper[0] = kSupport.second;
	upper[1] = xSupport.second;

	if ( splittingReweight() ) {
	pair<double,double> bounds = splittingReweight()->reweightBounds(generatedSplitting.index());
	int pos = 4;
	if ( generatedSplitting.index().emitterPDF().pdf() ) {
	++pos;
	}
	if ( generatedSplitting.index().spectatorPDF().pdf() ) {
	++pos;
	}
	lower[pos] = bounds.first;
	upper[pos] = bounds.second;
	}

	theSupport.first = lower;
	theSupport.second = upper;

	return theSupport;

	}

	void DipoleSplittingGenerator::startPresampling() {
	assert(!wrapping());
	splittingKernel()->startPresampling(generatedSplitting.index());
	presampling = true;
	}

	void DipoleSplittingGenerator::stopPresampling() {
	assert(!wrapping());
	splittingKernel()->stopPresampling(generatedSplitting.index());
	presampling = false;
	}

	bool DipoleSplittingGenerator::haveOverestimate() const {

	assert(!wrapping());
	assert(prepared);

	return
	generatedSplitting.splittingKinematics()->haveOverestimate() &&
	splittingKernel()->haveOverestimate(generatedSplitting);

	}

	bool DipoleSplittingGenerator::overestimate(const vector<double>& point) {

	assert(!wrapping());
	assert(prepared);
	assert(!presampling);
	assert(haveOverestimate());

	if ( ! generatedSplitting.splittingKinematics()->generateSplitting(point[0],point[1],point[2],
	generatedSplitting,
	*splittingKernel()) )
	return 0.;

	generatedSplitting.splittingKinematics()->prepareSplitting(generatedSplitting);

	return
	( generatedSplitting.splittingKinematics()->jacobianOverestimate() *
	splittingKernel()->overestimate(generatedSplitting) *
	(splittingReweight() ? splittingReweight()->evaluate(generatedSplitting) : 1.) );

	}

	double DipoleSplittingGenerator::invertOverestimateIntegral(double value) const {

	assert(!wrapping());
	assert(prepared);
	assert(!presampling);
	assert(haveOverestimate());

	return
	splittingKernel()->invertOverestimateIntegral(generatedSplitting,value);

	}

	double DipoleSplittingGenerator::evaluate(const vector<double>& point,Energy fixed) {

	assert(!wrapping());
	assert(prepared);
	assert(generator());

	DipoleSplittingInfo& split =
	( !presampling ? generatedSplitting : presampledSplitting );
	if(fixed>0.*GeV){
	split.fixedScale(fixed);
	}else{
	split.fixedScale(-1.*GeV);
	}

	split.continuesEvolving();

	size_t shift = 4;

	if ( presampling ) {

	split.scale(point[3] * generator()->maximumCMEnergy());

	if ( split.index().emitterPDF().pdf() &&
	split.index().spectatorPDF().pdf() ) {
	split.emitterX(point[4]);
	split.spectatorX(point[5]);
	shift += 2;
	}

	if ( split.index().emitterPDF().pdf() &&
	!split.index().spectatorPDF().pdf() ) {
	split.emitterX(point[4]);
	++shift;
	}

	if ( !split.index().emitterPDF().pdf() &&
	split.index().spectatorPDF().pdf() ) {
	split.spectatorX(point[4]);
	++shift;
	}

	if ( splittingReweight() )
	++shift;

	if ( splittingKernel()->nDimAdditional() )
	copy(point.begin()+shift,point.end(),split.splittingParameters().begin());

	split.hardPt(split.splittingKinematics()->ptMax(split.scale(),
	split.emitterX(),
	split.spectatorX(),
	split.index(),
	*splittingKernel()));

	}

	if ( ! split.splittingKinematics()->generateSplitting(point[0],point[1],point[2],split,*splittingKernel()) ) {
	split.lastValue(0.);
	return 0.;
	}

	split.splittingKinematics()->prepareSplitting(split);

	if ( split.stoppedEvolving() ) {
	split.lastValue(0.);
	return 0.;
	}

	double kernel = splittingKernel()->evaluate(split);
	if ( splittingReweight() ) {
	if ( !presampling )
	kernel *= splittingReweight()->evaluate(split);
	else
	kernel *= point[shift-1];
	}
	double jac = split.splittingKinematics()->jacobian();

	split.lastValue( abs(jac) * kernel );

	if ( isnan(split.lastValue()) \|\| isinf(split.lastValue()) ) {
	generator()->log() << "DipoleSplittingGenerator:evaluate(): problematic splitting kernel encountered for "
	<< splittingKernel()->name() << "\n" << flush;
	split.lastValue(0.0);
	}

	if ( kernel < 0. )
	return 0.;

	return split.lastValue();

	}

	void DipoleSplittingGenerator::doGenerate(Energy optCutoff) {

	assert(!wrapping());

	double res = 0.;

	Energy startPt = generatedSplitting.hardPt();
	double optKappaCutoff = 0.0;
	if ( optCutoff > splittingKinematics()->IRCutoff() ) {
	optKappaCutoff = splittingKinematics()->ptToRandom(optCutoff,
	generatedSplitting.scale(),
	generatedSplitting.emitterX(),
	generatedSplitting.spectatorX(),
	generatedSplitting.index(),
	*splittingKernel());
	}

	while (true) {
	try {
	if ( optKappaCutoff == 0.0 ) {
	res = theExponentialGenerator->generate();
	} else {
	res = theExponentialGenerator->generate(optKappaCutoff);
	}
	} catch (exsample::exponential_regenerate&) {
	generatedSplitting.hardPt(startPt);
	continue;
	} catch (exsample::hit_and_miss_maxtry&) {
	throw DipoleShowerHandler::RedoShower();
	} catch (exsample::selection_maxtry&) {
	throw DipoleShowerHandler::RedoShower();
	}
	break;
	}

	if ( res == 0. ) {
	generatedSplitting.lastPt(0.0*GeV);
	generatedSplitting.didStopEvolving();
	} else {

	generatedSplitting.continuesEvolving();

	if ( theMCCheck )
	theMCCheck->book(generatedSplitting.emitterX(),
	generatedSplitting.spectatorX(),
	generatedSplitting.scale(),
	startPt,
	generatedSplitting.lastPt(),
	generatedSplitting.lastZ(),
	1.);

	}

	}

	Energy DipoleSplittingGenerator::generate(const DipoleSplittingInfo& split,
	Energy optHardPt,
	Energy optCutoff) {

	fixParameters(split,optHardPt);

	if ( wrapping() ) {
	return theOtherGenerator->generateWrapped(generatedSplitting,optHardPt,optCutoff);
	}

	doGenerate(optCutoff);

	return generatedSplitting.lastPt();

	}

	double DipoleSplittingGenerator::unlops(const DipoleSplittingInfo& split,Energy down,Energy fixedScale){
	fixParameters(split);
	if ( wrapping() ) {
	return theOtherGenerator->wrappedUnlops( generatedSplitting, down,fixedScale);
	}
	return dounlops( split, down,fixedScale);
	}

	double DipoleSplittingGenerator::sudakov(const DipoleSplittingInfo& split,Energy down, bool fast){
	fixParameters(split);
	if ( wrapping() ) {
	return theOtherGenerator->wrappedSudakov( generatedSplitting, down,fast);
	}
	return dosudakov( split, down,fast);
	}



	double DipoleSplittingGenerator::dounlops(const DipoleSplittingInfo& ,Energy down,Energy fixedScale){
	assert(down > splittingKinematics()->IRCutoff());
	double optKappaCutoff = splittingKinematics()->ptToRandom(down,
	generatedSplitting.scale(),
	generatedSplitting.emitterX(),
	generatedSplitting.spectatorX(),
	generatedSplitting.index(),
	*splittingKernel());
	vector<double> RN;
	RN.resize(parameters.size());
	double res=0.;
	double resq=0.;
	double var=10.;
	double varx=10.;
	int k=0;
	while ((k<500.\|\|varx>0.1)&&k<5000){
	k+=1.;
	RN[0]= optKappaCutoff+(1-optKappaCutoff)*UseRandom::rnd();
	for (size_t rn=1;rn< parameters.size();rn++)RN[rn]=UseRandom::rnd();
	double tmp=(1-optKappaCutoff)*evaluate(RN,fixedScale);
	res+= tmp;
	resq+=pow(tmp,2.);
	if(k%50==0.){
	varx=sqrt((resq/pow(1.k,2)-pow(res,2)/pow(1.k,3)));
	}
	}
	return -res/(1.0*k);
	}















	double DipoleSplittingGenerator::dosudakov(const DipoleSplittingInfo& ,Energy down,bool fast){

	double est=splittingKernel()->estimate(generatedSplitting.scale(),down);;
	if(fast){
	//cout<<"\nDipoleSplittingGenerator::dosudakov estimate "<< est;
	Ptr<StandardEventHandler>::ptr eH =
	dynamic_ptr_cast<Ptr<StandardEventHandler>::ptr>(generator()->eventHandler());
	eH->estimatedDSigRD();
	return est;
	}

	double optKappaCutoff = splittingKinematics()->ptToRandom(down,
	generatedSplitting.scale(),
	generatedSplitting.emitterX(),
	generatedSplitting.spectatorX(),
	generatedSplitting.index(),
	*splittingKernel());
	vector<double> RN;
	RN.resize(parameters.size());
	double res=0.;
	double resq=0.;
	double var=10.;
	double varx=10.;
	int k=0;

	while (((k<500.\|\|var>0.05)&&k<50000)\|\|(fast&&k<50)){
	k+=1.;
	RN[0]= optKappaCutoff+(1-optKappaCutoff)*UseRandom::rnd();
	for (size_t rn=1;rn< parameters.size();rn++)RN[rn]=UseRandom::rnd();
	double tmp=(1-optKappaCutoff)*evaluate(RN);
	res+= tmp;
	resq+=pow(tmp,2.);
	varx=sqrt((resq/pow(1.k,2)-pow(res,2)/pow(1.k,3)));
	if(k%50==0.)var= exp(-(res)/(1.0k)+varx)-exp(-(res)/(1.0k)-varx);

	}
	- // cout<<"\n"<<splittingKernel()->name()<<" "<<generatedSplitting.scale()/GeV<<" "<<down/GeV<<" "<<log(est)<<" "<<-res/(1.0k)<<" "<<log(est)/(-res/(1.0k));
	+ //if (abs(est/exp(-res/(1.0*k))-1.)>2) {
	+ cout<<"\n"<<splittingKernel()->name()<<" "<<generatedSplitting.scale()/GeV<<" "<<down/GeV<<" "<<est<<" "<<exp(-res/(1.0k))<<" "<<est/exp(-res/(1.0k));
	+ //}
	+
	return exp(-res/(1.0*k));

	}


	double DipoleSplittingGenerator::wrappedUnlops(DipoleSplittingInfo& split,
	Energy down,Energy fixedScale) {

	assert(!wrapping());

	DipoleSplittingInfo backup = generatedSplitting;
	generatedSplitting = split;

	fixParameters(split);
	double res=dounlops( split, down,fixedScale);

	split = generatedSplitting;
	generatedSplitting = backup;

	return res;

	}

	double DipoleSplittingGenerator::wrappedSudakov(DipoleSplittingInfo& split,
	Energy down, bool fast) {

	assert(!wrapping());

	DipoleSplittingInfo backup = generatedSplitting;
	generatedSplitting = split;

	fixParameters(split);
	double res=dosudakov( split, down,fast);

	split = generatedSplitting;
	generatedSplitting = backup;

	return res;

	}

	Energy DipoleSplittingGenerator::generateWrapped(DipoleSplittingInfo& split,
	Energy optHardPt,
	Energy optCutoff) {

	assert(!wrapping());

	DipoleSplittingInfo backup = generatedSplitting;
	generatedSplitting = split;

	fixParameters(split,optHardPt);

	try {
	doGenerate(optCutoff);
	} catch (...) {
	split = generatedSplitting;
	generatedSplitting = backup;
	throw;
	}

	Energy pt = generatedSplitting.lastPt();

	split = generatedSplitting;
	generatedSplitting = backup;

	return pt;

	}

	void DipoleSplittingGenerator::completeSplitting(DipoleSplittingInfo& sp) const {
	pair<bool,bool> conf = sp.configuration();
	sp = generatedSplitting;
	sp.configuration(conf);
	}

	Ptr<DipoleSplittingKernel>::tptr DipoleSplittingGenerator::splittingKernel() const {
	if ( wrapping() )
	return theOtherGenerator->splittingKernel();
	return theSplittingKernel;
	}

	Ptr<DipoleSplittingReweight>::tptr DipoleSplittingGenerator::splittingReweight() const {
	if ( wrapping() )
	return theOtherGenerator->splittingReweight();
	return theSplittingReweight;
	}

	Ptr<DipoleSplittingKinematics>::tptr DipoleSplittingGenerator::splittingKinematics() const {
	if ( wrapping() )
	return theOtherGenerator->splittingKinematics();
	return theSplittingKernel->splittingKinematics();
	}

	void DipoleSplittingGenerator::splittingKernel(Ptr<DipoleSplittingKernel>::tptr sp) {
	theSplittingKernel = sp;
	if ( theSplittingKernel->mcCheck() )
	theMCCheck = theSplittingKernel->mcCheck();
	}

	void DipoleSplittingGenerator::splittingReweight(Ptr<DipoleSplittingReweight>::tptr sp) {
	theSplittingReweight = sp;
	}

	void DipoleSplittingGenerator::debugGenerator(ostream& os) const {

	os << "--- DipoleSplittingGenerator ---------------------------------------------------\n";

	os << " generating splittings using\n"
	<< " splittingKernel = " << splittingKernel()->name()
	<< " splittingKinematics = " << generatedSplitting.splittingKinematics()->name() << "\n"
	<< " to sample splittings of type:\n";

	os << generatedSplitting;

	os << "--------------------------------------------------------------------------------\n";

	}

	void DipoleSplittingGenerator::debugLastEvent(ostream& os) const {

	os << "--- DipoleSplittingGenerator ---------------------------------------------------\n";

	os << " last generated event:\n";

	os << generatedSplitting;

	os << "--------------------------------------------------------------------------------\n";

	}

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


	void DipoleSplittingGenerator::persistentOutput(PersistentOStream & os) const {
	os << theOtherGenerator << theSplittingKernel << theSplittingReweight << theMCCheck << theDoCompensate;
	}

	void DipoleSplittingGenerator::persistentInput(PersistentIStream & is, int) {
	is >> theOtherGenerator >> theSplittingKernel >> theSplittingReweight >> theMCCheck >> theDoCompensate;
	}

	ClassDescription<DipoleSplittingGenerator> DipoleSplittingGenerator::initDipoleSplittingGenerator;
	// Definition of the static class description member.

	void DipoleSplittingGenerator::Init() {

	static ClassDocumentation<DipoleSplittingGenerator> documentation
	("DipoleSplittingGenerator is used by the dipole shower "
	"to sample splittings from a given dipole splitting kernel.");


	static Reference<DipoleSplittingGenerator,DipoleSplittingKernel> interfaceSplittingKernel
	("SplittingKernel",
	"Set the splitting kernel to sample from.",
	&DipoleSplittingGenerator::theSplittingKernel, false, false, true, false, false);

	static Reference<DipoleSplittingGenerator,DipoleSplittingReweight> interfaceSplittingReweight
	("SplittingReweight",
	"Set the splitting reweight.",
	&DipoleSplittingGenerator::theSplittingReweight, false, false, true, true, false);

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

	interfaceMCCheck.rank(-1);

	}

	diff --git a/DipoleShower/Kernels/FIgx2ggxDipoleKernel.h b/DipoleShower/Kernels/FIgx2ggxDipoleKernel.h
	--- a/DipoleShower/Kernels/FIgx2ggxDipoleKernel.h
	+++ b/DipoleShower/Kernels/FIgx2ggxDipoleKernel.h
	@@ -1,188 +1,190 @@
	// -- C++ --
	#ifndef HERWIG_FIgx2ggxDipoleKernel_H
	#define HERWIG_FIgx2ggxDipoleKernel_H
	//
	// This is the declaration of the FIgx2ggxDipoleKernel class.
	//

	#include "DipoleSplittingKernel.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* \ingroup DipoleShower
	* \author Simon Platzer
	*
	* \brief FIgx2ggxDipoleKernel implements the g -> gg
	* splitting off a final-initial dipole
	*
	*/
	class FIgx2ggxDipoleKernel: public DipoleSplittingKernel {

	public:

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

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

	public:

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

	/**
	* 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;

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

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

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

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

	- virtual double estimate(Energy up,Energy down) const {assert(false);return 1.;}
	+ virtual double estimate(Energy up,Energy down) const {
	+ return exp(-1.theSymmetryFactor3.(pow(log(up/down),2)- 11./12. log(up/down)) 0.2/(2.3.14));
	+ }

	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;
	//@}


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


	private:
	/**
	* Symmetry factor for final state gluon splittings (should be 1/2).
	*/

	double theSymmetryFactor;

	/**
	* The static object used to initialize the description of this class.
	* Indicates that this is a concrete class with persistent data.
	*/
	static ClassDescription<FIgx2ggxDipoleKernel> initFIgx2ggxDipoleKernel;

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

	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

	/** This template specialization informs ThePEG about the name of
	* the FIgx2ggxDipoleKernel class and the shared object where it is defined. */
	template <>
	struct ClassTraits<Herwig::FIgx2ggxDipoleKernel>
	: public ClassTraitsBase<Herwig::FIgx2ggxDipoleKernel> {
	/** Return a platform-independent class name */
	static string className() { return "Herwig::FIgx2ggxDipoleKernel"; }
	/**
	* The name of a file containing the dynamic library where the class
	* FIgx2ggxDipoleKernel is implemented. It may also include several, space-separated,
	* libraries if the class FIgx2ggxDipoleKernel 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_FIgx2ggxDipoleKernel_H */
	diff --git a/DipoleShower/Kernels/FIgx2qqxDipoleKernel.h b/DipoleShower/Kernels/FIgx2qqxDipoleKernel.h
	--- a/DipoleShower/Kernels/FIgx2qqxDipoleKernel.h
	+++ b/DipoleShower/Kernels/FIgx2qqxDipoleKernel.h
	@@ -1,183 +1,185 @@
	// -- C++ --
	#ifndef HERWIG_FIgx2qqxDipoleKernel_H
	#define HERWIG_FIgx2qqxDipoleKernel_H
	//
	// This is the declaration of the FIgx2qqxDipoleKernel class.
	//

	#include "DipoleSplittingKernel.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* \ingroup DipoleShower
	* \author Simon Platzer
	*
	* \brief FIgx2qqxDipoleKernel implements the g -> qqbar
	* splitting off a final-initial dipole
	*
	*/
	class FIgx2qqxDipoleKernel: public DipoleSplittingKernel {

	public:

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

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

	public:

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

	/**
	* 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;

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

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

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

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

	- virtual double estimate(Energy up,Energy down) const {assert(false);return 1.;}
	+ virtual double estimate(Energy up,Energy down) const {
	+ return exp(-1.1./4.log(up/down)0.2/(2.3.14));
	+ }

	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;
	//@}


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


	private:

	/**
	* The static object used to initialize the description of this class.
	* Indicates that this is a concrete class with persistent data.
	*/
	static ClassDescription<FIgx2qqxDipoleKernel> initFIgx2qqxDipoleKernel;

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

	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

	/** This template specialization informs ThePEG about the name of
	* the FIgx2qqxDipoleKernel class and the shared object where it is defined. */
	template <>
	struct ClassTraits<Herwig::FIgx2qqxDipoleKernel>
	: public ClassTraitsBase<Herwig::FIgx2qqxDipoleKernel> {
	/** Return a platform-independent class name */
	static string className() { return "Herwig::FIgx2qqxDipoleKernel"; }
	/**
	* The name of a file containing the dynamic library where the class
	* FIgx2qqxDipoleKernel is implemented. It may also include several, space-separated,
	* libraries if the class FIgx2qqxDipoleKernel 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_FIgx2qqxDipoleKernel_H */
	diff --git a/DipoleShower/Kernels/FIqx2qgxDipoleKernel.h b/DipoleShower/Kernels/FIqx2qgxDipoleKernel.h
	--- a/DipoleShower/Kernels/FIqx2qgxDipoleKernel.h
	+++ b/DipoleShower/Kernels/FIqx2qgxDipoleKernel.h
	@@ -1,183 +1,185 @@
	// -- C++ --
	#ifndef HERWIG_FIqx2qgxDipoleKernel_H
	#define HERWIG_FIqx2qgxDipoleKernel_H
	//
	// This is the declaration of the FIqx2qgxDipoleKernel class.
	//

	#include "DipoleSplittingKernel.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* \ingroup DipoleShower
	* \author Simon Platzer
	*
	* \brief FIqx2qgxDipoleKernel implements the q -> qg
	* splitting off a final-initial dipole
	*
	*/
	class FIqx2qgxDipoleKernel: public DipoleSplittingKernel {

	public:

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

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

	public:

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

	/**
	* 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;

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

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

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

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

	- virtual double estimate(Energy up,Energy down) const {assert(false);return 1.;}
	-
	+ virtual double estimate(Energy up,Energy down) const {
	+ return exp(-1.4./3.(pow(log(up/down),2)-3./4.log(up/down))0.2/(2.*3.14));
	+ }
	+
	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;
	//@}


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


	private:

	/**
	* The static object used to initialize the description of this class.
	* Indicates that this is a concrete class with persistent data.
	*/
	static ClassDescription<FIqx2qgxDipoleKernel> initFIqx2qgxDipoleKernel;

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

	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

	/** This template specialization informs ThePEG about the name of
	* the FIqx2qgxDipoleKernel class and the shared object where it is defined. */
	template <>
	struct ClassTraits<Herwig::FIqx2qgxDipoleKernel>
	: public ClassTraitsBase<Herwig::FIqx2qgxDipoleKernel> {
	/** Return a platform-independent class name */
	static string className() { return "Herwig::FIqx2qgxDipoleKernel"; }
	/**
	* The name of a file containing the dynamic library where the class
	* FIqx2qgxDipoleKernel is implemented. It may also include several, space-separated,
	* libraries if the class FIqx2qgxDipoleKernel 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_FIqx2qgxDipoleKernel_H */
	diff --git a/DipoleShower/Kernels/IFgx2ggxDipoleKernel.h b/DipoleShower/Kernels/IFgx2ggxDipoleKernel.h
	--- a/DipoleShower/Kernels/IFgx2ggxDipoleKernel.h
	+++ b/DipoleShower/Kernels/IFgx2ggxDipoleKernel.h
	@@ -1,183 +1,185 @@
	// -- C++ --
	#ifndef HERWIG_IFgx2ggxDipoleKernel_H
	#define HERWIG_IFgx2ggxDipoleKernel_H
	//
	// This is the declaration of the IFgx2ggxDipoleKernel class.
	//

	#include "DipoleSplittingKernel.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* \ingroup DipoleShower
	* \author Simon Platzer
	*
	* \brief IFgx2ggxDipoleKernel implements the g -> gg
	* splitting off an initial-final dipole
	*
	*/
	class IFgx2ggxDipoleKernel: public DipoleSplittingKernel {

	public:

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

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

	public:

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

	/**
	* 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;

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

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

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

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

	- virtual double estimate(Energy up,Energy down) const {assert(false);return 1.;}
	+ virtual double estimate(Energy up,Energy down) const {
	+ return exp(-1.3.(pow(log(up/down),2)- 11./12.* log(up/down)) 0.2/(2.3.14));
	+ }

	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;
	//@}


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


	private:

	/**
	* The static object used to initialize the description of this class.
	* Indicates that this is a concrete class with persistent data.
	*/
	static ClassDescription<IFgx2ggxDipoleKernel> initIFgx2ggxDipoleKernel;

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

	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

	/** This template specialization informs ThePEG about the name of
	* the IFgx2ggxDipoleKernel class and the shared object where it is defined. */
	template <>
	struct ClassTraits<Herwig::IFgx2ggxDipoleKernel>
	: public ClassTraitsBase<Herwig::IFgx2ggxDipoleKernel> {
	/** Return a platform-independent class name */
	static string className() { return "Herwig::IFgx2ggxDipoleKernel"; }
	/**
	* The name of a file containing the dynamic library where the class
	* IFgx2ggxDipoleKernel is implemented. It may also include several, space-separated,
	* libraries if the class IFgx2ggxDipoleKernel 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_IFgx2ggxDipoleKernel_H */
	diff --git a/DipoleShower/Kernels/IFgx2qqxDipoleKernel.h b/DipoleShower/Kernels/IFgx2qqxDipoleKernel.h
	--- a/DipoleShower/Kernels/IFgx2qqxDipoleKernel.h
	+++ b/DipoleShower/Kernels/IFgx2qqxDipoleKernel.h
	@@ -1,183 +1,185 @@
	// -- C++ --
	#ifndef HERWIG_IFgx2qqxDipoleKernel_H
	#define HERWIG_IFgx2qqxDipoleKernel_H
	//
	// This is the declaration of the IFgx2qqxDipoleKernel class.
	//

	#include "DipoleSplittingKernel.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* \ingroup DipoleShower
	* \author Simon Platzer
	*
	* \brief IFgx2qqxDipoleKernel implements the g -> qq
	* splitting off an initial-final dipole
	*
	*/
	class IFgx2qqxDipoleKernel: public DipoleSplittingKernel {

	public:

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

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

	public:

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

	/**
	* 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;

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

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

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

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

	- virtual double estimate(Energy up,Energy down) const {assert(false);return 1.;}
	-
	+ virtual double estimate(Energy up,Energy down) const {
	+ return exp(-1.0.54./3.(pow(log(up/down),2)-3./4.log(up/down))0.2/(2.3.14));
	+ }
	+
	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;
	//@}


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


	private:

	/**
	* The static object used to initialize the description of this class.
	* Indicates that this is a concrete class with persistent data.
	*/
	static ClassDescription<IFgx2qqxDipoleKernel> initIFgx2qqxDipoleKernel;

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

	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

	/** This template specialization informs ThePEG about the name of
	* the IFgx2qqxDipoleKernel class and the shared object where it is defined. */
	template <>
	struct ClassTraits<Herwig::IFgx2qqxDipoleKernel>
	: public ClassTraitsBase<Herwig::IFgx2qqxDipoleKernel> {
	/** Return a platform-independent class name */
	static string className() { return "Herwig::IFgx2qqxDipoleKernel"; }
	/**
	* The name of a file containing the dynamic library where the class
	* IFgx2qqxDipoleKernel is implemented. It may also include several, space-separated,
	* libraries if the class IFgx2qqxDipoleKernel 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_IFgx2qqxDipoleKernel_H */
	diff --git a/DipoleShower/Kernels/IFqx2gqxDipoleKernel.h b/DipoleShower/Kernels/IFqx2gqxDipoleKernel.h
	--- a/DipoleShower/Kernels/IFqx2gqxDipoleKernel.h
	+++ b/DipoleShower/Kernels/IFqx2gqxDipoleKernel.h
	@@ -1,183 +1,185 @@
	// -- C++ --
	#ifndef HERWIG_IFqx2gqxDipoleKernel_H
	#define HERWIG_IFqx2gqxDipoleKernel_H
	//
	// This is the declaration of the IFqx2gqxDipoleKernel class.
	//

	#include "DipoleSplittingKernel.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* \ingroup DipoleShower
	* \author Simon Platzer
	*
	* \brief IFqx2gqxDipoleKernel implements the q -> gqbar
	* splitting off an initial-final dipole
	*
	*/
	class IFqx2gqxDipoleKernel: public DipoleSplittingKernel {

	public:

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

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

	public:

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

	/**
	* 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;

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

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

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

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

	- virtual double estimate(Energy up,Energy down) const {assert(false);return 1.;}
	-
	+ virtual double estimate(Energy up,Energy down) const {
	+ return exp(-1.1./4.log(up/down)0.2/(2.3.14));
	+ }
	+
	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;
	//@}


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


	private:

	/**
	* The static object used to initialize the description of this class.
	* Indicates that this is a concrete class with persistent data.
	*/
	static ClassDescription<IFqx2gqxDipoleKernel> initIFqx2gqxDipoleKernel;

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

	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

	/** This template specialization informs ThePEG about the name of
	* the IFqx2gqxDipoleKernel class and the shared object where it is defined. */
	template <>
	struct ClassTraits<Herwig::IFqx2gqxDipoleKernel>
	: public ClassTraitsBase<Herwig::IFqx2gqxDipoleKernel> {
	/** Return a platform-independent class name */
	static string className() { return "Herwig::IFqx2gqxDipoleKernel"; }
	/**
	* The name of a file containing the dynamic library where the class
	* IFqx2gqxDipoleKernel is implemented. It may also include several, space-separated,
	* libraries if the class IFqx2gqxDipoleKernel 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_IFqx2gqxDipoleKernel_H */
	diff --git a/DipoleShower/Kernels/IFqx2qgxDipoleKernel.h b/DipoleShower/Kernels/IFqx2qgxDipoleKernel.h
	--- a/DipoleShower/Kernels/IFqx2qgxDipoleKernel.h
	+++ b/DipoleShower/Kernels/IFqx2qgxDipoleKernel.h
	@@ -1,183 +1,185 @@
	// -- C++ --
	#ifndef HERWIG_IFqx2qgxDipoleKernel_H
	#define HERWIG_IFqx2qgxDipoleKernel_H
	//
	// This is the declaration of the IFqx2qgxDipoleKernel class.
	//

	#include "DipoleSplittingKernel.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* \ingroup DipoleShower
	* \author Simon Platzer
	*
	* \brief IFqx2qgxDipoleKernel implements the q -> qg
	* splitting off an initial-final dipole
	*
	*/
	class IFqx2qgxDipoleKernel: public DipoleSplittingKernel {

	public:

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

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

	public:

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

	/**
	* 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;

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

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

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

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

	- virtual double estimate(Energy up,Energy down) const {assert(false);return 1.;}
	-
	+ virtual double estimate(Energy up,Energy down) const {
	+ return exp(-1.0.54./3.(pow(log(up/down),2)-3./4.log(up/down))0.2/(2.3.14));
	+ }
	+
	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;
	//@}


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


	private:

	/**
	* The static object used to initialize the description of this class.
	* Indicates that this is a concrete class with persistent data.
	*/
	static ClassDescription<IFqx2qgxDipoleKernel> initIFqx2qgxDipoleKernel;

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

	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

	/** This template specialization informs ThePEG about the name of
	* the IFqx2qgxDipoleKernel class and the shared object where it is defined. */
	template <>
	struct ClassTraits<Herwig::IFqx2qgxDipoleKernel>
	: public ClassTraitsBase<Herwig::IFqx2qgxDipoleKernel> {
	/** Return a platform-independent class name */
	static string className() { return "Herwig::IFqx2qgxDipoleKernel"; }
	/**
	* The name of a file containing the dynamic library where the class
	* IFqx2qgxDipoleKernel is implemented. It may also include several, space-separated,
	* libraries if the class IFqx2qgxDipoleKernel 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_IFqx2qgxDipoleKernel_H */
	diff --git a/DipoleShower/Kernels/IIgx2ggxDipoleKernel.h b/DipoleShower/Kernels/IIgx2ggxDipoleKernel.h
	--- a/DipoleShower/Kernels/IIgx2ggxDipoleKernel.h
	+++ b/DipoleShower/Kernels/IIgx2ggxDipoleKernel.h
	@@ -1,185 +1,185 @@
	// -- C++ --
	#ifndef HERWIG_IIgx2ggxDipoleKernel_H
	#define HERWIG_IIgx2ggxDipoleKernel_H
	//
	// This is the declaration of the IIgx2ggxDipoleKernel class.
	//

	#include "DipoleSplittingKernel.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* \ingroup DipoleShower
	* \author Simon Platzer
	*
	* \brief IIgx2ggxDipoleKernel implements the g -> gg
	* splitting off an initial-initial dipole
	*
	*/
	class IIgx2ggxDipoleKernel: public DipoleSplittingKernel {

	public:

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

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

	public:

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

	/**
	* 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;

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

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

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

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

	virtual double estimate(Energy up,Energy down) const {
	- return exp(-1.3.(pow(log(up/down),2)- 11./12.* log(up/down)) 0.2/(2.3.14));
	+ return exp(-1.0.53.(pow(log(up/down),2)- 11./12. log(up/down)) 0.2/(2.3.14));
	}

	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;
	//@}


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


	private:

	/**
	* The static object used to initialize the description of this class.
	* Indicates that this is a concrete class with persistent data.
	*/
	static ClassDescription<IIgx2ggxDipoleKernel> initIIgx2ggxDipoleKernel;

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

	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

	/** This template specialization informs ThePEG about the name of
	* the IIgx2ggxDipoleKernel class and the shared object where it is defined. */
	template <>
	struct ClassTraits<Herwig::IIgx2ggxDipoleKernel>
	: public ClassTraitsBase<Herwig::IIgx2ggxDipoleKernel> {
	/** Return a platform-independent class name */
	static string className() { return "Herwig::IIgx2ggxDipoleKernel"; }
	/**
	* The name of a file containing the dynamic library where the class
	* IIgx2ggxDipoleKernel is implemented. It may also include several, space-separated,
	* libraries if the class IIgx2ggxDipoleKernel 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_IIgx2ggxDipoleKernel_H */
	diff --git a/Sampling/MonacoSampler.cc b/Sampling/MonacoSampler.cc
	--- a/Sampling/MonacoSampler.cc
	+++ b/Sampling/MonacoSampler.cc
	@@ -1,440 +1,440 @@
	// -- C++ --
	//
	// MonacoSampler.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2012 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 MonacoSampler class.
	//

	#include "MonacoSampler.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/EventRecord/Particle.h"
	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "ThePEG/Repository/Repository.h"

	#include "ThePEG/Interface/Parameter.h"

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

	#include "ThePEG/Handlers/StandardEventHandler.h"
	#include "ThePEG/Handlers/StandardXComb.h"

	#include <boost/progress.hpp>

	#include "MonacoSampler.h"
	#include "Herwig/Sampling/GeneralSampler.h"
	#include <ctime>

	using namespace Herwig;

	MonacoSampler::MonacoSampler()
	: BinSampler(),
	theAlpha(0.875),
	theGridDivisions(48),
	theIterationPoints(0) {}

	MonacoSampler::~MonacoSampler() {}

	IBPtr MonacoSampler::clone() const {
	return new_ptr(*this);
	}

	IBPtr MonacoSampler::fullclone() const {
	return new_ptr(*this);
	}

	double MonacoSampler::generate() {

	double w = 1.;
	double wref=1.;
	// cout<<"\npoint: ";
	std::valarray<int> upperb(dimension());
	for ( int k = 0; k < dimension(); ++k ) {
	double div = (1 - UseRandom::rnd()) * theGridDivisions;
	upperb[k] = static_cast<int>(div);
	double gupper, glower;
	if ( upperb[k] <= 0 ) {
	upperb[k] = 0;
	glower = 0.;
	gupper = theGrid(k,0);
	} else if (upperb[k] >= static_cast<int>(theGridDivisions)) {
	upperb[k] = theGridDivisions-1;
	glower = theGrid(k,theGridDivisions-2);
	gupper = theGrid(k,theGridDivisions-1);
	} else {
	glower = theGrid(k,upperb[k]-1);
	gupper = theGrid(k,upperb[k]);
	}
	double gdiff = gupper - glower;
	lastPoint()[k] = glower + (div-upperb[k])*gdiff;
	w = gdiff theGridDivisions;
	}
	double elapsed_secs;

	bool canImprove = sampler()->almostUnweighted();
	eventHandler()->resetDidEstimate();
	try {
	/**
	* For AlmostUnweighted Events we can speed up the calculation.
	* 1. Calculate a fast estimate (e.g. Born, analytic CKKW ...)
	* 2. If accepted reweight with full/estimate.
	* (Born+Virt, CKKW with sampled sudakovs)
	**/
	clock_t begin = clock();


	wref=eventHandler()->dSigDR(lastPoint(),canImprove) / nanobarn;
	clock_t end = clock();
	elapsed_secs = double(end - begin) / CLOCKS_PER_SEC;

	w *= wref;
	} catch (Veto&) {
	w = 0.0;
	} catch (...) {
	throw;
	}
	// only store numbers
	double wgt = w;
	if ( isnan(wgt) \|\| isinf(wgt) ) wgt = 0;
	// save results for later grid optimization
	theIterationPoints++;
	for ( int k = 0; k < dimension(); ++k ) {
	theGridData(k,upperb[k]) += wgt*wgt;
	}

	if (randomNumberString()!="")
	for ( size_t k = 0; k < lastPoint().size(); ++k ) {
	RandomNumberHistograms[RandomNumberIndex(id(),k)].first.book(lastPoint()[k],wgt);
	RandomNumberHistograms[RandomNumberIndex(id(),k)].second+=wgt;
	}

	if ( !weighted() && initialized() ) {
	double p = min(abs(w),referenceWeight())/referenceWeight();
	double sign = w >= 0. ? 1. : -1.;
	if ( p < 1 && UseRandom::rnd() > p )
	w = 0.;
	else
	w = sign*max(abs(w),referenceWeight());
	}



	//if (eventHandler()->didEstimate())
	//cout<<"\n"<<canImprove <<" "<<eventHandler()->didEstimate()<<" "<<weighted()<<" "<<initialized()<<" ---------------------------->"<<w;

	if(canImprove && eventHandler()->didEstimate()
	- && (w!=0.) && ((!weighted()&&initialized())\|\|UseRandom::rnd()<0.01)){
	+ && (w!=0.) && ((!weighted()&&initialized())\|\|UseRandom::rnd()<0.1)){


	double x=elapsed_secs;

	clock_t begin = clock();


	double wfull=eventHandler()->dSigDR(lastPoint(),false) / nanobarn;
	clock_t end = clock();
	elapsed_secs = double(end - begin) / CLOCKS_PER_SEC;
	cout<<"\nwfull/wref: "<<wfull/wref <<" time ratio "<<x/elapsed_secs;


	w*=abs(wfull/wref);

	}
	select(w);
	if ( w != 0.0 )
	accept();
	return w;
	}

	void MonacoSampler::saveGrid() const {
	XML::Element grid = toXML();
	grid.appendAttribute("process",id());
	sampler()->grids().append(grid);
	}

	void MonacoSampler::initialize(bool progress) {

	//read in grid
	bool haveGrid = false;
	list<XML::Element>::iterator git = sampler()->grids().children().begin();
	for ( ; git != sampler()->grids().children().end(); ++git ) {
	if ( git->type() != XML::ElementTypes::Element )
	continue;
	if ( git->name() != "Monaco" )
	continue;
	string proc;
	git->getFromAttribute("process",proc);
	if ( proc == id() ) {
	haveGrid = true;
	break;
	}
	}

	if ( haveGrid ) {
	fromXML(*git);
	sampler()->grids().erase(git);
	didReadGrids();
	} else {
	// flat grid
	theGrid.resize(dimension(),theGridDivisions);
	for (int k = 0; k < dimension(); k++)
	for (size_t l = 0; l < theGridDivisions; l++)
	theGrid(k,l) = (l+1)/static_cast<double>(theGridDivisions);
	theGridData = boost::numeric::ublas::zero_matrix<double>(dimension(),theGridDivisions);
	theIterationPoints = 0;
	}

	lastPoint().resize(dimension());
	if (randomNumberString()!="")
	for(size_t i=0;i<lastPoint().size();i++){
	RandomNumberHistograms[RandomNumberIndex(id(),i)] = make_pair( RandomNumberHistogram(),0.);
	}

	if ( initialized() ) {
	if ( !hasGrids() )
	throw Exception() << "MonacoSampler: Require existing grid when starting to run.\n"
	<< "Did you miss setting --setupfile?"
	<< Exception::abortnow;
	return;
	}

	if ( haveGrid ) {
	if ( !integrated() ) {
	runIteration(initialPoints(),progress);
	adapt();
	}
	isInitialized();
	return;
	}

	// if ( !sampler()->grids().children().empty() ) {
	// nIterations(1);
	// }
	unsigned long points = initialPoints();
	for ( unsigned long k = 0; k < nIterations(); ++k ) {
	runIteration(points,progress);
	if ( k < nIterations() - 1 ) {
	points = (unsigned long)(points*enhancementFactor());
	adapt();
	nextIteration();
	}
	}
	adapt();
	didReadGrids();
	isInitialized();
	}

	void MonacoSampler::adapt() {

	int dim = dimension();

	// refine grid
	std::valarray<double> gridcumul(dim);
	for (int k=0; k<dim; ++k) {
	double gridold = theGridData(k,0);
	double gridnew = theGridData(k,1);
	theGridData(k,0) = (gridold + gridnew) / 2.0;
	gridcumul[k] = theGridData(k,0);
	for (size_t l=1; l<theGridDivisions-1; ++l) {
	theGridData(k,l) = gridold + gridnew;
	gridold = gridnew;
	gridnew = theGridData(k,l+1);
	theGridData(k,l) = (theGridData(k,l) + gridnew) / 3.0;
	gridcumul[k] += theGridData(k,l);
	}
	theGridData(k,theGridDivisions-1) = (gridnew + gridold) / 2.0;
	gridcumul[k] += theGridData(k,theGridDivisions-1);
	}

	for (int k=0; k<dim; ++k) {
	double rc = 0.;
	std::valarray<double> ri(theGridDivisions);
	for (size_t l=0; l<theGridDivisions; ++l) {
	ri[l] = 0.;
	if ((theGridData(k,l) >= 0) && (gridcumul[k] != 0)) {
	theGridData(k,l) = max( 1.0e-30, theGridData(k,l) );
	double gpart = gridcumul[k] / theGridData(k,l);
	ri[l] = pow( (gpart - 1.0) / (gpart * log( gpart )), theAlpha);
	} else {
	ri[l] = pow( 1. / log( 1e30 ), theAlpha);
	}
	rc += ri[l];
	}
	rc /= theGridDivisions;

	double gridold = 0, gridnew = 0.;
	double deltar = 0.;
	unsigned int m = 0;
	std::valarray<double> theGridRowNew(theGridDivisions);

	for (size_t l = 0; l < theGridDivisions; ++l) {
	deltar += ri[l];
	gridold = gridnew;
	gridnew = theGrid(k,l);
	for (; deltar > rc; m++) {
	deltar -= rc;
	theGridRowNew[m] = gridnew - (gridnew - gridold) * deltar / ri[l];
	}
	}
	for (size_t l = 0; l < theGridDivisions-1; ++l) {
	theGrid(k,l) = theGridRowNew[l];
	}
	theGrid(k,theGridDivisions-1) = 1.0;
	}

	theGridData = boost::numeric::ublas::zero_matrix<double>(dimension(),theGridDivisions);
	theIterationPoints = 0;
	}

	void MonacoSampler::finalize(bool) {
	// save grid
	adapt();
	XML::Element grid = MonacoSampler::toXML();
	grid.appendAttribute("process",id());
	sampler()->grids().append(grid);

	if (randomNumberString()!="")
	for ( map<RandomNumberIndex,pair<RandomNumberHistogram,double> >::
	const_iterator b = RandomNumberHistograms.begin();
	b != RandomNumberHistograms.end(); ++b ) {
	b->second.first.dump(randomNumberString(), b->first.first,shortprocess(),b->first.second);
	}

	}

	void MonacoSampler::fromXML(const XML::Element& grid) {
	int dim = 0;
	grid.getFromAttribute("Dimension",dim);
	if ( dim != dimension() ) {
	throw std::runtime_error("[MonacoSampler] Number of dimensions in grid file does not match expectation.");
	}

	size_t griddivisions = 0;
	grid.getFromAttribute("GridDivisions",griddivisions);
	boost::numeric::ublas::matrix<double> tmpgrid(dim,griddivisions);

	pair<multimap<pair<int,string>,list<XML::Element>::iterator>::const_iterator,multimap<pair<int,string>,list<XML::Element>::iterator>::const_iterator> cit;
	cit = grid.findAll(XML::ElementTypes::Element,"GridVector");

	if ( cit.first->second == grid.children().end() )
	throw std::runtime_error("[MonacoSampler] Expected a GridVector element.");

	for (multimap<pair<int,string>,list<XML::Element>::iterator>::const_iterator iit=cit.first; iit!=cit.second; ++iit) {
	const XML::Element& gridvector = *iit->second;


	int k = 0;
	gridvector.getFromAttribute("Index",k);
	if ( k >= dim ) {
	throw std::runtime_error("[MonacoSampler] Index of grid dimension larger than grid size.");
	} else {
	list<XML::Element>::const_iterator git;
	git = gridvector.findFirst(XML::ElementTypes::ParsedCharacterData,"");
	if ( git == gridvector.children().end() )
	throw std::runtime_error("[MonacoSampler] Expected grid data.");
	istringstream bdata(git->content());

	for ( size_t l = 0; l < griddivisions; ++l ) {
	bdata >> tmpgrid(k,l);
	}
	}

	}

	// store back into main variable
	// if griddivisions do not match, rebin preserving bin density

	theGrid.resize(dim,theGridDivisions);
	theIterationPoints = 0;
	double divratio = griddivisions / static_cast<double>(theGridDivisions);
	for (int k = 0; k < dim; k++) {
	double xold = 0, xnew = 0, deltar = 0;
	size_t l = 0;
	for (size_t m = 0; m < griddivisions; m++) {
	deltar += 1;
	xold = xnew;
	xnew = tmpgrid(k,m);
	for (; deltar > divratio; l++) {
	deltar -= divratio;
	theGrid(k,l) = xnew - (xnew - xold) * deltar;
	}
	}
	theGrid(k,theGridDivisions-1) = 1.0;
	}

	theGridData = boost::numeric::ublas::zero_matrix<double>(dimension(),theGridDivisions);

	}

	XML::Element MonacoSampler::toXML() const {

	XML::Element grid(XML::ElementTypes::Element,"Monaco");
	grid.appendAttribute("Dimension",dimension());
	grid.appendAttribute("GridDivisions",theGridDivisions);

	for ( int k = 0; k < dimension(); ++k ) {
	XML::Element gridvector(XML::ElementTypes::Element,"GridVector");
	gridvector.appendAttribute("Index",k);
	ostringstream bdata;
	bdata << setprecision(20);
	for ( size_t l = 0; l < theGridDivisions; ++l )
	bdata << theGrid(k,l) << " ";

	XML::Element belem(XML::ElementTypes::ParsedCharacterData,bdata.str());
	gridvector.append(belem);
	grid.append(gridvector);
	}

	return grid;

	}

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


	void MonacoSampler::persistentOutput(PersistentOStream & os) const {
	BinSampler::put(os);
	os << theAlpha << theGridDivisions;
	}

	void MonacoSampler::persistentInput(PersistentIStream & is, int) {
	BinSampler::get(is);
	is >> theAlpha >> theGridDivisions;
	}


	// * Attention * The following static variable is needed for the type
	// description system in ThePEG. Please check that the template arguments
	// are correct (the class and its base class), and that the constructor
	// arguments are correct (the class name and the name of the dynamically
	// loadable library where the class implementation can be found).
	DescribeClass<MonacoSampler,BinSampler>
	describeHerwigMonacoSampler("Herwig::MonacoSampler", "HwSampling.so");

	void MonacoSampler::Init() {

	static ClassDocumentation<MonacoSampler> documentation
	("MonacoSampler samples XCombs bins. This implementation performs weighted MC integration using Monaco, an adapted Vegas algorithm.");

	static Parameter<MonacoSampler,double> interfaceAlpha
	("Alpha",
	"Rate of grid modification (0 for no modification).",
	&MonacoSampler::theAlpha, 0.875, 0.0, 0,
	false, false, Interface::lowerlim);

	static Parameter<MonacoSampler,size_t> interfaceGridDivisions
	("GridDivisions",
	"The number of divisions per grid dimension.",
	&MonacoSampler::theGridDivisions, 48, 1, 0,
	false, false, Interface::lowerlim);

	}

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 3:32 PM (1 d, 18 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3801233
Default Alt Text: (76 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions