No OneTemporary
Actions

Size

118 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,587 +1,608 @@
	// -- 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::resetVariations() {
	+ for ( map<string,double>::iterator w = currentWeights.begin();
	+ w != currentWeights.end(); ++w )
	+ w->second = 1.;
	+}
	+
	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) {

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

	DipoleSplittingInfo& split =
	( !presampling ? generatedSplitting : presampledSplitting );

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

	// multiply in the profile scales when relevant
	assert(ShowerHandler::currentHandler());
	if ( ShowerHandler::currentHandler()->firstInteraction() &&
	ShowerHandler::currentHandler()->profileScales() &&
	!presampling ) {
	Energy hard = ShowerHandler::currentHandler()->hardScale();
	if ( hard > ZERO )
	kernel *= ShowerHandler::currentHandler()->profileScales()->hardScaleProfile(hard,split.lastPt());
	}

	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) {
	+void DipoleSplittingGenerator::doGenerate(map<string,double>& variations,
	+ 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());
	}

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

	+ for ( map<string,double>::const_iterator w = currentWeights.begin();
	+ w != currentWeights.end(); ++w ) {
	+ map<string,double>::iterator v = variations.find(w->first);
	+ if ( v != variations.end() )
	+ v->second *= w->second;
	+ else
	+ variations[w->first] = w->second;
	+ }
	+
	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,
	+ map<string,double>& variations,
	Energy optHardPt,
	Energy optCutoff) {

	fixParameters(split,optHardPt);

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

	- doGenerate(optCutoff);
	+ doGenerate(variations,optCutoff);

	return generatedSplitting.lastPt();

	}

	Energy DipoleSplittingGenerator::generateWrapped(DipoleSplittingInfo& split,
	+ map<string,double>& variations,
	Energy optHardPt,
	Energy optCutoff) {

	assert(!wrapping());

	DipoleSplittingInfo backup = generatedSplitting;
	generatedSplitting = split;

	fixParameters(split,optHardPt);

	try {
	- doGenerate(optCutoff);
	+ doGenerate(variations,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/Base/DipoleSplittingGenerator.h b/DipoleShower/Base/DipoleSplittingGenerator.h
	--- a/DipoleShower/Base/DipoleSplittingGenerator.h
	+++ b/DipoleShower/Base/DipoleSplittingGenerator.h
	@@ -1,449 +1,476 @@
	// -- C++ --
	//
	// DipoleSplittingGenerator.h 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.
	//
	#ifndef HERWIG_DipoleSplittingGenerator_H
	#define HERWIG_DipoleSplittingGenerator_H
	//
	// This is the declaration of the DipoleSplittingGenerator class.
	//

	#include "ThePEG/Handlers/HandlerBase.h"

	#include "Herwig/DipoleShower/Kernels/DipoleSplittingKernel.h"
	#include "DipoleSplittingReweight.h"
	#include "Herwig/DipoleShower/Utility/DipoleMCCheck.h"
	#include "Herwig/Sampling/exsample/exponential_generator.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* \ingroup DipoleShower
	* \author Simon Platzer
	*
	* \brief DipoleSplittingGenerator is used by the dipole shower
	* to sample splittings from a given dipole splitting kernel.
	*
	* @see \ref DipoleSplittingGeneratorInterfaces "The interfaces"
	* defined for DipoleSplittingGenerator.
	*/
	class DipoleSplittingGenerator: public HandlerBase {

	public:

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

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

	public:

	/**
	* Return the dipole splitting kernel.
	*/
	Ptr<DipoleSplittingKernel>::tptr splittingKernel() const;

	/**
	* Return the dipole splitting reweight.
	*/
	Ptr<DipoleSplittingReweight>::tptr splittingReweight() const;

	/**
	* Return the dipole splitting kinematics.
	*/
	Ptr<DipoleSplittingKinematics>::tptr splittingKinematics() const;

	/**
	* Set the dipole splitting kernel.
	*/
	void splittingKernel(Ptr<DipoleSplittingKernel>::tptr sp);

	/**
	* Set the dipole splitting reweight.
	*/
	void splittingReweight(Ptr<DipoleSplittingReweight>::tptr sp);

	/**
	* Make a wrapper around another generator.
	*/
	void wrap(Ptr<DipoleSplittingGenerator>::ptr other);

	/**
	* Return true, if this is actually a wrapper around
	* another splitting generator.
	*/
	bool wrapping() const { return theOtherGenerator; }

	public:

	/**
	+ * Reset the current variations to one
	+ */
	+ void resetVariations();
	+
	+ /**
	* Prepare to fill the given splitting.
	*/
	void prepare(const DipoleSplittingInfo&);

	/**
	* Fix parameters from the fiven DipoleSplittingInfo
	* and generate the next splitting. Return the
	* pt selected for the next splitting.
	*/
	Energy generate(const DipoleSplittingInfo&,
	+ map<string,double>& variations,
	Energy optHardPt = ZERO,
	Energy optCutoff = ZERO);

	/**
	* Fix parameters from the fiven DipoleSplittingInfo
	* and generate the next splitting. Return the
	* pt selected for the next splitting when called
	* from a wrapping generator.
	*/
	Energy generateWrapped(DipoleSplittingInfo&,
	+ map<string,double>& variations,
	Energy optHardPt = ZERO,
	Energy optCutoff = ZERO);

	/**
	* Complete the given splitting.
	*/
	void completeSplitting(DipoleSplittingInfo&) const;

	/**
	* Return the last generated splitting
	*/
	const DipoleSplittingInfo& lastSplitting() const { return generatedSplitting; }

	public:

	/**
	* Print debug information on the splitting
	* handled.
	*/
	void debugGenerator(ostream&) const;

	/**
	* Print debug information on the last
	* generated event.
	*/
	void debugLastEvent(ostream&) const;

	protected:

	/**
	* Update parameters given a splitting.
	*/
	void fixParameters(const DipoleSplittingInfo&,
	Energy optHardPt = ZERO);

	/**
	* With the parameters previuosly supplied
	* through fixParameters generate the next
	* splitting.
	*/
	- void doGenerate(Energy optCutoff = ZERO);
	+ void doGenerate(map<string,double>& variations,
	+ Energy optCutoff = ZERO);

	public:

	/**
	* Return the number of random numbers
	* needed to sample this kernel.
	*/
	int nDim() const;

	/**
	* Flag, which variables are free variables.
	*/
	const vector<bool>& sampleFlags();

	/**
	* Return the support of the splitting kernel.
	* The lower bound on the first variable is
	* assumed to correspond to the cutoff on the
	* evolution variable.
	*/
	const pair<vector<double>,vector<double> >& support();

	/**
	* Return the parameter point associated to the splitting
	* previously supplied through fixParameters.
	*/
	const vector<double>& parameterPoint() const { return parameters; }

	/**
	* Indicate that presampling of this kernel
	* will be performed in the next calls to
	* evaluate until stopPresampling() is called.
	*/
	void startPresampling();

	/**
	* Indicate that presampling of this kernel
	* is done until startPresampling() is called.
	*/
	void stopPresampling();

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

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

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

	/**
	* Return true, if this splitting generator
	* is able to deliver an overestimate to the sampled
	* kernel.
	*/
	bool haveOverestimate() const;

	/**
	* Return an overestimate to the sampled kernel.
	*/
	bool overestimate(const vector<double>&);

	/**
	* Invert the integral over the overestimate to equal
	* the given value.
	*/
	double invertOverestimateIntegral(double) const;

	/**
	* Evalute the splitting kernel.
	*/
	double evaluate(const vector<double>&);

	/**
	+ * Indicate that a veto with the given kernel value and overestimate has occured.
	+ */
	+ void veto(const vector<double>&, double p, double r) {
	+ splittingKernel()->veto(generatedSplitting, p, r, currentWeights);
	+ }
	+
	+ /**
	+ * Indicate that an accept with the given kernel value and overestimate has occured.
	+ */
	+ void accept(const vector<double>&, double p, double r) {
	+ splittingKernel()->accept(generatedSplitting, p, r, currentWeights);
	+ }
	+
	+ /**
	* True, if sampler should apply compensation
	*/
	void doCompensate(bool yes = true) { theDoCompensate = yes; }

	public:

	/*@name Wrap to the exsample2 interface until this is finally cleaned up. /
	//@{

	inline const vector<bool>& variable_flags () {
	return sampleFlags();
	}

	inline size_t evolution_variable () const { return 0; }

	inline double evolution_cutoff () { return support().first[0]; }

	inline const vector<double>& parameter_point () const {
	return parameterPoint();
	}

	inline void start_presampling () {
	startPresampling();
	}

	inline void stop_presampling () {
	stopPresampling();
	}

	inline size_t dimension () const {
	return nDim();
	}

	inline unsigned long presampling_points () const {
	return presamplingPoints();
	}

	//@}

	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:

	/**
	* Pointer to another generator to wrap around.
	*/
	Ptr<DipoleSplittingGenerator>::ptr theOtherGenerator;

	/**
	* The dipole splitting kernel to sample
	* splitting from.
	*/
	Ptr<DipoleSplittingKernel>::ptr theSplittingKernel;

	/**
	* The dipole splitting reweight.
	*/
	Ptr<DipoleSplittingReweight>::ptr theSplittingReweight;

	/**
	* Pointer to the exponential generator
	*/
	exsample::exponential_generator<DipoleSplittingGenerator,UseRandom>*
	theExponentialGenerator;

	/**
	* The dipole splitting to be completed.
	*/
	DipoleSplittingInfo generatedSplitting;

	/**
	* A backup of the dipole splitting to be
	* completed, if this generator is presampled.
	*/
	DipoleSplittingInfo presampledSplitting;

	/**
	* True, if prepared to sample splittings
	* of a given kind.
	*/
	bool prepared;

	/**
	* Wether or not the kernel is currently
	* being presampled.
	*/
	bool presampling;

	/**
	* The parameter point.
	*/
	vector<double> parameters;

	/**
	* The sampling flags
	*/
	vector<bool> theFlags;

	/**
	* The support.
	*/
	pair<vector<double>,vector<double> > theSupport;

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

	/**
	* True, if sampler should apply compensation
	*/
	bool theDoCompensate;

	+ /**
	+ * The currently used weight map
	+ */
	+ map<string,double> currentWeights;
	+
	private:

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

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

	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

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

	#include <config.h>
	#include "DipoleShowerHandler.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Interface/RefVector.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Switch.h"

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

	// include theses to have complete types
	#include "Herwig/Shower/Base/Evolver.h"
	#include "Herwig/Shower/Base/ShowerParticle.h"
	#include "Herwig/PDF/MPIPDF.h"
	#include "Herwig/PDF/MinBiasPDF.h"
	#include "Herwig/Shower/Base/ShowerTree.h"
	#include "Herwig/Shower/Base/KinematicsReconstructor.h"
	#include "Herwig/Shower/Base/PartnerFinder.h"
	#include "Herwig/PDF/HwRemDecayer.h"

	#include "Herwig/DipoleShower/Utility/DipolePartonSplitter.h"

	#include "Herwig/MatrixElement/Matchbox/Base/SubtractedME.h"
	#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h"

	using namespace Herwig;

	bool DipoleShowerHandler::firstWarn = true;

	DipoleShowerHandler::DipoleShowerHandler()
	: ShowerHandler(), chainOrderVetoScales(true),
	nEmissions(0), discardNoEmissions(false), firstMCatNLOEmission(false),
	doFSR(true), doISR(true), realignmentScheme(0),
	verbosity(0), printEvent(0), nTries(0),
	didRadiate(false), didRealign(false),
	theRenormalizationScaleFreeze(1.*GeV),
	theFactorizationScaleFreeze(2.*GeV),
	isMCatNLOSEvent(false),
	isMCatNLOHEvent(false), theDoCompensate(false),
	theFreezeGrid(500000), maxPt(ZERO),
	muPt(ZERO) {}

	DipoleShowerHandler::~DipoleShowerHandler() {}

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

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

	tPPair DipoleShowerHandler::cascade(tSubProPtr sub, XCPtr,
	Energy optHardPt, Energy optCutoff) {

	useMe();

	prepareCascade(sub);

	+ currentWeights.clear();
	+
	if ( !doFSR && ! doISR )
	return sub->incoming();

	eventRecord().clear();
	eventRecord().prepare(sub,dynamic_ptr_cast<tStdXCombPtr>(lastXCombPtr()),pdfs());

	if ( eventRecord().outgoing().empty() && !doISR )
	return sub->incoming();
	if ( !eventRecord().incoming().first->coloured() &&
	!eventRecord().incoming().second->coloured() &&
	!doFSR )
	return sub->incoming();

	nTries = 0;

	while ( true ) {

	try {

	didRadiate = false;
	didRealign = false;

	isMCatNLOSEvent = false;
	isMCatNLOHEvent = false;

	if ( eventRecord().xcombPtr() ) {

	Ptr<SubtractedME>::tptr subme =
	dynamic_ptr_cast<Ptr<SubtractedME>::tptr>(eventRecord().xcombPtr()->matrixElement());
	Ptr<MatchboxMEBase>::tptr me =
	dynamic_ptr_cast<Ptr<MatchboxMEBase>::tptr>(eventRecord().xcombPtr()->matrixElement());
	Ptr<SubtractionDipole>::tptr dipme =
	dynamic_ptr_cast<Ptr<SubtractionDipole>::tptr>(eventRecord().xcombPtr()->matrixElement());


	if ( subme ) {
	if ( subme->showerApproximation() ) {
	// don't do this for POWHEG-type corrections
	if ( !subme->showerApproximation()->needsSplittingGenerator() ) {
	theShowerApproximation = subme->showerApproximation();
	if ( subme->realShowerSubtraction() )
	isMCatNLOHEvent = true;
	else if ( subme->virtualShowerSubtraction() )
	isMCatNLOSEvent = true;
	}
	}
	} else if ( me ) {
	if ( me->factory()->showerApproximation() ) {
	if ( !me->factory()->showerApproximation()->needsSplittingGenerator() ) {
	theShowerApproximation = me->factory()->showerApproximation();
	isMCatNLOSEvent = true;
	}
	}
	}

	string error = "Inconsistent hard emission set-up in DipoleShowerHandler::cascade. ";
	if (evolver()->hardEmissionMode()==1 \|\| evolver()->hardEmissionMode()==3 )
	throw Exception() << error
	<< "Cannot generate POWHEG corrections "
	<< "for particle decays using DipoleShowerHandler. "
	<< "Check value of Evolver:HardEmissionMode."
	<< Exception::runerror;
	if ( ( isMCatNLOSEvent \|\| isMCatNLOHEvent ) && evolver()->hardEmissionMode()==2)
	throw Exception() << error
	<< "Cannot generate POWHEG matching with MC@NLO shower approximation. "
	<< "Add 'set Evolver:HardEmissionMode 0' to input file."
	<< Exception::runerror;
	if (me && me->factory()->showerApproximation()){
	if(me->factory()->showerApproximation()->needsTruncatedShower())
	throw Exception() << error
	<< "No truncated shower needed with DipoleShowerHandler. Add "
	<< "'set MEMatching:TruncatedShower No' to input file."
	<< Exception::runerror;
	if (!( isMCatNLOSEvent \|\| isMCatNLOHEvent ) &&
	evolver()->hardEmissionMode()==0 && firstWarn){
	firstWarn=false;
	throw Exception() << error
	<< "Evolver:HardEmissionMode will be set to 'MatchboxPOWHEG'"
	<< Exception::warning;
	}
	}
	else if (subme && subme->factory()->showerApproximation()){
	if(subme->factory()->showerApproximation()->needsTruncatedShower())
	throw Exception() << error
	<< "No truncated shower needed with DipoleShowerHandler. Add "
	<< "'set MEMatching:TruncatedShower No' to input file."
	<< Exception::runerror;
	if (!( isMCatNLOSEvent \|\| isMCatNLOHEvent ) &&
	evolver()->hardEmissionMode()==0 && firstWarn){
	firstWarn=false;
	throw Exception() << error
	<< "Evolver:HardEmissionMode will be set to 'MatchboxPOWHEG'"
	<< Exception::warning;
	}
	}
	else if (dipme && evolver()->hardEmissionMode() == 0 && firstWarn){
	firstWarn=false;
	throw Exception() << error
	<< "Evolver:HardEmissionMode will be set to 'MatchboxPOWHEG'"
	<< Exception::warning;
	}
	else if (!dipme && evolver()->hardEmissionMode()==2 &&
	ShowerHandler::currentHandler()->firstInteraction())
	throw Exception() << error
	<< "POWHEG matching requested for LO events. Include "
	<< "'set Factory:ShowerApproximation MEMatching' in input file."
	<< Exception::runerror;

	}

	hardScales(lastXCombPtr()->lastCentralScale());

	if ( verbosity > 1 ) {
	generator()->log() << "DipoleShowerHandler starting off:\n";
	eventRecord().debugLastEvent(generator()->log());
	generator()->log() << flush;
	}

	unsigned int nEmitted = 0;

	if ( firstMCatNLOEmission ) {

	if ( !isMCatNLOHEvent )
	nEmissions = 1;
	else
	nEmissions = 0;

	}

	if ( !firstMCatNLOEmission ) {

	doCascade(nEmitted,optHardPt,optCutoff);

	if ( discardNoEmissions ) {
	if ( !didRadiate )
	throw Veto();
	if ( nEmissions )
	if ( nEmissions < nEmitted )
	throw Veto();
	}

	} else {

	if ( nEmissions == 1 )
	doCascade(nEmitted,optHardPt,optCutoff);

	}

	if ( intrinsicPtGenerator ) {
	if ( eventRecord().incoming().first->coloured() &&
	eventRecord().incoming().second->coloured() ) {
	SpinOneLorentzRotation rot =
	intrinsicPtGenerator->kick(eventRecord().incoming(),
	eventRecord().intermediates());
	eventRecord().transform(rot);
	}
	}

	didRealign = realign();

	constituentReshuffle();

	break;

	} catch (RedoShower&) {

	if ( ++nTries > maxtry() )
	throw ShowerTriesVeto(maxtry());

	eventRecord().clear();
	eventRecord().prepare(sub,dynamic_ptr_cast<tStdXCombPtr>(lastXCombPtr()),pdfs());

	continue;

	} catch (...) {
	throw;
	}

	}

	+ tEventPtr event = eventHandler()->currentEvent();
	+
	+ if ( event->optionalWeights().empty() ) {
	+ event->optionalWeights() = currentWeights;
	+ } else {
	+ map<string,double> combineMap;
	+ for ( map<string,double>::const_iterator w =
	+ event->optionalWeights().begin();
	+ w != event->optionalWeights().end(); ++w ) {
	+ for ( map<string,double>::const_iterator v = currentWeights.begin();
	+ v != currentWeights.end(); ++v ) {
	+ combineMap[w->first + "-" + v->first] = w->second*v->second;
	+ }
	+ }
	+ event->optionalWeights() = combineMap;
	+ }
	+
	return eventRecord().fillEventRecord(newStep(),firstInteraction(),didRealign);

	}

	void DipoleShowerHandler::constituentReshuffle() {

	if ( constituentReshuffler ) {
	constituentReshuffler->reshuffle(eventRecord().outgoing(),
	eventRecord().incoming(),
	eventRecord().intermediates());
	}

	}

	void DipoleShowerHandler::hardScales(Energy2 muf) {

	maxPt = generator()->maximumCMEnergy();

	if ( restrictPhasespace() ) {
	if ( !hardScaleIsMuF() \|\| !firstInteraction() ) {
	if ( !eventRecord().outgoing().empty() ) {
	for ( PList::const_iterator p = eventRecord().outgoing().begin();
	p != eventRecord().outgoing().end(); ++p )
	maxPt = min(maxPt,(**p).momentum().mt());
	} else {
	assert(!eventRecord().hard().empty());
	Lorentz5Momentum phard(ZERO,ZERO,ZERO,ZERO);
	for ( PList::const_iterator p = eventRecord().hard().begin();
	p != eventRecord().hard().end(); ++p )
	phard += (**p).momentum();
	Energy mhard = phard.m();
	maxPt = mhard;
	}
	maxPt *= hardScaleFactor();
	} else {
	maxPt = hardScaleFactor()*sqrt(muf);
	}
	muPt = maxPt;
	} else {
	muPt = hardScaleFactor()*sqrt(muf);
	}

	for ( list<DipoleChain>::iterator ch = eventRecord().chains().begin();
	ch != eventRecord().chains().end(); ++ch ) {

	Energy minVetoScale = -1.*GeV;

	for ( list<Dipole>::iterator dip = ch->dipoles().begin();
	dip != ch->dipoles().end(); ++dip ) {

	// max scale per config
	Energy maxFirst = 0.0*GeV;
	Energy maxSecond = 0.0*GeV;

	for ( vector<Ptr<DipoleSplittingKernel>::ptr>::iterator k =
	kernels.begin(); k != kernels.end(); ++k ) {

	pair<bool,bool> conf = make_pair(true,false);

	if ( (**k).canHandle(dip->index(conf)) ) {
	Energy scale =
	evolutionOrdering()->hardScale(dip->emitter(conf),dip->spectator(conf),
	dip->emitterX(conf),dip->spectatorX(conf),
	**k,dip->index(conf));
	maxFirst = max(maxFirst,scale);
	}

	conf = make_pair(false,true);

	if ( (**k).canHandle(dip->index(conf)) ) {
	Energy scale =
	evolutionOrdering()->hardScale(dip->emitter(conf),dip->spectator(conf),
	dip->emitterX(conf),dip->spectatorX(conf),
	**k,dip->index(conf));
	maxSecond = max(maxSecond,scale);
	}

	}

	if ( dip->leftParticle()->vetoScale() >= ZERO ) {
	maxFirst = min(maxFirst,sqrt(dip->leftParticle()->vetoScale()));
	if ( minVetoScale >= ZERO )
	minVetoScale = min(minVetoScale,sqrt(dip->leftParticle()->vetoScale()));
	else
	minVetoScale = sqrt(dip->leftParticle()->vetoScale());
	}

	if ( dip->rightParticle()->vetoScale() >= ZERO ) {
	maxSecond = min(maxSecond,sqrt(dip->rightParticle()->vetoScale()));
	if ( minVetoScale >= ZERO )
	minVetoScale = min(minVetoScale,sqrt(dip->rightParticle()->vetoScale()));
	else
	minVetoScale = sqrt(dip->rightParticle()->vetoScale());
	}

	maxFirst = min(maxPt,maxFirst);
	dip->emitterScale(make_pair(true,false),maxFirst);

	maxSecond = min(maxPt,maxSecond);
	dip->emitterScale(make_pair(false,true),maxSecond);

	}

	if ( !evolutionOrdering()->independentDipoles() &&
	chainOrderVetoScales &&
	minVetoScale >= ZERO ) {
	for ( list<Dipole>::iterator dip = ch->dipoles().begin();
	dip != ch->dipoles().end(); ++dip ) {
	dip->leftScale(min(dip->leftScale(),minVetoScale));
	dip->rightScale(min(dip->rightScale(),minVetoScale));
	}
	}

	}

	}

	Energy DipoleShowerHandler::getWinner(DipoleSplittingInfo& winner,
	const Dipole& dip,
	pair<bool,bool> conf,
	Energy optHardPt,
	Energy optCutoff) {
	return
	getWinner(winner,dip.index(conf),
	dip.emitterX(conf),dip.spectatorX(conf),
	conf,dip.emitter(conf),dip.spectator(conf),
	dip.emitterScale(conf),optHardPt,optCutoff);
	}

	Energy DipoleShowerHandler::getWinner(SubleadingSplittingInfo& winner,
	Energy optHardPt,
	Energy optCutoff) {
	return
	getWinner(winner,winner.index(),
	winner.emitterX(),winner.spectatorX(),
	winner.configuration(),
	winner.emitter(),winner.spectator(),
	winner.startScale(),optHardPt,optCutoff);
	}

	Energy DipoleShowerHandler::getWinner(DipoleSplittingInfo& winner,
	const DipoleIndex& index,
	double emitterX, double spectatorX,
	pair<bool,bool> conf,
	tPPtr emitter, tPPtr spectator,
	Energy startScale,
	Energy optHardPt,
	Energy optCutoff) {

	if ( !index.initialStateEmitter() &&
	!doFSR ) {
	winner.didStopEvolving();
	return 0.0*GeV;
	}

	if ( index.initialStateEmitter() &&
	!doISR ) {
	winner.didStopEvolving();
	return 0.0*GeV;
	}

	DipoleSplittingInfo candidate;
	candidate.index(index);
	candidate.configuration(conf);
	candidate.emitterX(emitterX);
	candidate.spectatorX(spectatorX);

	if ( generators().find(candidate.index()) == generators().end() )
	getGenerators(candidate.index());

	//
	// NOTE -- needs proper fixing at some point
	//
	// For some very strange reason, equal_range gives back
	// key ranges it hasn't been asked for. This particularly
	// happens e.g. for FI dipoles of the same kind, but different
	// PDF (hard vs MPI PDF). I can't see a reason for this,
	// as DipoleIndex properly implements comparison for equality
	// and (lexicographic) ordering; for the time being, we
	// use equal_range, extented by an explicit check for wether
	// the key is indeed what we wanted. See line after (*) comment
	// below.
	//

	pair<GeneratorMap::iterator,GeneratorMap::iterator> gens
	= generators().equal_range(candidate.index());

	Energy winnerScale = 0.0*GeV;
	GeneratorMap::iterator winnerGen = generators().end();

	for ( GeneratorMap::iterator gen = gens.first; gen != gens.second; ++gen ) {

	// (*) see NOTE above
	if ( !(gen->first == candidate.index()) )
	continue;

	if ( startScale <= gen->second->splittingKinematics()->IRCutoff() )
	continue;

	Energy dScale =
	gen->second->splittingKinematics()->dipoleScale(emitter->momentum(),
	spectator->momentum());

	// in very exceptional cases happening in DIS
	if ( isnan(dScale/GeV ) )
	throw RedoShower();

	candidate.scale(dScale);
	candidate.continuesEvolving();
	Energy hardScale = evolutionOrdering()->maxPt(startScale,candidate,*(gen->second->splittingKernel()));

	Energy maxPossible =
	gen->second->splittingKinematics()->ptMax(candidate.scale(),
	candidate.emitterX(), candidate.spectatorX(),
	candidate.index(),
	*gen->second->splittingKernel());

	Energy ircutoff =
	optCutoff < gen->second->splittingKinematics()->IRCutoff() ?
	gen->second->splittingKinematics()->IRCutoff() :
	optCutoff;

	if ( maxPossible <= ircutoff ) {
	continue;
	}

	if ( maxPossible >= hardScale )
	candidate.hardPt(hardScale);
	else {
	hardScale = maxPossible;
	candidate.hardPt(maxPossible);
	}

	- gen->second->generate(candidate,optHardPt,optCutoff);
	+ gen->second->generate(candidate,currentWeights,optHardPt,optCutoff);
	Energy nextScale = evolutionOrdering()->evolutionScale(gen->second->lastSplitting(),*(gen->second->splittingKernel()));

	if ( nextScale > winnerScale ) {
	winner.fill(candidate);
	gen->second->completeSplitting(winner);
	winnerGen = gen;
	winnerScale = nextScale;
	}

	}

	if ( winnerGen == generators().end() ) {
	winner.didStopEvolving();
	return 0.0*GeV;
	}

	if ( winner.stoppedEvolving() )
	return 0.0*GeV;

	return winnerScale;

	}

	void DipoleShowerHandler::doCascade(unsigned int& emDone,
	Energy optHardPt,
	Energy optCutoff) {

	if ( nEmissions )
	if ( emDone == nEmissions )
	return;

	DipoleSplittingInfo winner;
	DipoleSplittingInfo dipoleWinner;

	while ( eventRecord().haveChain() ) {

	if ( verbosity > 2 ) {
	generator()->log() << "DipoleShowerHandler selecting splittings for the chain:\n"
	<< eventRecord().currentChain() << flush;
	}

	list<Dipole>::iterator winnerDip = eventRecord().currentChain().dipoles().end();

	Energy winnerScale = 0.0*GeV;
	Energy nextLeftScale = 0.0*GeV;
	Energy nextRightScale = 0.0*GeV;

	for ( list<Dipole>::iterator dip = eventRecord().currentChain().dipoles().begin();
	dip != eventRecord().currentChain().dipoles().end(); ++dip ) {

	nextLeftScale = getWinner(dipoleWinner,*dip,make_pair(true,false),optHardPt,optCutoff);

	if ( nextLeftScale > winnerScale ) {
	winnerScale = nextLeftScale;
	winner = dipoleWinner;
	winnerDip = dip;
	}

	nextRightScale = getWinner(dipoleWinner,*dip,make_pair(false,true),optHardPt,optCutoff);

	if ( nextRightScale > winnerScale ) {
	winnerScale = nextRightScale;
	winner = dipoleWinner;
	winnerDip = dip;
	}

	if ( evolutionOrdering()->independentDipoles() ) {
	Energy dipScale = max(nextLeftScale,nextRightScale);
	if ( dip->leftScale() > dipScale )
	dip->leftScale(dipScale);
	if ( dip->rightScale() > dipScale )
	dip->rightScale(dipScale);
	}

	}

	if ( verbosity > 1 ) {
	if ( winnerDip != eventRecord().currentChain().dipoles().end() )
	generator()->log() << "DipoleShowerHandler selected the splitting:\n"
	<< winner << " for the dipole\n"
	<< (*winnerDip) << flush;
	else
	generator()->log() << "DipoleShowerHandler could not select a splitting above the IR cutoff\n"
	<< flush;
	}

	// pop the chain if no dipole did radiate
	if ( winnerDip == eventRecord().currentChain().dipoles().end() ) {
	if ( theEventReweight ) {
	double w = theEventReweight->weightNoEmission(eventRecord().incoming(),
	eventRecord().outgoing(),
	eventRecord().hard(),theGlobalAlphaS);
	Ptr<StandardEventHandler>::tptr eh =
	dynamic_ptr_cast<Ptr<StandardEventHandler>::tptr>(eventHandler());
	assert(eh);
	eh->reweight(w);
	}
	eventRecord().popChain();
	continue;
	}

	// otherwise perform the splitting

	didRadiate = true;

	isMCatNLOSEvent = false;
	isMCatNLOHEvent = false;

	pair<list<Dipole>::iterator,list<Dipole>::iterator> children;

	DipoleChain* firstChain = 0;
	DipoleChain* secondChain = 0;

	eventRecord().split(winnerDip,winner,children,firstChain,secondChain);

	assert(firstChain && secondChain);

	evolutionOrdering()->setEvolutionScale(winnerScale,winner,*firstChain,children);

	if ( !secondChain->dipoles().empty() )
	evolutionOrdering()->setEvolutionScale(winnerScale,winner,*secondChain,children);

	if ( verbosity > 1 ) {
	generator()->log() << "DipoleShowerHandler did split the last selected dipole into:\n"
	<< (children.first) << (children.second) << flush;
	}

	if ( verbosity > 2 ) {
	generator()->log() << "After splitting the last selected dipole, "
	<< "DipoleShowerHandler encountered the following chains:\n"
	<< (firstChain) << (secondChain) << flush;
	}

	if ( theEventReweight ) {
	double w = theEventReweight->weight(eventRecord().incoming(),
	eventRecord().outgoing(),
	eventRecord().hard(),theGlobalAlphaS);
	Ptr<StandardEventHandler>::tptr eh =
	dynamic_ptr_cast<Ptr<StandardEventHandler>::tptr>(eventHandler());
	assert(eh);
	eh->reweight(w);
	}

	if ( nEmissions )
	if ( ++emDone == nEmissions )
	return;

	}

	}

	bool DipoleShowerHandler::realign() {

	if ( !didRadiate && !intrinsicPtGenerator )
	return false;

	if ( eventRecord().incoming().first->coloured() \|\|
	eventRecord().incoming().second->coloured() ) {

	if ( eventRecord().incoming().first->momentum().perp2()/GeV2 < 1e-10 &&
	eventRecord().incoming().second->momentum().perp2()/GeV2 < 1e-10 )
	return false;

	pair<Lorentz5Momentum,Lorentz5Momentum> inMomenta
	(eventRecord().incoming().first->momentum(),
	eventRecord().incoming().second->momentum());

	SpinOneLorentzRotation transform((inMomenta.first+inMomenta.second).findBoostToCM());

	Axis dir = (transform * inMomenta.first).vect().unit();
	Axis rot (-dir.y(),dir.x(),0);
	double theta = dir.theta();

	if ( lastParticles().first->momentum().z() < ZERO )
	theta = -theta;

	transform.rotate(-theta,rot);

	inMomenta.first = transform*inMomenta.first;
	inMomenta.second = transform*inMomenta.second;

	assert(inMomenta.first.z() > ZERO &&
	inMomenta.second.z() < ZERO);

	Energy2 sHat =
	(eventRecord().incoming().first->momentum() +
	eventRecord().incoming().second->momentum()).m2();

	pair<Energy,Energy> masses(eventRecord().incoming().first->mass(),
	eventRecord().incoming().second->mass());
	pair<Energy,Energy> qs;

	if ( !eventRecord().incoming().first->coloured() ) {
	assert(masses.second == ZERO);
	qs.first = eventRecord().incoming().first->momentum().z();
	qs.second = (sHat-sqr(masses.first))/(2.*(qs.first+sqrt(sqr(masses.first)+sqr(qs.first))));
	} else if ( !eventRecord().incoming().second->coloured() ) {
	assert(masses.first == ZERO);
	qs.second = eventRecord().incoming().second->momentum().z();
	qs.first = (sHat-sqr(masses.second))/(2.*(qs.second+sqrt(sqr(masses.second)+sqr(qs.second))));
	} else {
	assert(masses.first == ZERO && masses.second == ZERO);
	if ( realignmentScheme == 0 ) {
	double yX = eventRecord().pX().rapidity();
	double yInt = (transform*eventRecord().pX()).rapidity();
	double dy = yX-yInt;
	qs.first = (sqrt(sHat)/2.)*exp(dy);
	qs.second = (sqrt(sHat)/2.)*exp(-dy);
	} else if ( realignmentScheme == 1 ) {
	Energy sS = sqrt((lastParticles().first->momentum() +
	lastParticles().second->momentum()).m2());
	qs.first = eventRecord().fractions().first * sS / 2.;
	qs.second = eventRecord().fractions().second * sS / 2.;
	}
	}

	double beta =
	(qs.first-qs.second) /
	( sqrt(sqr(masses.first)+sqr(qs.first)) +
	sqrt(sqr(masses.second)+sqr(qs.second)) );
	transform.boostZ(beta);

	Lorentz5Momentum tmp;

	if ( eventRecord().incoming().first->coloured() ) {
	tmp = eventRecord().incoming().first->momentum();
	tmp = transform * tmp;
	eventRecord().incoming().first->set5Momentum(tmp);
	}
	if ( eventRecord().incoming().second->coloured() ) {
	tmp = eventRecord().incoming().second->momentum();
	tmp = transform * tmp;
	eventRecord().incoming().second->set5Momentum(tmp);
	}
	eventRecord().transform(transform);
	return true;

	}

	return false;

	}

	void DipoleShowerHandler::resetAlphaS(Ptr<AlphaSBase>::tptr as) {

	for ( vector<Ptr<DipoleSplittingKernel>::ptr>::iterator k = kernels.begin();
	k != kernels.end(); ++k ) {
	(**k).alphaS(as);
	(**k).renormalizationScaleFreeze(theRenormalizationScaleFreeze);
	(**k).factorizationScaleFreeze(theFactorizationScaleFreeze);
	}

	// clear the generators to be rebuild
	// actually, there shouldn't be any generators
	// when this happens.
	generators().clear();

	}

	void DipoleShowerHandler::resetReweight(Ptr<DipoleSplittingReweight>::tptr rw) {
	for ( GeneratorMap::iterator k = generators().begin();
	k != generators().end(); ++k )
	k->second->splittingReweight(rw);
	}

	void DipoleShowerHandler::getGenerators(const DipoleIndex& ind,
	Ptr<DipoleSplittingReweight>::tptr rw) {

	bool gotone = false;

	for ( vector<Ptr<DipoleSplittingKernel>::ptr>::iterator k =
	kernels.begin(); k != kernels.end(); ++k ) {
	if ( (**k).canHandle(ind) ) {

	if ( verbosity > 0 ) {
	generator()->log() << "DipoleShowerHandler encountered the dipole configuration\n"
	<< ind << " in event number "
	<< eventHandler()->currentEvent()->number()
	<< "\nwhich can be handled by the splitting kernel '"
	<< (**k).name() << "'.\n" << flush;
	}

	gotone = true;

	Ptr<DipoleSplittingGenerator>::ptr nGenerator =
	new_ptr(DipoleSplittingGenerator());
	nGenerator->doCompensate(theDoCompensate);
	nGenerator->splittingKernel(*k);
	nGenerator->splittingKernel()->renormalizationScaleFactor(renormalizationScaleFactor());
	nGenerator->splittingKernel()->factorizationScaleFactor(factorizationScaleFactor());
	nGenerator->splittingKernel()->freezeGrid(theFreezeGrid);

	GeneratorMap::const_iterator equivalent = generators().end();

	for ( GeneratorMap::const_iterator eq = generators().begin();
	eq != generators().end(); ++eq ) {
	if ( !eq->second->wrapping() )
	if ( (*k).canHandleEquivalent(ind,(eq->second->splittingKernel()),eq->first) ) {

	equivalent = eq;

	if ( verbosity > 0 ) {
	generator()->log() << "The dipole configuration "
	<< ind
	<< " can equivalently be handled by the existing\n"
	<< "generator for configuration "
	<< eq->first << " using the kernel '"
	<< eq->second->splittingKernel()->name()
	<< "'\n" << flush;
	}

	break;

	}
	}

	if ( equivalent != generators().end() ) {
	nGenerator->wrap(equivalent->second);
	}

	DipoleSplittingInfo dummy;
	dummy.index(ind);
	nGenerator->splittingReweight(rw);
	nGenerator->prepare(dummy);

	generators().insert(make_pair(ind,nGenerator));


	}
	}

	if ( !gotone ) {
	generator()->logWarning(Exception()
	<< "DipoleShowerHandler could not "
	<< "find a splitting kernel which is able "
	<< "to handle splittings off the dipole "
	<< ind << ".\n"
	<< "Please check the input files."
	<< Exception::warning);
	}

	}

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

	void DipoleShowerHandler::doinit() {
	ShowerHandler::doinit();
	if ( theGlobalAlphaS )
	resetAlphaS(theGlobalAlphaS);
	}

	void DipoleShowerHandler::dofinish() {
	ShowerHandler::dofinish();
	}

	void DipoleShowerHandler::doinitrun() {
	ShowerHandler::doinitrun();
	}

	void DipoleShowerHandler::persistentOutput(PersistentOStream & os) const {
	os << kernels << theEvolutionOrdering
	<< constituentReshuffler << intrinsicPtGenerator
	<< theGlobalAlphaS << chainOrderVetoScales
	<< nEmissions << discardNoEmissions << firstMCatNLOEmission << doFSR << doISR
	<< realignmentScheme << verbosity << printEvent
	<< ounit(theRenormalizationScaleFreeze,GeV)
	<< ounit(theFactorizationScaleFreeze,GeV)
	<< isMCatNLOSEvent << isMCatNLOHEvent << theShowerApproximation
	<< theDoCompensate << theFreezeGrid
	<< theEventReweight << ounit(maxPt,GeV)
	<< ounit(muPt,GeV);
	}

	void DipoleShowerHandler::persistentInput(PersistentIStream & is, int) {
	is >> kernels >> theEvolutionOrdering
	>> constituentReshuffler >> intrinsicPtGenerator
	>> theGlobalAlphaS >> chainOrderVetoScales
	>> nEmissions >> discardNoEmissions >> firstMCatNLOEmission >> doFSR >> doISR
	>> realignmentScheme >> verbosity >> printEvent
	>> iunit(theRenormalizationScaleFreeze,GeV)
	>> iunit(theFactorizationScaleFreeze,GeV)
	>> isMCatNLOSEvent >> isMCatNLOHEvent >> theShowerApproximation
	>> theDoCompensate >> theFreezeGrid
	>> theEventReweight >> iunit(maxPt,GeV)
	>> iunit(muPt,GeV);
	}

	ClassDescription<DipoleShowerHandler> DipoleShowerHandler::initDipoleShowerHandler;
	// Definition of the static class description member.

	void DipoleShowerHandler::Init() {

	static ClassDocumentation<DipoleShowerHandler> documentation
	("The DipoleShowerHandler class manages the showering using "
	"the dipole shower algorithm.",
	"The shower evolution was performed using the algorithm described in "
	"\\cite{Platzer:2009jq} and \\cite{Platzer:2011bc}.",
	"%\\cite{Platzer:2009jq}\n"
	"\\bibitem{Platzer:2009jq}\n"
	"S.~Platzer and S.~Gieseke,\n"
	"``Coherent Parton Showers with Local Recoils,''\n"
	" JHEP {\\bf 1101}, 024 (2011)\n"
	"arXiv:0909.5593 [hep-ph].\n"
	"%%CITATION = ARXIV:0909.5593;%%\n"
	"%\\cite{Platzer:2011bc}\n"
	"\\bibitem{Platzer:2011bc}\n"
	"S.~Platzer and S.~Gieseke,\n"
	"``Dipole Showers and Automated NLO Matching in Herwig,''\n"
	"arXiv:1109.6256 [hep-ph].\n"
	"%%CITATION = ARXIV:1109.6256;%%");

	static RefVector<DipoleShowerHandler,DipoleSplittingKernel> interfaceKernels
	("Kernels",
	"Set the splitting kernels to be used by the dipole shower.",
	&DipoleShowerHandler::kernels, -1, false, false, true, false, false);


	static Reference<DipoleShowerHandler,DipoleEvolutionOrdering> interfaceEvolutionOrdering
	("EvolutionOrdering",
	"Set the evolution ordering to be used.",
	&DipoleShowerHandler::theEvolutionOrdering, false, false, true, false, false);


	static Reference<DipoleShowerHandler,ConstituentReshuffler> interfaceConstituentReshuffler
	("ConstituentReshuffler",
	"The object to be used to reshuffle partons to their constitutent mass shells.",
	&DipoleShowerHandler::constituentReshuffler, false, false, true, true, false);


	static Reference<DipoleShowerHandler,IntrinsicPtGenerator> interfaceIntrinsicPtGenerator
	("IntrinsicPtGenerator",
	"Set the object in charge to generate intrinsic pt for incoming partons.",
	&DipoleShowerHandler::intrinsicPtGenerator, false, false, true, true, false);

	static Reference<DipoleShowerHandler,AlphaSBase> interfaceGlobalAlphaS
	("GlobalAlphaS",
	"Set a global strong coupling for all splitting kernels.",
	&DipoleShowerHandler::theGlobalAlphaS, false, false, true, true, false);


	static Switch<DipoleShowerHandler,bool> interfaceDoFSR
	("DoFSR",
	"Switch on or off final state radiation.",
	&DipoleShowerHandler::doFSR, true, false, false);
	static SwitchOption interfaceDoFSROn
	(interfaceDoFSR,
	"On",
	"Switch on final state radiation.",
	true);
	static SwitchOption interfaceDoFSROff
	(interfaceDoFSR,
	"Off",
	"Switch off final state radiation.",
	false);

	static Switch<DipoleShowerHandler,bool> interfaceDoISR
	("DoISR",
	"Switch on or off initial state radiation.",
	&DipoleShowerHandler::doISR, true, false, false);
	static SwitchOption interfaceDoISROn
	(interfaceDoISR,
	"On",
	"Switch on initial state radiation.",
	true);
	static SwitchOption interfaceDoISROff
	(interfaceDoISR,
	"Off",
	"Switch off initial state radiation.",
	false);

	static Switch<DipoleShowerHandler,int> interfaceRealignmentScheme
	("RealignmentScheme",
	"The realignment scheme to use.",
	&DipoleShowerHandler::realignmentScheme, 0, false, false);
	static SwitchOption interfaceRealignmentSchemePreserveRapidity
	(interfaceRealignmentScheme,
	"PreserveRapidity",
	"Preserve the rapidity of non-coloured outgoing system.",
	0);
	static SwitchOption interfaceRealignmentSchemeEvolutionFractions
	(interfaceRealignmentScheme,
	"EvolutionFractions",
	"Use momentum fractions as generated by the evolution.",
	1);
	static SwitchOption interfaceRealignmentSchemeCollisionFrame
	(interfaceRealignmentScheme,
	"CollisionFrame",
	"Determine realignment from collision frame.",
	2);


	static Switch<DipoleShowerHandler,bool> interfaceChainOrderVetoScales
	("ChainOrderVetoScales",
	"[experimental] Switch on or off the chain ordering for veto scales.",
	&DipoleShowerHandler::chainOrderVetoScales, true, false, false);
	static SwitchOption interfaceChainOrderVetoScalesOn
	(interfaceChainOrderVetoScales,
	"On",
	"Switch on chain ordering for veto scales.",
	true);
	static SwitchOption interfaceChainOrderVetoScalesOff
	(interfaceChainOrderVetoScales,
	"Off",
	"Switch off chain ordering for veto scales.",
	false);

	interfaceChainOrderVetoScales.rank(-1);


	static Parameter<DipoleShowerHandler,unsigned int> interfaceNEmissions
	("NEmissions",
	"[debug option] Limit the number of emissions to be generated. Zero does not limit the number of emissions.",
	&DipoleShowerHandler::nEmissions, 0, 0, 0,
	false, false, Interface::lowerlim);

	interfaceNEmissions.rank(-1);


	static Switch<DipoleShowerHandler,bool> interfaceDiscardNoEmissions
	("DiscardNoEmissions",
	"[debug option] Discard events without radiation.",
	&DipoleShowerHandler::discardNoEmissions, false, false, false);
	static SwitchOption interfaceDiscardNoEmissionsOn
	(interfaceDiscardNoEmissions,
	"On",
	"Discard events without radiation.",
	true);
	static SwitchOption interfaceDiscardNoEmissionsOff
	(interfaceDiscardNoEmissions,
	"Off",
	"Do not discard events without radiation.",
	false);

	interfaceDiscardNoEmissions.rank(-1);

	static Switch<DipoleShowerHandler,bool> interfaceFirstMCatNLOEmission
	("FirstMCatNLOEmission",
	"[debug option] Only perform the first MC@NLO emission.",
	&DipoleShowerHandler::firstMCatNLOEmission, false, false, false);
	static SwitchOption interfaceFirstMCatNLOEmissionOn
	(interfaceFirstMCatNLOEmission,
	"On",
	"",
	true);
	static SwitchOption interfaceFirstMCatNLOEmissionOff
	(interfaceFirstMCatNLOEmission,
	"Off",
	"",
	false);

	interfaceFirstMCatNLOEmission.rank(-1);


	static Parameter<DipoleShowerHandler,int> interfaceVerbosity
	("Verbosity",
	"[debug option] Set the level of debug information provided.",
	&DipoleShowerHandler::verbosity, 0, 0, 0,
	false, false, Interface::lowerlim);

	interfaceVerbosity.rank(-1);


	static Parameter<DipoleShowerHandler,int> interfacePrintEvent
	("PrintEvent",
	"[debug option] The number of events for which debugging information should be provided.",
	&DipoleShowerHandler::printEvent, 0, 0, 0,
	false, false, Interface::lowerlim);

	interfacePrintEvent.rank(-1);

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

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

	static Switch<DipoleShowerHandler,bool> interfaceDoCompensate
	("DoCompensate",
	"",
	&DipoleShowerHandler::theDoCompensate, false, false, false);
	static SwitchOption interfaceDoCompensateYes
	(interfaceDoCompensate,
	"Yes",
	"",
	true);
	static SwitchOption interfaceDoCompensateNo
	(interfaceDoCompensate,
	"No",
	"",
	false);

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

	static Reference<DipoleShowerHandler,DipoleEventReweight> interfaceEventReweight
	("EventReweight",
	"",
	&DipoleShowerHandler::theEventReweight, false, false, true, true, false);

	}

	diff --git a/DipoleShower/DipoleShowerHandler.h b/DipoleShower/DipoleShowerHandler.h
	--- a/DipoleShower/DipoleShowerHandler.h
	+++ b/DipoleShower/DipoleShowerHandler.h
	@@ -1,504 +1,509 @@
	// -- C++ --
	//
	// DipoleShowerHandler.h 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.
	//
	#ifndef HERWIG_DipoleShowerHandler_H
	#define HERWIG_DipoleShowerHandler_H
	//
	// This is the declaration of the DipoleShowerHandler class.
	//

	#include "Herwig/Shower/ShowerHandler.h"

	#include "Herwig/DipoleShower/Base/DipoleSplittingInfo.h"
	#include "Herwig/DipoleShower/Base/DipoleSplittingReweight.h"
	#include "Herwig/DipoleShower/Kernels/DipoleSplittingKernel.h"
	#include "Herwig/DipoleShower/Base/DipoleSplittingGenerator.h"
	#include "Herwig/DipoleShower/Base/DipoleEventRecord.h"
	#include "Herwig/DipoleShower/Base/DipoleEvolutionOrdering.h"
	#include "Herwig/DipoleShower/Base/DipoleEventReweight.h"
	#include "Herwig/DipoleShower/Utility/ConstituentReshuffler.h"
	#include "Herwig/DipoleShower/Utility/IntrinsicPtGenerator.h"
	#include "Herwig/MatrixElement/Matchbox/Matching/ShowerApproximation.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* \ingroup DipoleShower
	* \author Simon Platzer
	*
	* \brief The DipoleShowerHandler class manages the showering using
	* the dipole shower algorithm.
	*
	* @see \ref DipoleShowerHandlerInterfaces "The interfaces"
	* defined for DipoleShowerHandler.
	*/
	class DipoleShowerHandler: public ShowerHandler {

	public:

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

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

	public:

	/**
	* Indicate a problem in the shower.
	*/
	struct RedoShower {};

	/**
	* Insert an additional splitting kernel.
	*/
	void addSplitting(Ptr<DipoleSplittingKernel>::ptr sp) {
	kernels.push_back(sp);
	}

	/**
	* Reset the alpha_s for all splitting kernels.
	*/
	void resetAlphaS(Ptr<AlphaSBase>::tptr);

	/**
	* Reset the splitting reweight for all splitting kernels.
	*/
	void resetReweight(Ptr<DipoleSplittingReweight>::tptr);

	/**
	* Return true, if the shower handler can generate a truncated
	* shower for POWHEG style events generated using Matchbox
	*/
	virtual bool canHandleMatchboxTrunc() const { return false; }

	/**
	* Return true, if this cascade handler will perform reshuffling from hard
	* process masses.
	*/
	virtual bool isReshuffling() const { return false; }

	/**
	* Return the relevant hard scale to be used in the profile scales
	*/
	virtual Energy hardScale() const {
	return muPt;
	}

	protected:

	typedef multimap<DipoleIndex,Ptr<DipoleSplittingGenerator>::ptr> GeneratorMap;

	/**
	* The main method which manages the showering of a subprocess.
	*/
	virtual tPPair cascade(tSubProPtr sub, XCPtr xcomb) {
	return cascade(sub,xcomb,ZERO,ZERO);
	}

	/**
	* The main method which manages the showering of a subprocess.
	*/
	tPPair cascade(tSubProPtr sub, XCPtr xcomb,
	Energy optHardPt, Energy optCutoff);

	/**
	* Build splitting generators for the given
	* dipole index.
	*/
	void getGenerators(const DipoleIndex&,
	Ptr<DipoleSplittingReweight>::tptr rw =
	Ptr<DipoleSplittingReweight>::tptr());

	/**
	* Setup the hard scales.
	*/
	void hardScales(Energy2 scale);

	/**
	* Return the evolution ordering
	*/
	Ptr<DipoleEvolutionOrdering>::tptr evolutionOrdering() const { return theEvolutionOrdering; }

	/**
	* Reshuffle to constituent mass shells
	*/
	void constituentReshuffle();

	/**
	* Access the generator map
	*/
	GeneratorMap& generators() { return theGenerators; }

	/**
	* Access the event record
	*/
	DipoleEventRecord& eventRecord() { return theEventRecord; }

	/**
	* Return the event record
	*/
	const DipoleEventRecord& eventRecord() const { return theEventRecord; }

	/**
	* Return the splitting kernels.
	*/
	const vector<Ptr<DipoleSplittingKernel>::ptr>& splittingKernels() const {
	return kernels;
	}

	/**
	* Realign the event such as to have the incoming partons along thre
	* beam axes.
	*/
	bool realign();

	private:

	/**
	* Perform the cascade.
	*/
	void doCascade(unsigned int& emDone,
	Energy optHardPt = ZERO,
	Energy optCutoff = ZERO);

	/**
	* Get the winning splitting for the
	* given dipole and configuration.
	*/
	Energy getWinner(DipoleSplittingInfo& winner,
	const Dipole& dip,
	pair<bool,bool> conf,
	Energy optHardPt = ZERO,
	Energy optCutoff = ZERO);

	/**
	* Get the winning splitting for the
	* given dipole and configuration.
	*/
	Energy getWinner(SubleadingSplittingInfo& winner,
	Energy optHardPt = ZERO,
	Energy optCutoff = ZERO);

	/**
	* Get the winning splitting for the
	* given dipole and configuration.
	*/
	Energy getWinner(DipoleSplittingInfo& winner,
	const DipoleIndex& index,
	double emitterX, double spectatorX,
	pair<bool,bool> conf,
	tPPtr emitter, tPPtr spectator,
	Energy startScale,
	Energy optHardPt = ZERO,
	Energy optCutoff = ZERO);

	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).


	protected:

	/** @name Standard Interfaced functions. */
	//@{
	/**
	* Initialize this object after the setup phase before saving an
	* EventGenerator to disk.
	* @throws InitException if object could not be initialized properly.
	*/
	virtual void doinit();

	/**
	* Initialize this object. Called in the run phase just before
	* a run begins.
	*/
	virtual void doinitrun();

	/**
	* Finalize this object. Called in the run phase just after a
	* run has ended. Used eg. to write out statistics.
	*/
	virtual void dofinish();
	//@}


	private:

	/**
	* The splitting kernels to be used.
	*/
	vector<Ptr<DipoleSplittingKernel>::ptr> kernels;

	/**
	* The evolution ordering considered
	*/
	Ptr<DipoleEvolutionOrdering>::ptr theEvolutionOrdering;

	/**
	* The ConstituentReshuffler to be used
	*/
	Ptr<ConstituentReshuffler>::ptr constituentReshuffler;

	/**
	* The intrinsic pt generator to be used.
	*/
	Ptr<IntrinsicPtGenerator>::ptr intrinsicPtGenerator;

	/**
	* A global alpha_s to be used for all splitting kernels.
	*/
	Ptr<AlphaSBase>::ptr theGlobalAlphaS;

	/**
	* Apply chain ordering to events from matrix
	* element corrections.
	*/
	bool chainOrderVetoScales;

	/**
	* Limit the number of emissions.
	* Limit applied if > 0.
	*/
	unsigned int nEmissions;

	/**
	* Discard events which did not radiate.
	*/
	bool discardNoEmissions;

	/**
	* Perform the first MC@NLO emission only.
	*/
	bool firstMCatNLOEmission;

	/**
	* Switch on or off final state radiation.
	*/
	bool doFSR;

	/**
	* Switch on or off initial state radiation.
	*/
	bool doISR;

	/**
	* The realignment scheme
	*/
	int realignmentScheme;

	private:

	/**
	* The verbosity level.
	* 0 - print no info
	* 1 - print diagnostic information on setting up
	* splitting generators etc.
	* 2 - print detailed event information for up to
	* printEvent events.
	* 3 - print dipole chains after each splitting.
	*/
	int verbosity;

	/**
	* See verbosity.
	*/
	int printEvent;

	private:

	/**
	* The splitting generators indexed by the dipole
	* indices they can work on.
	*/
	GeneratorMap theGenerators;

	/**
	* The evnt record used.
	*/
	DipoleEventRecord theEventRecord;

	/**
	* The number of shoer tries so far.
	*/
	unsigned int nTries;

	/**
	* Whether or not we did radiate anything
	*/
	bool didRadiate;

	/**
	* Whether or not we did realign the event
	*/
	bool didRealign;

	private:

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

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

	/**
	* True, if we are showering on a MC@NLO S event
	*/
	bool isMCatNLOSEvent;

	/**
	* True, if we are showering on a MC@NLO H event
	*/
	bool isMCatNLOHEvent;

	/**
	* The matching subtraction, if appropriate
	*/
	Ptr<ShowerApproximation>::tptr theShowerApproximation;

	/**
	* True, if sampler should apply compensation
	*/
	bool theDoCompensate;

	/**
	* Return the number of accepted points after which the grid should
	* be frozen
	*/
	unsigned long theFreezeGrid;

	/**
	* A pointer to the dipole event reweight object
	*/
	Ptr<DipoleEventReweight>::ptr theEventReweight;

	/**
	* True if no warnings have been issued yet
	*/
	static bool firstWarn;

	/**
	* The shower starting scale for the last event encountered
	*/
	Energy maxPt;

	/**
	* The shower hard scale for the last event encountered
	*/
	Energy muPt;

	+ /**
	+ * The shower variation weights
	+ */
	+ map<string,double> currentWeights;
	+
	private:

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

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

	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

	/** This template specialization informs ThePEG about the name of
	* the DipoleShowerHandler class and the shared object where it is defined. */
	template <>
	struct ClassTraits<Herwig::DipoleShowerHandler>
	: public ClassTraitsBase<Herwig::DipoleShowerHandler> {
	/** Return a platform-independent class name */
	static string className() { return "Herwig::DipoleShowerHandler"; }
	/**
	* The name of a file containing the dynamic library where the class
	* DipoleShowerHandler is implemented. It may also include several, space-separated,
	* libraries if the class DipoleShowerHandler 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_DipoleShowerHandler_H */
	diff --git a/DipoleShower/Kernels/DipoleSplittingKernel.h b/DipoleShower/Kernels/DipoleSplittingKernel.h
	--- a/DipoleShower/Kernels/DipoleSplittingKernel.h
	+++ b/DipoleShower/Kernels/DipoleSplittingKernel.h
	@@ -1,439 +1,451 @@
	// -- C++ --
	//
	// DipoleSplittingKernel.h 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.
	//
	#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/DipoleShower/Utility/PDFRatio.h"
	#include "Herwig/DipoleShower/Base/DipoleSplittingInfo.h"
	#include "Herwig/DipoleShower/Kinematics/DipoleSplittingKinematics.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&) {}

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

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

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

	/**
	+ * 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.; }

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

	/**
	* Return the number of accepted points after which the grid should
	* be frozen
	*/
	unsigned long theFreezeGrid;

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

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

	};

	}

	#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 */
	diff --git a/MatrixElement/Matchbox/Matching/ShowerApproximationKernel.h b/MatrixElement/Matchbox/Matching/ShowerApproximationKernel.h
	--- a/MatrixElement/Matchbox/Matching/ShowerApproximationKernel.h
	+++ b/MatrixElement/Matchbox/Matching/ShowerApproximationKernel.h
	@@ -1,462 +1,478 @@
	// -- C++ --
	//
	// ShowerApproximationKernel.h 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.
	//
	#ifndef Herwig_ShowerApproximationKernel_H
	#define Herwig_ShowerApproximationKernel_H
	//
	// This is the declaration of the ShowerApproximationKernel class.
	//

	#include "ThePEG/Handlers/HandlerBase.h"
	#include "ThePEG/Handlers/StandardXComb.h"
	#include "Herwig/MatrixElement/Matchbox/Matching/ShowerApproximation.h"
	#include "Herwig/MatrixElement/Matchbox/Phasespace/InvertedTildeKinematics.h"
	#include "Herwig/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h"
	#include "Herwig/Sampling/exsample/exponential_generator.h"

	namespace Herwig {

	using namespace ThePEG;

	class ShowerApproximationGenerator;

	/**
	* \ingroup Matchbox
	* \author Simon Platzer
	*
	* \brief ShowerApproximationKernel generates emissions according to a
	* shower approximation entering a NLO matching.
	*
	*/
	class ShowerApproximationKernel: public HandlerBase {

	public:

	/**
	* Exception to communicate sampler maxtry events.
	*/
	struct MaxTryException {};

	public:

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

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

	public:

	/**
	* Set the XComb object describing the Born process
	*/
	void setBornXComb(tStdXCombPtr xc) { theBornXComb = xc; }

	/**
	* Return the XComb object describing the Born process
	*/
	tcStdXCombPtr bornCXComb() const { return theBornXComb; }

	/**
	* Return the XComb object describing the Born process
	*/
	tStdXCombPtr bornXComb() const { return theBornXComb; }

	/**
	* Set the XComb object describing the real emission process
	*/
	void setRealXComb(tStdXCombPtr xc) { theRealXComb = xc; }

	/**
	* Return the XComb object describing the real emission process
	*/
	tcStdXCombPtr realCXComb() const { return theRealXComb; }

	/**
	* Return the XComb object describing the real emission process
	*/
	tStdXCombPtr realXComb() const { return theRealXComb; }

	/**
	* Set the tilde xcomb objects associated to the real xcomb
	*/
	void setTildeXCombs(const vector<StdXCombPtr>& xc) { theTildeXCombs = xc; }

	/**
	* Return the tilde xcomb objects associated to the real xcomb
	*/
	const vector<StdXCombPtr>& tildeXCombs() const { return theTildeXCombs; }

	/**
	* Set the dipole in charge for the emission
	*/
	void setDipole(Ptr<SubtractionDipole>::tptr dip) { theDipole = dip; }

	/**
	* Return the dipole in charge for the emission
	*/
	Ptr<SubtractionDipole>::tptr dipole() const { return theDipole; }

	/**
	* Set the shower approximation.
	*/
	void showerApproximation(Ptr<ShowerApproximation>::tptr app) { theShowerApproximation = app; }

	/**
	* Return the shower approximation.
	*/
	Ptr<ShowerApproximation>::tptr showerApproximation() const { return theShowerApproximation; }

	/**
	* Set the shower approximation generator.
	*/
	void showerApproximationGenerator(Ptr<ShowerApproximationGenerator>::tptr);

	/**
	* Return the shower approximation generator.
	*/
	Ptr<ShowerApproximationGenerator>::tptr showerApproximationGenerator() const;

	/**
	* Generate the next emission
	*/
	double generate();

	public:

	/**
	* Set a pt cut on the dipole to generate the radiation
	*/
	void ptCut(Energy pt) { dipole()->ptCut(pt); }

	/**
	* Return the number of random numbers
	* needed to sample this kernel.
	*/
	int nDim() const {
	return
	nDimBorn() +
	dipole()->nDimRadiation();
	}

	/**
	* Return the number of random numbers
	* needed to sample the Born process.
	*/
	int nDimBorn() const {
	return bornCXComb()->lastRandomNumbers().size();
	}

	/**
	* Flag, which variables are free variables.
	*/
	const vector<bool>& sampleFlags();

	/**
	* Return the support of the splitting kernel.
	* The lower bound on the first variable is
	* assumed to correspond to the cutoff on the
	* evolution variable.
	*/
	const pair<vector<double>,vector<double> >& support();

	/**
	* Return the parameter point associated to the splitting
	* previously supplied through fixParameters.
	*/
	const vector<double>& parameterPoint();

	/**
	* Indicate that presampling of this kernel
	* will be performed in the next calls to
	* evaluate until stopPresampling() is called.
	*/
	void startPresampling();

	/**
	* Indicate that presampling of this kernel
	* is done until startPresampling() is called.
	*/
	void stopPresampling();

	/**
	+ * Indicate that a veto with the given kernel value and overestimate has occured.
	+ */
	+ void veto(const vector<double>&, double, double) {
	+ /** use born and real xcombs in here to figure out what we need to reweight;
	+ it should have its kinematic variables completed at this step */
	+ }
	+
	+ /**
	+ * Indicate that an accept with the given kernel value and overestimate has occured.
	+ */
	+ void accept(const vector<double>&, double, double) {
	+ /** use born and real xcombs in here to figure out what we need to reweight;
	+ it should have its kinematic variables completed at this step */
	+ }
	+
	+ /**
	* Return true, if currently being presampled
	*/
	bool presampling() const { return thePresampling; }

	/**
	* 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 points to presample this
	* splitting generator.
	*/
	void presamplingPoints(unsigned long p) { thePresamplingPoints = p; }

	/**
	* Set the maximum number of trials
	* to generate a splitting.
	*/
	void maxtry(unsigned long p) { theMaxTry = p; }

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

	/**
	* Evalute the splitting kernel.
	*/
	double evaluate(const vector<double>&);

	/**
	* Return the index of the random number corresponding
	* to the evolution variable.
	*/
	int evolutionVariable() const {
	return
	nDimBorn() +
	(showerApproximation()->showerInvertedTildeKinematics() ?
	showerApproximation()->showerInvertedTildeKinematics()->evolutionVariable() :
	dipole()->invertedTildeKinematics()->evolutionVariable());
	}

	/**
	* Return the cutoff on the evolution
	* random number corresponding to the pt cut.
	*/
	double evolutionCutoff() const {
	return
	showerApproximation()->showerInvertedTildeKinematics() ?
	showerApproximation()->showerInvertedTildeKinematics()->evolutionCutoff() :
	dipole()->invertedTildeKinematics()->evolutionCutoff();
	}

	/**
	* True, if sampler should apply compensation
	*/
	void doCompensate(bool yes = true) { theDoCompensate = yes; }

	public:

	/*@name Wrap to the exsample2 interface until this is finally cleaned up. /
	//@{

	inline const vector<bool>& variable_flags () {
	return sampleFlags();
	}

	inline size_t evolution_variable () const { return evolutionVariable(); }

	inline double evolution_cutoff () const {
	return evolutionCutoff();
	}

	inline const vector<double>& parameter_point () {
	return parameterPoint();
	}

	inline void start_presampling () {
	startPresampling();
	}

	inline void stop_presampling () {
	stopPresampling();
	}

	inline size_t dimension () const {
	return nDim();
	}

	inline unsigned long presampling_points () const {
	return presamplingPoints();
	}

	//@}

	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 dipole in charge of the emission
	*/
	Ptr<SubtractionDipole>::ptr theDipole;

	/**
	* The shower approximation to consider
	*/
	Ptr<ShowerApproximation>::ptr theShowerApproximation;

	/**
	* The XComb off which radiation will be generated
	*/
	StdXCombPtr theBornXComb;

	/**
	* The XComb describing the process after radiation
	*/
	StdXCombPtr theRealXComb;

	/**
	* The tilde xcomb objects associated to the real xcomb
	*/
	vector<StdXCombPtr> theTildeXCombs;

	/**
	* True, if currently being presampled
	*/
	bool thePresampling;

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

	/**
	* The maximum number of trials
	* to generate a splitting.
	*/
	unsigned long theMaxTry;

	/**
	* Return the number of accepted points after which the grid should
	* be frozen
	*/
	unsigned long theFreezeGrid;

	/**
	* The sampling flags
	*/
	vector<bool> theFlags;

	/**
	* The support.
	*/
	pair<vector<double>,vector<double> > theSupport;

	/**
	* The shower approximation generator.
	*/
	Ptr<ShowerApproximationGenerator>::tptr theShowerApproximationGenerator;

	/**
	* The last parameter point
	*/
	vector<double> theLastParameterPoint;

	/**
	* The last random numbers used for Born sampling
	*/
	vector<double> theLastBornPoint;

	/**
	* Define the Sudakov sampler
	*/
	typedef
	exsample::exponential_generator<ShowerApproximationKernel,UseRandom>
	ExponentialGenerator;

	/**
	* Define a pointer to the Sudakov sampler
	*/
	typedef
	exsample::exponential_generator<ShowerApproximationKernel,UseRandom>*
	ExponentialGeneratorPtr;

	/**
	* The Sudakov sampler
	*/
	ExponentialGeneratorPtr sampler;

	/**
	* True, if sampler should apply compensation
	*/
	bool theDoCompensate;

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

	};

	}

	#endif /* Herwig_ShowerApproximationKernel_H */
	diff --git a/Sampling/exsample/exponential_generator.icc b/Sampling/exsample/exponential_generator.icc
	--- a/Sampling/exsample/exponential_generator.icc
	+++ b/Sampling/exsample/exponential_generator.icc
	@@ -1,380 +1,383 @@
	// -- C++ --
	//
	// exponential_generator.icc is part of ExSample -- A Library for Sampling Sudakov-Type Distributions
	//
	// Copyright (C) 2008-2011 Simon Platzer -- simon.plaetzer@desy.de
	//
	// ExSample is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	namespace exsample {

	template<class Function, class Random>
	void exponential_generator<Function,Random>::initialize() {
	adaption_info_.dimension = function_->dimension();
	adaption_info_.lower_left = function_->support().first;
	adaption_info_.upper_right = function_->support().second;
	if (adaption_info_.adapt.empty())
	adaption_info_.adapt = std::vector<bool>(adaption_info_.dimension,true);
	evolution_variable_ = function_->evolution_variable();
	evolution_cutoff_ = function_->evolution_cutoff();
	sample_variables_ = function_->variable_flags();
	sample_other_variables_ = sample_variables_;
	sample_other_variables_[evolution_variable_] = false;
	last_point_.resize(adaption_info_.dimension);
	parametric_selector_ = parametric_selector(&last_point_,sample_other_variables_);
	exponent_selector_ = parametric_selector(&last_point_,sample_variables_);
	missing_accessor_ = parametric_missing_accessor(&last_parameter_bin_);
	parametric_sampler_ = parametric_sampling_selector<rnd_generator<Random> >
	(&last_point_,&last_parameter_bin_,sample_other_variables_,rnd_gen_);
	if (initialized_) return;
	splits_ = 0;
	for ( std::size_t k = 0; k < adaption_info_.dimension; ++k ) {
	if ( sample_other_variables_[k] )
	continue;
	parameter_splits_[k].push_back(adaption_info_.lower_left[k]);
	parameter_splits_[k].push_back(adaption_info_.upper_right[k]);
	}
	root_cell_ =
	binary_tree<cell>(cell(adaption_info_.lower_left,
	adaption_info_.upper_right,
	sample_other_variables_,
	adaption_info_));
	root_cell_.value().info().explore(rnd_gen_,adaption_info_,function_);
	root_cell_.value().integral(root_cell_.value().info().volume() * root_cell_.value().info().overestimate());
	last_exponent_integrand_.resize(1);
	check_events_ = adaption_info_.presampling_points;
	initialized_ = true;
	}

	template<class Function, class Random>
	bool exponential_generator<Function,Random>::split () {
	if (adaption_info_.freeze_grid <= accepts_)
	return false;
	if (compensating_)
	return false;
	if (!(*last_cell_).info().bad(adaption_info_)) return false;
	bool dosplit = false;
	std::pair<std::size_t,double> sp =
	(*last_cell_).info().get_split(adaption_info_,dosplit);
	if (!dosplit) return false;
	if (!adaption_info_.adapt[sp.first]) return false;
	if (splits_ == parameter_hash_bits/2)
	return false;
	++splits_;
	last_cell_.node().split((*last_cell_).split(sp,rnd_gen_,function_,adaption_info_,sample_other_variables_));
	if ( !sample_other_variables_[sp.first] ) {
	if ( std::find(parameter_splits_[sp.first].begin(),parameter_splits_[sp.first].end(),sp.second)
	== parameter_splits_[sp.first].end() ) {
	parameter_splits_[sp.first].push_back(sp.second);
	std::sort(parameter_splits_[sp.first].begin(),parameter_splits_[sp.first].end());
	if ( sp.first == evolution_variable_ ) {
	last_exponent_integrand_.push_back(0.);
	}
	}
	}
	did_split_ = true;
	last_point_ = function_->parameter_point();
	root_cell_.tree_accumulate(parametric_selector_,integral_accessor_,std::plus<double>());
	exponents_.clear();
	get_exponent();
	return true;
	}

	template<class Function, class Random>
	void exponential_generator<Function,Random>::get_exponent () {
	last_parameter_bin_.reset();
	root_cell_.subtree_hash (exponent_selector_,last_parameter_bin_);
	last_exponent_ = exponents_.find(last_parameter_bin_);
	if (last_exponent_ != exponents_.end())
	return;
	exponents_[last_parameter_bin_] = linear_interpolator();
	last_exponent_ = exponents_.find(last_parameter_bin_);
	double old_evo = last_point_[evolution_variable_];
	std::vector<double>::iterator exp_it = last_exponent_integrand_.begin();
	for (std::vector<double>::iterator esp = parameter_splits_[evolution_variable_].begin();
	esp < boost::prior(parameter_splits_[evolution_variable_].end()); ++esp, ++exp_it) {
	last_point_[evolution_variable_] = (esp + boost::next(esp))/2.;
	*exp_it = root_cell_.accumulate(parametric_selector_,integral_accessor_,std::plus<double>());
	}
	exp_it = boost::prior(last_exponent_integrand_.end());
	double total = 0.;
	for (std::vector<double>::iterator esp = boost::prior(parameter_splits_[evolution_variable_].end());
	esp > parameter_splits_[evolution_variable_].begin(); --esp, --exp_it) {
	last_exponent_->second.set_interpolation(*esp,total);
	total += (exp_it) ((esp) - (boost::prior(esp)));
	}
	last_exponent_->second.set_interpolation(parameter_splits_[evolution_variable_].front(),total);
	last_point_[evolution_variable_] = old_evo;
	}

	template<class Function, class Random>
	std::set<std::vector<double> >
	exponential_generator<Function,Random>::parameter_points() {
	std::set<std::vector<double> > res;
	std::vector<double> pt(adaption_info_.dimension,0.);
	recursive_parameter_points(res,pt,0);
	return res;
	}

	template<class Function, class Random>
	void exponential_generator<Function,Random>::
	recursive_parameter_points(std::set<std::vector<double> >& res,
	std::vector<double>& pt,
	size_t current) {
	if ( current == adaption_info_.dimension ) {
	res.insert(pt);
	return;
	}
	if ( sample_variables_[current] ) {
	recursive_parameter_points(res,pt,current+1);
	return;
	}
	for ( std::vector<double>::const_iterator sp =
	parameter_splits_[current].begin();
	sp != boost::prior(parameter_splits_[current].end()); ++sp ) {
	pt[current] = (sp + boost::next(sp))/2.;
	recursive_parameter_points(res,pt,current+1);
	}
	}

	template<class Function, class Random>
	void exponential_generator<Function,Random>::compensate() {
	if (!did_split_ \|\| !docompensate_) {
	assert(did_split_ \|\| last_cell_ == root_cell_.begin());
	exponents_.clear();
	last_cell_->info().overestimate(last_value_,last_point_);
	last_cell_->integral(last_cell_->info().volume() * last_cell_->info().overestimate());
	last_point_ = function_->parameter_point();
	get_exponent();
	return;
	}
	std::vector<double> themaxpoint = last_point_;
	std::set<std::vector<double> > id_points
	= parameter_points();
	for ( std::set<std::vector<double> >::const_iterator id =
	id_points.begin(); id != id_points.end(); ++id ) {
	last_point_ = *id;
	get_exponent();
	}
	std::map<bit_container<parameter_hash_bits>,linear_interpolator >
	old_exponents = exponents_;
	double old_oe = last_cell_->info().overestimate();
	last_cell_->info().overestimate(last_value_,themaxpoint);
	last_cell_->integral(last_cell_->info().volume() * last_cell_->info().overestimate());
	exponents_.clear();
	for ( std::set<std::vector<double> >::const_iterator id =
	id_points.begin(); id != id_points.end(); ++id ) {
	last_point_ = *id;
	get_exponent();
	std::map<bit_container<parameter_hash_bits>,linear_interpolator >::iterator
	old_exp = old_exponents.find(last_parameter_bin_);
	std::map<bit_container<parameter_hash_bits>,linear_interpolator >::iterator
	new_exp = exponents_.find(last_parameter_bin_);
	assert(old_exp != old_exponents.end() && new_exp != exponents_.end());
	double old_norm = 1. - std::exp(-(old_exp->second)(adaption_info_.lower_left[evolution_variable_]));
	double new_norm = 1. - std::exp(-(new_exp->second)(adaption_info_.lower_left[evolution_variable_]));
	for (binary_tree<cell>::iterator it = root_cell_.begin();
	it != root_cell_.end(); ++it) {
	if ( !it->info().contains_parameter(last_point_,sample_variables_) )
	continue;
	double old_int = 0.;
	double new_int = 0.;
	for ( std::vector<double>::const_iterator sp = parameter_splits_[evolution_variable_].begin();
	sp != boost::prior(parameter_splits_[evolution_variable_].end()); ++sp ) {
	if ( *sp >= it->info().lower_left()[evolution_variable_] &&
	*sp < it->info().upper_right()[evolution_variable_] ) {
	double xl = *sp;
	double xxl = *boost::next(sp);
	double old_al =
	(old_exp->second.interpolation()[xxl] - old_exp->second.interpolation()[xl]) /
	(xxl-xl);
	double old_bl =
	(xxl * old_exp->second.interpolation()[xl] -
	xl * old_exp->second.interpolation()[xxl]) /
	(xxl-xl);
	double new_al =
	(new_exp->second.interpolation()[xxl] - new_exp->second.interpolation()[xl]) /
	(xxl-xl);
	double new_bl =
	(xxl * new_exp->second.interpolation()[xl] -
	xl * new_exp->second.interpolation()[xxl]) /
	(xxl-xl);
	if ( std::abs(old_al) > std::numeric_limits<double>::epsilon() ) {
	old_int += (exp(-(old_alxl+old_bl)) - exp(-(old_alxxl+old_bl)))/old_al;
	} else {
	old_int += (xxl-xl)*exp(-old_bl);
	}
	if ( std::abs(new_al) > std::numeric_limits<double>::epsilon() ) {
	new_int += (exp(-(new_alxl+new_bl)) - exp(-(new_alxxl+new_bl)))/new_al;
	} else {
	new_int += (xxl-xl)*exp(-new_bl);
	}
	}
	}
	double scaling;
	if (it != last_cell_) {
	if (old_int > std::numeric_limits<double>::epsilon() &&
	new_int > std::numeric_limits<double>::epsilon())
	scaling = ((old_norm * new_int) /
	(new_norm * old_int)) - 1.;
	else
	scaling = 0.;
	} else {
	if (old_int > std::numeric_limits<double>::epsilon() &&
	new_int > std::numeric_limits<double>::epsilon())
	scaling = ((last_value_ * old_norm * new_int) /
	(old_oe * new_norm * old_int)) - 1.;
	else
	scaling = 0.;
	}
	it->info().parametric_missing(last_parameter_bin_,
	it->info().parametric_missing(last_parameter_bin_) +
	static_cast<int>(round(scaling * it->info().attempted())));
	if (it->info().parametric_missing(last_parameter_bin_) != 0) {
	compensating_ = true;
	}
	}
	}
	last_point_ = function_->parameter_point();
	}

	template<class Function, class Random>
	double exponential_generator<Function,Random>::generate () {
	if (compensating_) {
	compensating_ = false;
	for (binary_tree<cell>::iterator it = root_cell_.begin();
	it != root_cell_.end(); ++it)
	if (it->info().parametric_compensating()) {
	compensating_ = true;
	break;
	}
	parametric_sampler_.compensate(compensating_);
	}
	last_point_ = function_->parameter_point();
	if (last_point_[evolution_variable_] < evolution_cutoff_) {
	return 0.;
	}
	unsigned long n_hit_miss = 0;
	unsigned long n_select = 0;
	double minus_log_r;
	root_cell_.tree_accumulate(parametric_selector_,integral_accessor_,std::plus<double>());
	get_exponent();
	while (true) {
	n_select = 0;
	minus_log_r = -std::log(rnd_gen_()) + last_exponent_->second(last_point_[evolution_variable_]);
	if (!last_exponent_->second.invertible(minus_log_r)) {
	return 0.;
	}
	try {
	last_point_[evolution_variable_] = last_exponent_->second.unique_inverse(minus_log_r);
	} catch (constant_interpolation& c) {
	last_point_[evolution_variable_] = rnd_gen_(c.range.first,c.range.second);
	}
	assert(!std::isnan(last_point_[evolution_variable_]) &&
	!std::isinf(last_point_[evolution_variable_]));
	if (last_point_[evolution_variable_] < evolution_cutoff_) {
	return 0.;
	}
	++attempts_;
	if (compensating_) {
	root_cell_.tree_accumulate(missing_accessor_,std::plus<int>());
	}
	if (parameter_splits_[evolution_variable_].size() > 2)
	root_cell_.tree_accumulate(parametric_selector_,integral_accessor_,std::plus<double>());
	if (did_split_)
	while ((last_cell_ = root_cell_.select(parametric_sampler_)) == root_cell_.end()) {
	root_cell_.tree_accumulate(missing_accessor_,std::plus<int>());
	if(++n_select > adaption_info_.maxtry)
	throw selection_maxtry();
	}
	else
	last_cell_ = root_cell_.begin();
	last_cell_->info().select(rnd_gen_,last_point_,sample_other_variables_);
	last_value_ = function_->evaluate(last_point_);
	assert(last_value_ >= 0.);
	last_cell_->info().selected(last_point_,last_value_,adaption_info_);
	if (last_value_ > last_cell_->info().overestimate()) {
	if ( std::abs(last_value_)/last_cell_->info().overestimate() > 2. ) {
	last_value_ =
	last_cell_->info().overestimate()*
	(1.+exp(2.*(2.-std::abs(last_value_)/last_cell_->info().overestimate())));
	}
	compensate();
	throw exponential_regenerate();
	}
	if (last_cell_->info().attempted() % check_events_ == 0) {
	if (split()) {
	throw exponential_regenerate();
	}
	}
	- if (last_value_/last_cell_->info().overestimate() > rnd_gen_())
	+ if (last_value_/last_cell_->info().overestimate() > rnd_gen_()) {
	+ function_->accept(last_point_,last_value_,last_cell_->info().overestimate());
	break;
	+ }
	+ function_->veto(last_point_,last_value_,last_cell_->info().overestimate());
	if(++n_hit_miss > adaption_info_.maxtry)
	throw hit_and_miss_maxtry();
	}
	if (last_value_ == 0.)
	return 0.;
	++accepts_;
	++check_events_;
	last_cell_->info().accept();
	return 1.;
	}


	template<class Function, class Random>
	template<class OStream>
	void exponential_generator<Function,Random>::put (OStream& os) const {
	os << check_events_; ostream_traits<OStream>::separator(os);
	adaption_info_.put(os);
	root_cell_.put(os);
	os << did_split_; ostream_traits<OStream>::separator(os);
	os << initialized_; ostream_traits<OStream>::separator(os);
	os << evolution_variable_; ostream_traits<OStream>::separator(os);
	os << evolution_cutoff_; ostream_traits<OStream>::separator(os);
	os << sample_variables_; ostream_traits<OStream>::separator(os);
	os << sample_other_variables_; ostream_traits<OStream>::separator(os);
	os << parameter_splits_; ostream_traits<OStream>::separator(os);
	// last_cell_ is selected new so we ignore it here
	os << last_point_; ostream_traits<OStream>::separator(os);
	os << last_value_; ostream_traits<OStream>::separator(os);
	last_parameter_bin_.put(os);
	os << exponents_.size(); ostream_traits<OStream>::separator(os);
	for ( std::map<bit_container<parameter_hash_bits>,linear_interpolator >::const_iterator
	ex = exponents_.begin(); ex != exponents_.end() ; ++ex ) {
	ex->first.put(os);
	ex->second.put(os);
	}
	os << last_exponent_integrand_; ostream_traits<OStream>::separator(os);
	os << compensating_; ostream_traits<OStream>::separator(os);
	os << attempts_; ostream_traits<OStream>::separator(os);
	os << accepts_; ostream_traits<OStream>::separator(os);
	os << splits_; ostream_traits<OStream>::separator(os);
	os << docompensate_; ostream_traits<OStream>::separator(os);
	}

	template<class Function, class Random>
	template<class IStream>
	void exponential_generator<Function,Random>::get (IStream& is) {
	is >> check_events_;
	adaption_info_.get(is);
	root_cell_.get(is);
	is >> did_split_ >> initialized_ >> evolution_variable_
	>> evolution_cutoff_ >> sample_variables_ >> sample_other_variables_
	>> parameter_splits_;
	// last_cell_ is selected new so we ignore it here
	is >> last_point_ >> last_value_;
	last_parameter_bin_.get(is);
	size_t dim; is >> dim;
	for ( size_t k = 0; k < dim ; ++k ) {
	bit_container<parameter_hash_bits> key;
	key.get(is);
	exponents_[key].get(is);
	}
	is >> last_exponent_integrand_;
	last_exponent_ = exponents_.find(last_parameter_bin_);
	is >> compensating_ >> attempts_ >> accepts_ >> splits_ >> docompensate_;
	}

	}

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 9:06 PM (1 d, 19 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3799248
Default Alt Text: (118 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions