No OneTemporary
Actions

Size

220 KB

Subscribers

None

View Options

	diff --git a/Shower/QTilde/QTildeShowerHandler.cc b/Shower/QTilde/QTildeShowerHandler.cc
	--- a/Shower/QTilde/QTildeShowerHandler.cc
	+++ b/Shower/QTilde/QTildeShowerHandler.cc
	@@ -1,3788 +1,3700 @@
	// -- C++ --
	//
	// QTildeShowerHandler.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 QTildeShowerHandler class.
	//

	#include "QTildeShowerHandler.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Interface/RefVector.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 "ThePEG/Utilities/EnumIO.h"
	#include "Herwig/Shower/QTilde/Base/ShowerParticle.h"
	#include "Herwig/PDF/MPIPDF.h"
	#include "Herwig/PDF/MinBiasPDF.h"
	#include "Herwig/Shower/QTilde/Base/ShowerTree.h"
	#include "Herwig/Shower/QTilde/Base/HardTree.h"
	#include "Herwig/Shower/QTilde/Base/KinematicsReconstructor.h"
	#include "Herwig/Shower/QTilde/Base/PartnerFinder.h"
	#include "Herwig/PDF/HwRemDecayer.h"
	#include "Herwig/Shower/QTilde/Base/ShowerVertex.h"
	#include "ThePEG/Repository/CurrentGenerator.h"
	#include "Herwig/MatrixElement/Matchbox/Base/SubtractedME.h"
	#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h"
	#include "ThePEG/PDF/PartonExtractor.h"
	#include "Herwig/Shower/RealEmissionProcess.h"

	using namespace Herwig;

	bool QTildeShowerHandler::_hardEmissionWarn = true;
	bool QTildeShowerHandler::_missingTruncWarn = true;

	QTildeShowerHandler::QTildeShowerHandler() :
	_maxtry(100), _meCorrMode(1), _reconOpt(0),
	_hardVetoReadOption(false),
	_iptrms(ZERO), _beta(0.), _gamma(ZERO), _iptmax(),
	_limitEmissions(0), _initialenhance(1.), _finalenhance(1.),
	_nReWeight(100), _reWeight(false),
	interaction_(ShowerInteraction::Both),
	_trunc_Mode(true), _hardEmission(1),
	_spinOpt(1), _softOpt(2), _hardPOWHEG(false), muPt(ZERO),
	_maxTryFSR(100000), _maxFailFSR(100), _fracFSR(0.001),
	_nFSR(0), _nFailedFSR(0)
	{}

	QTildeShowerHandler::~QTildeShowerHandler() {}

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

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

	void QTildeShowerHandler::persistentOutput(PersistentOStream & os) const {
	os << _model << _splittingGenerator << _maxtry
	<< _meCorrMode << _hardVetoReadOption
	<< _limitEmissions << _spinOpt << _softOpt << _hardPOWHEG
	<< ounit(_iptrms,GeV) << _beta << ounit(_gamma,GeV) << ounit(_iptmax,GeV)
	<< _vetoes << _fullShowerVetoes << _nReWeight << _reWeight
	<< _trunc_Mode << _hardEmission << _reconOpt
	<< ounit(muPt,GeV)
	<< oenum(interaction_) << _maxTryFSR << _maxFailFSR << _fracFSR;
	}

	void QTildeShowerHandler::persistentInput(PersistentIStream & is, int) {
	is >> _model >> _splittingGenerator >> _maxtry
	>> _meCorrMode >> _hardVetoReadOption
	>> _limitEmissions >> _spinOpt >> _softOpt >> _hardPOWHEG
	>> iunit(_iptrms,GeV) >> _beta >> iunit(_gamma,GeV) >> iunit(_iptmax,GeV)
	>> _vetoes >> _fullShowerVetoes >> _nReWeight >> _reWeight
	>> _trunc_Mode >> _hardEmission >> _reconOpt
	>> iunit(muPt,GeV)
	>> ienum(interaction_) >> _maxTryFSR >> _maxFailFSR >> _fracFSR;
	}


	// The following static variable is needed for the type
	// description system in ThePEG.
	DescribeClass<QTildeShowerHandler,ShowerHandler>
	describeHerwigQTildeShowerHandler("Herwig::QTildeShowerHandler", "HwShower.so");

	void QTildeShowerHandler::Init() {

	static ClassDocumentation<QTildeShowerHandler> documentation
	("TheQTildeShowerHandler class is the main class"
	" for the angular-ordered parton shower",
	"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<QTildeShowerHandler,SplittingGenerator>
	interfaceSplitGen("SplittingGenerator",
	"A reference to the SplittingGenerator object",
	&Herwig::QTildeShowerHandler::_splittingGenerator,
	false, false, true, false);

	static Reference<QTildeShowerHandler,ShowerModel> interfaceShowerModel
	("ShowerModel",
	"The pointer to the object which defines the shower evolution model.",
	&QTildeShowerHandler::_model, false, false, true, false, false);

	static Parameter<QTildeShowerHandler,unsigned int> interfaceMaxTry
	("MaxTry",
	"The maximum number of attempts to generate the shower from a"
	" particular ShowerTree",
	&QTildeShowerHandler::_maxtry, 100, 1, 100000,
	false, false, Interface::limited);

	static Parameter<QTildeShowerHandler,unsigned int> interfaceNReWeight
	("NReWeight",
	"The number of attempts for the shower when reweighting",
	&QTildeShowerHandler::_nReWeight, 100, 10, 10000,
	false, false, Interface::limited);

	static Switch<QTildeShowerHandler, unsigned int> ifaceMECorrMode
	("MECorrMode",
	"Choice of the ME Correction Mode",
	&QTildeShowerHandler::_meCorrMode, 1, false, false);
	static SwitchOption on
	(ifaceMECorrMode,"HardPlusSoft","hard+soft on", 1);
	static SwitchOption hard
	(ifaceMECorrMode,"Hard","only hard on", 2);
	static SwitchOption soft
	(ifaceMECorrMode,"Soft","only soft on", 3);

	static Switch<QTildeShowerHandler, bool> ifaceHardVetoReadOption
	("HardVetoReadOption",
	"Apply read-in scale veto to all collisions or just the primary one?",
	&QTildeShowerHandler::_hardVetoReadOption, false, false, false);
	static SwitchOption AllCollisions
	(ifaceHardVetoReadOption,
	"AllCollisions",
	"Read-in pT veto applied to primary and secondary collisions.",
	false);
	static SwitchOption PrimaryCollision
	(ifaceHardVetoReadOption,
	"PrimaryCollision",
	"Read-in pT veto applied to primary but not secondary collisions.",
	true);

	static Parameter<QTildeShowerHandler, Energy> ifaceiptrms
	("IntrinsicPtGaussian",
	"RMS of intrinsic pT of Gaussian distribution:\n"
	"2(1-Beta)exp(-sqr(intrinsicpT/RMS))/sqr(RMS)",
	&QTildeShowerHandler::_iptrms, GeV, ZERO, ZERO, 1000000.0*GeV,
	false, false, Interface::limited);

	static Parameter<QTildeShowerHandler, double> ifacebeta
	("IntrinsicPtBeta",
	"Proportion of inverse quadratic distribution in generating intrinsic pT.\n"
	"(1-Beta) is the proportion of Gaussian distribution",
	&QTildeShowerHandler::_beta, 0, 0, 1,
	false, false, Interface::limited);

	static Parameter<QTildeShowerHandler, Energy> ifacegamma
	("IntrinsicPtGamma",
	"Parameter for inverse quadratic:\n"
	"2BetaGamma/(sqr(Gamma)+sqr(intrinsicpT))",
	&QTildeShowerHandler::_gamma,GeV, ZERO, ZERO, 100000.0*GeV,
	false, false, Interface::limited);

	static Parameter<QTildeShowerHandler, Energy> ifaceiptmax
	("IntrinsicPtIptmax",
	"Upper bound on intrinsic pT for inverse quadratic",
	&QTildeShowerHandler::_iptmax,GeV, ZERO, ZERO, 100000.0*GeV,
	false, false, Interface::limited);

	static RefVector<QTildeShowerHandler,ShowerVeto> ifaceVetoes
	("Vetoes",
	"The vetoes to be checked during showering",
	&QTildeShowerHandler::_vetoes, -1,
	false,false,true,true,false);

	static RefVector<QTildeShowerHandler,FullShowerVeto> interfaceFullShowerVetoes
	("FullShowerVetoes",
	"The vetos to be appliede on the full final state of the shower",
	&QTildeShowerHandler::_fullShowerVetoes, -1, false, false, true, false, false);

	static Switch<QTildeShowerHandler,unsigned int> interfaceLimitEmissions
	("LimitEmissions",
	"Limit the number and type of emissions for testing",
	&QTildeShowerHandler::_limitEmissions, 0, false, false);
	static SwitchOption interfaceLimitEmissionsNoLimit
	(interfaceLimitEmissions,
	"NoLimit",
	"Allow an arbitrary number of emissions",
	0);
	static SwitchOption interfaceLimitEmissionsOneInitialStateEmission
	(interfaceLimitEmissions,
	"OneInitialStateEmission",
	"Allow one emission in the initial state and none in the final state",
	1);
	static SwitchOption interfaceLimitEmissionsOneFinalStateEmission
	(interfaceLimitEmissions,
	"OneFinalStateEmission",
	"Allow one emission in the final state and none in the initial state",
	2);
	static SwitchOption interfaceLimitEmissionsHardOnly
	(interfaceLimitEmissions,
	"HardOnly",
	"Only allow radiation from the hard ME correction",
	3);
	static SwitchOption interfaceLimitEmissionsOneEmission
	(interfaceLimitEmissions,
	"OneEmission",
	"Allow one emission in either the final state or initial state, but not both",
	4);

	static Switch<QTildeShowerHandler,bool> interfaceTruncMode
	("TruncatedShower", "Include the truncated shower?",
	&QTildeShowerHandler::_trunc_Mode, 1, false, false);
	static SwitchOption interfaceTruncMode0
	(interfaceTruncMode,"No","Truncated Shower is OFF", 0);
	static SwitchOption interfaceTruncMode1
	(interfaceTruncMode,"Yes","Truncated Shower is ON", 1);

	static Switch<QTildeShowerHandler,int> interfaceHardEmission
	("HardEmission",
	"Whether to use ME corrections or POWHEG for the hardest emission",
	&QTildeShowerHandler::_hardEmission, 0, false, false);
	static SwitchOption interfaceHardEmissionNone
	(interfaceHardEmission,
	"None",
	"No Corrections",
	0);
	static SwitchOption interfaceHardEmissionMECorrection
	(interfaceHardEmission,
	"MECorrection",
	"Old fashioned ME correction",
	1);
	static SwitchOption interfaceHardEmissionPOWHEG
	(interfaceHardEmission,
	"POWHEG",
	"Powheg style hard emission",
	2);

	static Switch<QTildeShowerHandler,ShowerInteraction> interfaceInteractions
	("Interactions",
	"The interactions to be used in the shower",
	&QTildeShowerHandler::interaction_, ShowerInteraction::Both, false, false);
	static SwitchOption interfaceInteractionsQCD
	(interfaceInteractions,
	"QCD",
	"Only QCD radiation",
	ShowerInteraction::QCD);
	static SwitchOption interfaceInteractionsQED
	(interfaceInteractions,
	"QED",
	"Only QEd radiation",
	ShowerInteraction::QED);
	static SwitchOption interfaceInteractionsQCDandQED
	(interfaceInteractions,
	"QCDandQED",
	"Both QED and QCD radiation",
	ShowerInteraction::Both);

	static Switch<QTildeShowerHandler,unsigned int> interfaceReconstructionOption
	("ReconstructionOption",
	"Treatment of the reconstruction of the transverse momentum of "
	"a branching from the evolution scale.",
	&QTildeShowerHandler::_reconOpt, 0, false, false);
	static SwitchOption interfaceReconstructionOptionCutOff
	(interfaceReconstructionOption,
	"CutOff",
	"Use the cut-off masses in the calculation",
	0);
	static SwitchOption interfaceReconstructionOptionOffShell
	(interfaceReconstructionOption,
	"OffShell",
	"Use the off-shell masses in the calculation veto the emission of the parent,"
	" no veto in generation of emissions from children",
	1);
	static SwitchOption interfaceReconstructionOptionOffShell2
	(interfaceReconstructionOption,
	"OffShell2",
	"Use the off-shell masses in the calculation veto the emissions from the children."
	" no veto in generation of emissions from children",
	2);
	static SwitchOption interfaceReconstructionOptionOffShell3
	(interfaceReconstructionOption,
	"OffShell3",
	"Use the off-shell masses in the calculation veto the emissions from the children."
	" veto in generation of emissions from children using cut-off for second parton",
	3);
	static SwitchOption interfaceReconstructionOptionOffShell4
	(interfaceReconstructionOption,
	"OffShell4",
	"Ass OffShell3 but with a restriction on the mass of final-state"
	" jets produced via backward evolution.",
	4);

	static Switch<QTildeShowerHandler,unsigned int> interfaceSpinCorrelations
	("SpinCorrelations",
	"Treatment of spin correlations in the parton shower",
	&QTildeShowerHandler::_spinOpt, 1, false, false);
	static SwitchOption interfaceSpinCorrelationsOff
	(interfaceSpinCorrelations,
	"No",
	"No spin correlations",
	0);
	static SwitchOption interfaceSpinCorrelationsSpin
	(interfaceSpinCorrelations,
	"Yes",
	"Include the azimuthal spin correlations only",
	1);

	static Switch<QTildeShowerHandler,unsigned int> interfaceSoftCorrelations
	("SoftCorrelations",
	"Option for the treatment of soft correlations in the parton shower",
	&QTildeShowerHandler::_softOpt, 2, false, false);
	static SwitchOption interfaceSoftCorrelationsNone
	(interfaceSoftCorrelations,
	"No",
	"No soft correlations",
	0);
	static SwitchOption interfaceSoftCorrelationsFull
	(interfaceSoftCorrelations,
	"Full",
	"Use the full eikonal",
	1);
	static SwitchOption interfaceSoftCorrelationsSingular
	(interfaceSoftCorrelations,
	"Singular",
	"Use original Webber-Marchisini form",
	2);

	static Switch<QTildeShowerHandler,bool> interfaceHardPOWHEG
	("HardPOWHEG",
	"Treatment of powheg emissions which are too hard to have a shower interpretation",
	&QTildeShowerHandler::_hardPOWHEG, false, false, false);
	static SwitchOption interfaceHardPOWHEGAsShower
	(interfaceHardPOWHEG,
	"AsShower",
	"Still interpret as shower emissions",
	false);
	static SwitchOption interfaceHardPOWHEGRealEmission
	(interfaceHardPOWHEG,
	"RealEmission",
	"Generate shower from the real emmission configuration",
	true);

	static Parameter<QTildeShowerHandler,unsigned int> interfaceMaxTryFSR
	("MaxTryFSR",
	"The maximum number of attempted FSR emissions in"
	" the generation of the FSR",
	&QTildeShowerHandler::_maxTryFSR, 100000, 10, 100000000,
	false, false, Interface::limited);

	static Parameter<QTildeShowerHandler,unsigned int> interfaceMaxFailFSR
	("MaxFailFSR",
	"Maximum number of failures generating the FSR",
	&QTildeShowerHandler::_maxFailFSR, 100, 1, 100000000,
	false, false, Interface::limited);

	static Parameter<QTildeShowerHandler,double> interfaceFSRFailureFraction
	("FSRFailureFraction",
	"Maximum fraction of events allowed to fail due to too many FSR emissions",
	&QTildeShowerHandler::_fracFSR, 0.001, 1e-10, 1,
	false, false, Interface::limited);

	}

	-void QTildeShowerHandler::cascade(tPVector finalstate) {
	-
	- hard_=ShowerTreePtr();
	- decay_.clear();
	- done_.clear();
	- tStdXCombPtr xcomb;
	- StepPtr pstep = newStep();
	-
	- unsigned int countFailures=0;
	- while (countFailures<maxtry()) {
	- try {
	- decay_.clear();
	- done_.clear();
	- DecayProcessMap decay;
	-
	-
	- PerturbativeProcessPtr hard=new_ptr(PerturbativeProcess());
	- PPtr p1,p2,p1tmp,p2tmp;
	- for(tPVector::iterator it=finalstate.begin();it!=finalstate.end();++it){
	- PPtr p=(*it);
	- while (!p->parents().empty()) {
	- if (p->parents().size()==2&&((!p1&&!p2)\|\|p1==p2)) {
	- p1=p->parents()[1];
	- p2=p->parents()[0];
	- while (!p2->parents().empty()) {
	- p2=p2->parents()[0];
	- }
	- while (!p1->parents().empty()) {
	- p1=p1->parents()[0];
	- }
	- }
	- p=p->parents()[0];
	- }
	- }
	-
	- assert(p1!=p2);
	- hard->incoming().push_back(make_pair(p1,PerturbativeProcessPtr()));
	- hard->incoming().push_back(make_pair(p2,PerturbativeProcessPtr()));
	- for(tPVector::iterator it=finalstate.begin();it!=finalstate.end();++it){
	- PPtr p=(*it);
	- PerturbativeProcessPtr newDecay=new_ptr(PerturbativeProcess());
	- if((decaysInShower((it).id())&&!(it).dataPtr()->stable())){
	- createDecayProcess(*it,hard,decay);
	- }else{
	- hard->outgoing().push_back(make_pair(*it,PerturbativeProcessPtr()));
	- }
	- }
	-
	- setCurrentHandler();
	- ShowerTree::constructTrees(hard_,decay_,hard,decay);
	-
	- done_.push_back(hard_);
	-
	-
	- if(decay_.empty()) return;
	-
	- while(!decay_.empty()) {
	- ShowerDecayMap::iterator dit = decay_.begin();
	- if (!dit->second->parent()) {
	- decay_.erase(dit);
	- continue;
	- }
	- ShowerTreePtr decayingTree = dit->second;
	- decay_.erase(dit);
	- QTildeShowerHandler::decay(decayingTree,decay_);
	- currentTree(decayingTree);
	- showerDecay(decayingTree);
	- done_.push_back(decayingTree);
	- decayingTree->updateAfterShower(decay_);
	-
	- }
	- break;
	- }
	- catch (KinematicsReconstructionVeto) {
	- resetWeights();
	- ++countFailures;
	- }
	- }
	-
	-
	- fillEventRecord();
	- setDidRunCascade(true);
	- hard_=ShowerTreePtr();
	- decay_.clear();
	- done_.clear();
	- return;
	-}
	-
	tPPair QTildeShowerHandler::cascade(tSubProPtr sub,
	XCPtr xcomb) {
	// use me for reference in tex file etc
	useMe();
	prepareCascade(sub);
	// set things up in the base class
	resetWeights();
	hard_=ShowerTreePtr();
	decay_.clear();
	done_.clear();
	// check if anything needs doing
	if ( !doFSR() && ! doISR() )
	return sub->incoming();
	// start of the try block for the whole showering process
	unsigned int countFailures=0;
	while (countFailures<maxtry()) {
	try {
	decay_.clear();
	done_.clear();
	PerturbativeProcessPtr hard;
	DecayProcessMap decay;
	splitHardProcess(firstInteraction() ? tagged() :
	tPVector(currentSubProcess()->outgoing().begin(),
	currentSubProcess()->outgoing().end()),
	hard,decay);
	ShowerTree::constructTrees(hard_,decay_,hard,decay);
	// 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
	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
	QTildeShowerHandler::decay(decayingTree,decay_);
	// now shower the decay
	showerDecay(decayingTree);
	done_.push_back(decayingTree);
	decayingTree->updateAfterShower(decay_);
	}
	// suceeded break out of the loop
	break;
	}
	catch (KinematicsReconstructionVeto) {
	resetWeights();
	++countFailures;
	}
	catch ( ... ) {
	hard_=ShowerTreePtr();
	decay_.clear();
	done_.clear();
	throw;
	}
	}
	// if loop exited because of too many tries, throw event away
	if (countFailures >= maxtry()) {
	resetWeights();
	hard_=ShowerTreePtr();
	decay_.clear();
	done_.clear();
	throw Exception() << "Too many tries for main while loop "
	<< "in QTildeShowerHandler::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(incomingBeams().first ) &&
	!isResolvedHadron(incomingBeams().second) )
	return incomingBeams();
	// 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));
	}

	void QTildeShowerHandler::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,doISR(),doFSR());
	}
	}

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

	void QTildeShowerHandler::doinit() {
	ShowerHandler::doinit();
	// interactions may have been changed through a setup file so we
	// clear it up here
	// calculate max no of FSR vetos
	_maxFailFSR = max(int(_maxFailFSR), int(_fracFSR*double(generator()->N())));
	// check on the reweighting
	for(unsigned int ix=0;ix<_fullShowerVetoes.size();++ix) {
	if(_fullShowerVetoes[ix]->behaviour()==1) {
	_reWeight = true;
	break;
	}
	}
	if(_reWeight && maximumTries()<_nReWeight) {
	throw Exception() << "Reweight being performed in the shower but the number of attempts for the"
	<< "shower is less than that for the reweighting.\n"
	<< "Maximum number of attempt for the shower "
	<< fullName() << ":MaxTry is " << maximumTries() << "\nand for reweighting is "
	<< fullName() << ":NReWeight is " << _nReWeight << "\n"
	<< "we recommend the number of attempts is 10 times the number for reweighting\n"
	<< Exception::runerror;
	}
	}

	void QTildeShowerHandler::generateIntrinsicpT(vector<ShowerProgenitorPtr> particlesToShower) {
	_intrinsic.clear();
	if ( !ipTon() \|\| !doISR() ) return;
	// don't do anything for the moment for secondary scatters
	if( !firstInteraction() ) return;
	// generate intrinsic pT
	for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
	// only consider initial-state particles
	if(particlesToShower[ix]->progenitor()->isFinalState()) continue;
	if(!particlesToShower[ix]->progenitor()->dataPtr()->coloured()) continue;
	Energy ipt;
	if(UseRandom::rnd() > _beta) {
	ipt=_iptrms*sqrt(-log(UseRandom::rnd()));
	}
	else {
	ipt=_gamma*sqrt(pow(1.+sqr(_iptmax/_gamma), UseRandom::rnd())-1.);
	}
	pair<Energy,double> pt = make_pair(ipt,UseRandom::rnd(Constants::twopi));
	_intrinsic[particlesToShower[ix]] = pt;
	}
	}

	void QTildeShowerHandler::setupMaximumScales(const vector<ShowerProgenitorPtr> & p,
	XCPtr xcomb) {
	// let POWHEG events radiate freely
	if(_hardEmission==2&&hardTree()) {
	vector<ShowerProgenitorPtr>::const_iterator ckt = p.begin();
	for (; ckt != p.end(); ckt++) (*ckt)->maxHardPt(Constants::MaxEnergy);
	return;
	}
	// return if no vetos
	if (!restrictPhasespace()) return;
	// find out if hard partonic subprocess.
	bool isPartonic(false);
	map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
	cit = _currenttree->incomingLines().begin();
	Lorentz5Momentum pcm;
	for(; cit!=currentTree()->incomingLines().end(); ++cit) {
	pcm += cit->first->progenitor()->momentum();
	isPartonic \|= cit->first->progenitor()->coloured();
	}
	// find minimum pt from hard process, the maximum pt from all outgoing
	// coloured lines (this is simpler and more general than
	// 2stu/(s^2+t^2+u^2)). Maximum scale for scattering processes will
	// be transverse mass.
	Energy ptmax = generator()->maximumCMEnergy();
	// general case calculate the scale
	if ( !hardScaleIsMuF() \|\| (hardVetoReadOption()&&!firstInteraction()) ) {
	// scattering process
	if(currentTree()->isHard()) {
	assert(xcomb);
	// coloured incoming particles
	if (isPartonic) {
	map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
	cjt = currentTree()->outgoingLines().begin();
	for(; cjt!=currentTree()->outgoingLines().end(); ++cjt) {
	if (cjt->first->progenitor()->coloured())
	ptmax = min(ptmax,cjt->first->progenitor()->momentum().mt());
	}
	}
	if (ptmax == generator()->maximumCMEnergy() ) ptmax = pcm.m();
	if(hardScaleIsMuF()&&hardVetoReadOption()&&
	!firstInteraction()) {
	ptmax=min(ptmax,sqrt(xcomb->lastShowerScale()));
	}
	}
	// decay, incoming() is the decaying particle.
	else {
	ptmax = currentTree()->incomingLines().begin()->first
	->progenitor()->momentum().mass();
	}
	}
	// hepeup.SCALUP is written into the lastXComb by the
	// LesHouchesReader itself - use this by user's choice.
	// Can be more general than this.
	else {
	if(currentTree()->isHard()) {
	assert(xcomb);
	ptmax = sqrt( xcomb->lastShowerScale() );
	}
	else {
	ptmax = currentTree()->incomingLines().begin()->first
	->progenitor()->momentum().mass();
	}
	}
	ptmax *= hardScaleFactor();
	// set maxHardPt for all progenitors. For partonic processes this
	// is now the max pt in the FS, for non-partonic processes or
	// processes with no coloured FS the invariant mass of the IS
	vector<ShowerProgenitorPtr>::const_iterator ckt = p.begin();
	for (; ckt != p.end(); ckt++) (*ckt)->maxHardPt(ptmax);
	}

	void QTildeShowerHandler::setupHardScales(const vector<ShowerProgenitorPtr> & p,
	XCPtr xcomb) {
	if ( hardScaleIsMuF() &&
	(!hardVetoReadOption() \|\| firstInteraction()) ) {
	Energy hardScale = ZERO;
	if(currentTree()->isHard()) {
	assert(xcomb);
	hardScale = sqrt( xcomb->lastShowerScale() );
	}
	else {
	hardScale = currentTree()->incomingLines().begin()->first
	->progenitor()->momentum().mass();
	}
	hardScale *= hardScaleFactor();
	vector<ShowerProgenitorPtr>::const_iterator ckt = p.begin();
	for (; ckt != p.end(); ckt++) (*ckt)->hardScale(hardScale);
	muPt = hardScale;
	}
	}

	void QTildeShowerHandler::showerHardProcess(ShowerTreePtr hard, XCPtr xcomb) {
	_hardme = HwMEBasePtr();
	// extract the matrix element
	tStdXCombPtr lastXC = dynamic_ptr_cast<tStdXCombPtr>(xcomb);
	if(lastXC) {
	_hardme = dynamic_ptr_cast<HwMEBasePtr>(lastXC->matrixElement());
	}
	_decayme = HwDecayerBasePtr();
	// set the current tree
	currentTree(hard);
	hardTree(HardTreePtr());
	// work out the type of event
	currentTree()->xcombPtr(dynamic_ptr_cast<StdXCombPtr>(xcomb));
	currentTree()->identifyEventType();
	checkFlags();
	// generate the showering
	doShowering(true,xcomb);
	}

	RealEmissionProcessPtr QTildeShowerHandler::hardMatrixElementCorrection(bool hard) {
	// set the initial enhancement factors for the soft correction
	_initialenhance = 1.;
	_finalenhance = 1.;
	// see if we can get the correction from the matrix element
	// or decayer
	RealEmissionProcessPtr real;
	if(hard) {
	if(_hardme&&_hardme->hasMECorrection()) {
	_hardme->initializeMECorrection(_currenttree->perturbativeProcess(),
	_initialenhance,_finalenhance);
	if(hardMEC())
	real =
	_hardme->applyHardMatrixElementCorrection(_currenttree->perturbativeProcess());
	}
	}
	else {
	if(_decayme&&_decayme->hasMECorrection()) {
	_decayme->initializeMECorrection(_currenttree->perturbativeProcess(),
	_initialenhance,_finalenhance);
	if(hardMEC())
	real = _decayme->applyHardMatrixElementCorrection(_currenttree->perturbativeProcess());
	}
	}
	return real;
	}

	ShowerParticleVector QTildeShowerHandler::createTimeLikeChildren(tShowerParticlePtr, IdList ids) {
	// Create the ShowerParticle objects for the two children of
	// the emitting particle; set the parent/child relationship
	// if same as definition create particles, otherwise create cc
	ShowerParticleVector children;
	for(unsigned int ix=0;ix<2;++ix) {
	children.push_back(new_ptr(ShowerParticle(ids[ix+1],true)));
	if(children[ix]->id()==_progenitor->id()&&!ids[ix+1]->stable())
	children[ix]->set5Momentum(Lorentz5Momentum(_progenitor->progenitor()->mass()));
	else
	children[ix]->set5Momentum(Lorentz5Momentum(ids[ix+1]->mass()));
	}
	return children;
	}

	bool QTildeShowerHandler::timeLikeShower(tShowerParticlePtr particle,
	ShowerInteraction type,
	Branching fb, bool first) {
	// don't do anything if not needed
	if(_limitEmissions == 1 \|\| hardOnly() \|\|
	( _limitEmissions == 2 && _nfs != 0) \|\|
	( _limitEmissions == 4 && _nfs + _nis != 0) ) {
	if(particle->spinInfo()) particle->spinInfo()->develop();
	return false;
	}
	// too many tries
	if(_nFSR>=_maxTryFSR) {
	++_nFailedFSR;
	// too many failed events
	if(_nFailedFSR>=_maxFailFSR)
	throw Exception() << "Too many events have failed due to too many shower emissions, in\n"
	<< "QTildeShowerHandler::timeLikeShower(). Terminating run\n"
	<< Exception::runerror;
	throw Exception() << "Too many attempted emissions in QTildeShowerHandler::timeLikeShower()\n"
	<< Exception::eventerror;
	}
	// generate the emission
	ShowerParticleVector children;
	int ntry=0;
	// generate the emission
	if(!fb.kinematics)
	fb = selectTimeLikeBranching(particle,type,HardBranchingPtr());
	// no emission, return
	if(!fb.kinematics) {
	if(particle->spinInfo()) particle->spinInfo()->develop();
	return false;
	}
	Branching fc[2];
	bool setupChildren = true;
	while (ntry<50) {
	fc[0] = Branching();
	fc[1] = Branching();
	++ntry;
	assert(fb.kinematics);
	// has emitted
	// Assign the shower kinematics to the emitting particle.
	if(setupChildren) {
	++_nFSR;
	particle->showerKinematics(fb.kinematics);
	// check highest pT
	if(fb.kinematics->pT()>progenitor()->highestpT())
	progenitor()->highestpT(fb.kinematics->pT());
	// create the children
	children = createTimeLikeChildren(particle,fb.ids);
	// update the children
	particle->showerKinematics()->
	updateChildren(particle, children,fb.type,_reconOpt>=3);
	// update number of emissions
	++_nfs;
	if(_limitEmissions!=0) {
	if(children[0]->spinInfo()) children[0]->spinInfo()->develop();
	if(children[1]->spinInfo()) children[1]->spinInfo()->develop();
	if(particle->spinInfo()) particle->spinInfo()->develop();
	return true;
	}
	setupChildren = false;
	}
	// select branchings for children
	fc[0] = selectTimeLikeBranching(children[0],type,HardBranchingPtr());
	fc[1] = selectTimeLikeBranching(children[1],type,HardBranchingPtr());
	// old default
	if(_reconOpt==0) {
	// shower the first particle
	if(fc[0].kinematics) timeLikeShower(children[0],type,fc[0],false);
	if(children[0]->spinInfo()) children[0]->spinInfo()->develop();
	// shower the second particle
	if(fc[1].kinematics) timeLikeShower(children[1],type,fc[1],false);
	if(children[1]->spinInfo()) children[1]->spinInfo()->develop();
	break;
	}
	// Herwig default
	else if(_reconOpt==1) {
	// shower the first particle
	if(fc[0].kinematics) timeLikeShower(children[0],type,fc[0],false);
	if(children[0]->spinInfo()) children[0]->spinInfo()->develop();
	// shower the second particle
	if(fc[1].kinematics) timeLikeShower(children[1],type,fc[1],false);
	if(children[1]->spinInfo()) children[1]->spinInfo()->develop();
	// branching has happened
	particle->showerKinematics()->updateParent(particle, children,fb.type);
	// clean up the vetoed emission
	if(particle->virtualMass()==ZERO) {
	particle->showerKinematics(ShoKinPtr());
	for(unsigned int ix=0;ix<children.size();++ix)
	particle->abandonChild(children[ix]);
	children.clear();
	if(particle->spinInfo()) particle->spinInfo()->decayVertex(VertexPtr());
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	// generate the new emission
	fb = selectTimeLikeBranching(particle,type,HardBranchingPtr());
	// no emission, return
	if(!fb.kinematics) {
	if(particle->spinInfo()) particle->spinInfo()->develop();
	return false;
	}
	setupChildren = true;
	continue;
	}
	else
	break;
	}
	// veto children
	else if(_reconOpt>=2) {
	// cut-off masses for the branching
	const vector<Energy> & virtualMasses = fb.sudakov->virtualMasses(fb.ids);
	// compute the masses of the children
	Energy masses[3];
	for(unsigned int ix=0;ix<2;++ix) {
	if(fc[ix].kinematics) {
	const vector<Energy> & vm = fc[ix].sudakov->virtualMasses(fc[ix].ids);
	Energy2 q2 =
	fc[ix].kinematics->z()(1.-fc[ix].kinematics->z())sqr(fc[ix].kinematics->scale());
	if(fc[ix].ids[0]->id()!=ParticleID::g) q2 += sqr(vm[0]);
	masses[ix+1] = sqrt(q2);
	}
	else {
	masses[ix+1] = virtualMasses[ix+1];
	}
	}
	masses[0] = fb.ids[0]->id()!=ParticleID::g ? virtualMasses[0] : ZERO;
	double z = fb.kinematics->z();
	Energy2 pt2 = z(1.-z)(z(1.-z)sqr(fb.kinematics->scale()) + sqr(masses[0]))
	- sqr(masses[1])(1.-z) - sqr(masses[2])z;
	if(pt2>=ZERO) {
	break;
	}
	else {
	// reset the scales for the children
	for(unsigned int ix=0;ix<2;++ix) {
	if(fc[ix].kinematics)
	children[ix]->vetoEmission(fc[ix].type,fc[ix].kinematics->scale());
	else
	children[ix]->vetoEmission(ShowerPartnerType::QCDColourLine,ZERO);
	children[ix]->virtualMass(ZERO);
	}
	}
	}
	};
	if(_reconOpt>=2) {
	// shower the first particle
	if(fc[0].kinematics) timeLikeShower(children[0],type,fc[0],false);
	if(children[0]->spinInfo()) children[0]->spinInfo()->develop();
	// shower the second particle
	if(fc[1].kinematics) timeLikeShower(children[1],type,fc[1],false);
	if(children[1]->spinInfo()) children[1]->spinInfo()->develop();
	// branching has happened
	particle->showerKinematics()->updateParent(particle, children,fb.type);
	}
	if(first&&!children.empty())
	particle->showerKinematics()->resetChildren(particle,children);
	if(particle->spinInfo()) particle->spinInfo()->develop();
	return true;
	}

	bool
	QTildeShowerHandler::spaceLikeShower(tShowerParticlePtr particle, PPtr beam,
	ShowerInteraction type) {
	//using the pdf's associated with the ShowerHandler assures, that
	//modified pdf's are used for the secondary interactions via
	//CascadeHandler::resetPDFs(...)
	tcPDFPtr pdf;
	if(firstPDF().particle() == _beam)
	pdf = firstPDF().pdf();
	if(secondPDF().particle() == _beam)
	pdf = secondPDF().pdf();
	Energy freeze = pdfFreezingScale();
	// don't do anything if not needed
	if(_limitEmissions == 2 \|\| hardOnly() \|\|
	( _limitEmissions == 1 && _nis != 0 ) \|\|
	( _limitEmissions == 4 && _nis + _nfs != 0 ) ) {
	if(particle->spinInfo()) particle->spinInfo()->develop();
	return false;
	}
	Branching bb;
	// generate branching
	while (true) {
	bb=_splittingGenerator->chooseBackwardBranching(*particle,beam,
	_initialenhance,
	_beam,type,
	pdf,freeze);
	// return if no emission
	if(!bb.kinematics) {
	if(particle->spinInfo()) particle->spinInfo()->develop();
	return false;
	}
	// if not vetoed break
	if(!spaceLikeVetoed(bb,particle)) break;
	// otherwise reset scale and continue
	particle->vetoEmission(bb.type,bb.kinematics->scale());
	if(particle->spinInfo()) particle->spinInfo()->decayVertex(VertexPtr());
	}
	// assign the splitting function and shower kinematics
	particle->showerKinematics(bb.kinematics);
	if(bb.kinematics->pT()>progenitor()->highestpT())
	progenitor()->highestpT(bb.kinematics->pT());
	// For the time being we are considering only 1->2 branching
	// particles as in Sudakov form factor
	tcPDPtr part[2]={bb.ids[0],bb.ids[2]};
	// Now create the actual particles, make the otherChild a final state
	// particle, while the newParent is not
	ShowerParticlePtr newParent = new_ptr(ShowerParticle(part[0],false));
	ShowerParticlePtr otherChild = new_ptr(ShowerParticle(part[1],true,true));
	ShowerParticleVector theChildren;
	theChildren.push_back(particle);
	theChildren.push_back(otherChild);
	//this updates the evolution scale
	particle->showerKinematics()->
	updateParent(newParent, theChildren,bb.type);
	// update the history if needed
	_currenttree->updateInitialStateShowerProduct(_progenitor,newParent);
	_currenttree->addInitialStateBranching(particle,newParent,otherChild);
	// for the reconstruction of kinematics, parent/child
	// relationships are according to the branching process:
	// now continue the shower
	++_nis;
	bool emitted = _limitEmissions==0 ?
	spaceLikeShower(newParent,beam,type) : false;
	if(newParent->spinInfo()) newParent->spinInfo()->develop();
	// now reconstruct the momentum
	if(!emitted) {
	if(_intrinsic.find(_progenitor)==_intrinsic.end()) {
	bb.kinematics->updateLast(newParent,ZERO,ZERO);
	}
	else {
	pair<Energy,double> kt=_intrinsic[_progenitor];
	bb.kinematics->updateLast(newParent,
	kt.first*cos(kt.second),
	kt.first*sin(kt.second));
	}
	}
	particle->showerKinematics()->
	updateChildren(newParent, theChildren,bb.type,_reconOpt>=4);
	if(_limitEmissions!=0) {
	if(particle->spinInfo()) particle->spinInfo()->develop();
	return true;
	}
	// perform the shower of the final-state particle
	timeLikeShower(otherChild,type,Branching(),true);
	updateHistory(otherChild);
	if(theChildren[1]->spinInfo()) theChildren[1]->spinInfo()->develop();
	// return the emitted
	if(particle->spinInfo()) particle->spinInfo()->develop();
	return true;
	}

	void QTildeShowerHandler::showerDecay(ShowerTreePtr decay) {
	// work out the type of event
	currentTree()->xcombPtr(StdXCombPtr());
	currentTree()->identifyEventType();
	_decayme = HwDecayerBasePtr();
	_hardme = HwMEBasePtr();
	// find the decayer
	// try the normal way if possible
	tDMPtr dm = decay->incomingLines().begin()->first->original() ->decayMode();
	if(!dm) dm = decay->incomingLines().begin()->first->copy() ->decayMode();
	if(!dm) dm = decay->incomingLines().begin()->first->progenitor()->decayMode();
	// otherwise make a string and look it up
	if(!dm) {
	string tag = decay->incomingLines().begin()->first->original()->dataPtr()->name()
	+ "->";
	OrderedParticles outgoing;
	for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
	it=decay->outgoingLines().begin();it!=decay->outgoingLines().end();++it) {
	if(abs(decay->incomingLines().begin()->first->original()->id()) == ParticleID::t &&
	abs(it->first->original()->id())==ParticleID::Wplus &&
	decay->treelinks().size() == 1) {
	ShowerTreePtr Wtree = decay->treelinks().begin()->first;
	for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
	it2=Wtree->outgoingLines().begin();it2!=Wtree->outgoingLines().end();++it2) {
	outgoing.insert(it2->first->original()->dataPtr());
	}
	}
	else {
	outgoing.insert(it->first->original()->dataPtr());
	}
	}
	for(OrderedParticles::const_iterator it=outgoing.begin(); it!=outgoing.end();++it) {
	if(it!=outgoing.begin()) tag += ",";
	tag +=(**it).name();
	}
	tag += ";";
	dm = findDecayMode(tag);
	}
	if(dm) _decayme = dynamic_ptr_cast<HwDecayerBasePtr>(dm->decayer());
	// set the ShowerTree to be showered
	currentTree(decay);
	decay->applyTransforms();
	hardTree(HardTreePtr());
	// generate the showering
	doShowering(false,XCPtr());
	// if no vetos
	// force calculation of spin correlations
	SpinPtr spInfo = decay->incomingLines().begin()->first->progenitor()->spinInfo();
	if(spInfo) {
	if(!spInfo->developed()) spInfo->needsUpdate();
	spInfo->develop();
	}
	}

	bool QTildeShowerHandler::spaceLikeDecayShower(tShowerParticlePtr particle,
	const ShowerParticle::EvolutionScales & maxScales,
	Energy minmass,ShowerInteraction type,
	Branching fb) {
	// too many tries
	if(_nFSR>=_maxTryFSR) {
	++_nFailedFSR;
	// too many failed events
	if(_nFailedFSR>=_maxFailFSR)
	throw Exception() << "Too many events have failed due to too many shower emissions, in\n"
	<< "QTildeShowerHandler::timeLikeShower(). Terminating run\n"
	<< Exception::runerror;
	throw Exception() << "Too many attempted emissions in QTildeShowerHandler::timeLikeShower()\n"
	<< Exception::eventerror;
	}
	// generate the emission
	ShowerParticleVector children;
	int ntry=0;
	// generate the emission
	if(!fb.kinematics)
	fb = selectSpaceLikeDecayBranching(particle,maxScales,minmass,type,
	HardBranchingPtr());
	// no emission, return
	if(!fb.kinematics) return false;
	Branching fc[2];
	bool setupChildren = true;
	while (ntry<50) {
	if(particle->virtualMass()==ZERO)
	particle->virtualMass(_progenitor->progenitor()->mass());
	fc[0] = Branching();
	fc[1] = Branching();
	++ntry;
	assert(fb.kinematics);
	// has emitted
	// Assign the shower kinematics to the emitting particle.
	if(setupChildren) {
	++_nFSR;
	// Assign the shower kinematics to the emitting particle.
	particle->showerKinematics(fb.kinematics);
	if(fb.kinematics->pT()>progenitor()->highestpT())
	progenitor()->highestpT(fb.kinematics->pT());
	// create the ShowerParticle objects for the two children
	children = createTimeLikeChildren(particle,fb.ids);
	// updateChildren the children
	particle->showerKinematics()->
	updateChildren(particle, children, fb.type,_reconOpt>=3);
	setupChildren = false;
	}
	// select branchings for children
	fc[0] = selectSpaceLikeDecayBranching(children[0],maxScales,minmass,
	type,HardBranchingPtr());
	fc[1] = selectTimeLikeBranching (children[1],type,HardBranchingPtr());
	// old default
	if(_reconOpt==0) {
	// shower the first particle
	_currenttree->updateInitialStateShowerProduct(_progenitor,children[0]);
	_currenttree->addInitialStateBranching(particle,children[0],children[1]);
	if(fc[0].kinematics) spaceLikeDecayShower(children[0],maxScales,minmass,type,Branching());
	// shower the second particle
	if(fc[1].kinematics) timeLikeShower(children[1],type,fc[1],true);
	updateHistory(children[1]);
	// branching has happened
	break;
	}
	// Herwig default
	else if(_reconOpt==1) {
	// shower the first particle
	_currenttree->updateInitialStateShowerProduct(_progenitor,children[0]);
	_currenttree->addInitialStateBranching(particle,children[0],children[1]);
	if(fc[0].kinematics) spaceLikeDecayShower(children[0],maxScales,minmass,type,Branching());
	// shower the second particle
	if(fc[1].kinematics) timeLikeShower(children[1],type,fc[1],true);
	updateHistory(children[1]);
	// branching has happened
	particle->showerKinematics()->updateParent(particle, children,fb.type);
	// clean up the vetoed emission
	if(particle->virtualMass()==ZERO) {
	particle->showerKinematics(ShoKinPtr());
	for(unsigned int ix=0;ix<children.size();++ix)
	particle->abandonChild(children[ix]);
	children.clear();
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	// generate the new emission
	fb = selectSpaceLikeDecayBranching(particle,maxScales,minmass,type,
	HardBranchingPtr());
	// no emission, return
	if(!fb.kinematics) {
	return false;
	}
	setupChildren = true;
	continue;
	}
	else
	break;
	}
	else if(_reconOpt>=2) {
	// cut-off masses for the branching
	const vector<Energy> & virtualMasses = fb.sudakov->virtualMasses(fb.ids);
	// compute the masses of the children
	Energy masses[3];
	// space-like children
	masses[1] = children[0]->virtualMass();
	// time-like child
	if(fc[1].kinematics) {
	const vector<Energy> & vm = fc[1].sudakov->virtualMasses(fc[1].ids);
	Energy2 q2 =
	fc[1].kinematics->z()(1.-fc[1].kinematics->z())sqr(fc[1].kinematics->scale());
	if(fc[1].ids[0]->id()!=ParticleID::g) q2 += sqr(vm[0]);
	masses[2] = sqrt(q2);
	}
	else {
	masses[2] = virtualMasses[2];
	}
	masses[0]=particle->virtualMass();
	double z = fb.kinematics->z();
	Energy2 pt2 = (1.-z)(zsqr(masses[0])-sqr(masses[1])-z/(1.-z)*sqr(masses[2]));
	if(pt2>=ZERO) {
	break;
	}
	else {
	// reset the scales for the children
	for(unsigned int ix=0;ix<2;++ix) {
	if(fc[ix].kinematics)
	children[ix]->vetoEmission(fc[ix].type,fc[ix].kinematics->scale());
	else {
	if(ix==0)
	children[ix]->vetoEmission(ShowerPartnerType::QCDColourLine,Constants::MaxEnergy);
	else
	children[ix]->vetoEmission(ShowerPartnerType::QCDColourLine,ZERO);
	}
	}
	children[0]->virtualMass(_progenitor->progenitor()->mass());
	children[1]->virtualMass(ZERO);
	}
	}
	};
	if(_reconOpt>=2) {
	// In the case of splittings which involves coloured particles,
	// set properly the colour flow of the branching.
	// update the history if needed
	_currenttree->updateInitialStateShowerProduct(_progenitor,children[0]);
	_currenttree->addInitialStateBranching(particle,children[0],children[1]);
	// shower the first particle
	if(fc[0].kinematics) spaceLikeDecayShower(children[0],maxScales,minmass,type,Branching());
	// shower the second particle
	if(fc[1].kinematics) timeLikeShower(children[1],type,fc[1],true);
	updateHistory(children[1]);
	// branching has happened
	particle->showerKinematics()->updateParent(particle, children,fb.type);
	}
	// branching has happened
	return true;
	}

	vector<ShowerProgenitorPtr> QTildeShowerHandler::setupShower(bool hard) {
	RealEmissionProcessPtr real;
	// generate hard me if needed
	if(_hardEmission==1) {
	real = hardMatrixElementCorrection(hard);
	if(real&&!real->outgoing().empty()) setupMECorrection(real);
	}
	// generate POWHEG hard emission if needed
	else if(_hardEmission==2)
	hardestEmission(hard);
	// set the initial colour partners
	setEvolutionPartners(hard,interaction_,false);
	// get the particles to be showered
	vector<ShowerProgenitorPtr> particlesToShower =
	currentTree()->extractProgenitors();
	// return the answer
	return particlesToShower;
	}

	void QTildeShowerHandler::setEvolutionPartners(bool hard,ShowerInteraction type,
	bool clear) {
	// match the particles in the ShowerTree and hardTree
	if(hardTree() && !hardTree()->connect(currentTree()))
	throw Exception() << "Can't match trees in "
	<< "QTildeShowerHandler::setEvolutionPartners()"
	<< Exception::eventerror;
	// extract the progenitors
	vector<ShowerParticlePtr> particles =
	currentTree()->extractProgenitorParticles();
	// clear the partners if needed
	if(clear) {
	for(unsigned int ix=0;ix<particles.size();++ix) {
	particles[ix]->partner(ShowerParticlePtr());
	particles[ix]->clearPartners();
	}
	}
	// sort out the colour partners
	if(hardTree()) {
	// find the partner
	for(unsigned int ix=0;ix<particles.size();++ix) {
	tShowerParticlePtr partner = hardTree()->particles()[particles[ix]]->branchingParticle()->partner();
	if(!partner) continue;
	for(map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
	it=hardTree()->particles().begin();
	it!=hardTree()->particles().end();++it) {
	if(it->second->branchingParticle()==partner) {
	particles[ix]->partner(it->first);
	break;
	}
	}
	if(!particles[ix]->partner())
	throw Exception() << "Can't match partners in "
	<< "QTildeShowerHandler::setEvolutionPartners()"
	<< Exception::eventerror;
	}
	}

	// Set the initial evolution scales
	showerModel()->partnerFinder()->
	setInitialEvolutionScales(particles,!hard,interaction_,!_hardtree);
	if(hardTree() && _hardPOWHEG) {
	bool tooHard=false;
	map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
	eit=hardTree()->particles().end();
	for(unsigned int ix=0;ix<particles.size();++ix) {
	map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
	mit = hardTree()->particles().find(particles[ix]);
	Energy hardScale(ZERO);
	ShowerPartnerType type(ShowerPartnerType::Undefined);
	// final-state
	if(particles[ix]->isFinalState()) {
	if(mit!= eit && !mit->second->children().empty()) {
	hardScale = mit->second->scale();
	type = mit->second->type();
	}
	}
	// initial-state
	else {
	if(mit!= eit && mit->second->parent()) {
	hardScale = mit->second->parent()->scale();
	type = mit->second->parent()->type();
	}
	}
	if(type!=ShowerPartnerType::Undefined) {
	if(type==ShowerPartnerType::QED) {
	tooHard \|= particles[ix]->scales().QED_noAO<hardScale;
	}
	else if(type==ShowerPartnerType::QCDColourLine) {
	tooHard \|= particles[ix]->scales().QCD_c_noAO<hardScale;
	}
	else if(type==ShowerPartnerType::QCDAntiColourLine) {
	tooHard \|= particles[ix]->scales().QCD_ac_noAO<hardScale;
	}
	}
	}
	if(tooHard) convertHardTree(hard,type);
	}
	}

	void QTildeShowerHandler::updateHistory(tShowerParticlePtr particle) {
	if(!particle->children().empty()) {
	ShowerParticleVector theChildren;
	for(unsigned int ix=0;ix<particle->children().size();++ix) {
	ShowerParticlePtr part = dynamic_ptr_cast<ShowerParticlePtr>
	(particle->children()[ix]);
	theChildren.push_back(part);
	}
	// update the history if needed
	if(particle==_currenttree->getFinalStateShowerProduct(_progenitor))
	_currenttree->updateFinalStateShowerProduct(_progenitor,
	particle,theChildren);
	_currenttree->addFinalStateBranching(particle,theChildren);
	for(unsigned int ix=0;ix<theChildren.size();++ix)
	updateHistory(theChildren[ix]);
	}
	}

	bool QTildeShowerHandler::startTimeLikeShower(ShowerInteraction type) {
	_nFSR = 0;
	// initialize basis vectors etc
	if(!progenitor()->progenitor()->partner()) return false;
	progenitor()->progenitor()->initializeFinalState();
	if(hardTree()) {
	map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
	eit=hardTree()->particles().end(),
	mit = hardTree()->particles().find(progenitor()->progenitor());
	if( mit != eit && !mit->second->children().empty() ) {
	bool output=truncatedTimeLikeShower(progenitor()->progenitor(),
	mit->second ,type,Branching(),true);
	if(output) updateHistory(progenitor()->progenitor());
	return output;
	}
	}
	// do the shower
	bool output = hardOnly() ? false :
	timeLikeShower(progenitor()->progenitor() ,type,Branching(),true) ;
	if(output) updateHistory(progenitor()->progenitor());
	return output;
	}

	bool QTildeShowerHandler::startSpaceLikeShower(PPtr parent, ShowerInteraction type) {
	// initialise the basis vectors
	if(!progenitor()->progenitor()->partner()) return false;
	progenitor()->progenitor()->initializeInitialState(parent);
	if(hardTree()) {
	map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
	eit =hardTree()->particles().end(),
	mit = hardTree()->particles().find(progenitor()->progenitor());
	if( mit != eit && mit->second->parent() ) {
	return truncatedSpaceLikeShower( progenitor()->progenitor(),
	parent, mit->second->parent(), type );
	}
	}
	// perform the shower
	return hardOnly() ? false :
	spaceLikeShower(progenitor()->progenitor(),parent,type);
	}

	bool QTildeShowerHandler::
	startSpaceLikeDecayShower(const ShowerParticle::EvolutionScales & maxScales,
	Energy minimumMass,ShowerInteraction type) {
	_nFSR = 0;
	// set up the particle basis vectors
	if(!progenitor()->progenitor()->partner()) return false;
	progenitor()->progenitor()->initializeDecay();
	if(hardTree()) {
	map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
	eit =hardTree()->particles().end(),
	mit = hardTree()->particles().find(progenitor()->progenitor());
	if( mit != eit && mit->second->parent() ) {
	HardBranchingPtr branch=mit->second;
	while(branch->parent()) branch=branch->parent();
	return truncatedSpaceLikeDecayShower(progenitor()->progenitor(),maxScales,
	minimumMass, branch ,type, Branching());
	}
	}
	// perform the shower
	return hardOnly() ? false :
	spaceLikeDecayShower(progenitor()->progenitor(),maxScales,minimumMass,type,Branching());
	}

	bool QTildeShowerHandler::timeLikeVetoed(const Branching & fb,
	ShowerParticlePtr particle) {
	// work out type of interaction
	ShowerInteraction type = convertInteraction(fb.type);
	// check whether emission was harder than largest pt of hard subprocess
	if ( restrictPhasespace() && fb.kinematics->pT() > _progenitor->maxHardPt() )
	return true;
	// soft matrix element correction veto
	if( softMEC()) {
	if(_hardme && _hardme->hasMECorrection()) {
	if(_hardme->softMatrixElementVeto(_progenitor,particle,fb))
	return true;
	}
	else if(_decayme && _decayme->hasMECorrection()) {
	if(_decayme->softMatrixElementVeto(_progenitor,particle,fb))
	return true;
	}
	}
	// veto on maximum pt
	if(fb.kinematics->pT()>_progenitor->maximumpT(type)) return true;
	// general vetos
	if (fb.kinematics && !_vetoes.empty()) {
	bool vetoed=false;
	for (vector<ShowerVetoPtr>::iterator v = _vetoes.begin();
	v != _vetoes.end(); ++v) {
	bool test = (**v).vetoTimeLike(_progenitor,particle,fb);
	switch((**v).vetoType()) {
	case ShowerVeto::Emission:
	vetoed \|= test;
	break;
	case ShowerVeto::Shower:
	if(test) throw VetoShower();
	break;
	case ShowerVeto::Event:
	if(test) throw Veto();
	break;
	}
	}
	if(vetoed) return true;
	}
	if ( firstInteraction() &&
	profileScales() ) {
	double weight =
	profileScales()->
	hardScaleProfile(_progenitor->hardScale(),fb.kinematics->pT());
	if ( UseRandom::rnd() > weight )
	return true;
	}
	return false;
	}

	bool QTildeShowerHandler::spaceLikeVetoed(const Branching & bb,
	ShowerParticlePtr particle) {
	// work out type of interaction
	ShowerInteraction type = convertInteraction(bb.type);
	// check whether emission was harder than largest pt of hard subprocess
	if (restrictPhasespace() && bb.kinematics->pT() > _progenitor->maxHardPt())
	return true;
	// apply the soft correction
	if( softMEC() && _hardme && _hardme->hasMECorrection() ) {
	if(_hardme->softMatrixElementVeto(_progenitor,particle,bb))
	return true;
	}
	// the more general vetos

	// check vs max pt for the shower
	if(bb.kinematics->pT()>_progenitor->maximumpT(type)) return true;

	if (!_vetoes.empty()) {
	bool vetoed=false;
	for (vector<ShowerVetoPtr>::iterator v = _vetoes.begin();
	v != _vetoes.end(); ++v) {
	bool test = (**v).vetoSpaceLike(_progenitor,particle,bb);
	switch ((**v).vetoType()) {
	case ShowerVeto::Emission:
	vetoed \|= test;
	break;
	case ShowerVeto::Shower:
	if(test) throw VetoShower();
	break;
	case ShowerVeto::Event:
	if(test) throw Veto();
	break;
	}
	}
	if (vetoed) return true;
	}
	if ( firstInteraction() &&
	profileScales() ) {
	double weight =
	profileScales()->
	hardScaleProfile(_progenitor->hardScale(),bb.kinematics->pT());
	if ( UseRandom::rnd() > weight )
	return true;
	}
	return false;
	}

	bool QTildeShowerHandler::spaceLikeDecayVetoed( const Branching & fb,
	ShowerParticlePtr particle) {
	// work out type of interaction
	ShowerInteraction type = convertInteraction(fb.type);
	// apply the soft correction
	if( softMEC() && _decayme && _decayme->hasMECorrection() ) {
	if(_decayme->softMatrixElementVeto(_progenitor,particle,fb))
	return true;
	}
	// veto on hardest pt in the shower
	if(fb.kinematics->pT()> _progenitor->maximumpT(type)) return true;
	// general vetos
	if (!_vetoes.empty()) {
	bool vetoed=false;
	for (vector<ShowerVetoPtr>::iterator v = _vetoes.begin();
	v != _vetoes.end(); ++v) {
	bool test = (**v).vetoSpaceLike(_progenitor,particle,fb);
	switch((**v).vetoType()) {
	case ShowerVeto::Emission:
	vetoed \|= test;
	break;
	case ShowerVeto::Shower:
	if(test) throw VetoShower();
	break;
	case ShowerVeto::Event:
	if(test) throw Veto();
	break;
	}
	if (vetoed) return true;
	}
	}
	return false;
	}

	void QTildeShowerHandler::hardestEmission(bool hard) {
	HardTreePtr ISRTree;
	// internal POWHEG in production or decay
	if( (( _hardme && _hardme->hasPOWHEGCorrection()!=0 ) \|\|
	( _decayme && _decayme->hasPOWHEGCorrection()!=0 ) ) ) {
	RealEmissionProcessPtr real;
	unsigned int type(0);
	// production
	if(_hardme) {
	assert(hard);
	real = _hardme->generateHardest( currentTree()->perturbativeProcess(),
	interaction_);
	type = _hardme->hasPOWHEGCorrection();
	}
	// decay
	else {
	assert(!hard);
	real = _decayme->generateHardest( currentTree()->perturbativeProcess() );
	type = _decayme->hasPOWHEGCorrection();
	}
	if(real) {
	// set up ther hard tree
	if(!real->outgoing().empty()) _hardtree = new_ptr(HardTree(real));
	// set up the vetos
	currentTree()->setVetoes(real->pT(),type);
	}
	// store initial state POWHEG radiation
	if(_hardtree && _hardme && _hardme->hasPOWHEGCorrection()==1)
	ISRTree = _hardtree;
	}
	else if (hard) {
	// Get minimum pT cutoff used in shower approximation
	Energy maxpt = 1.*GeV;

	if ( currentTree()->showerApproximation() ) {
	int colouredIn = 0;
	int colouredOut = 0;
	for( map< ShowerProgenitorPtr, tShowerParticlePtr >::iterator it
	= currentTree()->outgoingLines().begin();
	it != currentTree()->outgoingLines().end(); ++it ) {
	if( it->second->coloured() ) ++colouredOut;
	}
	for( map< ShowerProgenitorPtr, ShowerParticlePtr >::iterator it
	= currentTree()->incomingLines().begin();
	it != currentTree()->incomingLines().end(); ++it ) {
	if( it->second->coloured() ) ++colouredIn;
	}
	if ( currentTree()->showerApproximation()->ffPtCut() == currentTree()->showerApproximation()->fiPtCut() &&
	currentTree()->showerApproximation()->ffPtCut() == currentTree()->showerApproximation()->iiPtCut() )
	maxpt = currentTree()->showerApproximation()->ffPtCut();
	else if ( colouredIn == 2 && colouredOut == 0 )
	maxpt = currentTree()->showerApproximation()->iiPtCut();
	else if ( colouredIn == 0 && colouredOut > 1 )
	maxpt = currentTree()->showerApproximation()->ffPtCut();
	else if ( colouredIn == 2 && colouredOut == 1 )
	maxpt = min(currentTree()->showerApproximation()->iiPtCut(), currentTree()->showerApproximation()->fiPtCut());
	else if ( colouredIn == 1 && colouredOut > 1 )
	maxpt = min(currentTree()->showerApproximation()->ffPtCut(), currentTree()->showerApproximation()->fiPtCut());
	else
	maxpt = min(min(currentTree()->showerApproximation()->iiPtCut(), currentTree()->showerApproximation()->fiPtCut()),
	currentTree()->showerApproximation()->ffPtCut());
	}

	// Generate hardtree from born and real emission subprocesses
	_hardtree = generateCKKW(currentTree());

	// Find transverse momentum of hardest emission
	if (_hardtree){
	for(set<HardBranchingPtr>::iterator it=_hardtree->branchings().begin();
	it!=_hardtree->branchings().end();++it) {
	if ((it)->parent() && (it)->status()==HardBranching::Incoming)
	maxpt=(*it)->branchingParticle()->momentum().perp();
	if ((it)->children().size()==2 && (it)->status()==HardBranching::Outgoing){
	if ((*it)->branchingParticle()->id()!=21 &&
	abs((*it)->branchingParticle()->id())>5 ){
	if ((*it)->children()[0]->branchingParticle()->id()==21 \|\|
	abs((*it)->children()[0]->branchingParticle()->id())<6)
	maxpt=(*it)->children()[0]->branchingParticle()->momentum().perp();
	else if ((*it)->children()[1]->branchingParticle()->id()==21 \|\|
	abs((*it)->children()[1]->branchingParticle()->id())<6)
	maxpt=(*it)->children()[1]->branchingParticle()->momentum().perp();
	}
	else {
	if ( abs((*it)->branchingParticle()->id())<6){
	if (abs((*it)->children()[0]->branchingParticle()->id())<6)
	maxpt = (*it)->children()[1]->branchingParticle()->momentum().perp();
	else
	maxpt = (*it)->children()[0]->branchingParticle()->momentum().perp();
	}
	else maxpt = (*it)->children()[1]->branchingParticle()->momentum().perp();
	}
	}
	}
	}


	// Hardest (pt) emission should be the first powheg emission.
	maxpt=min(sqrt(lastXCombPtr()->lastShowerScale()),maxpt);

	// set maximum pT for subsequent emissions from S events
	if ( currentTree()->isPowhegSEvent() ) {
	for( map< ShowerProgenitorPtr, tShowerParticlePtr >::iterator it
	= currentTree()->outgoingLines().begin();
	it != currentTree()->outgoingLines().end(); ++it ) {
	if( ! it->second->coloured() ) continue;
	it->first->maximumpT(maxpt, ShowerInteraction::QCD );
	}
	for( map< ShowerProgenitorPtr, ShowerParticlePtr >::iterator it
	= currentTree()->incomingLines().begin();
	it != currentTree()->incomingLines().end(); ++it ) {
	if( ! it->second->coloured() ) continue;
	it->first->maximumpT(maxpt, ShowerInteraction::QCD );
	}
	}
	}
	else
	_hardtree = generateCKKW(currentTree());

	// if hard me doesn't have a FSR powheg
	// correction use decay powheg correction
	if (_hardme && _hardme->hasPOWHEGCorrection()<2) {
	addFSRUsingDecayPOWHEG(ISRTree);
	}
	// connect the trees
	if(_hardtree) {
	connectTrees(currentTree(),_hardtree,hard);
	}
	}

	void QTildeShowerHandler::addFSRUsingDecayPOWHEG(HardTreePtr ISRTree) {
	// check for intermediate colour singlet resonance
	const ParticleVector inter = _hardme->subProcess()->intermediates();
	if (inter.size()!=1 \|\| inter[0]->momentum().m2()/GeV2 < 0 \|\|
	inter[0]->dataPtr()->iColour()!=PDT::Colour0) {
	return;
	}

	// ignore cases where outgoing particles are not coloured
	map<ShowerProgenitorPtr, tShowerParticlePtr > out = currentTree()->outgoingLines();
	if (out.size() != 2 \|\|
	out. begin()->second->dataPtr()->iColour()==PDT::Colour0 \|\|
	out.rbegin()->second->dataPtr()->iColour()==PDT::Colour0) {
	return;
	}

	// look up decay mode
	tDMPtr dm;
	string tag;
	string inParticle = inter[0]->dataPtr()->name() + "->";
	vector<string> outParticles;
	outParticles.push_back(out.begin ()->first->progenitor()->dataPtr()->name());
	outParticles.push_back(out.rbegin()->first->progenitor()->dataPtr()->name());
	for (int it=0; it<2; ++it){
	tag = inParticle + outParticles[it] + "," + outParticles[(it+1)%2] + ";";
	dm = generator()->findDecayMode(tag);
	if(dm) break;
	}

	// get the decayer
	HwDecayerBasePtr decayer;
	if(dm) decayer = dynamic_ptr_cast<HwDecayerBasePtr>(dm->decayer());
	// check if decayer has a FSR POWHEG correction
	if (!decayer \|\| decayer->hasPOWHEGCorrection()<2) {
	return;
	}
	// generate the hardest emission
	// create RealEmissionProcess
	PPtr in = new_ptr(*inter[0]);
	RealEmissionProcessPtr newProcess(new_ptr(RealEmissionProcess()));
	newProcess->bornIncoming().push_back(in);
	newProcess->bornOutgoing().push_back(out.begin ()->first->progenitor());
	newProcess->bornOutgoing().push_back(out.rbegin()->first->progenitor());
	// generate the FSR
	newProcess = decayer->generateHardest(newProcess);
	HardTreePtr FSRTree;
	if(newProcess) {
	// set up ther hard tree
	if(!newProcess->outgoing().empty()) FSRTree = new_ptr(HardTree(newProcess));
	// set up the vetos
	currentTree()->setVetoes(newProcess->pT(),2);
	}
	if(!FSRTree) return;

	// if there is no ISRTree make _hardtree from FSRTree
	if (!ISRTree){
	vector<HardBranchingPtr> inBranch,hardBranch;
	for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
	cit =currentTree()->incomingLines().begin();
	cit!=currentTree()->incomingLines().end();++cit ) {
	inBranch.push_back(new_ptr(HardBranching(cit->second,SudakovPtr(),
	HardBranchingPtr(),
	HardBranching::Incoming)));
	inBranch.back()->beam(cit->first->original()->parents()[0]);
	hardBranch.push_back(inBranch.back());
	}
	if(inBranch[0]->branchingParticle()->dataPtr()->coloured()) {
	inBranch[0]->colourPartner(inBranch[1]);
	inBranch[1]->colourPartner(inBranch[0]);
	}
	for(set<HardBranchingPtr>::iterator it=FSRTree->branchings().begin();
	it!=FSRTree->branchings().end();++it) {
	if((**it).branchingParticle()->id()!=in->id())
	hardBranch.push_back(*it);
	}
	hardBranch[2]->colourPartner(hardBranch[3]);
	hardBranch[3]->colourPartner(hardBranch[2]);
	HardTreePtr newTree = new_ptr(HardTree(hardBranch,inBranch,
	ShowerInteraction::QCD));
	_hardtree = newTree;
	}

	// Otherwise modify the ISRTree to include the emission in FSRTree
	else {
	vector<tShowerParticlePtr> FSROut, ISROut;
	set<HardBranchingPtr>::iterator itFSR, itISR;
	// get outgoing particles
	for(itFSR =FSRTree->branchings().begin();
	itFSR!=FSRTree->branchings().end();++itFSR){
	if ((**itFSR).status()==HardBranching::Outgoing)
	FSROut.push_back((*itFSR)->branchingParticle());
	}
	for(itISR =ISRTree->branchings().begin();
	itISR!=ISRTree->branchings().end();++itISR){
	if ((**itISR).status()==HardBranching::Outgoing)
	ISROut.push_back((*itISR)->branchingParticle());
	}

	// find COM frame formed by outgoing particles
	LorentzRotation eventFrameFSR, eventFrameISR;
	eventFrameFSR = ((FSROut[0]->momentum()+FSROut[1]->momentum()).findBoostToCM());
	eventFrameISR = ((ISROut[0]->momentum()+ISROut[1]->momentum()).findBoostToCM());

	// find rotation between ISR and FSR frames
	int j=0;
	if (ISROut[0]->id()!=FSROut[0]->id()) j=1;
	eventFrameISR.rotateZ( (eventFrameFSR*FSROut[0]->momentum()).phi()-
	(eventFrameISR*ISROut[j]->momentum()).phi() );
	eventFrameISR.rotateY( (eventFrameFSR*FSROut[0]->momentum()).theta()-
	(eventFrameISR*ISROut[j]->momentum()).theta() );
	eventFrameISR.invert();

	for (itFSR=FSRTree->branchings().begin();
	itFSR!=FSRTree->branchings().end();++itFSR){
	if ((**itFSR).branchingParticle()->id()==in->id()) continue;
	for (itISR =ISRTree->branchings().begin();
	itISR!=ISRTree->branchings().end();++itISR){
	if ((**itISR).status()==HardBranching::Incoming) continue;
	if ((**itFSR).branchingParticle()->id()==
	(**itISR).branchingParticle()->id()){
	// rotate FSRTree particle to ISRTree event frame
	(*itISR).branchingParticle()->setMomentum(eventFrameISR
	eventFrameFSR*
	(**itFSR).branchingParticle()->momentum());
	(**itISR).branchingParticle()->rescaleMass();
	// add the children of the FSRTree particles to the ISRTree
	if(!(**itFSR).children().empty()){
	(itISR).addChild((itFSR).children()[0]);
	(itISR).addChild((itFSR).children()[1]);
	// rotate momenta to ISRTree event frame
	(*itISR).children()[0]->branchingParticle()->setMomentum(eventFrameISR
	eventFrameFSR*
	(**itFSR).children()[0]->branchingParticle()->momentum());
	(*itISR).children()[1]->branchingParticle()->setMomentum(eventFrameISR
	eventFrameFSR*
	(**itFSR).children()[1]->branchingParticle()->momentum());
	}
	}
	}
	}
	_hardtree = ISRTree;
	}
	}

	bool QTildeShowerHandler::truncatedTimeLikeShower(tShowerParticlePtr particle,
	HardBranchingPtr branch,
	ShowerInteraction type,
	Branching fb, bool first) {
	// select a branching if we don't have one
	if(!fb.kinematics)
	fb = selectTimeLikeBranching(particle,type,branch);
	// must be an emission, the forced one it not a truncated one
	assert(fb.kinematics);
	ShowerParticleVector children;
	int ntry=0;
	Branching fc[2];
	bool setupChildren = true;
	while (ntry<50) {
	if(!fc[0].hard) fc[0] = Branching();
	if(!fc[1].hard) fc[1] = Branching();
	++ntry;
	// Assign the shower kinematics to the emitting particle.
	if(setupChildren) {
	++_nFSR;
	// Assign the shower kinematics to the emitting particle.
	particle->showerKinematics(fb.kinematics);
	if(fb.kinematics->pT()>progenitor()->highestpT())
	progenitor()->highestpT(fb.kinematics->pT());
	// create the children
	children = createTimeLikeChildren(particle,fb.ids);
	// update the children
	particle->showerKinematics()->
	updateChildren(particle, children,fb.type,_reconOpt>=3);
	setupChildren = false;
	}
	// select branchings for children
	if(!fc[0].kinematics) {
	// select branching for first particle
	if(!fb.hard && fb.iout ==1 )
	fc[0] = selectTimeLikeBranching(children[0],type,branch);
	else if(fb.hard && !branch->children()[0]->children().empty() )
	fc[0] = selectTimeLikeBranching(children[0],type,branch->children()[0]);
	else
	fc[0] = selectTimeLikeBranching(children[0],type,HardBranchingPtr());
	}
	// select branching for the second particle
	if(!fc[1].kinematics) {
	// select branching for first particle
	if(!fb.hard && fb.iout ==2 )
	fc[1] = selectTimeLikeBranching(children[1],type,branch);
	else if(fb.hard && !branch->children()[1]->children().empty() )
	fc[1] = selectTimeLikeBranching(children[1],type,branch->children()[1]);
	else
	fc[1] = selectTimeLikeBranching(children[1],type,HardBranchingPtr());
	}
	// old default
	if(_reconOpt==0 \|\| (_reconOpt==1 && fb.hard) ) {
	// shower the first particle
	if(fc[0].kinematics) {
	// the parent has truncated emission and following line
	if(!fb.hard && fb.iout == 1)
	truncatedTimeLikeShower(children[0],branch,type,fc[0],false);
	// hard emission and subsquent hard emissions
	else if(fb.hard && !branch->children()[0]->children().empty() )
	truncatedTimeLikeShower(children[0],branch->children()[0],type,fc[0],false);
	// normal shower
	else
	timeLikeShower(children[0],type,fc[0],false);
	}
	if(children[0]->spinInfo()) children[0]->spinInfo()->develop();
	// shower the second particle
	if(fc[1].kinematics) {
	// the parent has truncated emission and following line
	if(!fb.hard && fb.iout == 2)
	truncatedTimeLikeShower(children[1],branch,type,fc[1],false);
	// hard emission and subsquent hard emissions
	else if(fb.hard && !branch->children()[1]->children().empty() )
	truncatedTimeLikeShower(children[1],branch->children()[1],type,fc[1],false);
	else
	timeLikeShower(children[1],type,fc[1],false);
	}
	if(children[1]->spinInfo()) children[1]->spinInfo()->develop();
	// branching has happened
	particle->showerKinematics()->updateParent(particle, children,fb.type);
	break;
	}
	// H7 default
	else if(_reconOpt==1) {
	// shower the first particle
	if(fc[0].kinematics) {
	// the parent has truncated emission and following line
	if(!fb.hard && fb.iout == 1)
	truncatedTimeLikeShower(children[0],branch,type,fc[0],false);
	else
	timeLikeShower(children[0],type,fc[0],false);
	}
	if(children[0]->spinInfo()) children[0]->spinInfo()->develop();
	// shower the second particle
	if(fc[1].kinematics) {
	// the parent has truncated emission and following line
	if(!fb.hard && fb.iout == 2)
	truncatedTimeLikeShower(children[1],branch,type,fc[1],false);
	else
	timeLikeShower(children[1],type,fc[1],false);
	}
	if(children[1]->spinInfo()) children[1]->spinInfo()->develop();
	// branching has happened
	particle->showerKinematics()->updateParent(particle, children,fb.type);
	// clean up the vetoed emission
	if(particle->virtualMass()==ZERO) {
	particle->showerKinematics(ShoKinPtr());
	for(unsigned int ix=0;ix<children.size();++ix)
	particle->abandonChild(children[ix]);
	children.clear();
	if(particle->spinInfo()) particle->spinInfo()->decayVertex(VertexPtr());
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	// generate the new emission
	fb = selectTimeLikeBranching(particle,type,branch);
	// must be at least hard emission
	assert(fb.kinematics);
	setupChildren = true;
	continue;
	}
	else
	break;
	}
	else if(_reconOpt>=2) {
	// cut-off masses for the branching
	const vector<Energy> & virtualMasses = fb.sudakov->virtualMasses(fb.ids);
	// compute the masses of the children
	Energy masses[3];
	for(unsigned int ix=0;ix<2;++ix) {
	if(fc[ix].kinematics) {
	const vector<Energy> & vm = fc[ix].sudakov->virtualMasses(fc[ix].ids);
	Energy2 q2 =
	fc[ix].kinematics->z()(1.-fc[ix].kinematics->z())sqr(fc[ix].kinematics->scale());
	if(fc[ix].ids[0]->id()!=ParticleID::g) q2 += sqr(vm[0]);
	masses[ix+1] = sqrt(q2);
	}
	else {
	masses[ix+1] = virtualMasses[ix+1];
	}
	}
	masses[0] = fb.ids[0]->id()!=ParticleID::g ? virtualMasses[0] : ZERO;
	double z = fb.kinematics->z();
	Energy2 pt2 = z(1.-z)(z(1.-z)sqr(fb.kinematics->scale()) + sqr(masses[0]))
	- sqr(masses[1])(1.-z) - sqr(masses[2])z;
	if(pt2>=ZERO) {
	break;
	}
	// if only the hard emission have to accept it
	else if ((fc[0].hard && !fc[1].kinematics) \|\|
	(fc[1].hard && !fc[0].kinematics) ) {
	break;
	}
	else {
	// reset the scales for the children
	for(unsigned int ix=0;ix<2;++ix) {
	if(fc[ix].hard) continue;
	if(fc[ix].kinematics && ! fc[ix].hard )
	children[ix]->vetoEmission(fc[ix].type,fc[ix].kinematics->scale());
	else
	children[ix]->vetoEmission(ShowerPartnerType::QCDColourLine,ZERO);
	children[ix]->virtualMass(ZERO);
	}
	}
	}
	};
	if(_reconOpt>=2) {
	// shower the first particle
	if(fc[0].kinematics) {
	// the parent has truncated emission and following line
	if(!fb.hard && fb.iout == 1)
	truncatedTimeLikeShower(children[0],branch,type,fc[0],false);
	// hard emission and subsquent hard emissions
	else if(fb.hard && !branch->children()[0]->children().empty() )
	truncatedTimeLikeShower(children[0],branch->children()[0],type,fc[0],false);
	// normal shower
	else
	timeLikeShower(children[0],type,fc[0],false);
	}
	if(children[0]->spinInfo()) children[0]->spinInfo()->develop();
	// shower the second particle
	if(fc[1].kinematics) {
	// the parent has truncated emission and following line
	if(!fb.hard && fb.iout == 2)
	truncatedTimeLikeShower(children[1],branch,type,fc[1],false);
	// hard emission and subsquent hard emissions
	else if(fb.hard && !branch->children()[1]->children().empty() )
	truncatedTimeLikeShower(children[1],branch->children()[1],type,fc[1],false);
	else
	timeLikeShower(children[1],type,fc[1],false);
	}
	if(children[1]->spinInfo()) children[1]->spinInfo()->develop();
	// branching has happened
	particle->showerKinematics()->updateParent(particle, children,fb.type);
	}
	if(first&&!children.empty())
	particle->showerKinematics()->resetChildren(particle,children);
	if(particle->spinInfo()) particle->spinInfo()->develop();
	return true;
	}

	bool QTildeShowerHandler::truncatedSpaceLikeShower(tShowerParticlePtr particle, PPtr beam,
	HardBranchingPtr branch,
	ShowerInteraction type) {
	tcPDFPtr pdf;
	if(firstPDF().particle() == beamParticle())
	pdf = firstPDF().pdf();
	if(secondPDF().particle() == beamParticle())
	pdf = secondPDF().pdf();
	Energy freeze = pdfFreezingScale();
	Branching bb;
	// parameters of the force branching
	double z(0.);
	HardBranchingPtr timelike;
	for( unsigned int ix = 0; ix < branch->children().size(); ++ix ) {
	if( branch->children()[ix]->status() ==HardBranching::Outgoing) {
	timelike = branch->children()[ix];
	}
	if( branch->children()[ix]->status() ==HardBranching::Incoming )
	z = branch->children()[ix]->z();
	}
	// generate truncated branching
	tcPDPtr part[2];
	if(z>=0.&&z<=1.) {
	while (true) {
	if( !isTruncatedShowerON() \|\| hardOnly() ) break;
	bb = splittingGenerator()->chooseBackwardBranching( *particle,
	beam, 1., beamParticle(),
	type , pdf,freeze);
	if( !bb.kinematics \|\| bb.kinematics->scale() < branch->scale() ) {
	bb = Branching();
	break;
	}
	// particles as in Sudakov form factor
	part[0] = bb.ids[0];
	part[1] = bb.ids[2];
	double zsplit = bb.kinematics->z();
	// apply the vetos for the truncated shower
	// if doesn't carry most of momentum
	ShowerInteraction type2 = convertInteraction(bb.type);
	if(type2==branch->sudakov()->interactionType() &&
	zsplit < 0.5) {
	particle->vetoEmission(bb.type,bb.kinematics->scale());
	continue;
	}
	// others
	if( part[0]->id() != particle->id() \|\| // if particle changes type
	bb.kinematics->pT() > progenitor()->maximumpT(type2) \|\| // pt veto
	bb.kinematics->scale() < branch->scale()) { // angular ordering veto
	particle->vetoEmission(bb.type,bb.kinematics->scale());
	continue;
	}
	// and those from the base class
	if(spaceLikeVetoed(bb,particle)) {
	particle->vetoEmission(bb.type,bb.kinematics->scale());
	continue;
	}
	break;
	}
	}
	if( !bb.kinematics ) {
	//do the hard emission
	ShoKinPtr kinematics =
	branch->sudakov()->createInitialStateBranching( branch->scale(), z, branch->phi(),
	branch->children()[0]->pT() );
	// assign the splitting function and shower kinematics
	particle->showerKinematics( kinematics );
	if(kinematics->pT()>progenitor()->highestpT())
	progenitor()->highestpT(kinematics->pT());
	// For the time being we are considering only 1->2 branching
	// Now create the actual particles, make the otherChild a final state
	// particle, while the newParent is not
	ShowerParticlePtr newParent =
	new_ptr( ShowerParticle( branch->branchingParticle()->dataPtr(), false ) );
	ShowerParticlePtr otherChild =
	new_ptr( ShowerParticle( timelike->branchingParticle()->dataPtr(),
	true, true ) );
	ShowerParticleVector theChildren;
	theChildren.push_back( particle );
	theChildren.push_back( otherChild );
	particle->showerKinematics()->
	updateParent( newParent, theChildren, branch->type());
	// update the history if needed
	currentTree()->updateInitialStateShowerProduct( progenitor(), newParent );
	currentTree()->addInitialStateBranching( particle, newParent, otherChild );
	// for the reconstruction of kinematics, parent/child
	// relationships are according to the branching process:
	// now continue the shower
	bool emitted=false;
	if(!hardOnly()) {
	if( branch->parent() ) {
	emitted = truncatedSpaceLikeShower( newParent, beam, branch->parent() , type);
	}
	else {
	emitted = spaceLikeShower( newParent, beam , type);
	}
	}
	if( !emitted ) {
	if( intrinsicpT().find( progenitor() ) == intrinsicpT().end() ) {
	kinematics->updateLast( newParent, ZERO, ZERO );
	}
	else {
	pair<Energy,double> kt = intrinsicpT()[progenitor()];
	kinematics->updateLast( newParent,
	kt.first*cos( kt.second ),
	kt.first*sin( kt.second ) );
	}
	}
	particle->showerKinematics()->
	updateChildren( newParent, theChildren,bb.type,false);
	if(hardOnly()) return true;
	// perform the shower of the final-state particle
	if( timelike->children().empty() ) {
	timeLikeShower( otherChild , type,Branching(),true);
	}
	else {
	truncatedTimeLikeShower( otherChild, timelike , type,Branching(), true);
	}
	updateHistory(otherChild);
	// return the emitted
	return true;
	}
	// assign the splitting function and shower kinematics
	particle->showerKinematics( bb.kinematics );
	if(bb.kinematics->pT()>progenitor()->highestpT())
	progenitor()->highestpT(bb.kinematics->pT());
	// For the time being we are considering only 1->2 branching
	// Now create the actual particles, make the otherChild a final state
	// particle, while the newParent is not
	ShowerParticlePtr newParent = new_ptr( ShowerParticle( part[0], false ) );
	ShowerParticlePtr otherChild = new_ptr( ShowerParticle( part[1], true, true ) );
	ShowerParticleVector theChildren;
	theChildren.push_back( particle );
	theChildren.push_back( otherChild );
	particle->showerKinematics()->
	updateParent( newParent, theChildren, bb.type);
	// update the history if needed
	currentTree()->updateInitialStateShowerProduct( progenitor(), newParent );
	currentTree()->addInitialStateBranching( particle, newParent, otherChild );
	// for the reconstruction of kinematics, parent/child
	// relationships are according to the branching process:
	// now continue the shower
	bool emitted = truncatedSpaceLikeShower( newParent, beam, branch,type);
	// now reconstruct the momentum
	if( !emitted ) {
	if( intrinsicpT().find( progenitor() ) == intrinsicpT().end() ) {
	bb.kinematics->updateLast( newParent, ZERO, ZERO );
	}
	else {
	pair<Energy,double> kt = intrinsicpT()[ progenitor() ];
	bb.kinematics->updateLast( newParent,
	kt.first*cos( kt.second ),
	kt.first*sin( kt.second ) );
	}
	}
	particle->showerKinematics()->
	updateChildren( newParent, theChildren, bb.type,false);
	// perform the shower of the final-state particle
	timeLikeShower( otherChild , type,Branching(),true);
	updateHistory(otherChild);
	// return the emitted
	return true;
	}

	bool QTildeShowerHandler::
	truncatedSpaceLikeDecayShower(tShowerParticlePtr particle,
	const ShowerParticle::EvolutionScales & maxScales,
	Energy minmass, HardBranchingPtr branch,
	ShowerInteraction type, Branching fb) {
	// select a branching if we don't have one
	if(!fb.kinematics)
	fb = selectSpaceLikeDecayBranching(particle,maxScales,minmass,type,branch);
	// must be an emission, the forced one it not a truncated one
	assert(fb.kinematics);
	ShowerParticleVector children;
	int ntry=0;
	Branching fc[2];
	bool setupChildren = true;
	while (ntry<50) {
	if(!fc[0].hard) fc[0] = Branching();
	if(!fc[1].hard) fc[1] = Branching();
	++ntry;
	if(setupChildren) {
	++_nFSR;
	// Assign the shower kinematics to the emitting particle.
	particle->showerKinematics(fb.kinematics);
	if(fb.kinematics->pT()>progenitor()->highestpT())
	progenitor()->highestpT(fb.kinematics->pT());
	// create the ShowerParticle objects for the two children
	children = createTimeLikeChildren(particle,fb.ids);
	// updateChildren the children
	particle->showerKinematics()->
	updateChildren(particle, children, fb.type,_reconOpt>=3);
	setupChildren = false;
	}
	// select branchings for children
	if(!fc[0].kinematics) {
	if(children[0]->id()==particle->id()) {
	// select branching for first particle
	if(!fb.hard)
	fc[0] = selectSpaceLikeDecayBranching(children[0],maxScales,minmass,type,branch);
	else if(fb.hard && ! branch->children()[0]->children().empty() )
	fc[0] = selectSpaceLikeDecayBranching(children[0],maxScales,minmass,type,
	branch->children()[0]);
	else
	fc[0] = selectSpaceLikeDecayBranching(children[0],maxScales,minmass,type,
	HardBranchingPtr());
	}
	else {
	// select branching for first particle
	if(fb.hard && !branch->children()[0]->children().empty() )
	fc[0] = selectTimeLikeBranching(children[0],type,branch->children()[0]);
	else
	fc[0] = selectTimeLikeBranching(children[0],type,HardBranchingPtr());
	}
	}
	// select branching for the second particle
	if(!fc[1].kinematics) {
	if(children[1]->id()==particle->id()) {
	// select branching for first particle
	if(!fb.hard)
	fc[1] = selectSpaceLikeDecayBranching(children[1],maxScales,minmass,type,branch);
	else if(fb.hard && ! branch->children()[1]->children().empty() )
	fc[1] = selectSpaceLikeDecayBranching(children[1],maxScales,minmass,type,
	branch->children()[1]);
	else
	fc[1] = selectSpaceLikeDecayBranching(children[1],maxScales,minmass,type,
	HardBranchingPtr());
	}
	else {
	if(fb.hard && !branch->children()[1]->children().empty() )
	fc[1] = selectTimeLikeBranching(children[1],type,branch->children()[1]);
	else
	fc[1] = selectTimeLikeBranching(children[1],type,HardBranchingPtr());
	}
	}
	// old default
	if(_reconOpt==0 \|\| (_reconOpt==1 && fb.hard) ) {
	// update the history if needed
	currentTree()->updateInitialStateShowerProduct(progenitor(),children[0]);
	currentTree()->addInitialStateBranching(particle,children[0],children[1]);
	// shower the first particle
	if(fc[0].kinematics) {
	if(children[0]->id()==particle->id()) {
	if(!fb.hard)
	truncatedSpaceLikeDecayShower( children[0],maxScales,minmass,
	branch,type,fc[0]);
	else if(fb.hard && ! branch->children()[0]->children().empty() )
	truncatedSpaceLikeDecayShower( children[0],maxScales,minmass,
	branch->children()[0],type,fc[0]);
	else
	spaceLikeDecayShower( children[0],maxScales,minmass,type,fc[0]);
	}
	else {
	if(fb.hard && !branch->children()[0]->children().empty() )
	truncatedTimeLikeShower(children[0],branch->children()[0],type,fc[0],false);
	// normal shower
	else
	timeLikeShower(children[0],type,fc[0],false);
	}
	}
	// shower the second particle
	if(fc[1].kinematics) {
	if(children[0]->id()==particle->id()) {
	if(!fb.hard)
	truncatedSpaceLikeDecayShower( children[0],maxScales,minmass,
	branch,type,fc[1]);
	else if(fb.hard && ! branch->children()[0]->children().empty() )
	truncatedSpaceLikeDecayShower( children[0],maxScales,minmass,
	branch->children()[0],type,fc[1]);
	else
	spaceLikeDecayShower( children[0],maxScales,minmass,type,fc[1]);
	}
	else {
	if(fb.hard && !branch->children()[0]->children().empty() )
	truncatedTimeLikeShower(children[0],branch->children()[0],type,fc[1],false);
	// normal shower
	else
	timeLikeShower(children[0],type,fc[1],false);
	}
	}
	updateHistory(children[1]);
	// branching has happened
	break;
	}
	// H7 default
	else if(_reconOpt==1) {
	// update the history if needed
	currentTree()->updateInitialStateShowerProduct(progenitor(),children[0]);
	currentTree()->addInitialStateBranching(particle,children[0],children[1]);
	// shower the first particle
	if(fc[0].kinematics) {
	if(children[0]->id()==particle->id()) {
	if(!fb.hard)
	truncatedSpaceLikeDecayShower( children[0],maxScales,minmass,
	branch,type,fc[0]);
	else if(fb.hard && ! branch->children()[0]->children().empty() )
	truncatedSpaceLikeDecayShower( children[0],maxScales,minmass,
	branch->children()[0],type,fc[0]);
	else
	spaceLikeDecayShower( children[0],maxScales,minmass,type,fc[0]);
	}
	else {
	if(fb.hard && !branch->children()[0]->children().empty() )
	truncatedTimeLikeShower(children[0],branch->children()[0],type,fc[0],false);
	// normal shower
	else
	timeLikeShower(children[0],type,fc[0],false);
	}
	}
	// shower the second particle
	if(fc[1].kinematics) {
	if(children[0]->id()==particle->id()) {
	if(!fb.hard)
	truncatedSpaceLikeDecayShower( children[0],maxScales,minmass,
	branch,type,fc[1]);
	else if(fb.hard && ! branch->children()[0]->children().empty() )
	truncatedSpaceLikeDecayShower( children[0],maxScales,minmass,
	branch->children()[0],type,fc[1]);
	else
	spaceLikeDecayShower( children[0],maxScales,minmass,type,fc[1]);
	}
	else {
	if(fb.hard && !branch->children()[0]->children().empty() )
	truncatedTimeLikeShower(children[0],branch->children()[0],type,fc[1],false);
	// normal shower
	else
	timeLikeShower(children[0],type,fc[1],false);
	}
	}
	// clean up the vetoed emission
	if(particle->virtualMass()==ZERO) {
	particle->showerKinematics(ShoKinPtr());
	for(unsigned int ix=0;ix<children.size();++ix)
	particle->abandonChild(children[ix]);
	children.clear();
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	// generate the new emission
	fb = selectSpaceLikeDecayBranching(particle,maxScales,minmass,type,branch);
	// must be at least hard emission
	assert(fb.kinematics);
	setupChildren = true;
	continue;
	}
	else {
	updateHistory(children[1]);
	break;
	}
	}
	else if(_reconOpt>=2) {
	// cut-off masses for the branching
	const vector<Energy> & virtualMasses = fb.sudakov->virtualMasses(fb.ids);
	// compute the masses of the children
	Energy masses[3];
	// space-like children
	masses[1] = children[0]->virtualMass();
	// time-like child
	if(fc[1].kinematics) {
	const vector<Energy> & vm = fc[1].sudakov->virtualMasses(fc[1].ids);
	Energy2 q2 =
	fc[1].kinematics->z()(1.-fc[1].kinematics->z())sqr(fc[1].kinematics->scale());
	if(fc[1].ids[0]->id()!=ParticleID::g) q2 += sqr(vm[0]);
	masses[2] = sqrt(q2);
	}
	else {
	masses[2] = virtualMasses[2];
	}
	masses[0]=particle->virtualMass();
	double z = fb.kinematics->z();
	Energy2 pt2 = (1.-z)(zsqr(masses[0])-sqr(masses[1])-z/(1.-z)*sqr(masses[2]));
	if(pt2>=ZERO) {
	break;
	}
	else {
	// reset the scales for the children
	for(unsigned int ix=0;ix<2;++ix) {
	if(fc[ix].kinematics)
	children[ix]->vetoEmission(fc[ix].type,fc[ix].kinematics->scale());
	else {
	if(ix==0)
	children[ix]->vetoEmission(ShowerPartnerType::QCDColourLine,Constants::MaxEnergy);
	else
	children[ix]->vetoEmission(ShowerPartnerType::QCDColourLine,ZERO);
	}
	}
	children[0]->virtualMass(_progenitor->progenitor()->mass());
	children[1]->virtualMass(ZERO);
	}
	}
	};
	if(_reconOpt>=2) {
	// update the history if needed
	currentTree()->updateInitialStateShowerProduct(progenitor(),children[0]);
	currentTree()->addInitialStateBranching(particle,children[0],children[1]);
	// shower the first particle
	if(fc[0].kinematics) {
	if(children[0]->id()==particle->id()) {
	if(!fb.hard)
	truncatedSpaceLikeDecayShower( children[0],maxScales,minmass,
	branch,type,fc[0]);
	else if(fb.hard && ! branch->children()[0]->children().empty() )
	truncatedSpaceLikeDecayShower( children[0],maxScales,minmass,
	branch->children()[0],type,fc[0]);
	else
	spaceLikeDecayShower( children[0],maxScales,minmass,type,fc[0]);
	}
	else {
	if(fb.hard && !branch->children()[0]->children().empty() )
	truncatedTimeLikeShower(children[0],branch->children()[0],type,fc[0],false);
	// normal shower
	else
	timeLikeShower(children[0],type,fc[0],false);
	}
	}
	// shower the second particle
	if(fc[1].kinematics) {
	if(children[0]->id()==particle->id()) {
	if(!fb.hard)
	truncatedSpaceLikeDecayShower( children[0],maxScales,minmass,
	branch,type,fc[1]);
	else if(fb.hard && ! branch->children()[0]->children().empty() )
	truncatedSpaceLikeDecayShower( children[0],maxScales,minmass,
	branch->children()[0],type,fc[1]);
	else
	spaceLikeDecayShower( children[0],maxScales,minmass,type,fc[1]);
	}
	else {
	if(fb.hard && !branch->children()[0]->children().empty() )
	truncatedTimeLikeShower(children[0],branch->children()[0],type,fc[1],false);
	// normal shower
	else
	timeLikeShower(children[0],type,fc[1],false);
	}
	}
	updateHistory(children[1]);
	}
	return true;
	}

	void QTildeShowerHandler::connectTrees(ShowerTreePtr showerTree,
	HardTreePtr hardTree, bool hard ) {
	ShowerParticleVector particles;
	// find the Sudakovs
	for(set<HardBranchingPtr>::iterator cit=hardTree->branchings().begin();
	cit!=hardTree->branchings().end();++cit) {
	// Sudakovs for ISR
	if((cit).parent()&&(cit).status()==HardBranching::Incoming) {
	++_nis;
	vector<long> br(3);
	br[0] = (**cit).parent()->branchingParticle()->id();
	br[1] = (**cit). branchingParticle()->id();
	br[2] = (*cit).parent()->children()[0]==cit ?
	(**cit).parent()->children()[1]->branchingParticle()->id() :
	(**cit).parent()->children()[0]->branchingParticle()->id();
	BranchingList branchings = splittingGenerator()->initialStateBranchings();
	if(br[1]<0&&br[0]==br[1]) {
	br[0] = abs(br[0]);
	br[1] = abs(br[1]);
	}
	else if(br[1]<0) {
	br[1] = -br[1];
	br[2] = -br[2];
	}
	long index = abs(br[1]);
	SudakovPtr sudakov;
	for(BranchingList::const_iterator cjt = branchings.lower_bound(index);
	cjt != branchings.upper_bound(index); ++cjt ) {
	IdList ids = cjt->second.particles;
	if(ids[0]->id()==br[0]&&ids[1]->id()==br[1]&&ids[2]->id()==br[2]) {
	sudakov=cjt->second.sudakov;
	break;
	}
	}
	if(!sudakov) throw Exception() << "Can't find Sudakov for the hard emission in "
	<< "QTildeShowerHandler::connectTrees() for ISR"
	<< Exception::runerror;
	(**cit).parent()->sudakov(sudakov);
	}
	// Sudakovs for FSR
	else if(!(**cit).children().empty()) {
	++_nfs;
	vector<long> br(3);
	br[0] = (**cit) .branchingParticle()->id();
	br[1] = (**cit).children()[0]->branchingParticle()->id();
	br[2] = (**cit).children()[1]->branchingParticle()->id();
	BranchingList branchings = splittingGenerator()->finalStateBranchings();
	if(br[0]<0) {
	br[0] = abs(br[0]);
	br[1] = abs(br[1]);
	br[2] = abs(br[2]);
	}
	long index = br[0];
	SudakovPtr sudakov;
	for(BranchingList::const_iterator cjt = branchings.lower_bound(index);
	cjt != branchings.upper_bound(index); ++cjt ) {
	IdList ids = cjt->second.particles;
	if(ids[0]->id()==br[0]&&ids[1]->id()==br[1]&&ids[2]->id()==br[2]) {
	sudakov=cjt->second.sudakov;
	break;
	}
	}
	if(!sudakov) {
	throw Exception() << "Can't find Sudakov for the hard emission in "
	<< "QTildeShowerHandler::connectTrees()"
	<< Exception::runerror;
	}
	(**cit).sudakov(sudakov);
	}
	}
	// calculate the evolution scale
	for(set<HardBranchingPtr>::iterator cit=hardTree->branchings().begin();
	cit!=hardTree->branchings().end();++cit) {
	particles.push_back((*cit)->branchingParticle());
	}
	showerModel()->partnerFinder()->
	setInitialEvolutionScales(particles,!hard,interaction_,true);
	hardTree->partnersSet(true);
	// inverse reconstruction
	if(hard) {
	showerModel()->kinematicsReconstructor()->
	deconstructHardJets(hardTree,interaction_);
	}
	else
	showerModel()->kinematicsReconstructor()->
	deconstructDecayJets(hardTree,interaction_);
	// now reset the momenta of the showering particles
	vector<ShowerProgenitorPtr> particlesToShower=showerTree->extractProgenitors();
	// match them
	map<ShowerProgenitorPtr,HardBranchingPtr> partners;
	for(set<HardBranchingPtr>::const_iterator bit=hardTree->branchings().begin();
	bit!=hardTree->branchings().end();++bit) {
	Energy2 dmin( 1e30*GeV2 );
	ShowerProgenitorPtr partner;
	for(vector<ShowerProgenitorPtr>::const_iterator pit=particlesToShower.begin();
	pit!=particlesToShower.end();++pit) {
	if(partners.find(*pit)!=partners.end()) continue;
	if( (bit).branchingParticle()->id() != (pit).progenitor()->id() ) continue;
	if( (**bit).branchingParticle()->isFinalState() !=
	(**pit).progenitor()->isFinalState() ) continue;
	if( (**pit).progenitor()->isFinalState() ) {
	Energy2 dtest =
	sqr( (pit).progenitor()->momentum().x() - (bit).showerMomentum().x() ) +
	sqr( (pit).progenitor()->momentum().y() - (bit).showerMomentum().y() ) +
	sqr( (pit).progenitor()->momentum().z() - (bit).showerMomentum().z() ) +
	sqr( (pit).progenitor()->momentum().t() - (bit).showerMomentum().t() );
	// add mass difference for identical particles (e.g. Z0 Z0 production)
	dtest += 1e10sqr((pit).progenitor()->momentum().m()-(*bit).showerMomentum().m());
	if( dtest < dmin ) {
	partner = *pit;
	dmin = dtest;
	}
	}
	else {
	// ensure directions are right
	if((pit).progenitor()->momentum().z()/(bit).showerMomentum().z()>ZERO) {
	partner = *pit;
	break;
	}
	}
	}
	if(!partner) throw Exception() << "Failed to match shower and hard trees in QTildeShowerHandler::hardestEmission"
	<< Exception::eventerror;
	partners[partner] = *bit;
	}
	for(vector<ShowerProgenitorPtr>::const_iterator pit=particlesToShower.begin();
	pit!=particlesToShower.end();++pit) {
	HardBranchingPtr partner = partners[*pit];
	if((**pit).progenitor()->dataPtr()->stable()) {
	(**pit).progenitor()->set5Momentum(partner->showerMomentum());
	(**pit).copy()->set5Momentum(partner->showerMomentum());
	}
	else {
	Lorentz5Momentum oldMomentum = (**pit).progenitor()->momentum();
	Lorentz5Momentum newMomentum = partner->showerMomentum();
	LorentzRotation boost( oldMomentum.findBoostToCM(),oldMomentum.e()/oldMomentum.mass());
	(**pit).progenitor()->transform(boost);
	(**pit).copy() ->transform(boost);
	boost = LorentzRotation(-newMomentum.findBoostToCM(),newMomentum.e()/newMomentum.mass());
	(**pit).progenitor()->transform(boost);
	(**pit).copy() ->transform(boost);
	}
	}
	// correction boosts for daughter trees
	for(map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator
	tit = showerTree->treelinks().begin();
	tit != showerTree->treelinks().end();++tit) {
	ShowerTreePtr decayTree = tit->first;
	map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
	cit = decayTree->incomingLines().begin();
	// reset the momentum of the decay particle
	Lorentz5Momentum oldMomentum = cit->first->progenitor()->momentum();
	Lorentz5Momentum newMomentum = tit->second.second->momentum();
	LorentzRotation boost( oldMomentum.findBoostToCM(),oldMomentum.e()/oldMomentum.mass());
	decayTree->transform(boost,true);
	boost = LorentzRotation(-newMomentum.findBoostToCM(),newMomentum.e()/newMomentum.mass());
	decayTree->transform(boost,true);
	}
	}

	void QTildeShowerHandler::doShowering(bool hard,XCPtr xcomb) {
	// zero number of emissions
	_nis = _nfs = 0;
	// if MC@NLO H event and limited emissions
	// indicate both final and initial state emission
	if ( currentTree()->isMCatNLOHEvent() && _limitEmissions != 0 ) {
	_nis = _nfs = 1;
	}
	// extract particles to shower
	vector<ShowerProgenitorPtr> particlesToShower(setupShower(hard));
	// check if we should shower
	bool colCharge = false;
	for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
	if(particlesToShower[ix]->progenitor()->dataPtr()->coloured() \|\|
	particlesToShower[ix]->progenitor()->dataPtr()->charged()) {
	colCharge = true;
	break;
	}
	}
	if(!colCharge) {
	_currenttree->hasShowered(true);
	return;
	}
	// setup the maximum scales for the shower
	if (restrictPhasespace()) setupMaximumScales(particlesToShower,xcomb);
	// set the hard scales for the profiles
	setupHardScales(particlesToShower,xcomb);
	// specific stuff for hard processes and decays
	Energy minmass(ZERO), mIn(ZERO);
	// hard process generate the intrinsic p_T once and for all
	if(hard) {
	generateIntrinsicpT(particlesToShower);
	}
	// decay compute the minimum mass of the final-state
	else {
	for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
	if(particlesToShower[ix]->progenitor()->isFinalState()) {
	if(particlesToShower[ix]->progenitor()->dataPtr()->stable())
	minmass += particlesToShower[ix]->progenitor()->dataPtr()->constituentMass();
	else
	minmass += particlesToShower[ix]->progenitor()->mass();
	}
	else {
	mIn = particlesToShower[ix]->progenitor()->mass();
	}
	}
	// throw exception if decay can't happen
	if ( minmass > mIn ) {
	throw Exception() << "QTildeShowerHandler.cc: Mass of decaying particle is "
	<< "below constituent masses of decay products."
	<< Exception::eventerror;
	}
	}
	// setup for reweighted
	bool reWeighting = _reWeight && hard && ShowerHandler::currentHandler()->firstInteraction();
	double eventWeight=0.;
	unsigned int nTryReWeight(0);
	// create random particle vector (only need to do once)
	vector<ShowerProgenitorPtr> tmp;
	unsigned int nColouredIncoming = 0;
	while(particlesToShower.size()>0){
	unsigned int xx=UseRandom::irnd(particlesToShower.size());
	tmp.push_back(particlesToShower[xx]);
	particlesToShower.erase(particlesToShower.begin()+xx);
	}
	particlesToShower=tmp;
	for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
	if(!particlesToShower[ix]->progenitor()->isFinalState() &&
	particlesToShower[ix]->progenitor()->coloured()) ++nColouredIncoming;
	}
	bool switchRecon = hard && nColouredIncoming !=1;
	// main shower loop
	unsigned int ntry(0);
	bool reconstructed = false;
	do {
	// clear results of last attempt if needed
	if(ntry!=0) {
	currentTree()->clear();
	setEvolutionPartners(hard,interaction_,true);
	_nis = _nfs = 0;
	// if MC@NLO H event and limited emissions
	// indicate both final and initial state emission
	if ( currentTree()->isMCatNLOHEvent() && _limitEmissions != 0 ) {
	_nis = _nfs = 1;
	}
	for(unsigned int ix=0; ix<particlesToShower.size();++ix) {
	SpinPtr spin = particlesToShower[ix]->progenitor()->spinInfo();
	if(spin && spin->decayVertex() &&
	dynamic_ptr_cast<tcSVertexPtr>(spin->decayVertex())) {
	spin->decayVertex(VertexPtr());
	}
	}
	}
	// loop over particles
	for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
	// extract the progenitor
	progenitor(particlesToShower[ix]);
	// final-state radiation
	if(progenitor()->progenitor()->isFinalState()) {
	if(!doFSR()) continue;
	// perform shower
	progenitor()->hasEmitted(startTimeLikeShower(interaction_));
	}
	// initial-state radiation
	else {
	if(!doISR()) continue;
	// hard process
	if(hard) {
	// get the PDF
	setBeamParticle(_progenitor->beam());
	assert(beamParticle());
	// perform the shower
	// set the beam particle
	tPPtr beamparticle=progenitor()->original();
	if(!beamparticle->parents().empty())
	beamparticle=beamparticle->parents()[0];
	// generate the shower
	progenitor()->hasEmitted(startSpaceLikeShower(beamparticle,
	interaction_));
	}
	// decay
	else {
	// skip colour and electrically neutral particles
	if(!progenitor()->progenitor()->dataPtr()->coloured() &&
	!progenitor()->progenitor()->dataPtr()->charged()) {
	progenitor()->hasEmitted(false);
	continue;
	}
	// perform shower
	// set the scales correctly. The current scale is the maximum scale for
	// emission not the starting scale
	ShowerParticle::EvolutionScales maxScales(progenitor()->progenitor()->scales());
	progenitor()->progenitor()->scales() = ShowerParticle::EvolutionScales();
	if(progenitor()->progenitor()->dataPtr()->charged()) {
	progenitor()->progenitor()->scales().QED = progenitor()->progenitor()->mass();
	progenitor()->progenitor()->scales().QED_noAO = progenitor()->progenitor()->mass();
	}
	if(progenitor()->progenitor()->hasColour()) {
	progenitor()->progenitor()->scales().QCD_c = progenitor()->progenitor()->mass();
	progenitor()->progenitor()->scales().QCD_c_noAO = progenitor()->progenitor()->mass();
	}
	if(progenitor()->progenitor()->hasAntiColour()) {
	progenitor()->progenitor()->scales().QCD_ac = progenitor()->progenitor()->mass();
	progenitor()->progenitor()->scales().QCD_ac_noAO = progenitor()->progenitor()->mass();
	}
	// perform the shower
	progenitor()->hasEmitted(startSpaceLikeDecayShower(maxScales,minmass,
	interaction_));
	}
	}
	}
	// do the kinematic reconstruction, checking if it worked
	reconstructed = hard ?
	showerModel()->kinematicsReconstructor()->
	reconstructHardJets (currentTree(),intrinsicpT(),interaction_,
	switchRecon && ntry>maximumTries()/2) :
	showerModel()->kinematicsReconstructor()->
	reconstructDecayJets(currentTree(),interaction_);
	if(!reconstructed) continue;
	// apply vetos on the full shower
	for(vector<FullShowerVetoPtr>::const_iterator it=_fullShowerVetoes.begin();
	it!=_fullShowerVetoes.end();++it) {
	int veto = (**it).applyVeto(currentTree());
	if(veto<0) continue;
	// veto the shower
	if(veto==0) {
	reconstructed = false;
	break;
	}
	// veto the shower and reweight
	else if(veto==1) {
	reconstructed = false;
	break;
	}
	// veto the event
	else if(veto==2) {
	throw Veto();
	}
	}
	if(reWeighting) {
	if(reconstructed) eventWeight += 1.;
	reconstructed=false;
	++nTryReWeight;
	if(nTryReWeight==_nReWeight) {
	reWeighting = false;
	if(eventWeight==0.) throw Veto();
	}
	}
	}
	while(!reconstructed&&maximumTries()>++ntry);
	// check if failed to generate the shower
	if(ntry==maximumTries()) {
	if(hard)
	throw ShowerHandler::ShowerTriesVeto(ntry);
	else
	throw Exception() << "Failed to generate the shower after "
	<< ntry << " attempts in QTildeShowerHandler::showerDecay()"
	<< Exception::eventerror;
	}
	// handle the weights and apply any reweighting required
	if(nTryReWeight>0) {
	tStdEHPtr seh = dynamic_ptr_cast<tStdEHPtr>(generator()->currentEventHandler());
	static bool first = true;
	if(seh) {
	seh->reweight(eventWeight/double(nTryReWeight));
	}
	else if(first) {
	generator()->log() << "Reweighting the shower only works with internal Herwig7 processes"
	<< "Presumably you are showering Les Houches Events. These will not be"
	<< "reweighted\n";
	first = false;
	}
	}
	// tree has now showered
	_currenttree->hasShowered(true);
	hardTree(HardTreePtr());
	}

	void QTildeShowerHandler:: convertHardTree(bool hard,ShowerInteraction type) {
	map<ColinePtr,ColinePtr> cmap;
	// incoming particles
	for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
	cit=currentTree()->incomingLines().begin();cit!=currentTree()->incomingLines().end();++cit) {
	map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
	mit = hardTree()->particles().find(cit->first->progenitor());
	// put the colour lines in the map
	ShowerParticlePtr oldParticle = cit->first->progenitor();
	ShowerParticlePtr newParticle = mit->second->branchingParticle();
	ColinePtr cLine = oldParticle-> colourLine();
	ColinePtr aLine = oldParticle->antiColourLine();
	if(newParticle->colourLine() &&
	cmap.find(newParticle-> colourLine())==cmap.end())
	cmap[newParticle-> colourLine()] = cLine;
	if(newParticle->antiColourLine() &&
	cmap.find(newParticle->antiColourLine())==cmap.end())
	cmap[newParticle->antiColourLine()] = aLine;
	// check whether or not particle emits
	bool emission = mit->second->parent();
	if(emission) {
	if(newParticle->colourLine()) {
	ColinePtr ctemp = newParticle-> colourLine();
	ctemp->removeColoured(newParticle);
	}
	if(newParticle->antiColourLine()) {
	ColinePtr ctemp = newParticle->antiColourLine();
	ctemp->removeAntiColoured(newParticle);
	}
	newParticle = mit->second->parent()->branchingParticle();
	}
	// get the new colour lines
	ColinePtr newCLine,newALine;
	// sort out colour lines
	if(newParticle->colourLine()) {
	ColinePtr ctemp = newParticle-> colourLine();
	ctemp->removeColoured(newParticle);
	if(cmap.find(ctemp)!=cmap.end()) {
	newCLine = cmap[ctemp];
	}
	else {
	newCLine = new_ptr(ColourLine());
	cmap[ctemp] = newCLine;
	}
	}
	// and anticolour lines
	if(newParticle->antiColourLine()) {
	ColinePtr ctemp = newParticle->antiColourLine();
	ctemp->removeAntiColoured(newParticle);
	if(cmap.find(ctemp)!=cmap.end()) {
	newALine = cmap[ctemp];
	}
	else {
	newALine = new_ptr(ColourLine());
	cmap[ctemp] = newALine;
	}
	}
	// remove colour lines from old particle
	if(aLine) {
	aLine->removeAntiColoured(cit->first->copy());
	aLine->removeAntiColoured(cit->first->progenitor());
	}
	if(cLine) {
	cLine->removeColoured(cit->first->copy());
	cLine->removeColoured(cit->first->progenitor());
	}
	// add particle to colour lines
	if(newCLine) newCLine->addColoured (newParticle);
	if(newALine) newALine->addAntiColoured(newParticle);
	// insert new particles
	cit->first->copy(newParticle);
	ShowerParticlePtr sp(new_ptr(ShowerParticle(*newParticle,1,false)));
	cit->first->progenitor(sp);
	currentTree()->incomingLines()[cit->first]=sp;
	cit->first->perturbative(!emission);
	// and the emitted particle if needed
	if(emission) {
	ShowerParticlePtr newOut = mit->second->parent()->children()[1]->branchingParticle();
	if(newOut->colourLine()) {
	ColinePtr ctemp = newOut-> colourLine();
	ctemp->removeColoured(newOut);
	assert(cmap.find(ctemp)!=cmap.end());
	cmap[ctemp]->addColoured (newOut);
	}
	if(newOut->antiColourLine()) {
	ColinePtr ctemp = newOut->antiColourLine();
	ctemp->removeAntiColoured(newOut);
	assert(cmap.find(ctemp)!=cmap.end());
	cmap[ctemp]->addAntiColoured(newOut);
	}
	ShowerParticlePtr sout=new_ptr(ShowerParticle(*newOut,1,true));
	ShowerProgenitorPtr out=new_ptr(ShowerProgenitor(cit->first->original(),newOut,sout));
	out->perturbative(false);
	currentTree()->outgoingLines().insert(make_pair(out,sout));
	}
	if(hard) {
	// sort out the value of x
	if(mit->second->beam()->momentum().z()>ZERO) {
	sp->x(newParticle->momentum(). plus()/mit->second->beam()->momentum(). plus());
	}
	else {
	sp->x(newParticle->momentum().minus()/mit->second->beam()->momentum().minus());
	}
	}
	}
	// outgoing particles
	for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
	cit=currentTree()->outgoingLines().begin();cit!=currentTree()->outgoingLines().end();++cit) {
	map<tShowerTreePtr,pair<tShowerProgenitorPtr,
	tShowerParticlePtr> >::const_iterator tit;
	for(tit = currentTree()->treelinks().begin();
	tit != currentTree()->treelinks().end();++tit) {
	if(tit->second.first && tit->second.second==cit->first->progenitor())
	break;
	}
	map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
	mit = hardTree()->particles().find(cit->first->progenitor());
	if(mit==hardTree()->particles().end()) continue;
	// put the colour lines in the map
	ShowerParticlePtr oldParticle = cit->first->progenitor();
	ShowerParticlePtr newParticle = mit->second->branchingParticle();
	ShowerParticlePtr newOut;
	ColinePtr cLine = oldParticle-> colourLine();
	ColinePtr aLine = oldParticle->antiColourLine();
	if(newParticle->colourLine() &&
	cmap.find(newParticle-> colourLine())==cmap.end())
	cmap[newParticle-> colourLine()] = cLine;
	if(newParticle->antiColourLine() &&
	cmap.find(newParticle->antiColourLine())==cmap.end())
	cmap[newParticle->antiColourLine()] = aLine;
	// check whether or not particle emits
	bool emission = !mit->second->children().empty();
	if(emission) {
	if(newParticle->colourLine()) {
	ColinePtr ctemp = newParticle-> colourLine();
	ctemp->removeColoured(newParticle);
	}
	if(newParticle->antiColourLine()) {
	ColinePtr ctemp = newParticle->antiColourLine();
	ctemp->removeAntiColoured(newParticle);
	}
	newParticle = mit->second->children()[0]->branchingParticle();
	newOut = mit->second->children()[1]->branchingParticle();
	if(newParticle->id()!=oldParticle->id()&&newParticle->id()==newOut->id())
	swap(newParticle,newOut);
	}
	// get the new colour lines
	ColinePtr newCLine,newALine;
	// sort out colour lines
	if(newParticle->colourLine()) {
	ColinePtr ctemp = newParticle-> colourLine();
	ctemp->removeColoured(newParticle);
	if(cmap.find(ctemp)!=cmap.end()) {
	newCLine = cmap[ctemp];
	}
	else {
	newCLine = new_ptr(ColourLine());
	cmap[ctemp] = newCLine;
	}
	}
	// and anticolour lines
	if(newParticle->antiColourLine()) {
	ColinePtr ctemp = newParticle->antiColourLine();
	ctemp->removeAntiColoured(newParticle);
	if(cmap.find(ctemp)!=cmap.end()) {
	newALine = cmap[ctemp];
	}
	else {
	newALine = new_ptr(ColourLine());
	cmap[ctemp] = newALine;
	}
	}
	// remove colour lines from old particle
	if(aLine) {
	aLine->removeAntiColoured(cit->first->copy());
	aLine->removeAntiColoured(cit->first->progenitor());
	}
	if(cLine) {
	cLine->removeColoured(cit->first->copy());
	cLine->removeColoured(cit->first->progenitor());
	}
	// special for unstable particles
	if(newParticle->id()==oldParticle->id() &&
	(tit!=currentTree()->treelinks().end()\|\|!oldParticle->dataPtr()->stable())) {
	Lorentz5Momentum oldMomentum = oldParticle->momentum();
	Lorentz5Momentum newMomentum = newParticle->momentum();
	LorentzRotation boost( oldMomentum.findBoostToCM(),oldMomentum.e()/oldMomentum.mass());
	if(tit!=currentTree()->treelinks().end()) tit->first->transform(boost,false);
	oldParticle->transform(boost);
	boost = LorentzRotation(-newMomentum.findBoostToCM(),newMomentum.e()/newMomentum.mass());
	oldParticle->transform(boost);
	if(tit!=currentTree()->treelinks().end()) tit->first->transform(boost,false);
	newParticle=oldParticle;
	}
	// add particle to colour lines
	if(newCLine) newCLine->addColoured (newParticle);
	if(newALine) newALine->addAntiColoured(newParticle);
	// insert new particles
	cit->first->copy(newParticle);
	ShowerParticlePtr sp(new_ptr(ShowerParticle(*newParticle,1,true)));
	cit->first->progenitor(sp);
	currentTree()->outgoingLines()[cit->first]=sp;
	cit->first->perturbative(!emission);
	// and the emitted particle if needed
	if(emission) {
	if(newOut->colourLine()) {
	ColinePtr ctemp = newOut-> colourLine();
	ctemp->removeColoured(newOut);
	assert(cmap.find(ctemp)!=cmap.end());
	cmap[ctemp]->addColoured (newOut);
	}
	if(newOut->antiColourLine()) {
	ColinePtr ctemp = newOut->antiColourLine();
	ctemp->removeAntiColoured(newOut);
	assert(cmap.find(ctemp)!=cmap.end());
	cmap[ctemp]->addAntiColoured(newOut);
	}
	ShowerParticlePtr sout=new_ptr(ShowerParticle(*newOut,1,true));
	ShowerProgenitorPtr out=new_ptr(ShowerProgenitor(cit->first->original(),newOut,sout));
	out->perturbative(false);
	currentTree()->outgoingLines().insert(make_pair(out,sout));
	}
	// update any decay products
	if(tit!=currentTree()->treelinks().end())
	currentTree()->updateLink(tit->first,make_pair(cit->first,sp));
	}
	// reset the tree
	currentTree()->resetShowerProducts();
	// reextract the particles and set the colour partners
	vector<ShowerParticlePtr> particles =
	currentTree()->extractProgenitorParticles();
	// clear the partners
	for(unsigned int ix=0;ix<particles.size();++ix) {
	particles[ix]->partner(ShowerParticlePtr());
	particles[ix]->clearPartners();
	}
	// clear the tree
	hardTree(HardTreePtr());
	// Set the initial evolution scales
	showerModel()->partnerFinder()->
	setInitialEvolutionScales(particles,!hard,type,!_hardtree);
	}

	Branching QTildeShowerHandler::selectTimeLikeBranching(tShowerParticlePtr particle,
	ShowerInteraction type,
	HardBranchingPtr branch) {
	Branching fb;
	unsigned int iout=0;
	while (true) {
	// break if doing truncated shower and no truncated shower needed
	if(branch && (!isTruncatedShowerON()\|\|hardOnly())) break;
	fb=_splittingGenerator->chooseForwardBranching(*particle,_finalenhance,type);
	// no emission break
	if(!fb.kinematics) break;
	// special for truncated shower
	if(branch) {
	// check haven't evolved too far
	if(fb.kinematics->scale() < branch->scale()) {
	fb=Branching();
	break;
	}
	// find the truncated line
	iout=0;
	if(fb.ids[1]->id()!=fb.ids[2]->id()) {
	if(fb.ids[1]->id()==particle->id()) iout=1;
	else if (fb.ids[2]->id()==particle->id()) iout=2;
	}
	else if(fb.ids[1]->id()==particle->id()) {
	if(fb.kinematics->z()>0.5) iout=1;
	else iout=2;
	}
	// apply the vetos for the truncated shower
	// no flavour changing branchings
	if(iout==0) {
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	continue;
	}
	double zsplit = iout==1 ? fb.kinematics->z() : 1-fb.kinematics->z();
	// only if same interaction for forced branching
	ShowerInteraction type2 = convertInteraction(fb.type);
	// and evolution
	if(type2==branch->sudakov()->interactionType()) {
	if(zsplit < 0.5 \|\| // hardest line veto
	fb.kinematics->scale()*zsplit < branch->scale() ) { // angular ordering veto
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	continue;
	}
	}
	// pt veto
	if(fb.kinematics->pT() > progenitor()->maximumpT(type2)) {
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	continue;
	}
	}
	// standard vetos for all emissions
	if(timeLikeVetoed(fb,particle)) {
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	if(particle->spinInfo()) particle->spinInfo()->decayVertex(VertexPtr());
	continue;
	}
	// special for already decayed particles
	// don't allow flavour changing branchings
	bool vetoDecay = false;
	for(map<tShowerTreePtr,pair<tShowerProgenitorPtr,
	tShowerParticlePtr> >::const_iterator tit = currentTree()->treelinks().begin();
	tit != currentTree()->treelinks().end();++tit) {
	if(tit->second.first == progenitor()) {
	map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
	it = currentTree()->outgoingLines().find(progenitor());
	if(it!=currentTree()->outgoingLines().end() && particle == it->second &&
	fb.ids[0]!=fb.ids[1] && fb.ids[1]!=fb.ids[2]) {
	vetoDecay = true;
	break;
	}
	}
	}
	if(vetoDecay) {
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	if(particle->spinInfo()) particle->spinInfo()->decayVertex(VertexPtr());
	continue;
	}
	break;
	}
	// normal case
	if(!branch) {
	if(fb.kinematics) fb.hard = false;
	return fb;
	}
	// truncated emission
	if(fb.kinematics) {
	fb.hard = false;
	fb.iout = iout;
	return fb;
	}
	// otherwise need to return the hard emission
	// construct the kinematics for the hard emission
	ShoKinPtr showerKin=
	branch->sudakov()->createFinalStateBranching(branch->scale(),
	branch->children()[0]->z(),
	branch->phi(),
	branch->children()[0]->pT());
	IdList idlist(3);
	idlist[0] = particle->dataPtr();
	idlist[1] = branch->children()[0]->branchingParticle()->dataPtr();
	idlist[2] = branch->children()[1]->branchingParticle()->dataPtr();
	fb = Branching( showerKin, idlist, branch->sudakov(),branch->type() );
	fb.hard = true;
	fb.iout=0;
	// return it
	return fb;
	}

	Branching QTildeShowerHandler::selectSpaceLikeDecayBranching(tShowerParticlePtr particle,
	const ShowerParticle::EvolutionScales & maxScales,
	Energy minmass,ShowerInteraction type,
	HardBranchingPtr branch) {
	Branching fb;
	unsigned int iout=0;
	while (true) {
	// break if doing truncated shower and no truncated shower needed
	if(branch && (!isTruncatedShowerON()\|\|hardOnly())) break;
	// select branching
	fb=_splittingGenerator->chooseDecayBranching(*particle,maxScales,minmass,
	_initialenhance,type);
	// return if no radiation
	if(!fb.kinematics) break;
	// special for truncated shower
	if(branch) {
	// check haven't evolved too far
	if(fb.kinematics->scale() < branch->scale()) {
	fb=Branching();
	break;
	}
	// find the truncated line
	iout=0;
	if(fb.ids[1]->id()!=fb.ids[2]->id()) {
	if(fb.ids[1]->id()==particle->id()) iout=1;
	else if (fb.ids[2]->id()==particle->id()) iout=2;
	}
	else if(fb.ids[1]->id()==particle->id()) {
	if(fb.kinematics->z()>0.5) iout=1;
	else iout=2;
	}
	// apply the vetos for the truncated shower
	// no flavour changing branchings
	if(iout==0) {
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	continue;
	}
	ShowerInteraction type2 = convertInteraction(fb.type);
	double zsplit = iout==1 ? fb.kinematics->z() : 1-fb.kinematics->z();
	if(type2==branch->sudakov()->interactionType()) {
	if(zsplit < 0.5 \|\| // hardest line veto
	fb.kinematics->scale()*zsplit < branch->scale() ) { // angular ordering veto
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	continue;
	}
	}
	// pt veto
	if(fb.kinematics->pT() > progenitor()->maximumpT(type2)) {
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	continue;
	}
	}
	// if not vetoed break
	if(spaceLikeDecayVetoed(fb,particle)) {
	// otherwise reset scale and continue
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	continue;
	}
	break;
	}
	// normal case
	if(!branch) {
	if(fb.kinematics) fb.hard = false;
	return fb;
	}
	// truncated emission
	if(fb.kinematics) {
	fb.hard = false;
	fb.iout = iout;
	return fb;
	}
	// otherwise need to return the hard emission
	// construct the kinematics for the hard emission
	ShoKinPtr showerKin=
	branch->sudakov()->createDecayBranching(branch->scale(),
	branch->children()[0]->z(),
	branch->phi(),
	branch->children()[0]->pT());
	IdList idlist(3);
	idlist[0] = particle->dataPtr();
	idlist[1] = branch->children()[0]->branchingParticle()->dataPtr();
	idlist[2] = branch->children()[1]->branchingParticle()->dataPtr();
	// create the branching
	fb = Branching( showerKin, idlist, branch->sudakov(),ShowerPartnerType::QCDColourLine );
	fb.hard=true;
	fb.iout=0;
	// return it
	return fb;
	}

	void QTildeShowerHandler::checkFlags() {
	string error = "Inconsistent hard emission set-up in QTildeShowerHandler::showerHardProcess(). ";
	if ( ( currentTree()->isMCatNLOSEvent() \|\| currentTree()->isMCatNLOHEvent() ) ) {
	if (_hardEmission ==2 )
	throw Exception() << error
	<< "Cannot generate POWHEG matching with MC@NLO shower "
	<< "approximation. Add 'set QTildeShowerHandler:HardEmission 0' to input file."
	<< Exception::runerror;
	if ( canHandleMatchboxTrunc() )
	throw Exception() << error
	<< "Cannot use truncated qtilde shower with MC@NLO shower "
	<< "approximation. Set LHCGenerator:EventHandler"
	<< ":CascadeHandler to '/Herwig/Shower/ShowerHandler' or "
	<< "'/Herwig/Shower/Dipole/DipoleShowerHandler'."
	<< Exception::runerror;
	}
	else if ( ((currentTree()->isPowhegSEvent() \|\| currentTree()->isPowhegHEvent()) ) &&
	_hardEmission != 2){
	if ( canHandleMatchboxTrunc())
	throw Exception() << error
	<< "Unmatched events requested for POWHEG shower "
	<< "approximation. Set QTildeShowerHandler:HardEmission to "
	<< "'POWHEG'."
	<< Exception::runerror;
	else if (_hardEmissionWarn) {
	_hardEmissionWarn = false;
	_hardEmission=2;
	throw Exception() << error
	<< "Unmatched events requested for POWHEG shower "
	<< "approximation. Changing QTildeShowerHandler:HardEmission from "
	<< _hardEmission << " to 2"
	<< Exception::warning;
	}
	}

	if ( currentTree()->isPowhegSEvent() \|\| currentTree()->isPowhegHEvent()) {
	if (currentTree()->showerApproximation()->needsTruncatedShower() &&
	!canHandleMatchboxTrunc() )
	throw Exception() << error
	<< "Current shower handler cannot generate truncated shower. "
	<< "Set Generator:EventHandler:CascadeHandler to "
	<< "'/Herwig/Shower/PowhegShowerHandler'."
	<< Exception::runerror;
	}
	else if ( currentTree()->truncatedShower() && _missingTruncWarn) {
	_missingTruncWarn=false;
	throw Exception() << "Warning: POWHEG shower approximation used without "
	<< "truncated shower. Set Generator:EventHandler:"
	<< "CascadeHandler to '/Herwig/Shower/PowhegShowerHandler' and "
	<< "'MEMatching:TruncatedShower Yes'."
	<< Exception::warning;
	}
	// else if ( !dipme && _hardEmissionMode > 1 &&
	// firstInteraction())
	// throw Exception() << error
	// << "POWHEG matching requested for LO events. Include "
	// << "'set Factory:ShowerApproximation MEMatching' in input file."
	// << Exception::runerror;
	}


	tPPair QTildeShowerHandler::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 QTildeShowerHandler::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 == incomingBeams().first \|\|
	parent == incomingBeams().second ) return seed;
	else return findFirstParton(parent);
	}

	void QTildeShowerHandler::decay(ShowerTreePtr tree, ShowerDecayMap & decay) {
	// must be one incoming particle
	assert(tree->incomingLines().size()==1);
	// apply any transforms
	tree->applyTransforms();
	// if already decayed return
	if(!tree->outgoingLines().empty()) return;
	// now we need to replace the particle with a new copy after the shower
	// find particle after the shower
	map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator
	tit = tree->parent()->treelinks().find(tree);
	assert(tit!=tree->parent()->treelinks().end());
	ShowerParticlePtr newparent=tit->second.second;
	PerturbativeProcessPtr newProcess = new_ptr(PerturbativeProcess());
	newProcess->incoming().push_back(make_pair(newparent,PerturbativeProcessPtr()));
	DecayProcessMap decayMap;
	ShowerHandler::decay(newProcess,decayMap);
	ShowerTree::constructTrees(tree,decay,newProcess,decayMap);
	}

	namespace {

	ShowerProgenitorPtr
	findFinalStateLine(ShowerTreePtr tree, long id, Lorentz5Momentum momentum) {
	map<ShowerProgenitorPtr,tShowerParticlePtr>::iterator partner;
	Energy2 dmin(1e30*GeV2);
	for(map<ShowerProgenitorPtr,tShowerParticlePtr>::iterator
	cit =tree->outgoingLines().begin(); cit!=tree->outgoingLines().end(); ++cit) {
	if(cit->second->id()!=id) continue;
	Energy2 test =
	sqr(cit->second->momentum().x()-momentum.x())+
	sqr(cit->second->momentum().y()-momentum.y())+
	sqr(cit->second->momentum().z()-momentum.z())+
	sqr(cit->second->momentum().t()-momentum.t());
	if(test<dmin) {
	dmin = test;
	partner = cit;
	}
	}
	return partner->first;
	}

	ShowerProgenitorPtr
	findInitialStateLine(ShowerTreePtr tree, long id, Lorentz5Momentum momentum) {
	map<ShowerProgenitorPtr,ShowerParticlePtr>::iterator partner;
	Energy2 dmin(1e30*GeV2);
	for(map<ShowerProgenitorPtr,ShowerParticlePtr>::iterator
	cit =tree->incomingLines().begin(); cit!=tree->incomingLines().end(); ++cit) {
	if(cit->second->id()!=id) continue;
	Energy2 test =
	sqr(cit->second->momentum().x()-momentum.x())+
	sqr(cit->second->momentum().y()-momentum.y())+
	sqr(cit->second->momentum().z()-momentum.z())+
	sqr(cit->second->momentum().t()-momentum.t());
	if(test<dmin) {
	dmin = test;
	partner = cit;
	}
	}
	return partner->first;
	}

	void fixSpectatorColours(PPtr newSpect,ShowerProgenitorPtr oldSpect,
	ColinePair & cline,ColinePair & aline, bool reconnect) {
	cline.first = oldSpect->progenitor()->colourLine();
	cline.second = newSpect->colourLine();
	aline.first = oldSpect->progenitor()->antiColourLine();
	aline.second = newSpect->antiColourLine();
	if(!reconnect) return;
	if(cline.first) {
	cline.first ->removeColoured(oldSpect->copy());
	cline.first ->removeColoured(oldSpect->progenitor());
	cline.second->removeColoured(newSpect);
	cline.first ->addColoured(newSpect);
	}
	if(aline.first) {
	aline.first ->removeAntiColoured(oldSpect->copy());
	aline.first ->removeAntiColoured(oldSpect->progenitor());
	aline.second->removeAntiColoured(newSpect);
	aline.first ->addAntiColoured(newSpect);
	}
	}

	void fixInitialStateEmitter(ShowerTreePtr tree, PPtr newEmit,PPtr emitted, ShowerProgenitorPtr emitter,
	ColinePair cline,ColinePair aline,double x) {
	// sort out the colours
	if(emitted->dataPtr()->iColour()==PDT::Colour8) {
	// emitter
	if(cline.first && cline.first == emitter->progenitor()->antiColourLine() &&
	cline.second !=newEmit->antiColourLine()) {
	// sort out not radiating line
	ColinePtr col = emitter->progenitor()->colourLine();
	if(col) {
	col->removeColoured(emitter->copy());
	col->removeColoured(emitter->progenitor());
	newEmit->colourLine()->removeColoured(newEmit);
	col->addColoured(newEmit);
	}
	}
	else if(aline.first && aline.first == emitter->progenitor()->colourLine() &&
	aline.second !=newEmit->colourLine()) {
	// sort out not radiating line
	ColinePtr anti = emitter->progenitor()->antiColourLine();
	if(anti) {
	anti->removeAntiColoured(emitter->copy());
	anti->removeAntiColoured(emitter->progenitor());
	newEmit->colourLine()->removeAntiColoured(newEmit);
	anti->addAntiColoured(newEmit);
	}
	}
	else
	assert(false);
	// emitted
	if(cline.first && cline.second==emitted->colourLine()) {
	cline.second->removeColoured(emitted);
	cline.first->addColoured(emitted);
	}
	else if(aline.first && aline.second==emitted->antiColourLine()) {
	aline.second->removeAntiColoured(emitted);
	aline.first->addAntiColoured(emitted);
	}
	else
	assert(false);
	}
	else {
	if(emitter->progenitor()->antiColourLine() ) {
	ColinePtr col = emitter->progenitor()->antiColourLine();
	col->removeAntiColoured(emitter->copy());
	col->removeAntiColoured(emitter->progenitor());
	if(newEmit->antiColourLine()) {
	newEmit->antiColourLine()->removeAntiColoured(newEmit);
	col->addAntiColoured(newEmit);
	}
	else if (emitted->colourLine()) {
	emitted->colourLine()->removeColoured(emitted);
	col->addColoured(emitted);
	}
	else
	assert(false);
	}
	if(emitter->progenitor()->colourLine() ) {
	ColinePtr col = emitter->progenitor()->colourLine();
	col->removeColoured(emitter->copy());
	col->removeColoured(emitter->progenitor());
	if(newEmit->colourLine()) {
	newEmit->colourLine()->removeColoured(newEmit);
	col->addColoured(newEmit);
	}
	else if (emitted->antiColourLine()) {
	emitted->antiColourLine()->removeAntiColoured(emitted);
	col->addAntiColoured(emitted);
	}
	else
	assert(false);
	}
	}
	// update the emitter
	emitter->copy(newEmit);
	ShowerParticlePtr sp = new_ptr(ShowerParticle(*newEmit,1,false));
	sp->x(x);
	emitter->progenitor(sp);
	tree->incomingLines()[emitter]=sp;
	emitter->perturbative(false);
	// add emitted
	sp=new_ptr(ShowerParticle(*emitted,1,true));
	ShowerProgenitorPtr gluon=new_ptr(ShowerProgenitor(emitter->original(),emitted,sp));
	gluon->perturbative(false);
	tree->outgoingLines().insert(make_pair(gluon,sp));
	}

	void fixFinalStateEmitter(ShowerTreePtr tree, PPtr newEmit,PPtr emitted, ShowerProgenitorPtr emitter,
	ColinePair cline,ColinePair aline) {
	map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator tit;
	// special case if decayed
	for(tit = tree->treelinks().begin(); tit != tree->treelinks().end();++tit) {
	if(tit->second.first && tit->second.second==emitter->progenitor())
	break;
	}
	// sort out the colour lines
	if(cline.first && cline.first == emitter->progenitor()->antiColourLine() &&
	cline.second !=newEmit->antiColourLine()) {
	// sort out not radiating line
	ColinePtr col = emitter->progenitor()->colourLine();
	if(col) {
	col->removeColoured(emitter->copy());
	col->removeColoured(emitter->progenitor());
	newEmit->colourLine()->removeColoured(newEmit);
	col->addColoured(newEmit);
	}
	}
	else if(aline.first && aline.first == emitter->progenitor()->colourLine() &&
	aline.second !=newEmit->colourLine()) {
	// sort out not radiating line
	ColinePtr anti = emitter->progenitor()->antiColourLine();
	if(anti) {
	anti->removeAntiColoured(emitter->copy());
	anti->removeAntiColoured(emitter->progenitor());
	newEmit->colourLine()->removeAntiColoured(newEmit);
	anti->addAntiColoured(newEmit);
	}
	}
	else
	assert(false);
	// update the emitter
	emitter->copy(newEmit);
	ShowerParticlePtr sp = new_ptr(ShowerParticle(*newEmit,1,true));
	emitter->progenitor(sp);
	tree->outgoingLines()[emitter]=sp;
	emitter->perturbative(false);
	// update for decaying particles
	if(tit!=tree->treelinks().end())
	tree->updateLink(tit->first,make_pair(emitter,sp));
	// add the emitted particle
	// sort out the colour
	if(cline.first && cline.second==emitted->antiColourLine()) {
	cline.second->removeAntiColoured(emitted);
	cline.first->addAntiColoured(emitted);
	}
	else if(aline.first && aline.second==emitted->colourLine()) {
	aline.second->removeColoured(emitted);
	aline.first->addColoured(emitted);
	}
	else
	assert(false);
	sp=new_ptr(ShowerParticle(*emitted,1,true));
	ShowerProgenitorPtr gluon=new_ptr(ShowerProgenitor(emitter->original(),
	emitted,sp));
	gluon->perturbative(false);
	tree->outgoingLines().insert(make_pair(gluon,sp));
	}

	}

	void QTildeShowerHandler::setupMECorrection(RealEmissionProcessPtr real) {
	assert(real);
	currentTree()->hardMatrixElementCorrection(true);
	// II emission
	if(real->emitter() < real->incoming().size() &&
	real->spectator() < real->incoming().size()) {
	// recoiling system
	for( map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
	cjt= currentTree()->outgoingLines().begin();
	cjt != currentTree()->outgoingLines().end();++cjt ) {
	cjt->first->progenitor()->transform(real->transformation());
	cjt->first->copy()->transform(real->transformation());
	}
	// the the radiating system
	ShowerProgenitorPtr emitter,spectator;
	unsigned int iemit = real->emitter();
	unsigned int ispect = real->spectator();
	int ig = int(real->emitted())-int(real->incoming().size());
	emitter = findInitialStateLine(currentTree(),
	real->bornIncoming()[iemit]->id(),
	real->bornIncoming()[iemit]->momentum());
	spectator = findInitialStateLine(currentTree(),
	real->bornIncoming()[ispect]->id(),
	real->bornIncoming()[ispect]->momentum());
	// sort out the colours
	ColinePair cline,aline;
	fixSpectatorColours(real->incoming()[ispect],spectator,cline,aline,true);
	// update the spectator
	spectator->copy(real->incoming()[ispect]);
	ShowerParticlePtr sp(new_ptr(ShowerParticle(*real->incoming()[ispect],1,false)));
	sp->x(ispect ==0 ? real->x().first :real->x().second);
	spectator->progenitor(sp);
	currentTree()->incomingLines()[spectator]=sp;
	spectator->perturbative(true);
	// now for the emitter
	fixInitialStateEmitter(currentTree(),real->incoming()[iemit],real->outgoing()[ig],
	emitter,cline,aline,iemit ==0 ? real->x().first :real->x().second);
	}
	// FF emission
	else if(real->emitter() >= real->incoming().size() &&
	real->spectator() >= real->incoming().size()) {
	assert(real->outgoing()[real->emitted()-real->incoming().size()]->id()==ParticleID::g);
	// find the emitter and spectator in the shower tree
	ShowerProgenitorPtr emitter,spectator;
	int iemit = int(real->emitter())-int(real->incoming().size());
	emitter = findFinalStateLine(currentTree(),
	real->bornOutgoing()[iemit]->id(),
	real->bornOutgoing()[iemit]->momentum());
	int ispect = int(real->spectator())-int(real->incoming().size());
	spectator = findFinalStateLine(currentTree(),
	real->bornOutgoing()[ispect]->id(),
	real->bornOutgoing()[ispect]->momentum());
	map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator tit;
	// first the spectator
	// special case if decayed
	for(tit = currentTree()->treelinks().begin(); tit != currentTree()->treelinks().end();++tit) {
	if(tit->second.first && tit->second.second==spectator->progenitor())
	break;
	}
	// sort out the colours
	ColinePair cline,aline;
	fixSpectatorColours(real->outgoing()[ispect],spectator,cline,aline,true);
	// update the spectator
	spectator->copy(real->outgoing()[ispect]);
	ShowerParticlePtr sp(new_ptr(ShowerParticle(*real->outgoing()[ispect],1,true)));
	spectator->progenitor(sp);
	currentTree()->outgoingLines()[spectator]=sp;
	spectator->perturbative(true);
	// update for decaying particles
	if(tit!=currentTree()->treelinks().end())
	currentTree()->updateLink(tit->first,make_pair(spectator,sp));
	// now the emitting particle
	int ig = int(real->emitted())-int(real->incoming().size());
	fixFinalStateEmitter(currentTree(),real->outgoing()[iemit],
	real->outgoing()[ig],
	emitter,cline,aline);
	}
	// IF emission
	else {
	// scattering process
	if(real->incoming().size()==2) {
	ShowerProgenitorPtr emitter,spectator;
	unsigned int iemit = real->emitter();
	unsigned int ispect = real->spectator();
	int ig = int(real->emitted())-int(real->incoming().size());
	ColinePair cline,aline;
	// incoming spectator
	if(ispect<2) {
	spectator = findInitialStateLine(currentTree(),
	real->bornIncoming()[ispect]->id(),
	real->bornIncoming()[ispect]->momentum());
	fixSpectatorColours(real->incoming()[ispect],spectator,cline,aline,true);
	// update the spectator
	spectator->copy(real->incoming()[ispect]);
	ShowerParticlePtr sp(new_ptr(ShowerParticle(*real->incoming()[ispect],1,false)));
	sp->x(ispect ==0 ? real->x().first :real->x().second);
	spectator->progenitor(sp);
	currentTree()->incomingLines()[spectator]=sp;
	spectator->perturbative(true);
	}
	// outgoing spectator
	else {
	spectator = findFinalStateLine(currentTree(),
	real->bornOutgoing()[ispect-real->incoming().size()]->id(),
	real->bornOutgoing()[ispect-real->incoming().size()]->momentum());
	// special case if decayed
	map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator tit;
	for(tit = currentTree()->treelinks().begin(); tit != currentTree()->treelinks().end();++tit) {
	if(tit->second.first && tit->second.second==spectator->progenitor())
	break;
	}
	fixSpectatorColours(real->outgoing()[ispect-real->incoming().size()],spectator,cline,aline,true);
	// update the spectator
	spectator->copy(real->outgoing()[ispect-real->incoming().size()]);
	ShowerParticlePtr sp(new_ptr(ShowerParticle(*real->outgoing()[ispect-real->incoming().size()],1,true)));
	spectator->progenitor(sp);
	currentTree()->outgoingLines()[spectator]=sp;
	spectator->perturbative(true);
	// update for decaying particles
	if(tit!=currentTree()->treelinks().end())
	currentTree()->updateLink(tit->first,make_pair(spectator,sp));
	}
	// incoming emitter
	if(iemit<2) {
	emitter = findInitialStateLine(currentTree(),
	real->bornIncoming()[iemit]->id(),
	real->bornIncoming()[iemit]->momentum());
	fixInitialStateEmitter(currentTree(),real->incoming()[iemit],real->outgoing()[ig],
	emitter,aline,cline,iemit ==0 ? real->x().first :real->x().second);
	}
	// outgoing emitter
	else {
	emitter = findFinalStateLine(currentTree(),
	real->bornOutgoing()[iemit-real->incoming().size()]->id(),
	real->bornOutgoing()[iemit-real->incoming().size()]->momentum());
	fixFinalStateEmitter(currentTree(),real->outgoing()[iemit-real->incoming().size()],
	real->outgoing()[ig],emitter,aline,cline);
	}
	}
	// decay process
	else {
	assert(real->spectator()==0);
	unsigned int iemit = real->emitter()-real->incoming().size();
	int ig = int(real->emitted())-int(real->incoming().size());
	ColinePair cline,aline;
	// incoming spectator
	ShowerProgenitorPtr spectator = findInitialStateLine(currentTree(),
	real->bornIncoming()[0]->id(),
	real->bornIncoming()[0]->momentum());
	fixSpectatorColours(real->incoming()[0],spectator,cline,aline,false);
	// find the emitter
	ShowerProgenitorPtr emitter =
	findFinalStateLine(currentTree(),
	real->bornOutgoing()[iemit]->id(),
	real->bornOutgoing()[iemit]->momentum());
	// recoiling system
	for( map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
	cjt= currentTree()->outgoingLines().begin();
	cjt != currentTree()->outgoingLines().end();++cjt ) {
	if(cjt->first==emitter) continue;
	cjt->first->progenitor()->transform(real->transformation());
	cjt->first->copy()->transform(real->transformation());
	}
	// sort out the emitter
	fixFinalStateEmitter(currentTree(),real->outgoing()[iemit],
	real->outgoing()[ig],emitter,aline,cline);
	}
	}
	// clean up the shower tree
	_currenttree->resetShowerProducts();
	}
	diff --git a/Shower/QTilde/QTildeShowerHandler.h b/Shower/QTilde/QTildeShowerHandler.h
	--- a/Shower/QTilde/QTildeShowerHandler.h
	+++ b/Shower/QTilde/QTildeShowerHandler.h
	@@ -1,859 +1,852 @@
	// -- C++ --
	#ifndef Herwig_QTildeShowerHandler_H
	#define Herwig_QTildeShowerHandler_H
	//
	// This is the declaration of the QTildeShowerHandler class.
	//

	#include "QTildeShowerHandler.fh"
	#include "Herwig/Shower/ShowerHandler.h"
	#include "Herwig/Shower/QTilde/Base/ShowerModel.h"
	#include "Herwig/Shower/QTilde/SplittingFunctions/SplittingGenerator.h"
	#include "Herwig/Shower/QTilde/Base/ShowerTree.h"
	#include "Herwig/Shower/QTilde/Base/ShowerProgenitor.fh"
	#include "Herwig/Shower/QTilde/Base/HardTree.h"
	#include "Herwig/Shower/QTilde/Base/Branching.h"
	#include "Herwig/Shower/QTilde/Base/ShowerVeto.h"
	#include "Herwig/Shower/QTilde/Base/FullShowerVeto.h"
	#include "Herwig/MatrixElement/HwMEBase.h"
	#include "Herwig/Decay/HwDecayerBase.h"
	#include "Herwig/MatrixElement/Matchbox/Matching/ShowerApproximation.h"
	#include "Herwig/Shower/RealEmissionProcess.h"
	#include "Herwig/Utilities/Statistic.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* The QTildeShowerHandler class.
	*
	* @see \ref QTildeShowerHandlerInterfaces "The interfaces"
	* defined for QTildeShowerHandler.
	*/
	class QTildeShowerHandler: public ShowerHandler {

	public:

	/**
	* Pointer to an XComb object
	*/
	typedef Ptr<XComb>::pointer XCPtr;

	public:

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

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

	public:

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

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

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

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

	/**
	* Members to perform the shower
	*/
	//@{
	/**
	* Perform the shower of the hard process
	*/
	virtual void showerHardProcess(ShowerTreePtr,XCPtr);

	/**
	* Perform the shower of a decay
	*/
	virtual void showerDecay(ShowerTreePtr);
	//@}

	/**
	* Access to the flags and shower variables
	*/
	//@{
	/**
	* Get the ShowerModel
	*/
	ShowerModelPtr showerModel() const {return _model;}

	/**
	* Get the SplittingGenerator
	*/
	tSplittingGeneratorPtr splittingGenerator() const { return _splittingGenerator; }

	/**
	* Mode for hard emissions
	*/
	int hardEmission() const {return _hardEmission;}
	//@}

	/**
	* Connect the Hard and Shower trees
	*/
	virtual void connectTrees(ShowerTreePtr showerTree, HardTreePtr hardTree, bool hard );

	/**
	* Access to switches for spin correlations
	*/
	//@{
	/**
	* Spin Correlations
	*/
	unsigned int spinCorrelations() const {
	return _spinOpt;
	}

	/**
	* Soft correlations
	*/
	unsigned int softCorrelations() const {
	return _softOpt;
	}

	/**
	* Any correlations
	*/
	bool correlations() const {
	return _spinOpt!=0\|\|_softOpt!=0;
	}
	//@}

	protected:

	/**
	* Perform the shower
	*/
	void doShowering(bool hard,XCPtr);

	/**
	* Generate the hard matrix element correction
	*/
	virtual RealEmissionProcessPtr hardMatrixElementCorrection(bool);

	/**
	* Generate the hardest emission
	*/
	virtual void hardestEmission(bool hard);

	/**
	* Set up for applying a matrix element correction
	*/
	void setupMECorrection(RealEmissionProcessPtr real);

	/**
	* Extract the particles to be showered, set the evolution scales
	* and apply the hard matrix element correction
	* @param hard Whether this is a hard process or decay
	* @return The particles to be showered
	*/
	virtual vector<ShowerProgenitorPtr> setupShower(bool hard);

	/**
	* set the colour partners
	*/
	virtual void setEvolutionPartners(bool hard,ShowerInteraction,
	bool clear);

	/**
	* Methods to perform the evolution of an individual particle, including
	* recursive calling on the products
	*/
	//@{
	/**
	* It does the forward evolution of the time-like input particle
	* (and recursively for all its radiation products).
	* accepting only emissions which conforms to the showerVariables
	* and soft matrix element correction.
	* If at least one emission has occurred then the method returns true.
	* @param particle The particle to be showered
	*/
	virtual bool timeLikeShower(tShowerParticlePtr particle, ShowerInteraction,
	Branching fb, bool first);

	/**
	* It does the backward evolution of the space-like input particle
	* (and recursively for all its time-like radiation products).
	* accepting only emissions which conforms to the showerVariables.
	* If at least one emission has occurred then the method returns true
	* @param particle The particle to be showered
	* @param beam The beam particle
	*/
	virtual bool spaceLikeShower(tShowerParticlePtr particle,PPtr beam,
	ShowerInteraction);

	/**
	* If does the forward evolution of the input on-shell particle
	* involved in a decay
	* (and recursively for all its time-like radiation products).
	* accepting only emissions which conforms to the showerVariables.
	* @param particle The particle to be showered
	* @param maxscale The maximum scale for the shower.
	* @param minimumMass The minimum mass of the final-state system
	*/
	virtual bool
	spaceLikeDecayShower(tShowerParticlePtr particle,
	const ShowerParticle::EvolutionScales & maxScales,
	Energy minimumMass,ShowerInteraction,
	Branching fb);

	/**
	* Truncated shower from a time-like particle
	*/
	virtual bool truncatedTimeLikeShower(tShowerParticlePtr particle,
	HardBranchingPtr branch,
	ShowerInteraction type,
	Branching fb, bool first);

	/**
	* Truncated shower from a space-like particle
	*/
	virtual bool truncatedSpaceLikeShower(tShowerParticlePtr particle,PPtr beam,
	HardBranchingPtr branch,
	ShowerInteraction type);

	/**
	* Truncated shower from a time-like particle
	*/
	virtual bool truncatedSpaceLikeDecayShower(tShowerParticlePtr particle,
	const ShowerParticle::EvolutionScales & maxScales,
	Energy minimumMass, HardBranchingPtr branch,
	ShowerInteraction type, Branching fb);
	//@}

	/**
	* Switches for matrix element corrections
	*/
	//@{
	/**
	* Any ME correction?
	*/
	bool MECOn() const {
	return _hardEmission == 1;
	}

	/**
	* Any hard ME correction?
	*/
	bool hardMEC() const {
	return _hardEmission == 1 && (_meCorrMode == 1 \|\| _meCorrMode == 2);
	}

	/**
	* Any soft ME correction?
	*/
	bool softMEC() const {
	return _hardEmission == 1 && (_meCorrMode == 1 \|\| _meCorrMode > 2);
	}
	//@}

	/**
	* Is the truncated shower on?
	*/
	bool isTruncatedShowerON() const {return _trunc_Mode;}

	/**
	* Switch for intrinsic pT
	*/
	//@{
	/**
	* Any intrinsic pT?
	*/
	bool ipTon() const {
	return _iptrms != ZERO \|\| ( _beta == 1.0 && _gamma != ZERO && _iptmax !=ZERO );
	}
	//@}

	/*@name Additional shower vetoes /
	//@{
	/**
	* Insert a veto.
	*/
	void addVeto (ShowerVetoPtr v) { _vetoes.push_back(v); }

	/**
	* Remove a veto.
	*/
	void removeVeto (ShowerVetoPtr v) {
	vector<ShowerVetoPtr>::iterator vit = find(_vetoes.begin(),_vetoes.end(),v);
	if (vit != _vetoes.end())
	_vetoes.erase(vit);
	}

	//@}

	/**
	* Switches for vetoing hard emissions
	*/
	//@{
	/**
	* Returns true if the hard veto read-in is to be applied to only
	* the primary collision and false otherwise.
	*/
	bool hardVetoReadOption() const {return _hardVetoReadOption;}
	//@}

	/**
	* Enhancement factors for radiation needed to generate the soft matrix
	* element correction.
	*/
	//@{
	/**
	* Access the enhancement factor for initial-state radiation
	*/
	double initialStateRadiationEnhancementFactor() const { return _initialenhance; }

	/**
	* Access the enhancement factor for final-state radiation
	*/
	double finalStateRadiationEnhancementFactor() const { return _finalenhance; }

	/**
	* Set the enhancement factor for initial-state radiation
	*/
	void initialStateRadiationEnhancementFactor(double in) { _initialenhance=in; }

	/**
	* Set the enhancement factor for final-state radiation
	*/
	void finalStateRadiationEnhancementFactor(double in) { _finalenhance=in; }
	//@}

	/**
	* Access to set/get the HardTree currently beinging showered
	*/
	//@{
	/**
	* The HardTree currently being showered
	*/
	tHardTreePtr hardTree() {return _hardtree;}

	/**
	* The HardTree currently being showered
	*/
	void hardTree(tHardTreePtr in) {_hardtree = in;}
	//@}

	/**
	* Access/set the beam particle for the current initial-state shower
	*/
	//@{
	/**
	* Get the beam particle data
	*/
	Ptr<BeamParticleData>::const_pointer beamParticle() const { return _beam; }

	/**
	* Set the beam particle data
	*/
	void setBeamParticle(Ptr<BeamParticleData>::const_pointer in) { _beam=in; }
	//@}

	/**
	* Set/Get the current tree being evolver for inheriting classes
	*/
	//@{
	/**
	* Get the tree
	*/
	tShowerTreePtr currentTree() { return _currenttree; }

	/**
	* Set the tree
	*/
	void currentTree(tShowerTreePtr tree) { _currenttree=tree; }

	//@}

	/**
	* Access the maximum number of attempts to generate the shower
	*/
	unsigned int maximumTries() const { return _maxtry; }

	/**
	* Set/Get the ShowerProgenitor for the current shower
	*/
	//@{
	/**
	* Access the progenitor
	*/
	ShowerProgenitorPtr progenitor() { return _progenitor; }

	/**
	* Set the progenitor
	*/
	void progenitor(ShowerProgenitorPtr in) { _progenitor=in; }
	//@}

	/**
	* Calculate the intrinsic \f$p_T\f$.
	*/
	virtual void generateIntrinsicpT(vector<ShowerProgenitorPtr>);

	/**
	* Access to the intrinsic \f$p_T\f$ for inheriting classes
	*/
	map<tShowerProgenitorPtr,pair<Energy,double> > & intrinsicpT() { return _intrinsic; }

	/**
	* find the maximally allowed pt acc to the hard process.
	*/
	void setupMaximumScales(const vector<ShowerProgenitorPtr> &,XCPtr);

	/**
	* find the relevant hard scales for profile scales.
	*/
	void setupHardScales(const vector<ShowerProgenitorPtr> &,XCPtr);

	/**
	* Convert the HardTree into an extra shower emission
	*/
	void convertHardTree(bool hard,ShowerInteraction type);

	protected:

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

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

	protected:

	/**
	* Start the shower of a timelike particle
	*/
	virtual bool startTimeLikeShower(ShowerInteraction);

	/**
	* Update of the time-like stuff
	*/
	void updateHistory(tShowerParticlePtr particle);

	/**
	* Start the shower of a spacelike particle
	*/
	virtual bool startSpaceLikeShower(PPtr,ShowerInteraction);

	/**
	* Start the shower of a spacelike particle
	*/
	virtual bool
	startSpaceLikeDecayShower(const ShowerParticle::EvolutionScales & maxScales,
	Energy minimumMass,ShowerInteraction);

	/**
	* Select the branching for the next time-like emission
	*/
	Branching selectTimeLikeBranching(tShowerParticlePtr particle,
	ShowerInteraction type,
	HardBranchingPtr branch);

	/**
	* Select the branching for the next space-like emission in a decay
	*/
	Branching selectSpaceLikeDecayBranching(tShowerParticlePtr particle,
	const ShowerParticle::EvolutionScales & maxScales,
	Energy minmass,ShowerInteraction type,
	HardBranchingPtr branch);
	/**
	* Create the timelike child of a branching
	*/
	ShowerParticleVector createTimeLikeChildren(tShowerParticlePtr particle,
	IdList ids);

	/**
	* Vetos for the timelike shower
	*/
	virtual bool timeLikeVetoed(const Branching &,ShowerParticlePtr);

	/**
	* Vetos for the spacelike shower
	*/
	virtual bool spaceLikeVetoed(const Branching &,ShowerParticlePtr);

	/**
	* Vetos for the spacelike shower
	*/
	virtual bool spaceLikeDecayVetoed(const Branching &,ShowerParticlePtr);

	/**
	* Only generate the hard emission, for testing only.
	*/
	bool hardOnly() const {return _limitEmissions==3;}

	/**
	* Check the flags
	*/
	void checkFlags();

	/**
	*
	*/
	void addFSRUsingDecayPOWHEG(HardTreePtr ISRTree);

	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:

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

	/**
	- * If cascade was called before, but there are still decaying partialces
	- * in the DecayInShower list to be showered. Also possible with another shower
	- * like combining Dipole and Qtilde.
	- */
	- virtual void cascade(tPVector) ;
	-
	- /**
	* Decay a ShowerTree
	*/
	void decay(ShowerTreePtr tree, ShowerDecayMap & decay);

	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:

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

	private:

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

	private:


	/**
	* Stuff from the ShowerHandler
	*/
	//@{

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

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

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

	private :

	/**
	* Pointer to the model for the shower evolution model
	*/
	ShowerModelPtr _model;

	/**
	* Pointer to the splitting generator
	*/
	SplittingGeneratorPtr _splittingGenerator;

	/**
	* Maximum number of tries to generate the shower of a particular tree
	*/
	unsigned int _maxtry;

	/**
	* Matrix element correction switch
	*/
	unsigned int _meCorrMode;

	/**
	* Control of the reconstruction option
	*/
	unsigned int _reconOpt;

	/**
	* If hard veto pT scale is being read-in this determines
	* whether the read-in value is applied to primary and
	* secondary (MPI) scatters or just the primary one, with
	* the usual computation of the veto being performed for
	* the secondary (MPI) scatters.
	*/
	bool _hardVetoReadOption;

	/**
	* rms intrinsic pT of Gaussian distribution
	*/
	Energy _iptrms;

	/**
	* Proportion of inverse quadratic intrinsic pT distribution
	*/
	double _beta;

	/**
	* Parameter for inverse quadratic: 2BetaGamma/(sqr(Gamma)+sqr(intrinsicpT))
	*/
	Energy _gamma;

	/**
	* Upper bound on intrinsic pT for inverse quadratic
	*/
	Energy _iptmax;

	/**
	* Limit the number of emissions for testing
	*/
	unsigned int _limitEmissions;

	/**
	* The progenitor of the current shower
	*/
	ShowerProgenitorPtr _progenitor;

	/**
	* Matrix element
	*/
	HwMEBasePtr _hardme;

	/**
	* Decayer
	*/
	HwDecayerBasePtr _decayme;

	/**
	* The ShowerTree currently being showered
	*/
	ShowerTreePtr _currenttree;

	/**
	* The HardTree currently being showered
	*/
	HardTreePtr _hardtree;

	/**
	* Radiation enhancement factors for use with the veto algorithm
	* if needed by the soft matrix element correction
	*/
	//@{
	/**
	* Enhancement factor for initial-state radiation
	*/
	double _initialenhance;

	/**
	* Enhancement factor for final-state radiation
	*/
	double _finalenhance;
	//@}

	/**
	* The beam particle data for the current initial-state shower
	*/
	Ptr<BeamParticleData>::const_pointer _beam;

	/**
	* Storage of the intrinsic \f$p_t\f$ of the particles
	*/
	map<tShowerProgenitorPtr,pair<Energy,double> > _intrinsic;

	/**
	* Vetoes
	*/
	vector<ShowerVetoPtr> _vetoes;

	/**
	* Full Shower Vetoes
	*/
	vector<FullShowerVetoPtr> _fullShowerVetoes;

	/**
	* Number of iterations for reweighting
	*/
	unsigned int _nReWeight;

	/**
	* Whether or not we are reweighting
	*/
	bool _reWeight;

	/**
	* number of IS emissions
	*/
	unsigned int _nis;

	/**
	* Number of FS emissions
	*/
	unsigned int _nfs;

	/**
	* The option for wqhich interactions to use
	*/
	ShowerInteraction interaction_;

	/**
	* Truncated shower switch
	*/
	bool _trunc_Mode;

	/**
	* Count of the number of truncated emissions
	*/
	unsigned int _truncEmissions;

	/**
	* Mode for the hard emissions
	*/
	int _hardEmission;

	/**
	* Option to include spin correlations
	*/
	unsigned int _spinOpt;

	/**
	* Option for the kernal for soft correlations
	*/
	unsigned int _softOpt;

	/**
	* Option for hard radiation in POWHEG events
	*/
	bool _hardPOWHEG;

	/**
	* True if no warnings about incorrect hard emission
	* mode setting have been issued yet
	*/
	static bool _hardEmissionWarn;

	/**
	* True if no warnings about missing truncated shower
	* have been issued yet
	*/
	static bool _missingTruncWarn;

	/**
	* The relevant hard scale to be used in the profile scales
	*/
	Energy muPt;

	/**
	* Maximum number of emission attempts for FSR
	*/
	unsigned int _maxTryFSR;

	/**
	* Maximum number of failures for FSR generation
	*/
	unsigned int _maxFailFSR;

	/**
	* Failure fraction for FSR generation
	*/
	double _fracFSR;

	/**
	* Counter for number of FSR emissions
	*/
	unsigned int _nFSR;

	/**
	* Counter for the number of failed events due to FSR emissions
	*/
	unsigned int _nFailedFSR;

	};

	}

	#endif /* HERWIG_QTildeShowerHandler_H */
	diff --git a/Shower/ShowerHandler.cc b/Shower/ShowerHandler.cc
	--- a/Shower/ShowerHandler.cc
	+++ b/Shower/ShowerHandler.cc
	@@ -1,1096 +1,1081 @@
	// -- 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 "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/QTilde/Base/ShowerTree.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() {}

	tShowerHandlerPtr ShowerHandler::currentHandler_ = tShowerHandlerPtr();

	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::_vmin2 = vMin_;
	ShowerTree::_spaceTime = includeSpaceTime_;
	}

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

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

	ShowerHandler::ShowerHandler() :
	maxtry_(10),maxtryMPI_(10),maxtryDP_(10),maxtryDecay_(100),
	factorizationScaleFactor_(1.0),
	renormalizationScaleFactor_(1.0),
	hardScaleFactor_(1.0),
	restrictPhasespace_(true), maxPtIsMuF_(false),
	pdfFreezingScale_(2.5*GeV),
	doFSR_(true), doISR_(true),
	splitHardProcess_(true),
	includeSpaceTime_(false), vMin_(0.1*GeV2),
	reweight_(1.0) {
	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::_vmin2 = vMin_;
	ShowerTree::_spaceTime = includeSpaceTime_;
	}

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

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

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

	void ShowerHandler::Init() {

	static ClassDocumentation<ShowerHandler> documentation
	("Main driver class for the showering.");

	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 Parameter<ShowerHandler,unsigned int> interfaceMaxTryDecay
	("MaxTryDecay",
	"The maximum number of attempts to generate a decay",
	&ShowerHandler::maxtryDecay_, 200, 10, 0,
	false, false, Interface::lowerlim);

	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 Command<ShowerHandler> interfaceAddVariation
	("AddVariation",
	"Add a shower variation.",
	&ShowerHandler::doAddVariation, false);

	static Switch<ShowerHandler,bool> interfaceDoFSR
	("DoFSR",
	"Switch on or off final state radiation.",
	&ShowerHandler::doFSR_, true, false, false);
	static SwitchOption interfaceDoFSROn
	(interfaceDoFSR,
	"Yes",
	"Switch on final state radiation.",
	true);
	static SwitchOption interfaceDoFSROff
	(interfaceDoFSR,
	"No",
	"Switch off final state radiation.",
	false);

	static Switch<ShowerHandler,bool> interfaceDoISR
	("DoISR",
	"Switch on or off initial state radiation.",
	&ShowerHandler::doISR_, true, false, false);
	static SwitchOption interfaceDoISROn
	(interfaceDoISR,
	"Yes",
	"Switch on initial state radiation.",
	true);
	static SwitchOption interfaceDoISROff
	(interfaceDoISR,
	"No",
	"Switch off initial state radiation.",
	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);
	}

	Energy ShowerHandler::hardScale() const {
	assert(false);
	}

	void ShowerHandler::cascade() {
	useMe();
	// Initialise the weights in the event object
	// so that any variations are output regardless of
	// whether showering occurs for the given event
	initializeWeights();
	// get the PDF's from ThePEG (if locally overridden use the local versions)
	tcPDFPtr first = PDFA_ ? tcPDFPtr(PDFA_) : firstPDF().pdf();
	tcPDFPtr second = PDFB_ ? tcPDFPtr(PDFB_) : secondPDF().pdf();
	resetPDFs(make_pair(first,second));
	-
	- if(didRunCascade()){
	- tPVector finalstate=eventHandler()->currentStep()->getFinalState();
	- bool issecondshower=false;
	- for(tPVector::iterator it=finalstate.begin();it!=finalstate.end();++it){
	- if((decaysInShower((it).id())&&!(it).dataPtr()->stable())){
	- issecondshower=true;
	- break;
	- }
	- }
	- if (issecondshower) cascade(finalstate);
	- return;
	- }
	-
	-
	// set the PDFs for the remnant
	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();
	remnantDecayer()->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
	setCurrentHandler();
	// first shower the hard process
	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);
	// perform the reweighting for the hard process shower
	combineWeights();
	// unset the current ShowerHandler
	unSetCurrentHandler();
	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
	remnantDecayer()->initialize(remnants, incoming_, *currentStep(), pdfFreezingScale());
	// do the first forcedSplitting
	try {
	remnantDecayer()->doSplit(incomingPartons, make_pair(rempdfs_.first,rempdfs_.second), true);
	}
	catch (ExtraScatterVeto) {
	throw Exception() << "Remnant extraction failed in "
	<< "ShowerHandler::cascade() from primary interaction"
	<< Exception::eventerror;
	}
	// perform the reweighting for the hard process shower
	combineWeights();
	// if no MPI return
	if( !isMPIOn() ) {
	remnantDecayer()->finalize();
	// boost back to lab if needed
	if(btotal) boostCollision(true);
	// unset the current ShowerHandler
	unSetCurrentHandler();
	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
	remnantDecayer()->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;
	}
	// perform the reweighting for the additional hard scatter shower
	combineWeights();
	}
	// 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
	remnantDecayer()->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;
	// perform the reweighting for the MPI process shower
	combineWeights();
	}
	// finalize the remnants
	remnantDecayer()->finalize(getMPIHandler()->colourDisrupt(),
	getMPIHandler()->softMultiplicity());
	// boost back to lab if needed
	if(btotal) boostCollision(true);
	// unset the current ShowerHandler
	unSetCurrentHandler();
	getMPIHandler()->clean();
	}

	void ShowerHandler::initializeWeights() {
	if ( !showerVariations().empty() ) {

	tEventPtr event = eventHandler()->currentEvent();

	for ( map<string,ShowerVariation>::const_iterator var =
	showerVariations().begin();
	var != showerVariations().end(); ++var ) {

	// Check that this is behaving as intended
	//map<string,double>::iterator wi = event->optionalWeights().find(var->first);
	//assert(wi == event->optionalWeights().end() );

	event->optionalWeights()[var->first] = 1.0;
	currentWeights_[var->first] = 1.0;
	}
	}
	reweight_ = 1.0;
	}

	void ShowerHandler::resetWeights() {
	for ( map<string,double>::iterator w = currentWeights_.begin();
	w != currentWeights_.end(); ++w ) {
	w->second = 1.0;
	}
	reweight_ = 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 {
	assert(false && "Weight name unknown.");
	//event->optionalWeights()[w->first] = w->second;
	}
	}
	if ( reweight_ != 1.0 ) {
	Ptr<StandardEventHandler>::tptr eh =
	dynamic_ptr_cast<Ptr<StandardEventHandler>::tptr>(eventHandler());
	if ( !eh ) {
	throw Exception() << "ShowerHandler::combineWeights() : Cross section reweighting "
	<< "through the shower is currently only available with standard "
	<< "event generators" << Exception::runerror;
	}
	eh->reweight(reweight_);
	}
	}

	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, XCPtr) {
	assert(false);
	}

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

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

	namespace {

	bool decayProduct(tSubProPtr subProcess,
	tPPtr particle) {
	// must be time-like and not incoming
	if(particle->momentum().m2()<=ZERO\|\|
	particle == subProcess->incoming().first\|\|
	particle == subProcess->incoming().second) return false;
	// if only 1 outgoing and this is it
	if(subProcess->outgoing().size()==1 &&
	subProcess->outgoing()[0]==particle) return true;
	// must not be the s-channel intermediate otherwise
	if(find(subProcess->incoming().first->children().begin(),
	subProcess->incoming().first->children().end(),particle)!=
	subProcess->incoming().first->children().end()&&
	find(subProcess->incoming().second->children().begin(),
	subProcess->incoming().second->children().end(),particle)!=
	subProcess->incoming().second->children().end()&&
	subProcess->incoming().first ->children().size()==1&&
	subProcess->incoming().second->children().size()==1)
	return false;
	// if non-coloured this is enough
	if(!particle->dataPtr()->coloured()) return true;
	// if coloured must be unstable
	if(particle->dataPtr()->stable()) return false;
	// must not have same particle type as a child
	int id = particle->id();
	for(unsigned int ix=0;ix<particle->children().size();++ix)
	if(particle->children()[ix]->id()==id) return false;
	// otherwise its a decaying particle
	return true;
	}

	PPtr findParent(PPtr original, bool & isHard,
	set<PPtr> outgoingset,
	tSubProPtr subProcess) {
	PPtr parent=original;
	isHard \|=(outgoingset.find(original) != outgoingset.end());
	if(!original->parents().empty()) {
	PPtr orig=original->parents()[0];
	if(CurrentGenerator::current().currentEventHandler()->currentStep()->
	find(orig)&&decayProduct(subProcess,orig)) {
	parent=findParent(orig,isHard,outgoingset,subProcess);
	}
	}
	return parent;
	}
	}

	void ShowerHandler::findDecayProducts(PPtr in,PerturbativeProcessPtr hard,
	DecayProcessMap decay) const {
	ParticleVector children=in->children();
	for(ParticleVector::const_iterator it=children.begin(); it!=children.end();++it) {
	// if decayed or should be decayed in shower make the PerturbaitveProcess
	bool radiates = false;
	if(!(**it).children().empty()) {
	// remove d,u,s,c,b quarks and leptons other than on-shell taus
	if( StandardQCDPartonMatcher::Check((**it).id()) \|\|
	( LeptonMatcher::Check((it).id()) && !(abs((it).id())==ParticleID::tauminus &&
	abs((it).mass()-(it).dataPtr()->mass())<MeV))) {
	radiates = true;
	}
	else {
	bool foundParticle(false),foundGauge(false);
	for(unsigned int iy=0;iy<(**it).children().size();++iy) {
	if((it).children()[iy]->id()==(it).id()) {
	foundParticle = true;
	}
	else if((**it).children()[iy]->id()==ParticleID::g \|\|
	(**it).children()[iy]->id()==ParticleID::gamma) {
	foundGauge = true;
	}
	}
	radiates = foundParticle && foundGauge;
	}
	}
	if(radiates) {
	findDecayProducts(*it,hard,decay);
	}
	else if(!(**it).children().empty()\|\|
	(decaysInShower((it).id())&&!(it).dataPtr()->stable())) {
	createDecayProcess(in,hard,decay);
	}
	else {
	hard->outgoing().push_back(make_pair(*it,PerturbativeProcessPtr()));
	}
	}
	}

	void ShowerHandler::splitHardProcess(tPVector tagged, PerturbativeProcessPtr & hard,
	DecayProcessMap & decay) const {
	// temporary storage of the particles
	set<PPtr> hardParticles;
	// tagged particles in a set
	set<PPtr> outgoingset(tagged.begin(),tagged.end());
	bool isHard=false;
	// loop over the tagged particles
	for (tParticleVector::const_iterator taggedP = tagged.begin();
	taggedP != tagged.end(); ++taggedP) {
	// skip remnants
	if (eventHandler()->currentCollision()&&
	eventHandler()->currentCollision()->isRemnant(*taggedP)) continue;
	// find the parent and whether its a decaying particle
	bool isDecayProd=false;
	// check if hard
	isHard \|=(outgoingset.find(*taggedP) != outgoingset.end());
	if(splitHardProcess_) {
	tPPtr parent = *taggedP;
	// check if from s channel decaying colourless particle
	while(parent&&!parent->parents().empty()&&!isDecayProd) {
	parent = parent->parents()[0];
	if(parent == subProcess_->incoming().first \|\|
	parent == subProcess_->incoming().second ) break;
	isDecayProd = decayProduct(subProcess_,parent);
	}
	if (isDecayProd)
	hardParticles.insert(findParent(parent,isHard,outgoingset,subProcess_));
	}
	if (!isDecayProd)
	hardParticles.insert(*taggedP);
	}
	// there must be something to shower
	if(hardParticles.empty())
	throw Exception() << "No particles to shower in "
	<< "ShowerHandler::splitHardProcess()"
	<< Exception::eventerror;
	// must be a hard process
	if(!isHard)
	throw Exception() << "Starting on decay not yet implemented in "
	<< "ShowerHandler::splitHardProcess()"
	<< Exception::runerror;
	// create the hard process
	hard = new_ptr(PerturbativeProcess());
	// incoming particles
	hard->incoming().push_back(make_pair(subProcess_->incoming().first ,PerturbativeProcessPtr()));
	hard->incoming().push_back(make_pair(subProcess_->incoming().second,PerturbativeProcessPtr()));
	// outgoing particles
	for(set<PPtr>::const_iterator it=hardParticles.begin();it!=hardParticles.end();++it) {
	// if decayed or should be decayed in shower make the tree
	PPtr orig = *it;
	bool radiates = false;
	if(!orig->children().empty()) {
	// remove d,u,s,c,b quarks and leptons other than on-shell taus
	if( StandardQCDPartonMatcher::Check(orig->id()) \|\|
	( LeptonMatcher::Check(orig->id()) &&
	!(abs(orig->id())==ParticleID::tauminus && abs(orig->mass()-orig->dataPtr()->mass())<MeV))) {
	radiates = true;
	}
	else {
	bool foundParticle(false),foundGauge(false);
	for(unsigned int iy=0;iy<orig->children().size();++iy) {
	if(orig->children()[iy]->id()==orig->id()) {
	foundParticle = true;
	}
	else if(orig->children()[iy]->id()==ParticleID::g \|\|
	orig->children()[iy]->id()==ParticleID::gamma) {
	foundGauge = true;
	}
	}
	radiates = foundParticle && foundGauge;
	}
	}
	if(radiates) {
	findDecayProducts(orig,hard,decay);
	}
	else if(!(**it).children().empty()\|\|
	(decaysInShower((it).id())&&!(it).dataPtr()->stable())) {
	createDecayProcess(*it,hard,decay);
	}
	else {
	hard->outgoing().push_back(make_pair(*it,PerturbativeProcessPtr()));
	}
	}
	}

	void ShowerHandler::createDecayProcess(PPtr in,PerturbativeProcessPtr hard, DecayProcessMap & decay) const {
	// there must be an incoming particle
	assert(in);
	// create the new process and connect with the parent
	PerturbativeProcessPtr newDecay=new_ptr(PerturbativeProcess());
	newDecay->incoming().push_back(make_pair(in,hard));
	Energy width=in->dataPtr()->generateWidth(in->mass());
	decay.insert(make_pair(width,newDecay));
	hard->outgoing().push_back(make_pair(in,newDecay));
	// we need to deal with the decay products if decayed
	ParticleVector children = in->children();
	if(!children.empty()) {
	for(ParticleVector::const_iterator it = children.begin();
	it!= children.end(); ++it) {
	// if decayed or should be decayed in shower make the tree
	in->abandonChild(*it);
	bool radiates = false;
	if(!(**it).children().empty()) {
	if(StandardQCDPartonMatcher::Check((**it).id())\|\|
	(LeptonMatcher::Check((it).id())&& !(abs((it).id())==ParticleID::tauminus &&
	abs((it).mass()-(it).dataPtr()->mass())<MeV))) {
	radiates = true;
	}
	else {
	bool foundParticle(false),foundGauge(false);
	for(unsigned int iy=0;iy<(**it).children().size();++iy) {
	if((it).children()[iy]->id()==(it).id()) {
	foundParticle = true;
	}
	else if((**it).children()[iy]->id()==ParticleID::g \|\|
	(**it).children()[iy]->id()==ParticleID::gamma) {
	foundGauge = true;
	}
	}
	radiates = foundParticle && foundGauge;
	}
	// finally assume all non-decaying particles are in this class
	// pr 27/11/15 not sure about this bit
	// if(!radiates) {
	// radiates = !decaysInShower((**it).id());
	// }
	}
	if(radiates) {
	findDecayProducts(*it,newDecay,decay);
	}
	else if(!(**it).children().empty()\|\|
	(decaysInShower((it).id())&&!(it).dataPtr()->stable())) {
	createDecayProcess(*it,newDecay,decay);
	}
	else {
	newDecay->outgoing().push_back(make_pair(*it,PerturbativeProcessPtr()));
	}
	}
	}
	}

	tDMPtr ShowerHandler::decay(PerturbativeProcessPtr process,
	DecayProcessMap & decayMap) const {
	PPtr parent = process->incoming()[0].first;
	assert(parent);
	if(parent->spinInfo()) parent->spinInfo()->decay(true);
	unsigned int ntry = 0;
	ParticleVector children;
	tDMPtr dm = DMPtr();
	while (true) {
	// exit if fails
	if (++ntry>=maxtryDecay_)
	throw Exception() << "Failed to perform decay in ShowerHandler::decay()"
	<< " after " << maxtryDecay_
	<< " attempts for " << parent->PDGName()
	<< Exception::eventerror;
	// select decay mode
	dm = parent->data().selectMode(*parent);

	if(!dm)
	throw Exception() << "Failed to select decay mode in ShowerHandler::decay()"
	<< "for " << parent->PDGName()
	<< Exception::eventerror;
	if(!dm->decayer())
	throw Exception() << "No Decayer for selected decay mode "
	<< " in ShowerHandler::decay()"
	<< Exception::runerror;
	// start of try block
	try {
	children = dm->decayer()->decay(dm, parent);
	// if no children have another go
	if(children.empty()) continue;
	// set up parent
	parent->decayMode(dm);
	// add children
	for (unsigned int i = 0, N = children.size(); i < N; ++i ) {
	children[i]->setLabVertex(parent->labDecayVertex());
	//parent->addChild(children[i]);
	}
	// if succeeded break out of loop
	break;
	}
	catch(Veto) {
	}
	}
	assert(!children.empty());
	for(ParticleVector::const_iterator it = children.begin();
	it!= children.end(); ++it) {
	if(!(**it).children().empty()\|\|
	(decaysInShower((it).id())&&!(it).dataPtr()->stable())) {
	createDecayProcess(*it,process,decayMap);
	}
	else {
	process->outgoing().push_back(make_pair(*it,PerturbativeProcessPtr()));
	}
	}
	return dm;
	}


	// Note: The tag must be constructed from an ordered particle container.
	tDMPtr ShowerHandler::findDecayMode(const string & tag) const {
	static map<string,DMPtr> cache;
	map<string,DMPtr>::const_iterator pos = cache.find(tag);

	if ( pos != cache.end() )
	return pos->second;

	tDMPtr dm = CurrentGenerator::current().findDecayMode(tag);
	cache[tag] = dm;
	return dm;
	}


	/**
	* Operator for the particle ordering
	* @param p1 The first ParticleData object
	* @param p2 The second ParticleData object
	*/

	bool ShowerHandler::ParticleOrdering::operator() (tcPDPtr p1, tcPDPtr p2) {
	return abs(p1->id()) > abs(p2->id()) \|\|
	( abs(p1->id()) == abs(p2->id()) && p1->id() > p2->id() ) \|\|
	( p1->id() == p2->id() && p1->fullName() > p2->fullName() );
	}

	diff --git a/Shower/ShowerHandler.h b/Shower/ShowerHandler.h
	--- a/Shower/ShowerHandler.h
	+++ b/Shower/ShowerHandler.h
	@@ -1,814 +1,803 @@
	// -- 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 "ShowerVariation.h"
	#include "Herwig/PDF/HwRemDecayer.fh"
	#include "ThePEG/EventRecord/RemnantParticle.fh"
	#include "UEBase.h"
	#include "PerturbativeProcess.h"
	#include "Herwig/MatrixElement/Matchbox/Matching/HardScaleProfile.h"
	#include "ShowerHandler.fh"

	namespace Herwig {
	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;

	public:

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

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

	public:

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

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

	/**
	* Get the PDF freezing scale
	*/
	Energy pdfFreezingScale() const { return pdfFreezingScale_; }

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

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

	/**
	* Split the hard process into production and decays
	* @param tagged The tagged particles from the StepHandler
	* @param hard The hard perturbative process
	* @param decay The decay particles
	*/
	void splitHardProcess(tPVector tagged, PerturbativeProcessPtr & hard,
	DecayProcessMap & decay) const;

	/**
	* Decay a particle
	*/
	tDMPtr decay(PerturbativeProcessPtr,
	DecayProcessMap & decay) const;


	/**
	* Cached lookup of decay modes.
	* Generator::findDecayMode() is not efficient.
	*/
	tDMPtr findDecayMode(const string & tag) const;


	/**
	* A struct to order the particles in the same way as in the DecayMode's
	*/

	struct ParticleOrdering {

	bool operator() (tcPDPtr p1, tcPDPtr p2);

	};


	/**
	* A container for ordered particles required
	* for constructing tags for decay mode lookup.
	*/
	typedef multiset<tcPDPtr,ParticleOrdering> OrderedParticles;


	public:

	/**
	* @name Switches for initial- and final-state radiation
	*/
	//@{
	/**
	* Switch for any radiation
	*/
	bool doRadiation() const {return doFSR_ \|\| doISR_;}

	/**
	* Switch on or off final state radiation.
	*/
	bool doFSR() const { return doFSR_;}

	/**
	* Switch on or off initial state radiation.
	*/
	bool doISR() const { return doISR_;}
	//@}

	public:

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

	/**
	* The factorization scale factor.
	*/
	double factorizationScaleFactor() const {
	return factorizationScaleFactor_;
	}

	/**
	* The renormalization scale factor.
	*/
	double renFac() const {
	return renormalizationScaleFactor_;
	}

	/**
	* The factorization scale factor.
	*/
	double facFac() const {
	return factorizationScaleFactor_;
	}

	/**
	* The renormalization scale factor.
	*/
	double renormalizationScaleFactor() const {
	return renormalizationScaleFactor_;
	}

	/**
	* The scale factor for the hard scale
	*/
	double hardScaleFactor() const {
	return hardScaleFactor_;
	}
	/**
	* 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;
	//@}

	public:

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

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

	/**
	* Access the current Weights
	*/
	map<string,double>& currentWeights() {
	return currentWeights_;
	}

	/**
	* Return the current Weights
	*/
	const map<string,double>& currentWeights() const {
	return currentWeights_;
	}

	/**
	* Change the current reweighting factor
	*/
	void reweight(double w) {
	reweight_ = w;
	}

	/**
	* Return the current reweighting factor
	*/
	double reweight() const {
	return reweight_;
	}

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

	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 Functions to perform the cascade
	*/
	//@{
	/**
	* The main method which manages the showering of a subprocess.
	*/
	virtual tPPair cascade(tSubProPtr sub, XCPtr xcomb);

	/**
	- * If cascade was called before, but there are still decaying partialces
	- * in the DecayInShower list to be showered. Also possible with another shower
	- * like combining Dipole and Qtilde.
	- */
	- virtual void cascade(tPVector) {
	- throw Exception()
	- <<"\nError: Cascade Handler was called multiple times."
	- << Exception::eventerror;
	- };
	-
	- /**
	* Set up for the cascade
	*/
	void prepareCascade(tSubProPtr sub) {
	current_ = currentStep();
	subProcess_ = sub;
	}

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

	protected:

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

	/**
	* Unset the current handler
	*/
	void unSetCurrentHandler() {
	currentHandler_ = tShowerHandlerPtr();
	}
	//@}

	protected:

	/**
	* @name Members relating to the underlying event and MPI
	*/
	//@{
	/**
	* Return true if multiple parton interactions are switched on
	* and can be used for this beam setup.
	*/
	bool isMPIOn() const {
	return MPIHandler_ && MPIHandler_->beamOK();
	}

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

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

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

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

	public:

	/**
	* Check if a particle decays in the shower
	* @param id The PDG code for the particle
	*/
	bool decaysInShower(long id) const {
	return ( particlesDecayInShower_.find( abs(id) ) !=
	particlesDecayInShower_.end() );
	}

	protected:

	/**
	* Members to handle splitting up of hard process and decays
	*/
	//@{
	/**
	* Find decay products from the hard process and create decay processes
	* @param parent The parent particle
	* @param hard The hard process
	* @param decay The decay processes
	*/
	void findDecayProducts(PPtr parent, PerturbativeProcessPtr hard, DecayProcessMap decay) const;

	/**
	* Find decay products from the hard process and create decay processes
	* @param parent The parent particle
	* @param hard The parent hard process
	* @param decay The decay processes
	*/
	void createDecayProcess(PPtr parent,PerturbativeProcessPtr hard, DecayProcessMap & decay) const;
	//@}

	/**
	* @name Functions to return information relevant to the process being showered
	*/
	//@{
	/**
	* Return the currently used SubProcess.
	*/
	tSubProPtr currentSubProcess() const {
	assert(subProcess_);
	return subProcess_;
	}

	/**
	* Access to the incoming beam particles
	*/
	tPPair incomingBeams() const {
	return incoming_;
	}
	//@}

	protected:

	/**
	* Weight handling for shower variations
	*/
	//@
	/**
	* Combine the variation weights which have been encountered
	*/
	void combineWeights();

	/**
	* Initialise the weights in currentEvent()
	*/
	void initializeWeights();

	/**
	* Reset the current weights
	*/
	void resetWeights();
	//@}

	protected:

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

	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:

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

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

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

	private:

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

	private:

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

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

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

	private:

	/**
	* Maximum tries for various stages of the showering process
	*/
	//@{
	/**
	* 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_;

	/**
	* Maximum number of attempts to generate a decay
	*/
	unsigned int maxtryDecay_;
	//@}

	private:

	/**
	* Factors for the various scales
	*/
	//@{
	/**
	* The factorization scale factor.
	*/
	double factorizationScaleFactor_;

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

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

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

	private:

	/**
	* Storage of information about the current event
	*/
	//@{
	/**
	* The incoming beam particles for the current collision
	*/
	tPPair incoming_;

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

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

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

	private:

	/**
	* PDFs to be used for the various stages and related parameters
	*/
	//@{
	/**
	* The PDF freezing scale
	*/
	Energy pdfFreezingScale_;

	/**
	* PDFs to be used for the various stages and related parameters
	*/
	//@{
	/**
	* 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 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_;
	//@}

	private:

	/**
	* @name Parameters for initial- and final-state radiation
	*/
	//@{
	/**
	* Switch on or off final state radiation.
	*/
	bool doFSR_;

	/**
	* Switch on or off initial state radiation.
	*/
	bool doISR_;
	//@}

	private:

	/**
	* @name Parameters for particle decays
	*/
	//@{
	/**
	* Whether or not to split into hard and decay trees
	*/
	bool splitHardProcess_;

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

	private:

	/**
	* Parameters for the space-time model
	*/
	//@{
	/**
	* Whether or not to include spa-cetime distances in the shower
	*/
	bool includeSpaceTime_;

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

	private:

	/**
	* Parameters relevant for reweight and variations
	*/
	//@{
	/**
	* The shower variations
	*/
	map<string,ShowerVariation> showerVariations_;

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

	/**
	* A reweighting factor applied by the showering
	*/
	double reweight_;

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

	}

	#endif /* HERWIG_ShowerHandler_H */

File Metadata

Mime Type: text/x-diff
Expires: Sat, May 3, 6:34 AM (1 d, 17 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4972230
Default Alt Text: (220 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions