No OneTemporary
Actions

Size

84 KB

Subscribers

None

View Options

	diff --git a/DipoleShower/DipoleShowerHandler.cc b/DipoleShower/DipoleShowerHandler.cc
	--- a/DipoleShower/DipoleShowerHandler.cc
	+++ b/DipoleShower/DipoleShowerHandler.cc
	@@ -1,1098 +1,1103 @@
	// -- C++ --
	//
	// DipoleShowerHandler.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2007 The Herwig Collaboration
	//
	// Herwig is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the DipoleShowerHandler class.
	//

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

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

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

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

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

	using namespace Herwig;

	bool DipoleShowerHandler::firstWarn = true;

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

	DipoleShowerHandler::~DipoleShowerHandler() {}

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

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

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

	useMe();

	prepareCascade(sub);
	+ resetWeights();

	if ( !doFSR && ! doISR )
	return sub->incoming();

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

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

	nTries = 0;

	while ( true ) {

	try {

	didRadiate = false;
	didRealign = false;

	isMCatNLOSEvent = false;
	isMCatNLOHEvent = false;

	if ( eventRecord().xcombPtr() ) {

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


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

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

	}

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

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

	unsigned int nEmitted = 0;

	if ( firstMCatNLOEmission ) {

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

	}

	if ( !firstMCatNLOEmission ) {

	doCascade(nEmitted,optHardPt,optCutoff);

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

	} else {

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

	}

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

	didRealign = realign();

	constituentReshuffle();

	break;

	} catch (RedoShower&) {

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

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

	continue;

	} catch (...) {
	throw;
	}

	}

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

	}

	void DipoleShowerHandler::constituentReshuffle() {

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

	}

	void DipoleShowerHandler::hardScales(Energy2 muf) {

	maxPt = generator()->maximumCMEnergy();

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

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

	Energy minVetoScale = -1.*GeV;

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

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

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

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

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

	conf = make_pair(false,true);

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

	}

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

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

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

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

	}

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

	}

	}

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

	if ( maxPossible <= ircutoff ) {
	continue;
	}

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

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

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

	}

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

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

	return winnerScale;

	}

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

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

	DipoleSplittingInfo winner;
	DipoleSplittingInfo dipoleWinner;

	while ( eventRecord().haveChain() ) {

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

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

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

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

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

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

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

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

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

	}

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

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

	// otherwise perform the splitting

	didRadiate = true;

	isMCatNLOSEvent = false;
	isMCatNLOHEvent = false;

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

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

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

	assert(firstChain && secondChain);

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

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

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

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

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

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

	}

	}

	bool DipoleShowerHandler::realign() {

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

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

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

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

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

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

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

	transform.rotate(-theta,rot);

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

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

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

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

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

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

	Lorentz5Momentum tmp;

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

	}

	return false;

	}

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

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

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

	}

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

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

	bool gotone = false;

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

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

	gotone = true;

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

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

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

	equivalent = eq;

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

	break;

	}
	}

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

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

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


	}
	}

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

	}

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

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

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

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

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

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

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

	void DipoleShowerHandler::Init() {

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

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


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


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


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

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


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

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

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


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

	interfaceChainOrderVetoScales.rank(-1);


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

	interfaceNEmissions.rank(-1);


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

	interfaceDiscardNoEmissions.rank(-1);

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

	interfaceFirstMCatNLOEmission.rank(-1);


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

	interfaceVerbosity.rank(-1);


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

	interfacePrintEvent.rank(-1);

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

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

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

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

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

	}

	diff --git a/Shower/ShowerHandler.cc b/Shower/ShowerHandler.cc
	--- a/Shower/ShowerHandler.cc
	+++ b/Shower/ShowerHandler.cc
	@@ -1,900 +1,903 @@
	// -- C++ --
	//
	// ShowerHandler.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 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 ShowerHandler class.
	//

	#include "ShowerHandler.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/ParVector.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Interface/Command.h"
	#include "ThePEG/PDF/PartonExtractor.h"
	#include "ThePEG/PDF/PartonBinInstance.h"
	#include "Herwig/PDT/StandardMatchers.h"
	#include "ThePEG/Cuts/Cuts.h"
	#include "ThePEG/Handlers/StandardXComb.h"
	#include "ThePEG/Utilities/Throw.h"
	#include "ThePEG/Utilities/StringUtils.h"
	#include "Herwig/Shower/Base/Evolver.h"
	#include "Herwig/Shower/Base/ShowerParticle.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "Herwig/Utilities/EnumParticles.h"
	#include "Herwig/PDF/MPIPDF.h"
	#include "Herwig/PDF/MinBiasPDF.h"
	#include "ThePEG/Handlers/EventHandler.h"
	#include "Herwig/Shower/Base/ShowerTree.h"
	#include "Herwig/Shower/Base/HardTree.h"
	#include "Herwig/Shower/Base/KinematicsReconstructor.h"
	#include "Herwig/Shower/Base/PartnerFinder.h"
	#include "Herwig/PDF/HwRemDecayer.h"
	#include <cassert>
	#include "ThePEG/Utilities/DescribeClass.h"

	using namespace Herwig;

	DescribeClass<ShowerHandler,CascadeHandler>
	describeShowerHandler ("Herwig::ShowerHandler","HwShower.so");

	ShowerHandler::~ShowerHandler() {}

	ShowerHandler * ShowerHandler::currentHandler_ = 0;

	void ShowerHandler::doinit() {
	CascadeHandler::doinit();
	// copy particles to decay before showering from input vector to the
	// set used in the simulation
	if ( particlesDecayInShower_.empty() )
	particlesDecayInShower_.insert(inputparticlesDecayInShower_.begin(),
	inputparticlesDecayInShower_.end());
	ShowerTree::_decayInShower = particlesDecayInShower_;
	ShowerTree::_vmin2 = vMin_;
	ShowerTree::_spaceTime = includeSpaceTime_;
	}

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

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

	ShowerHandler::ShowerHandler() :
	pdfFreezingScale_(2.5*GeV),
	maxtry_(10),maxtryMPI_(10),maxtryDP_(10),
	includeSpaceTime_(false), vMin_(0.1*GeV2), subProcess_(),
	factorizationScaleFactor_(1.0),
	renormalizationScaleFactor_(1.0),
	hardScaleFactor_(1.0), scaleFactorOption_(0),
	restrictPhasespace_(true), maxPtIsMuF_(false),
	splitHardProcess_(true) {
	inputparticlesDecayInShower_.push_back( 6 ); // top
	inputparticlesDecayInShower_.push_back( 23 ); // Z0
	inputparticlesDecayInShower_.push_back( 24 ); // W+/-
	inputparticlesDecayInShower_.push_back( 25 ); // h0
	}

	void ShowerHandler::doinitrun(){
	CascadeHandler::doinitrun();
	//can't use isMPIOn here, because the EventHandler is not set at that stage
	if(MPIHandler_){
	MPIHandler_->initialize();
	if(MPIHandler_->softInt())
	remDec_->initSoftInteractions(MPIHandler_->Ptmin(), MPIHandler_->beta());
	}
	ShowerTree::_decayInShower = particlesDecayInShower_;
	ShowerTree::_vmin2 = vMin_;
	ShowerTree::_spaceTime = includeSpaceTime_;
	}

	void ShowerHandler::dofinish(){
	CascadeHandler::dofinish();
	if(MPIHandler_) MPIHandler_->finalize();
	}

	void ShowerHandler::persistentOutput(PersistentOStream & os) const {
	os << evolver_ << remDec_ << ounit(pdfFreezingScale_,GeV) << maxtry_
	<< maxtryMPI_ << maxtryDP_ << inputparticlesDecayInShower_
	<< particlesDecayInShower_ << MPIHandler_ << PDFA_ << PDFB_
	<< PDFARemnant_ << PDFBRemnant_
	<< includeSpaceTime_ << ounit(vMin_,GeV2)
	<< factorizationScaleFactor_ << renormalizationScaleFactor_
	<< hardScaleFactor_ << scaleFactorOption_
	<< restrictPhasespace_ << maxPtIsMuF_ << hardScaleProfile_
	<< splitHardProcess_ << showerVariations_;
	}

	void ShowerHandler::persistentInput(PersistentIStream & is, int) {
	is >> evolver_ >> remDec_ >> iunit(pdfFreezingScale_,GeV) >> maxtry_
	>> maxtryMPI_ >> maxtryDP_ >> inputparticlesDecayInShower_
	>> particlesDecayInShower_ >> MPIHandler_ >> PDFA_ >> PDFB_
	>> PDFARemnant_ >> PDFBRemnant_
	>> includeSpaceTime_ >> iunit(vMin_,GeV2)
	>> factorizationScaleFactor_ >> renormalizationScaleFactor_
	>> hardScaleFactor_ >> scaleFactorOption_
	>> restrictPhasespace_ >> maxPtIsMuF_ >> hardScaleProfile_
	>> splitHardProcess_ >> showerVariations_;
	}

	void ShowerHandler::Init() {

	static ClassDocumentation<ShowerHandler> documentation
	("Main driver class for the showering.",
	"The Shower evolution was performed using an algorithm described in "
	"\\cite{Marchesini:1983bm,Marchesini:1987cf,Gieseke:2003rz,Bahr:2008pv}.",
	"%\\cite{Marchesini:1983bm}\n"
	"\\bibitem{Marchesini:1983bm}\n"
	" G.~Marchesini and B.~R.~Webber,\n"
	" ``Simulation Of QCD Jets Including Soft Gluon Interference,''\n"
	" Nucl.\\ Phys.\\ B {\\bf 238}, 1 (1984).\n"
	" %%CITATION = NUPHA,B238,1;%%\n"
	"%\\cite{Marchesini:1987cf}\n"
	"\\bibitem{Marchesini:1987cf}\n"
	" G.~Marchesini and B.~R.~Webber,\n"
	" ``Monte Carlo Simulation of General Hard Processes with Coherent QCD\n"
	" Radiation,''\n"
	" Nucl.\\ Phys.\\ B {\\bf 310}, 461 (1988).\n"
	" %%CITATION = NUPHA,B310,461;%%\n"
	"%\\cite{Gieseke:2003rz}\n"
	"\\bibitem{Gieseke:2003rz}\n"
	" S.~Gieseke, P.~Stephens and B.~Webber,\n"
	" ``New formalism for QCD parton showers,''\n"
	" JHEP {\\bf 0312}, 045 (2003)\n"
	" [arXiv:hep-ph/0310083].\n"
	" %%CITATION = JHEPA,0312,045;%%\n"
	);

	static Reference<ShowerHandler,Evolver>
	interfaceEvolver("Evolver",
	"A reference to the Evolver object",
	&Herwig::ShowerHandler::evolver_,
	false, false, true, false);

	static Reference<ShowerHandler,HwRemDecayer>
	interfaceRemDecayer("RemDecayer",
	"A reference to the Remnant Decayer object",
	&Herwig::ShowerHandler::remDec_,
	false, false, true, false);

	static Parameter<ShowerHandler,Energy> interfacePDFFreezingScale
	("PDFFreezingScale",
	"The PDF freezing scale",
	&ShowerHandler::pdfFreezingScale_, GeV, 2.5GeV, 2.0GeV, 10.0*GeV,
	false, false, Interface::limited);

	static Parameter<ShowerHandler,unsigned int> interfaceMaxTry
	("MaxTry",
	"The maximum number of attempts for the main showering loop",
	&ShowerHandler::maxtry_, 10, 1, 100,
	false, false, Interface::limited);

	static Parameter<ShowerHandler,unsigned int> interfaceMaxTryMPI
	("MaxTryMPI",
	"The maximum number of regeneration attempts for an additional scattering",
	&ShowerHandler::maxtryMPI_, 10, 0, 100,
	false, false, Interface::limited);

	static Parameter<ShowerHandler,unsigned int> interfaceMaxTryDP
	("MaxTryDP",
	"The maximum number of regeneration attempts for an additional hard scattering",
	&ShowerHandler::maxtryDP_, 10, 0, 100,
	false, false, Interface::limited);

	static ParVector<ShowerHandler,long> interfaceDecayInShower
	("DecayInShower",
	"PDG codes of the particles to be decayed in the shower",
	&ShowerHandler::inputparticlesDecayInShower_, -1, 0l, -10000000l, 10000000l,
	false, false, Interface::limited);

	static Reference<ShowerHandler,UEBase> interfaceMPIHandler
	("MPIHandler",
	"The object that administers all additional scatterings.",
	&ShowerHandler::MPIHandler_, false, false, true, true);

	static Reference<ShowerHandler,PDFBase> interfacePDFA
	("PDFA",
	"The PDF for beam particle A. Overrides the particle's own PDF setting."
	"By default used for both the shower and forced splitting in the remnant",
	&ShowerHandler::PDFA_, false, false, true, true, false);

	static Reference<ShowerHandler,PDFBase> interfacePDFB
	("PDFB",
	"The PDF for beam particle B. Overrides the particle's own PDF setting."
	"By default used for both the shower and forced splitting in the remnant",
	&ShowerHandler::PDFB_, false, false, true, true, false);

	static Reference<ShowerHandler,PDFBase> interfacePDFARemnant
	("PDFARemnant",
	"The PDF for beam particle A used to generate forced splittings of the remnant."
	" This overrides both the particle's own PDF setting and the value set by PDFA if used.",
	&ShowerHandler::PDFARemnant_, false, false, true, true, false);

	static Reference<ShowerHandler,PDFBase> interfacePDFBRemnant
	("PDFBRemnant",
	"The PDF for beam particle B used to generate forced splittings of the remnant."
	" This overrides both the particle's own PDF setting and the value set by PDFB if used.",
	&ShowerHandler::PDFBRemnant_, false, false, true, true, false);


	static Switch<ShowerHandler,bool> interfaceIncludeSpaceTime
	("IncludeSpaceTime",
	"Whether to include the model for the calculation of space-time distances",
	&ShowerHandler::includeSpaceTime_, false, false, false);
	static SwitchOption interfaceIncludeSpaceTimeYes
	(interfaceIncludeSpaceTime,
	"Yes",
	"Include the model",
	true);
	static SwitchOption interfaceIncludeSpaceTimeNo
	(interfaceIncludeSpaceTime,
	"No",
	"Only include the displacement from the particle-s lifetime for decaying particles",
	false);

	static Parameter<ShowerHandler,Energy2> interfaceMinimumVirtuality
	("MinimumVirtuality",
	"The minimum virtuality for the space-time model",
	&ShowerHandler::vMin_, GeV2, 0.1GeV2, 0.0GeV2, 1000.0*GeV2,
	false, false, Interface::limited);

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

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

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

	static Switch<ShowerHandler,int> interfaceScaleFactorOption
	("ScaleFactorOption",
	"Where to apply scale factors.",
	&ShowerHandler::scaleFactorOption_, 0, false, false);
	static SwitchOption interfaceScaleFactorOptionAll
	(interfaceScaleFactorOption,
	"All",
	"Apply in first and secondary scatterings.",
	0);
	static SwitchOption interfaceScaleFactorOptionHard
	(interfaceScaleFactorOption,
	"Hard",
	"Only apply for the hard process.",
	1);
	static SwitchOption interfaceScaleFactorOptionSecondary
	(interfaceScaleFactorOption,
	"Secondary",
	"Only apply to secondary scatterings.",
	2);

	static Reference<ShowerHandler,HardScaleProfile> interfaceHardScaleProfile
	("HardScaleProfile",
	"The hard scale profile to use.",
	&ShowerHandler::hardScaleProfile_, false, false, true, true, false);

	static Switch<ShowerHandler,bool> interfaceMaxPtIsMuF
	("MaxPtIsMuF",
	"",
	&ShowerHandler::maxPtIsMuF_, false, false, false);
	static SwitchOption interfaceMaxPtIsMuFYes
	(interfaceMaxPtIsMuF,
	"Yes",
	"",
	true);
	static SwitchOption interfaceMaxPtIsMuFNo
	(interfaceMaxPtIsMuF,
	"No",
	"",
	false);

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


	static Switch<ShowerHandler,bool> interfaceSplitHardProcess
	("SplitHardProcess",
	"Whether or not to try and split the hard process into production and decay processes",
	&ShowerHandler::splitHardProcess_, true, false, false);
	static SwitchOption interfaceSplitHardProcessYes
	(interfaceSplitHardProcess,
	"Yes",
	"Split the hard process",
	true);
	static SwitchOption interfaceSplitHardProcessNo
	(interfaceSplitHardProcess,
	"No",
	"Don't split the hard process",
	false);

	static Command<ShowerHandler> interfaceAddVariation
	("AddVariation",
	"Add a shower variation.",
	&ShowerHandler::doAddVariation, false);

	}

	Energy ShowerHandler::hardScale() const {
	return evolver_->hardScale();
	}

	void ShowerHandler::cascade() {

	tcPDFPtr first = firstPDF().pdf();
	tcPDFPtr second = secondPDF().pdf();

	if ( PDFA_ ) first = PDFA_;
	if ( PDFB_ ) second = PDFB_;

	resetPDFs(make_pair(first,second));

	if( ! rempdfs_.first)
	rempdfs_.first = PDFARemnant_ ? PDFPtr(PDFARemnant_) : const_ptr_cast<PDFPtr>(first);
	if( ! rempdfs_.second)
	rempdfs_.second = PDFBRemnant_ ? PDFPtr(PDFBRemnant_) : const_ptr_cast<PDFPtr>(second);

	// get the incoming partons
	tPPair incomingPartons =
	eventHandler()->currentCollision()->primarySubProcess()->incoming();
	// and the parton bins
	PBIPair incomingBins =
	make_pair(lastExtractor()->partonBinInstance(incomingPartons.first),
	lastExtractor()->partonBinInstance(incomingPartons.second));
	// and the incoming hadrons
	tPPair incomingHadrons =
	eventHandler()->currentCollision()->incoming();
	remDec_->setHadronContent(incomingHadrons);
	// check if incoming hadron == incoming parton
	// and get the incoming hadron if exists or parton otherwise
	incoming_ = make_pair(incomingBins.first ?
	incomingBins.first ->particle() : incomingPartons.first,
	incomingBins.second ?
	incomingBins.second->particle() : incomingPartons.second);
	// check the collision is of the beam particles
	// and if not boost collision to the right frame
	// i.e. the hadron-hadron CMF of the collision
	bool btotal(false);
	LorentzRotation rtotal;
	if(incoming_.first != incomingHadrons.first \|\|
	incoming_.second != incomingHadrons.second ) {
	btotal = true;
	boostCollision(false);
	}
	// set the current ShowerHandler
	currentHandler_ = this;
	// first shower the hard process
	useMe();
	try {
	SubProPtr sub = eventHandler()->currentCollision()->primarySubProcess();
	incomingPartons = cascade(sub,lastXCombPtr());
	}
	catch(ShowerTriesVeto &veto){
	throw Exception() << "Failed to generate the shower after "
	<< veto.tries
	<< " attempts in ShowerHandler::cascade()"
	<< Exception::eventerror;
	}
	if(showerHardProcessVeto()) throw Veto();
	// if a non-hadron collision return (both incoming non-hadronic)
	if( ( !incomingBins.first\|\|
	!isResolvedHadron(incomingBins.first ->particle()))&&
	( !incomingBins.second\|\|
	!isResolvedHadron(incomingBins.second->particle()))) {
	// boost back to lab if needed
	if(btotal) boostCollision(true);
	// unset the current ShowerHandler
	currentHandler_ = 0;
	return;
	}
	// get the remnants for hadronic collision
	pair<tRemPPtr,tRemPPtr> remnants(getRemnants(incomingBins));
	// set the starting scale of the forced splitting to the PDF freezing scale
	remDec_->initialize(remnants, incoming_, *currentStep(), pdfFreezingScale());
	// do the first forcedSplitting
	try {
	remDec_->doSplit(incomingPartons, make_pair(rempdfs_.first,rempdfs_.second), true);
	}
	catch (ExtraScatterVeto) {
	throw Exception() << "Remnant extraction failed in "
	<< "ShowerHandler::cascade() from primary interaction"
	<< Exception::eventerror;
	}
	// if no MPI return
	if( !isMPIOn() ) {
	remDec_->finalize();
	// boost back to lab if needed
	if(btotal) boostCollision(true);
	// unset the current ShowerHandler
	currentHandler_ = 0;
	return;
	}
	// generate the multiple scatters use modified pdf's now:
	setMPIPDFs();
	// additional "hard" processes
	unsigned int tries(0);
	// This is the loop over additional hard scatters (most of the time
	// only one, but who knows...)
	for(unsigned int i=1; i <= getMPIHandler()->additionalHardProcs(); i++){
	//counter for regeneration
	unsigned int multSecond = 0;
	// generate the additional scatters
	while( multSecond < getMPIHandler()->multiplicity(i) ) {
	// generate the hard scatter
	tStdXCombPtr lastXC = getMPIHandler()->generate(i);
	SubProPtr sub = lastXC->construct();
	// add to the Step
	newStep()->addSubProcess(sub);
	// increment the counters
	tries++;
	multSecond++;
	if(tries == maxtryDP_)
	throw Exception() << "Failed to establish the requested number "
	<< "of additional hard processes. If this error "
	<< "occurs often, your selection of additional "
	<< "scatter is probably unphysical"
	<< Exception::eventerror;
	// Generate the shower. If not possible veto the event
	try {
	incomingPartons = cascade(sub,lastXC);
	}
	catch(ShowerTriesVeto &veto){
	throw Exception() << "Failed to generate the shower of "
	<< "a secondary hard process after "
	<< veto.tries
	<< " attempts in Evolver::showerHardProcess()"
	<< Exception::eventerror;
	}
	try {
	// do the forcedSplitting
	remDec_->doSplit(incomingPartons, make_pair(remmpipdfs_.first,remmpipdfs_.second), false);
	}
	catch(ExtraScatterVeto){
	//remove all particles associated with the subprocess
	newStep()->removeParticle(incomingPartons.first);
	newStep()->removeParticle(incomingPartons.second);
	//remove the subprocess from the list
	newStep()->removeSubProcess(sub);
	//regenerate the scattering
	multSecond--;
	continue;
	}
	// connect with the remnants but don't set Remnant colour,
	// because that causes problems due to the multiple colour lines.
	if ( !remnants.first ->extract(incomingPartons.first , false) \|\|
	!remnants.second->extract(incomingPartons.second, false) )
	throw Exception() << "Remnant extraction failed in "
	<< "ShowerHandler::cascade() for additional scatter"
	<< Exception::runerror;
	}
	}
	// the underlying event processes
	unsigned int ptveto(1), veto(0);
	unsigned int max(getMPIHandler()->multiplicity());
	for(unsigned int i=0; i<max; i++) {
	// check how often this scattering has been regenerated
	if(veto > maxtryMPI_) break;
	//generate PSpoint
	tStdXCombPtr lastXC = getMPIHandler()->generate();
	SubProPtr sub = lastXC->construct();
	//If Algorithm=1 additional scatters of the signal type
	// with pt > ptmin have to be vetoed
	//with probability 1/(m+1), where m is the number of occurances in this event
	if( getMPIHandler()->Algorithm() == 1 ){
	//get the pT
	Energy pt = sub->outgoing().front()->momentum().perp();
	if(pt > getMPIHandler()->PtForVeto() && UseRandom::rnd() < 1./(ptveto+1) ){
	ptveto++;
	i--;
	continue;
	}
	}
	// add to the SubProcess to the step
	newStep()->addSubProcess(sub);
	// Run the Shower. If not possible veto the scattering
	try {
	incomingPartons = cascade(sub,lastXC);
	}
	// discard this extra scattering, but try the next one
	catch(ShowerTriesVeto) {
	newStep()->removeSubProcess(sub);
	//regenerate the scattering
	veto++;
	i--;
	continue;
	}
	try{
	//do the forcedSplitting
	remDec_->doSplit(incomingPartons, make_pair(remmpipdfs_.first,remmpipdfs_.second), false);
	}
	catch (ExtraScatterVeto) {
	//remove all particles associated with the subprocess
	newStep()->removeParticle(incomingPartons.first);
	newStep()->removeParticle(incomingPartons.second);
	//remove the subprocess from the list
	newStep()->removeSubProcess(sub);
	//regenerate the scattering
	veto++;
	i--;
	continue;
	}
	//connect with the remnants but don't set Remnant colour,
	//because that causes problems due to the multiple colour lines.
	if ( !remnants.first ->extract(incomingPartons.first , false) \|\|
	!remnants.second->extract(incomingPartons.second, false) )
	throw Exception() << "Remnant extraction failed in "
	<< "ShowerHandler::cascade() for MPI hard scattering"
	<< Exception::runerror;
	//reset veto counter
	veto = 0;
	}
	// finalize the remnants
	remDec_->finalize(getMPIHandler()->colourDisrupt(),
	getMPIHandler()->softMultiplicity());
	// boost back to lab if needed
	if(btotal) boostCollision(true);
	// unset the current ShowerHandler
	currentHandler_ = 0;
	getMPIHandler()->clean();
	}

	void ShowerHandler::fillEventRecord() {
	// create a new step
	StepPtr pstep = newStep();
	assert(!done_.empty());
	assert(done_[0]->isHard());
	// insert the steps
	for(unsigned int ix=0;ix<done_.size();++ix) {
	done_[ix]->fillEventRecord(pstep,
	evolver_->isISRadiationON(),
	evolver_->isFSRadiationON());
	}
	}

	void ShowerHandler::prepareCascade(tSubProPtr sub) {
	current_ = currentStep();
	subProcess_ = sub;
	+}
	+
	+void ShowerHandler::resetWeights() {
	for ( map<string,double>::iterator w = currentWeights_.begin();
	w != currentWeights_.end(); ++w ) {
	w->second = 1.0;
	}
	-}
	+}

	void ShowerHandler::combineWeights() {
	tEventPtr event = eventHandler()->currentEvent();
	for ( map<string,double>::const_iterator w =
	currentWeights_.begin(); w != currentWeights_.end(); ++w ) {
	map<string,double>::iterator ew = event->optionalWeights().find(w->first);
	if ( ew != event->optionalWeights().end() )
	ew->second *= w->second;
	else
	event->optionalWeights()[w->first] = w->second;
	}
	}

	string ShowerHandler::ShowerVariation::fromInFile(const string& in) {
	// pretty simple for the moment, just to try
	// TODO make this better
	istringstream read(in);
	string where;
	read >> renormalizationScaleFactor >> factorizationScaleFactor >> where;
	if ( !read )
	return "something went wrong with: " + in;
	if ( where == "FirstInteraction" \|\| where == "AllInteractions" )
	firstInteraction = true;
	else
	firstInteraction = false;
	if ( where == "SecondaryInteractions" \|\| where == "AllInteractions" )
	secondaryInteractions = true;
	else
	secondaryInteractions = false;
	return "";
	}

	void ShowerHandler::ShowerVariation::put(PersistentOStream& os) const {
	os << renormalizationScaleFactor << factorizationScaleFactor
	<< firstInteraction << secondaryInteractions;
	}

	void ShowerHandler::ShowerVariation::get(PersistentIStream& is) {
	is >> renormalizationScaleFactor >> factorizationScaleFactor
	>> firstInteraction >> secondaryInteractions;
	}

	string ShowerHandler::doAddVariation(string in) {
	if ( in.empty() )
	return "expecting a name and a variation specification";
	string name = StringUtils::car(in);
	ShowerVariation var;
	string res = var.fromInFile(StringUtils::cdr(in));
	if ( res.empty() ) {
	if ( !var.firstInteraction && !var.secondaryInteractions ) {
	// TODO what about decay showers?
	return "variation does not apply to any shower";
	}
	if ( var.renormalizationScaleFactor == 1.0 &&
	var.factorizationScaleFactor == 1.0 ) {
	return "variation does not vary anything";
	}
	/*
	Repository::clog() << "adding a variation with tag '" << name << "' using\nxir = "
	<< var.renormalizationScaleFactor
	<< " xif = "
	<< var.factorizationScaleFactor
	<< "\napplying to:\n"
	<< "first interaction = " << var.firstInteraction << " "
	<< "secondary interactions = " << var.secondaryInteractions << "\n"
	<< flush;
	*/
	showerVariations()[name] = var;
	}
	return res;
	}

	tPPair ShowerHandler::cascade(tSubProPtr sub,
	XCPtr xcomb) {
	prepareCascade(sub);
	// set the scale variation factors; needs to go after prepareCascade
	// to trigger possible different variations for hard and secondary
	// scatters
	evolver()->renormalizationScaleFactor(renormalizationScaleFactor());
	evolver()->factorizationScaleFactor(factorizationScaleFactor());
	evolver()->restrictPhasespace(restrictPhasespace());
	evolver()->hardScaleIsMuF(hardScaleIsMuF());
	// start of the try block for the whole showering process
	unsigned int countFailures=0;
	while (countFailures<maxtry_) {
	try {
	decay_.clear();
	done_.clear();
	ShowerTree::constructTrees(currentSubProcess(),hard_,decay_,
	firstInteraction() ? tagged() :
	tPVector(currentSubProcess()->outgoing().begin(),
	currentSubProcess()->outgoing().end()),
	splitHardProcess_);
	// if no hard process
	if(!hard_) throw Exception() << "Shower starting with a decay"
	<< "is not implemented"
	<< Exception::runerror;
	// perform the shower for the hard process
	evolver_->showerHardProcess(hard_,xcomb);
	done_.push_back(hard_);
	hard_->updateAfterShower(decay_);
	// if no decaying particles to shower break out of the loop
	if(decay_.empty()) break;
	// shower the decay products
	while(!decay_.empty()) {
	// find particle whose production process has been showered
	ShowerDecayMap::iterator dit = decay_.begin();
	while(!dit->second->parent()->hasShowered() && dit!=decay_.end()) ++dit;
	assert(dit!=decay_.end());
	// get the particle
	ShowerTreePtr decayingTree = dit->second;
	// remove it from the multimap
	decay_.erase(dit);
	// make sure the particle has been decayed
	decayingTree->decay(decay_);
	// now shower the decay
	evolver_->showerDecay(decayingTree);
	done_.push_back(decayingTree);
	decayingTree->updateAfterShower(decay_);
	}
	// suceeded break out of the loop
	break;
	}
	catch (KinematicsReconstructionVeto) {
	++countFailures;
	}
	}
	// if loop exited because of too many tries, throw event away
	if (countFailures >= maxtry_) {
	hard_=ShowerTreePtr();
	decay_.clear();
	done_.clear();
	throw Exception() << "Too many tries for main while loop "
	<< "in ShowerHandler::cascade()."
	<< Exception::eventerror;
	}
	//enter the particles in the event record
	fillEventRecord();
	// clear storage
	hard_=ShowerTreePtr();
	decay_.clear();
	done_.clear();
	// non hadronic case return
	if (!isResolvedHadron(incoming_.first ) &&
	!isResolvedHadron(incoming_.second) )
	return incoming_;
	// remake the remnants (needs to be after the colours are sorted
	// out in the insertion into the event record)
	if ( firstInteraction() ) return remakeRemnant(sub->incoming());
	//Return the new pair of incoming partons. remakeRemnant is not
	//necessary here, because the secondary interactions are not yet
	//connected to the remnants.
	return make_pair(findFirstParton(sub->incoming().first ),
	findFirstParton(sub->incoming().second));
	}

	ShowerHandler::RemPair
	ShowerHandler::getRemnants(PBIPair incomingBins) {
	RemPair remnants;
	// first beam particle
	if(incomingBins.first&&!incomingBins.first->remnants().empty()) {
	remnants.first =
	dynamic_ptr_cast<tRemPPtr>(incomingBins.first->remnants()[0] );
	if(remnants.first) {
	ParticleVector children=remnants.first->children();
	for(unsigned int ix=0;ix<children.size();++ix) {
	if(children[ix]->dataPtr()==remnants.first->dataPtr())
	remnants.first = dynamic_ptr_cast<RemPPtr>(children[ix]);
	}
	//remove existing colour lines from the remnants
	if(remnants.first->colourLine())
	remnants.first->colourLine()->removeColoured(remnants.first);
	if(remnants.first->antiColourLine())
	remnants.first->antiColourLine()->removeAntiColoured(remnants.first);
	}
	}
	// seconnd beam particle
	if(incomingBins.second&&!incomingBins. second->remnants().empty()) {
	remnants.second =
	dynamic_ptr_cast<tRemPPtr>(incomingBins.second->remnants()[0] );
	if(remnants.second) {
	ParticleVector children=remnants.second->children();
	for(unsigned int ix=0;ix<children.size();++ix) {
	if(children[ix]->dataPtr()==remnants.second->dataPtr())
	remnants.second = dynamic_ptr_cast<RemPPtr>(children[ix]);
	}
	//remove existing colour lines from the remnants
	if(remnants.second->colourLine())
	remnants.second->colourLine()->removeColoured(remnants.second);
	if(remnants.second->antiColourLine())
	remnants.second->antiColourLine()->removeAntiColoured(remnants.second);
	}
	}
	assert(remnants.first \|\| remnants.second);
	return remnants;
	}

	tPPair ShowerHandler::remakeRemnant(tPPair oldp){
	// get the parton extractor
	PartonExtractor & pex = *lastExtractor();
	// get the new partons
	tPPair newp = make_pair(findFirstParton(oldp.first ),
	findFirstParton(oldp.second));
	// if the same do nothing
	if(newp == oldp) return oldp;
	// Creates the new remnants and returns the new PartonBinInstances
	// ATTENTION Broken here for very strange configuration
	PBIPair newbins = pex.newRemnants(oldp, newp, newStep());
	newStep()->addIntermediate(newp.first);
	newStep()->addIntermediate(newp.second);
	// return the new partons
	return newp;
	}

	PPtr ShowerHandler::findFirstParton(tPPtr seed) const{
	if(seed->parents().empty()) return seed;
	tPPtr parent = seed->parents()[0];
	//if no parent there this is a loose end which will
	//be connected to the remnant soon.
	if(!parent \|\| parent == incoming_.first \|\|
	parent == incoming_.second ) return seed;
	else return findFirstParton(parent);
	}


	namespace {

	void addChildren(tPPtr in,set<tPPtr> & particles) {
	particles.insert(in);
	for(unsigned int ix=0;ix<in->children().size();++ix)
	addChildren(in->children()[ix],particles);
	}
	}

	void ShowerHandler::boostCollision(bool boost) {
	// calculate boost from lab to rest
	if(!boost) {
	Lorentz5Momentum ptotal=incoming_.first ->momentum()+incoming_.second->momentum();
	boost_ = LorentzRotation(-ptotal.boostVector());
	Axis axis((boost_*incoming_.first ->momentum()).vect().unit());
	if(axis.perp2()>0.) {
	double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
	boost_.rotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
	}
	}
	// first call performs the boost and second inverse
	// get the particles to be boosted
	set<tPPtr> particles;
	addChildren(incoming_.first,particles);
	addChildren(incoming_.second,particles);
	// apply the boost
	for(set<tPPtr>::const_iterator cit=particles.begin();
	cit!=particles.end();++cit) {
	(*cit)->transform(boost_);
	}
	if(!boost) boost_.invert();
	}

	void ShowerHandler::setMPIPDFs() {

	if ( !mpipdfs_.first ) {
	// first have to check for MinBiasPDF
	tcMinBiasPDFPtr first = dynamic_ptr_cast<tcMinBiasPDFPtr>(firstPDF().pdf());
	if(first)
	mpipdfs_.first = new_ptr(MPIPDF(first->originalPDF()));
	else
	mpipdfs_.first = new_ptr(MPIPDF(firstPDF().pdf()));
	}

	if ( !mpipdfs_.second ) {
	tcMinBiasPDFPtr second = dynamic_ptr_cast<tcMinBiasPDFPtr>(secondPDF().pdf());
	if(second)
	mpipdfs_.second = new_ptr(MPIPDF(second->originalPDF()));
	else
	mpipdfs_.second = new_ptr(MPIPDF(secondPDF().pdf()));
	}

	if( !remmpipdfs_.first ) {
	tcMinBiasPDFPtr first = dynamic_ptr_cast<tcMinBiasPDFPtr>(rempdfs_.first);
	if(first)
	remmpipdfs_.first = new_ptr(MPIPDF(first->originalPDF()));
	else
	remmpipdfs_.first = new_ptr(MPIPDF(rempdfs_.first));
	}

	if( !remmpipdfs_.second ) {
	tcMinBiasPDFPtr second = dynamic_ptr_cast<tcMinBiasPDFPtr>(rempdfs_.second);
	if(second)
	remmpipdfs_.second = new_ptr(MPIPDF(second->originalPDF()));
	else
	remmpipdfs_.second = new_ptr(MPIPDF(rempdfs_.second));
	}
	// reset the PDFs stored in the base class
	resetPDFs(mpipdfs_);

	}

	bool ShowerHandler::isResolvedHadron(tPPtr particle) {
	if(!HadronMatcher::Check(particle->data())) return false;
	for(unsigned int ix=0;ix<particle->children().size();++ix) {
	if(particle->children()[ix]->id()==ParticleID::Remnant) return true;
	}
	return false;
	}

	HardTreePtr ShowerHandler::generateCKKW(ShowerTreePtr ) const {
	return HardTreePtr();
	}

	diff --git a/Shower/ShowerHandler.h b/Shower/ShowerHandler.h
	--- a/Shower/ShowerHandler.h
	+++ b/Shower/ShowerHandler.h
	@@ -1,676 +1,681 @@
	// -- C++ --
	//
	// ShowerHandler.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 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_ShowerHandler_H
	#define HERWIG_ShowerHandler_H
	//
	// This is the declaration of the ShowerHandler class.
	//

	#include "ThePEG/Handlers/EventHandler.h"
	#include "ThePEG/Handlers/CascadeHandler.h"
	#include "Herwig/Shower/UEBase.h"
	#include "Herwig/Shower/Base/Evolver.fh"
	#include "Herwig/Shower/Base/ShowerParticle.fh"
	#include "Herwig/Shower/Base/ShowerTree.fh"
	#include "Herwig/Shower/Base/HardTree.fh"
	#include "Herwig/PDF/HwRemDecayer.fh"
	#include "ThePEG/EventRecord/RemnantParticle.fh"
	#include "ShowerHandler.fh"
	#include "Herwig/MatrixElement/Matchbox/Matching/HardScaleProfile.h"

	namespace Herwig {

	/**
	* Typedef for the ShowerTree for the decays
	*/
	typedef multimap<Energy,ShowerTreePtr,std::greater<Energy> > ShowerDecayMap;

	using namespace ThePEG;

	/** \ingroup Shower
	*
	* This class is the main driver of the shower: it is responsible for
	* the proper handling of all other specific collaborating classes
	* and for the storing of the produced particles in the event record.
	*
	* @see \ref ShowerHandlerInterfaces "The interfaces"
	*
	* @see ThePEG::CascadeHandler
	* @see MPIHandler
	* @see HwRemDecayer
	*/
	class ShowerHandler: public CascadeHandler {

	public:

	/** Typedef for a pair of ThePEG::RemnantParticle pointers. */
	typedef pair<tRemPPtr, tRemPPtr> RemPair;

	/**
	* The default constructor.
	*/
	ShowerHandler();

	/**
	* Destructor
	*/
	virtual ~ShowerHandler();

	public:

	/**
	* The main method which manages the multiple interactions and starts
	* the shower by calling cascade(sub, lastXC).
	*/
	virtual void cascade();

	/**
	* Hook to allow vetoing of event after showering hard sub-process
	* as in e.g. MLM merging.
	*/
	virtual bool showerHardProcessVeto() { return false; };

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

	public:

	/*@name Methods related to PDF freezing /
	//@{
	/**
	* Get the PDF freezing scale
	*/
	Energy pdfFreezingScale() const {return pdfFreezingScale_;}
	//@}

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

	public:

	/** @name Functions to access information. */
	//@{

	/**
	* Return true if currently the primary subprocess is showered.
	*/
	bool firstInteraction() const {
	return ( subProcess_ ==
	eventHandler()->currentCollision()->primarySubProcess() );
	}

	/**
	* Return the currently used SubProcess.
	*/
	tSubProPtr currentSubProcess() const {
	assert(subProcess_);
	return subProcess_;
	}

	/**
	* Return true if multiple parton interactions are switched on
	* and can be used for this beam setup.
	*/
	bool isMPIOn() const {
	return MPIHandler_ && MPIHandler_->beamOK();
	}

	/**
	* Return the remnant decayer.
	*/
	tHwRemDecPtr remnantDecayer() const { return remDec_; }
	//@}

	/**
	* Access to the Evolver
	*/
	tEvolverPtr evolver() const {return evolver_;}

	/**
	* Generate hard emissions for CKKW etc
	*/
	virtual HardTreePtr generateCKKW(ShowerTreePtr tree) const;

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

	/**
	* The factorization scale factor.
	*/
	double factorizationScaleFactor() const {
	if ( scaleFactorOption_ == 0 \|\| !subProcess_ )
	return factorizationScaleFactor_;
	if ( scaleFactorOption_ == 1 )
	return firstInteraction() ? factorizationScaleFactor_ : 1.0;
	if ( scaleFactorOption_ == 2 )
	return !firstInteraction() ? factorizationScaleFactor_ : 1.0;
	return 1.0;
	}

	/**
	* The renormalization scale factor.
	*/
	double renormalizationScaleFactor() const {
	if ( scaleFactorOption_ == 0 \|\| !subProcess_ )
	return renormalizationScaleFactor_;
	if ( scaleFactorOption_ == 1 )
	return firstInteraction() ? renormalizationScaleFactor_ : 1.0;
	if ( scaleFactorOption_ == 2 )
	return !firstInteraction() ? renormalizationScaleFactor_ : 1.0;
	return 1.0;
	}

	/**
	* The scale factor for the hard scale
	*/
	double hardScaleFactor() const {
	if ( scaleFactorOption_ == 0 \|\| !subProcess_ )
	return hardScaleFactor_;
	if ( scaleFactorOption_ == 1 )
	return firstInteraction() ? hardScaleFactor_ : 1.0;
	if ( scaleFactorOption_ == 2 )
	return !firstInteraction() ? hardScaleFactor_ : 1.0;
	return 1.0;
	}

	/**
	* The option on when to apply the scale factors
	*/
	int scaleFactorOption() const { return scaleFactorOption_; }

	/**
	* Return true, if the phase space restrictions of the dipole shower should
	* be applied.
	*/
	bool restrictPhasespace() const { return restrictPhasespace_; }

	/**
	* Return profile scales
	*/
	Ptr<HardScaleProfile>::tptr profileScales() const { return hardScaleProfile_; }

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

	/**
	* Return true if maximum pt should be deduced from the factorization scale
	*/
	bool hardScaleIsMuF() const { return maxPtIsMuF_; }

	/**
	* A struct identifying a shower variation
	*/
	struct ShowerVariation {

	/**
	* Vary the renormalization scale by the given factor.
	*/
	double renormalizationScaleFactor;

	/**
	* Vary the factorization scale by the given factor.
	*/
	double factorizationScaleFactor;

	/**
	* Apply the variation to the first interaction
	*/
	bool firstInteraction;

	/**
	* Apply the variation to the secondary interactions
	*/
	bool secondaryInteractions;

	/**
	* Default constructor
	*/
	ShowerVariation()
	: renormalizationScaleFactor(1.0),
	factorizationScaleFactor(1.0),
	firstInteraction(true),
	secondaryInteractions(false) {}

	/**
	* Parse from in file command
	*/
	string fromInFile(const string&);

	/**
	* Put to persistent stream
	*/
	void put(PersistentOStream& os) const;

	/**
	* Get from persistent stream
	*/
	void get(PersistentIStream& is);

	};

	/**
	* Access the shower variations
	*/
	map<string,ShowerVariation>& showerVariations() {
	return showerVariations_;
	}

	/**
	* Return the shower variations
	*/
	const map<string,ShowerVariation>& showerVariations() const {
	return showerVariations_;
	}

	protected:

	/**
	* The shower variation weights
	*/
	map<string,double> currentWeights_;

	/**
	* Combine the variation weights which have been encountered
	*/
	void combineWeights();

	+ /**
	+ * Reset the current weights
	+ */
	+ void resetWeights();
	+
	protected:

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

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

	protected:

	/**
	* Prepare to shower the given subprocess
	*/
	void prepareCascade(tSubProPtr sub);

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

	/**
	* Return the maximum number of attempts for showering
	* a given subprocess.
	*/
	unsigned int maxtry() const { return maxtry_; }

	/**
	* At the end of the Showering, transform ShowerParticle objects
	* into ThePEG particles and fill the event record with them.
	* Notice that the parent/child relationships and the
	* transformation from ShowerColourLine objects into ThePEG
	* ColourLine ones must be properly handled.
	*/
	void fillEventRecord();

	/**
	* Find the parton extracted from the incoming particle after ISR
	*/
	PPtr findFirstParton(tPPtr seed) const;

	/**
	* Fix Remnant connections after ISR
	*/
	tPPair remakeRemnant(tPPair oldp);

	/**
	* Get the remnants from the ThePEG::PartonBinInstance es and
	* do some checks.
	*/
	RemPair getRemnants(PBIPair incbins);

	/**
	* Make the remnant after the shower
	*/
	void makeRemnants();

	/**
	* Reset the PDF's after the hard collision has been showered
	*/
	void setMPIPDFs();

	/**
	* Boost all the particles in the collision so that the collision always occurs
	* in the rest frame with the incoming particles along the z axis
	*/
	void boostCollision(bool boost);

	/**
	* Is a beam particle where hadronic structure is resolved
	*/
	bool isResolvedHadron(tPPtr);

	protected:

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

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

	/**
	* Called at the end of the run phase.
	*/
	virtual void dofinish();
	//@}

	private:

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

	private:

	/**
	* Access function for the MPIHandler, it should only be called after
	* checking with isMPIOn.
	*/
	tUEBasePtr getMPIHandler() const {
	assert(MPIHandler_);
	return MPIHandler_;
	}

	private:

	/**
	* a MPIHandler to administer the creation of several (semihard)
	* partonic interactions.
	*/
	UEBasePtr MPIHandler_;

	/**
	* Pointer to the evolver
	*/
	EvolverPtr evolver_;

	/**
	* Pointer to the HwRemDecayer
	*/
	HwRemDecPtr remDec_;

	/**
	* The PDF for beam particle A. Overrides the particle's own PDF setting.
	*/
	PDFPtr PDFA_;

	/**
	* The PDF for beam particle B. Overrides the particle's own PDF setting.
	*/
	PDFPtr PDFB_;

	/**
	* The PDF for beam particle A for remnant splitting. Overrides the particle's own PDF setting.
	*/
	PDFPtr PDFARemnant_;

	/**
	* The PDF for beam particle B for remnant splitting. Overrides the particle's own PDF setting.
	*/
	PDFPtr PDFBRemnant_;

	/**
	* The PDF freezing scale
	*/
	Energy pdfFreezingScale_;

	/**
	* Maximum number of attempts for the
	* main showering loop
	*/
	unsigned int maxtry_;

	/**
	* Maximum number of attempts for the regeneration of an additional
	* scattering, before the number of scatters is reduced.
	*/
	unsigned int maxtryMPI_;

	/**
	* Maximum number of attempts for the regeneration of an additional
	* hard scattering, before this event is vetoed.
	*/
	unsigned int maxtryDP_;

	/**
	* PDG codes of the particles which decay during showering
	* this is fast storage for use during running
	*/
	set<long> particlesDecayInShower_;

	/**
	* PDG codes of the particles which decay during showering
	* this is a vector that is interfaced so they can be changed
	*/
	vector<long> inputparticlesDecayInShower_;

	/**
	* Whether or not to include spa-cetime distances in the shower
	*/
	bool includeSpaceTime_;

	/**
	* The minimum virtuality for the space-time model
	*/
	Energy2 vMin_;

	/**
	* The ShowerTree for the hard process
	*/
	ShowerTreePtr hard_;

	/**
	* The incoming beam particles for the current collision
	*/
	tPPair incoming_;

	/**
	* The ShowerTree for the decays
	*/
	ShowerDecayMap decay_;

	/**
	* The ShowerTrees for which the initial shower
	*/
	vector<ShowerTreePtr> done_;

	/**
	* Const pointer to the current step
	*/
	tcStepPtr current_;

	/**
	* Const pointer to the currently handeled ThePEG::SubProcess
	*/
	tSubProPtr subProcess_;

	/**
	* pointer to "this", the current ShowerHandler.
	*/
	static ShowerHandler * currentHandler_;

	/**
	* Boost to get back to the lab
	*/
	LorentzRotation boost_;

	/**
	* The MPI PDF's to be used for secondary scatters.
	*/
	pair <PDFPtr, PDFPtr> mpipdfs_;

	/**
	* The MPI PDF's to be used for secondary scatters.
	*/
	pair <PDFPtr, PDFPtr> rempdfs_;

	/**
	* The MPI PDF's to be used for secondary scatters.
	*/
	pair <PDFPtr, PDFPtr> remmpipdfs_;

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

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

	/**
	* The scale factor for the hard scale
	*/
	double hardScaleFactor_;

	/**
	* The option on when to apply the scale factors
	*/
	int scaleFactorOption_;

	/**
	* True, if the phase space restrictions of the dipole shower should
	* be applied.
	*/
	bool restrictPhasespace_;

	/**
	* True if maximum pt should be deduced from the factorization scale
	*/
	bool maxPtIsMuF_;

	/**
	* The profile scales
	*/
	Ptr<HardScaleProfile>::ptr hardScaleProfile_;

	/**
	* Whether or not to split into hard and decay trees
	*/
	bool splitHardProcess_;

	/**
	* The shower variations
	*/
	map<string,ShowerVariation> showerVariations_;

	/**
	* Command to add a shower variation
	*/
	string doAddVariation(string);

	public:

	/**
	* struct that is used to catch exceptions which are thrown
	* due to energy conservation issues of additional scatters
	*/
	struct ExtraScatterVeto {};

	/**
	* struct that is used to catch exceptions which are thrown
	* due to fact that the Shower has been invoked more than
	* a defined threshold on a certain configuration
	*/
	struct ShowerTriesVeto {
	/** variable to store the number of attempts */
	const int tries;

	/** constructor */
	ShowerTriesVeto(int t) : tries(t) {}
	};

	/**
	* pointer to "this", the current ShowerHandler.
	*/
	static const ShowerHandler * currentHandler() {
	assert(currentHandler_);
	return currentHandler_;
	}

	protected:

	/**
	* Set the current handler
	*/
	void setCurrentHandler() {
	currentHandler_ = this;
	}

	};

	inline PersistentOStream& operator<<(PersistentOStream& os, const ShowerHandler::ShowerVariation& var) {
	var.put(os); return os;
	}

	inline PersistentIStream& operator>>(PersistentIStream& is, ShowerHandler::ShowerVariation& var) {
	var.get(is); return is;
	}

	}

	#endif /* HERWIG_ShowerHandler_H */

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 4:51 PM (1 d, 10 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3799224
Default Alt Text: (84 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions