No OneTemporary
Actions

Size

222 KB

Subscribers

None

View Options

	diff --git a/DipoleShower/Base/DipoleChainOrdering.cc b/DipoleShower/Base/DipoleChainOrdering.cc
	--- a/DipoleShower/Base/DipoleChainOrdering.cc
	+++ b/DipoleShower/Base/DipoleChainOrdering.cc
	@@ -1,146 +1,146 @@
	// -- C++ --
	//
	// DipoleChainOrdering.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 DipoleChainOrdering class.
	//

	#include "DipoleChainOrdering.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Switch.h"

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

	using namespace Herwig;

	DipoleChainOrdering::DipoleChainOrdering()
	: DipoleEvolutionOrdering(), virtualityOrdering(false) {}

	DipoleChainOrdering::~DipoleChainOrdering() {}

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

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

	Energy DipoleChainOrdering::hardScale(tPPtr emitter, tPPtr spectator,
	double emitterX, double spectatorX,
	const DipoleSplittingKernel& split,
	const DipoleIndex& index) const {

	Energy scale =
	split.splittingKinematics()->dipoleScale(emitter->momentum(),
	spectator->momentum());

	return
	virtualityOrdering ?
	- split.splittingKinematics()->QMax(scale,emitterX,spectatorX,index) :
	+ split.splittingKinematics()->QMax(scale,emitterX,spectatorX,index,split) :
	split.splittingKinematics()->ptMax(scale,emitterX,spectatorX,index,split);

	}


	void DipoleChainOrdering::setEvolutionScale(Energy scale,
	const DipoleSplittingInfo&,
	DipoleChain& chain,
	pair<list<Dipole>::iterator,list<Dipole>::iterator>) const {

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

	if ( dip->emitterScale(make_pair(true,false)) > scale )
	dip->emitterScale(make_pair(true,false),scale);

	if ( dip->emitterScale(make_pair(false,true)) > scale )
	dip->emitterScale(make_pair(false,true),scale);

	}

	}

	void DipoleChainOrdering::setEvolutionScale(Energy scale,
	const DipoleSplittingInfo&,
	DipoleChain& chain,
	list<Dipole>::iterator) const {

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

	if ( dip->emitterScale(make_pair(true,false)) > scale )
	dip->emitterScale(make_pair(true,false),scale);

	if ( dip->emitterScale(make_pair(false,true)) > scale )
	dip->emitterScale(make_pair(false,true),scale);

	}

	}

	Energy DipoleChainOrdering::evolutionScale(const DipoleSplittingInfo& split,
	const DipoleSplittingKernel& spkernel) const {
	return
	virtualityOrdering ?
	spkernel.splittingKinematics()->QFromPt(split.lastPt(),split) :
	split.lastPt();
	}

	Energy DipoleChainOrdering::maxPt(Energy scale,
	const DipoleSplittingInfo& split,
	const DipoleSplittingKernel& spkernel) const {
	return
	virtualityOrdering ?
	spkernel.splittingKinematics()->ptMax(scale,split.emitterX(),split.spectatorX(),split.index(),spkernel) :
	scale;
	}

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


	void DipoleChainOrdering::persistentOutput(PersistentOStream & os) const {
	os << virtualityOrdering;
	}

	void DipoleChainOrdering::persistentInput(PersistentIStream & is, int) {
	is >> virtualityOrdering;
	}

	ClassDescription<DipoleChainOrdering> DipoleChainOrdering::initDipoleChainOrdering;
	// Definition of the static class description member.

	void DipoleChainOrdering::Init() {

	static ClassDocumentation<DipoleChainOrdering> documentation
	("DipoleChainOrdering performs ordering on "
	"complete colour singlet dipole chains.");


	static Switch<DipoleChainOrdering,bool> interfaceVirtualityOrdering
	("Ordering",
	"[experimental] Switch between virtuality and pt ordering.",
	&DipoleChainOrdering::virtualityOrdering, false, false, false);
	static SwitchOption interfaceVirtualityOrderingPt
	(interfaceVirtualityOrdering,
	"Pt",
	"Perform pt ordering",
	false);
	static SwitchOption interfaceVirtualityOrderingVirtuality
	(interfaceVirtualityOrdering,
	"Virtuality",
	"Perform virtuality ordering",
	true);

	interfaceVirtualityOrdering.rank(-1);

	}

	diff --git a/DipoleShower/Base/DipoleSplittingGenerator.cc b/DipoleShower/Base/DipoleSplittingGenerator.cc
	--- a/DipoleShower/Base/DipoleSplittingGenerator.cc
	+++ b/DipoleShower/Base/DipoleSplittingGenerator.cc
	@@ -1,553 +1,548 @@
	// -- 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 "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) {
	if ( ShowerHandler::currentHandler() )
	setGenerator(ShowerHandler::currentHandler()->generator());
	}

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

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

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

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

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

	generatedSplitting = sp;

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

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

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

	presampledSplitting = generatedSplitting;

	prepared = true;

	parameters.resize(nDim());

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

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

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

	}

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

	assert(generator());

	assert(!presampling);
	assert(prepared);

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

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

	- Energy maxPossible =
	- generatedSplitting.splittingKinematics()->ptMax(sp.scale(),
	+ generatedSplitting.hardPt(sp.hardPt());
	+
	+ parameters[0] = splittingKinematics()->ptToRandom(generatedSplitting.hardPt(),
	+ sp.scale(),
	sp.emitterX(), sp.spectatorX(),
	generatedSplitting.index(),
	*splittingKernel());

	- if ( maxPossible >= sp.hardPt() )
	- generatedSplitting.hardPt(sp.hardPt());
	- else
	- generatedSplitting.hardPt(maxPossible);
	-
	- parameters[0] = splittingKinematics()->ptToRandom(generatedSplitting.hardPt(),
	- sp.scale(),generatedSplitting.index());
	-
	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) )
	+ 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) ) {
	+ if ( ! split.splittingKinematics()->generateSplitting(point[0],point[1],point[2],split,*splittingKernel()) ) {
	split.lastValue(0.);
	return 0.;
	}

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

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

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

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

	if ( kernel < 0. )
	return 0.;

	return split.lastValue();

	}

	void DipoleSplittingGenerator::doGenerate() {

	assert(!wrapping());

	double res = 0.;

	Energy startPt = generatedSplitting.hardPt();

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

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

	generatedSplitting.continuesEvolving();

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

	}

	}

	Energy DipoleSplittingGenerator::generate(const DipoleSplittingInfo& split) {

	fixParameters(split);

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

	doGenerate();

	return generatedSplitting.lastPt();

	}

	Energy DipoleSplittingGenerator::generateWrapped(DipoleSplittingInfo& split) {

	assert(!wrapping());

	DipoleSplittingInfo backup = generatedSplitting;
	generatedSplitting = split;

	fixParameters(split);

	try {
	doGenerate();
	} 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;
	}

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

	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/DipoleShowerHandler.cc b/DipoleShower/DipoleShowerHandler.cc
	--- a/DipoleShower/DipoleShowerHandler.cc
	+++ b/DipoleShower/DipoleShowerHandler.cc
	@@ -1,992 +1,1008 @@
	// -- 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 "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;

	DipoleShowerHandler::DipoleShowerHandler()
	: ShowerHandler(), chainOrderVetoScales(true),
	nEmissions(0), discardNoEmissions(false), firstMCatNLOEmission(false),
	doFSR(true), doISR(true), realignmentScheme(0),
	hardFirstEmission(false),
	verbosity(0), printEvent(0), nTries(0),
	didRadiate(false), didRealign(false),
	theRenormalizationScaleFreeze(1.*GeV),
	theFactorizationScaleFreeze(1.*GeV),
	theFactorizationScaleFactor(1.0),
	theRenormalizationScaleFactor(1.0),
	theHardScaleFactor(1.0),
	isMCatNLOSEvent(false),
	isMCatNLOHEvent(false) {}

	DipoleShowerHandler::~DipoleShowerHandler() {}

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

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

	tPPair DipoleShowerHandler::cascade(tSubProPtr sub, XCPtr) {

	prepareCascade(sub);

	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;

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

	if ( subme ) {
	if ( subme->showerApproximation() ) {
	theShowerApproximation = subme->showerApproximation();
	if ( subme->realShowerSubtraction() )
	isMCatNLOHEvent = true;
	else if ( subme->virtualShowerSubtraction() )
	isMCatNLOSEvent = true;
	}
	} else if ( me ) {
	if ( me->factory()->showerApproximation() ) {
	theShowerApproximation = me->factory()->showerApproximation();
	isMCatNLOSEvent = true;
	}
	}

	if ( theShowerApproximation ) {
	theHardScaleFactor = theShowerApproximation->hardScaleFactor();
	}

	hardScales();

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

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

	} else {

	if ( nEmissions == 1 )
	doCascade(nEmitted);

	}

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

	}

	return eventRecord().fillEventRecord(newStep(),firstInteraction(),didRealign);

	}

	void DipoleShowerHandler::constituentReshuffle() {

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

	}

	void DipoleShowerHandler::hardScales() {

	Energy maxPt = generator()->maximumCMEnergy();

	bool restrictPhasespace = !hardFirstEmission;

	if ( !restrictPhasespace ) {
	assert(eventRecord().xcombPtr());
	Ptr<SubtractedME>::tptr subme =
	dynamic_ptr_cast<Ptr<SubtractedME>::tptr>(eventRecord().xcombPtr()->matrixElement());
	if ( subme )
	if ( subme->realShowerSubtraction() )
	restrictPhasespace = true;
	}

	if ( (eventRecord().incoming().first->coloured() \|\|
	eventRecord().incoming().second->coloured()) &&
	restrictPhasespace ) {
	if ( !eventRecord().outgoing().empty() ) {
	for ( PList::const_iterator p = eventRecord().outgoing().begin();
	p != eventRecord().outgoing().end(); ++p )
	maxPt = min(maxPt,(**p).momentum().perp());
	} 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 *= theHardScaleFactor;
	}

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

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

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

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

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

	tPPtr emitter = dip.emitter(conf);
	tPPtr spectator = dip.spectator(conf);
	Energy startScale = dip.emitterScale(conf);
	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()));
	- candidate.hardPt(hardScale);
	+ Energy maxPossible =
	+ gen->second->splittingKinematics()->ptMax(candidate.scale(),
	+ candidate.emitterX(), candidate.spectatorX(),
	+ candidate.index(),
	+ *gen->second->splittingKernel());
	+
	+ if ( maxPossible <= gen->second->splittingKinematics()->IRCutoff() ) {
	+ continue;
	+ }
	+
	+ if ( maxPossible >= hardScale )
	+ candidate.hardPt(hardScale);
	+ else
	+ candidate.hardPt(maxPossible);

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

	if ( isMCatNLOSEvent ) {
	- assert(theShowerApproximation);
	- if ( theShowerApproximation->restrictPhasespace() &&
	- theShowerApproximation->profileScales() ) {
	- while ( UseRandom::rnd() > theShowerApproximation->hardScaleProfile(hardScale,nextScale) ) {
	- candidate.continuesEvolving();
	- Energy nextHardScale = evolutionOrdering()->maxPt(nextScale,candidate,*(gen->second->splittingKernel()));
	- candidate.hardPt(nextHardScale);
	- gen->second->generate(candidate);
	- nextScale = evolutionOrdering()->evolutionScale(gen->second->lastSplitting(),*(gen->second->splittingKernel()));
	- if ( nextScale == ZERO \|\| candidate.stoppedEvolving() )
	- break;
	+ if ( eventRecord().incoming().first->coloured() \|\|
	+ eventRecord().incoming().second->coloured() ) {
	+ assert(theShowerApproximation);
	+ if ( theShowerApproximation->restrictPhasespace() &&
	+ theShowerApproximation->profileScales() ) {
	+ while ( UseRandom::rnd() > theShowerApproximation->hardScaleProfile(hardScale,nextScale) ) {
	+ candidate.continuesEvolving();
	+ Energy nextHardScale = evolutionOrdering()->maxPt(nextScale,candidate,*(gen->second->splittingKernel()));
	+ candidate.hardPt(nextHardScale);
	+ gen->second->generate(candidate);
	+ nextScale = evolutionOrdering()->evolutionScale(gen->second->lastSplitting(),*(gen->second->splittingKernel()));
	+ if ( nextScale == ZERO \|\| candidate.stoppedEvolving() )
	+ break;
	+ }
	}
	}
	}

	if ( nextScale > winnerScale ) {
	winner = 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) {

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

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

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

	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() ) {
	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 ( 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->splittingKernel(*k);
	nGenerator->splittingKernel()->renormalizationScaleFactor(theRenormalizationScaleFactor);
	nGenerator->splittingKernel()->factorizationScaleFactor(theFactorizationScaleFactor);

	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 << hardFirstEmission << verbosity << printEvent
	<< ounit(theRenormalizationScaleFreeze,GeV)
	<< ounit(theFactorizationScaleFreeze,GeV)
	<< theFactorizationScaleFactor << theRenormalizationScaleFactor
	<< theHardScaleFactor
	<< isMCatNLOSEvent << isMCatNLOHEvent << theShowerApproximation;
	}

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

	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,bool> interfaceHardFirstEmission
	("HardFirstEmission",
	"Switch on or off hard first emission.",
	&DipoleShowerHandler::hardFirstEmission, false, false, false);
	static SwitchOption interfaceHardFirstEmissionOn
	(interfaceHardFirstEmission,
	"On",
	"Switch on hard first emission.",
	true);
	static SwitchOption interfaceHardFirstEmissionOff
	(interfaceHardFirstEmission,
	"Off",
	"Switch off hard first emission.",
	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,double> interfaceFactorizationScaleFactor
	("FactorizationScaleFactor",
	"The factorization scale factor.",
	&DipoleShowerHandler::theFactorizationScaleFactor, 1.0, 0.0, 0,
	false, false, Interface::lowerlim);

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

	static Parameter<DipoleShowerHandler,double> interfaceHardScaleFactor
	("HardScaleFactor",
	"The hard scale factor.",
	&DipoleShowerHandler::theHardScaleFactor, 1.0, 0.0, 0,
	false, false, Interface::lowerlim);

	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, 1.0GeV, 0.0GeV, 0*GeV,
	false, false, Interface::lowerlim);

	}

	diff --git a/DipoleShower/Kinematics/DipoleSplittingKinematics.cc b/DipoleShower/Kinematics/DipoleSplittingKinematics.cc
	--- a/DipoleShower/Kinematics/DipoleSplittingKinematics.cc
	+++ b/DipoleShower/Kinematics/DipoleSplittingKinematics.cc
	@@ -1,167 +1,238 @@
	// -- C++ --
	//
	// DipoleSplittingKinematics.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.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the DipoleSplittingKinematics class.
	//

	#include "DipoleSplittingKinematics.h"
	#include "Herwig++/DipoleShower/Base/DipoleSplittingInfo.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Interface/Parameter.h"

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

	+#include "Herwig++/MatrixElement/Matchbox/Phasespace/RandomHelpers.h"
	+
	#include <limits>

	using namespace Herwig;

	DipoleSplittingKinematics::DipoleSplittingKinematics()
	: HandlerBase(), theIRCutoff(1.0*GeV),
	theXMin(1.e-5), theJacobian(0.0),
	theLastPt(0.0*GeV), theLastZ(0.0), theLastPhi(0.0),
	theLastEmitterZ(1.0), theLastSpectatorZ(1.0),
	theLastSplittingParameters() {}

	DipoleSplittingKinematics::~DipoleSplittingKinematics() {}

	void DipoleSplittingKinematics::persistentOutput(PersistentOStream & os) const {
	os << ounit(theIRCutoff,GeV) << theXMin << theMCCheck;
	}

	void DipoleSplittingKinematics::persistentInput(PersistentIStream & is, int) {
	is >> iunit(theIRCutoff,GeV) >> theXMin >> theMCCheck;
	}

	void DipoleSplittingKinematics::prepareSplitting(DipoleSplittingInfo& dInfo) {

	dInfo.splittingKinematics(this);

	if ( lastPt() > IRCutoff() )
	dInfo.lastPt(lastPt());
	else {
	dInfo.lastPt(0.0*GeV);
	dInfo.didStopEvolving();
	}

	dInfo.lastZ(lastZ());
	dInfo.lastPhi(lastPhi());
	dInfo.lastEmitterZ(lastEmitterZ());
	dInfo.lastSpectatorZ(lastSpectatorZ());
	dInfo.splittingParameters().resize(lastSplittingParameters().size());
	copy(lastSplittingParameters().begin(),lastSplittingParameters().end(),
	dInfo.splittingParameters().begin());

	}

	+Energy DipoleSplittingKinematics::generatePt(double r, Energy dScale,
	+ double emX, double specX,
	+ const DipoleIndex& dIndex,
	+ const DipoleSplittingKernel& split,
	+ double& weight) const {
	+
	+ Energy maxPt = ptMax(dScale,emX,specX,dIndex,split);
	+ if ( maxPt <= IRCutoff() ) {
	+ weight = 0.0;
	+ return ZERO;
	+ }
	+
	+ weight *= log(sqr(maxPt/IRCutoff()));
	+
	+ return IRCutoff()*pow(maxPt/IRCutoff(),r);
	+
	+}
	+
	+double DipoleSplittingKinematics::ptToRandom(Energy pt, Energy dScale,
	+ double emX, double specX,
	+ const DipoleIndex& dIndex,
	+ const DipoleSplittingKernel& split) const {
	+
	+ Energy maxPt = ptMax(dScale,emX,specX,dIndex,split);
	+ assert(pt >= IRCutoff() && pt <= maxPt);
	+
	+ return log(pt/IRCutoff())/log(maxPt/IRCutoff());
	+
	+}
	+
	+double DipoleSplittingKinematics::generateZ(double r, Energy pt, int sampling,
	+ const DipoleSplittingInfo& dInfo,
	+ const DipoleSplittingKernel& split,
	+ double& weight) const {
	+
	+ pair<double,double> zLims = zBoundaries(pt,dInfo,split);
	+
	+ using namespace RandomHelpers;
	+
	+ if ( sampling == FlatZ ) {
	+ pair<double,double> kw = generate(flat(zLims.first,zLims.second),r);
	+ weight *= kw.second;
	+ return kw.first;
	+ }
	+
	+ if ( sampling == OneOverZ ) {
	+ pair<double,double> kw = generate(inverse(0.0,zLims.first,zLims.second),r);
	+ weight *= kw.second;
	+ return kw.first;
	+ }
	+
	+ if ( sampling == OneOverOneMinusZ ) {
	+ pair<double,double> kw = generate(inverse(1.0,zLims.first,zLims.second),r);
	+ weight *= kw.second;
	+ return kw.first;
	+ }
	+
	+ if ( sampling == OneOverOneMinusZ ) {
	+ pair<double,double> kw = generate(inverse(0.0,zLims.first,zLims.second) +
	+ inverse(1.0,zLims.first,zLims.second),r);
	+ weight *= kw.second;
	+ return kw.first;
	+ }
	+
	+ weight = 0.0;
	+ return 0.0;
	+
	+}
	+
	Lorentz5Momentum DipoleSplittingKinematics::getKt(const Lorentz5Momentum& p1,
	const Lorentz5Momentum& p2,
	Energy pt,
	double phi,
	bool spacelike) const {

	Lorentz5Momentum P;
	if ( !spacelike )
	P = p1 + p2;
	else
	P = p1 - p2;

	Energy2 Q2 = abs(P.m2());

	Lorentz5Momentum Q =
	!spacelike ?
	Lorentz5Momentum(ZERO,ZERO,ZERO,sqrt(Q2),sqrt(Q2)) :
	Lorentz5Momentum(ZERO,ZERO,sqrt(Q2),ZERO,-sqrt(Q2));

	if ( spacelike && Q.z() < P.z() )
	Q.setZ(-Q.z());

	bool boost =
	abs((P-Q).vect().mag2()/GeV2) > 1e-10 \|\|
	abs((P-Q).t()/GeV) > 1e-5;

	Lorentz5Momentum inFrame1;
	if ( boost )
	inFrame1 = p1 + ((Pp1-Qp1)/(PQ-Q.mass2()))(P-Q);
	else
	inFrame1 = p1;

	Energy ptx = inFrame1.x();
	Energy pty = inFrame1.y();
	Energy q = 2.*inFrame1.z();

	Energy Qp = sqrt(4.*(sqr(ptx)+sqr(pty))+sqr(q));
	Energy Qy = sqrt(4.*sqr(pty)+sqr(q));

	double cPhi = cos(phi);
	double sPhi = sqrt(1.-sqr(cPhi));
	if ( phi > Constants::pi )
	sPhi = -sPhi;

	Lorentz5Momentum kt;

	if ( !spacelike ) {
	kt.setT(ZERO);
	kt.setX(ptQycPhi/Qp);
	kt.setY(-pt(4ptxptycPhi/Qp+q*sPhi)/Qy);
	kt.setZ(2.pt(-ptxqcPhi/Qp + pty*sPhi)/Qy);
	} else {
	kt.setT(2.pt(ptxqcPhi+ptyQpsPhi)/(q*Qy));
	kt.setX(pt(QpqcPhi+4.ptxptysPhi)/(q*Qy));
	kt.setY(ptQysPhi/q);
	kt.setZ(ZERO);
	}

	if ( boost )
	kt = kt + ((Pkt-Qkt)/(PQ-Q.mass2()))(P-Q);
	kt.setMass(-pt);
	kt.rescaleRho();

	return kt;

	}

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

	AbstractClassDescription<DipoleSplittingKinematics> DipoleSplittingKinematics::initDipoleSplittingKinematics;
	// Definition of the static class description member.

	void DipoleSplittingKinematics::Init() {

	static ClassDocumentation<DipoleSplittingKinematics> documentation
	("DipoleSplittingKinematics is the base class for dipole splittings "
	"as performed in the dipole shower.");


	static Parameter<DipoleSplittingKinematics,Energy> interfaceIRCutoff
	("IRCutoff",
	"The IR cutoff to be used by this splitting kinematics.",
	&DipoleSplittingKinematics::theIRCutoff, GeV, 1.0GeV, 0.0GeV, 0*GeV,
	false, false, Interface::lowerlim);


	static Parameter<DipoleSplittingKinematics,double> interfaceXMin
	("XMin",
	"The minimum momentum fraction for incoming partons",
	&DipoleSplittingKinematics::theXMin, 1.0e-5, 0.0, 1.0,
	false, false, Interface::limited);


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

	interfaceMCCheck.rank(-1);

	}

	diff --git a/DipoleShower/Kinematics/DipoleSplittingKinematics.h b/DipoleShower/Kinematics/DipoleSplittingKinematics.h
	--- a/DipoleShower/Kinematics/DipoleSplittingKinematics.h
	+++ b/DipoleShower/Kinematics/DipoleSplittingKinematics.h
	@@ -1,476 +1,507 @@
	// -- C++ --
	//
	// DipoleSplittingKinematics.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_DipoleSplittingKinematics_H
	#define HERWIG_DipoleSplittingKinematics_H
	//
	// This is the declaration of the DipoleSplittingKinematics class.
	//

	#include "ThePEG/Handlers/HandlerBase.h"
	#include "ThePEG/Vectors/Lorentz5Vector.h"

	#include "Herwig++/DipoleShower/Utility/DipoleMCCheck.h"

	namespace Herwig {

	using namespace ThePEG;

	class DipoleIndex;
	class DipoleSplittingInfo;
	class DipoleSplittingKernel;

	/**
	* \ingroup DipoleShower
	* \author Simon Platzer
	*
	* \brief DipoleSplittingKinematics is the base class for dipole splittings
	* as performed in the dipole shower.
	*
	* @see \ref DipoleSplittingKinematicsInterfaces "The interfaces"
	* defined for DipoleSplittingKinematics.
	*/
	class DipoleSplittingKinematics: public HandlerBase {

	public:

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

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

	public:

	/**
	* Return the boundaries in between the evolution
	* variable random number is to be sampled; the lower
	* cuoff is assumed to correspond to the infrared cutoff.
	*/
	- virtual pair<double,double> kappaSupport(const DipoleSplittingInfo& dIndex) const = 0;
	+ virtual pair<double,double> kappaSupport(const DipoleSplittingInfo&) const {
	+ return make_pair(0.0,1.0);
	+ }

	/**
	* Return the boundaries in between the momentum
	* fraction random number is to be sampled.
	*/
	- virtual pair<double,double> xiSupport(const DipoleSplittingInfo& dIndex) const = 0;
	+ virtual pair<double,double> xiSupport(const DipoleSplittingInfo&) const {
	+ return make_pair(0.0,1.0);
	+ }

	/**
	* Return the dipole scale associated to the
	* given pair of emitter and spectator. This
	* should be the invariant mass or absolute value
	* final/final or initial/initial and the absolute
	* value of the momentum transfer for intial/final or
	* final/initial dipoles.
	*/
	virtual Energy dipoleScale(const Lorentz5Momentum& pEmitter,
	- const Lorentz5Momentum& pSpectator) const = 0;
	+ const Lorentz5Momentum& pSpectator) const {
	+ return sqrt(2.pEmitterpSpectator);
	+ }

	/**
	* Return the maximum pt for the given dipole scale.
	*/
	virtual Energy ptMax(Energy dScale,
	double emX, double specX,
	const DipoleIndex& dIndex,
	const DipoleSplittingKernel& split) const = 0;

	/**
	* Return the maximum virtuality for the given dipole scale.
	*/
	virtual Energy QMax(Energy dScale,
	double emX, double specX,
	- const DipoleIndex& dIndex) const = 0;
	+ const DipoleIndex& dIndex,
	+ const DipoleSplittingKernel& split) const = 0;

	/**
	* Return the pt given a virtuality.
	*/
	virtual Energy PtFromQ(Energy scale, const DipoleSplittingInfo&) const = 0;

	/**
	* Return the virtuality given a pt.
	*/
	virtual Energy QFromPt(Energy scale, const DipoleSplittingInfo&) const = 0;

	/**
	* Return the infrared cutoff.
	*/
	virtual Energy IRCutoff() const { return theIRCutoff; }

	/**
	* Return the minimum momentum fraction for
	* incoming partons
	*/
	double xMin() const { return theXMin; }

	/**
	+ * Generate a pt
	+ */
	+ Energy generatePt(double r, Energy dScale,
	+ double emX, double specX,
	+ const DipoleIndex& dIndex,
	+ const DipoleSplittingKernel& split,
	+ double& weight) const;
	+
	+ /**
	* Return the random number associated to
	* the given pt.
	*/
	virtual double ptToRandom(Energy pt, Energy dScale,
	- const DipoleIndex& dIndex) const = 0;
	+ double emX, double specX,
	+ const DipoleIndex& dIndex,
	+ const DipoleSplittingKernel& split) const;
	+
	+ /**
	+ * Return the boundaries on the momentum fraction
	+ */
	+ virtual pair<double,double> zBoundaries(Energy pt,
	+ const DipoleSplittingInfo& dInfo,
	+ const DipoleSplittingKernel& split) const = 0;
	+
	+ /**
	+ * Enumerate the variants of sampling z
	+ */
	+ enum ZSamplingOptions {
	+
	+ FlatZ = 0,
	+ OneOverZ,
	+ OneOverOneMinusZ,
	+ OneOverZOneMinusZ
	+
	+ };
	+
	+ /**
	+ * Generate a z value flat
	+ */
	+ double generateZ(double r, Energy pt, int sampling,
	+ const DipoleSplittingInfo& dInfo,
	+ const DipoleSplittingKernel& split,
	+ double& weight) const;

	/**
	* Generate splitting variables given three random numbers
	* and the momentum fractions of the emitter and spectator.
	* Return true on success.
	*/
	virtual bool generateSplitting(double kappa, double xi, double phi,
	- DipoleSplittingInfo& dIndex) = 0;
	-
	- /**
	- * For the splitting products present in the given dipole splitting
	- * info object calculate the kinematics parameters and return the
	- * propagator factor.
	- */
	- virtual InvEnergy2 setKinematics(DipoleSplittingInfo&) const = 0;
	-
	- /**
	- * For the splitting parameters given in the dipole splitting info
	- * object, calculate the phasespace Jacobian times the propagator
	- * factor.
	- */
	- virtual double jacobianTimesPropagator(const DipoleSplittingInfo&,
	- Energy) const = 0;
	+ DipoleSplittingInfo& dIndex,
	+ const DipoleSplittingKernel&) = 0;

	/**
	* Get the splitting phasespace weight associated to
	* the last call to generateSplitting. This is taken to
	* be the single particle phasespace times 16 \pi^2 divided
	* by the relevant propagator invariant.
	*/
	double jacobian() const { return theJacobian; }

	/**
	* Return true, if this splitting kinematics
	* class is capable of delivering an overestimate
	* to the jacobian.
	*/
	virtual bool haveOverestimate() const { return false; }

	/**
	* Return the overestimated jacobian for the
	* last generated parameters.
	*/
	virtual double jacobianOverestimate() const { return -1.; }

	/**
	* Return the last generated pt
	*/
	Energy lastPt() const { return theLastPt; }

	/**
	* Return the last generated momentum fraction.
	*/
	double lastZ() const { return theLastZ; }

	/**
	* Return the last generated azimuthal angle.
	*/
	double lastPhi() const { return theLastPhi; }

	/**
	* Return the momentum fraction, by which the emitter's
	* momentum fraction should be divided after the splitting.
	*/
	double lastEmitterZ() const { return theLastEmitterZ; }

	/**
	* Return the momentum fraction, by which the spectator's
	* momentum fraction should be divided after the splitting.
	*/
	double lastSpectatorZ() const { return theLastSpectatorZ; }

	/**
	* Return any additional parameters needed to
	* evaluate the splitting kernel or to generate the
	* full splitting.
	*/
	const vector<double>& lastSplittingParameters() const { return theLastSplittingParameters; }

	/**
	* Complete a DipoleSplittingInfo object with
	* the parameters generated by the last call to
	* generateSplitting()
	*/
	void prepareSplitting(DipoleSplittingInfo& dInfo);

	public:

	/**
	* Generate the full kinematics given emitter and
	* spectator momentum and a previously completeted
	* DipoleSplittingInfo object.
	*/
	virtual void generateKinematics(const Lorentz5Momentum& pEmitter,
	const Lorentz5Momentum& pSpectator,
	const DipoleSplittingInfo& dInfo) = 0;

	/**
	* Return the emitter's momentum after the splitting.
	*/
	const Lorentz5Momentum& lastEmitterMomentum() const { return theEmitterMomentum; }

	/**
	* Return the spectator's momentum after the splitting.
	*/
	const Lorentz5Momentum& lastSpectatorMomentum() const { return theSpectatorMomentum; }

	/**
	* Return the emission's momentum.
	*/
	const Lorentz5Momentum& lastEmissionMomentum() const { return theEmissionMomentum; }

	/*
	* Return true, if there is a transformation which should
	* be applied to all other final state particles except the ones
	* involved in the splitting after having performed the splitting
	*/
	virtual bool doesTransform () const { return false; }

	/*
	* perform the transformation, if existing
	*/
	virtual Lorentz5Momentum transform (const Lorentz5Momentum& p) const { return p; }

	protected:

	/**
	* Calculate a transverse momentum for the given momenta,
	* invariant pt and azimuth.
	*/
	Lorentz5Momentum getKt(const Lorentz5Momentum& p1,
	const Lorentz5Momentum& p2,
	Energy pt,
	double phi,
	bool spacelike = false) const;

	/**
	* Set the splitting phasespace weight associated to
	* the last call to generateSplitting. This is taken to
	* be the single particle phasespace times 16 \pi^2 divided
	* by the relevant propagator invariant.
	*/
	void jacobian(double w) { theJacobian = w; }

	/**
	* Set the last generated pt
	*/
	void lastPt(Energy p) { theLastPt = p; }

	/**
	* Set the last generated momentum fraction.
	*/
	void lastZ(double z) { theLastZ = z; }

	/**
	* Set the last generated azimuthal angle.
	*/
	void lastPhi(double p) { theLastPhi = p; }

	/**
	* Set the momentum fraction, by which the emitter's
	* momentum fraction should be divided after the splitting.
	*/
	void lastEmitterZ(double z) { theLastEmitterZ = z; }

	/**
	* Set the momentum fraction, by which the spectator's
	* momentum fraction should be divided after the splitting.
	*/
	void lastSpectatorZ(double z) { theLastSpectatorZ = z; }

	/**
	* Access any additional parameters needed to
	* evaluate the splitting kernel or to generate the
	* full splitting.
	*/
	vector<double>& splittingParameters() { return theLastSplittingParameters; }

	/**
	* Set the emitter's momentum after the splitting.
	*/
	void emitterMomentum(const Lorentz5Momentum& p) { theEmitterMomentum = p; }

	/**
	* Set the spectator's momentum after the splitting.
	*/
	void spectatorMomentum(const Lorentz5Momentum& p) { theSpectatorMomentum = p; }

	/**
	* Set the emission's momentum.
	*/
	void emissionMomentum(const Lorentz5Momentum& p) { theEmissionMomentum = p; }

	public:

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

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


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

	private:

	/**
	* The infrared cutoff associated to this
	* splitting kinematics.
	*/
	Energy theIRCutoff;

	/**
	* The minimum momentum fraction for
	* incoming partons
	*/
	double theXMin;

	/**
	* The last calculated splitting phase space weight.
	*/
	double theJacobian;

	/**
	* The last generated pt
	*/
	Energy theLastPt;

	/**
	* The last generated momentum fraction.
	*/
	double theLastZ;

	/**
	* The last generated azimuthal angle.
	*/
	double theLastPhi;

	/**
	* The momentum fraction, by which the emitter's
	* momentum fraction should be divided after the splitting.
	*/
	double theLastEmitterZ;

	/**
	* The momentum fraction, by which the spectator's
	* momentum fraction should be divided after the splitting.
	*/
	double theLastSpectatorZ;

	/**
	* Any additional parameters needed to
	* evaluate the splitting kernel or to generate the
	* full splitting.
	*/
	vector<double> theLastSplittingParameters;

	/**
	* The emitter's momentum after the splitting.
	*/
	Lorentz5Momentum theEmitterMomentum;

	/**
	* The emission's momentum after the splitting.
	*/
	Lorentz5Momentum theEmissionMomentum;

	/**
	* The spectator's momentum after the splitting.
	*/
	Lorentz5Momentum theSpectatorMomentum;

	protected:

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

	private:

	/**
	* The static object used to initialize the description of this class.
	* Indicates that this is an abstract class.
	*/
	static AbstractClassDescription<DipoleSplittingKinematics> initDipoleSplittingKinematics;

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

	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

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

	#include "FFLightKinematics.h"
	#include "ThePEG/Interface/ClassDocumentation.h"

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

	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "Herwig++/DipoleShower/Base/DipoleSplittingInfo.h"

	using namespace Herwig;

	FFLightKinematics::FFLightKinematics()
	: DipoleSplittingKinematics() {}

	FFLightKinematics::~FFLightKinematics() {}

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

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

	-pair<double,double> FFLightKinematics::kappaSupport(const DipoleSplittingInfo&) const {
	- return make_pair(0.0,1.0);
	-}
	-
	-pair<double,double> FFLightKinematics::xiSupport(const DipoleSplittingInfo& split) const {
	- double c = sqrt(1.-4.*sqr(IRCutoff()/generator()->maximumCMEnergy()));
	- if ( split.index().emitterData()->id() == ParticleID::g ) {
	- if ( split.emissionData()->id() != ParticleID::g )
	- return make_pair(0.5(1.-c),0.5(1.+c));
	- double b = log((1.+c)/(1.-c));
	- return make_pair(-b,b);
	- }
	- return make_pair(-log(0.5(1.+c)),-log(0.5(1.-c)));
	-}
	-
	-Energy FFLightKinematics::dipoleScale(const Lorentz5Momentum& pEmitter,
	- const Lorentz5Momentum& pSpectator) const {
	- return sqrt(2.(pEmitterpSpectator));
	-}
	-
	Energy FFLightKinematics::ptMax(Energy dScale,
	double, double,
	const DipoleIndex&,
	const DipoleSplittingKernel&) const {
	return dScale/2.;
	}

	Energy FFLightKinematics::QMax(Energy dScale,
	double, double,
	- const DipoleIndex&) const {
	+ const DipoleIndex&,
	+ const DipoleSplittingKernel&) const {
	return dScale;
	}

	Energy FFLightKinematics::PtFromQ(Energy scale, const DipoleSplittingInfo& split) const {
	double z = split.lastZ();
	return scalesqrt(z(1.-z));
	}

	Energy FFLightKinematics::QFromPt(Energy scale, const DipoleSplittingInfo& split) const {
	double z = split.lastZ();
	return scale/sqrt(z*(1.-z));
	}

	-double FFLightKinematics::ptToRandom(Energy pt, Energy,
	- const DipoleIndex&) const {
	- return log(pt/IRCutoff()) / log(0.5 * generator()->maximumCMEnergy()/IRCutoff());
	+pair<double,double> FFLightKinematics::zBoundaries(Energy pt,
	+ const DipoleSplittingInfo& dInfo,
	+ const DipoleSplittingKernel&) const {
	+ double s = sqrt(1.-sqr(pt/dInfo.hardPt()));
	+ return make_pair(0.5(1.-s),0.5(1.+s));
	}

	bool FFLightKinematics::generateSplitting(double kappa, double xi, double rphi,
	- DipoleSplittingInfo& info) {
	+ DipoleSplittingInfo& info,
	+ const DipoleSplittingKernel& split) {

	- Energy pt = IRCutoff() * pow(0.5 * generator()->maximumCMEnergy()/IRCutoff(),kappa);
	+ double weight = 1.0;

	- if ( pt > info.hardPt() ) {
	+ Energy pt = generatePt(kappa,info.scale(),
	+ info.emitterX(),info.spectatorX(),
	+ info.index(),split,
	+ weight);
	+
	+ if ( pt < IRCutoff() \|\| pt > info.hardPt() ) {
	jacobian(0.0);
	return false;
	}

	- double z;
	- double mapZJacobian;
	+ double z = 0.0;

	if ( info.index().emitterData()->id() == ParticleID::g ) {
	if ( info.emissionData()->id() != ParticleID::g ) {
	- z = xi;
	- mapZJacobian = 1.;
	+ z = generateZ(xi,pt,FlatZ,
	+ info,split,weight);
	} else {
	- z = exp(xi)/(1.+exp(xi));
	- mapZJacobian = z*(1.-z);
	+ z = generateZ(xi,pt,OneOverZOneMinusZ,
	+ info,split,weight);
	}
	} else {
	- z = 1.-exp(-xi);
	- mapZJacobian = 1.-z;
	+ z = generateZ(xi,pt,OneOverOneMinusZ,
	+ info,split,weight);
	}

	- double s = z*(1.-z);
	- double zp = 0.5*(1.+sqrt(1.-sqr(pt/info.hardPt())));
	- double zm = 0.5*(1.-sqrt(1.-sqr(pt/info.hardPt())));
	+ double y = sqr(pt/info.scale())/(z*(1.-z));

	- if ( pt < IRCutoff() \|\|
	- pt > info.hardPt() \|\|
	- z > zp \|\| z < zm ) {
	+ if ( z < 0.0 \|\| z > 1.0 \|\|
	+ y < 0.0 \|\| y > 1.0 ) {
	jacobian(0.0);
	return false;
	}

	double phi = 2.Constants::pirphi;

	- jacobian( 2. * mapZJacobian * (1. - sqr(pt) / (s * sqr(info.scale())) ) *
	- log(0.5 * generator()->maximumCMEnergy()/IRCutoff()) );
	+ jacobian(weight*(1.-y));

	lastPt(pt);
	lastZ(z);
	lastPhi(phi);

	if ( theMCCheck )
	theMCCheck->book(1.,1.,info.scale(),info.hardPt(),pt,z,jacobian());

	return true;

	}

	-InvEnergy2 FFLightKinematics::setKinematics(DipoleSplittingInfo& split) const {
	-
	- Lorentz5Momentum emitter = split.splitEmitter()->momentum();
	- Lorentz5Momentum emission = split.emission()->momentum();
	- Lorentz5Momentum spectator = split.splitSpectator()->momentum();
	-
	- split.splittingKinematics(const_cast<FFLightKinematics*>(this));
	-
	- Energy2 scale = 2.(emissionemitter + emissionspectator + emitterspectator);
	- split.scale(sqrt(scale));
	-
	- double y = 2.emissionemitter / scale;
	- double z = emitterspectator / (emitterspectator + emission*spectator);
	-
	- split.lastPt(split.scale() * sqrt(yz(1.-z)));
	- split.lastZ(z);
	-
	- split.hardPt(split.lastPt());
	-
	- if ( split.hardPt() > IRCutoff() ) {
	- split.continuesEvolving();
	- } else {
	- split.didStopEvolving();
	- }
	-
	- return 1./(2.(emitteremission));
	-
	-}
	-
	-double FFLightKinematics::
	-jacobianTimesPropagator(const DipoleSplittingInfo& split,
	- Energy scale) const {
	-
	- Energy pt = split.lastPt();
	- double z = split.lastZ();
	- double s = z*(1.-z);
	- double zp = 0.5*(1.+sqrt(1.-sqr(pt/split.hardPt())));
	- double zm = 0.5*(1.-sqrt(1.-sqr(pt/split.hardPt())));
	-
	- if ( pt < IRCutoff() \|\|
	- pt > split.hardPt() \|\|
	- z > zp \|\| z < zm ) {
	- return 0.;
	- }
	-
	- return (2.scale/pt)(1.-sqr(pt)/(s*sqr(scale)));
	-
	-}
	-
	-
	void FFLightKinematics::generateKinematics(const Lorentz5Momentum& pEmitter,
	const Lorentz5Momentum& pSpectator,
	const DipoleSplittingInfo& dInfo) {

	double z = dInfo.lastZ();
	Energy pt = dInfo.lastPt();
	double y = sqr(pt / (pEmitter+pSpectator).m()) / (z*(1.-z));

	Lorentz5Momentum kt =
	getKt(pEmitter, pSpectator, pt, dInfo.lastPhi());

	Lorentz5Momentum em = zpEmitter + y(1.-z)*pSpectator + kt;
	em.setMass(0.*GeV);
	em.rescaleEnergy();

	Lorentz5Momentum emm = (1.-z)pEmitter + zy*pSpectator - kt;
	emm.setMass(0.*GeV);
	emm.rescaleEnergy();

	Lorentz5Momentum spe = (1.-y)*pSpectator;
	spe.setMass(0.*GeV);
	spe.rescaleEnergy();

	emitterMomentum(em);
	emissionMomentum(emm);
	spectatorMomentum(spe);

	}

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


	void FFLightKinematics::persistentOutput(PersistentOStream & ) const {
	}

	void FFLightKinematics::persistentInput(PersistentIStream & , int) {
	}

	ClassDescription<FFLightKinematics> FFLightKinematics::initFFLightKinematics;
	// Definition of the static class description member.

	void FFLightKinematics::Init() {

	static ClassDocumentation<FFLightKinematics> documentation
	("FFLightKinematics implements massless splittings "
	"off a final-final dipole.");

	}

	diff --git a/DipoleShower/Kinematics/FFLightKinematics.h b/DipoleShower/Kinematics/FFLightKinematics.h
	--- a/DipoleShower/Kinematics/FFLightKinematics.h
	+++ b/DipoleShower/Kinematics/FFLightKinematics.h
	@@ -1,237 +1,200 @@
	// -- C++ --
	//
	// FFLightKinematics.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_FFLightKinematics_H
	#define HERWIG_FFLightKinematics_H
	//
	// This is the declaration of the FFLightKinematics class.
	//

	#include "DipoleSplittingKinematics.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* \ingroup DipoleShower
	* \author Simon Platzer
	*
	* \brief FFLightKinematics implements massless splittings
	* off a final-final dipole.
	*
	*/
	class FFLightKinematics: public DipoleSplittingKinematics {

	public:

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

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

	public:

	/**
	- * Return the boundaries in between the evolution
	- * variable random number is to be sampled; the lower
	- * cuoff is assumed to correspond to the infrared cutoff.
	- */
	- virtual pair<double,double> kappaSupport(const DipoleSplittingInfo& dIndex) const;
	-
	- /**
	- * Return the boundaries in between the momentum
	- * fraction random number is to be sampled.
	- */
	- virtual pair<double,double> xiSupport(const DipoleSplittingInfo& dIndex) const;
	-
	- /**
	- * Return the dipole scale associated to the
	- * given pair of emitter and spectator. This
	- * should be the invariant mass or absolute value
	- * final/final or initial/initial and the absolute
	- * value of the momentum transfer for intial/final or
	- * final/initial dipoles.
	- */
	- virtual Energy dipoleScale(const Lorentz5Momentum& pEmitter,
	- const Lorentz5Momentum& pSpectator) const;
	-
	- /**
	* Return the maximum pt for the given dipole scale.
	*/
	virtual Energy ptMax(Energy dScale,
	double emX, double specX,
	const DipoleIndex& dIndex,
	const DipoleSplittingKernel&) const;

	/**
	* Return the maximum virtuality for the given dipole scale.
	*/
	virtual Energy QMax(Energy dScale,
	double emX, double specX,
	- const DipoleIndex& dIndex) const;
	+ const DipoleIndex& dIndex,
	+ const DipoleSplittingKernel&) const;

	/**
	* Return the pt given a virtuality.
	*/
	virtual Energy PtFromQ(Energy scale, const DipoleSplittingInfo&) const;

	/**
	* Return the virtuality given a pt.
	*/
	virtual Energy QFromPt(Energy scale, const DipoleSplittingInfo&) const;

	/**
	- * Return the random number associated to
	- * the given pt.
	+ * Return the boundaries on the momentum fraction
	*/
	- virtual double ptToRandom(Energy pt, Energy dScale,
	- const DipoleIndex& dIndex) const;
	+ virtual pair<double,double> zBoundaries(Energy pt,
	+ const DipoleSplittingInfo& dInfo,
	+ const DipoleSplittingKernel& split) const;

	/**
	* Generate splitting variables given three random numbers
	* and the momentum fractions of the emitter and spectator.
	* Return true on success.
	*/
	virtual bool generateSplitting(double kappa, double xi, double phi,
	- DipoleSplittingInfo& dIndex);
	-
	- /**
	- * For the splitting products present in the given dipole splitting
	- * info object calculate the kinematics parameters and return the
	- * propagator factor.
	- */
	- virtual InvEnergy2 setKinematics(DipoleSplittingInfo&) const;
	-
	- /**
	- * For the splitting parameters given in the dipole splitting info
	- * object, calculate the phasespace Jacobian times the propagator
	- * factor.
	- */
	- virtual double jacobianTimesPropagator(const DipoleSplittingInfo&,
	- Energy) const;
	+ DipoleSplittingInfo& dIndex,
	+ const DipoleSplittingKernel&);

	/**
	* Generate the full kinematics given emitter and
	* spectator momentum and a previously completeted
	* DipoleSplittingInfo object.
	*/
	virtual void generateKinematics(const Lorentz5Momentum& pEmitter,
	const Lorentz5Momentum& pSpectator,
	const DipoleSplittingInfo& dInfo);

	public:

	/** @name Functions used by the persistent I/O system. */
	//@{
	/**
	* Function used to write out object persistently.
	* @param os the persistent output stream written to.
	*/
	void persistentOutput(PersistentOStream & os) const;

	/**
	* Function used to read in object persistently.
	* @param is the persistent input stream read from.
	* @param version the version number of the object when written.
	*/
	void persistentInput(PersistentIStream & is, int version);
	//@}

	/**
	* The standard Init function used to initialize the interfaces.
	* Called exactly once for each class by the class description system
	* before the main function starts or
	* when this class is dynamically loaded.
	*/
	static void Init();

	protected:

	/** @name Clone Methods. */
	//@{
	/**
	* Make a simple clone of this object.
	* @return a pointer to the new object.
	*/
	virtual IBPtr clone() const;

	/** Make a clone of this object, possibly modifying the cloned object
	* to make it sane.
	* @return a pointer to the new object.
	*/
	virtual IBPtr fullclone() const;
	//@}


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


	private:

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

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

	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

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

	#include "FFMassiveKinematics.h"
	#include "ThePEG/Interface/ClassDocumentation.h"

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

	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "Herwig++/DipoleShower/Base/DipoleSplittingInfo.h"
	#include "Herwig++/DipoleShower/Kernels/DipoleSplittingKernel.h"

	// TODO: remove after verification
	// only for checking for NaN or inf
	#include <gsl/gsl_math.h>

	using namespace Herwig;

	FFMassiveKinematics::FFMassiveKinematics()
	: DipoleSplittingKinematics() {}

	FFMassiveKinematics::~FFMassiveKinematics() {}

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

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

	pair<double,double> FFMassiveKinematics::kappaSupport(const DipoleSplittingInfo&) const {
	return make_pair(0.0,1.0);
	}

	pair<double,double> FFMassiveKinematics::xiSupport(const DipoleSplittingInfo& split) const {
	double c = sqrt(1.-4.*sqr(IRCutoff()/generator()->maximumCMEnergy()));
	if ( split.index().emitterData()->id() == ParticleID::g ) {
	if ( split.emissionData()->id() != ParticleID::g )
	return make_pair(0.5(1.-c),0.5(1.+c));
	double b = log((1.+c)/(1.-c));
	return make_pair(-b,b);
	}
	return make_pair(-log(0.5(1.+c)),-log(0.5(1.-c)));
	}

	Energy FFMassiveKinematics::dipoleScale(const Lorentz5Momentum& pEmitter,
	const Lorentz5Momentum& pSpectator) const {
	return (pEmitter+pSpectator).m();
	}

	Energy FFMassiveKinematics::ptMax(Energy dScale,
	double, double,
	const DipoleIndex& ind,
	const DipoleSplittingKernel& split) const {
	double mui2 = sqr( split.emitter(ind)->mass() / dScale );
	double mu2 = sqr( split.emission(ind)->mass() / dScale );
	double muj2 = sqr( split.spectator(ind)->mass() / dScale );

	// stolen from generateSplitting
	Energy ptmax = rootOfKallen( mui2, mu2, sqr(1.-sqrt(muj2)) ) /
	( 2.-2.sqrt(muj2) ) dScale;

	return ptmax > 0.GeV ? ptmax : 0.GeV;
	}

	Energy FFMassiveKinematics::QMax(Energy dScale,
	double, double,
	- const DipoleIndex& ind) const {
	+ const DipoleIndex& ind,
	+ const DipoleSplittingKernel&) const {
	assert(false && "implementation missing");
	double Muj = ind.spectatorData()->mass() / dScale;
	return dScale * ( 1.-2.*Muj+sqr(Muj) );
	}

	// relict
	Energy FFMassiveKinematics::PtFromQ(Energy, const DipoleSplittingInfo&) const {
	assert(false && "implementation missing");
	return 0.*GeV;
	}
	Energy FFMassiveKinematics::QFromPt(Energy, const DipoleSplittingInfo&) const {
	assert(false && "implementation missing");
	return 0.*GeV;
	}

	double FFMassiveKinematics::ptToRandom(Energy pt, Energy,
	- const DipoleIndex&) const {
	+ double,double,
	+ const DipoleIndex&,
	+ const DipoleSplittingKernel&) const {
	return log(pt/IRCutoff()) / log(0.5 * generator()->maximumCMEnergy()/IRCutoff());
	}

	// own kinematics close to Dinsdale,Ternick,Weinzierl
	bool FFMassiveKinematics::generateSplitting(double kappa, double xi, double rphi,
	- DipoleSplittingInfo& info) {
	+ DipoleSplittingInfo& info,
	+ const DipoleSplittingKernel&) {

	Energy pt = IRCutoff() * pow(0.5 * generator()->maximumCMEnergy()/IRCutoff(),kappa);

	if ( pt > info.hardPt() \|\| pt < IRCutoff() ) {
	jacobian(0.0);
	return false;
	}

	double z;
	double mapZJacobian;

	if ( info.index().emitterData()->id() == ParticleID::g ) {
	if ( info.emissionData()->id() != ParticleID::g ) {
	z = xi;
	mapZJacobian = 1.;
	} else {
	z = exp(xi)/(1.+exp(xi));
	mapZJacobian = z*(1.-z);
	}
	} else {
	z = 1.-exp(-xi);
	mapZJacobian = 1.-z;
	}

	// masses
	double mui2 = sqr( info.emitterData()->mass() / info.scale() );
	double mu2 = sqr( info.emissionData()->mass() / info.scale() );
	double muj2 = sqr( info.spectatorData()->mass() / info.scale() );
	double Mui2 = 0.;
	if ( info.emitterData()->id() + info.emissionData()->id() == 0 ) Mui2 = 0.; // gluon
	else Mui2 = mui2; // (anti)quark
	double Muj2 = muj2;

	// new: 2011-08-31
	// 2011-11-08: this does happen
	if( sqrt(mui2)+sqrt(mu2)+sqrt(muj2) > 1. ){
	jacobian(0.0);
	return false;
	}

	double bar = 1.-mui2-mu2-muj2;
	double y = ( sqr( pt / info.scale() ) + sqr(1.-z)mui2 + zz*mu2 ) /
	(z(1.-z)bar);

	// do this here for simplicity
	Energy ptmax1 = rootOfKallen( mui2, mu2, sqr(1.-sqrt(muj2)) ) /
	( 2.-2.sqrt(muj2) ) info.scale();
	Energy auxHardPt = ptmax1 > info.hardPt() ? info.hardPt() : ptmax1;
	// 2011-11-09
	assert(ptmax1>=info.hardPt());

	// phasespace constraint to incorporate ptMax
	double zp1 = ( 1.+mui2-mu2+muj2-2.*sqrt(muj2) +
	rootOfKallen(mui2,mu2,sqr(1-sqrt(muj2))) *
	sqrt( 1.-sqr(pt/auxHardPt) ) ) /
	( 2.*sqr(1.-sqrt(muj2)) );
	double zm1 = ( 1.+mui2-mu2+muj2-2.*sqrt(muj2) -
	rootOfKallen(mui2,mu2,sqr(1-sqrt(muj2))) *
	sqrt( 1.-sqr(pt/auxHardPt) ) ) /
	( 2.*sqr(1.-sqrt(muj2)) );
	if ( z > zp1 \|\| z < zm1 \|\|
	pt > auxHardPt) {
	jacobian(0.0);
	return false;
	}

	// kinematic phasespace boundaries for (y,z)
	// same as in Dittmaier hep-ph/9904440v2 (equivalent to CS)
	double ym = 2.sqrt(mui2)sqrt(mu2)/bar;
	double yp = 1. - 2.sqrt(muj2)(1.-sqrt(muj2))/bar;
	if ( y < ym \|\| y > yp ) {
	jacobian(0.0);
	return false;
	}

	double zm = ( (2.mui2+bary)(1.-y) - sqrt(yy-ymym)sqrt(sqr(2.muj2+bar-bary)-4.*muj2) ) /
	( 2.(1.-y)(mui2+mu2+bar*y) );
	double zp = ( (2.mui2+bary)(1.-y) + sqrt(yy-ymym)sqrt(sqr(2.muj2+bar-bary)-4.*muj2) ) /
	( 2.(1.-y)(mui2+mu2+bar*y) );

	if ( z < zm \|\| z > zp ) {
	jacobian(0.0);
	return false;
	}

	double phi = 2.Constants::pirphi;

	jacobian( 2. * mapZJacobian * (1.-y) *
	log(0.5 * generator()->maximumCMEnergy()/IRCutoff()) *
	bar / rootOfKallen(1.,Mui2,Muj2) );

	lastPt(pt);
	lastZ(z);
	lastPhi(phi);

	if ( theMCCheck )
	theMCCheck->book(1.,1.,info.scale(),info.hardPt(),pt,z,jacobian());

	return true;

	}

	-InvEnergy2 FFMassiveKinematics::setKinematics(DipoleSplittingInfo& split) const {
	-
	- // masses
	- double mui2 = sqr( split.emitterData()->mass() / split.scale() );
	- double mu2 = sqr( split.emissionData()->mass() / split.scale() );
	- double muj2 = sqr( split.spectatorData()->mass() / split.scale() );
	- double Mui2 = 0.;
	- if ( split.emitterData()->id() + split.emissionData()->id() == 0 ) Mui2 = 0.; // gluon
	- else Mui2 = mui2; // (anti)quark
	-
	- split.splittingKinematics(const_cast<FFMassiveKinematics*>(this));
	-
	- Lorentz5Momentum emitter = split.splitEmitter()->momentum();
	- Lorentz5Momentum emission = split.emission()->momentum();
	- Lorentz5Momentum spectator = split.splitSpectator()->momentum();
	-
	- Energy2 scale = (emitter+emission+spectator).m2();
	- split.scale(sqrt(scale));
	-
	- double y = 2.emissionemitter / scale / (1.-mui2-mu2-muj2);
	- double z = emitterspectator / (emitterspectator + emission*spectator);
	-
	- split.lastPt( split.scale() * sqrt( y * (1.-mui2-mu2-muj2) * z(1.-z) - sqr(1.-z)mui2 - sqr(z)*mu2 ) );
	- split.lastZ(z);
	-
	- split.hardPt(split.lastPt());
	-
	- if ( split.hardPt() > IRCutoff() ) {
	- split.continuesEvolving();
	- } else {
	- split.didStopEvolving();
	- }
	-
	- return 1./((emitter+emission).m2()-Mui2*sqr(split.scale()));
	-
	-}
	-
	-double FFMassiveKinematics::
	-jacobianTimesPropagator(const DipoleSplittingInfo&,
	- Energy) const {
	- assert(false && "implementation missing");
	- return 0.;
	-}
	-
	-
	// kinematics close to Dinsdale,Ternick,Weinzierl
	// revised 2011-08-22
	// revised 2011-11-06
	void FFMassiveKinematics::generateKinematics(const Lorentz5Momentum& pEmitter,
	const Lorentz5Momentum& pSpectator,
	const DipoleSplittingInfo& dInfo) {

	double z = dInfo.lastZ();
	Energy pt = dInfo.lastPt();

	// masses
	double mui2 = sqr( dInfo.emitterData()->mass() / dInfo.scale() );
	double mu2 = sqr( dInfo.emissionData()->mass() / dInfo.scale() );
	double muj2 = sqr( dInfo.spectatorData()->mass() / dInfo.scale() );

	double y = ( sqr( pt / dInfo.scale() ) + sqr(1.-z)mui2 + zz*mu2 ) /
	(z(1.-z)(1.-mui2-mu2-muj2));

	Energy2 sbar = sqr(dInfo.scale()) *(1.-mui2-mu2-muj2);

	// CMF: particle energies
	Energy Ei = ( sbar(1.-(1.-z)(1.-y)) + 2.sqr(dInfo.scale())mui2 ) / (2.*dInfo.scale());
	Energy E = ( sbar(1.- z (1.-y)) + 2.sqr(dInfo.scale())mu2 ) / (2.*dInfo.scale());
	Energy Ej = ( sbar(1.- y ) + 2.sqr(dInfo.scale())muj2 ) / (2.dInfo.scale());
	// CMF: momenta in z-direction (axis of pEmitter & pSpectator)
	Energy qi3 = (2.EiEj-z(1.-y)sbar ) / 2./sqrt(EjEj-sqr(dInfo.scale())muj2);
	Energy q3 = (2.E Ej-(1.-z)(1.-y)sbar) / 2./sqrt(EjEj-sqr(dInfo.scale())muj2);
	Energy qj3 = sqrt( sqr(Ej) - sqr(dInfo.scale())*muj2 );

	// get z axis in the dipole's CMF which is parallel to pSpectator
	Boost toCMF = (pEmitter+pSpectator).findBoostToCM();
	Lorentz5Momentum pjAux = pSpectator; pjAux.boost(toCMF);
	ThreeVector<double> pjAxis = pjAux.vect().unit();

	// set the momenta in this special reference frame
	// note that pt might in some cases differ from the physical pt!
	// phi is defined exactly as in getKt
	Energy ptResc = sqrt( sqr(Ei)-sqr(dInfo.scale())*mui2-sqr(qi3) );
	Lorentz5Momentum em ( ptResccos(dInfo.lastPhi()), -ptRescsin(dInfo.lastPhi()), qi3, Ei );
	Lorentz5Momentum emm ( -ptResccos(dInfo.lastPhi()), ptRescsin(dInfo.lastPhi()), q3, E );
	Lorentz5Momentum spe ( 0.GeV, 0.GeV, qj3, Ej );

	// output the mismatch between pt and physical pt
	// ofstream output1("ptDiffOnPtAxis-uub-m.dat",ofstream::app);
	// ofstream output2("ptDiffOnCosAxis-uub-m.dat",ofstream::app);
	// if( abs(dInfo.spectatorData()->id())==5 && dInfo.emitterData()->id()+dInfo.emissionData()->id()==0 &&
	// abs(dInfo.emitterData()->id())==1 ) {
	// output1 << pt/dInfo.scale() << " " << abs(ptResc-pt)/(ptResc+pt) << " " << endl;
	// output2 << em.vect().unit()*emm.vect().unit() << " " << abs(ptResc-pt)/(ptResc+pt) << " " << endl;
	// }
	// output1.close(); output2.close();

	// rotate back
	em.rotateUz (pjAxis);
	emm.rotateUz(pjAxis);
	spe.rotateUz(pjAxis);

	// boost back
	em.boost (-toCMF);
	emm.boost(-toCMF);
	spe.boost(-toCMF);

	// mass shells, rescale energy
	em.setMass(dInfo.scale()*sqrt(mui2));
	em.rescaleEnergy();
	emm.setMass(dInfo.scale()*sqrt(mu2));
	emm.rescaleEnergy();
	spe.setMass(dInfo.scale()*sqrt(muj2));
	spe.rescaleEnergy();

	// book
	emitterMomentum(em);
	emissionMomentum(emm);
	spectatorMomentum(spe);

	// TODO: remove
	// 2011-11-09: never occurred
	if(em.t()/GeV>=0. && emm.t()/GeV>=0. && spe.t()/GeV>=0.);
	else cout << "FFMassiveKinematics::generateKinematics momenta corrupt" << endl;

	// 2011-11-03 LEP run with full masses:
	// x,y,t no problem
	// z order > 5.e-7 happend 41 times in 10000 runs
	// maximum was 2.5e-6

	// double order=2.4e-6;
	// Lorentz5Momentum pDif=em+emm+spe-pEmitter-pSpectator, pSum=em+emm+spe+pEmitter+pSpectator;
	// if(abs(pDif.x()/(pSum.x()==ZERO ? 1.*GeV : pSum.x())) <order \|\| (pDif-pSum).x()==ZERO);
	// else cout << "FFMassiveKinematics::generateKinematics momenta corrupt: x " <<
	// abs(pDif.x()/(pSum.x()==ZERO ? 1.*GeV : pSum.x())) << endl <<
	// " " << em/GeV << " " << emm/GeV << " " << spe/GeV << endl <<
	// " " << pEmitter/GeV << " " << pSpectator/GeV << endl;
	// if(abs(pDif.y()/(pSum.y()==ZERO ? 1.*GeV : pSum.y())) <order \|\| (pDif-pSum).y()==ZERO);
	// else cout << "FFMassiveKinematics::generateKinematics momenta corrupt: y " <<
	// abs(pDif.y()/(pSum.y()==ZERO ? 1.*GeV : pSum.y())) << endl;
	// if(abs(pDif.z()/(pSum.z()==ZERO ? 1.*GeV : pSum.z())) <order \| (pDif-pSum).z()==ZERO);
	// else cout << "FFMassiveKinematics::generateKinematics momenta corrupt: z " <<
	// abs(pDif.z()/(pSum.z()==ZERO ? 1.*GeV : pSum.z())) << endl <<
	// " " << em/GeV << " " << emm/GeV << " " << spe/GeV << endl <<
	// " " << pEmitter/GeV << " " << pSpectator/GeV << endl;
	// if(abs(pDif.t()/(pSum.t()==ZERO ? 1.*GeV : pSum.t())) <order \|\| (pDif-pSum).t()==ZERO);
	// else cout << "FFMassiveKinematics::generateKinematics momenta corrupt: t " <<
	// // abs(pDif.t()/(pSum.t()==ZERO ? 1.*GeV : pSum.t())) << endl;
	// cout << endl;

	}

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


	void FFMassiveKinematics::persistentOutput(PersistentOStream & ) const {
	}

	void FFMassiveKinematics::persistentInput(PersistentIStream & , int) {
	}

	ClassDescription<FFMassiveKinematics> FFMassiveKinematics::initFFMassiveKinematics;
	// Definition of the static class description member.

	void FFMassiveKinematics::Init() {

	static ClassDocumentation<FFMassiveKinematics> documentation
	("FFMassiveKinematics implements massive splittings "
	"off a final-final dipole.");

	}

	diff --git a/DipoleShower/Kinematics/FFMassiveKinematics.h b/DipoleShower/Kinematics/FFMassiveKinematics.h
	--- a/DipoleShower/Kinematics/FFMassiveKinematics.h
	+++ b/DipoleShower/Kinematics/FFMassiveKinematics.h
	@@ -1,263 +1,261 @@
	// -- C++ --
	//
	// FFMassiveKinematics.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_FFMassiveKinematics_H
	#define HERWIG_FFMassiveKinematics_H
	//
	// This is the declaration of the FFMassiveKinematics class.
	//

	#include "DipoleSplittingKinematics.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* \ingroup DipoleShower
	* \author Simon Platzer
	*
	* \brief FFMassiveKinematics implements massive splittings
	* off a final-final dipole.
	*
	*/
	class FFMassiveKinematics: public DipoleSplittingKinematics {

	public:

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

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

	public:

	/**
	* Return the boundaries in between the evolution
	* variable random number is to be sampled; the lower
	* cuoff is assumed to correspond to the infrared cutoff.
	*/
	virtual pair<double,double> kappaSupport(const DipoleSplittingInfo& dIndex) const;

	/**
	* Return the boundaries in between the momentum
	* fraction random number is to be sampled.
	*/
	virtual pair<double,double> xiSupport(const DipoleSplittingInfo& dIndex) const;

	/**
	+ * Return the boundaries on the momentum fraction
	+ */
	+ virtual pair<double,double> zBoundaries(Energy,
	+ const DipoleSplittingInfo&,
	+ const DipoleSplittingKernel&) const {
	+ return make_pair(0.0,1.0);
	+ }
	+
	+ /**
	* Return the dipole scale associated to the
	* given pair of emitter and spectator. This
	* should be the invariant mass or absolute value
	* final/final or initial/initial and the absolute
	* value of the momentum transfer for intial/final or
	* final/initial dipoles.
	*/
	virtual Energy dipoleScale(const Lorentz5Momentum& pEmitter,
	const Lorentz5Momentum& pSpectator) const;

	/**
	* Return the maximum pt for the given dipole scale.
	*/
	virtual Energy ptMax(Energy dScale,
	double emX, double specX,
	const DipoleIndex& ind,
	const DipoleSplittingKernel& split) const;

	/**
	* Return the maximum virtuality for the given dipole scale.
	*/
	virtual Energy QMax(Energy dScale,
	double emX, double specX,
	- const DipoleIndex& dIndex) const;
	+ const DipoleIndex& dIndex,
	+ const DipoleSplittingKernel& split) const;

	/**
	* Return the pt given a virtuality.
	*/
	virtual Energy PtFromQ(Energy scale, const DipoleSplittingInfo&) const;

	/**
	* Return the virtuality given a pt.
	*/
	virtual Energy QFromPt(Energy scale, const DipoleSplittingInfo&) const;

	/**
	* Return the random number associated to
	* the given pt.
	*/
	virtual double ptToRandom(Energy pt, Energy dScale,
	- const DipoleIndex& dIndex) const;
	+ double emX, double specX,
	+ const DipoleIndex& dIndex,
	+ const DipoleSplittingKernel& split) const;

	/**
	* Generate splitting variables given three random numbers
	* and the momentum fractions of the emitter and spectator.
	* Return true on success.
	*/
	virtual bool generateSplitting(double kappa, double xi, double phi,
	- DipoleSplittingInfo& dIndex);
	-
	- /**
	- * For the splitting products present in the given dipole splitting
	- * info object calculate the kinematics parameters and return the
	- * propagator factor.
	- */
	- virtual InvEnergy2 setKinematics(DipoleSplittingInfo&) const;
	-
	- /**
	- * For the splitting parameters given in the dipole splitting info
	- * object, calculate the phasespace Jacobian times the propagator
	- * factor.
	- */
	- virtual double jacobianTimesPropagator(const DipoleSplittingInfo&,
	- Energy) const;
	+ DipoleSplittingInfo& dIndex,
	+ const DipoleSplittingKernel& split);

	/**
	* Generate the full kinematics given emitter and
	* spectator momentum and a previously completeted
	* DipoleSplittingInfo object.
	*/
	virtual void generateKinematics(const Lorentz5Momentum& pEmitter,
	const Lorentz5Momentum& pSpectator,
	const DipoleSplittingInfo& dInfo);

	public:

	/**
	* Triangular / Kallen function
	*/
	template <class T>
	inline double rootOfKallen (T a, T b, T c) const {
	return sqrt( aa + bb + cc - 2.( ab+ac+b*c ) ); }

	/**
	* Perform a rotation on both momenta such that the first one will
	* point along the (positive) z axis. Rotate back to the original
	* reference frame by applying rotateUz(returnedVector) to each momentum.
	*/
	ThreeVector<double> rotateToZ (Lorentz5Momentum& pTarget, Lorentz5Momentum& p1){
	ThreeVector<double> oldAxis = pTarget.vect().unit();
	double ct = oldAxis.z(); double st = sqrt( 1.-sqr(ct) ); // cos,sin(theta)
	double cp = oldAxis.x()/st; double sp = oldAxis.y()/st; // cos,sin(phi)
	pTarget.setZ( pTarget.vect().mag() ); pTarget.setX( 0.GeV ); pTarget.setY( 0.GeV );
	Lorentz5Momentum p1old = p1;
	p1.setX( spp1old.x() - cpp1old.y() );
	p1.setY( ctcpp1old.x() + ctspp1old.y() - st*p1old.z() );
	p1.setZ( stcpp1old.x() + stspp1old.y() + ct*p1old.z() );
	return oldAxis;
	}

	public:

	/** @name Functions used by the persistent I/O system. */
	//@{
	/**
	* Function used to write out object persistently.
	* @param os the persistent output stream written to.
	*/
	void persistentOutput(PersistentOStream & os) const;

	/**
	* Function used to read in object persistently.
	* @param is the persistent input stream read from.
	* @param version the version number of the object when written.
	*/
	void persistentInput(PersistentIStream & is, int version);
	//@}

	/**
	* The standard Init function used to initialize the interfaces.
	* Called exactly once for each class by the class description system
	* before the main function starts or
	* when this class is dynamically loaded.
	*/
	static void Init();

	protected:

	/** @name Clone Methods. */
	//@{
	/**
	* Make a simple clone of this object.
	* @return a pointer to the new object.
	*/
	virtual IBPtr clone() const;

	/** Make a clone of this object, possibly modifying the cloned object
	* to make it sane.
	* @return a pointer to the new object.
	*/
	virtual IBPtr fullclone() const;
	//@}


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


	private:

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

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

	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

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

	#include "FILightKinematics.h"
	#include "ThePEG/Interface/ClassDocumentation.h"


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

	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "Herwig++/DipoleShower/Base/DipoleSplittingInfo.h"

	using namespace Herwig;

	FILightKinematics::FILightKinematics()
	: DipoleSplittingKinematics() {}

	FILightKinematics::~FILightKinematics() {}

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

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

	-pair<double,double> FILightKinematics::kappaSupport(const DipoleSplittingInfo&) const {
	- return make_pair(0.0,1.0);
	-}
	-
	-pair<double,double> FILightKinematics::xiSupport(const DipoleSplittingInfo& split) const {
	- double c = sqrt(1.-4.*sqr(IRCutoff()/generator()->maximumCMEnergy()));
	- if ( split.index().emitterData()->id() == ParticleID::g ) {
	- if ( split.emissionData()->id() != ParticleID::g )
	- return make_pair(0.5(1.-c),0.5(1.+c));
	- double b = log((1.+c)/(1.-c));
	- return make_pair(-b,b);
	- }
	- return make_pair(-log(0.5(1.+c)),-log(0.5(1.-c)));
	-}
	-
	-Energy FILightKinematics::dipoleScale(const Lorentz5Momentum& pEmitter,
	- const Lorentz5Momentum& pSpectator) const {
	- return sqrt(2.(pEmitterpSpectator));
	-}
	-
	Energy FILightKinematics::ptMax(Energy dScale,
	double, double specX,
	const DipoleIndex&,
	const DipoleSplittingKernel&) const {
	return dScale * sqrt((1.-specX)/specX) /2.;
	}

	Energy FILightKinematics::QMax(Energy dScale,
	double, double specX,
	- const DipoleIndex&) const {
	+ const DipoleIndex&,
	+ const DipoleSplittingKernel&) const {
	return dScale * sqrt((1.-specX)/specX);
	}

	Energy FILightKinematics::PtFromQ(Energy scale, const DipoleSplittingInfo& split) const {
	double z = split.lastZ();
	return scalesqrt(z(1.-z));
	}

	Energy FILightKinematics::QFromPt(Energy scale, const DipoleSplittingInfo& split) const {
	double z = split.lastZ();
	return scale/sqrt(z*(1.-z));
	}

	-
	-double FILightKinematics::ptToRandom(Energy pt, Energy,
	- const DipoleIndex&) const {
	- return log(pt/IRCutoff()) / log(0.5 * generator()->maximumCMEnergy()/IRCutoff());
	+pair<double,double> FILightKinematics::zBoundaries(Energy pt,
	+ const DipoleSplittingInfo& dInfo,
	+ const DipoleSplittingKernel&) const {
	+ double s = sqrt(1.-sqr(pt/dInfo.hardPt()));
	+ return make_pair(0.5(1.-s),0.5(1.+s));
	}

	bool FILightKinematics::generateSplitting(double kappa, double xi, double rphi,
	- DipoleSplittingInfo& info) {
	+ DipoleSplittingInfo& info,
	+ const DipoleSplittingKernel& split) {

	if ( info.spectatorX() < xMin() ) {
	jacobian(0.0);
	return false;
	}

	- Energy pt = IRCutoff() * pow(0.5 * generator()->maximumCMEnergy()/IRCutoff(),kappa);
	+ double weight = 1.0;

	- if ( pt > info.hardPt() ) {
	+ Energy pt = generatePt(kappa,info.scale(),
	+ info.emitterX(),info.spectatorX(),
	+ info.index(),split,
	+ weight);
	+
	+ if ( pt < IRCutoff() \|\| pt > info.hardPt() ) {
	jacobian(0.0);
	return false;
	}

	- double z;
	- double mapZJacobian;
	+ double z = 0.0;

	if ( info.index().emitterData()->id() == ParticleID::g ) {
	if ( info.emissionData()->id() != ParticleID::g ) {
	- z = xi;
	- mapZJacobian = 1.;
	+ z = generateZ(xi,pt,FlatZ,
	+ info,split,weight);
	} else {
	- z = exp(xi)/(1.+exp(xi));
	- mapZJacobian = z*(1.-z);
	+ z = generateZ(xi,pt,OneOverZOneMinusZ,
	+ info,split,weight);
	}
	} else {
	- z = 1.-exp(-xi);
	- mapZJacobian = 1.-z;
	+ z = generateZ(xi,pt,OneOverOneMinusZ,
	+ info,split,weight);
	}

	- double s = z*(1.-z);
	+ double x = 1./(1.+sqr(pt/info.scale())/(z*(1.-z)));

	- double xs = info.spectatorX();
	-
	- double x = 1. / ( 1. + sqr(pt/info.scale()) / s );
	-
	- double zp = 0.5*(1.+sqrt(1.-sqr(pt/info.hardPt())));
	- double zm = 0.5*(1.-sqrt(1.-sqr(pt/info.hardPt())));
	-
	- if ( pt < IRCutoff() \|\|
	- pt > info.hardPt() \|\|
	- z > zp \|\| z < zm \|\|
	- x < xs ) {
	+ if ( z < 0.0 \|\| z > 1.0 \|\|
	+ x < info.spectatorX() \|\| x > 1.0 ) {
	jacobian(0.0);
	return false;
	}

	double phi = 2.Constants::pirphi;

	- jacobian( 2. * mapZJacobian * log(0.5 * generator()->maximumCMEnergy()/IRCutoff()));
	+ jacobian(weight);

	lastPt(pt);
	lastZ(z);
	lastPhi(phi);
	lastSpectatorZ(x);

	if ( theMCCheck )
	theMCCheck->book(1.,info.spectatorX(),info.scale(),info.hardPt(),pt,z,jacobian());

	return true;

	}

	-InvEnergy2 FILightKinematics::setKinematics(DipoleSplittingInfo& split) const {
	-
	- Lorentz5Momentum emitter = split.splitEmitter()->momentum();
	- Lorentz5Momentum emission = split.emission()->momentum();
	- Lorentz5Momentum spectator = split.splitSpectator()->momentum();
	-
	- split.splittingKinematics(const_cast<FILightKinematics*>(this));
	-
	- Energy2 scale = 2.(- emissionemitter + emissionspectator + emitterspectator);
	- split.scale(sqrt(scale));
	-
	- double x =
	- scale / (2.(emitterspectator + emission*spectator));
	- double z = emitterspectator / (emitterspectator + emission*spectator);
	-
	- split.lastPt(split.scale() * sqrt(z(1.-z)(1.-x)/x));
	- split.lastZ(z);
	-
	- split.hardPt(split.lastPt());
	-
	- if ( split.hardPt() > IRCutoff() ) {
	- split.continuesEvolving();
	- } else {
	- split.didStopEvolving();
	- }
	-
	- return 1./(2.x(emitter*emission));
	-
	-}
	-
	-double FILightKinematics::
	-jacobianTimesPropagator(const DipoleSplittingInfo&,
	- Energy) const {
	- assert(false && "implementation missing");
	- return 0.;
	-}
	-
	-
	void FILightKinematics::generateKinematics(const Lorentz5Momentum& pEmitter,
	const Lorentz5Momentum& pSpectator,
	const DipoleSplittingInfo& dInfo) {

	Energy pt = dInfo.lastPt();
	double z = dInfo.lastZ();
	+ double x = 1./(1.+sqr(pt/dInfo.scale())/(z*(1.-z)));

	Lorentz5Momentum kt =
	getKt (pSpectator, pEmitter, pt, dInfo.lastPhi(),true);

	- double ratio = sqr(pt/(-(pEmitter-pSpectator).m()));
	- double xInv = (1.+ratio/(z*(1.-z)));
	-
	- Lorentz5Momentum em = zpEmitter + (ratio/z)pSpectator + kt;
	+ Lorentz5Momentum em = zpEmitter + (1.-z)((1.-x)/x)*pSpectator + kt;
	em.setMass(0.*GeV);
	em.rescaleEnergy();

	- Lorentz5Momentum emm = (1.-z)pEmitter + (ratio/(1.-z))pSpectator - kt;
	+ Lorentz5Momentum emm = (1.-z)pEmitter + z((1.-x)/x)*pSpectator - kt;
	emm.setMass(0.*GeV);
	emm.rescaleEnergy();

	- Lorentz5Momentum spe = xInv*pSpectator;
	+ Lorentz5Momentum spe = (1./x)*pSpectator;
	spe.setMass(0.*GeV);
	spe.rescaleEnergy();

	emitterMomentum(em);
	emissionMomentum(emm);
	spectatorMomentum(spe);

	}

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


	void FILightKinematics::persistentOutput(PersistentOStream & ) const {
	}

	void FILightKinematics::persistentInput(PersistentIStream & , int) {
	}

	ClassDescription<FILightKinematics> FILightKinematics::initFILightKinematics;
	// Definition of the static class description member.

	void FILightKinematics::Init() {

	static ClassDocumentation<FILightKinematics> documentation
	("FILightKinematics implements massless splittings "
	"off a final-initial dipole.");

	}

	diff --git a/DipoleShower/Kinematics/FILightKinematics.h b/DipoleShower/Kinematics/FILightKinematics.h
	--- a/DipoleShower/Kinematics/FILightKinematics.h
	+++ b/DipoleShower/Kinematics/FILightKinematics.h
	@@ -1,237 +1,200 @@
	// -- C++ --
	//
	// FILightKinematics.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_FILightKinematics_H
	#define HERWIG_FILightKinematics_H
	//
	// This is the declaration of the FILightKinematics class.
	//

	#include "DipoleSplittingKinematics.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* \ingroup DipoleShower
	* \author Simon Platzer
	*
	* \brief FILightKinematics implements massless splittings
	* off a final-initial dipole.
	*
	*/
	class FILightKinematics: public DipoleSplittingKinematics {

	public:

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

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

	public:

	/**
	- * Return the boundaries in between the evolution
	- * variable random number is to be sampled; the lower
	- * cuoff is assumed to correspond to the infrared cutoff.
	- */
	- virtual pair<double,double> kappaSupport(const DipoleSplittingInfo& dIndex) const;
	-
	- /**
	- * Return the boundaries in between the momentum
	- * fraction random number is to be sampled.
	- */
	- virtual pair<double,double> xiSupport(const DipoleSplittingInfo& dIndex) const;
	-
	- /**
	- * Return the dipole scale associated to the
	- * given pair of emitter and spectator. This
	- * should be the invariant mass or absolute value
	- * final/final or initial/initial and the absolute
	- * value of the momentum transfer for intial/final or
	- * final/initial dipoles.
	- */
	- virtual Energy dipoleScale(const Lorentz5Momentum& pEmitter,
	- const Lorentz5Momentum& pSpectator) const;
	-
	- /**
	* Return the maximum pt for the given dipole scale.
	*/
	virtual Energy ptMax(Energy dScale,
	double emX, double specX,
	const DipoleIndex& dIndex,
	const DipoleSplittingKernel&) const;

	/**
	* Return the maximum virtuality for the given dipole scale.
	*/
	virtual Energy QMax(Energy dScale,
	double emX, double specX,
	- const DipoleIndex& dIndex) const;
	+ const DipoleIndex& dIndex,
	+ const DipoleSplittingKernel&) const;

	/**
	* Return the pt given a virtuality.
	*/
	virtual Energy PtFromQ(Energy scale, const DipoleSplittingInfo&) const;

	/**
	* Return the virtuality given a pt.
	*/
	virtual Energy QFromPt(Energy scale, const DipoleSplittingInfo&) const;

	/**
	- * Return the random number associated to
	- * the given pt.
	+ * Return the boundaries on the momentum fraction
	*/
	- virtual double ptToRandom(Energy pt, Energy dScale,
	- const DipoleIndex& dIndex) const;
	+ virtual pair<double,double> zBoundaries(Energy pt,
	+ const DipoleSplittingInfo& dInfo,
	+ const DipoleSplittingKernel& split) const;

	/**
	* Generate splitting variables given three random numbers
	* and the momentum fractions of the emitter and spectator.
	* Return true on success.
	*/
	virtual bool generateSplitting(double kappa, double xi, double phi,
	- DipoleSplittingInfo& dIndex);
	-
	- /**
	- * For the splitting products present in the given dipole splitting
	- * info object calculate the kinematics parameters and return the
	- * propagator factor.
	- */
	- virtual InvEnergy2 setKinematics(DipoleSplittingInfo&) const;
	-
	- /**
	- * For the splitting parameters given in the dipole splitting info
	- * object, calculate the phasespace Jacobian times the propagator
	- * factor.
	- */
	- virtual double jacobianTimesPropagator(const DipoleSplittingInfo&,
	- Energy) const;
	+ DipoleSplittingInfo& dIndex,
	+ const DipoleSplittingKernel&);

	/**
	* Generate the full kinematics given emitter and
	* spectator momentum and a previously completeted
	* DipoleSplittingInfo object.
	*/
	virtual void generateKinematics(const Lorentz5Momentum& pEmitter,
	const Lorentz5Momentum& pSpectator,
	const DipoleSplittingInfo& dInfo);

	public:

	/** @name Functions used by the persistent I/O system. */
	//@{
	/**
	* Function used to write out object persistently.
	* @param os the persistent output stream written to.
	*/
	void persistentOutput(PersistentOStream & os) const;

	/**
	* Function used to read in object persistently.
	* @param is the persistent input stream read from.
	* @param version the version number of the object when written.
	*/
	void persistentInput(PersistentIStream & is, int version);
	//@}

	/**
	* The standard Init function used to initialize the interfaces.
	* Called exactly once for each class by the class description system
	* before the main function starts or
	* when this class is dynamically loaded.
	*/
	static void Init();

	protected:

	/** @name Clone Methods. */
	//@{
	/**
	* Make a simple clone of this object.
	* @return a pointer to the new object.
	*/
	virtual IBPtr clone() const;

	/** Make a clone of this object, possibly modifying the cloned object
	* to make it sane.
	* @return a pointer to the new object.
	*/
	virtual IBPtr fullclone() const;
	//@}


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


	private:

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

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

	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

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

	#include "FIMassiveKinematics.h"
	#include "ThePEG/Interface/ClassDocumentation.h"


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

	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "Herwig++/DipoleShower/Base/DipoleSplittingInfo.h"
	#include "Herwig++/DipoleShower/Kernels/DipoleSplittingKernel.h"

	using namespace Herwig;

	FIMassiveKinematics::FIMassiveKinematics()
	: DipoleSplittingKinematics() {}

	FIMassiveKinematics::~FIMassiveKinematics() {}

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

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

	pair<double,double> FIMassiveKinematics::kappaSupport(const DipoleSplittingInfo&) const {
	return make_pair(0.0,1.0);
	}

	pair<double,double> FIMassiveKinematics::xiSupport(const DipoleSplittingInfo& split) const {
	double c = sqrt(1.-4.*sqr(IRCutoff()/generator()->maximumCMEnergy()));
	if ( split.index().emitterData()->id() == ParticleID::g ) {
	if ( split.emissionData()->id() != ParticleID::g )
	return make_pair(0.5(1.-c),0.5(1.+c));
	double b = log((1.+c)/(1.-c));
	return make_pair(-b,b);
	}
	return make_pair(-log(0.5(1.+c)),-log(0.5(1.-c)));
	}

	// sbar
	Energy FIMassiveKinematics::dipoleScale(const Lorentz5Momentum& pEmitter,
	const Lorentz5Momentum& pSpectator) const {
	return sqrt(2.(pEmitterpSpectator));
	}

	Energy FIMassiveKinematics::ptMax(Energy dScale,
	double, double specX,
	const DipoleIndex& ind,
	const DipoleSplittingKernel& split) const {
	Energy2 mi2 = sqr(split.emitter(ind)->mass());
	Energy2 m2 = sqr(split.emission(ind)->mass());
	// Energy2 Mi2 = split.emitter(int)->id() + split.emission(int)->id() == 0 ?
	// 0.*GeV2 : mi2;
	Energy2 Mi2 = mi2 == m2 ? 0.*GeV2 : mi2;

	// s^star/x
	Energy2 s = sqr(dScale) * (1.-specX)/specX + Mi2;

	Energy ptmax = .5 * sqrt(s) * rootOfKallen( s/s, mi2/s, m2/s );
	return ptmax > 0.GeV ? ptmax : 0.GeV;
	}

	// what is this? in FF it is given by y+dScale = sqrt( 2qiq / bar )->max
	Energy FIMassiveKinematics::QMax(Energy dScale,
	double, double specX,
	- const DipoleIndex&) const {
	+ const DipoleIndex&,
	+ const DipoleSplittingKernel&) const {
	assert(false && "implementation missing");
	cout << "FIMassiveKinematics::QMax called.\n" << flush;
	// this is sqrt( 2qi*q ) -> max;
	return dScale * sqrt((1.-specX)/specX);
	}

	Energy FIMassiveKinematics::PtFromQ(Energy, const DipoleSplittingInfo&) const {
	assert(false && "add this");
	return 0.0*GeV;
	}

	Energy FIMassiveKinematics::QFromPt(Energy, const DipoleSplittingInfo&) const {
	assert(false && "add this");
	return 0.0*GeV;
	}


	double FIMassiveKinematics::ptToRandom(Energy pt, Energy,
	- const DipoleIndex&) const {
	+ double,double,
	+ const DipoleIndex&,
	+ const DipoleSplittingKernel&) const {
	return log(pt/IRCutoff()) / log(0.5 * generator()->maximumCMEnergy()/IRCutoff());
	}

	bool FIMassiveKinematics::generateSplitting(double kappa, double xi, double rphi,
	- DipoleSplittingInfo& info) {
	+ DipoleSplittingInfo& info,
	+ const DipoleSplittingKernel&) {

	if ( info.spectatorX() < xMin() ) {
	jacobian(0.0);
	return false;
	}

	Energy pt = IRCutoff() * pow(0.5 * generator()->maximumCMEnergy()/IRCutoff(),kappa);

	if ( pt > info.hardPt() \|\| pt < IRCutoff() ) {
	jacobian(0.0);
	return false;
	}

	double z;
	double mapZJacobian;

	if ( info.index().emitterData()->id() == ParticleID::g ) {
	if ( info.emissionData()->id() != ParticleID::g ) {
	z = xi;
	mapZJacobian = 1.;
	} else {
	z = exp(xi)/(1.+exp(xi));
	mapZJacobian = z*(1.-z);
	}
	} else {
	z = 1.-exp(-xi);
	mapZJacobian = 1.-z;
	}

	// double s = z*(1.-z);

	// double xs = info.spectatorX();

	// double x = 1. / ( 1. + sqr(pt/info.scale()) / s );

	// double zp = 0.5*(1.+sqrt(1.-sqr(pt/info.hardPt())));
	// double zm = 0.5*(1.-sqrt(1.-sqr(pt/info.hardPt())));

	Energy2 mi2 = sqr(info.emitterData()->mass());
	Energy2 m2 = sqr(info.emissionData()->mass());
	Energy2 Mi2 = info.emitterData()->id()+info.emissionData()->id() == 0 ?
	0.*GeV2 : mi2;

	// s^star/x
	Energy2 s = sqr(info.scale()) * (1.-info.spectatorX())/info.spectatorX() + Mi2;

	double xs = info.spectatorX();
	double x = 1. / ( 1. +
	( sqr(pt) + (1.-z)mi2 + zm2 - z(1.-z)Mi2 ) /
	( z(1.-z)s ) );

	double zm1 = .5( 1.+(mi2-m2)/s - rootOfKallen(s/s,mi2/s,m2/s)
	sqrt( 1.-sqr(pt/info.hardPt()) ) );
	double zp1 = .5( 1.+(mi2-m2)/s + rootOfKallen(s/s,mi2/s,m2/s)
	sqrt( 1.-sqr(pt/info.hardPt()) ) );

	if ( // pt < IRCutoff() \|\|
	// pt > info.hardPt() \|\|
	z > zp1 \|\| z < zm1 \|\|
	x < xs ) {
	jacobian(0.0);
	return false;
	}

	// additional purely kinematic constraints
	double mui2 = x*mi2/sqr(info.scale());
	double mu2 = x*m2/sqr(info.scale());
	double Mui2 = x*Mi2/sqr(info.scale());
	double xp = 1. + Mui2 - sqr(sqrt(mui2)+sqrt(mu2));
	double zm = .5*( 1.-x+Mui2+mui2-mui2 -
	sqrt( sqr(1.-x+Mui2-mui2-mu2)-4.mui2mu2 ) ) /
	(1.-x+Mui2);
	double zp = .5*( 1.-x+Mui2+mui2-mui2 +
	sqrt( sqr(1.-x+Mui2-mui2-mu2)-4.mui2mu2 ) ) /
	(1.-x+Mui2);

	if (x > xp \|\|
	z > zp \|\| z < zm ) {
	jacobian(0.0);
	return false;
	}

	double phi = 2.Constants::pirphi;

	jacobian( 2. * mapZJacobian * log(0.5 * generator()->maximumCMEnergy()/IRCutoff()));

	lastPt(pt);
	lastZ(z);
	lastPhi(phi);
	lastSpectatorZ(x);

	if ( theMCCheck )
	theMCCheck->book(1.,info.spectatorX(),info.scale(),info.hardPt(),pt,z,jacobian());

	return true;

	}

	-InvEnergy2 FIMassiveKinematics::setKinematics(DipoleSplittingInfo& split) const {
	-
	- Energy2 mi2 = sqr(split.emitterData()->mass());
	- Energy2 m2 = sqr(split.emissionData()->mass());
	- Energy2 Mi2 = split.emitterData()->id() + split.emissionData()->id() == 0 ?
	- 0.*GeV2 : mi2;
	-
	- Lorentz5Momentum emitter = split.splitEmitter()->momentum();
	- Lorentz5Momentum emission = split.emission()->momentum();
	- Lorentz5Momentum spectator = split.splitSpectator()->momentum();
	-
	- split.splittingKinematics(const_cast<FIMassiveKinematics*>(this));
	-
	- // sbar is the relevant scale
	- // Energy2 scale = 2.(- emissionemitter + emissionspectator + emitterspectator);
	- Energy2 scale = Mi2 - (emitter+emission-spectator).m2();
	- split.scale(sqrt(scale));
	-
	- double x =
	- scale / (2.(emitterspectator + emission*spectator));
	- double z = emitterspectator / (emitterspectator + emission*spectator);
	-
	- split.lastPt( sqrt( z(1.-z)(1.-x)/x*scale -
	- ((1.-z)mi2+zm2-z(1.-z)Mi2) ) );
	- split.lastZ(z);
	-
	- split.hardPt(split.lastPt());
	-
	- if ( split.hardPt() > IRCutoff() ) {
	- split.continuesEvolving();
	- } else {
	- split.didStopEvolving();
	- }
	-
	- return 1./(2.x((emitter+emission).m2()-Mi2));
	-
	-}
	-
	-double FIMassiveKinematics::
	-jacobianTimesPropagator(const DipoleSplittingInfo&,
	- Energy) const {
	- assert(false && "implementation missing");
	- return 0.;
	-}
	-
	void FIMassiveKinematics::generateKinematics(const Lorentz5Momentum& pEmitter,
	const Lorentz5Momentum& pSpectator,
	const DipoleSplittingInfo& dInfo) {

	Energy pt = dInfo.lastPt();
	double z = dInfo.lastZ();

	Lorentz5Momentum kt =
	getKt (pSpectator, pEmitter, pt, dInfo.lastPhi(),true);

	Energy2 mi2 = sqr(dInfo.emitterData()->mass());
	Energy2 m2 = sqr(dInfo.emissionData()->mass());
	Energy2 Mi2 = dInfo.emitterData()->id() + dInfo.emissionData()->id() == 0 ?
	0.*GeV2 : mi2;

	double xInv = ( 1. +
	(ptpt+(1.-z)mi2+zm2-z(1.-z)*Mi2) /
	(z(1.-z)sqr(dInfo.scale())) );

	Lorentz5Momentum em = z*pEmitter +
	(sqr(pt)+mi2-zzMi2)/(zsqr(dInfo.scale()))pSpectator + kt;
	em.setMass(sqrt(mi2));
	em.rescaleEnergy();

	Lorentz5Momentum emm = (1.-z)*pEmitter +
	(ptpt+m2-sqr(1.-z)Mi2)/((1.-z)sqr(dInfo.scale()))pSpectator - kt;
	emm.setMass(sqrt(m2));
	emm.rescaleEnergy();

	Lorentz5Momentum spe = xInv*pSpectator;
	spe.setMass(ZERO);
	spe.rescaleEnergy();

	emitterMomentum(em);
	emissionMomentum(emm);
	spectatorMomentum(spe);

	}

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


	void FIMassiveKinematics::persistentOutput(PersistentOStream & ) const {
	}

	void FIMassiveKinematics::persistentInput(PersistentIStream & , int) {
	}

	ClassDescription<FIMassiveKinematics> FIMassiveKinematics::initFIMassiveKinematics;
	// Definition of the static class description member.

	void FIMassiveKinematics::Init() {

	static ClassDocumentation<FIMassiveKinematics> documentation
	("FIMassiveKinematics implements massless splittings "
	"off a final-initial dipole.");

	}

	diff --git a/DipoleShower/Kinematics/FIMassiveKinematics.h b/DipoleShower/Kinematics/FIMassiveKinematics.h
	--- a/DipoleShower/Kinematics/FIMassiveKinematics.h
	+++ b/DipoleShower/Kinematics/FIMassiveKinematics.h
	@@ -1,246 +1,244 @@
	// -- C++ --
	//
	// FILightKinematics.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_FILightKinematics_H
	#define HERWIG_FILightKinematics_H
	//
	// This is the declaration of the FILightKinematics class.
	//

	#include "DipoleSplittingKinematics.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* \ingroup DipoleShower
	* \author Simon Platzer, Martin Stoll
	*
	* \brief FIMassiveKinematics implements massless splittings
	* off a final-initial dipole.
	*
	*/
	class FIMassiveKinematics: public DipoleSplittingKinematics {

	public:

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

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

	public:

	/**
	* Return the boundaries in between the evolution
	* variable random number is to be sampled; the lower
	* cuoff is assumed to correspond to the infrared cutoff.
	*/
	virtual pair<double,double> kappaSupport(const DipoleSplittingInfo& dIndex) const;

	/**
	* Return the boundaries in between the momentum
	* fraction random number is to be sampled.
	*/
	virtual pair<double,double> xiSupport(const DipoleSplittingInfo& dIndex) const;

	/**
	+ * Return the boundaries on the momentum fraction
	+ */
	+ virtual pair<double,double> zBoundaries(Energy,
	+ const DipoleSplittingInfo&,
	+ const DipoleSplittingKernel&) const {
	+ return make_pair(0.0,1.0);
	+ }
	+
	+ /**
	* Return the dipole scale associated to the
	* given pair of emitter and spectator. This
	* should be the invariant mass or absolute value
	* final/final or initial/initial and the absolute
	* value of the momentum transfer for intial/final or
	* final/initial dipoles.
	*/
	virtual Energy dipoleScale(const Lorentz5Momentum& pEmitter,
	const Lorentz5Momentum& pSpectator) const;

	/**
	* Return the maximum pt for the given dipole scale.
	*/
	virtual Energy ptMax(Energy dScale,
	double emX, double specX,
	const DipoleIndex& dIndex,
	const DipoleSplittingKernel&) const;

	/**
	* Return the maximum virtuality for the given dipole scale.
	*/
	virtual Energy QMax(Energy dScale,
	double emX, double specX,
	- const DipoleIndex& dIndex) const;
	+ const DipoleIndex& dIndex,
	+ const DipoleSplittingKernel&) const;

	/**
	* Return the pt given a virtuality.
	*/
	virtual Energy PtFromQ(Energy scale, const DipoleSplittingInfo&) const;

	/**
	* Return the virtuality given a pt.
	*/
	virtual Energy QFromPt(Energy scale, const DipoleSplittingInfo&) const;

	/**
	* Return the random number associated to
	* the given pt.
	*/
	virtual double ptToRandom(Energy pt, Energy dScale,
	- const DipoleIndex& dIndex) const;
	+ double emX, double specX,
	+ const DipoleIndex& dIndex,
	+ const DipoleSplittingKernel&) const;

	/**
	* Generate splitting variables given three random numbers
	* and the momentum fractions of the emitter and spectator.
	* Return true on success.
	*/
	virtual bool generateSplitting(double kappa, double xi, double phi,
	- DipoleSplittingInfo& dIndex);
	-
	- /**
	- * For the splitting products present in the given dipole splitting
	- * info object calculate the kinematics parameters and return the
	- * propagator factor.
	- */
	- virtual InvEnergy2 setKinematics(DipoleSplittingInfo&) const;
	-
	- /**
	- * For the splitting parameters given in the dipole splitting info
	- * object, calculate the phasespace Jacobian times the propagator
	- * factor.
	- */
	- virtual double jacobianTimesPropagator(const DipoleSplittingInfo&,
	- Energy) const;
	+ DipoleSplittingInfo& dIndex,
	+ const DipoleSplittingKernel&);

	/**
	* Generate the full kinematics given emitter and
	* spectator momentum and a previously completeted
	* DipoleSplittingInfo object.
	*/
	virtual void generateKinematics(const Lorentz5Momentum& pEmitter,
	const Lorentz5Momentum& pSpectator,
	const DipoleSplittingInfo& dInfo);

	public:

	/**
	* Triangular / Kallen function
	*/
	template <class T>
	inline double rootOfKallen (T a, T b, T c) const {
	return sqrt( aa + bb + cc - 2.( ab+ac+b*c ) ); }

	public:

	/** @name Functions used by the persistent I/O system. */
	//@{
	/**
	* Function used to write out object persistently.
	* @param os the persistent output stream written to.
	*/
	void persistentOutput(PersistentOStream & os) const;

	/**
	* Function used to read in object persistently.
	* @param is the persistent input stream read from.
	* @param version the version number of the object when written.
	*/
	void persistentInput(PersistentIStream & is, int version);
	//@}

	/**
	* The standard Init function used to initialize the interfaces.
	* Called exactly once for each class by the class description system
	* before the main function starts or
	* when this class is dynamically loaded.
	*/
	static void Init();

	protected:

	/** @name Clone Methods. */
	//@{
	/**
	* Make a simple clone of this object.
	* @return a pointer to the new object.
	*/
	virtual IBPtr clone() const;

	/** Make a clone of this object, possibly modifying the cloned object
	* to make it sane.
	* @return a pointer to the new object.
	*/
	virtual IBPtr fullclone() const;
	//@}


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


	private:

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

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

	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

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

	#include "IFLightKinematics.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Switch.h"

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

	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "Herwig++/DipoleShower/Base/DipoleSplittingInfo.h"

	using namespace Herwig;

	IFLightKinematics::IFLightKinematics()
	: DipoleSplittingKinematics(), theCollinearScheme(false) {}

	IFLightKinematics::~IFLightKinematics() {}

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

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

	-pair<double,double> IFLightKinematics::kappaSupport(const DipoleSplittingInfo&) const {
	- return make_pair(0.0,1.0);
	-}
	-
	-pair<double,double> IFLightKinematics::xiSupport(const DipoleSplittingInfo& split) const {
	-
	- double c = sqrt(1.-4.*sqr(IRCutoff()/generator()->maximumCMEnergy()));
	-
	- if ( split.index().emitterData()->id() == ParticleID::g ) {
	- if ( split.emitterData()->id() == ParticleID::g ) {
	- double b = log((1.+c)/(1.-c));
	- return make_pair(-b,b);
	- } else {
	- return make_pair(log(0.5(1.-c)),log(0.5(1.+c)));
	- }
	- }
	-
	- if ( split.index().emitterData()->id() != ParticleID::g &&
	- split.emitterData()->id() != ParticleID::g ) {
	- return make_pair(-log(0.5(1.+c)),-log(0.5(1.-c)));
	- }
	-
	- return make_pair(0.5(1.-c),0.5(1.+c));
	-
	-}
	-
	-Energy IFLightKinematics::dipoleScale(const Lorentz5Momentum& pEmitter,
	- const Lorentz5Momentum& pSpectator) const {
	- return sqrt(2.(pEmitterpSpectator));
	-}
	-
	Energy IFLightKinematics::ptMax(Energy dScale,
	double emX, double,
	const DipoleIndex&,
	const DipoleSplittingKernel&) const {
	return dScale * sqrt((1.-emX)/emX) /2.;
	}

	Energy IFLightKinematics::QMax(Energy,
	double, double,
	- const DipoleIndex&) const {
	+ const DipoleIndex&,
	+ const DipoleSplittingKernel&) const {
	assert(false && "add this");
	return 0.0*GeV;
	}

	-Energy IFLightKinematics::PtFromQ(Energy, const DipoleSplittingInfo&) const {
	- assert(false && "add this");
	- return 0.0*GeV;
	+Energy IFLightKinematics::PtFromQ(Energy scale, const DipoleSplittingInfo& split) const {
	+ double z = split.lastZ();
	+ return scale*sqrt(1.-z);
	}

	-Energy IFLightKinematics::QFromPt(Energy, const DipoleSplittingInfo&) const {
	- assert(false && "add this");
	- return 0.0*GeV;
	+Energy IFLightKinematics::QFromPt(Energy scale, const DipoleSplittingInfo& split) const {
	+ double z = split.lastZ();
	+ return scale/sqrt(1.-z);
	}

	-
	-double IFLightKinematics::ptToRandom(Energy pt, Energy,
	- const DipoleIndex&) const {
	- return log(pt/IRCutoff()) / log(0.5 * generator()->maximumCMEnergy()/IRCutoff());
	+pair<double,double> IFLightKinematics::zBoundaries(Energy pt,
	+ const DipoleSplittingInfo& dInfo,
	+ const DipoleSplittingKernel&) const {
	+ double x = dInfo.emitterX();
	+ double s = sqrt(1.-sqr(pt/dInfo.hardPt()));
	+ return make_pair(0.5(1.+x-(1.-x)s),0.5(1.+x+(1.-x)s));
	}

	+
	bool IFLightKinematics::generateSplitting(double kappa, double xi, double rphi,
	- DipoleSplittingInfo& info) {
	+ DipoleSplittingInfo& info,
	+ const DipoleSplittingKernel& split) {

	if ( info.emitterX() < xMin() ) {
	jacobian(0.0);
	return false;
	}

	- Energy pt = IRCutoff() * pow(0.5 * generator()->maximumCMEnergy()/IRCutoff(),kappa);
	+ double weight = 1.0;

	- if ( sqr(pt) > sqr(info.hardPt())/(1.+4.*sqr(info.hardPt()/info.scale())) ) {
	+ Energy pt = generatePt(kappa,info.scale(),
	+ info.emitterX(),info.spectatorX(),
	+ info.index(),split,
	+ weight);
	+
	+ if ( pt < IRCutoff() \|\| pt > info.hardPt() ) {
	jacobian(0.0);
	return false;
	}

	- double z = 0.;
	- double mapZJacobian = 0.;
	+ double z = 0.0;

	if ( info.index().emitterData()->id() == ParticleID::g ) {
	if ( info.emitterData()->id() == ParticleID::g ) {
	- z = exp(xi)/(1.+exp(xi));
	- mapZJacobian = z*(1.-z);
	+ z = generateZ(xi,pt,OneOverZOneMinusZ,
	+ info,split,weight);
	} else {
	- z = exp(xi);
	- mapZJacobian = z;
	+ z = generateZ(xi,pt,OneOverZ,
	+ info,split,weight);
	}
	}

	if ( info.index().emitterData()->id() != ParticleID::g ) {
	if ( info.emitterData()->id() != ParticleID::g ) {
	- z = 1.-exp(-xi);
	- mapZJacobian = 1.-z;
	+ z = generateZ(xi,pt,OneOverOneMinusZ,
	+ info,split,weight);
	} else {
	- z = xi;
	- mapZJacobian = 1.;
	+ z = generateZ(xi,pt,FlatZ,
	+ info,split,weight);
	}
	}

	double ratio = sqr(pt/info.scale());

	- double x = ( z*(1.-z) - ratio ) / ( 1. - z - ratio );
	- double u = ratio/(1.-z);
	-
	- if ( x < 0. \|\| x > 1. \|\| u > 1. ) {
	+ double rho = 1. - 4.ratioz*(1.-z)/(1.-z+ratio);
	+ if ( rho < 0.0 ) {
	jacobian(0.0);
	return false;
	}

	- double xe = info.emitterX();
	+ double x = 0.5((1.-z+ratio)/ratio)(1.+sqrt(rho));
	+ double u = 0.5((1.-z+ratio)/(1.-z))(1.+sqrt(rho));

	- double zpx = 0.5*( 1.+ xe +
	- (1.-xe)sqrt(1.-sqr(2.pt/info.scale())/(1.-xe) ) );
	- double zmx = 0.5*( 1.+ xe -
	- (1.-xe)sqrt(1.-sqr(2.pt/info.scale())/(1.-xe) ) );
	-
	- double xq = sqr(pt/info.hardPt());
	-
	- double zpq = 0.5*( 1.+ xq +
	- (1.-xq)sqrt(1.-sqr(2.pt/info.scale())/(1.-xq) ) );
	- double zmq = 0.5*( 1.+ xq -
	- (1.-xq)sqrt(1.-sqr(2.pt/info.scale())/(1.-xq) ) );
	-
	- double zp = min(zpx,zpq);
	- double zm = max(zmx,zmq);
	-
	- if ( pt < IRCutoff() \|\|
	- pt > info.hardPt() \|\|
	- z > zp \|\| z < zm \|\|
	- x < xe ) {
	+ if ( x < info.emitterX() \|\| x > 1. \|\|
	+ u < 0. \|\| u > 1. ) {
	jacobian(0.0);
	return false;
	}

	double phi = 2.Constants::pirphi;

	- jacobian(2. * mapZJacobian * (1.-z)/(z(1.-z)-ratio) log(0.5 * generator()->maximumCMEnergy()/IRCutoff()));
	+ jacobian(weight*(1./z));

	lastPt(pt);
	lastZ(z);
	lastPhi(phi);
	lastEmitterZ(x);

	if ( theMCCheck )
	theMCCheck->book(info.emitterX(),1.,info.scale(),info.hardPt(),pt,z,jacobian());

	return true;

	}

	-InvEnergy2 IFLightKinematics::setKinematics(DipoleSplittingInfo& split) const {
	-
	- Lorentz5Momentum emitter = split.splitEmitter()->momentum();
	- Lorentz5Momentum emission = split.emission()->momentum();
	- Lorentz5Momentum spectator = split.splitSpectator()->momentum();
	-
	- split.splittingKinematics(const_cast<IFLightKinematics*>(this));
	-
	- Energy2 scale = 2.(emissionemitter - emissionspectator + emitterspectator);
	- split.scale(sqrt(scale));
	-
	- double x =
	- scale / (2.(emitteremission + emitter*spectator));
	- double u = emitteremission / (emitteremission + emitter*spectator);
	-
	- split.lastPt(split.scale() * sqrt(u(1.-u)(1.-x)));
	- split.lastZ(u+x*(1.-u));
	-
	- split.hardPt(split.lastPt());
	-
	- if ( split.hardPt() > IRCutoff() ) {
	- split.continuesEvolving();
	- } else {
	- split.didStopEvolving();
	- }
	-
	- return 1./(2.x(emitter*emission));
	-
	-}
	-
	-double IFLightKinematics::
	-jacobianTimesPropagator(const DipoleSplittingInfo&,
	- Energy) const {
	- assert(false && "implementation missing");
	- return 0.;
	-}
	-
	void IFLightKinematics::generateKinematics(const Lorentz5Momentum& pEmitter,
	const Lorentz5Momentum& pSpectator,
	const DipoleSplittingInfo& dInfo) {

	Energy pt = dInfo.lastPt();
	- double ratio = sqr(pt/-(pEmitter-pSpectator).m());
	double z = dInfo.lastZ();
	- double x = (z*(1.-z)-ratio)/(1.-z-ratio);
	- double u = ratio / (1.-z);

	- pt = -(pEmitter-pSpectator).m()sqrt(u(1.-u)*(1.-x)/x);
	+ double ratio = sqr(pt)/(2.pEmitterpSpectator);
	+ double rho = 1. - 4.ratioz*(1.-z)/(1.-z+ratio);
	+
	+ double x = 0.5((1.-z+ratio)/ratio)(1.+sqrt(rho));
	+ double u = 0.5((1.-z+ratio)/(1.-z))(1.+sqrt(rho));

	Lorentz5Momentum kt =
	getKt (pEmitter, pSpectator, pt, dInfo.lastPhi(),true);

	Lorentz5Momentum em;
	Lorentz5Momentum emm;
	Lorentz5Momentum spe;

	if ( !theCollinearScheme &&
	x > u && (1.-x)/(x-u) < 1. ) {

	em =
	((1.-u)/(x-u))pEmitter + ((u/x)(1.-x)/(x-u))*pSpectator - kt/(x-u);
	em.setMass(0.*GeV);
	em.rescaleEnergy();

	emm =
	((1.-x)/(x-u))pEmitter + ((u/x)(1.-u)/(x-u))*pSpectator - kt/(x-u);
	emm.setMass(0.*GeV);
	emm.rescaleEnergy();

	spe =
	(1.-u/x)*pSpectator;
	spe.setMass(0.*GeV);
	spe.rescaleEnergy();

	} else {

	em = (1./x)*pEmitter;

	emm = ((1.-x)(1.-u)/x)pEmitter + u*pSpectator + kt;
	emm.setMass(0.*GeV);
	emm.rescaleEnergy();

	spe = ((1.-x)u/x)pEmitter + (1.-u)*pSpectator - kt;
	spe.setMass(0.*GeV);
	spe.rescaleEnergy();

	}

	emitterMomentum(em);
	emissionMomentum(emm);
	spectatorMomentum(spe);

	}

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


	void IFLightKinematics::persistentOutput(PersistentOStream & os) const {
	os << theCollinearScheme;
	}

	void IFLightKinematics::persistentInput(PersistentIStream & is, int) {
	is >> theCollinearScheme;
	}

	ClassDescription<IFLightKinematics> IFLightKinematics::initIFLightKinematics;
	// Definition of the static class description member.

	void IFLightKinematics::Init() {

	static ClassDocumentation<IFLightKinematics> documentation
	("IFLightKinematics implements massless splittings "
	"off a initial-final dipole.");


	static Switch<IFLightKinematics,bool> interfaceCollinearScheme
	("CollinearScheme",
	"[experimental] Switch on or off the collinear scheme",
	&IFLightKinematics::theCollinearScheme, false, false, false);
	static SwitchOption interfaceCollinearSchemeOn
	(interfaceCollinearScheme,
	"On",
	"Switch on the collinear scheme.",
	true);
	static SwitchOption interfaceCollinearSchemeOff
	(interfaceCollinearScheme,
	"Off",
	"Switch off the collinear scheme",
	false);

	interfaceCollinearScheme.rank(-1);

	}

	diff --git a/DipoleShower/Kinematics/IFLightKinematics.h b/DipoleShower/Kinematics/IFLightKinematics.h
	--- a/DipoleShower/Kinematics/IFLightKinematics.h
	+++ b/DipoleShower/Kinematics/IFLightKinematics.h
	@@ -1,246 +1,209 @@
	// -- C++ --
	//
	// IFLightKinematics.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_IFLightKinematics_H
	#define HERWIG_IFLightKinematics_H
	//
	// This is the declaration of the IFLightKinematics class.
	//

	#include "DipoleSplittingKinematics.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* \ingroup DipoleShower
	* \author Simon Platzer
	*
	* \brief IFLightKinematics implements massless splittings
	* off an initial-final dipole.
	*
	* @see \ref IFLightKinematicsInterfaces "The interfaces"
	* defined for IFLightKinematics.
	*/
	class IFLightKinematics: public DipoleSplittingKinematics {

	public:

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

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

	public:

	/**
	- * Return the boundaries in between the evolution
	- * variable random number is to be sampled; the lower
	- * cuoff is assumed to correspond to the infrared cutoff.
	- */
	- virtual pair<double,double> kappaSupport(const DipoleSplittingInfo& dIndex) const;
	-
	- /**
	- * Return the boundaries in between the momentum
	- * fraction random number is to be sampled.
	- */
	- virtual pair<double,double> xiSupport(const DipoleSplittingInfo& dIndex) const;
	-
	- /**
	- * Return the dipole scale associated to the
	- * given pair of emitter and spectator. This
	- * should be the invariant mass or absolute value
	- * final/final or initial/initial and the absolute
	- * value of the momentum transfer for intial/final or
	- * final/initial dipoles.
	- */
	- virtual Energy dipoleScale(const Lorentz5Momentum& pEmitter,
	- const Lorentz5Momentum& pSpectator) const;
	-
	- /**
	* Return the maximum pt for the given dipole scale.
	*/
	virtual Energy ptMax(Energy dScale,
	double emX, double specX,
	const DipoleIndex&,
	const DipoleSplittingKernel&) const;

	/**
	* Return the maximum virtuality for the given dipole scale.
	*/
	virtual Energy QMax(Energy dScale,
	double emX, double specX,
	- const DipoleIndex& dIndex) const;
	+ const DipoleIndex& dIndex,
	+ const DipoleSplittingKernel&) const;

	/**
	* Return the pt given a virtuality.
	*/
	virtual Energy PtFromQ(Energy scale, const DipoleSplittingInfo&) const;

	/**
	* Return the virtuality given a pt.
	*/
	virtual Energy QFromPt(Energy scale, const DipoleSplittingInfo&) const;

	/**
	- * Return the random number associated to
	- * the given pt.
	+ * Return the boundaries on the momentum fraction
	*/
	- virtual double ptToRandom(Energy pt, Energy dScale,
	- const DipoleIndex& dIndex) const;
	+ virtual pair<double,double> zBoundaries(Energy pt,
	+ const DipoleSplittingInfo& dInfo,
	+ const DipoleSplittingKernel& split) const;

	/**
	* Generate splitting variables given three random numbers
	* and the momentum fractions of the emitter and spectator.
	* Return true on success.
	*/
	virtual bool generateSplitting(double kappa, double xi, double phi,
	- DipoleSplittingInfo& dIndex);
	-
	- /**
	- * For the splitting products present in the given dipole splitting
	- * info object calculate the kinematics parameters and return the
	- * propagator factor.
	- */
	- virtual InvEnergy2 setKinematics(DipoleSplittingInfo&) const;
	-
	- /**
	- * For the splitting parameters given in the dipole splitting info
	- * object, calculate the phasespace Jacobian times the propagator
	- * factor.
	- */
	- virtual double jacobianTimesPropagator(const DipoleSplittingInfo&,
	- Energy) const;
	+ DipoleSplittingInfo& dIndex,
	+ const DipoleSplittingKernel&);

	/**
	* Generate the full kinematics given emitter and
	* spectator momentum and a previously completeted
	* DipoleSplittingInfo object.
	*/
	virtual void generateKinematics(const Lorentz5Momentum& pEmitter,
	const Lorentz5Momentum& pSpectator,
	const DipoleSplittingInfo& dInfo);

	public:

	/** @name Functions used by the persistent I/O system. */
	//@{
	/**
	* Function used to write out object persistently.
	* @param os the persistent output stream written to.
	*/
	void persistentOutput(PersistentOStream & os) const;

	/**
	* Function used to read in object persistently.
	* @param is the persistent input stream read from.
	* @param version the version number of the object when written.
	*/
	void persistentInput(PersistentIStream & is, int version);
	//@}

	/**
	* The standard Init function used to initialize the interfaces.
	* Called exactly once for each class by the class description system
	* before the main function starts or
	* when this class is dynamically loaded.
	*/
	static void Init();

	protected:

	/** @name Clone Methods. */
	//@{
	/**
	* Make a simple clone of this object.
	* @return a pointer to the new object.
	*/
	virtual IBPtr clone() const;

	/** Make a clone of this object, possibly modifying the cloned object
	* to make it sane.
	* @return a pointer to the new object.
	*/
	virtual IBPtr fullclone() const;
	//@}


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


	private:

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

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

	private:

	/**
	* Wether or not to choose the `collinear' scheme
	*/
	bool theCollinearScheme;

	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

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

	#include "IFMassiveKinematics.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Switch.h"

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

	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "Herwig++/DipoleShower/Base/DipoleSplittingInfo.h"

	using namespace Herwig;

	IFMassiveKinematics::IFMassiveKinematics()
	: DipoleSplittingKinematics(), theCollinearScheme(false) {}

	IFMassiveKinematics::~IFMassiveKinematics() {}

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

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

	pair<double,double> IFMassiveKinematics::kappaSupport(const DipoleSplittingInfo&) const {
	return make_pair(0.0,1.0);
	}

	pair<double,double> IFMassiveKinematics::xiSupport(const DipoleSplittingInfo& split) const {

	double c = sqrt(1.-4.*sqr(IRCutoff()/generator()->maximumCMEnergy()));

	if ( split.index().emitterData()->id() == ParticleID::g ) {
	if ( split.emitterData()->id() == ParticleID::g ) {
	double b = log((1.+c)/(1.-c));
	return make_pair(-b,b);
	} else {
	return make_pair(log(0.5(1.-c)),log(0.5(1.+c)));
	}
	}

	if ( split.index().emitterData()->id() != ParticleID::g &&
	split.emitterData()->id() != ParticleID::g ) {
	return make_pair(-log(0.5(1.+c)),-log(0.5(1.-c)));
	}

	return make_pair(0.5(1.-c),0.5(1.+c));

	}

	// sbar
	Energy IFMassiveKinematics::dipoleScale(const Lorentz5Momentum& pEmitter,
	const Lorentz5Momentum& pSpectator) const {
	return sqrt(2.(pEmitterpSpectator));
	}

	Energy IFMassiveKinematics::ptMax(Energy dScale,
	double emX, double,
	const DipoleIndex&,
	const DipoleSplittingKernel&) const {
	return dScale * sqrt(1.-emX) /2.;
	}

	Energy IFMassiveKinematics::QMax(Energy,
	- double, double,
	- const DipoleIndex&) const {
	+ double, double,
	+ const DipoleIndex&,
	+ const DipoleSplittingKernel&) const {
	assert(false && "add this");
	return 0.0*GeV;
	}

	Energy IFMassiveKinematics::PtFromQ(Energy, const DipoleSplittingInfo&) const {
	assert(false && "add this");
	return 0.0*GeV;
	}

	Energy IFMassiveKinematics::QFromPt(Energy, const DipoleSplittingInfo&) const {
	assert(false && "add this");
	return 0.0*GeV;
	}


	double IFMassiveKinematics::ptToRandom(Energy pt, Energy,
	- const DipoleIndex&) const {
	+ double,double,
	+ const DipoleIndex&,
	+ const DipoleSplittingKernel&) const {
	return log(pt/IRCutoff()) / log(0.5 * generator()->maximumCMEnergy()/IRCutoff());
	}

	bool IFMassiveKinematics::generateSplitting(double kappa, double xi, double rphi,
	- DipoleSplittingInfo& info) {
	+ DipoleSplittingInfo& info,
	+ const DipoleSplittingKernel&) {

	if ( info.emitterX() < xMin() ) {
	jacobian(0.0);
	return false;
	}

	Energy pt = IRCutoff() * pow(0.5 * generator()->maximumCMEnergy()/IRCutoff(),kappa);

	if ( pt > info.hardPt() ) {
	jacobian(0.0);
	return false;
	}

	double z = 0.;
	double mapZJacobian = 0.;

	if ( info.index().emitterData()->id() == ParticleID::g ) {
	if ( info.emitterData()->id() == ParticleID::g ) {
	z = exp(xi)/(1.+exp(xi));
	mapZJacobian = z*(1.-z);
	} else {
	z = exp(xi);
	mapZJacobian = z;
	}
	}

	if ( info.index().emitterData()->id() != ParticleID::g ) {
	if ( info.emitterData()->id() != ParticleID::g ) {
	z = 1.-exp(-xi);
	mapZJacobian = 1.-z;
	} else {
	z = xi;
	mapZJacobian = 1.;
	}
	}

	double ratio = sqr(pt/info.scale());

	double alpha = 1. - 2.*sqr(info.spectatorData()->mass()/info.scale());

	double x = alpha == 1. ? ( z*(1.-z) - ratio ) / ( 1. - z - ratio ) :
	( sqr(alpha)ratio + 2.z - alpha*(1.+z) +
	alphasqrt( sqr(1.-z+alpharatio) - 4.ratio(1.-z) ) ) /
	(2.*(1.-alpha));
	double u = ( 1.-z + alpha*ratio -
	sqrt( sqr(1.-z+alpharatio) - 4.ratio*(1.-z) ) ) /
	(2.*(1.-z));
	// double x = ( z*(1.-z) - ratio ) / ( 1. - z - ratio );
	// double u = ratio/(1.-z);

	double up = (1.-x) /
	( 1.-x + x*sqr(info.spectatorData()->mass()/info.scale()) );

	if ( x < 0. \|\| x > 1. \|\| u > up ) {
	jacobian(0.0);
	return false;
	}

	double xe = info.emitterX();

	double zp = 0.5( alpha + xe - (alpha-1.)xe +
	alpha(1.-xe)sqrt(1.-sqr(pt/info.hardPt()) ) );
	double zm = 0.5( alpha + xe - (alpha-1.)xe +
	alpha(1.-xe)sqrt(1.-sqr(pt/info.hardPt()) ) );

	if ( pt < IRCutoff() \|\|
	pt > info.hardPt() \|\|
	z > zp \|\| z < zm \|\|
	x < xe ) {
	jacobian(0.0);
	return false;
	}

	double phi = 2.Constants::pirphi;

	jacobian(2. * mapZJacobian / x * log(0.5 * generator()->maximumCMEnergy()/IRCutoff()) * (1.-u)/(1.-2.u+uu*alpha) );

	lastPt(pt);
	lastZ(z);
	lastPhi(phi);
	lastEmitterZ(x);

	if ( theMCCheck )
	theMCCheck->book(info.emitterX(),1.,info.scale(),info.hardPt(),pt,z,jacobian());

	return true;

	}

	-InvEnergy2 IFMassiveKinematics::setKinematics(DipoleSplittingInfo& split) const {
	-
	- Lorentz5Momentum emitter = split.splitEmitter()->momentum();
	- Lorentz5Momentum emission = split.emission()->momentum();
	- Lorentz5Momentum spectator = split.splitSpectator()->momentum();
	-
	- split.splittingKinematics(const_cast<IFMassiveKinematics*>(this));
	-
	- // sbar
	- Energy2 scale = 2.(emissionemitter - emissionspectator + emitterspectator);
	- split.scale(sqrt(scale));
	-
	- double x =
	- scale / (2.(emitteremission + emitter*spectator));
	- double u = emitteremission / (emitteremission + emitter*spectator);
	-
	- split.lastPt(split.scale() * sqrt(u(1.-u)(1.-x)));
	- split.lastZ( x+u(1.-x)
	- ( 1. - 2.*sqr(split.spectatorData()->mass())/scale ) );
	-
	- split.hardPt(split.lastPt());
	-
	- if ( split.hardPt() > IRCutoff() ) {
	- split.continuesEvolving();
	- } else {
	- split.didStopEvolving();
	- }
	-
	- return 1./(2.x(emitter*emission));
	-
	-}
	-
	-double IFMassiveKinematics::
	-jacobianTimesPropagator(const DipoleSplittingInfo&,
	- Energy) const {
	- assert(false && "implementation missing");
	- return 0.;
	-}
	-
	void IFMassiveKinematics::generateKinematics(const Lorentz5Momentum& pEmitter,
	const Lorentz5Momentum& pSpectator,
	const DipoleSplittingInfo& dInfo) {
	Energy2 sbar = 2.pEmitterpSpectator;
	Energy pt = dInfo.lastPt();
	double ratio = pt*pt/sbar;
	double z = dInfo.lastZ();
	double x = (z*(1.-z)-ratio)/(1.-z-ratio);
	double u = ratio / (1.-z);

	pt = sqrt(sbaru(1.-u)*(1.-x)/x);

	Lorentz5Momentum kt =
	getKt (pEmitter, pSpectator, pt, dInfo.lastPhi(),true);

	Lorentz5Momentum em;
	Lorentz5Momentum emm;
	Lorentz5Momentum spe;

	Energy2 mj2 = dInfo.spectatorData()->mass()*dInfo.spectatorData()->mass();
	double alpha = 1. - 2.*mj2/sbar;

	// TODO: adjust phasespace boundary condition
	if ( !theCollinearScheme &&
	x > u && (1.-x)/(x-u) < 1. ) {

	double fkt = sqrt(sqr(x-u)+4.xu*mj2/sbar);

	// em =
	// ((1.-u)/(x-u))pEmitter + ((u/x)(1.-x)/(x-u))*pSpectator - kt/(x-u);
	Energy2 fa = (sbar(x+u-2.xz)+2.mj2xu) / sqrt(sqr(x-u)+4.xu*mj2/sbar);
	double a = (-sbar+fa) / (2.x(sbar-mj2));
	double ap = (sbar+alphafa) / (2.x*(sbar-mj2));
	em = appEmitter + apSpectator - fkt*kt;
	em.setMass(ZERO);
	em.rescaleEnergy();

	// emm =
	// ((1.-x)/(x-u))pEmitter + ((u/x)(1.-u)/(x-u))*pSpectator - kt/(x-u);
	Energy2 fb = abs(sbar(u(1.-u)-x(1.-x))+2.mj2xu) / sqrt(sqr(x-u)+4.xu*mj2/sbar);
	double b = (-sbar(1.-x-u)+fb) / (2.x*(sbar-mj2));
	double bp = (sbar(1.-x-u)+alphafb) / (2.x(sbar-mj2));
	emm = bppEmitter + bpSpectator + fkt*kt;
	emm.setMass(ZERO);
	emm.rescaleEnergy();

	// spe =
	// (1.-u/x)*pSpectator;
	Energy2 fc = sqrt(sqr(sbar(x-u))+4.sbarmj2x*u);
	double c = (sbar(x-u)-2.xmj2+fc) / (2.x*(sbar-mj2));
	double cp = (-sbar(x-u)+2.xmj2+alphafc) / (2.x(sbar-mj2));
	spe = cppEmitter + cpSpectator;
	spe.setMass(dInfo.spectatorData()->mass());
	spe.rescaleEnergy();

	} else {

	em = (1./x)*pEmitter;
	em.setMass(ZERO);
	em.rescaleEnergy();

	// emm = ((1.-x)(1.-u)/x)pEmitter + u*pSpectator + kt;
	emm = (ptpt-uumj2)/(usbar)*pEmitter +
	u*pSpectator + kt;
	emm.setMass(ZERO);
	emm.rescaleEnergy();

	// spe = ((1.-x)u/x)pEmitter + (1.-u)*pSpectator - kt;
	spe = (ptpt+mj2-sqr(1.-u)mj2)/((1.-u)sbar)pEmitter +
	(1.-u)*pSpectator - kt;
	spe.setMass(dInfo.spectatorData()->mass());
	spe.rescaleEnergy();

	}

	emitterMomentum(em);
	emissionMomentum(emm);
	spectatorMomentum(spe);

	}

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


	void IFMassiveKinematics::persistentOutput(PersistentOStream & os) const {
	os << theCollinearScheme;
	}

	void IFMassiveKinematics::persistentInput(PersistentIStream & is, int) {
	is >> theCollinearScheme;
	}

	ClassDescription<IFMassiveKinematics> IFMassiveKinematics::initIFMassiveKinematics;
	// Definition of the static class description member.

	void IFMassiveKinematics::Init() {

	static ClassDocumentation<IFMassiveKinematics> documentation
	("IFMassiveKinematics implements massless splittings "
	"off a initial-final dipole.");


	static Switch<IFMassiveKinematics,bool> interfaceCollinearScheme
	("CollinearScheme",
	"[experimental] Switch on or off the collinear scheme",
	&IFMassiveKinematics::theCollinearScheme, false, false, false);
	static SwitchOption interfaceCollinearSchemeOn
	(interfaceCollinearScheme,
	"On",
	"Switch on the collinear scheme.",
	true);
	static SwitchOption interfaceCollinearSchemeOff
	(interfaceCollinearScheme,
	"Off",
	"Switch off the collinear scheme",
	false);

	interfaceCollinearScheme.rank(-1);

	}

	diff --git a/DipoleShower/Kinematics/IFMassiveKinematics.h b/DipoleShower/Kinematics/IFMassiveKinematics.h
	--- a/DipoleShower/Kinematics/IFMassiveKinematics.h
	+++ b/DipoleShower/Kinematics/IFMassiveKinematics.h
	@@ -1,246 +1,244 @@
	// -- C++ --
	//
	// IFLightKinematics.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_IFLightKinematics_H
	#define HERWIG_IFLightKinematics_H
	//
	// This is the declaration of the IFLightKinematics class.
	//

	#include "DipoleSplittingKinematics.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* \ingroup DipoleShower
	* \author Simon Platzer, Martin Stoll
	*
	* \brief IFMassiveKinematics implements massless splittings
	* off an initial-final dipole.
	*
	* @see \ref IFMassiveKinematicsInterfaces "The interfaces"
	* defined for IFMassiveKinematics.
	*/
	class IFMassiveKinematics: public DipoleSplittingKinematics {

	public:

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

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

	public:

	/**
	* Return the boundaries in between the evolution
	* variable random number is to be sampled; the lower
	* cuoff is assumed to correspond to the infrared cutoff.
	*/
	virtual pair<double,double> kappaSupport(const DipoleSplittingInfo& dIndex) const;

	/**
	* Return the boundaries in between the momentum
	* fraction random number is to be sampled.
	*/
	virtual pair<double,double> xiSupport(const DipoleSplittingInfo& dIndex) const;

	/**
	+ * Return the boundaries on the momentum fraction
	+ */
	+ virtual pair<double,double> zBoundaries(Energy,
	+ const DipoleSplittingInfo&,
	+ const DipoleSplittingKernel&) const {
	+ return make_pair(0.0,1.0);
	+ }
	+
	+ /**
	* Return the dipole scale associated to the
	* given pair of emitter and spectator. This
	* should be the invariant mass or absolute value
	* final/final or initial/initial and the absolute
	* value of the momentum transfer for intial/final or
	* final/initial dipoles.
	*/
	virtual Energy dipoleScale(const Lorentz5Momentum& pEmitter,
	const Lorentz5Momentum& pSpectator) const;

	/**
	* Return the maximum pt for the given dipole scale.
	*/
	virtual Energy ptMax(Energy dScale,
	double emX, double specX,
	const DipoleIndex&,
	const DipoleSplittingKernel&) const;

	/**
	* Return the maximum virtuality for the given dipole scale.
	*/
	virtual Energy QMax(Energy dScale,
	double emX, double specX,
	- const DipoleIndex& dIndex) const;
	+ const DipoleIndex& dIndex,
	+ const DipoleSplittingKernel&) const;

	/**
	* Return the pt given a virtuality.
	*/
	virtual Energy PtFromQ(Energy scale, const DipoleSplittingInfo&) const;

	/**
	* Return the virtuality given a pt.
	*/
	virtual Energy QFromPt(Energy scale, const DipoleSplittingInfo&) const;

	/**
	* Return the random number associated to
	* the given pt.
	*/
	virtual double ptToRandom(Energy pt, Energy dScale,
	- const DipoleIndex& dIndex) const;
	+ double emX, double specX,
	+ const DipoleIndex& dIndex,
	+ const DipoleSplittingKernel&) const;

	/**
	* Generate splitting variables given three random numbers
	* and the momentum fractions of the emitter and spectator.
	* Return true on success.
	*/
	virtual bool generateSplitting(double kappa, double xi, double phi,
	- DipoleSplittingInfo& dIndex);
	-
	- /**
	- * For the splitting products present in the given dipole splitting
	- * info object calculate the kinematics parameters and return the
	- * propagator factor.
	- */
	- virtual InvEnergy2 setKinematics(DipoleSplittingInfo&) const;
	-
	- /**
	- * For the splitting parameters given in the dipole splitting info
	- * object, calculate the phasespace Jacobian times the propagator
	- * factor.
	- */
	- virtual double jacobianTimesPropagator(const DipoleSplittingInfo&,
	- Energy) const;
	+ DipoleSplittingInfo& dIndex,
	+ const DipoleSplittingKernel&);

	/**
	* Generate the full kinematics given emitter and
	* spectator momentum and a previously completeted
	* DipoleSplittingInfo object.
	*/
	virtual void generateKinematics(const Lorentz5Momentum& pEmitter,
	const Lorentz5Momentum& pSpectator,
	const DipoleSplittingInfo& dInfo);

	public:

	/** @name Functions used by the persistent I/O system. */
	//@{
	/**
	* Function used to write out object persistently.
	* @param os the persistent output stream written to.
	*/
	void persistentOutput(PersistentOStream & os) const;

	/**
	* Function used to read in object persistently.
	* @param is the persistent input stream read from.
	* @param version the version number of the object when written.
	*/
	void persistentInput(PersistentIStream & is, int version);
	//@}

	/**
	* The standard Init function used to initialize the interfaces.
	* Called exactly once for each class by the class description system
	* before the main function starts or
	* when this class is dynamically loaded.
	*/
	static void Init();

	protected:

	/** @name Clone Methods. */
	//@{
	/**
	* Make a simple clone of this object.
	* @return a pointer to the new object.
	*/
	virtual IBPtr clone() const;

	/** Make a clone of this object, possibly modifying the cloned object
	* to make it sane.
	* @return a pointer to the new object.
	*/
	virtual IBPtr fullclone() const;
	//@}


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


	private:

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

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

	private:

	/**
	* Wether or not to choose the `collinear' scheme
	*/
	bool theCollinearScheme;

	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

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

	#include "IILightKinematics.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Switch.h"

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

	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "Herwig++/DipoleShower/Base/DipoleSplittingInfo.h"

	using namespace Herwig;

	IILightKinematics::IILightKinematics()
	: DipoleSplittingKinematics(), theCollinearScheme(false), didCollinear(false) {}

	IILightKinematics::~IILightKinematics() {}

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

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

	-pair<double,double> IILightKinematics::kappaSupport(const DipoleSplittingInfo&) const {
	- return make_pair(0.0,1.0);
	-}
	-
	-pair<double,double> IILightKinematics::xiSupport(const DipoleSplittingInfo& split) const {
	-
	- double c = sqrt(1.-4.*sqr(IRCutoff()/generator()->maximumCMEnergy()));
	-
	- if ( split.index().emitterData()->id() == ParticleID::g ) {
	- if ( split.emitterData()->id() == ParticleID::g ) {
	- double b = log((1.+c)/(1.-c));
	- return make_pair(-b,b);
	- } else {
	- return make_pair(log(0.5(1.-c)),log(0.5(1.+c)));
	- }
	- }
	-
	- if ( split.index().emitterData()->id() != ParticleID::g &&
	- split.emitterData()->id() != ParticleID::g ) {
	- return make_pair(-log(0.5(1.+c)),-log(0.5(1.-c)));
	- }
	-
	- return make_pair(0.5(1.-c),0.5(1.+c));
	-
	-}
	-
	-Energy IILightKinematics::dipoleScale(const Lorentz5Momentum& pEmitter,
	- const Lorentz5Momentum& pSpectator) const {
	- return sqrt(2.(pEmitterpSpectator));
	-}
	-
	Energy IILightKinematics::ptMax(Energy dScale,
	double emX, double specX,
	const DipoleIndex&,
	const DipoleSplittingKernel&) const {
	double tau = emX*specX;
	return (1.-tau) * dScale / (2.*sqrt(tau));
	}

	Energy IILightKinematics::QMax(Energy,
	double, double,
	- const DipoleIndex&) const {
	+ const DipoleIndex&,
	+ const DipoleSplittingKernel&) const {
	assert(false && "add this");
	return 0.0*GeV;
	}

	-Energy IILightKinematics::PtFromQ(Energy, const DipoleSplittingInfo&) const {
	- assert(false && "add this");
	- return 0.0*GeV;
	+Energy IILightKinematics::PtFromQ(Energy scale, const DipoleSplittingInfo& split) const {
	+ double z = split.lastZ();
	+ return scale*sqrt(1.-z);
	}

	-Energy IILightKinematics::QFromPt(Energy, const DipoleSplittingInfo&) const {
	- assert(false && "add this");
	- return 0.0*GeV;
	+Energy IILightKinematics::QFromPt(Energy scale, const DipoleSplittingInfo& split) const {
	+ double z = split.lastZ();
	+ return scale/sqrt(1.-z);
	}

	-
	-double IILightKinematics::ptToRandom(Energy pt, Energy,
	- const DipoleIndex&) const {
	- return log(pt/IRCutoff()) / log(0.5 * generator()->maximumCMEnergy()/IRCutoff());
	+pair<double,double> IILightKinematics::zBoundaries(Energy pt,
	+ const DipoleSplittingInfo& dInfo,
	+ const DipoleSplittingKernel&) const {
	+ double x = dInfo.emitterX()*dInfo.spectatorX();
	+ double s = sqrt(1.-sqr(pt/dInfo.hardPt()));
	+ return make_pair(0.5(1.+x-(1.-x)s),0.5(1.+x+(1.-x)s));
	}

	bool IILightKinematics::generateSplitting(double kappa, double xi, double rphi,
	- DipoleSplittingInfo& info) {
	+ DipoleSplittingInfo& info,
	+ const DipoleSplittingKernel& split) {

	if ( info.emitterX() < xMin() \|\|
	info.spectatorX() < xMin() ) {
	jacobian(0.0);
	return false;
	}

	- Energy pt = IRCutoff() * pow(0.5 * generator()->maximumCMEnergy()/IRCutoff(),kappa);
	+ double weight = 1.0;

	- double r = sqr(info.hardPt()/info.scale());
	- if ( sqr(pt) > sqr(info.hardPt())r((2.+1./r)-2.*sqrt(1.+1./r)) ) {
	+ Energy pt = generatePt(kappa,info.scale(),
	+ info.emitterX(),info.spectatorX(),
	+ info.index(),split,
	+ weight);
	+
	+ if ( pt < IRCutoff() \|\| pt > info.hardPt() ) {
	jacobian(0.0);
	return false;
	}

	- double z = 0.;
	- double mapZJacobian = 0.;
	+ double z = 0.0;

	if ( info.index().emitterData()->id() == ParticleID::g ) {
	if ( info.emitterData()->id() == ParticleID::g ) {
	- z = exp(xi)/(1.+exp(xi));
	- mapZJacobian = z*(1.-z);
	+ z = generateZ(xi,pt,OneOverZOneMinusZ,
	+ info,split,weight);
	} else {
	- z = exp(xi);
	- mapZJacobian = z;
	+ z = generateZ(xi,pt,OneOverZ,
	+ info,split,weight);
	}
	}

	if ( info.index().emitterData()->id() != ParticleID::g ) {
	if ( info.emitterData()->id() != ParticleID::g ) {
	- z = 1.-exp(-xi);
	- mapZJacobian = (1.-z);
	+ z = generateZ(xi,pt,OneOverOneMinusZ,
	+ info,split,weight);
	} else {
	- z = xi;
	- mapZJacobian = 1.;
	+ z = generateZ(xi,pt,FlatZ,
	+ info,split,weight);
	}
	}

	double ratio = sqr(pt/info.scale());

	- double x = ( z*(1.-z) - ratio ) / ( 1. - z );
	- double v = ratio / (1.-z);
	+ double x = z*(1.-z)/(1.-z+ratio);
	+ double v = ratio*z /(1.-z+ratio);

	- if ( x < 0. \|\| x > 1. \|\| v > 1. \|\| v > 1.-x ) {
	- jacobian(0.0);
	- return false;
	- }
	-
	- double tau = info.emitterX()*info.spectatorX();
	-
	- double zpx = 0.5*( 1.+ tau +
	- (1.-tau)sqrt(1.-sqr(2.pt/((1.-tau)*info.scale()))) );
	- double zmx = 0.5*( 1.+ tau -
	- (1.-tau)sqrt(1.-sqr(2.pt/((1.-tau)*info.scale()))) );
	-
	- double xq = sqr(pt/info.hardPt());
	-
	- double zpq = 0.5*( 1.+ xq +
	- (1.-xq)sqrt(1.-sqr(2.pt/((1.-xq)*info.scale()))) );
	- double zmq = 0.5*( 1.+ xq -
	- (1.-xq)sqrt(1.-sqr(2.pt/((1.-xq)*info.scale()))) );
	-
	- double zp = min(zpx,zpq);
	- double zm = max(zmx,zmq);
	-
	- if ( pt < IRCutoff() \|\|
	- pt > info.hardPt() \|\|
	- z > zp \|\| z < zm ) {
	+ if ( x < 0. \|\| x > 1. \|\| v < 0. \|\| v > 1.-x ) {
	jacobian(0.0);
	return false;
	}

	if ( !theCollinearScheme &&
	(1.-v-x)/(v+x) < 1. ) {
	if ( (x+v) < info.emitterX() \|\|
	x/(x+v) < info.spectatorX() ) {
	jacobian(0.0);
	return false;
	}
	} else {
	if ( x < info.emitterX() ) {
	jacobian(0.0);
	return false;
	}
	}

	double phi = 2.Constants::pirphi;

	- jacobian(2. * mapZJacobian * (1.-z)/(z(1.-z)-ratio) log(0.5 * generator()->maximumCMEnergy()/IRCutoff()));
	+ jacobian(weight*(1./z));

	lastPt(pt);
	lastZ(z);
	lastPhi(phi);

	if ( !theCollinearScheme &&
	(1.-v-x)/(v+x) < 1. ) {
	lastEmitterZ(x+v);
	lastSpectatorZ(x/(x+v));
	} else {
	lastEmitterZ(x);
	lastSpectatorZ(1.);
	}

	if ( theMCCheck )
	theMCCheck->book(info.emitterX(),info.spectatorX(),info.scale(),info.hardPt(),pt,z,jacobian());

	return true;

	}

	-InvEnergy2 IILightKinematics::setKinematics(DipoleSplittingInfo& split) const {
	-
	- Lorentz5Momentum emitter = split.splitEmitter()->momentum();
	- Lorentz5Momentum emission = split.emission()->momentum();
	- Lorentz5Momentum spectator = split.splitSpectator()->momentum();
	-
	- split.splittingKinematics(const_cast<IILightKinematics*>(this));
	-
	- Energy2 scale = 2.(-emissionemitter - emissionspectator + emitterspectator);
	- split.scale(sqrt(scale));
	-
	- double x = scale/(2.(emitterspectator));
	- double v = (emitteremission)/(emitterspectator);
	-
	- split.lastZ(v+x);
	- split.lastPt(split.scale() * sqrt(v*(1.-x-v)));
	-
	- split.hardPt(split.lastPt());
	-
	- if ( split.hardPt() > IRCutoff() ) {
	- split.continuesEvolving();
	- } else {
	- split.didStopEvolving();
	- }
	-
	- return 1./(2.x(emitter*emission));
	-
	-}
	-
	-double IILightKinematics::
	-jacobianTimesPropagator(const DipoleSplittingInfo&,
	- Energy) const {
	- assert(false && "implementation missing");
	- return 0.;
	-}
	-
	void IILightKinematics::generateKinematics(const Lorentz5Momentum& pEmitter,
	const Lorentz5Momentum& pSpectator,
	const DipoleSplittingInfo& dInfo) {

	Energy pt = dInfo.lastPt();
	double z = dInfo.lastZ();

	- double ratio = sqr(pt/(pEmitter+pSpectator).m());
	+ double ratio = sqr(pt)/(2.pEmitterpSpectator);

	- double x = ( z*(1.-z) - ratio ) / ( 1. - z );
	- double v = ratio / (1.-z);
	-
	- pt = sqrt(v(1.-x-v)/x) (pEmitter+pSpectator).m();
	+ double x = z*(1.-z)/(1.-z+ratio);
	+ double v = ratio*z /(1.-z+ratio);

	Lorentz5Momentum kt =
	getKt (pEmitter, pSpectator, pt, dInfo.lastPhi());

	if ( !theCollinearScheme &&
	(1.-v-x)/(v+x) < 1. ) {

	Lorentz5Momentum em =
	(1./(v+x))pEmitter+(v(1.-v-x)/(x(x+v)))pSpectator+kt/(x+v);
	em.setMass(0.*GeV);
	em.rescaleEnergy();

	Lorentz5Momentum emm =
	((1.-v-x)/(v+x))pEmitter+(v/(x(x+v)))*pSpectator+kt/(x+v);
	emm.setMass(0.*GeV);
	emm.rescaleEnergy();

	Lorentz5Momentum spe =
	(1.+v/x)*pSpectator;
	spe.setMass(0.*GeV);
	spe.rescaleEnergy();

	emitterMomentum(em);
	emissionMomentum(emm);
	spectatorMomentum(spe);

	didCollinear = false;

	} else {

	Lorentz5Momentum em =
	(1./x)*pEmitter;
	em.setMass(0.*GeV);
	em.rescaleEnergy();

	Lorentz5Momentum emm =
	((1.-x-v)/x)pEmitter+vpSpectator+kt;
	emm.setMass(0.*GeV);
	emm.rescaleEnergy();

	Lorentz5Momentum spe =
	pSpectator;

	emitterMomentum(em);
	emissionMomentum(emm);
	spectatorMomentum(spe);

	K = em + spe - emm;
	K2 = K.m2();

	Ktilde = pEmitter + pSpectator;
	KplusKtilde = K + Ktilde;

	KplusKtilde2 = KplusKtilde.m2();

	didCollinear = true;

	}

	}

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


	void IILightKinematics::persistentOutput(PersistentOStream & os) const {
	os << theCollinearScheme;
	}

	void IILightKinematics::persistentInput(PersistentIStream & is, int) {
	is >> theCollinearScheme;
	}

	ClassDescription<IILightKinematics> IILightKinematics::initIILightKinematics;
	// Definition of the static class description member.

	void IILightKinematics::Init() {

	static ClassDocumentation<IILightKinematics> documentation
	("IILightKinematics implements massless splittings "
	"off an initial-initial dipole.");


	static Switch<IILightKinematics,bool> interfaceCollinearScheme
	("CollinearScheme",
	"[experimental] Switch on or off the collinear scheme",
	&IILightKinematics::theCollinearScheme, false, false, false);
	static SwitchOption interfaceCollinearSchemeOn
	(interfaceCollinearScheme,
	"On",
	"Switch on the collinear scheme.",
	true);
	static SwitchOption interfaceCollinearSchemeOff
	(interfaceCollinearScheme,
	"Off",
	"Switch off the collinear scheme",
	false);

	interfaceCollinearScheme.rank(-1);

	}

	diff --git a/DipoleShower/Kinematics/IILightKinematics.h b/DipoleShower/Kinematics/IILightKinematics.h
	--- a/DipoleShower/Kinematics/IILightKinematics.h
	+++ b/DipoleShower/Kinematics/IILightKinematics.h
	@@ -1,269 +1,232 @@
	// -- C++ --
	//
	// IILightKinematics.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_IILightKinematics_H
	#define HERWIG_IILightKinematics_H
	//
	// This is the declaration of the IILightKinematics class.
	//

	#include "DipoleSplittingKinematics.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* \ingroup DipoleShower
	* \author Simon Platzer
	*
	* \brief IILightKinematics implements massless splittings
	* off an initial-initial dipole.
	*
	* @see \ref IILightKinematicsInterfaces "The interfaces"
	* defined for IILightKinematics.
	*/
	class IILightKinematics: public DipoleSplittingKinematics {

	public:

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

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

	public:

	/**
	- * Return the boundaries in between the evolution
	- * variable random number is to be sampled; the lower
	- * cuoff is assumed to correspond to the infrared cutoff.
	- */
	- virtual pair<double,double> kappaSupport(const DipoleSplittingInfo& dIndex) const;
	-
	- /**
	- * Return the boundaries in between the momentum
	- * fraction random number is to be sampled.
	- */
	- virtual pair<double,double> xiSupport(const DipoleSplittingInfo& dIndex) const;
	-
	- /**
	- * Return the dipole scale associated to the
	- * given pair of emitter and spectator. This
	- * should be the invariant mass or absolute value
	- * final/final or initial/initial and the absolute
	- * value of the momentum transfer for intial/final or
	- * final/initial dipoles.
	- */
	- virtual Energy dipoleScale(const Lorentz5Momentum& pEmitter,
	- const Lorentz5Momentum& pSpectator) const;
	-
	- /**
	* Return the maximum pt for the given dipole scale.
	*/
	virtual Energy ptMax(Energy dScale,
	double emX, double specX,
	const DipoleIndex&,
	const DipoleSplittingKernel&) const;

	/**
	* Return the maximum virtuality for the given dipole scale.
	*/
	virtual Energy QMax(Energy dScale,
	double emX, double specX,
	- const DipoleIndex& dIndex) const;
	+ const DipoleIndex& dIndex,
	+ const DipoleSplittingKernel&) const;

	/**
	* Return the pt given a virtuality.
	*/
	virtual Energy PtFromQ(Energy scale, const DipoleSplittingInfo&) const;

	/**
	* Return the virtuality given a pt.
	*/
	virtual Energy QFromPt(Energy scale, const DipoleSplittingInfo&) const;

	/**
	- * Return the random number associated to
	- * the given pt.
	+ * Return the boundaries on the momentum fraction
	*/
	- virtual double ptToRandom(Energy pt, Energy dScale,
	- const DipoleIndex& dIndex) const;
	+ virtual pair<double,double> zBoundaries(Energy pt,
	+ const DipoleSplittingInfo& dInfo,
	+ const DipoleSplittingKernel& split) const;

	/**
	* Generate splitting variables given three random numbers
	* and the momentum fractions of the emitter and spectator.
	* Return true on success.
	*/
	virtual bool generateSplitting(double kappa, double xi, double phi,
	- DipoleSplittingInfo& dIndex);
	-
	- /**
	- * For the splitting products present in the given dipole splitting
	- * info object calculate the kinematics parameters and return the
	- * propagator factor.
	- */
	- virtual InvEnergy2 setKinematics(DipoleSplittingInfo&) const;
	-
	- /**
	- * For the splitting parameters given in the dipole splitting info
	- * object, calculate the phasespace Jacobian times the propagator
	- * factor.
	- */
	- virtual double jacobianTimesPropagator(const DipoleSplittingInfo&,
	- Energy) const;
	+ DipoleSplittingInfo& dIndex,
	+ const DipoleSplittingKernel&);

	/**
	* Generate the full kinematics given emitter and
	* spectator momentum and a previously completeted
	* DipoleSplittingInfo object.
	*/
	virtual void generateKinematics(const Lorentz5Momentum& pEmitter,
	const Lorentz5Momentum& pSpectator,
	const DipoleSplittingInfo& dInfo);

	/*
	* Return true, if there is a transformation which should
	* be applied to all other final state particles except the ones
	* involved in the splitting after having performed the splitting
	*/
	virtual bool doesTransform () const { return theCollinearScheme \|\| didCollinear; }

	/*
	* perform the transformation, if existing
	*/
	virtual Lorentz5Momentum transform (const Lorentz5Momentum& p) const {
	if ( !theCollinearScheme && !didCollinear ) return p;
	return p-(2.(KplusKtildep)/KplusKtilde2)KplusKtilde+(2.(Ktildep)/K2)K;
	}

	public:

	/** @name Functions used by the persistent I/O system. */
	//@{
	/**
	* Function used to write out object persistently.
	* @param os the persistent output stream written to.
	*/
	void persistentOutput(PersistentOStream & os) const;

	/**
	* Function used to read in object persistently.
	* @param is the persistent input stream read from.
	* @param version the version number of the object when written.
	*/
	void persistentInput(PersistentIStream & is, int version);
	//@}

	/**
	* The standard Init function used to initialize the interfaces.
	* Called exactly once for each class by the class description system
	* before the main function starts or
	* when this class is dynamically loaded.
	*/
	static void Init();

	protected:

	/** @name Clone Methods. */
	//@{
	/**
	* Make a simple clone of this object.
	* @return a pointer to the new object.
	*/
	virtual IBPtr clone() const;

	/** Make a clone of this object, possibly modifying the cloned object
	* to make it sane.
	* @return a pointer to the new object.
	*/
	virtual IBPtr fullclone() const;
	//@}


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


	private:

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

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

	private:

	/**
	* Wether or not to choose the `collinear' scheme
	*/
	bool theCollinearScheme;

	bool didCollinear;

	Lorentz5Momentum K;
	Energy2 K2;
	Lorentz5Momentum Ktilde;
	Lorentz5Momentum KplusKtilde;
	Energy2 KplusKtilde2;

	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

	/** This template specialization informs ThePEG about the name of
	* the IILightKinematics class and the shared object where it is defined. */
	template <>
	struct ClassTraits<Herwig::IILightKinematics>
	: public ClassTraitsBase<Herwig::IILightKinematics> {
	/** Return a platform-independent class name */
	static string className() { return "Herwig::IILightKinematics"; }
	/**
	* The name of a file containing the dynamic library where the class
	* IILightKinematics is implemented. It may also include several, space-separated,
	* libraries if the class IILightKinematics 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_IILightKinematics_H */
	diff --git a/MatrixElement/Matchbox/Matching/ShowerApproximation.cc b/MatrixElement/Matchbox/Matching/ShowerApproximation.cc
	--- a/MatrixElement/Matchbox/Matching/ShowerApproximation.cc
	+++ b/MatrixElement/Matchbox/Matching/ShowerApproximation.cc
	@@ -1,547 +1,556 @@
	// -- C++ --
	//
	// ShowerApproximation.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2012 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the ShowerApproximation class.
	//

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

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

	#include "Herwig++/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h"

	using namespace Herwig;

	ShowerApproximation::ShowerApproximation()
	: HandlerBase(), theBelowCutoff(false),
	theFFPtCut(1.0*GeV), theFFScreeningScale(ZERO),
	theFIPtCut(1.0*GeV), theFIScreeningScale(ZERO),
	theIIPtCut(1.0*GeV), theIIScreeningScale(ZERO),
	theRestrictPhasespace(true), theHardScaleFactor(1.0),
	theRenormalizationScaleFactor(1.0), theFactorizationScaleFactor(1.0),
	theExtrapolationX(0.65),
	theRealEmissionScaleInSubtraction(showerScale),
	theBornScaleInSubtraction(showerScale),
	theEmissionScaleInSubtraction(showerScale),
	theRealEmissionScaleInSplitting(showerScale),
	theBornScaleInSplitting(showerScale),
	theEmissionScaleInSplitting(showerScale),
	theRenormalizationScaleFreeze(1.*GeV),
	theFactorizationScaleFreeze(1.*GeV),
	theProfileScales(true),
	theProfileRho(0.3) {}

	ShowerApproximation::~ShowerApproximation() {}

	void ShowerApproximation::setDipole(Ptr<SubtractionDipole>::tcptr dip) { theDipole = dip; }

	Ptr<SubtractionDipole>::tcptr ShowerApproximation::dipole() const { return theDipole; }

	bool ShowerApproximation::isAboveCutoff() const {

	if ( dipole()->bornEmitter() > 1 &&
	dipole()->bornSpectator() > 1 ) {
	return dipole()->lastPt() >= ffPtCut();
	} else if ( ( dipole()->bornEmitter() > 1 &&
	dipole()->bornSpectator() < 2 ) \|\|
	( dipole()->bornEmitter() < 2 &&
	dipole()->bornSpectator() > 1 ) ) {
	return dipole()->lastPt() >= fiPtCut();
	} else {
	assert(dipole()->bornEmitter() < 2 &&
	dipole()->bornSpectator() < 2);
	return dipole()->lastPt() >= iiPtCut();
	}

	return true;

	}

	Energy ShowerApproximation::hardScale() const {
	+ if ( !bornCXComb()->mePartonData()[0]->coloured() &&
	+ !bornCXComb()->mePartonData()[1]->coloured() ) {
	+ Energy maxPt = (bornCXComb()->meMomenta()[0] + bornCXComb()->meMomenta()[1]).m();
	+ maxPt *= hardScaleFactor();
	+ return maxPt;
	+ }
	Energy maxPt = generator()->maximumCMEnergy();
	vector<Lorentz5Momentum>::const_iterator p =
	bornCXComb()->meMomenta().begin() + 2;
	cPDVector::const_iterator pp =
	bornCXComb()->mePartonData().begin() + 2;
	for ( ; p != bornCXComb()->meMomenta().end(); ++p, ++pp )
	if ( (**pp).coloured() )
	maxPt = min(maxPt,p->perp());
	if ( maxPt == generator()->maximumCMEnergy() )
	maxPt = (bornCXComb()->meMomenta()[0] + bornCXComb()->meMomenta()[1]).m();
	maxPt *= hardScaleFactor();
	return maxPt;
	}

	double ShowerApproximation::hardScaleProfile(Energy hard, Energy soft) const {
	+ if ( !bornCXComb()->mePartonData()[0]->coloured() &&
	+ !bornCXComb()->mePartonData()[1]->coloured() )
	+ return 1;
	double x = soft/hard;
	if ( theProfileScales ) {
	if ( x > 1. ) {
	return 0.;
	} else if ( x <= 1. && x > 1. - theProfileRho ) {
	return sqr(1.-x)/(2.*sqr(theProfileRho));
	} else if ( x <= 1. - theProfileRho &&
	x > 1. - 2.*theProfileRho ) {
	return 1. - sqr(1.-2.theProfileRho-x)/(2.sqr(theProfileRho));
	} else {
	return 1.;
	}
	}
	if ( x <= 1. )
	return 1.;
	return 0.;
	}

	bool ShowerApproximation::isInShowerPhasespace() const {

	if ( !isAboveCutoff() )
	return false;
	if ( !restrictPhasespace() )
	return true;

	return dipole()->lastPt() <= hardScale();

	}

	Energy2 ShowerApproximation::showerEmissionScale() const {

	Energy2 mur = sqr(dipole()->lastPt());

	if ( dipole()->bornEmitter() > 1 &&
	dipole()->bornSpectator() > 1 ) {
	return mur + sqr(ffScreeningScale());
	} else if ( ( dipole()->bornEmitter() > 1 &&
	dipole()->bornSpectator() < 2 ) \|\|
	( dipole()->bornEmitter() < 2 &&
	dipole()->bornSpectator() > 1 ) ) {
	return mur + sqr(fiScreeningScale());
	} else {
	assert(dipole()->bornEmitter() < 2 &&
	dipole()->bornSpectator() < 2);
	return mur + sqr(iiScreeningScale());
	}

	return mur;

	}

	Energy2 ShowerApproximation::bornRenormalizationScale() const {
	return
	sqr(dipole()->underlyingBornME()->renormalizationScaleFactor()) *
	dipole()->underlyingBornME()->renormalizationScale();
	}

	Energy2 ShowerApproximation::bornFactorizationScale() const {
	return
	sqr(dipole()->underlyingBornME()->factorizationScaleFactor()) *
	dipole()->underlyingBornME()->factorizationScale();
	}

	Energy2 ShowerApproximation::realRenormalizationScale() const {
	return
	sqr(dipole()->realEmissionME()->renormalizationScaleFactor()) *
	dipole()->realEmissionME()->renormalizationScale();
	}

	Energy2 ShowerApproximation::realFactorizationScale() const {
	return
	sqr(dipole()->realEmissionME()->factorizationScaleFactor()) *
	dipole()->realEmissionME()->factorizationScale();
	}

	double ShowerApproximation::bornPDFWeight(Energy2 muf) const {
	if ( !bornCXComb()->mePartonData()[0]->coloured() &&
	!bornCXComb()->mePartonData()[1]->coloured() )
	return 1.;
	if ( muf < sqr(theFactorizationScaleFreeze) )
	muf = sqr(theFactorizationScaleFreeze);
	double pdfweight = 1.;
	if ( bornCXComb()->mePartonData()[0]->coloured() &&
	dipole()->underlyingBornME()->havePDFWeight1() )
	pdfweight *= dipole()->underlyingBornME()->pdf1(muf,theExtrapolationX);
	if ( bornCXComb()->mePartonData()[1]->coloured() &&
	dipole()->underlyingBornME()->havePDFWeight2() )
	pdfweight *= dipole()->underlyingBornME()->pdf2(muf,theExtrapolationX);
	return pdfweight;
	}

	double ShowerApproximation::realPDFWeight(Energy2 muf) const {
	if ( !realCXComb()->mePartonData()[0]->coloured() &&
	!realCXComb()->mePartonData()[1]->coloured() )
	return 1.;
	if ( muf < sqr(theFactorizationScaleFreeze) )
	muf = sqr(theFactorizationScaleFreeze);
	double pdfweight = 1.;
	if ( realCXComb()->mePartonData()[0]->coloured() &&
	dipole()->realEmissionME()->havePDFWeight1() )
	pdfweight *= dipole()->realEmissionME()->pdf1(muf,theExtrapolationX);
	if ( realCXComb()->mePartonData()[1]->coloured() &&
	dipole()->realEmissionME()->havePDFWeight2() )
	pdfweight *= dipole()->realEmissionME()->pdf2(muf,theExtrapolationX);
	return pdfweight;
	}

	double ShowerApproximation::scaleWeight(int rScale, int bScale, int eScale) const {

	double emissionAlpha = 1.;
	Energy2 emissionScale = ZERO;
	Energy2 showerscale = ZERO;
	if ( eScale == showerScale \|\| bScale == showerScale \|\| eScale == showerScale ) {
	showerscale = showerRenormalizationScale();
	if ( showerscale < sqr(theRenormalizationScaleFreeze) )
	showerscale = sqr(theFactorizationScaleFreeze);
	}
	if ( eScale == showerScale ) {
	emissionAlpha = SM().alphaS(showerscale);
	emissionScale = showerFactorizationScale();
	} else if ( eScale == realScale ) {
	emissionAlpha = dipole()->realEmissionME()->lastXComb().lastAlphaS();
	emissionScale = dipole()->realEmissionME()->lastScale();
	} else if ( eScale == bornScale ) {
	emissionAlpha = dipole()->underlyingBornME()->lastXComb().lastAlphaS();
	emissionScale = dipole()->underlyingBornME()->lastScale();
	}
	double emissionPDF = realPDFWeight(emissionScale);

	double couplingFactor = 1.;
	if ( bScale != rScale ) {
	double bornAlpha = 1.;
	if ( bScale == showerScale ) {
	bornAlpha = SM().alphaS(showerscale);
	} else if ( bScale == realScale ) {
	bornAlpha = dipole()->realEmissionME()->lastXComb().lastAlphaS();
	} else if ( bScale == bornScale ) {
	bornAlpha = dipole()->underlyingBornME()->lastXComb().lastAlphaS();
	}
	double realAlpha = 1.;
	if ( rScale == showerScale ) {
	realAlpha = SM().alphaS(showerscale);
	} else if ( rScale == realScale ) {
	realAlpha = dipole()->realEmissionME()->lastXComb().lastAlphaS();
	} else if ( rScale == bornScale ) {
	realAlpha = dipole()->underlyingBornME()->lastXComb().lastAlphaS();
	}
	couplingFactor *=
	pow(realAlpha/bornAlpha,(double)(dipole()->underlyingBornME()->orderInAlphaS()));
	}

	Energy2 hardScale = ZERO;
	if ( bScale == showerScale ) {
	hardScale = showerFactorizationScale();
	} else if ( bScale == realScale ) {
	hardScale = dipole()->realEmissionME()->lastScale();
	} else if ( bScale == bornScale ) {
	hardScale = dipole()->underlyingBornME()->lastScale();
	}
	double bornPDF = bornPDFWeight(hardScale);

	if ( bornPDF < 1e-8 )
	bornPDF = 0.0;

	if ( emissionPDF < 1e-8 )
	emissionPDF = 0.0;

	if ( emissionPDF == 0.0 \|\| bornPDF == 0.0 )
	return 0.0;

	return
	emissionAlpha * emissionPDF *
	couplingFactor / bornPDF;

	}

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


	void ShowerApproximation::persistentOutput(PersistentOStream & os) const {
	os << theBornXComb << theRealXComb << theTildeXCombs << theDipole << theBelowCutoff
	<< ounit(theFFPtCut,GeV) << ounit(theFFScreeningScale,GeV)
	<< ounit(theFIPtCut,GeV) << ounit(theFIScreeningScale,GeV)
	<< ounit(theIIPtCut,GeV) << ounit(theIIScreeningScale,GeV)
	<< theRestrictPhasespace << theHardScaleFactor
	<< theRenormalizationScaleFactor << theFactorizationScaleFactor
	<< theExtrapolationX
	<< theRealEmissionScaleInSubtraction << theBornScaleInSubtraction
	<< theEmissionScaleInSubtraction << theRealEmissionScaleInSplitting
	<< theBornScaleInSplitting << theEmissionScaleInSplitting
	<< ounit(theRenormalizationScaleFreeze,GeV)
	<< ounit(theFactorizationScaleFreeze,GeV)
	<< theProfileScales << theProfileRho;
	}

	void ShowerApproximation::persistentInput(PersistentIStream & is, int) {
	is >> theBornXComb >> theRealXComb >> theTildeXCombs >> theDipole >> theBelowCutoff
	>> iunit(theFFPtCut,GeV) >> iunit(theFFScreeningScale,GeV)
	>> iunit(theFIPtCut,GeV) >> iunit(theFIScreeningScale,GeV)
	>> iunit(theIIPtCut,GeV) >> iunit(theIIScreeningScale,GeV)
	>> theRestrictPhasespace >> theHardScaleFactor
	>> theRenormalizationScaleFactor >> theFactorizationScaleFactor
	>> theExtrapolationX
	>> theRealEmissionScaleInSubtraction >> theBornScaleInSubtraction
	>> theEmissionScaleInSubtraction >> theRealEmissionScaleInSplitting
	>> theBornScaleInSplitting >> theEmissionScaleInSplitting
	>> iunit(theRenormalizationScaleFreeze,GeV)
	>> iunit(theFactorizationScaleFreeze,GeV)
	>> theProfileScales >> theProfileRho;
	}


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

	void ShowerApproximation::Init() {

	static ClassDocumentation<ShowerApproximation> documentation
	("ShowerApproximation describes the shower emission to be used "
	"in NLO matching.");

	static Parameter<ShowerApproximation,Energy> interfaceFFPtCut
	("FFPtCut",
	"Set the pt infrared cutoff",
	&ShowerApproximation::theFFPtCut, GeV, 1.0GeV, 0.0GeV, 0*GeV,
	false, false, Interface::lowerlim);

	static Parameter<ShowerApproximation,Energy> interfaceFIPtCut
	("FIPtCut",
	"Set the pt infrared cutoff",
	&ShowerApproximation::theFIPtCut, GeV, 1.0GeV, 0.0GeV, 0*GeV,
	false, false, Interface::lowerlim);

	static Parameter<ShowerApproximation,Energy> interfaceIIPtCut
	("IIPtCut",
	"Set the pt infrared cutoff",
	&ShowerApproximation::theIIPtCut, GeV, 1.0GeV, 0.0GeV, 0*GeV,
	false, false, Interface::lowerlim);

	static Parameter<ShowerApproximation,Energy> interfaceFFScreeningScale
	("FFScreeningScale",
	"Set the screening scale",
	&ShowerApproximation::theFFScreeningScale, GeV, 0.0GeV, 0.0GeV, 0*GeV,
	false, false, Interface::lowerlim);

	static Parameter<ShowerApproximation,Energy> interfaceFIScreeningScale
	("FIScreeningScale",
	"Set the screening scale",
	&ShowerApproximation::theFIScreeningScale, GeV, 0.0GeV, 0.0GeV, 0*GeV,
	false, false, Interface::lowerlim);

	static Parameter<ShowerApproximation,Energy> interfaceIIScreeningScale
	("IIScreeningScale",
	"Set the screening scale",
	&ShowerApproximation::theIIScreeningScale, GeV, 0.0GeV, 0.0GeV, 0*GeV,
	false, false, Interface::lowerlim);

	static Switch<ShowerApproximation,bool> interfaceRestrictPhasespace
	("RestrictPhasespace",
	"Switch on or off phasespace restrictions",
	&ShowerApproximation::theRestrictPhasespace, true, false, false);
	static SwitchOption interfaceRestrictPhasespaceOn
	(interfaceRestrictPhasespace,
	"On",
	"Perform phasespace restrictions",
	true);
	static SwitchOption interfaceRestrictPhasespaceOff
	(interfaceRestrictPhasespace,
	"Off",
	"Do not perform phasespace restrictions",
	false);

	static Parameter<ShowerApproximation,double> interfaceHardScaleFactor
	("HardScaleFactor",
	"The hard scale factor.",
	&ShowerApproximation::theHardScaleFactor, 1.0, 0.0, 0,
	false, false, Interface::lowerlim);

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

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

	static Parameter<ShowerApproximation,double> interfaceExtrapolationX
	("ExtrapolationX",
	"The x from which on extrapolation should be performed.",
	&ShowerApproximation::theExtrapolationX, 0.65, 0.0, 1.0,
	false, false, Interface::limited);

	static Switch<ShowerApproximation,int> interfaceRealEmissionScaleInSubtraction
	("RealEmissionScaleInSubtraction",
	"Set the scale choice for the real emission cross section in the matching subtraction.",
	&ShowerApproximation::theRealEmissionScaleInSubtraction, showerScale, false, false);
	static SwitchOption interfaceRealEmissionScaleInSubtractionRealScale
	(interfaceRealEmissionScaleInSubtraction,
	"RealScale",
	"Use the real emission scale.",
	realScale);
	static SwitchOption interfaceRealEmissionScaleInSubtractionBornScale
	(interfaceRealEmissionScaleInSubtraction,
	"BornScale",
	"Use the Born scale.",
	bornScale);
	static SwitchOption interfaceRealEmissionScaleInSubtractionShowerScale
	(interfaceRealEmissionScaleInSubtraction,
	"ShowerScale",
	"Use the shower scale",
	showerScale);

	static Switch<ShowerApproximation,int> interfaceBornScaleInSubtraction
	("BornScaleInSubtraction",
	"Set the scale choice for the Born cross section in the matching subtraction.",
	&ShowerApproximation::theBornScaleInSubtraction, showerScale, false, false);
	static SwitchOption interfaceBornScaleInSubtractionRealScale
	(interfaceBornScaleInSubtraction,
	"RealScale",
	"Use the real emission scale.",
	realScale);
	static SwitchOption interfaceBornScaleInSubtractionBornScale
	(interfaceBornScaleInSubtraction,
	"BornScale",
	"Use the Born scale.",
	bornScale);
	static SwitchOption interfaceBornScaleInSubtractionShowerScale
	(interfaceBornScaleInSubtraction,
	"ShowerScale",
	"Use the shower scale",
	showerScale);

	static Switch<ShowerApproximation,int> interfaceEmissionScaleInSubtraction
	("EmissionScaleInSubtraction",
	"Set the scale choice for the emission in the matching subtraction.",
	&ShowerApproximation::theEmissionScaleInSubtraction, showerScale, false, false);
	static SwitchOption interfaceEmissionScaleInSubtractionRealScale
	(interfaceEmissionScaleInSubtraction,
	"RealScale",
	"Use the real emission scale.",
	realScale);
	static SwitchOption interfaceEmissionScaleInSubtractionEmissionScale
	(interfaceEmissionScaleInSubtraction,
	"BornScale",
	"Use the Born scale.",
	bornScale);
	static SwitchOption interfaceEmissionScaleInSubtractionShowerScale
	(interfaceEmissionScaleInSubtraction,
	"ShowerScale",
	"Use the shower scale",
	showerScale);

	static Switch<ShowerApproximation,int> interfaceRealEmissionScaleInSplitting
	("RealEmissionScaleInSplitting",
	"Set the scale choice for the real emission cross section in the splitting.",
	&ShowerApproximation::theRealEmissionScaleInSplitting, showerScale, false, false);
	static SwitchOption interfaceRealEmissionScaleInSplittingRealScale
	(interfaceRealEmissionScaleInSplitting,
	"RealScale",
	"Use the real emission scale.",
	realScale);
	static SwitchOption interfaceRealEmissionScaleInSplittingBornScale
	(interfaceRealEmissionScaleInSplitting,
	"BornScale",
	"Use the Born scale.",
	bornScale);
	static SwitchOption interfaceRealEmissionScaleInSplittingShowerScale
	(interfaceRealEmissionScaleInSplitting,
	"ShowerScale",
	"Use the shower scale",
	showerScale);

	static Switch<ShowerApproximation,int> interfaceBornScaleInSplitting
	("BornScaleInSplitting",
	"Set the scale choice for the Born cross section in the splitting.",
	&ShowerApproximation::theBornScaleInSplitting, showerScale, false, false);
	static SwitchOption interfaceBornScaleInSplittingRealScale
	(interfaceBornScaleInSplitting,
	"RealScale",
	"Use the real emission scale.",
	realScale);
	static SwitchOption interfaceBornScaleInSplittingBornScale
	(interfaceBornScaleInSplitting,
	"BornScale",
	"Use the Born scale.",
	bornScale);
	static SwitchOption interfaceBornScaleInSplittingShowerScale
	(interfaceBornScaleInSplitting,
	"ShowerScale",
	"Use the shower scale",
	showerScale);

	static Switch<ShowerApproximation,int> interfaceEmissionScaleInSplitting
	("EmissionScaleInSplitting",
	"Set the scale choice for the emission in the splitting.",
	&ShowerApproximation::theEmissionScaleInSplitting, showerScale, false, false);
	static SwitchOption interfaceEmissionScaleInSplittingRealScale
	(interfaceEmissionScaleInSplitting,
	"RealScale",
	"Use the real emission scale.",
	realScale);
	static SwitchOption interfaceEmissionScaleInSplittingEmissionScale
	(interfaceEmissionScaleInSplitting,
	"BornScale",
	"Use the Born scale.",
	bornScale);
	static SwitchOption interfaceEmissionScaleInSplittingShowerScale
	(interfaceEmissionScaleInSplitting,
	"ShowerScale",
	"Use the shower scale",
	showerScale);

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

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

	static Switch<ShowerApproximation,bool> interfaceProfileScales
	("ProfileScales",
	"Switch on or off use of profile scales.",
	&ShowerApproximation::theProfileScales, true, false, false);
	static SwitchOption interfaceProfileScalesYes
	(interfaceProfileScales,
	"Yes",
	"Use profile scales.",
	true);
	static SwitchOption interfaceProfileScalesNo
	(interfaceProfileScales,
	"No",
	"Use a hard cutoff.",
	false);

	static Parameter<ShowerApproximation,double> interfaceProfileRho
	("ProfileRho",
	"The rho parameter of the profile scales.",
	&ShowerApproximation::theProfileRho, 0.3, 0.0, 1.0,
	false, false, Interface::limited);

	}

File Metadata

Mime Type: text/x-diff
Expires: Wed, May 14, 11:09 AM (17 h, 46 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 5092916
Default Alt Text: (222 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions