No OneTemporary
Actions

Size

101 KB

Subscribers

None

View Options

	diff --git a/Shower/Base/SudakovFormFactor.cc b/Shower/Base/SudakovFormFactor.cc
	--- a/Shower/Base/SudakovFormFactor.cc
	+++ b/Shower/Base/SudakovFormFactor.cc
	@@ -1,672 +1,675 @@
	// -- C++ --
	//
	// SudakovFormFactor.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 SudakovFormFactor class.
	//

	#include "SudakovFormFactor.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ShowerKinematics.h"
	#include "ShowerParticle.h"
	#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
	#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
	#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
	#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
	#include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h"
	#include "ThePEG/Utilities/DescribeClass.h"

	using namespace Herwig;

	DescribeAbstractClass<SudakovFormFactor,Interfaced>
	describeSudakovFormFactor ("Herwig::SudakovFormFactor","");

	void SudakovFormFactor::persistentOutput(PersistentOStream & os) const {
	os << splittingFn_ << alpha_ << pdfmax_ << particles_ << pdffactor_
	<< a_ << b_ << ounit(c_,GeV) << ounit(kinCutoffScale_,GeV) << cutOffOption_
	<< ounit(vgcut_,GeV) << ounit(vqcut_,GeV)
	<< ounit(pTmin_,GeV) << ounit(pT2min_,GeV2);
	}

	void SudakovFormFactor::persistentInput(PersistentIStream & is, int) {
	is >> splittingFn_ >> alpha_ >> pdfmax_ >> particles_ >> pdffactor_
	>> a_ >> b_ >> iunit(c_,GeV) >> iunit(kinCutoffScale_,GeV) >> cutOffOption_
	>> iunit(vgcut_,GeV) >> iunit(vqcut_,GeV)
	>> iunit(pTmin_,GeV) >> iunit(pT2min_,GeV2);
	}

	void SudakovFormFactor::Init() {

	static ClassDocumentation<SudakovFormFactor> documentation
	("The SudakovFormFactor class is the base class for the implementation of Sudakov"
	" form factors in Herwig++");

	static Reference<SudakovFormFactor,SplittingFunction>
	interfaceSplittingFunction("SplittingFunction",
	"A reference to the SplittingFunction object",
	&Herwig::SudakovFormFactor::splittingFn_,
	false, false, true, false);

	static Reference<SudakovFormFactor,ShowerAlpha>
	interfaceAlpha("Alpha",
	"A reference to the Alpha object",
	&Herwig::SudakovFormFactor::alpha_,
	false, false, true, false);

	static Parameter<SudakovFormFactor,double> interfacePDFmax
	("PDFmax",
	"Maximum value of PDF weight. ",
	&SudakovFormFactor::pdfmax_, 35.0, 1.0, 100000.0,
	false, false, Interface::limited);

	static Switch<SudakovFormFactor,unsigned int> interfacePDFFactor
	("PDFFactor",
	"Include additional factors in the overestimate for the PDFs",
	&SudakovFormFactor::pdffactor_, 0, false, false);
	static SwitchOption interfacePDFFactorOff
	(interfacePDFFactor,
	"Off",
	"Don't include any factors",
	0);
	static SwitchOption interfacePDFFactorOverZ
	(interfacePDFFactor,
	"OverZ",
	"Include an additional factor of 1/z",
	1);
	static SwitchOption interfacePDFFactorOverOneMinusZ
	(interfacePDFFactor,
	"OverOneMinusZ",
	"Include an additional factor of 1/(1-z)",
	2);
	static SwitchOption interfacePDFFactorOverZOneMinusZ
	(interfacePDFFactor,
	"OverZOneMinusZ",
	"Include an additional factor of 1/z/(1-z)",
	3);


	static Switch<SudakovFormFactor,unsigned int> interfaceCutOffOption
	("CutOffOption",
	"The type of cut-off to use to end the shower",
	&SudakovFormFactor::cutOffOption_, 0, false, false);
	static SwitchOption interfaceCutOffOptionDefault
	(interfaceCutOffOption,
	"Default",
	"Use the standard Herwig++ cut-off on virtualities with the minimum"
	" virtuality depending on the mass of the branching particle",
	0);
	static SwitchOption interfaceCutOffOptionFORTRAN
	(interfaceCutOffOption,
	"FORTRAN",
	"Use a FORTRAN-like cut-off on virtualities",
	1);
	static SwitchOption interfaceCutOffOptionpT
	(interfaceCutOffOption,
	"pT",
	"Use a cut on the minimum allowed pT",
	2);

	static Parameter<SudakovFormFactor,double> interfaceaParameter
	("aParameter",
	"The a parameter for the kinematic cut-off",
	&SudakovFormFactor::a_, 0.3, -10.0, 10.0,
	false, false, Interface::limited);

	static Parameter<SudakovFormFactor,double> interfacebParameter
	("bParameter",
	"The b parameter for the kinematic cut-off",
	&SudakovFormFactor::b_, 2.3, -10.0, 10.0,
	false, false, Interface::limited);

	static Parameter<SudakovFormFactor,Energy> interfacecParameter
	("cParameter",
	"The c parameter for the kinematic cut-off",
	&SudakovFormFactor::c_, GeV, 0.3GeV, 0.1GeV, 10.0*GeV,
	false, false, Interface::limited);

	static Parameter<SudakovFormFactor,Energy>
	interfaceKinScale ("cutoffKinScale",
	"kinematic cutoff scale for the parton shower phase"
	" space (unit [GeV])",
	&SudakovFormFactor::kinCutoffScale_, GeV,
	2.3GeV, 0.001GeV, 10.0*GeV,false,false,false);

	static Parameter<SudakovFormFactor,Energy> interfaceGluonVirtualityCut
	("GluonVirtualityCut",
	"For the FORTRAN cut-off option the minimum virtuality of the gluon",
	&SudakovFormFactor::vgcut_, GeV, 0.85GeV, 0.1GeV, 10.0*GeV,
	false, false, Interface::limited);

	static Parameter<SudakovFormFactor,Energy> interfaceQuarkVirtualityCut
	("QuarkVirtualityCut",
	"For the FORTRAN cut-off option the minimum virtuality added to"
	" the mass for particles other than the gluon",
	&SudakovFormFactor::vqcut_, GeV, 0.85GeV, 0.1GeV, 10.0*GeV,
	false, false, Interface::limited);

	static Parameter<SudakovFormFactor,Energy> interfacepTmin
	("pTmin",
	"The minimum pT if using a cut-off on the pT",
	&SudakovFormFactor::pTmin_, GeV, 1.0GeV, ZERO, 10.0GeV,
	false, false, Interface::limited);
	}

	bool SudakovFormFactor::
	PDFVeto(const Energy2 t, const double x,
	const tcPDPtr parton0, const tcPDPtr parton1,
	Ptr<BeamParticleData>::transient_const_pointer beam) const {
	assert(pdf_);

	Energy2 theScale = t;
	if (theScale < sqr(freeze_)) theScale = sqr(freeze_);

	double newpdf(0.0), oldpdf(0.0);
	//different treatment of MPI ISR is done via CascadeHandler::resetPDFs()
	newpdf=pdf_->xfx(beam,parton0,theScale,x/z());
	oldpdf=pdf_->xfx(beam,parton1,theScale,x);


	if(newpdf<=0.) return true;
	if(oldpdf<=0.) return false;
	double ratio = newpdf/oldpdf;

	double maxpdf(pdfmax_);
	switch (pdffactor_) {
	case 1:
	maxpdf /= z();
	break;
	case 2:
	maxpdf /= 1.-z();
	break;
	case 3:
	maxpdf /= (z()*(1.-z()));
	break;
	}

	// ratio / PDFMax must be a probability <= 1.0
	if (ratio > maxpdf) {
	generator()->log() << "PDFVeto warning: Ratio > " << name()
	<< ":PDFmax (by a factor of "
	<< ratio/maxpdf <<") for "
	<< parton0->PDGName() << " to "
	<< parton1->PDGName() << "\n";
	}
	return ratio < UseRandom::rnd()*maxpdf;
	}

	void SudakovFormFactor::addSplitting(const IdList & in) {
	bool add=true;
	for(unsigned int ix=0;ix<particles_.size();++ix) {
	if(particles_[ix].size()==in.size()) {
	bool match=true;
	for(unsigned int iy=0;iy<in.size();++iy) {
	if(particles_[ix][iy]!=in[iy]) {
	match=false;
	break;
	}
	}
	if(match) {
	add=false;
	break;
	}
	}
	}
	if(add) particles_.push_back(in);
	}

	namespace {

	LorentzRotation boostToShower(const vector<Lorentz5Momentum> & basis,
	ShowerKinematics::Frame frame,
	Lorentz5Momentum & porig) {
	LorentzRotation output;
	if(frame==ShowerKinematics::BackToBack) {
	// we are doing the evolution in the back-to-back frame
	// work out the boostvector
	Boost boostv(-(basis[0]+basis[1]).boostVector());
	// momentum of the parton
	Lorentz5Momentum ptest(basis[0]);
	// construct the Lorentz boost
	output = LorentzRotation(boostv);
	ptest *= output;
	Axis axis(ptest.vect().unit());
	// now rotate so along the z axis as needed for the splitting functions
	if(axis.perp2()>0.) {
	double sinth(sqrt(1.-sqr(axis.z())));
	output.rotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
	}
	else if(axis.z()<0.) {
	output.rotate(Constants::pi,Axis(1.,0.,0.));
	}
	porig = output*basis[0];
	porig.setX(ZERO);
	porig.setY(ZERO);
	}
	else {
	output = LorentzRotation(-basis[0].boostVector());
	porig = output*basis[0];
	porig.setX(ZERO);
	porig.setY(ZERO);
	porig.setZ(ZERO);
	}
	return output;
	}

	RhoDMatrix bosonMapping(ShowerParticle & particle,
	const Lorentz5Momentum & porig,
	VectorSpinPtr vspin,
	const LorentzRotation & rot) {
	// rotate the original basis
	vector<LorentzPolarizationVector> sbasis;
	for(unsigned int ix=0;ix<3;++ix) {
	sbasis.push_back(vspin->getProductionBasisState(ix));
	sbasis.back().transform(rot);
	}
	// splitting basis
	vector<LorentzPolarizationVector> fbasis;
	bool massless(particle.id()==ParticleID::g\|\|particle.id()==ParticleID::gamma);
	VectorWaveFunction wave(porig,particle.dataPtr(),outgoing);
	for(unsigned int ix=0;ix<3;++ix) {
	if(massless&&ix==1) {
	fbasis.push_back(LorentzPolarizationVector());
	}
	else {
	wave.reset(ix);
	fbasis.push_back(wave.wave());
	}
	}
	// work out the mapping
	RhoDMatrix mapping=RhoDMatrix(PDT::Spin1,false);
	for(unsigned int ix=0;ix<3;++ix) {
	for(unsigned int iy=0;iy<3;++iy) {
	mapping(ix,iy)= sbasis[iy].dot(fbasis[ix].conjugate());
	if(particle.id()<0)
	mapping(ix,iy)=conj(mapping(ix,iy));
	}
	}
	// \todo need to fix this
	mapping = RhoDMatrix(PDT::Spin1,false);
	if(massless) {
	mapping(0,0) = 1.;
	mapping(2,2) = 1.;
	}
	else {
	mapping(0,0) = 1.;
	mapping(1,1) = 1.;
	mapping(2,2) = 1.;
	}
	return mapping;
	}

	RhoDMatrix fermionMapping(ShowerParticle & particle,
	const Lorentz5Momentum & porig,
	FermionSpinPtr fspin,
	const LorentzRotation & rot) {
	// extract the original basis states
	vector<LorentzSpinor<SqrtEnergy> > sbasis;
	for(unsigned int ix=0;ix<2;++ix) {
	sbasis.push_back(fspin->getProductionBasisState(ix));
	sbasis.back().transform(rot);
	}
	// calculate the states in the splitting basis
	vector<LorentzSpinor<SqrtEnergy> > fbasis;
	SpinorWaveFunction wave(porig,particle.dataPtr(),
	particle.id()>0 ? incoming : outgoing);
	for(unsigned int ix=0;ix<2;++ix) {
	wave.reset(ix);
	fbasis.push_back(wave.dimensionedWave());
	}
	RhoDMatrix mapping=RhoDMatrix(PDT::Spin1Half,false);
	for(unsigned int ix=0;ix<2;++ix) {
	if(fbasis[0].s2()==SqrtEnergy()) {
	mapping(ix,0) = sbasis[ix].s3()/fbasis[0].s3();
	mapping(ix,1) = sbasis[ix].s2()/fbasis[1].s2();
	}
	else {
	mapping(ix,0) = sbasis[ix].s2()/fbasis[0].s2();
	mapping(ix,1) = sbasis[ix].s3()/fbasis[1].s3();
	}
	}
	return mapping;
	}

	FermionSpinPtr createFermionSpinInfo(ShowerParticle & particle,
	const Lorentz5Momentum & porig,
	const LorentzRotation & rot,
	Helicity::Direction dir) {
	// calculate the splitting basis for the branching
	// and rotate back to construct the basis states
	LorentzRotation rinv = rot.inverse();
	SpinorWaveFunction wave;
	if(particle.id()>0)
	wave=SpinorWaveFunction(porig,particle.dataPtr(),incoming);
	else
	wave=SpinorWaveFunction(porig,particle.dataPtr(),outgoing);
	FermionSpinPtr fspin = new_ptr(FermionSpinInfo(particle.momentum(),dir==outgoing));
	for(unsigned int ix=0;ix<2;++ix) {
	wave.reset(ix);
	LorentzSpinor<SqrtEnergy> basis = wave.dimensionedWave();
	basis.transform(rinv);
	fspin->setBasisState(ix,basis);
	fspin->setDecayState(ix,basis);
	}
	particle.spinInfo(fspin);
	return fspin;
	}

	VectorSpinPtr createVectorSpinInfo(ShowerParticle & particle,
	const Lorentz5Momentum & porig,
	const LorentzRotation & rot,
	Helicity::Direction dir) {
	// calculate the splitting basis for the branching
	// and rotate back to construct the basis states
	LorentzRotation rinv = rot.inverse();
	bool massless(particle.id()==ParticleID::g\|\|particle.id()==ParticleID::gamma);
	VectorWaveFunction wave(porig,particle.dataPtr(),dir);
	VectorSpinPtr vspin = new_ptr(VectorSpinInfo(particle.momentum(),dir==outgoing));
	for(unsigned int ix=0;ix<3;++ix) {
	LorentzPolarizationVector basis;
	if(massless&&ix==1) {
	basis = LorentzPolarizationVector();
	}
	else {
	wave.reset(ix);
	basis = wave.wave();
	}
	basis *= rinv;
	vspin->setBasisState(ix,basis);
	vspin->setDecayState(ix,basis);
	}
	particle.spinInfo(vspin);
	vspin-> DMatrix() = RhoDMatrix(PDT::Spin1);
	vspin->rhoMatrix() = RhoDMatrix(PDT::Spin1);
	if(massless) {
	vspin-> DMatrix()(0,0) = 0.5;
	vspin->rhoMatrix()(0,0) = 0.5;
	vspin-> DMatrix()(2,2) = 0.5;
	vspin->rhoMatrix()(2,2) = 0.5;
	}
	return vspin;
	}
	}

	bool SudakovFormFactor::getMapping(SpinPtr & output, RhoDMatrix & mapping,
	ShowerParticle & particle,ShoKinPtr showerkin) {
	// if the particle is not from the hard process
	if(!particle.perturbative()) {
	// mapping is the identity
	output=particle.spinInfo();
	mapping=RhoDMatrix(particle.dataPtr()->iSpin());
	if(output) {
	return false;
	}
	else {
	assert(particle.parents().size()==1);
	Lorentz5Momentum porig;
	LorentzRotation rot = boostToShower(showerkin->getBasis(),showerkin->frame(),porig);
	if(particle.dataPtr()->iSpin()==PDT::Spin0) {
	assert(false);
	}
	else if(particle.dataPtr()->iSpin()==PDT::Spin1Half) {
	output = createFermionSpinInfo(particle,porig,rot,outgoing);
	}
	else if(particle.dataPtr()->iSpin()==PDT::Spin1) {
	output = createVectorSpinInfo(particle,porig,rot,outgoing);
	}
	else {
	assert(false);
	}
	return false;
	}
	}
	// if particle is final-state and is from the hard process
	else if(particle.isFinalState()) {
	assert(particle.perturbative()==1 \|\| particle.perturbative()==2);
	// get transform to shower frame
	Lorentz5Momentum porig;
	LorentzRotation rot = boostToShower(showerkin->getBasis(),showerkin->frame(),porig);
	// the rest depends on the spin of the particle
	PDT::Spin spin(particle.dataPtr()->iSpin());
	mapping=RhoDMatrix(spin,false);
	// do the spin dependent bit
	if(spin==PDT::Spin0) {
	- ScalarWaveFunction::constructSpinInfo(&particle,outgoing,true);
	+ ScalarSpinPtr sspin=dynamic_ptr_cast<ScalarSpinPtr>(particle.spinInfo());
	+ if(!sspin) {
	+ ScalarWaveFunction::constructSpinInfo(&particle,outgoing,true);
	+ }
	output=particle.spinInfo();
	return false;
	}
	else if(spin==PDT::Spin1Half) {
	FermionSpinPtr fspin=dynamic_ptr_cast<FermionSpinPtr>(particle.spinInfo());
	// spin info exists get information from it
	if(fspin) {
	output=fspin;
	mapping = fermionMapping(particle,porig,fspin,rot);
	return true;
	}
	// spin info does not exist create it
	else {
	output = createFermionSpinInfo(particle,porig,rot,outgoing);
	return false;
	}
	}
	else if(spin==PDT::Spin1) {
	VectorSpinPtr vspin=dynamic_ptr_cast<VectorSpinPtr>(particle.spinInfo());
	// spin info exists get information from it
	if(vspin) {
	output=vspin;
	mapping = bosonMapping(particle,porig,vspin,rot);
	return true;
	}
	else {
	output = createVectorSpinInfo(particle,porig,rot,outgoing);
	return false;
	}
	}
	// not scalar/fermion/vector
	else
	assert(false);
	}
	// incoming to hard process
	else if(particle.perturbative()==1 && !particle.isFinalState()) {
	// get the basis vectors
	// get transform to shower frame
	Lorentz5Momentum porig;
	LorentzRotation rot = boostToShower(showerkin->getBasis(),showerkin->frame(),porig);
	porig *= particle.x();
	// the rest depends on the spin of the particle
	PDT::Spin spin(particle.dataPtr()->iSpin());
	mapping=RhoDMatrix(spin);
	// do the spin dependent bit
	if(spin==PDT::Spin0) {
	cerr << "testing spin 0 not yet implemented " << endl;
	assert(false);
	}
	// spin-1/2
	else if(spin==PDT::Spin1Half) {
	FermionSpinPtr fspin=dynamic_ptr_cast<FermionSpinPtr>(particle.spinInfo());
	// spin info exists get information from it
	if(fspin) {
	output=fspin;
	mapping = fermionMapping(particle,porig,fspin,rot);
	return true;
	}
	// spin info does not exist create it
	else {
	output = createFermionSpinInfo(particle,porig,rot,incoming);
	return false;
	}
	}
	// spin-1
	else if(spin==PDT::Spin1) {
	VectorSpinPtr vspin=dynamic_ptr_cast<VectorSpinPtr>(particle.spinInfo());
	// spinInfo exists map it
	if(vspin) {
	output=vspin;
	mapping = bosonMapping(particle,porig,vspin,rot);
	return true;
	}
	// create the spininfo
	else {
	output = createVectorSpinInfo(particle,porig,rot,incoming);
	return false;
	}
	}
	assert(false);
	}
	// incoming to decay
	else if(particle.perturbative() == 2 && !particle.isFinalState()) {
	// get the basis vectors
	Lorentz5Momentum porig;
	LorentzRotation rot=boostToShower(showerkin->getBasis(),
	showerkin->frame(),porig);
	// the rest depends on the spin of the particle
	PDT::Spin spin(particle.dataPtr()->iSpin());
	mapping=RhoDMatrix(spin);
	// do the spin dependent bit
	if(spin==PDT::Spin0) {
	cerr << "testing spin 0 not yet implemented " << endl;
	assert(false);
	}
	// spin-1/2
	else if(spin==PDT::Spin1Half) {
	// FermionSpinPtr fspin=dynamic_ptr_cast<FermionSpinPtr>(particle.spinInfo());
	// // spin info exists get information from it
	// if(fspin) {
	// output=fspin;
	// mapping = fermionMapping(particle,porig,fspin,rot);
	// return true;
	// // spin info does not exist create it
	// else {
	// output = createFermionSpinInfo(particle,porig,rot,incoming);
	// return false;
	// }
	// }
	assert(false);
	}
	// // spin-1
	// else if(spin==PDT::Spin1) {
	// VectorSpinPtr vspin=dynamic_ptr_cast<VectorSpinPtr>(particle.spinInfo());
	// // spinInfo exists map it
	// if(vspin) {
	// output=vspin;
	// mapping = bosonMapping(particle,porig,vspin,rot);
	// return true;
	// }
	// // create the spininfo
	// else {
	// output = createVectorSpinInfo(particle,porig,rot,incoming);
	// return false;
	// }
	// }
	// assert(false);
	assert(false);
	}
	else
	assert(false);
	return true;
	}

	void SudakovFormFactor::removeSplitting(const IdList & in) {
	for(vector<IdList>::iterator it=particles_.begin();
	it!=particles_.end();++it) {
	if(it->size()==in.size()) {
	bool match=true;
	for(unsigned int iy=0;iy<in.size();++iy) {
	if((*it)[iy]!=in[iy]) {
	match=false;
	break;
	}
	}
	if(match) {
	vector<IdList>::iterator itemp=it;
	--itemp;
	particles_.erase(it);
	it = itemp;
	}
	}
	}
	}

	Energy2 SudakovFormFactor::guesst(Energy2 t1,unsigned int iopt,
	const IdList &ids,
	double enhance,bool ident) const {
	unsigned int pdfopt = iopt!=1 ? 0 : pdffactor_;
	double c =
	1./((splittingFn_->integOverP(zlimits_.second,ids,pdfopt) -
	splittingFn_->integOverP(zlimits_.first ,ids,pdfopt))*
	alpha_->overestimateValue()/Constants::twopi*enhance);
	assert(iopt<=2);
	if(iopt==1) {
	c/=pdfmax_;
	if(ident) c*=0.5;
	}
	else if(iopt==2) c*=-1.;
	if(splittingFn_->interactionOrder()==1) {
	double r = UseRandom::rnd();
	if(iopt!=2 \|\| c*log(r)<log(Constants::MaxEnergy2/t1)) {
	return t1*pow(r,c);
	}
	else
	return Constants::MaxEnergy2;
	}
	else {
	assert(false && "Units are dubious here.");
	int nm(splittingFn()->interactionOrder()-1);
	c/=Math::powi(alpha_->overestimateValue()/Constants::twopi,nm);
	return t1 / pow (1. - nmclog(UseRandom::rnd())
	* Math::powi(t1*UnitRemoval::InvE2,nm)
	,1./double(nm));
	}
	}

	double SudakovFormFactor::guessz (unsigned int iopt, const IdList &ids) const {
	unsigned int pdfopt = iopt!=1 ? 0 : pdffactor_;
	double lower = splittingFn_->integOverP(zlimits_.first,ids,pdfopt);
	return splittingFn_->invIntegOverP
	(lower + UseRandom::rnd()*(splittingFn_->integOverP(zlimits_.second,ids,pdfopt) -
	lower),ids,pdfopt);
	}

	void SudakovFormFactor::doinit() {
	Interfaced::doinit();
	pT2min_ = cutOffOption()==2 ? sqr(pTmin_) : ZERO;
	}

	const vector<Energy> & SudakovFormFactor::virtualMasses(const IdList & ids) {
	static vector<Energy> output;
	output.clear();
	if(cutOffOption() == 0) {
	for(unsigned int ix=0;ix<ids.size();++ix)
	output.push_back(getParticleData(ids[ix])->mass());
	Energy kinCutoff=
	kinematicCutOff(kinScale(),*std::max_element(output.begin(),output.end()));
	for(unsigned int ix=0;ix<output.size();++ix)
	output[ix]=max(kinCutoff,output[ix]);
	}
	else if(cutOffOption() == 1) {
	for(unsigned int ix=0;ix<ids.size();++ix) {
	output.push_back(getParticleData(ids[ix])->mass());
	output.back() += ids[ix]==ParticleID::g ? vgCut() : vqCut();
	}
	}
	else if(cutOffOption() == 2) {
	for(unsigned int ix=0;ix<ids.size();++ix)
	output.push_back(getParticleData(ids[ix])->mass());
	}
	else {
	throw Exception() << "Unknown option for the cut-off"
	<< " in SudakovFormFactor::virtualMasses()"
	<< Exception::runerror;
	}
	return output;
	}
	diff --git a/Shower/Default/Decay_QTildeShowerKinematics1to2.cc b/Shower/Default/Decay_QTildeShowerKinematics1to2.cc
	--- a/Shower/Default/Decay_QTildeShowerKinematics1to2.cc
	+++ b/Shower/Default/Decay_QTildeShowerKinematics1to2.cc
	@@ -1,100 +1,113 @@
	// -- C++ --
	//
	// Decay_QTildeShowerKinematics1to2.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 Decay_QTildeShowerKinematics1to2 class.
	//

	#include "Decay_QTildeShowerKinematics1to2.h"
	#include "ThePEG/PDT/EnumParticles.h"
	#include "Herwig++/Shower/SplittingFunctions/SplittingFunction.h"
	#include "Herwig++/Shower/Base/ShowerParticle.h"
	#include <cassert>
	+#include "Herwig++/Shower/ShowerHandler.h"
	+#include "Herwig++/Shower/Base/Evolver.h"
	+#include "Herwig++/Shower/Base/PartnerFinder.h"
	+#include "Herwig++/Shower/Base/ShowerModel.h"
	+#include "Herwig++/Shower/Base/KinematicsReconstructor.h"
	+#include "Herwig++/Shower/Base/ShowerVertex.h"

	using namespace Herwig;

	void Decay_QTildeShowerKinematics1to2::
	updateChildren(const tShowerParticlePtr parent,
	const ShowerParticleVector & children,
	ShowerPartnerType::Type partnerType) const {
	assert(children.size() == 2);
	// calculate the scales
	splittingFn()->evaluateDecayScales(partnerType,scale(),z(),parent,
	children[0],children[1]);
	// determine alphas of children according to interpretation of z
	const ShowerParticle::Parameters & params = parent->showerParameters();
	ShowerParticle::Parameters & child0 = children[0]->showerParameters();
	ShowerParticle::Parameters & child1 = children[1]->showerParameters();
	child0.alpha = z() * params.alpha;
	child1.alpha = (1.-z()) * params.alpha;

	child0.ptx = pT() * cos(phi()) + z()* params.ptx;
	child0.pty = pT() * sin(phi()) + z()* params.pty;
	child0.pt = sqrt( sqr(child0.ptx) + sqr(child0.pty) );

	child1.ptx = -pT() * cos(phi()) + (1.-z()) * params.ptx;
	child1.pty = -pT() * sin(phi()) + (1.-z()) * params.pty;
	child1.pt = sqrt( sqr(child1.ptx) + sqr(child1.pty) );

	// set up the colour connections
	splittingFn()->colourConnection(parent,children[0],children[1],partnerType,false);
	// make the products children of the parent
	parent->addChild(children[0]);
	parent->addChild(children[1]);
	+ if(! ShowerHandler::currentHandler()->evolver()->correlations()) return;
	+ SpinPtr pspin(parent->spinInfo());
	+ // set the momenta of the children
	+ ShowerParticleVector::const_iterator pit;
	+ for(pit=children.begin();pit!=children.end();++pit) {
	+ setMomentum(*pit,true);
	+ }
	}

	void Decay_QTildeShowerKinematics1to2::
	reconstructParent( const tShowerParticlePtr, const ParticleVector &) const {
	throw Exception() << "Decay_QTildeShowerKinematics1to2::updateParent not implemented"
	<< Exception::abortnow;
	}

	void Decay_QTildeShowerKinematics1to2::
	reconstructLast(const tShowerParticlePtr last, Energy mass) const {
	// set beta component and consequently all missing data from that,
	// using the nominal (i.e. PDT) mass.
	Energy theMass = mass>ZERO ? mass : last->data().constituentMass();
	last->showerParameters().beta=
	(sqr(theMass) + sqr(last->showerParameters().pt)
	- sqr( last->showerParameters().alpha )*pVector().m2())
	/ ( 2.last->showerParameters().alphap_dot_n() );
	// set that new momentum
	last->set5Momentum( sudakov2Momentum( last->showerParameters().alpha,
	last->showerParameters().beta,
	last->showerParameters().ptx,
	last->showerParameters().pty) );
	}

	void Decay_QTildeShowerKinematics1to2::initialize(ShowerParticle & particle,PPtr) {
	Lorentz5Momentum p, n, ppartner, pcm;
	Frame frame;
	assert(particle.perturbative()!=1);
	// this is for the initial decaying particle
	if(particle.perturbative()==2) {
	p = particle.momentum();
	ShowerParticlePtr partner=particle.partner();
	Lorentz5Momentum ppartner(partner->momentum());
	// reomved to make inverse recon work properly
	//if(partner->thePEGBase()) ppartner=partner->thePEGBase()->momentum();
	pcm=ppartner;
	Boost boost(p.findBoostToCM());
	pcm.boost(boost);
	n = Lorentz5Momentum( ZERO,0.5p.mass()pcm.vect().unit());
	n.boost( -boost);
	frame = Rest;
	}
	else {
	tShoKinPtr kin=dynamic_ptr_cast<ShowerParticlePtr>(particle.parents()[0])
	->showerKinematics();
	p = kin->getBasis()[0];
	n = kin->getBasis()[1];
	frame = kin->frame();
	}
	setBasis(p,n,frame);
	}
	diff --git a/Shower/Default/FS_QTildeShowerKinematics1to2.cc b/Shower/Default/FS_QTildeShowerKinematics1to2.cc
	--- a/Shower/Default/FS_QTildeShowerKinematics1to2.cc
	+++ b/Shower/Default/FS_QTildeShowerKinematics1to2.cc
	@@ -1,213 +1,217 @@
	// -- C++ --
	//
	// FS_QTildeShowerKinematics1to2.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 FS_QTildeShowerKinematics1to2 class.
	//

	#include "FS_QTildeShowerKinematics1to2.h"
	#include "ThePEG/PDT/EnumParticles.h"
	#include "Herwig++/Shower/SplittingFunctions/SplittingFunction.h"
	#include "Herwig++/Shower/Base/ShowerParticle.h"
	#include "ThePEG/Utilities/Debug.h"
	#include "Herwig++/Shower/ShowerHandler.h"
	#include "Herwig++/Shower/Base/Evolver.h"
	#include "Herwig++/Shower/Base/PartnerFinder.h"
	#include "Herwig++/Shower/Base/ShowerModel.h"
	#include "Herwig++/Shower/Base/KinematicsReconstructor.h"
	#include "Herwig++/Shower/Base/ShowerVertex.h"

	using namespace Herwig;

	void FS_QTildeShowerKinematics1to2::
	updateParameters(tShowerParticlePtr theParent,
	tShowerParticlePtr theChild0,
	tShowerParticlePtr theChild1,
	bool setAlpha) const {
	const ShowerParticle::Parameters & parent = theParent->showerParameters();
	ShowerParticle::Parameters & child0 = theChild0->showerParameters();
	ShowerParticle::Parameters & child1 = theChild1->showerParameters();
	// determine alphas of children according to interpretation of z
	if ( setAlpha ) {
	child0.alpha = z() * parent.alpha;
	child1.alpha = (1.-z()) * parent.alpha;
	}
	// set the values
	double cphi = cos(phi());
	double sphi = sin(phi());

	child0.ptx = pT() * cphi + z() * parent.ptx;
	child0.pty = pT() * sphi + z() * parent.pty;
	child0.pt = sqrt( sqr(child0.ptx) + sqr(child0.pty) );

	child1.ptx = -pT() * cphi + (1.-z())* parent.ptx;
	child1.pty = -pT() * sphi + (1.-z())* parent.pty;
	child1.pt = sqrt( sqr(child1.ptx) + sqr(child1.pty) );

	}

	void FS_QTildeShowerKinematics1to2::
	updateChildren(const tShowerParticlePtr parent,
	const ShowerParticleVector & children,
	ShowerPartnerType::Type partnerType) const {
	assert(children.size()==2);
	// calculate the scales
	splittingFn()->evaluateFinalStateScales(partnerType,scale(),z(),parent,
	children[0],children[1]);
	// update the parameters
	updateParameters(parent, children[0], children[1], true);
	// set up the colour connections
	splittingFn()->colourConnection(parent,children[0],children[1],partnerType,false);
	// make the products children of the parent
	parent->addChild(children[0]);
	parent->addChild(children[1]);
	// sort out the helicity stuff
	if(! ShowerHandler::currentHandler()->evolver()->correlations()) return;
	SpinPtr pspin(parent->spinInfo());
	- if(!pspin) return;
	+ // set the momenta of the children
	+ ShowerParticleVector::const_iterator pit;
	+ for(pit=children.begin();pit!=children.end();++pit) {
	+ setMomentum(*pit,true);
	+ }
	+ if(!pspin \|\| !ShowerHandler::currentHandler()->evolver()->spinCorrelations() ) return;
	Energy2 t = sqr(scale())z()(1.-z());
	IdList ids;
	ids.push_back(parent->id());
	ids.push_back(children[0]->id());
	ids.push_back(children[1]->id());
	// create the vertex
	SVertexPtr vertex(new_ptr(ShowerVertex()));
	// set the matrix element
	vertex->ME(splittingFn()->matrixElement(z(),t,ids,phi()));
	// set the incoming particle for the vertex
	parent->spinInfo()->decayVertex(vertex);
	- ShowerParticleVector::const_iterator pit;
	for(pit=children.begin();pit!=children.end();++pit) {
	// construct the spin info for the children
	constructSpinInfo(*pit,true);
	// connect the spinInfo object to the vertex
	(*pit)->spinInfo()->productionVertex(vertex);
	}
	}

	void FS_QTildeShowerKinematics1to2::
	reconstructParent(const tShowerParticlePtr parent,
	const ParticleVector & children ) const {
	assert(children.size() == 2);
	ShowerParticlePtr c1 = dynamic_ptr_cast<ShowerParticlePtr>(children[0]);
	ShowerParticlePtr c2 = dynamic_ptr_cast<ShowerParticlePtr>(children[1]);
	parent->showerParameters().beta=
	c1->showerParameters().beta + c2->showerParameters().beta;
	parent->set5Momentum( c1->momentum() + c2->momentum() );
	}

	void FS_QTildeShowerKinematics1to2::reconstructLast(const tShowerParticlePtr theLast,
	Energy mass) const {
	// set beta component and consequently all missing data from that,
	// using the nominal (i.e. PDT) mass.
	Energy theMass = mass > ZERO ? mass : theLast->data().constituentMass();
	ShowerParticle::Parameters & last = theLast->showerParameters();
	last.beta = ( sqr(theMass) + sqr(last.pt) - sqr(last.alpha) * pVector().m2() )
	/ ( 2. * last.alpha * p_dot_n() );
	// set that new momentum
	theLast->set5Momentum(sudakov2Momentum( last.alpha, last.beta,
	last.ptx, last.pty) );
	}

	void FS_QTildeShowerKinematics1to2::initialize(ShowerParticle & particle,PPtr) {
	// set the basis vectors
	Lorentz5Momentum p,n;
	Frame frame;
	if(particle.perturbative()!=0) {
	// find the partner and its momentum
	ShowerParticlePtr partner=particle.partner();
	Lorentz5Momentum ppartner(partner->momentum());
	// momentum of the emitting particle
	p = particle.momentum();
	Lorentz5Momentum pcm;
	// if the partner is a final-state particle then the reference
	// vector is along the partner in the rest frame of the pair
	if(partner->isFinalState()) {
	Boost boost=(p + ppartner).findBoostToCM();
	pcm = ppartner;
	pcm.boost(boost);
	n = Lorentz5Momentum(ZERO,pcm.vect());
	n.boost( -boost);
	}
	else if(!partner->isFinalState()) {
	// if the partner is an initial-state particle then the reference
	// vector is along the partner which should be massless
	if(particle.perturbative()==1)
	{n = Lorentz5Momentum(ZERO,ppartner.vect());}
	// if the partner is an initial-state decaying particle then the reference
	// vector is along the backwards direction in rest frame of decaying particle
	else {
	Boost boost=ppartner.findBoostToCM();
	pcm = p;
	pcm.boost(boost);
	n = Lorentz5Momentum( ZERO, -pcm.vect());
	n.boost( -boost);
	}
	}
	frame = BackToBack;
	}
	else if(particle.initiatesTLS()) {
	tShoKinPtr kin=dynamic_ptr_cast<ShowerParticlePtr>
	(particle.parents()[0]->children()[0])->showerKinematics();
	p = kin->getBasis()[0];
	n = kin->getBasis()[1];
	frame = kin->frame();
	}
	else {
	tShoKinPtr kin=dynamic_ptr_cast<ShowerParticlePtr>(particle.parents()[0])
	->showerKinematics();
	p = kin->getBasis()[0];
	n = kin->getBasis()[1];
	frame = kin->frame();
	}
	// set the basis vectors
	setBasis(p,n,frame);
	}

	void FS_QTildeShowerKinematics1to2::updateParent(const tShowerParticlePtr parent,
	const ShowerParticleVector & children,
	ShowerPartnerType::Type) const {
	IdList ids(3);
	ids[0] = parent->id();
	ids[1] = children[0]->id();
	ids[2] = children[1]->id();
	const vector<Energy> & virtualMasses = SudakovFormFactor()->virtualMasses(ids);
	if(children[0]->children().empty()) children[0]->virtualMass(virtualMasses[1]);
	if(children[1]->children().empty()) children[1]->virtualMass(virtualMasses[2]);
	// compute the new pT of the branching
	Energy2 pt2=sqr(z()(1.-z()))sqr(scale())
	- sqr(children[0]->virtualMass())*(1.-z())
	- sqr(children[1]->virtualMass())* z() ;
	if(ids[0]!=ParticleID::g) pt2 += z()(1.-z())sqr(virtualMasses[0]);
	Energy2 q2 =
	sqr(children[0]->virtualMass())/z() +
	sqr(children[1]->virtualMass())/(1.-z()) +
	pt2/z()/(1.-z());
	if(pt2<ZERO) {
	parent->virtualMass(ZERO);
	}
	else {
	parent->virtualMass(sqrt(q2));
	pT(sqrt(pt2));
	}
	}

	void FS_QTildeShowerKinematics1to2::
	resetChildren(const tShowerParticlePtr parent,
	const ShowerParticleVector & children) const {
	updateParameters(parent, children[0], children[1], false);
	for(unsigned int ix=0;ix<children.size();++ix) {
	if(children[ix]->children().empty()) continue;
	ShowerParticleVector newChildren;
	for(unsigned int iy=0;iy<children[ix]->children().size();++iy)
	newChildren.push_back(dynamic_ptr_cast<ShowerParticlePtr>
	(children[ix]->children()[iy]));
	children[ix]->showerKinematics()->resetChildren(children[ix],newChildren);
	}
	}
	diff --git a/Shower/Default/IS_QTildeShowerKinematics1to2.cc b/Shower/Default/IS_QTildeShowerKinematics1to2.cc
	--- a/Shower/Default/IS_QTildeShowerKinematics1to2.cc
	+++ b/Shower/Default/IS_QTildeShowerKinematics1to2.cc
	@@ -1,190 +1,192 @@
	// -- C++ --
	//
	// IS_QTildeShowerKinematics1to2.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 IS_QTildeShowerKinematics1to2 class.
	//

	#include "IS_QTildeShowerKinematics1to2.h"
	#include "ThePEG/PDT/EnumParticles.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "Herwig++/Shower/Base/ShowerParticle.h"
	#include "ThePEG/Utilities/Debug.h"
	#include "Herwig++/Shower/ShowerHandler.h"
	#include "Herwig++/Shower/Base/Evolver.h"
	#include "Herwig++/Shower/Base/PartnerFinder.h"
	#include "Herwig++/Shower/Base/ShowerModel.h"
	#include "Herwig++/Shower/Base/KinematicsReconstructor.h"
	#include "Herwig++/Shower/Base/ShowerVertex.h"
	#include <cassert>

	using namespace Herwig;

	void IS_QTildeShowerKinematics1to2::
	updateChildren( const tShowerParticlePtr theParent,
	const ShowerParticleVector & theChildren,
	ShowerPartnerType::Type) const {
	const ShowerParticle::Parameters & parent = theParent->showerParameters();
	ShowerParticle::Parameters & child0 = theChildren[0]->showerParameters();
	ShowerParticle::Parameters & child1 = theChildren[1]->showerParameters();
	double cphi = cos(phi());
	double sphi = sin(phi());

	child1.alpha = (1.-z()) * parent.alpha;

	child1.ptx = (1.-z()) * parent.ptx - cphi * pT();
	child1.pty = (1.-z()) * parent.pty - sphi * pT();
	child1.pt = sqrt( sqr(child1.ptx) + sqr(child1.pty) );
	// space-like child
	child0.alpha = parent.alpha - child1.alpha;
	child0.beta = parent.beta - child1.beta;
	child0.ptx = parent.ptx - child1.ptx;
	child0.pty = parent.pty - child1.pty;
	}


	void IS_QTildeShowerKinematics1to2::
	updateParent(const tShowerParticlePtr parent,
	const ShowerParticleVector & children,
	ShowerPartnerType::Type partnerType) const {
	// calculate the scales
	splittingFn()->evaluateInitialStateScales(partnerType,scale(),z(),parent,
	children[0],children[1]);
	// set proper colour connections
	splittingFn()->colourConnection(parent,children[0],children[1],
	partnerType,true);
	// set proper parent/child relationships
	parent->addChild(children[0]);
	parent->addChild(children[1]);
	parent->x(children[0]->x()/z());
	// sort out the helicity stuff
	if(! ShowerHandler::currentHandler()->evolver()->correlations()) return;
	+ // construct the spin info for parent and timelike child
	+ // temporary assignment of shower parameters to calculate correlations
	+ parent->showerParameters().alpha = parent->x();
	+ children[1]->showerParameters().alpha = (1.-z()) * parent->x();
	+ children[1]->showerParameters().ptx = - cos(phi()) * pT();
	+ children[1]->showerParameters().pty = - sin(phi()) * pT();
	+ children[1]->showerParameters().pt = pT();
	+ setMomentum(parent,false);
	+ setMomentum(children[1],true);
	SpinPtr pspin(children[0]->spinInfo());
	- if(!pspin) return;
	+ if(!pspin \|\| !ShowerHandler::currentHandler()->evolver()->spinCorrelations() ) return;
	// compute the matrix element for spin correlations
	IdList ids;
	ids.push_back(parent->id());
	ids.push_back(children[0]->id());
	ids.push_back(children[1]->id());
	Energy2 t = (1.-z())*sqr(scale())/z();
	// create the vertex
	SVertexPtr vertex(new_ptr(ShowerVertex()));
	// set the matrix element
	vertex->ME(splittingFn()->matrixElement(z(),t,ids,phi()));
	// set the incoming particle for the vertex
	// (in reality the first child as going backwards)
	pspin->decayVertex(vertex);
	- // construct the spin info for parent and timelike child
	- // temporary assignment of shower parameters to calculate correlations
	- parent->showerParameters().alpha = parent->x();
	- children[1]->showerParameters().alpha = (1.-z()) * parent->x();
	- children[1]->showerParameters().ptx = - cos(phi()) * pT();
	- children[1]->showerParameters().pty = - sin(phi()) * pT();
	- children[1]->showerParameters().pt = pT();
	// construct the spin infos
	constructSpinInfo(parent,false);
	constructSpinInfo(children[1],true);
	// connect the spinInfo objects to the vertex
	parent ->spinInfo()->productionVertex(vertex);
	children[1]->spinInfo()->productionVertex(vertex);
	}

	void IS_QTildeShowerKinematics1to2::
	reconstructParent(const tShowerParticlePtr theParent,
	const ParticleVector & theChildren ) const {
	PPtr c1 = theChildren[0];
	ShowerParticlePtr c2 = dynamic_ptr_cast<ShowerParticlePtr>(theChildren[1]);
	ShowerParticle::Parameters & c2param = c2->showerParameters();
	// get shower variables from 1st child in order to keep notation
	// parent->(c1, c2) clean even though the splitting was initiated
	// from c1. The name updateParent is still referring to the
	// timelike branching though.
	// on-shell child
	c2param.beta = 0.5*( sqr(c2->data().constituentMass()) + sqr(c2param.pt) )
	/ ( c2param.alpha * p_dot_n() );
	c2->set5Momentum( sudakov2Momentum(c2param.alpha, c2param.beta,
	c2param.ptx , c2param.pty) );
	// spacelike child
	Lorentz5Momentum pc1(theParent->momentum() - c2->momentum());
	pc1.rescaleMass();
	c1->set5Momentum(pc1);
	}

	void IS_QTildeShowerKinematics1to2::
	updateLast( const tShowerParticlePtr theLast,Energy px,Energy py) const {
	if(theLast->isFinalState()) return;
	ShowerParticle::Parameters & last = theLast->showerParameters();
	Energy2 pt2 = sqr(px) + sqr(py);
	last.alpha = theLast->x();
	last.beta = 0.5 * pt2 / last.alpha / p_dot_n();
	last.ptx = ZERO;
	last.pty = ZERO;
	last.pt = ZERO;
	// momentum
	Lorentz5Momentum ntemp = Lorentz5Momentum(ZERO,-pVector().vect());
	double beta = 0.5 * pt2 / last.alpha / (pVector() * ntemp);
	Lorentz5Momentum plast =
	Lorentz5Momentum( (pVector().z()>ZERO ? px : -px), py, ZERO, ZERO)
	+ theLast->x() * pVector() + beta * ntemp;
	plast.rescaleMass();
	theLast->set5Momentum(plast);
	}

	void IS_QTildeShowerKinematics1to2::initialize(ShowerParticle & particle, PPtr parent) {
	// For the time being we are considering only 1->2 branching
	Lorentz5Momentum p, n, pthis, pcm;
	assert(particle.perturbative()!=2);
	Frame frame;
	if(particle.perturbative()==1) {
	// find the partner and its momentum
	ShowerParticlePtr partner=particle.partner();
	assert(partner);
	if(partner->isFinalState()) {
	Lorentz5Momentum pa = -particle.momentum()+partner->momentum();
	Lorentz5Momentum pb = particle.momentum();
	Energy scale=parent->momentum().t();
	Lorentz5Momentum pbasis(ZERO,parent->momentum().vect().unit()*scale);
	Axis axis(pa.vect().unit());
	LorentzRotation rot;
	double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
	if(axis.perp2()>1e-20) {
	rot.setRotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
	rot.rotateX(Constants::pi);
	}
	if(abs(1.-pa.e()/pa.vect().mag())>1e-6) rot.boostZ( pa.e()/pa.vect().mag());
	pb *= rot;


	if(pb.perp2()/GeV2>1e-20) {
	Boost trans = -1./pb.e()*pb.vect();
	trans.setZ(0.);
	rot.boost(trans);
	}
	pbasis *=rot;
	rot.invert();
	n = rot*Lorentz5Momentum(ZERO,-pbasis.vect());
	p = rot*Lorentz5Momentum(ZERO, pbasis.vect());
	}
	else {
	pcm = parent->momentum();
	p = Lorentz5Momentum(ZERO, pcm.vect());
	n = Lorentz5Momentum(ZERO, -pcm.vect());
	}
	frame = BackToBack;
	}
	else {
	p = dynamic_ptr_cast<ShowerParticlePtr>(particle.children()[0])
	->showerKinematics()->getBasis()[0];
	n = dynamic_ptr_cast<ShowerParticlePtr>(particle.children()[0])
	->showerKinematics()->getBasis()[1];
	frame = dynamic_ptr_cast<ShowerParticlePtr>(particle.children()[0])
	->showerKinematics()->frame();
	}
	setBasis(p,n,frame);
	}
	diff --git a/Shower/Default/QTildeShowerKinematics1to2.cc b/Shower/Default/QTildeShowerKinematics1to2.cc
	--- a/Shower/Default/QTildeShowerKinematics1to2.cc
	+++ b/Shower/Default/QTildeShowerKinematics1to2.cc
	@@ -1,158 +1,148 @@
	// -- C++ --
	//
	// QTildeShowerKinematics1to2.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 QTildeShowerKinematics1to2 class.
	//

	#include "QTildeShowerKinematics1to2.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "Herwig++/Shower/Base/ShowerParticle.h"
	#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
	#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
	#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
	#include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h"
	#include "ThePEG/Helicity/LorentzSpinorBar.h"

	using namespace Herwig;
	using namespace ThePEG::Helicity;

	vector<Lorentz5Momentum> QTildeShowerKinematics1to2::getBasis() const {
	vector<Lorentz5Momentum> dum;
	dum.push_back( _pVector );
	dum.push_back( _nVector );
	return dum;
	}

	void QTildeShowerKinematics1to2::setBasis(const Lorentz5Momentum &p,
	const Lorentz5Momentum & n,
	Frame inframe) {
	_pVector=p;
	_nVector=n;
	frame(inframe);
	}

	Lorentz5Momentum QTildeShowerKinematics1to2::
	sudakov2Momentum(double alpha, double beta, Energy px, Energy py) const {
	if(isnan(beta)\|\|isinf(beta))
	throw Exception() << "beta infinite in "
	<< "QTildeShowerKinematics1to2::sudakov2Momentum()"
	<< Exception::eventerror;
	- Lorentz5Momentum dq;
	+ Boost beta_bb;
	if(frame()==BackToBack) {
	- const Boost beta_bb = -(_pVector + _nVector).boostVector();
	- Lorentz5Momentum p_bb = _pVector;
	- Lorentz5Momentum n_bb = _nVector;
	- p_bb.boost( beta_bb );
	- n_bb.boost( beta_bb );
	- // set first in b2b frame along z-axis (assuming that p and n are
	- // b2b as checked above)
	- dq=Lorentz5Momentum(ZERO, ZERO, (alpha - beta)*p_bb.vect().mag(),
	- alphap_bb.t() + betan_bb.t());
	- // add transverse components
	- dq.setX(px);
	- dq.setY(py);
	- // rotate to have z-axis parallel to p
	- // this rotation changed by PR to a different rotation with the same effect
	- // but different azimuthal angle to make implementing spin correlations easier
	- // dq.rotateUz( unitVector(p_bb.vect()) );
	- Axis axis(p_bb.vect().unit());
	- LorentzRotation rot;
	- if(axis.perp2()>0.) {
	- double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
	- rot.rotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
	- }
	- else if(axis.z()<0.) {
	- rot.rotate(Constants::pi,Axis(1.,0.,0.));
	- }
	- dq.transform(rot);
	- // boost back
	- dq.boost( -beta_bb );
	- dq.rescaleMass();
	- // return the momentum
	+ beta_bb = -(_pVector + _nVector).boostVector();
	}
	else if(frame()==Rest) {
	- const Boost beta_bb = -pVector().boostVector();
	- Lorentz5Momentum p_bb = pVector();
	- Lorentz5Momentum n_bb = nVector();
	- p_bb.boost( beta_bb );
	- n_bb.boost( beta_bb );
	- // set first in b2b frame along z-axis (assuming that p and n are
	- // b2b as checked above)
	- dq=Lorentz5Momentum (ZERO, ZERO, 0.5betapVector().mass(),
	- alphapVector().mass() + 0.5beta*pVector().mass());
	- // add transverse components
	- dq.setX(px);
	- dq.setY(py);
	- // changed to be same as other case
	-// dq.rotateUz( unitVector(n_bb.vect()) );
	- Axis axis(n_bb.vect().unit());
	- LorentzRotation rot;
	- if(axis.perp2()>0.) {
	- double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
	- rot.rotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
	- }
	- else if(axis.z()<0.) {
	- rot.rotate(Constants::pi,Axis(1.,0.,0.));
	- }
	- dq.transform(rot);
	- // boost back
	- dq.boost( -beta_bb );
	- dq.rescaleMass();
	+ beta_bb = -pVector().boostVector();
	}
	else
	assert(false);
	+ Lorentz5Momentum p_bb = pVector();
	+ Lorentz5Momentum n_bb = nVector();
	+ p_bb.boost( beta_bb );
	+ n_bb.boost( beta_bb );
	+ // momentum without transverse components
	+ Lorentz5Momentum dq;
	+ if(frame()==BackToBack) {
	+ dq=Lorentz5Momentum(ZERO, ZERO, (alpha - beta)*p_bb.vect().mag(),
	+ alphap_bb.t() + betan_bb.t());
	+ }
	+ else if(frame()==Rest) {
	+ dq=Lorentz5Momentum (ZERO, ZERO, 0.5betapVector().mass(),
	+ alphapVector().mass() + 0.5beta*pVector().mass());
	+ }
	+ else
	+ assert(false);
	+ // add transverse components
	+ dq.setX(px);
	+ dq.setY(py);
	+ // rotate to have z-axis parallel to p/n
	+ Axis axis;
	+ if(frame()==BackToBack) {
	+ axis = p_bb.vect().unit();
	+ }
	+ else if(frame()==Rest) {
	+ axis = n_bb.vect().unit();
	+ }
	+ else
	+ assert(false);
	+ LorentzRotation rot;
	+ if(axis.perp2()>0.) {
	+ double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
	+ rot.rotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
	+ }
	+ else if(axis.z()<0.) {
	+ rot.rotate(Constants::pi,Axis(1.,0.,0.));
	+ }
	+ dq.transform(rot);
	+ // boost back
	+ dq.boost( -beta_bb );
	+ dq.rescaleMass();
	return dq;
	}

	-void QTildeShowerKinematics1to2::constructSpinInfo(tShowerParticlePtr particle,
	- bool timeLike) const {
	+void QTildeShowerKinematics1to2::setMomentum(tShowerParticlePtr particle,
	+ bool timeLike) const {
	Energy mass = particle->data().mass();
	// calculate the momentum of the assuming on-shell
	Energy2 pt2 = sqr(particle->showerParameters().pt);
	double alpha = timeLike ? particle->showerParameters().alpha : particle->x();
	double beta = 0.5(sqr(mass) + pt2 - sqr(alpha)pVector().m2())/(alpha*p_dot_n());
	Lorentz5Momentum porig=sudakov2Momentum(alpha,beta,
	particle->showerParameters().ptx,
	particle->showerParameters().pty);
	porig.setMass(mass);
	+ particle->set5Momentum(porig);
	+}
	+
	+void QTildeShowerKinematics1to2::constructSpinInfo(tShowerParticlePtr particle,
	+ bool timeLike) const {
	// now construct the required spininfo and calculate the basis states
	PDT::Spin spin(particle->dataPtr()->iSpin());
	- particle->set5Momentum(porig);
	if(spin==PDT::Spin0) {
	ScalarWaveFunction::constructSpinInfo(particle,outgoing,timeLike);
	}
	// calculate the basis states and construct the SpinInfo for a spin-1/2 particle
	else if(spin==PDT::Spin1Half) {
	// outgoing particle
	if(particle->id()>0) {
	vector<LorentzSpinorBar<SqrtEnergy> > stemp;
	SpinorBarWaveFunction::calculateWaveFunctions(stemp,particle,outgoing);
	SpinorBarWaveFunction::constructSpinInfo(stemp,particle,outgoing,timeLike);
	}
	// outgoing antiparticle
	else {
	vector<LorentzSpinor<SqrtEnergy> > stemp;
	SpinorWaveFunction::calculateWaveFunctions(stemp,particle,outgoing);
	SpinorWaveFunction::constructSpinInfo(stemp,particle,outgoing,timeLike);
	}
	}
	// calculate the basis states and construct the SpinInfo for a spin-1 particle
	else if(spin==PDT::Spin1) {
	bool massless(particle->id()==ParticleID::g\|\|particle->id()==ParticleID::gamma);
	vector<Helicity::LorentzPolarizationVector> vtemp;
	VectorWaveFunction::calculateWaveFunctions(vtemp,particle,outgoing,massless);
	VectorWaveFunction::constructSpinInfo(vtemp,particle,outgoing,timeLike,massless);
	}
	else {
	throw Exception() << "Spins higher than 1 are not yet implemented in "
	<< "FS_QtildaShowerKinematics1to2::constructVertex() "
	<< Exception::runerror;
	}
	}
	diff --git a/Shower/Default/QTildeShowerKinematics1to2.h b/Shower/Default/QTildeShowerKinematics1to2.h
	--- a/Shower/Default/QTildeShowerKinematics1to2.h
	+++ b/Shower/Default/QTildeShowerKinematics1to2.h
	@@ -1,110 +1,115 @@
	// -- C++ --
	//
	// QTildeShowerKinematics1to2.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_QTildeShowerKinematics1to2_H
	#define HERWIG_QTildeShowerKinematics1to2_H
	//
	// This is the declaration of the QTildeShowerKinematics1to2 class.
	//

	#include "Herwig++/Shower/Base/ShowerKinematics.h"
	#include "ThePEG/Vectors/Lorentz5Vector.h"
	#include "QTildeShowerKinematics1to2.fh"

	namespace Herwig {

	using namespace ThePEG;

	/** \ingroup Shower
	*
	* This abstract class describes the common features for initial and final
	* state radiation kinematics for \f$1\to2\f$ branchings and for
	* the choice of \f$\tilde{q}\f$ as evolution variable.
	*
	* @see ShowerKinematics
	* @see IS_QTildeShowerKinematics1to2
	* @see FS_QTildeShowerKinematics1to2
	* @see KinematicsReconstructor
	*/
	class QTildeShowerKinematics1to2: public ShowerKinematics {

	public:

	/**
	* Implementation of the virtual function returning a set of basis vectors, specific to
	* the type of evolution. This function will be used by the
	* ForwardShowerEvolver in order to access \f$p\f$
	* and \f$n\f$.
	*/
	virtual vector<Lorentz5Momentum> getBasis() const;

	/**
	* Access to the \f$p\f$ vector used to describe the kinematics.
	*/
	const Lorentz5Momentum & pVector() const {return _pVector;}

	/**
	* Access to the \f$n\f$ vector used to describe the kinematics.
	*/
	const Lorentz5Momentum & nVector() const {return _nVector;}

	/**
	* Dot product of thew basis vectors
	*/
	Energy2 p_dot_n() const {return _pVector*_nVector;}

	/**
	* Converts a Sudakov parametrization of a momentum w.r.t. the given
	* basis \f$p\f$ and \f$n\f$ into a 5 momentum.
	* @param alpha The \f$\alpha\f$ parameter of the Sudakov parameterisation
	* @param beta The \f$\beta\f$ parameter of the Sudakov parameterisation
	* @param px The \f$x\f$-component of the transverse momentum in the Sudakov
	* parameterisation
	* @param py The \f$x\f$-component of the transverse momentum in the Sudakov
	* parameterisation
	*/
	Lorentz5Momentum sudakov2Momentum(double alpha, double beta,
	Energy px, Energy py) const;

	protected:

	/**
	* Set the basis vectors
	*/
	void setBasis(const Lorentz5Momentum &p, const Lorentz5Momentum & n, Frame frame);

	/**
	+ * Set a preliminary momentum for the particle
	+ */
	+ void setMomentum(tShowerParticlePtr,bool timelike) const;
	+
	+ /**
	* Construct the spin info object for a shower particle
	*/
	void constructSpinInfo(tShowerParticlePtr,bool timelike) const;

	private:

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

	private:

	/**
	* The \f$p\f$ reference vector
	*/
	Lorentz5Momentum _pVector;

	/**
	* The \f$n\f$ reference vector
	*/
	Lorentz5Momentum _nVector;

	};

	}

	#endif /* HERWIG_QTildeShowerKinematics1to2_H */
	diff --git a/Shower/Default/QTildeSudakov.cc b/Shower/Default/QTildeSudakov.cc
	--- a/Shower/Default/QTildeSudakov.cc
	+++ b/Shower/Default/QTildeSudakov.cc
	@@ -1,585 +1,868 @@
	// -- C++ --
	//
	// QTildeSudakov.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 QTildeSudakov class.
	//

	#include "QTildeSudakov.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/PDT/ParticleData.h"
	#include "ThePEG/EventRecord/Event.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/PDT/EnumParticles.h"
	#include "Herwig++/Shower/Default/FS_QTildeShowerKinematics1to2.h"
	#include "Herwig++/Shower/Default/IS_QTildeShowerKinematics1to2.h"
	#include "Herwig++/Shower/Default/Decay_QTildeShowerKinematics1to2.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "Herwig++/Shower/Base/ShowerVertex.h"
	#include "Herwig++/Shower/Base/ShowerParticle.h"
	#include "Herwig++/Shower/ShowerHandler.h"
	#include "Herwig++/Shower/Base/Evolver.h"
	#include "Herwig++/Shower/Base/PartnerFinder.h"
	#include "Herwig++/Shower/Base/ShowerModel.h"
	#include "Herwig++/Shower/Base/KinematicsReconstructor.h"

	using namespace Herwig;

	DescribeNoPIOClass<QTildeSudakov,Herwig::SudakovFormFactor>
	describeQTildeSudakov ("Herwig::QTildeSudakov","HwShower.so");

	void QTildeSudakov::Init() {

	static ClassDocumentation<QTildeSudakov> documentation
	("The QTildeSudakov class implements the Sudakov form factor for ordering it"
	" qtilde");
	}

	bool QTildeSudakov::guessTimeLike(Energy2 &t,Energy2 tmin,double enhance) {
	Energy2 told = t;
	// calculate limits on z and if lower>upper return
	if(!computeTimeLikeLimits(t)) return false;
	// guess values of t and z
	t = guesst(told,0,ids_,enhance,ids_[1]==ids_[2]);
	z(guessz(0,ids_));
	// actual values for z-limits
	if(!computeTimeLikeLimits(t)) return false;
	if(t<tmin) {
	t=-1.0*GeV2;
	return false;
	}
	else
	return true;
	}

	bool QTildeSudakov::guessSpaceLike(Energy2 &t, Energy2 tmin, const double x,
	double enhance) {
	Energy2 told = t;
	// calculate limits on z if lower>upper return
	if(!computeSpaceLikeLimits(t,x)) return false;
	// guess values of t and z
	t = guesst(told,1,ids_,enhance,ids_[1]==ids_[2]);
	z(guessz(1,ids_));
	// actual values for z-limits
	if(!computeSpaceLikeLimits(t,x)) return false;
	if(t<tmin) {
	t=-1.0*GeV2;
	return false;
	}
	else
	return true;
	}

	bool QTildeSudakov::PSVeto(const Energy2 t) {
	// still inside PS, return true if outside
	// check vs overestimated limits
	if(z() < zLimits().first \|\| z() > zLimits().second) return true;
	// compute the pts
	Energy2 pt2=sqr(z()(1.-z()))t-masssquared_[1](1.-z())-masssquared_[2]z();
	if(ids_[0]!=ParticleID::g) pt2+=z()(1.-z())masssquared_[0];
	// if pt2<0 veto
	if(pt2<pT2min()) return true;
	// otherwise calculate pt and return
	pT(sqrt(pt2));
	return false;
	}

	ShoKinPtr QTildeSudakov::generateNextTimeBranching(const Energy startingScale,
	const IdList &ids,const bool cc,
	double enhance) {
	// First reset the internal kinematics variables that can
	// have been eventually set in the previous call to the method.
	q_ = ZERO;
	z(0.);
	phi(0.);
	// perform initialization
	Energy2 tmax(sqr(startingScale)),tmin;
	initialize(ids,tmin,cc);
	// check max > min
	if(tmax<=tmin) return ShoKinPtr();
	// calculate next value of t using veto algorithm
	Energy2 t(tmax);
	do {
	if(!guessTimeLike(t,tmin,enhance)) break;
	}
	while(PSVeto(t) \|\| SplittingFnVeto(z()(1.-z())t,ids,true) \|\|
	alphaSVeto(sqr(z()(1.-z()))t));
	q_ = t > ZERO ? Energy(sqrt(t)) : -1.*MeV;
	if(q_ < ZERO) return ShoKinPtr();
	// return the ShowerKinematics object
	return createFinalStateBranching(q_,z(),phi(),pT());
	}

	ShoKinPtr QTildeSudakov::
	generateNextSpaceBranching(const Energy startingQ,
	const IdList &ids,
	double x,bool cc,
	double enhance,
	Ptr<BeamParticleData>::transient_const_pointer beam) {
	// First reset the internal kinematics variables that can
	// have been eventually set in the previous call to the method.
	q_ = ZERO;
	z(0.);
	phi(0.);
	// perform the initialization
	Energy2 tmax(sqr(startingQ)),tmin;
	initialize(ids,tmin,cc);
	// check max > min
	if(tmax<=tmin) return ShoKinPtr();
	// extract the partons which are needed for the PDF veto
	// Different order, incoming parton is id = 1, outgoing are id=0,2
	tcPDPtr parton0 = getParticleData(ids[0]);
	tcPDPtr parton1 = getParticleData(ids[1]);
	if(cc) {
	if(parton0->CC()) parton0 = parton0->CC();
	if(parton1->CC()) parton1 = parton1->CC();
	}
	// calculate next value of t using veto algorithm
	Energy2 t(tmax),pt2(ZERO);
	do {
	if(!guessSpaceLike(t,tmin,x,enhance)) break;
	pt2=sqr(1.-z())t-z()masssquared_[2];
	}
	while(z() > zLimits().second \|\|
	SplittingFnVeto((1.-z())*t/z(),ids,true) \|\|
	alphaSVeto(sqr(1.-z())*t) \|\|
	PDFVeto(t,x,parton0,parton1,beam) \|\| pt2 < pT2min() );
	if(t > ZERO && zLimits().first < zLimits().second) q_ = sqrt(t);
	else return ShoKinPtr();
	pT(sqrt(pt2));
	// create the ShowerKinematics and return it
	return createInitialStateBranching(q_,z(),phi(),pT());
	}

	void QTildeSudakov::initialize(const IdList & ids, Energy2 & tmin,const bool cc) {
	ids_=ids;
	if(cc) {
	for(unsigned int ix=0;ix<ids.size();++ix) {
	if(getParticleData(ids[ix])->CC()) ids_[ix]*=-1;
	}
	}
	tmin = cutOffOption() != 2 ? ZERO : 4.*pT2min();
	masses_ = virtualMasses(ids);
	masssquared_.clear();
	for(unsigned int ix=0;ix<masses_.size();++ix) {
	masssquared_.push_back(sqr(masses_[ix]));
	if(ix>0) tmin=max(masssquared_[ix],tmin);
	}
	}

	ShoKinPtr QTildeSudakov::generateNextDecayBranching(const Energy startingScale,
	const Energy stoppingScale,
	const Energy minmass,
	const IdList &ids,
	const bool cc,
	double enhance) {
	// First reset the internal kinematics variables that can
	// have been eventually set in the previous call to this method.
	q_ = Constants::MaxEnergy;
	z(0.);
	phi(0.);
	// perform initialisation
	Energy2 tmax(sqr(stoppingScale)),tmin;
	initialize(ids,tmin,cc);
	tmin=sqr(startingScale);
	// check some branching possible
	if(tmax<=tmin) return ShoKinPtr();
	// perform the evolution
	Energy2 t(tmin),pt2(-MeV2);
	do {
	if(!guessDecay(t,tmax,minmass,enhance)) break;
	pt2 = sqr(1.-z())(t-masssquared_[0])-z()masssquared_[2];
	}
	while(SplittingFnVeto((1.-z())*t/z(),ids,true)\|\|
	alphaSVeto(sqr(1.-z())*t) \|\|
	pt2<pT2min() \|\|
	t*(1.-z())>masssquared_[0]-sqr(minmass));
	if(t > ZERO) {
	q_ = sqrt(t);
	pT(sqrt(pt2));
	}
	else return ShoKinPtr();
	phi(0.);
	// create the ShowerKinematics object
	return createDecayBranching(q_,z(),phi(),pT());
	}

	bool QTildeSudakov::guessDecay(Energy2 &t,Energy2 tmax, Energy minmass,
	double enhance) {
	// previous scale
	Energy2 told = t;
	// overestimated limits on z
	if(tmax<masssquared_[0]) {
	t=-1.0*GeV2;
	return false;
	}
	Energy2 tm2 = tmax-masssquared_[0];
	Energy tm = sqrt(tm2);
	pair<double,double> limits=make_pair(sqr(minmass/masses_[0]),
	1.-sqrt(masssquared_[2]+pT2min()+
	0.25*sqr(masssquared_[2])/tm2)/tm
	+0.5*masssquared_[2]/tm2);
	zLimits(limits);
	if(zLimits().second<zLimits().first) {
	t=-1.0*GeV2;
	return false;
	}
	// guess values of t and z
	t = guesst(told,2,ids_,enhance,ids_[1]==ids_[2]);
	z(guessz(2,ids_));
	// actual values for z-limits
	if(t<masssquared_[0]) {
	t=-1.0*GeV2;
	return false;
	}
	tm2 = t-masssquared_[0];
	tm = sqrt(tm2);
	limits=make_pair(sqr(minmass/masses_[0]),
	1.-sqrt(masssquared_[2]+pT2min()+
	0.25*sqr(masssquared_[2])/tm2)/tm
	+0.5*masssquared_[2]/tm2);
	zLimits(limits);
	if(t>tmax\|\|zLimits().second<zLimits().first) {
	t=-1.0*GeV2;
	return false;
	}
	else
	return true;
	}

	bool QTildeSudakov::computeTimeLikeLimits(Energy2 & t) {
	if (t < 1e-20 * GeV2) {
	t=-1.*GeV2;
	return false;
	}
	// special case for gluon radiating
	pair<double,double> limits;
	if(ids_[0]==ParticleID::g\|\|ids_[0]==ParticleID::gamma) {
	// no emission possible
	if(t<16.*masssquared_[1]) {
	t=-1.*GeV2;
	return false;
	}
	// overestimate of the limits
	limits.first = 0.5(1.-sqrt(1.-4.sqrt((masssquared_[1]+pT2min())/t)));
	limits.second = 1.-limits.first;
	}
	// special case for radiated particle is gluon
	else if(ids_[2]==ParticleID::g\|\|ids_[2]==ParticleID::gamma) {
	limits.first = sqrt((masssquared_[1]+pT2min())/t);
	limits.second = 1.-sqrt((masssquared_[2]+pT2min())/t);
	}
	else if(ids_[1]==ParticleID::g\|\|ids_[1]==ParticleID::gamma) {
	limits.second = sqrt((masssquared_[2]+pT2min())/t);
	limits.first = 1.-sqrt((masssquared_[1]+pT2min())/t);
	}
	else {
	limits.first = (masssquared_[1]+pT2min())/t;
	limits.second = 1.-(masssquared_[2]+pT2min())/t;
	}
	if(limits.first>=limits.second) {
	t=-1.*GeV2;
	return false;
	}
	zLimits(limits);
	return true;
	}

	bool QTildeSudakov::computeSpaceLikeLimits(Energy2 & t, double x) {
	if (t < 1e-20 * GeV2) {
	t=-1.*GeV2;
	return false;
	}
	pair<double,double> limits;
	// compute the limits
	limits.first = x;
	double yy = 1.+0.5*masssquared_[2]/t;
	limits.second = yy - sqrt(sqr(yy)-1.+pT2min()/t);
	// return false if lower>upper
	zLimits(limits);
	if(limits.second<limits.first) {
	t=-1.*GeV2;
	return false;
	}
	else
	return true;
	}

	+namespace {
	+
	+tShowerParticlePtr findCorrelationPartner(ShowerParticle & particle,
	+ bool forward,
	+ ShowerInteraction::Type inter) {
	+ tPPtr child = &particle;
	+ tPPtr mother;
	+ if(forward) {
	+ mother = !particle.parents().empty() ? particle.parents()[0] : tPPtr();
	+ }
	+ else {
	+ mother = particle.children().size()==2 ? &particle : tPPtr();
	+ }
	+ tShowerParticlePtr partner;
	+ while(mother) {
	+ tPPtr otherChild;
	+ if(forward) {
	+ for (unsigned int ix=0;ix<mother->children().size();++ix) {
	+ if(mother->children()[ix]!=child) {
	+ otherChild = mother->children()[ix];
	+ break;
	+ }
	+ }
	+ }
	+ else {
	+ otherChild = mother->children()[1];
	+ }
	+ tShowerParticlePtr other = dynamic_ptr_cast<tShowerParticlePtr>(otherChild);
	+ if((inter==ShowerInteraction::QCD && otherChild->dataPtr()->coloured()) \|\|
	+ (inter==ShowerInteraction::QED && otherChild->dataPtr()->charged())) {
	+ partner = other;
	+ break;
	+ }
	+ if(forward && !other->isFinalState()) {
	+ partner = dynamic_ptr_cast<tShowerParticlePtr>(mother);
	+ break;
	+ }
	+ child = mother;
	+ if(forward) {
	+ mother = ! mother->parents().empty() ? mother->parents()[0] : tPPtr();
	+ }
	+ else {
	+ if(mother->children()[0]->children().size()!=2)
	+ break;
	+ mother = mother->children()[0];
	+ }
	+ }
	+ if(!partner) {
	+ if(forward) {
	+ partner = dynamic_ptr_cast<tShowerParticlePtr>( child)->partner();
	+ }
	+ else {
	+ if(mother) {
	+ tShowerParticlePtr parent;
	+ if(!mother->children().empty()) {
	+ parent = dynamic_ptr_cast<tShowerParticlePtr>(mother->children()[0]);
	+ }
	+ if(!parent) {
	+ parent = dynamic_ptr_cast<tShowerParticlePtr>(mother);
	+ }
	+ partner = parent->partner();
	+ }
	+ else {
	+ partner = dynamic_ptr_cast<tShowerParticlePtr>(&particle)->partner();
	+ }
	+ }
	+ }
	+ return partner;
	+}
	+
	+}
	+
	double QTildeSudakov::generatePhiForward(ShowerParticle & particle,
	const IdList & ids,
	ShoKinPtr kinematics) {
	// no correlations, return flat phi
	if(! ShowerHandler::currentHandler()->evolver()->correlations())
	return Constants::twopi*UseRandom::rnd();
	- // get the spin density matrix and the mapping
	- RhoDMatrix mapping;
	- SpinPtr inspin;
	- bool needMapping = getMapping(inspin,mapping,particle,kinematics);
	- // set the decayed flag
	- inspin->decay();
	- // get the spin density matrix
	- RhoDMatrix rho=inspin->rhoMatrix();
	- // map to the shower basis if needed
	- if(needMapping) {
	- RhoDMatrix rhop(rho.iSpin(),false);
	- for(int ixa=0;ixa<rho.iSpin();++ixa) {
	- for(int ixb=0;ixb<rho.iSpin();++ixb) {
	- for(int iya=0;iya<rho.iSpin();++iya) {
	- for(int iyb=0;iyb<rho.iSpin();++iyb) {
	- rhop(ixa,ixb) += rho(iya,iyb)mapping(iya,ixa)conj(mapping(iyb,ixb));
	+ // get the kinematic variables
	+ double z = kinematics->z();
	+ Energy2 t = z(1.-z)sqr(kinematics->scale());
	+ Energy pT = kinematics->pT();
	+ // if soft correlations
	+ Energy2 pipj,pik;
	+ bool canBeSoft[2] = {ids[1]==ParticleID::g \|\| ids[1]==ParticleID::gamma,
	+ ids[2]==ParticleID::g \|\| ids[2]==ParticleID::gamma };
	+ vector<Energy2> pjk(3,ZERO);
	+ Energy2 m12(ZERO),m22(ZERO);
	+ InvEnergy2 aziMax(ZERO);
	+ bool softAllowed = ShowerHandler::currentHandler()->evolver()->softCorrelations()&&
	+ (canBeSoft[0] \|\| canBeSoft[1]);
	+ bool swapOrder;
	+ if(softAllowed) {
	+ // find the partner for the soft correlations
	+ tShowerParticlePtr partner=findCorrelationPartner(particle,true,splittingFn()->interactionType());
	+ // remember we want the softer gluon
	+ swapOrder = !canBeSoft[1] \|\| (canBeSoft[0] && canBeSoft[1] && z < 0.5);
	+ double zFact = !swapOrder ? (1.-z) : z;
	+ // compute the transforms to the shower reference frame
	+ // first the boost
	+ vector<Lorentz5Momentum> basis = kinematics->getBasis();
	+ Lorentz5Momentum pVect = basis[0];
	+ Lorentz5Momentum nVect = basis[1];
	+ Boost beta_bb;
	+ if(kinematics->frame()==ShowerKinematics::BackToBack) {
	+ beta_bb = -(pVect + nVect).boostVector();
	+ }
	+ else if(kinematics->frame()==ShowerKinematics::Rest) {
	+ beta_bb = -pVect.boostVector();
	+ }
	+ else
	+ assert(false);
	+ pVect.boost(beta_bb);
	+ nVect.boost(beta_bb);
	+ Axis axis;
	+ if(kinematics->frame()==ShowerKinematics::BackToBack) {
	+ axis = pVect.vect().unit();
	+ }
	+ else if(kinematics->frame()==ShowerKinematics::Rest) {
	+ axis = nVect.vect().unit();
	+ }
	+ else
	+ assert(false);
	+ // and then the rotation
	+ LorentzRotation rot;
	+ if(axis.perp2()>0.) {
	+ double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
	+ rot.rotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
	+ }
	+ else if(axis.z()<0.) {
	+ rot.rotate(Constants::pi,Axis(1.,0.,0.));
	+ }
	+ rot.invert();
	+ pVect *= rot;
	+ nVect *= rot;
	+ // shower parameters
	+ Energy2 pn = pVect*nVect;
	+ Energy2 m2 = pVect.m2();
	+ double alpha0 = particle.showerParameters().alpha;
	+ double beta0 = 0.5/alpha0/pn*
	+ (sqr(particle.dataPtr()->mass())-sqr(alpha0)*m2+sqr(particle.showerParameters().pt));
	+ Lorentz5Momentum qperp0(particle.showerParameters().ptx,
	+ particle.showerParameters().pty,ZERO,ZERO);
	+ assert(partner);
	+ Lorentz5Momentum pj = partner->momentum();
	+ pj.boost(beta_bb);
	+ pj *= rot;
	+ // compute the two phi independent dot products
	+ pik = 0.5zFact(sqr(alpha0)m2 - sqr(particle.showerParameters().pt) + 2.alpha0beta0pn )
	+ +0.5*sqr(pT)/zFact;
	+ Energy2 dot1 = pj*pVect;
	+ Energy2 dot2 = pj*nVect;
	+ Energy2 dot3 = pj*qperp0;
	+ pipj = alpha0dot1+beta0dot2+dot3;
	+ // compute the constants for the phi dependent dot product
	+ pjk[0] = zFact(alpha0dot1+dot3-0.5dot2/pn(alpha0*m2-sqr(particle.showerParameters().pt)/alpha0))
	+ +0.5sqr(pT)dot2/pn/zFact/alpha0;
	+ pjk[1] = (pj.x() - dot2/alpha0/pnqperp0.x())pT;
	+ pjk[2] = (pj.y() - dot2/alpha0/pnqperp0.y())pT;
	+ m12 = sqr(particle.dataPtr()->mass());
	+ m22 = sqr(partner->dataPtr()->mass());
	+ if(swapOrder) {
	+ pjk[1] *= -1.;
	+ pjk[2] *= -1.;
	+ }
	+ Energy2 mag = sqrt(sqr(pjk[1])+sqr(pjk[2]));
	+ aziMax = -m12/sqr(pik) -m22/sqr(pjk[0]+mag) +2.*pipj/pik/(pjk[0]-mag);
	+ }
	+ // if spin correlations
	+ vector<pair<int,Complex> > wgts;
	+ if(ShowerHandler::currentHandler()->evolver()->spinCorrelations()) {
	+ // get the spin density matrix and the mapping
	+ RhoDMatrix mapping;
	+ SpinPtr inspin;
	+ bool needMapping = getMapping(inspin,mapping,particle,kinematics);
	+ // set the decayed flag
	+ inspin->decay();
	+ // get the spin density matrix
	+ RhoDMatrix rho=inspin->rhoMatrix();
	+ // map to the shower basis if needed
	+ if(needMapping) {
	+ RhoDMatrix rhop(rho.iSpin(),false);
	+ for(int ixa=0;ixa<rho.iSpin();++ixa) {
	+ for(int ixb=0;ixb<rho.iSpin();++ixb) {
	+ for(int iya=0;iya<rho.iSpin();++iya) {
	+ for(int iyb=0;iyb<rho.iSpin();++iyb) {
	+ rhop(ixa,ixb) += rho(iya,iyb)mapping(iya,ixa)conj(mapping(iyb,ixb));
	+ }
	}
	}
	}
	+ rhop.normalize();
	+ rho = rhop;
	}
	- rhop.normalize();
	- rho = rhop;
	+ // calculate the weights
	+ wgts = splittingFn()->generatePhiForward(z,t,ids,rho);
	}
	- // get the kinematic variables
	- double z = kinematics->z();
	- Energy2 t = z(1.-z)sqr(kinematics->scale());
	+ else {
	+ wgts = vector<pair<int,Complex> >(1,make_pair(0,1.));
	+ }
	// generate the azimuthal angle
	- double phi;
	- Complex wgt;
	- vector<pair<int,Complex> >
	- wgts = splittingFn()->generatePhiForward(z,t,ids,rho);
	+ double phi,wgt;
	static const Complex ii(0.,1.);
	+ unsigned int ntry(0);
	+ double phiMax(0.),wgtMax(0.);
	do {
	phi = Constants::twopi*UseRandom::rnd();
	- wgt = 0.;
	+ // first the spin correlations bit (gives 1 if correlations off)
	+ Complex spinWgt = 0.;
	for(unsigned int ix=0;ix<wgts.size();++ix) {
	if(wgts[ix].first==0)
	- wgt += wgts[ix].second;
	+ spinWgt += wgts[ix].second;
	else
	- wgt += exp(double(wgts[ix].first)iiphi)*wgts[ix].second;
	+ spinWgt += exp(double(wgts[ix].first)iiphi)*wgts[ix].second;
	}
	- if(wgt.real()-1.>1e-10) {
	- cerr << "Forward weight problem " << wgt << " " << wgt.real()-1.
	- << " " << ids[0] << " " << ids[1] << " " << ids[2] << " " << " " << z << " " << phi << "\n";
	- cerr << "Weights \n";
	+ wgt = spinWgt.real();
	+ if(wgt-1.>1e-10) {
	+ generator()->log() << "Forward spin weight problem " << wgt << " " << wgt-1.
	+ << " " << ids[0] << " " << ids[1] << " " << ids[2] << " " << " " << phi << "\n";
	+ generator()->log() << "Weights \n";
	for(unsigned int ix=0;ix<wgts.size();++ix)
	- cerr << wgts[ix].first << " " << wgts[ix].second << "\n";
	+ generator()->log() << wgts[ix].first << " " << wgts[ix].second << "\n";
	}
	+ // soft correlations bit
	+ double aziWgt = 1.;
	+ if(softAllowed) {
	+ Energy2 dot = pjk[0]+pjk[1]cos(phi)+pjk[2]sin(phi);
	+ aziWgt = (-m12/sqr(pik) -m22/sqr(dot) +2.*pipj/pik/dot)/aziMax;
	+ if(aziWgt-1.>1e-10\|\|aziWgt<-1e-10) {
	+ generator()->log() << "Forward soft weight problem " << aziWgt << " " << aziWgt-1.
	+ << " " << ids[0] << " " << ids[1] << " " << ids[2] << " " << " " << phi << "\n";
	+ }
	+ }
	+ wgt *= aziWgt;
	+ if(wgt>wgtMax) {
	+ phiMax = phi;
	+ wgtMax = wgt;
	+ }
	+ ++ntry;
	}
	- while(wgt.real()<UseRandom::rnd());
	+ while(wgt<UseRandom::rnd()&&ntry<10000);
	+ if(ntry==10000) {
	+ phi = phiMax;
	+ generator()->log() << "Too many tries to generate phi in forward evolution\n";
	+ }
	// return the azimuthal angle
	return phi;
	}

	double QTildeSudakov::generatePhiBackward(ShowerParticle & particle,
	const IdList & ids,
	ShoKinPtr kinematics) {
	// no correlations, return flat phi
	if(! ShowerHandler::currentHandler()->evolver()->correlations())
	return Constants::twopi*UseRandom::rnd();
	- // get the spin density matrix and the mapping
	- RhoDMatrix mapping;
	- SpinPtr inspin;
	- bool needMapping = getMapping(inspin,mapping,particle,kinematics);
	- // set the decayed flag (counterintuitive but going backward)
	- inspin->decay();
	- // get the spin density matrix
	- RhoDMatrix rho=inspin->DMatrix();
	- // map to the shower basis if needed
	- if(needMapping) {
	- RhoDMatrix rhop(rho.iSpin(),false);
	- for(int ixa=0;ixa<rho.iSpin();++ixa) {
	- for(int ixb=0;ixb<rho.iSpin();++ixb) {
	- for(int iya=0;iya<rho.iSpin();++iya) {
	- for(int iyb=0;iyb<rho.iSpin();++iyb) {
	- rhop(ixa,ixb) += rho(iya,iyb)mapping(iya,ixa)conj(mapping(iyb,ixb));
	- }
	- }
	- }
	- }
	- rhop.normalize();
	- rho = rhop;
	- }
	// get the kinematic variables
	double z = kinematics->z();
	Energy2 t = (1.-z)*sqr(kinematics->scale())/z;
	+ Energy pT = kinematics->pT();
	+ // if soft correlations
	+ bool softAllowed = ShowerHandler::currentHandler()->evolver()->softCorrelations()&&
	+ (ids[2]==ParticleID::g \|\| ids[2]==ParticleID::gamma);
	+ Energy2 pipj,pik,m12(ZERO),m22(ZERO);
	+ vector<Energy2> pjk(3,ZERO);
	+ InvEnergy2 aziMax(ZERO);
	+ if(softAllowed) {
	+ // find the partner for the soft correlations
	+ tShowerParticlePtr partner=findCorrelationPartner(particle,false,splittingFn()->interactionType());
	+ double zFact = (1.-z);
	+ // compute the transforms to the shower reference frame
	+ // first the boost
	+ vector<Lorentz5Momentum> basis = kinematics->getBasis();
	+ Lorentz5Momentum pVect = basis[0];
	+ Lorentz5Momentum nVect = basis[1];
	+ assert(kinematics->frame()==ShowerKinematics::BackToBack);
	+ Boost beta_bb = -(pVect + nVect).boostVector();
	+ pVect.boost(beta_bb);
	+ nVect.boost(beta_bb);
	+ Axis axis = pVect.vect().unit();
	+ // and then the rotation
	+ LorentzRotation rot;
	+ if(axis.perp2()>0.) {
	+ double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
	+ rot.rotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
	+ }
	+ else if(axis.z()<0.) {
	+ rot.rotate(Constants::pi,Axis(1.,0.,0.));
	+ }
	+ rot.invert();
	+ pVect *= rot;
	+ nVect *= rot;
	+ // shower parameters
	+ Energy2 pn = pVect*nVect;
	+ Energy2 m2 = pVect.m2();
	+ double alpha0 = particle.x();
	+ double beta0 = -0.5/alpha0/pnsqr(alpha0)m2;
	+ Lorentz5Momentum pj = partner->momentum();
	+ pj.boost(beta_bb);
	+ pj *= rot;
	+ double beta2 = 0.5(1.-zFact)(sqr(alpha0zFact/(1.-zFact))m2+sqr(pT))/alpha0/zFact/pn;
	+ // compute the two phi independent dot products
	+ Energy2 dot1 = pj*pVect;
	+ Energy2 dot2 = pj*nVect;
	+ pipj = alpha0dot1+beta0dot2;
	+ pik = alpha0(alpha0zFact/(1.-zFact)m2+pn(beta2+zFact/(1.-zFact)*beta0));
	+ // compute the constants for the phi dependent dot product
	+ pjk[0] = alpha0zFact/(1.-zFact)dot1+beta2*dot2;
	+ pjk[1] = pj.x()*pT;
	+ pjk[2] = pj.y()*pT;
	+ m12 = ZERO;
	+ m22 = sqr(partner->dataPtr()->mass());
	+ Energy2 mag = sqrt(sqr(pjk[1])+sqr(pjk[2]));
	+ aziMax = -m12/sqr(pik) -m22/sqr(pjk[0]+mag) +2.*pipj/pik/(pjk[0]-mag);
	+ }
	+ // if spin correlations
	+ vector<pair<int,Complex> > wgts;
	+ if(ShowerHandler::currentHandler()->evolver()->spinCorrelations()) {
	+ // get the spin density matrix and the mapping
	+ RhoDMatrix mapping;
	+ SpinPtr inspin;
	+ bool needMapping = getMapping(inspin,mapping,particle,kinematics);
	+ // set the decayed flag (counterintuitive but going backward)
	+ inspin->decay();
	+ // get the spin density matrix
	+ RhoDMatrix rho=inspin->DMatrix();
	+ // map to the shower basis if needed
	+ if(needMapping) {
	+ RhoDMatrix rhop(rho.iSpin(),false);
	+ for(int ixa=0;ixa<rho.iSpin();++ixa) {
	+ for(int ixb=0;ixb<rho.iSpin();++ixb) {
	+ for(int iya=0;iya<rho.iSpin();++iya) {
	+ for(int iyb=0;iyb<rho.iSpin();++iyb) {
	+ rhop(ixa,ixb) += rho(iya,iyb)mapping(iya,ixa)conj(mapping(iyb,ixb));
	+ }
	+ }
	+ }
	+ }
	+ rhop.normalize();
	+ rho = rhop;
	+ }
	+ wgts = splittingFn()->generatePhiBackward(z,t,ids,rho);
	+ }
	+ else {
	+ wgts = vector<pair<int,Complex> >(1,make_pair(0,1.));
	+ }
	// generate the azimuthal angle
	- double phi;
	- Complex wgt;
	- vector<pair<int,Complex> >
	- wgts = splittingFn()->generatePhiBackward(z,t,ids,rho);
	+ double phi,wgt;
	static const Complex ii(0.,1.);
	+ unsigned int ntry(0);
	+ double phiMax(0.),wgtMax(0.);
	+
	+
	do {
	phi = Constants::twopi*UseRandom::rnd();
	- wgt = 0.;
	+ Complex spinWgt = 0.;
	for(unsigned int ix=0;ix<wgts.size();++ix) {
	if(wgts[ix].first==0)
	- wgt += wgts[ix].second;
	+ spinWgt += wgts[ix].second;
	else
	- wgt += exp(double(wgts[ix].first)iiphi)*wgts[ix].second;
	+ spinWgt += exp(double(wgts[ix].first)iiphi)*wgts[ix].second;
	}
	- if(wgt.real()-1.>1e-10) {
	- cerr << "Backward weight problem " << wgt << " " << wgt.real()-1.
	+ wgt = spinWgt.real();
	+ if(wgt-1.>1e-10) {
	+ generator()->log() << "Backward weight problem " << wgt << " " << wgt-1.
	<< " " << ids[0] << " " << ids[1] << " " << ids[2] << " " << " " << z << " " << phi << "\n";
	- cerr << "Weights \n";
	+ generator()->log() << "Weights \n";
	for(unsigned int ix=0;ix<wgts.size();++ix)
	- cerr << wgts[ix].first << " " << wgts[ix].second << "\n";
	+ generator()->log() << wgts[ix].first << " " << wgts[ix].second << "\n";
	}
	+ // soft correlations bit
	+ double aziWgt = 1.;
	+ if(softAllowed) {
	+ Energy2 dot = pjk[0]+pjk[1]cos(phi)+pjk[2]sin(phi);
	+ aziWgt = (-m12/sqr(pik) -m22/sqr(dot) +2.*pipj/pik/dot)/aziMax;
	+ if(aziWgt-1.>1e-10\|\|aziWgt<-1e-10) {
	+ generator()->log() << "Forward soft weight problem " << aziWgt << " " << aziWgt-1.
	+ << " " << ids[0] << " " << ids[1] << " " << ids[2] << " " << " " << phi << "\n";
	+ }
	+ }
	+ wgt *= aziWgt;
	+ if(wgt>wgtMax) {
	+ phiMax = phi;
	+ wgtMax = wgt;
	+ }
	+ ++ntry;
	}
	- while(wgt.real()<UseRandom::rnd());
	+ while(wgt<UseRandom::rnd()&&ntry<10000);
	+ if(ntry==10000) {
	+ phi = phiMax;
	+ generator()->log() << "Too many tries to generate phi in backward evolution\n";
	+ }
	// return the azimuthal angle
	return phi;
	}

	double QTildeSudakov::generatePhiDecay(ShowerParticle & particle,
	const IdList & ids,
	ShoKinPtr kinematics) {
	- return Constants::twopi*UseRandom::rnd();
	- // cerr << particle.isFinalState() << " " << particle << "\n";
	- // cerr << particle.spinInfo() << "\n";
	- // if(particle.spinInfo()) {
	- // cerr << "testing spin info " << particle.spinInfo()->productionVertex() << " "
	- // << particle.spinInfo()->decayVertex() << " "
	- // << particle.spinInfo()->developed() << " " << particle.spinInfo()->decayed() << "\n";
	- // if(particle.spinInfo()->productionVertex()) {
	- // cerr << "production "
	- // << particle.spinInfo()->productionVertex()->incoming().size() << " "
	- // << particle.spinInfo()->productionVertex()->outgoing().size() << "\n";
	- // }
	- // if(particle.spinInfo()->decayVertex()) {
	- // cerr << "decay "
	- // << particle.spinInfo()->decayVertex()->incoming().size() << " "
	- // << particle.spinInfo()->decayVertex()->outgoing().size() << "\n";
	- // }
	- // }
	-
	- // cerr << "testing variables " << kinematics->z() << " " << kinematics->scale()/GeV << " "
	- // << particle.mass()/GeV << "\n";;
	-
	- // // no correlations, return flat phi
	- // if(! ShowerHandler::currentHandler()->evolver()->correlations())
	- // return Constants::twopi*UseRandom::rnd();
	- // // get the spin density matrix and the mapping
	- // RhoDMatrix mapping;
	- // SpinPtr inspin;
	- // bool needMapping = getMapping(inspin,mapping,particle,kinematics);
	- // assert(false);
	- //
	- // // set the decayed flag
	- // inspin->decay();
	- // // get the spin density matrix
	- // RhoDMatrix rho=inspin->rhoMatrix();
	- // // map to the shower basis if needed
	- // if(needMapping) {
	- // RhoDMatrix rhop(rho.iSpin(),false);
	- // for(int ixa=0;ixa<rho.iSpin();++ixa) {
	- // for(int ixb=0;ixb<rho.iSpin();++ixb) {
	- // for(int iya=0;iya<rho.iSpin();++iya) {
	- // for(int iyb=0;iyb<rho.iSpin();++iyb) {
	- // rhop(ixa,ixb) += rho(iya,iyb)mapping(iya,ixa)conj(mapping(iyb,ixb));
	- // }
	- // }
	- // }
	- // }
	- // rhop.normalize();
	- // rho = rhop;
	- // }
	- // // get the kinematic variables
	- // double z = kinematics->z();
	- // Energy2 t = z(1.-z)sqr(kinematics->scale());
	- // // generate the azimuthal angle
	- // double phi;
	- // Complex wgt;
	- // vector<pair<int,Complex> >
	- // wgts = splittingFn()->generatePhiForward(z,t,ids,rho);
	- // static const Complex ii(0.,1.);
	- // do {
	- // phi = Constants::twopi*UseRandom::rnd();
	- // wgt = 0.;
	- // for(unsigned int ix=0;ix<wgts.size();++ix) {
	- // if(wgts[ix].first==0)
	- // wgt += wgts[ix].second;
	- // else
	- // wgt += exp(double(wgts[ix].first)iiphi)*wgts[ix].second;
	- // }
	- // if(wgt.real()-1.>1e-10) {
	- // cerr << "Forward weight problem " << wgt << " " << wgt.real()-1.
	- // << " " << ids[0] << " " << ids[1] << " " << ids[2] << " " << " " << z << " " << phi << "\n";
	- // cerr << "Weights \n";
	- // for(unsigned int ix=0;ix<wgts.size();++ix)
	- // cerr << wgts[ix].first << " " << wgts[ix].second << "\n";
	- // }
	- // }
	- // while(wgt.real()<UseRandom::rnd());
	- // // compute the matrix element for spin correlations
	- // DecayMatrixElement me(splittingFn()->matrixElement(z,t,ids,phi));
	- // // create the vertex
	- // SVertexPtr Svertex(new_ptr(ShowerVertex()));
	- // // set the matrix element
	- // Svertex->ME().reset(me);
	- // // set the incoming particle for the vertex
	- // inspin->decayVertex(Svertex);
	- // // return the azimuthal angle
	- // return phi;
	+ // only soft correlations in this case
	+ bool canBeSoft[2] = {ids[1]==ParticleID::g \|\| ids[1]==ParticleID::gamma,
	+ ids[2]==ParticleID::g \|\| ids[2]==ParticleID::gamma };
	+ // no correlations, return flat phi
	+ if( ! (ShowerHandler::currentHandler()->evolver()->softCorrelations()&&
	+ (canBeSoft[0] \|\| canBeSoft[1]) ) )
	+ return Constants::twopi*UseRandom::rnd();
	+ // get the kinematic variables
	+ double z = kinematics->z();
	+ Energy pT = kinematics->pT();
	+ // if soft correlations
	+ // find the partner for the soft correlations
	+ tShowerParticlePtr partner = findCorrelationPartner(particle,true,splittingFn()->interactionType());
	+ // remember we want the softer gluon
	+ bool swapOrder = !canBeSoft[1] \|\| (canBeSoft[0] && canBeSoft[1] && z < 0.5);
	+ double zFact = !swapOrder ? (1.-z) : z;
	+ vector<Lorentz5Momentum> basis = kinematics->getBasis();
	+ // compute the transforms to the shower reference frame
	+ // first the boost
	+ Lorentz5Momentum pVect = basis[0];
	+ Lorentz5Momentum nVect = basis[1];
	+ assert(kinematics->frame()==ShowerKinematics::Rest);
	+ Boost beta_bb = -pVect.boostVector();
	+ pVect.boost(beta_bb);
	+ nVect.boost(beta_bb);
	+ Axis axis = nVect.vect().unit();
	+ // and then the rotation
	+ LorentzRotation rot;
	+ if(axis.perp2()>0.) {
	+ double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
	+ rot.rotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
	+ }
	+ else if(axis.z()<0.) {
	+ rot.rotate(Constants::pi,Axis(1.,0.,0.));
	+ }
	+ rot.invert();
	+ pVect *= rot;
	+ nVect *= rot;
	+ // shower parameters
	+ Energy2 pn = pVect*nVect;
	+ Energy2 m2 = pVect.m2();
	+ double alpha0 = particle.showerParameters().alpha;
	+ double beta0 = 0.5/alpha0/pn*
	+ (sqr(particle.dataPtr()->mass())-sqr(alpha0)*m2+sqr(particle.showerParameters().pt));
	+ Lorentz5Momentum qperp0(particle.showerParameters().ptx,
	+ particle.showerParameters().pty,ZERO,ZERO);
	+ Lorentz5Momentum pj = partner->momentum();
	+ pj.boost(beta_bb);
	+ pj *= rot;
	+ // compute the two phi independent dot products
	+ Energy2 pik = 0.5zFact(sqr(alpha0)m2 - sqr(particle.showerParameters().pt) + 2.alpha0beta0pn )
	+ +0.5*sqr(pT)/zFact;
	+ Energy2 dot1 = pj*pVect;
	+ Energy2 dot2 = pj*nVect;
	+ Energy2 dot3 = pj*qperp0;
	+ Energy2 pipj = alpha0dot1+beta0dot2+dot3;
	+ // compute the constants for the phi dependent dot product
	+ vector<Energy2> pjk(3,ZERO);
	+ pjk[0] = zFact(alpha0dot1+dot3-0.5dot2/pn(alpha0*m2-sqr(particle.showerParameters().pt)/alpha0))
	+ +0.5sqr(pT)dot2/pn/zFact/alpha0;
	+ pjk[1] = (pj.x() - dot2/alpha0/pnqperp0.x())pT;
	+ pjk[2] = (pj.y() - dot2/alpha0/pnqperp0.y())pT;
	+ Energy2 m12 = sqr(particle.dataPtr()->mass());
	+ Energy2 m22 = sqr(partner->dataPtr()->mass());
	+ if(swapOrder) {
	+ pjk[1] *= -1.;
	+ pjk[2] *= -1.;
	+ }
	+ Energy2 mag = sqrt(sqr(pjk[1])+sqr(pjk[2]));
	+ InvEnergy2 aziMax = -m12/sqr(pik) -m22/sqr(pjk[0]+mag) +2.*pipj/pik/(pjk[0]-mag);
	+ // generate the azimuthal angle
	+ double phi,wgt;
	+ unsigned int ntry(0);
	+ double phiMax(0.),wgtMax(0.);
	+ do {
	+ phi = Constants::twopi*UseRandom::rnd();
	+ Energy2 dot = pjk[0]+pjk[1]cos(phi)+pjk[2]sin(phi);
	+ wgt = (-m12/sqr(pik) -m22/sqr(dot) +2.*pipj/pik/dot)/aziMax;
	+ if(wgt-1.>1e-10\|\|wgt<-1e-10) {
	+ generator()->log() << "Decay soft weight problem " << wgt << " " << wgt-1.
	+ << " " << ids[0] << " " << ids[1] << " " << ids[2] << " " << " " << phi << "\n";
	+ }
	+ if(wgt>wgtMax) {
	+ phiMax = phi;
	+ wgtMax = wgt;
	+ }
	+ ++ntry;
	+ }
	+ while(wgt<UseRandom::rnd()&&ntry<10000);
	+ if(ntry==10000) {
	+ phi = phiMax;
	+ generator()->log() << "Too many tries to generate phi\n";
	+ }
	+ // return the azimuthal angle
	+ return phi;
	}


	Energy QTildeSudakov::calculateScale(double zin, Energy pt, IdList ids,
	unsigned int iopt) {
	Energy2 tmin;
	initialize(ids,tmin,false);
	// final-state branching
	if(iopt==0) {
	Energy2 scale=(sqr(pt)+masssquared_[1](1.-zin)+masssquared_[2]zin);
	if(ids[0]!=ParticleID::g) scale -= zin(1.-zin)masssquared_[0];
	scale /= sqr(zin*(1-zin));
	return scale<=ZERO ? sqrt(tmin) : sqrt(scale);
	}
	else if(iopt==1) {
	Energy2 scale=(sqr(pt)+zin*masssquared_[2])/sqr(1.-zin);
	return scale<=ZERO ? sqrt(tmin) : sqrt(scale);
	}
	else if(iopt==2) {
	Energy2 scale = (sqr(pt)+zin*masssquared_[2])/sqr(1.-zin)+masssquared_[0];
	return scale<=ZERO ? sqrt(tmin) : sqrt(scale);
	}
	else {
	throw Exception() << "Unknown option in QTildeSudakov::calculateScale() "
	<< "iopt = " << iopt << Exception::runerror;
	}
	}

	ShoKinPtr QTildeSudakov::createFinalStateBranching(Energy scale,double z,
	double phi, Energy pt) {
	ShoKinPtr showerKin = new_ptr(FS_QTildeShowerKinematics1to2());
	showerKin->scale(scale);
	showerKin->z(z);
	showerKin->phi(phi);
	showerKin->pT(pt);
	showerKin->SudakovFormFactor(this);
	return showerKin;
	}

	ShoKinPtr QTildeSudakov::createInitialStateBranching(Energy scale,double z,
	double phi, Energy pt) {
	ShoKinPtr showerKin = new_ptr(IS_QTildeShowerKinematics1to2());
	showerKin->scale(scale);
	showerKin->z(z);
	showerKin->phi(phi);
	showerKin->pT(pt);
	showerKin->SudakovFormFactor(this);
	return showerKin;
	}

	ShoKinPtr QTildeSudakov::createDecayBranching(Energy scale,double z,
	double phi, Energy pt) {
	ShoKinPtr showerKin = new_ptr(Decay_QTildeShowerKinematics1to2());
	showerKin->scale(scale);
	showerKin->z(z);
	showerKin->phi(phi);
	showerKin->pT(pt);
	showerKin->SudakovFormFactor(this);
	return showerKin;
	}
	diff --git a/src/defaults/Shower.in b/src/defaults/Shower.in
	--- a/src/defaults/Shower.in
	+++ b/src/defaults/Shower.in
	@@ -1,268 +1,268 @@
	############################################################
	# Setup of default parton shower
	#
	# Useful switches for users are marked near the top of
	# this file.
	#
	# Don't edit this file directly, but reset the switches
	# in your own input files!
	############################################################
	library HwMPI.so
	library HwShower.so

	mkdir /Herwig/Shower
	cd /Herwig/Shower

	create Herwig::ShowerHandler ShowerHandler
	newdef ShowerHandler:MPIHandler /Herwig/UnderlyingEvent/MPIHandler
	newdef ShowerHandler:RemDecayer /Herwig/Partons/RemnantDecayer

	#####################################
	# initial setup, don't change these!
	#####################################
	create Herwig::SplittingGenerator SplittingGenerator
	create Herwig::ShowerAlphaQCD AlphaQCD
	create Herwig::ShowerAlphaQED AlphaQED
	create Herwig::Evolver Evolver
	create Herwig::QTildeModel ShowerModel
	create Herwig::QTildeFinder PartnerFinder
	create Herwig::QTildeReconstructor KinematicsReconstructor
	newdef KinematicsReconstructor:ReconstructionOption Colour

	newdef /Herwig/Partons/RemnantDecayer:AlphaS AlphaQCD
	newdef /Herwig/Partons/RemnantDecayer:AlphaEM AlphaQED

	newdef ShowerHandler:Evolver Evolver
	newdef ShowerModel:PartnerFinder PartnerFinder
	newdef ShowerModel:KinematicsReconstructor KinematicsReconstructor
	newdef Evolver:ShowerModel ShowerModel
	newdef Evolver:SplittingGenerator SplittingGenerator
	newdef Evolver:HardScaleFactor 1.0
	newdef Evolver:Interactions BothAtOnce
	-newdef Evolver:SpinCorrelations Spin
	+newdef Evolver:SpinCorrelations Both

	##################################################################
	# Intrinsic pT
	#
	# Recommended:
	# 1.9 GeV for Tevatron W/Z production.
	# 2.1 GeV for LHC W/Z production at 10 TeV
	# 2.2 GeV for LHC W/Z production at 14 TeV
	#
	# Set all parameters to 0 to disable
	##################################################################
	newdef Evolver:IntrinsicPtGaussian 1.9*GeV
	newdef Evolver:IntrinsicPtBeta 0
	newdef Evolver:IntrinsicPtGamma 0*GeV
	newdef Evolver:IntrinsicPtIptmax 0*GeV
	#############################################################
	# Main control switches for the parton shower.
	#############################################################
	newdef SplittingGenerator:ISR Yes
	newdef SplittingGenerator:FSR Yes
	#############################################################


	#############################################################
	# End of interesting user servicable section.
	#
	# Anything that follows below should only be touched if you
	# know what you're doing.
	#
	# Really.
	#############################################################
	#
	# a few default values
	newdef Evolver:MECorrMode 1
	newdef AlphaQCD:ScaleFactor 1.0
	newdef AlphaQCD:RenormalizationScaleFactor 1.0
	newdef AlphaQCD:NPAlphaS 2
	newdef AlphaQCD:Qmin 0.935
	newdef AlphaQCD:NumberOfLoops 3
	newdef AlphaQCD:InputOption 1
	newdef AlphaQCD:AlphaMZ 0.120
	#
	#
	# Lets set up all the splittings
	create Herwig::HalfHalfOneSplitFn QtoQGammaSplitFn
	set QtoQGammaSplitFn:InteractionType QED
	set QtoQGammaSplitFn:ColourStructure ChargedChargedNeutral
	set QtoQGammaSplitFn:AngularOrdered Yes

	create Herwig::HalfHalfOneSplitFn QtoQGSplitFn
	newdef QtoQGSplitFn:InteractionType QCD
	newdef QtoQGSplitFn:ColourStructure TripletTripletOctet
	set QtoQGSplitFn:AngularOrdered Yes

	create Herwig::OneOneOneSplitFn GtoGGSplitFn
	newdef GtoGGSplitFn:InteractionType QCD
	newdef GtoGGSplitFn:ColourStructure OctetOctetOctet
	set GtoGGSplitFn:AngularOrdered Yes

	create Herwig::OneHalfHalfSplitFn GtoQQbarSplitFn
	newdef GtoQQbarSplitFn:InteractionType QCD
	newdef GtoQQbarSplitFn:ColourStructure OctetTripletTriplet
	#set GtoQQbarSplitFn:AngularOrdered No
	set GtoQQbarSplitFn:AngularOrdered Yes

	create Herwig::OneHalfHalfSplitFn GammatoQQbarSplitFn
	newdef GammatoQQbarSplitFn:InteractionType QED
	newdef GammatoQQbarSplitFn:ColourStructure NeutralChargedCharged
	#set GammatoQQbarSplitFn:AngularOrdered No
	set GammatoQQbarSplitFn:AngularOrdered Yes

	create Herwig::HalfOneHalfSplitFn QtoGQSplitFn
	newdef QtoGQSplitFn:InteractionType QCD
	newdef QtoGQSplitFn:ColourStructure TripletOctetTriplet
	set QtoGQSplitFn:AngularOrdered Yes

	create Herwig::HalfOneHalfSplitFn QtoGammaQSplitFn
	newdef QtoGammaQSplitFn:InteractionType QED
	newdef QtoGammaQSplitFn:ColourStructure ChargedNeutralCharged
	set QtoGammaQSplitFn:AngularOrdered Yes

	#
	# Now the Sudakovs
	create Herwig::QTildeSudakov SudakovCommon
	newdef SudakovCommon:Alpha AlphaQCD
	newdef SudakovCommon:cutoffKinScale 2.650*GeV
	newdef SudakovCommon:PDFmax 1.0

	cp SudakovCommon QtoQGSudakov
	newdef QtoQGSudakov:SplittingFunction QtoQGSplitFn
	newdef QtoQGSudakov:PDFmax 1.9

	cp SudakovCommon QtoQGammaSudakov
	set QtoQGammaSudakov:SplittingFunction QtoQGammaSplitFn
	set QtoQGammaSudakov:Alpha AlphaQED
	set QtoQGammaSudakov:PDFmax 1.9

	cp QtoQGammaSudakov LtoLGammaSudakov

	cp SudakovCommon GtoGGSudakov
	newdef GtoGGSudakov:SplittingFunction GtoGGSplitFn
	newdef GtoGGSudakov:PDFmax 2.0

	cp SudakovCommon GtoQQbarSudakov
	newdef GtoQQbarSudakov:SplittingFunction GtoQQbarSplitFn
	newdef GtoQQbarSudakov:PDFmax 120.0

	cp SudakovCommon GammatoQQbarSudakov
	newdef GammatoQQbarSudakov:SplittingFunction GammatoQQbarSplitFn
	newdef GammatoQQbarSudakov:PDFmax 120.0

	cp SudakovCommon GtobbbarSudakov
	newdef GtobbbarSudakov:SplittingFunction GtoQQbarSplitFn
	newdef GtobbbarSudakov:PDFmax 40000.0

	cp SudakovCommon GtoccbarSudakov
	newdef GtoccbarSudakov:SplittingFunction GtoQQbarSplitFn
	newdef GtoccbarSudakov:PDFmax 2000.0

	cp SudakovCommon QtoGQSudakov
	newdef QtoGQSudakov:SplittingFunction QtoGQSplitFn

	cp SudakovCommon QtoGammaQSudakov
	newdef QtoGammaQSudakov:SplittingFunction QtoGammaQSplitFn

	cp SudakovCommon utoGuSudakov
	newdef utoGuSudakov:SplittingFunction QtoGQSplitFn
	newdef utoGuSudakov:PDFFactor OverOneMinusZ
	newdef utoGuSudakov:PDFmax 5.0

	cp SudakovCommon dtoGdSudakov
	newdef dtoGdSudakov:SplittingFunction QtoGQSplitFn
	newdef dtoGdSudakov:PDFFactor OverOneMinusZ


	#
	# Now add the final splittings
	#
	do SplittingGenerator:AddFinalSplitting u->u,g; QtoQGSudakov
	do SplittingGenerator:AddFinalSplitting d->d,g; QtoQGSudakov
	do SplittingGenerator:AddFinalSplitting s->s,g; QtoQGSudakov
	do SplittingGenerator:AddFinalSplitting c->c,g; QtoQGSudakov
	do SplittingGenerator:AddFinalSplitting b->b,g; QtoQGSudakov
	do SplittingGenerator:AddFinalSplitting t->t,g; QtoQGSudakov
	#
	do SplittingGenerator:AddFinalSplitting g->g,g; GtoGGSudakov
	#
	do SplittingGenerator:AddFinalSplitting g->u,ubar; GtoQQbarSudakov
	do SplittingGenerator:AddFinalSplitting g->d,dbar; GtoQQbarSudakov
	do SplittingGenerator:AddFinalSplitting g->s,sbar; GtoQQbarSudakov
	do SplittingGenerator:AddFinalSplitting g->c,cbar; GtoccbarSudakov
	do SplittingGenerator:AddFinalSplitting g->b,bbar; GtobbbarSudakov
	do SplittingGenerator:AddFinalSplitting g->t,tbar; GtoQQbarSudakov
	#
	do SplittingGenerator:AddFinalSplitting gamma->u,ubar; GammatoQQbarSudakov
	do SplittingGenerator:AddFinalSplitting gamma->d,dbar; GammatoQQbarSudakov
	do SplittingGenerator:AddFinalSplitting gamma->s,sbar; GammatoQQbarSudakov
	do SplittingGenerator:AddFinalSplitting gamma->c,cbar; GammatoQQbarSudakov
	do SplittingGenerator:AddFinalSplitting gamma->b,bbar; GammatoQQbarSudakov
	do SplittingGenerator:AddFinalSplitting gamma->t,tbar; GammatoQQbarSudakov
	do SplittingGenerator:AddFinalSplitting gamma->e-,e+; GammatoQQbarSudakov
	do SplittingGenerator:AddFinalSplitting gamma->mu-,mu+; GammatoQQbarSudakov
	do SplittingGenerator:AddFinalSplitting gamma->tau-,tau+; GammatoQQbarSudakov
	#
	do SplittingGenerator:AddFinalSplitting u->u,gamma; QtoQGammaSudakov
	do SplittingGenerator:AddFinalSplitting d->d,gamma; QtoQGammaSudakov
	do SplittingGenerator:AddFinalSplitting s->s,gamma; QtoQGammaSudakov
	do SplittingGenerator:AddFinalSplitting c->c,gamma; QtoQGammaSudakov
	do SplittingGenerator:AddFinalSplitting b->b,gamma; QtoQGammaSudakov
	do SplittingGenerator:AddFinalSplitting t->t,gamma; QtoQGammaSudakov

	do SplittingGenerator:AddFinalSplitting e-->e-,gamma; LtoLGammaSudakov
	do SplittingGenerator:AddFinalSplitting mu-->mu-,gamma; LtoLGammaSudakov
	do SplittingGenerator:AddFinalSplitting tau-->tau-,gamma; LtoLGammaSudakov
	#
	# Now lets add the initial splittings. Remember the form a->b,c; means
	# that particle a is the particle given and we backward branch to
	# particle b which is initial state and particle c which is final state
	#
	do SplittingGenerator:AddInitialSplitting u->u,g; QtoQGSudakov
	do SplittingGenerator:AddInitialSplitting d->d,g; QtoQGSudakov
	do SplittingGenerator:AddInitialSplitting s->s,g; QtoQGSudakov
	do SplittingGenerator:AddInitialSplitting c->c,g; QtoQGSudakov
	do SplittingGenerator:AddInitialSplitting b->b,g; QtoQGSudakov
	do SplittingGenerator:AddInitialSplitting u->u,gamma; QtoQGammaSudakov
	do SplittingGenerator:AddInitialSplitting d->d,gamma; QtoQGammaSudakov
	do SplittingGenerator:AddInitialSplitting s->s,gamma; QtoQGammaSudakov
	do SplittingGenerator:AddInitialSplitting c->c,gamma; QtoQGammaSudakov
	do SplittingGenerator:AddInitialSplitting b->b,gamma; QtoQGammaSudakov
	do SplittingGenerator:AddInitialSplitting t->t,gamma; QtoQGammaSudakov

	do SplittingGenerator:AddInitialSplitting g->g,g; GtoGGSudakov
	#
	do SplittingGenerator:AddInitialSplitting g->d,dbar; GtoQQbarSudakov
	do SplittingGenerator:AddInitialSplitting g->u,ubar; GtoQQbarSudakov
	do SplittingGenerator:AddInitialSplitting g->s,sbar; GtoQQbarSudakov
	do SplittingGenerator:AddInitialSplitting g->c,cbar; GtoccbarSudakov
	do SplittingGenerator:AddInitialSplitting g->b,bbar; GtobbbarSudakov
	#
	do SplittingGenerator:AddInitialSplitting gamma->d,dbar; GammatoQQbarSudakov
	do SplittingGenerator:AddInitialSplitting gamma->u,ubar; GammatoQQbarSudakov
	do SplittingGenerator:AddInitialSplitting gamma->s,sbar; GammatoQQbarSudakov
	do SplittingGenerator:AddInitialSplitting gamma->c,cbar; GammatoQQbarSudakov
	do SplittingGenerator:AddInitialSplitting gamma->b,bbar; GammatoQQbarSudakov
	#
	do SplittingGenerator:AddInitialSplitting d->g,d; dtoGdSudakov
	do SplittingGenerator:AddInitialSplitting u->g,u; utoGuSudakov
	do SplittingGenerator:AddInitialSplitting s->g,s; QtoGQSudakov
	do SplittingGenerator:AddInitialSplitting c->g,c; QtoGQSudakov
	do SplittingGenerator:AddInitialSplitting b->g,b; QtoGQSudakov
	do SplittingGenerator:AddInitialSplitting dbar->g,dbar; dtoGdSudakov
	do SplittingGenerator:AddInitialSplitting ubar->g,ubar; utoGuSudakov
	do SplittingGenerator:AddInitialSplitting sbar->g,sbar; QtoGQSudakov
	do SplittingGenerator:AddInitialSplitting cbar->g,cbar; QtoGQSudakov
	do SplittingGenerator:AddInitialSplitting bbar->g,bbar; QtoGQSudakov
	#
	do SplittingGenerator:AddInitialSplitting d->gamma,d; QtoGammaQSudakov
	do SplittingGenerator:AddInitialSplitting u->gamma,u; QtoGammaQSudakov
	do SplittingGenerator:AddInitialSplitting s->gamma,s; QtoGammaQSudakov
	do SplittingGenerator:AddInitialSplitting c->gamma,c; QtoGammaQSudakov
	do SplittingGenerator:AddInitialSplitting b->gamma,b; QtoGammaQSudakov
	do SplittingGenerator:AddInitialSplitting dbar->gamma,dbar; QtoGammaQSudakov
	do SplittingGenerator:AddInitialSplitting ubar->gamma,ubar; QtoGammaQSudakov
	do SplittingGenerator:AddInitialSplitting sbar->gamma,sbar; QtoGammaQSudakov
	do SplittingGenerator:AddInitialSplitting cbar->gamma,cbar; QtoGammaQSudakov
	do SplittingGenerator:AddInitialSplitting bbar->gamma,bbar; QtoGammaQSudakov

File Metadata

Mime Type: text/x-diff
Expires: Mon, Jan 20, 8:56 PM (22 h, 58 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4242353
Default Alt Text: (101 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions