No OneTemporary
Actions

Size

77 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,685 +1,714 @@
	// -- 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"
	#include "Herwig/Shower/ShowerHandler.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)
	<< theFactorizationScaleFactor << theRenormalizationScaleFactor;
	}

	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)
	>> theFactorizationScaleFactor >> theRenormalizationScaleFactor;
	}

	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::alphaSVeto(Energy2 pt2) const {
	pt2 *= sqr(renormalizationScaleFactor());
	return UseRandom::rnd() > ThePEG::Math::powi(alpha_->ratio(pt2),
	splittingFn_->interactionOrder());
	}

	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 * sqr(factorizationScaleFactor());
	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()>1e-10) {
	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()==complex<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 {
	Lorentz5Momentum porig;
	LorentzRotation rot = boostToShower(showerkin->getBasis(),showerkin->frame(),porig);
	Helicity::Direction dir = particle.isFinalState() ? outgoing : incoming;
	if(particle.dataPtr()->iSpin()==PDT::Spin0) {
	assert(false);
	}
	else if(particle.dataPtr()->iSpin()==PDT::Spin1Half) {
	output = createFermionSpinInfo(particle,porig,rot,dir);
	}
	else if(particle.dataPtr()->iSpin()==PDT::Spin1) {
	output = createVectorSpinInfo(particle,porig,rot,dir);
	}
	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) {
	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_;
	//symmetry of FS gluon splitting
	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;
	}
	+
	+RhoDMatrix SudakovFormFactor::extractRhoMatrix(ShowerParticle & particle,
	+ ShoKinPtr kinematics,bool forward) {
	+ // 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 = forward ? inspin->rhoMatrix() : 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;
	+ }
	+ return rho;
	+}
	+
	diff --git a/Shower/Base/SudakovFormFactor.h b/Shower/Base/SudakovFormFactor.h
	--- a/Shower/Base/SudakovFormFactor.h
	+++ b/Shower/Base/SudakovFormFactor.h
	@@ -1,701 +1,704 @@
	// -- C++ --
	//
	// SudakovFormFactor.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_SudakovFormFactor_H
	#define HERWIG_SudakovFormFactor_H
	//
	// This is the declaration of the SudakovFormFactor class.
	//

	#include "ThePEG/Interface/Interfaced.h"
	#include "Herwig/Shower/SplittingFunctions/SplittingFunction.h"
	#include "Herwig/Shower/Couplings/ShowerAlpha.h"
	#include "Herwig/Shower/SplittingFunctions/SplittingGenerator.fh"
	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/PDF/BeamParticleData.h"
	#include "ThePEG/EventRecord/RhoDMatrix.h"
	#include "ThePEG/EventRecord/SpinInfo.h"
	#include "ShowerKinematics.fh"
	#include "SudakovFormFactor.fh"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* A typedef for the BeamParticleData
	*/
	typedef Ptr<BeamParticleData>::transient_const_pointer tcBeamPtr;

	/** \ingroup Shower
	*
	* This is the definition of the Sudakov form factor class. In general this
	* is the base class for the implementation of Sudakov form factors in Herwig.
	* The methods generateNextTimeBranching(), generateNextDecayBranching() and
	* generateNextSpaceBranching need to be implemented in classes inheriting from this
	* one.
	*
	* In addition a number of methods are implemented to assist with the calculation
	* of the form factor using the veto algorithm in classes inheriting from this one.
	*
	* In general the Sudakov form-factor, for final-state radiation, is given
	* by
	* \f[\Delta_{ba}(\tilde{q}_{i+1},\tilde{q}_i)=
	* \exp\left\{
	* -\int^{\tilde{q}^2_i}_{\tilde{q}^2_{i+1}}
	* \frac{{\rm d}\tilde{q}^2}{\tilde{q}^2}
	* \int\frac{\alpha_S(z,\tilde{q})}{2\pi}
	* P_{ba}(z,\tilde{q})\Theta(p_T)
	* \right\}.
	* \f]
	* We can solve this to obtain the next value of the scale \f$\tilde{q}_{i+1}\f$
	* given the previous value \f$\tilde{q}_i\f$
	* in the following way. First we obtain a simplified form of the integrand
	* which is greater than or equal to the true integrand for all values of
	* \f$\tilde{q}\f$.
	*
	* In practice it is easiest to obtain this over estimate in pieces. The ShowerAlpha
	* object contains an over estimate for \f$\alpha_S\f$, the splitting function
	* contains both an over estimate of the spltting function and its integral
	* which is needed to compute the over estimate of the \f$\tilde{q}\f$ integrand,
	* together with an over estimate of the limit of the \f$z\f$ integral.
	*
	* This gives an overestimate of the integrand
	* \f[g(\tilde{q}^2) = \frac{c}{\tilde{q}^2}, \f]
	* where because the over estimates are chosen to be independent of \f$\tilde{q}\f$ the
	* parameter
	* \f[c = \frac{\alpha_{\rm over}}{2\pi}\int^{z_1}_{z_0}P_{\rm over}(z),\f]
	* is a constant independent of \f$\tilde{q}\f$.
	*
	* The guesst() member can then be used to generate generate the value of
	* \f$\tilde{q}^2\f$ according to this result. This is done by solving the Sudakov
	* form factor, with the over estimates, is equal to a random number
	* \f$r\f$ in the interval \f$[0,1]\f$. This gives
	* \f[\tilde{q}^2_{i+1}=G^{-1}\left[G(\tilde{q}^2_i)+\ln r\right],\f]
	* where \f$G(\tilde{q}^2)=c\ln(\tilde{q}^2)\f$ is the infinite integral
	* of \f$g(\tilde{q}^2)\f$ and \f$G^{-1}(x)=\exp\left(\frac{x}c\right)\f$
	* is its inverse.
	* It this case we therefore obtain
	* \f[\tilde{q}^2_{i+1}=\tilde{q}^2_ir^{\frac1c}.\f]
	* The value of \f$z\f$ can then be calculated in a similar way
	* \f[z = I^{-1}\left[I(z_0)+r\left(I(z_1)-I(z_0)\right)\right],\f]
	* using the guessz() member,
	* where \f$I=\int P(z){\rm d}z\f$ and \f$I^{-1}\f$ is its inverse.
	*
	* The veto algorithm then uses rejection using the ratio of the
	* true value to the overestimated one to obtain the original distribution.
	* This is accomplished using the
	* - alphaSVeto() member for the \f$\alpha_S\f$ veto
	* - SplittingFnVeto() member for the veto on the value of the splitting function.
	* in general there must also be a chech that the emission is in the allowed
	* phase space but this is left to the inheriting classes as it will depend
	* on the ordering variable.
	*
	* The Sudakov form factor for the initial-scale shower is different because
	* it must include the PDF which guides the backward evolution.
	* It is given by
	* \f[\Delta_{ba}(\tilde{q}_{i+1},\tilde{q}_i)=
	* \exp\left\{
	* -\int^{\tilde{q}^2_i}_{\tilde{q}^2_{i+1}}
	* \frac{{\rm d}\tilde{q}^2}{\tilde{q}^2}
	* \int\frac{\alpha_S(z,\tilde{q})}{2\pi}
	* P_{ba}(z,\tilde{q})\frac{x'f_a(\frac{x}z,\tilde{q}^2)}{xf_b(x,\tilde{q^2})}
	* \right\},
	* \f]
	* where \f$x\f$ is the fraction of the beam momentum the parton \f$b\f$ had before
	* the backward evolution.
	* This can be solve in the same way as for the final-state branching but the constant
	* becomes
	* \f[c = \frac{\alpha_{\rm over}}{2\pi}\int^{z_1}_{z_0}P_{\rm over}(z)PDF_{\rm max},\f]
	* where
	* \f[PDF_{\rm max}=\max\frac{x'f_a(\frac{x}z,\tilde{q}^2)}{xf_b(x,\tilde{q^2})},\f]
	* which can be set using an interface.
	* In addition the PDFVeto() member then is needed to implement the relevant veto.
	*
	* @see SplittingGenerator
	* @see SplittingFunction
	* @see ShowerAlpha
	* @see \ref SudakovFormFactorInterfaces "The interfaces"
	* defined for SudakovFormFactor.
	*/
	class SudakovFormFactor: public Interfaced {

	/**
	* The SplittingGenerator is a friend to insert the particles in the
	* branchings at initialisation
	*/
	friend class SplittingGenerator;

	public:

	/**
	* The default constructor.
	*/
	SudakovFormFactor() : pdfmax_(35.0), pdffactor_(0),
	cutOffOption_(0), a_(0.3), b_(2.3), c_(0.3*GeV),
	kinCutoffScale_( 2.3GeV ), vgcut_(0.85GeV),
	vqcut_(0.85GeV), pTmin_(1.GeV), pT2min_(ZERO),
	z_( 0.0 ),phi_(0.0), pT_(),
	theFactorizationScaleFactor(1.0),
	theRenormalizationScaleFactor(1.0) {}

	/**
	* Members to generate the scale of the next branching
	*/
	//@{
	/**
	* Return the scale of the next time-like branching. If there is no
	* branching then it returns ZERO.
	* @param startingScale starting scale for the evolution
	* @param ids The PDG codes of the particles in the splitting
	* @param cc Whether this is the charge conjugate of the branching
	* @param enhance The radiation enhancement factor
	* @param maxQ2 The maximum \f$Q^2\f$ for the emission
	*/
	virtual ShoKinPtr generateNextTimeBranching(const Energy startingScale,
	const IdList &ids,const bool cc,
	double enhance, Energy2 maxQ2)=0;

	/**
	* Return the scale of the next space-like decay branching. If there is no
	* branching then it returns ZERO.
	* @param startingScale starting scale for the evolution
	* @param stoppingScale stopping scale for the evolution
	* @param minmass The minimum mass allowed for the spake-like particle.
	* @param ids The PDG codes of the particles in the splitting
	* @param cc Whether this is the charge conjugate of the branching
	* defined.
	* @param enhance The radiation enhancement factor
	*/
	virtual ShoKinPtr generateNextDecayBranching(const Energy startingScale,
	const Energy stoppingScale,
	const Energy minmass,
	const IdList &ids,
	const bool cc,
	double enhance)=0;

	/**
	* Return the scale of the next space-like branching. If there is no
	* branching then it returns ZERO.
	* @param startingScale starting scale for the evolution
	* @param ids The PDG codes of the particles in the splitting
	* @param x The fraction of the beam momentum
	* @param cc Whether this is the charge conjugate of the branching
	* defined.
	* @param beam The beam particle
	* @param enhance The radiation enhancement factor
	*/
	virtual ShoKinPtr generateNextSpaceBranching(const Energy startingScale,
	const IdList &ids,double x,
	const bool cc,double enhance,
	tcBeamPtr beam)=0;
	//@}

	/**
	* Generate the azimuthal angle of the branching for forward evolution
	* @param particle The branching particle
	* @param ids The PDG codes of the particles in the branchings
	* @param The Shower kinematics
	*/
	virtual double generatePhiForward(ShowerParticle & particle,const IdList & ids,
	ShoKinPtr kinematics)=0;

	/**
	* Generate the azimuthal angle of the branching for backward evolution
	* @param particle The branching particle
	* @param ids The PDG codes of the particles in the branchings
	* @param The Shower kinematics
	*/
	virtual double generatePhiBackward(ShowerParticle & particle,const IdList & ids,
	ShoKinPtr kinematics)=0;

	/**
	* Generate the azimuthal angle of the branching for ISR in decays
	* @param particle The branching particle
	* @param ids The PDG codes of the particles in the branchings
	* @param The Shower kinematics
	*/
	virtual double generatePhiDecay(ShowerParticle & particle,const IdList & ids,
	ShoKinPtr kinematics)=0;

	/**
	* Methods to provide public access to the private member variables
	*/
	//@{
	/**
	* Return the pointer to the SplittingFunction object.
	*/
	tSplittingFnPtr splittingFn() const { return splittingFn_; }

	/**
	* Return the pointer to the ShowerAlpha object.
	*/
	tShowerAlphaPtr alpha() const { return alpha_; }

	/**
	* The type of interaction
	*/
	inline ShowerInteraction::Type interactionType() const
	{return splittingFn_->interactionType();}
	//@}

	public:

	/**
	* Methods to access the kinematic variables for the branching
	*/
	//@{
	/**
	* The energy fraction
	*/
	double z() const { return z_; }

	/**
	* The azimuthal angle
	*/
	double phi() const { return phi_; }

	/**
	* The transverse momentum
	*/
	Energy pT() const { return pT_; }
	//@}

	/**
	* Access the maximum weight for the PDF veto
	*/
	double pdfMax() const { return pdfmax_;}

	/**
	* Method to return the evolution scale given the
	* transverse momentum, \f$p_T\f$ and \f$z\f$.
	*/
	virtual Energy calculateScale(double z, Energy pt, IdList ids,unsigned int iopt)=0;

	/**
	* Method to create the ShowerKinematics object for a final-state branching
	*/
	virtual ShoKinPtr createFinalStateBranching(Energy scale,double z,
	double phi, Energy pt)=0;

	/**
	* Method to create the ShowerKinematics object for an initial-state branching
	*/
	virtual ShoKinPtr createInitialStateBranching(Energy scale,double z,
	double phi, Energy pt)=0;

	/**
	* Method to create the ShowerKinematics object for a decay branching
	*/
	virtual ShoKinPtr createDecayBranching(Energy scale,double z,
	double phi, Energy pt)=0;

	public:

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

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

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

	protected:

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

	protected:

	/**
	* Methods to implement the veto algorithm to generate the scale of
	* the next branching
	*/
	//@{
	/**
	* Value of the energy fraction for the veto algorithm
	* @param iopt The option for calculating z
	* @param ids The PDG codes of the particles in the splitting
	* - 0 is final-state
	* - 1 is initial-state for the hard process
	* - 2 is initial-state for particle decays
	*/
	double guessz (unsigned int iopt, const IdList &ids) const;

	/**
	* Value of the scale for the veto algorithm
	* @param t1 The starting valoe of the scale
	* @param iopt The option for calculating t
	* @param ids The PDG codes of the particles in the splitting
	* - 0 is final-state
	* - 1 is initial-state for the hard process
	* - 2 is initial-state for particle decays
	* @param enhance The radiation enhancement factor
	* @param identical Whether or not the outgoing particles are identical
	*/
	Energy2 guesst (Energy2 t1,unsigned int iopt, const IdList &ids,
	double enhance, bool identical) const;

	/**
	* Veto on the PDF for the initial-state shower
	* @param t The scale
	* @param x The fraction of the beam momentum
	* @param parton0 Pointer to the particleData for the
	* new parent (this is the particle we evolved back to)
	* @param parton1 Pointer to the particleData for the
	* original particle
	* @param beam The BeamParticleData object
	*/
	bool PDFVeto(const Energy2 t, const double x,
	const tcPDPtr parton0, const tcPDPtr parton1,
	tcBeamPtr beam) const;

	/**
	* The veto on the splitting function.
	* @param t The scale
	* @param ids The PDG codes of the particles in the splitting
	* @param mass Whether or not to use the massive splitting functions
	* @return true if vetoed
	*/
	bool SplittingFnVeto(const Energy2 t,
	const IdList &ids,
	const bool mass) const {
	return UseRandom::rnd()>splittingFn_->ratioP(z_, t, ids,mass);
	}

	/**
	* The veto on the coupling constant
	* @param pt2 The value of ther transverse momentum squared, \f$p_T^2\f$.
	* @return true if vetoed
	*/
	bool alphaSVeto(Energy2 pt2) const;
	//@}

	/**
	* Methods to set the kinematic variables for the branching
	*/
	//@{
	/**
	* The energy fraction
	*/
	void z(double in) { z_=in; }

	/**
	* The azimuthal angle
	*/
	void phi(double in) { phi_=in; }

	/**
	* The transverse momentum
	*/
	void pT(Energy in) { pT_=in; }
	//@}

	/**
	* Set/Get the limits on the energy fraction for the splitting
	*/
	//@{
	/**
	* Get the limits
	*/
	pair<double,double> zLimits() const { return zlimits_;}

	/**
	* Set the limits
	*/
	void zLimits(pair<double,double> in) { zlimits_=in; }
	//@}

	/**
	* Set the particles in the splittings
	*/
	void addSplitting(const IdList &);

	/**
	* Delete the particles in the splittings
	*/
	void removeSplitting(const IdList &);

	/**
	* Access the potential branchings
	*/
	const vector<IdList> & particles() const { return particles_; }

	/**
	* For a particle which came from the hard process get the spin density and
	* the mapping required to the basis used in the Shower
	* @param rho The \f$\rho\f$ matrix
	* @param mapping The mapping
	* @param particle The particle
	* @param showerkin The ShowerKinematics object
	*/
	bool getMapping(SpinPtr &, RhoDMatrix & map,
	ShowerParticle & particle,ShoKinPtr showerkin);

	+ RhoDMatrix extractRhoMatrix(ShowerParticle & particle,
	+ ShoKinPtr kinematics, bool forward);
	+
	public:

	/**
	* @name Methods for the cut-off
	*/
	//@{
	/**
	* The option being used
	*/
	unsigned int cutOffOption() const { return cutOffOption_; }

	/**
	* The kinematic scale
	*/
	Energy kinScale() const {return kinCutoffScale_;}

	/**
	* The virtuality cut-off on the gluon \f$Q_g=\frac{\delta-am_q}{b}\f$
	* @param scale The scale \f$\delta\f$
	* @param mq The quark mass \f$m_q\f$.
	*/
	Energy kinematicCutOff(Energy scale, Energy mq) const
	{return max((scale -a_*mq)/b_,c_);}

	/**
	* The virtualilty cut-off for gluons
	*/
	Energy vgCut() const { return vgcut_; }

	/**
	* The virtuality cut-off for everything else
	*/
	Energy vqCut() const { return vqcut_; }

	/**
	* The minimum \f$p_T\f$ for the branching
	*/
	Energy pTmin() const { return pTmin_; }

	/**
	* The square of the minimum \f$p_T\f$
	*/
	Energy2 pT2min() const { return pT2min_; }

	/**
	* Calculate the virtual masses for a branchings
	*/
	const vector<Energy> & virtualMasses(const IdList & ids);
	//@}

	/**
	* Set the PDF
	*/
	void setPDF(tcPDFPtr pdf, Energy scale) {
	pdf_ = pdf;
	freeze_ = scale;
	}

	private:

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

	private:

	/**
	* Pointer to the splitting function for this Sudakov form factor
	*/
	SplittingFnPtr splittingFn_;

	/**
	* Pointer to the coupling for this Sudakov form factor
	*/
	ShowerAlphaPtr alpha_;

	/**
	* Maximum value of the PDF weight
	*/
	double pdfmax_;

	/**
	* List of the particles this Sudakov is used for to aid in setting up
	* interpolation tables if needed
	*/
	vector<IdList> particles_;

	/**
	* Option for the inclusion of a factor \f$1/(1-z)\f$ in the PDF estimate
	*/
	unsigned pdffactor_;

	private:

	/**
	* Option for the type of cut-off to be applied
	*/
	unsigned int cutOffOption_;

	/**
	* Parameters for the default Herwig cut-off option, i.e. the parameters for
	* the \f$Q_g=\max(\frac{\delta-am_q}{b},c)\f$ kinematic cut-off
	*/
	//@{
	/**
	* The \f$a\f$ parameter
	*/
	double a_;

	/**
	* The \f$b\f$ parameter
	*/
	double b_;

	/**
	* The \f$c\f$ parameter
	*/
	Energy c_;

	/**
	* Kinematic cutoff used in the parton shower phase space.
	*/
	Energy kinCutoffScale_;
	//@}

	/**
	* Parameters for the FORTRAN-like cut-off
	*/
	//@{
	/**
	* The virtualilty cut-off for gluons
	*/
	Energy vgcut_;

	/**
	* The virtuality cut-off for everything else
	*/
	Energy vqcut_;
	//@}

	/**
	* Parameters for the \f$p_T\f$ cut-off
	*/
	//@{
	/**
	* The minimum \f$p_T\f$ for the branching
	*/
	Energy pTmin_;

	/**
	* The square of the minimum \f$p_T\f$
	*/
	Energy2 pT2min_;
	//@}

	private:

	/**
	* Member variables to keep the shower kinematics information
	* generated by a call to generateNextTimeBranching or generateNextSpaceBranching
	*/
	//@{
	/**
	* The energy fraction
	*/
	double z_;

	/**
	* The azimuthal angle
	*/
	double phi_;

	/**
	* The transverse momentum
	*/
	Energy pT_;
	//@}

	/**
	* The limits of \f$z\f$ in the splitting
	*/
	pair<double,double> zlimits_;

	/**
	* Stuff for the PDFs
	*/
	//@{
	/**
	* PDf
	*/
	tcPDFPtr pdf_;

	/**
	* Freezing scale
	*/
	Energy freeze_;
	//@}

	public:

	/**
	* Get the factorization scale factor
	*/
	double factorizationScaleFactor() const { return theFactorizationScaleFactor; }

	/**
	* Set the factorization scale factor
	*/
	void factorizationScaleFactor(double f) { theFactorizationScaleFactor = f; }

	/**
	* Get the renormalization scale factor
	*/
	double renormalizationScaleFactor() const { return theRenormalizationScaleFactor; }

	/**
	* Set the renormalization scale factor
	*/
	void renormalizationScaleFactor(double f) { theRenormalizationScaleFactor = f; }

	private:

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

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

	};

	}

	#endif /* HERWIG_SudakovFormFactor_H */
	diff --git a/Shower/Default/QTildeSudakov.cc b/Shower/Default/QTildeSudakov.cc
	--- a/Shower/Default/QTildeSudakov.cc
	+++ b/Shower/Default/QTildeSudakov.cc
	@@ -1,981 +1,939 @@
	// -- 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/Repository/CurrentGenerator.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,
	const Energy2 maxQ2) {
	// still inside PS, return true if outside
	// check vs overestimated limits
	if(z() < zLimits().first \|\| z() > zLimits().second) return true;
	Energy2 q2 = z()(1.-z())t;
	if(ids_[0]!=ParticleID::g &&
	ids_[0]!=ParticleID::gamma ) q2 += masssquared_[0];
	if(q2>maxQ2) return true;
	// compute the pts
	Energy2 pt2 = z()(1.-z())q2 - masssquared_[1](1.-z()) - masssquared_[2]z();
	// 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, Energy2 maxQ2) {
	// 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,maxQ2) \|\| SplittingFnVeto(z()(1.-z())t,ids,true) \|\|
	alphaSVeto(splittingFn()->angularOrdered() ? sqr(z()(1.-z()))t : 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(splittingFn()->angularOrdered() ? sqr(1.-z())t : (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(splittingFn()->angularOrdered() ? sqr(1.-z())t : (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]+pT2min())) {
	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;
	tShowerParticlePtr mother;
	if(forward) {
	mother = !particle.parents().empty() ?
	dynamic_ptr_cast<tShowerParticlePtr>(particle.parents()[0]) : tShowerParticlePtr();
	}
	else {
	mother = particle.children().size()==2 ?
	dynamic_ptr_cast<tShowerParticlePtr>(&particle) : tShowerParticlePtr();
	}
	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() ?
	dynamic_ptr_cast<tShowerParticlePtr>(mother->parents()[0]) : tShowerParticlePtr();
	}
	else {
	if(mother->children()[0]->children().size()!=2)
	break;
	tShowerParticlePtr mtemp =
	dynamic_ptr_cast<tShowerParticlePtr>(mother->children()[0]);
	if(!mtemp)
	break;
	else
	mother=mtemp;
	}
	}
	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;
	}

	pair<double,double> softPhiMin(double phi0, double phi1, double A, double B, double C, double D) {
	double c01 = cos(phi0 - phi1);
	double s01 = sin(phi0 - phi1);
	double s012(sqr(s01)), c012(sqr(c01));
	double A2(AA), B2(BB), C2(CC), D2(DD);
	if(abs(B/A)<1e-10 && abs(D/C)<1e-10) return make_pair(phi0,phi0+Constants::pi);
	double root = sqr(B2)C2D2sqr(s012) + 2.AB2BC2CDc01s012 + 2.AB2BCD2Dc01*s012
	+ 4.A2B2C2D2c012 - A2B2C2D2s012 - A2B2sqr(D2)s012 - sqr(B2)sqr(C2)s012
	- sqr(B2)C2D2s012 - 4.A2ABCD2Dc01 - 4.AB2BC2CDc01 + sqr(A2)sqr(D2)
	+ 2.A2B2C2D2 + sqr(B2)*sqr(C2);
	if(root<0.) return make_pair(phi0,phi0+Constants::pi);
	root = sqrt(root);
	double denom = (-2.ABCDc01 + A2D2 + B2*C2);
	double denom2 = (-BCc01 + A*D);
	double num = B2CD*s012;
	return make_pair(atan2(Bs01(-C*(num + root) / denom + D) / denom2, -(num + root ) / denom) + phi0,
	atan2(Bs01(-C*(num - root) / denom + D) / denom2, -(num - root ) / denom) + phi0);
	}

	}

	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 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);
	vector<Energy> Ek(3,ZERO);
	Energy Ei,Ej;
	Energy2 m12(ZERO),m22(ZERO);
	InvEnergy2 aziMax(ZERO);
	bool softAllowed = ShowerHandler::currentHandler()->evolver()->softCorrelations()&&
	(canBeSoft[0] \|\| canBeSoft[1]);
	if(softAllowed) {
	// 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;
	// compute the transforms to the shower reference frame
	// first the boost
	vector<Lorentz5Momentum> basis = kinematics->getBasis();
	Lorentz5Momentum pVect = basis[0], 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, 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.;
	}
	Ek[0] = zFact(alpha0pVect.t()-0.5nVect.t()/pn(alpha0*m2-sqr(particle.showerParameters().pt)/alpha0))
	+0.5sqr(pT)nVect.t()/pn/zFact/alpha0;
	Ek[1] = -nVect.t()/alpha0/pnqperp0.x()pT;
	Ek[2] = -nVect.t()/alpha0/pnqperp0.y()pT;
	if(swapOrder) {
	Ek[1] *= -1.;
	Ek[2] *= -1.;
	}
	Energy mag2=sqrt(sqr(Ek[1])+sqr(Ek[2]));
	Ei = alpha0pVect.t()+beta0nVect.t();
	Ej = pj.t();
	double phi0 = atan2(-pjk[2],-pjk[1]);
	if(phi0<0.) phi0 += Constants::twopi;
	double phi1 = atan2(-Ek[2],-Ek[1]);
	if(phi1<0.) phi1 += Constants::twopi;
	double xi_min = pik/Ei/(Ek[0]+mag2), xi_max = pik/Ei/(Ek[0]-mag2), xi_ij = pipj/Ei/Ej;
	if(xi_min>xi_max) swap(xi_min,xi_max);
	if(xi_min>xi_ij) softAllowed = false;
	Energy2 mag = sqrt(sqr(pjk[1])+sqr(pjk[2]));
	if(ShowerHandler::currentHandler()->evolver()->softCorrelations()==1) {
	aziMax = -m12/sqr(pik) -m22/sqr(pjk[0]+mag) +2.*pipj/pik/(pjk[0]-mag);
	}
	else if(ShowerHandler::currentHandler()->evolver()->softCorrelations()==2) {
	double A = (pipjEk[0]- Ejpik)/Ej/sqr(Ej);
	double B = -sqrt(sqr(pipj)*(sqr(Ek[1])+sqr(Ek[2])))/Ej/sqr(Ej);
	double C = pjk[0]/sqr(Ej);
	double D = -sqrt(sqr(pjk[1])+sqr(pjk[2]))/sqr(Ej);
	pair<double,double> minima = softPhiMin(phi0,phi1,A,B,C,D);
	aziMax = 0.5/pik/(Ek[0]-mag2)(Ei-m12(Ek[0]-mag2)/pik + max(Ej(A+Bcos(minima.first -phi1))/(C+D*cos(minima.first -phi0)),
	Ej(A+Bcos(minima.second-phi1))/(C+D*cos(minima.second-phi0))));
	}
	else
	assert(false);
	}
	// 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;
	- }
	+ RhoDMatrix rho = extractRhoMatrix(particle,kinematics,true);
	// calculate the weights
	wgts = splittingFn()->generatePhiForward(z,t,ids,rho);
	}
	else {
	wgts = vector<pair<int,Complex> >(1,make_pair(0,1.));
	}
	// generate the azimuthal angle
	double phi,wgt;
	static const Complex ii(0.,1.);
	unsigned int ntry(0);
	double phiMax(0.),wgtMax(0.);
	do {
	phi = Constants::twopi*UseRandom::rnd();
	// 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)
	spinWgt += wgts[ix].second;
	else
	spinWgt += exp(double(wgts[ix].first)iiphi)*wgts[ix].second;
	}
	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)
	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);
	Energy Eg = Ek[0]+Ek[1]cos(phi)+Ek[2]sin(phi);
	if(pipjEg>pikEj) {
	if(ShowerHandler::currentHandler()->evolver()->softCorrelations()==1) {
	aziWgt = (-m12/sqr(pik) -m22/sqr(dot) +2.*pipj/pik/dot)/aziMax;
	}
	else if(ShowerHandler::currentHandler()->evolver()->softCorrelations()==2) {
	aziWgt = max(ZERO,0.5/pik/Eg(Ei-m12Eg/pik + (pipjEg - Ejpik)/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";
	}
	}
	else {
	aziWgt = 0.;
	}
	}
	wgt *= aziWgt;
	if(wgt>wgtMax) {
	phiMax = phi;
	wgtMax = wgt;
	}
	++ntry;
	}
	while(wgt<UseRandom::rnd()&&ntry<10000);
	if(ntry==10000) {
	generator()->log() << "Too many tries to generate phi in forward evolution\n";
	phi = phiMax;
	}
	// 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 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);
	Energy Ei,Ej,Ek;
	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]));
	if(ShowerHandler::currentHandler()->evolver()->softCorrelations()==1) {
	aziMax = -m12/sqr(pik) -m22/sqr(pjk[0]+mag) +2.*pipj/pik/(pjk[0]-mag);
	}
	else if(ShowerHandler::currentHandler()->evolver()->softCorrelations()==2) {
	Ek = alpha0zFact/(1.-zFact)pVect.t()+beta2*nVect.t();
	Ei = alpha0pVect.t()+beta0nVect.t();
	Ej = pj.t();
	if(pipjEk> Ejpik) {
	aziMax = 0.5/pik/Ek(Ei-m12Ek/pik + (pipjEk- Ejpik)/(pjk[0]-mag));
	}
	else {
	aziMax = 0.5/pik/Ek(Ei-m12Ek/pik);
	}
	}
	else {
	assert(ShowerHandler::currentHandler()->evolver()->softCorrelations()==0);
	}
	}
	// 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;
	- }
	+ RhoDMatrix rho = extractRhoMatrix(particle,kinematics,false);
	+ // get the weights
	wgts = splittingFn()->generatePhiBackward(z,t,ids,rho);
	}
	else {
	wgts = vector<pair<int,Complex> >(1,make_pair(0,1.));
	}
	// generate the azimuthal angle
	double phi,wgt;
	static const Complex ii(0.,1.);
	unsigned int ntry(0);
	double phiMax(0.),wgtMax(0.);
	do {
	phi = Constants::twopi*UseRandom::rnd();
	Complex spinWgt = 0.;
	for(unsigned int ix=0;ix<wgts.size();++ix) {
	if(wgts[ix].first==0)
	spinWgt += wgts[ix].second;
	else
	spinWgt += exp(double(wgts[ix].first)iiphi)*wgts[ix].second;
	}
	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";
	generator()->log() << "Weights \n";
	for(unsigned int ix=0;ix<wgts.size();++ix)
	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);
	if(ShowerHandler::currentHandler()->evolver()->softCorrelations()==1) {
	aziWgt = (-m12/sqr(pik) -m22/sqr(dot) +2.*pipj/pik/dot)/aziMax;
	}
	else if(ShowerHandler::currentHandler()->evolver()->softCorrelations()==2) {
	aziWgt = max(ZERO,0.5/pik/Ek(Ei-m12Ek/pik + pipjEk/dot - Ejpik/dot)/aziMax);
	}
	if(aziWgt-1.>1e-10\|\|aziWgt<-1e-10) {
	generator()->log() << "Backward 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<UseRandom::rnd()&&ntry<10000);
	if(ntry==10000) {
	generator()->log() << "Too many tries to generate phi in backward evolution\n";
	phi = phiMax;
	}
	// return the azimuthal angle
	return phi;
	}

	double QTildeSudakov::generatePhiDecay(ShowerParticle & particle,
	const IdList & ids,
	ShoKinPtr kinematics) {
	// only soft correlations in this case
	// no correlations, return flat phi
	if( !(ShowerHandler::currentHandler()->evolver()->softCorrelations() &&
	(ids[2]==ParticleID::g \|\| ids[2]==ParticleID::gamma )))
	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());
	double zFact(1.-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());
	Energy2 mag = sqrt(sqr(pjk[1])+sqr(pjk[2]));
	InvEnergy2 aziMax;
	vector<Energy> Ek(3,ZERO);
	Energy Ei,Ej;
	if(ShowerHandler::currentHandler()->evolver()->softCorrelations()==1) {
	aziMax = -m12/sqr(pik) -m22/sqr(pjk[0]+mag) +2.*pipj/pik/(pjk[0]-mag);
	}
	else if(ShowerHandler::currentHandler()->evolver()->softCorrelations()==2) {
	Ek[0] = zFact(alpha0pVect.t()+-0.5nVect.t()/pn(alpha0*m2-sqr(particle.showerParameters().pt)/alpha0))
	+0.5sqr(pT)nVect.t()/pn/zFact/alpha0;
	Ek[1] = -nVect.t()/alpha0/pnqperp0.x()pT;
	Ek[2] = -nVect.t()/alpha0/pnqperp0.y()pT;
	Energy mag2=sqrt(sqr(Ek[1])+sqr(Ek[2]));
	Ei = alpha0pVect.t()+beta0nVect.t();
	Ej = pj.t();
	aziMax = 0.5/pik/(Ek[0]-mag2)(Ei-m12(Ek[0]-mag2)/pik + pipj(Ek[0]+mag2)/(pjk[0]-mag) - Ejpik/(pjk[0]-mag) );
	}
	else
	assert(ShowerHandler::currentHandler()->evolver()->softCorrelations()==0);
	// generate the azimuthal angle
	double phi,wgt(0.);
	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);
	if(ShowerHandler::currentHandler()->evolver()->softCorrelations()==1) {
	wgt = (-m12/sqr(pik) -m22/sqr(dot) +2.*pipj/pik/dot)/aziMax;
	}
	else if(ShowerHandler::currentHandler()->evolver()->softCorrelations()==2) {
	if(qperp0.m2()==ZERO) {
	wgt = 1.;
	}
	else {
	Energy Eg = Ek[0]+Ek[1]cos(phi)+Ek[2]sin(phi);
	wgt = max(ZERO,0.5/pik/Eg(Ei-m12Eg/pik + (pipjEg - Ejpik)/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;
	}

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 8:31 PM (1 d, 5 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3806046
Default Alt Text: (77 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions