No OneTemporary
Actions

Size

67 KB

Subscribers

None

View Options

	diff --git a/Shower/QTilde/Base/PartnerFinder.cc b/Shower/QTilde/Base/PartnerFinder.cc
	--- a/Shower/QTilde/Base/PartnerFinder.cc
	+++ b/Shower/QTilde/Base/PartnerFinder.cc
	@@ -1,699 +1,699 @@
	// -- C++ --
	//
	// PartnerFinder.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2019 The Herwig Collaboration
	//
	// Herwig is licenced under version 3 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 PartnerFinder class.
	//

	#include "PartnerFinder.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "Herwig/Shower/QTilde/Base/ShowerParticle.h"
	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Utilities/Debug.h"
	#include "ThePEG/Utilities/DescribeClass.h"

	using namespace Herwig;

	DescribeClass<PartnerFinder,Interfaced>
	describePartnerFinder ("Herwig::PartnerFinder","HwShower.so");

	// some useful functions to avoid using #define
	namespace {

	// return bool if final-state particle
	inline bool FS(const tShowerParticlePtr a) {
	return a->isFinalState();
	}

	// return colour line pointer
	inline Ptr<ThePEG::ColourLine>::transient_pointer
	CL(const tShowerParticlePtr a, unsigned int index=0) {
	return a->colourInfo()->colourLines().empty() ? ThePEG::tColinePtr() :
	const_ptr_cast<ThePEG::tColinePtr>(a->colourInfo()->colourLines()[index]);
	}

	// return colour line size
	inline size_t
	CLSIZE(const tShowerParticlePtr a) {
	return a->colourInfo()->colourLines().size();
	}

	inline Ptr<ThePEG::ColourLine>::transient_pointer
	ACL(const tShowerParticlePtr a, unsigned int index=0) {
	return a->colourInfo()->antiColourLines().empty() ? ThePEG::tColinePtr() :
	const_ptr_cast<ThePEG::tColinePtr>(a->colourInfo()->antiColourLines()[index]);
	}

	inline size_t
	ACLSIZE(const tShowerParticlePtr a) {
	return a->colourInfo()->antiColourLines().size();
	}
	}

	void PartnerFinder::persistentOutput(PersistentOStream & os) const {
	os << partnerMethod_ << QEDPartner_ << scaleChoice_;
	}

	void PartnerFinder::persistentInput(PersistentIStream & is, int) {
	is >> partnerMethod_ >> QEDPartner_ >> scaleChoice_;
	}

	void PartnerFinder::Init() {

	static ClassDocumentation<PartnerFinder> documentation
	("This class is responsible for finding the partners for each interaction types ",
	"and within the evolution scale range specified by the ShowerVariables ",
	"then to determine the initial evolution scales for each pair of partners.");

	static Switch<PartnerFinder,int> interfacePartnerMethod
	("PartnerMethod",
	"Choice of partner finding method for gluon evolution.",
	&PartnerFinder::partnerMethod_, 0, false, false);
	static SwitchOption interfacePartnerMethodRandom
	(interfacePartnerMethod,
	"Random",
	"Choose partners of a gluon randomly.",
	0);
	static SwitchOption interfacePartnerMethodMaximum
	(interfacePartnerMethod,
	"Maximum",
	"Choose partner of gluon with largest angle.",
	1);

	static Switch<PartnerFinder,int> interfaceQEDPartner
	("QEDPartner",
	"Control of which particles to use as the partner for QED radiation",
	&PartnerFinder::QEDPartner_, 0, false, false);
	static SwitchOption interfaceQEDPartnerAll
	(interfaceQEDPartner,
	"All",
	"Consider all possible choices which give a positive contribution"
	" in the soft limit.",
	0);
	static SwitchOption interfaceQEDPartnerIIandFF
	(interfaceQEDPartner,
	"IIandFF",
	"Only allow initial-initial or final-final combinations",
	1);
	static SwitchOption interfaceQEDPartnerIF
	(interfaceQEDPartner,
	"IF",
	"Only allow initial-final combinations",
	2);

	static Switch<PartnerFinder,int> interfaceScaleChoice
	("ScaleChoice",
	"The choice of the evolution scales",
	&PartnerFinder::scaleChoice_, 0, false, false);
	static SwitchOption interfaceScaleChoicePartner
	(interfaceScaleChoice,
	"Partner",
	"Scale of all interactions is that of the evolution partner",
	0);
	static SwitchOption interfaceScaleChoiceDifferent
	(interfaceScaleChoice,
	"Different",
	"Allow each interaction to have different scales",
	1);
	}

	void PartnerFinder::setInitialEvolutionScales(const ShowerParticleVector &particles,
	const bool isDecayCase,
	ShowerInteraction type,
	const bool setPartners) {
	// clear the existing partners
	for(ShowerParticleVector::const_iterator cit = particles.begin();
	cit != particles.end(); ++cit) (*cit)->clearPartners();
	// set them
	if(type==ShowerInteraction::QCD) {
	setInitialQCDEvolutionScales(particles,isDecayCase,setPartners);
	}
	else if(type==ShowerInteraction::QED) {
	setInitialQEDEvolutionScales(particles,isDecayCase,setPartners);
	}
	else if(type==ShowerInteraction::EW) {
	setInitialEWEvolutionScales(particles,isDecayCase,false);
	}
	else if(type==ShowerInteraction::QEDQCD) {
	setInitialQCDEvolutionScales(particles,isDecayCase,setPartners);
	setInitialQEDEvolutionScales(particles,isDecayCase,false);
	}
	else if(type==ShowerInteraction::ALL) {
	setInitialQCDEvolutionScales(particles,isDecayCase,setPartners);
	setInitialQEDEvolutionScales(particles,isDecayCase,false);
	setInitialEWEvolutionScales(particles,isDecayCase,false);
	}
	else
	assert(false);
	// \todo EW scales here
	// print out for debugging
	if(Debug::level>=10) {
	for(ShowerParticleVector::const_iterator cit = particles.begin();
	cit != particles.end(); ++cit) {
	generator()->log() << "Particle: " << **cit << "\n";
	if(!(**cit).partner()) continue;
	generator()->log() << "Primary partner: " << (*cit).partner() << "\n";
	for(vector<ShowerParticle::EvolutionPartner>::const_iterator it= (**cit).partners().begin();
	it!=(**cit).partners().end();++it) {
	generator()->log() << static_cast<long>(it->type) << " "
	<< it->weight << " "
	<< it->scale/GeV << " "
	<< *(it->partner)
	<< "\n";
	}
	}
	generator()->log() << flush;
	}
	}

	void PartnerFinder::setInitialQCDEvolutionScales(const ShowerParticleVector &particles,
	const bool isDecayCase,
	const bool setPartners) {
	// Loop over particles and consider only coloured particles which don't
	// have already their colour partner fixed and that don't have children
	// (the latter requirement is relaxed in the case isDecayCase is true).
	// Build a map which has as key one of these particles (i.e. a pointer
	// to a ShowerParticle object) and as a corresponding value the vector
	// of all its possible normal candidate colour partners, defined as follows:
	// --- have colour, and no children (this is not required in the case
	// isDecayCase is true);
	// --- if both are initial (incoming) state particles, then the (non-null) colourLine()
	// of one of them must match the (non-null) antiColourLine() of the other.
	// --- if one is an initial (incoming) state particle and the other is
	// a final (outgoing) state particle, then both must have the
	// same (non-null) colourLine() or the same (non-null) antiColourLine();
	// Notice that this definition exclude the special case of baryon-violating
	// processes (as in R-parity Susy), which will show up as particles
	// without candidate colour partners, and that we will be treated a part later
	// (this means that no modifications of the following loop is needed!

	for ( const auto & sp : particles ) {
	// Skip colourless particles
	if(!sp->data().coloured()) continue;
	// find the partners
	auto partners = findQCDPartners(sp,particles);
	// must have a partner
	if(partners.empty()) {
	throw Exception() << "`Failed to make colour connections in "
	<< "PartnerFinder::setQCDInitialEvolutionScales"
	<< *sp
	<< Exception::eventerror;
	}
	// Calculate the evolution scales for all possible pairs of of particles
	vector<pair<Energy,Energy> > scales;
	int position = -1;
	for(size_t ix=0; ix< partners.size(); ++ix) {
	scales.push_back(calculateInitialEvolutionScales(ShowerPPair(sp, partners[ix].second),isDecayCase));
	if (!setPartners && partners[ix].second) position = ix;
	}
	assert(setPartners \|\| position >= 0);

	// set partners if required
	if (setPartners) {
	// In the case of more than one candidate colour partners,
	// there are now two approaches to choosing the partner. The
	// first method is based on two assumptions:
	// 1) the choice of which is the colour partner is done
	// randomly between the available candidates.
	// 2) the choice of which is the colour partner is done
	// independently from each particle: in other words,
	// if for a particle "i" its selected colour partner is
	// the particle "j", then the colour partner of "j"
	// does not have to be necessarily "i".
	// The second method always chooses the furthest partner
	// for hard gluons and gluinos.
	// random choice
	if( partnerMethod_ == 0 ) {
	// random choice of partner
	position = UseRandom::irnd(partners.size());
	}
	// take the one with largest angle
	else if (partnerMethod_ == 1 ) {
	if (sp->perturbative() == 1 &&
	sp->dataPtr()->iColour()==PDT::Colour8 ) {
	assert(partners.size()==2);
	// Determine largest angle
	double maxAngle(0.);
	for(unsigned int ix=0;ix<partners.size();++ix) {
	double angle = sp->momentum().vect().
	angle(partners[ix].second->momentum().vect());
	if(angle>maxAngle) {
	maxAngle = angle;
	position = ix;
	}
	}
	}
	else position = UseRandom::irnd(partners.size());
	}
	else assert(false);
	// set the evolution partner
	sp->partner(partners[position].second);
	}

	// primary partner set, set the others and do the scale
	for(size_t ix=0; ix<partners.size(); ++ix) {
	sp->addPartner(ShowerParticle::EvolutionPartner(partners[ix].second,1.,partners[ix].first,
	scales[ix].first));
	}

	// set scales for all interactions to that of the partner, default
	Energy scale = scales[position].first;
	for(unsigned int ix=0;ix<partners.size();++ix) {
	if(partners[ix].first==ShowerPartnerType::QCDColourLine) {
	sp->scales().QCD_c =
	sp->scales().QCD_c_noAO =
	(scaleChoice_==0 ? scale : scales[ix].first);
	}
	else if(partners[ix].first==ShowerPartnerType::QCDAntiColourLine) {
	sp->scales().QCD_ac =
	sp->scales().QCD_ac_noAO =
	(scaleChoice_==0 ? scale : scales[ix].first);
	}
	else assert(false);
	}
	}
	}

	void PartnerFinder::setInitialQEDEvolutionScales(const ShowerParticleVector &particles,
	const bool isDecayCase,
	const bool setPartners) {
	// loop over all the particles
	for(const auto & sp : particles) {
	// not charged or photon continue
	- if(!sp->dataPtr()->charged()) continue;
	+ if(!sp->dataPtr()->charged() && sp->dataPtr()->id()!=ParticleID::gamma) continue;
	// find the potential partners
	vector<pair<double,tShowerParticlePtr> > partners = findQEDPartners(sp,particles,isDecayCase);
	if(partners.empty()) {
	throw Exception() << "Failed to find partner in "
	<< "PartnerFinder::setQEDInitialEvolutionScales"
	<< *sp << Exception::eventerror;
	}
	// calculate the probabilities
	double prob(0.);
	for(unsigned int ix=0;ix<partners.size();++ix) prob += partners[ix].first;
	// normalise
	for(unsigned int ix=0;ix<partners.size();++ix) partners[ix].first /= prob;
	// set the partner if required
	int position(-1);
	// use QCD partner if set
	if(!setPartners&&sp->partner()) {
	for(unsigned int ix=0;ix<partners.size();++ix) {
	if(sp->partner()==partners[ix].second) {
	position = ix;
	break;
	}
	}
	}
	// set the partner
	if(setPartners\|\|!sp->partner()\|\|position<0) {
	prob = UseRandom::rnd();
	for(unsigned int ix=0;ix<partners.size();++ix) {
	if(partners[ix].first>prob) {
	position = ix;
	break;
	}
	prob -= partners[ix].first;
	}
	if(position>=0&&(setPartners\|\|!sp->partner())) {
	sp->partner(partners[position].second);
	}
	}
	// must have a partner
	if(position<0) throw Exception() << "Failed to find partner in "
	<< "PartnerFinder::setQEDInitialEvolutionScales"
	<< *sp << Exception::eventerror;
	// Calculate the evolution scales for all possible pairs of of particles
	vector<pair<Energy,Energy> > scales;
	for(unsigned int ix=0;ix< partners.size();++ix) {
	scales.push_back(calculateInitialEvolutionScales(ShowerPPair(sp,partners[ix].second),
	isDecayCase));
	}
	// store all the possible partners
	for(unsigned int ix=0;ix<partners.size();++ix) {
	sp->addPartner(ShowerParticle::EvolutionPartner(partners[ix].second,
	partners[ix].first,
	ShowerPartnerType::QED,
	scales[ix].first));
	}
	// set scales
	sp->scales().QED = scales[position].first;
	sp->scales().QED_noAO = scales[position].first;
	}
	}

	pair<Energy,Energy> PartnerFinder::
	calculateInitialEvolutionScales(const ShowerPPair &particlePair,
	const bool isDecayCase, int key) {
	bool FS1=FS(particlePair.first),FS2= FS(particlePair.second);
	if(FS1 && FS2){
	return calculateFinalFinalScales(particlePair.first->momentum(),particlePair.second->momentum(), key);
	}
	else if(FS1 && !FS2) {
	pair<Energy,Energy> rval = calculateInitialFinalScales(particlePair.second->momentum(),
	particlePair.first->momentum(),
	isDecayCase);
	return { rval.second, rval.first };
	}
	else if(!FS1 &&FS2)
	return calculateInitialFinalScales(particlePair.first->momentum(),particlePair.second->momentum(),isDecayCase);
	else
	return calculateInitialInitialScales(particlePair.first->momentum(),particlePair.second->momentum());
	}

	vector< pair<ShowerPartnerType, tShowerParticlePtr> >
	PartnerFinder::findQCDPartners(tShowerParticlePtr particle,
	const ShowerParticleVector &particles) {
	vector< pair<ShowerPartnerType, tShowerParticlePtr> > partners;
	for(const auto & sp : particles) {
	if(!sp->data().coloured() \|\| particle==sp) continue;
	// one initial-state and one final-state particle
	if(FS(particle) != FS(sp)) {
	// loop over all the colours of both particles
	for(size_t ix=0; ix<CLSIZE(particle); ++ix) {
	for(size_t jx=0; jx<CLSIZE(sp); ++jx) {
	if((CL(particle,ix) && CL(particle,ix)==CL(sp,jx))) {
	partners.push_back({ ShowerPartnerType:: QCDColourLine, sp });
	}
	}
	}
	//loop over all the anti-colours of both particles
	for(size_t ix=0; ix<ACLSIZE(particle); ++ix) {
	for(size_t jx=0; jx<ACLSIZE(sp); ++jx) {
	if((ACL(particle,ix) && ACL(particle,ix)==ACL(sp,jx))) {
	partners.push_back({ ShowerPartnerType::QCDAntiColourLine, sp });
	}
	}
	}
	}
	// two initial-state or two final-state particles
	else {
	//loop over the colours of the first particle and the anti-colours of the other
	for(size_t ix=0; ix<CLSIZE(particle); ++ix){
	for(size_t jx=0; jx<ACLSIZE(sp); ++jx){
	if(CL(particle,ix) && CL(particle,ix)==ACL(sp,jx)) {
	partners.push_back({ ShowerPartnerType:: QCDColourLine, sp });
	}
	}
	}
	//loop over the anti-colours of the first particle and the colours of the other
	for(size_t ix=0; ix<ACLSIZE(particle); ++ix){
	for(size_t jx=0; jx<CLSIZE(sp); jx++){
	if(ACL(particle,ix) && ACL(particle,ix)==CL(sp,jx)) {
	partners.push_back({ ShowerPartnerType::QCDAntiColourLine, sp });
	}
	}
	}
	}
	}
	// if we haven't found any partners look for RPV
	if (partners.empty()) {
	// special for RPV
	tColinePtr col = CL(particle);
	if(FS(particle)&&col&&col->sourceNeighbours().first) {
	tColinePair cpair = col->sourceNeighbours();
	for(const auto & sp : particles) {
	if(( FS(sp) && ( CL(sp) == cpair.first \|\| CL(sp) == cpair.second))\|\|
	(!FS(sp) && (ACL(sp) == cpair.first \|\| ACL(sp) == cpair.second ))) {
	partners.push_back({ ShowerPartnerType:: QCDColourLine, sp });
	}
	}
	}
	else if(col&&col->sinkNeighbours().first) {
	tColinePair cpair = col->sinkNeighbours();
	for(const auto & sp : particles) {
	if(( FS(sp) && (ACL(sp) == cpair.first \|\| ACL(sp) == cpair.second))\|\|
	(!FS(sp) && ( CL(sp) == cpair.first \|\| CL(sp) == cpair.second))) {
	partners.push_back({ ShowerPartnerType:: QCDColourLine, sp });
	}
	}
	}
	col = ACL(particle);
	if(FS(particle)&&col&&col->sinkNeighbours().first) {
	tColinePair cpair = col->sinkNeighbours();
	for(const auto & sp : particles) {
	if(( FS(sp) && (ACL(sp) == cpair.first \|\| ACL(sp) == cpair.second))\|\|
	(!FS(sp) && ( CL(sp) == cpair.first \|\| CL(sp) == cpair.second ))) {
	partners.push_back({ ShowerPartnerType::QCDAntiColourLine, sp });
	}
	}
	}
	else if(col&&col->sourceNeighbours().first) {
	tColinePair cpair = col->sourceNeighbours();
	for(const auto & sp : particles) {
	if(( FS(sp) && ( CL(sp) == cpair.first \|\| CL(sp) == cpair.second))\|\|
	(!FS(sp) && (ACL(sp) == cpair.first \|\|ACL(sp) == cpair.second))) {
	partners.push_back({ ShowerPartnerType::QCDAntiColourLine, sp });
	}
	}
	}
	}
	// return the partners
	return partners;
	}

	vector< pair<double, tShowerParticlePtr> >
	PartnerFinder::findQEDPartners(tShowerParticlePtr particle,
	const ShowerParticleVector & particles,
	const bool isDecayCase) {
	vector< pair<double, tShowerParticlePtr> > partners;
	const double pcharge =
	particle->id()==ParticleID::gamma ? 1 : double(particle->data().iCharge());
	vector< pair<double, tShowerParticlePtr> > photons;
	for (const auto & sp : particles) {
	if (particle == sp) continue;
	if (sp->id()==ParticleID::gamma) photons.push_back(make_pair(1.,sp));
	if (!sp->data().charged() ) continue;
	double charge = pcharge*double((sp)->data().iCharge());
	if ( FS(particle) != FS(sp) ) charge *=-1.;
	if ( QEDPartner_ != 0 && !isDecayCase ) {
	// only include II and FF as requested
	if ( QEDPartner_ == 1 && FS(particle) != FS(sp) )
	continue;
	// only include IF is requested
	else if (QEDPartner_ == 2 && FS(particle) == FS(sp) )
	continue;
	}
	if (particle->id()==ParticleID::gamma) charge = -abs(charge);
	// only keep positive dipoles
	if (charge<0.) partners.push_back({ -charge, sp });
	}
	if (particle->id()==ParticleID::gamma && partners.empty())
	return photons;
	return partners;
	}

	void PartnerFinder::setInitialEWEvolutionScales(const ShowerParticleVector &particles,
	const bool isDecayCase,
	const bool setPartners) {
	// loop over all the particles
	for(ShowerParticleVector::const_iterator cit = particles.begin();
	cit != particles.end(); ++cit) {
	// if not weakly interacting continue
	if( !weaklyInteracting( (**cit).dataPtr()))
	continue;
	// find the potential partners
	vector<pair<double,tShowerParticlePtr> > partners = findEWPartners(*cit,particles,isDecayCase);
	if(partners.empty()) {
	throw Exception() << "Failed to find partner in "
	<< "PartnerFinder::setEWInitialEvolutionScales"
	<< (**cit) << Exception::eventerror;
	}
	// calculate the probabilities
	double prob(0.);
	for(unsigned int ix=0;ix<partners.size();++ix) prob += partners[ix].first;
	// normalise
	for(unsigned int ix=0;ix<partners.size();++ix) partners[ix].first /= prob;
	// set the partner if required
	int position(-1);
	// use QCD partner if set
	if(!setPartners&&(*cit)->partner()) {
	for(unsigned int ix=0;ix<partners.size();++ix) {
	if((*cit)->partner()==partners[ix].second) {
	position = ix;
	break;
	}
	}
	}
	// set the partner
	if(setPartners\|\|!(*cit)->partner()\|\|position<0) {
	prob = UseRandom::rnd();
	for(unsigned int ix=0;ix<partners.size();++ix) {
	if(partners[ix].first>prob) {
	position = ix;
	break;
	}
	prob -= partners[ix].first;
	}
	if(position>=0&&(setPartners\|\|!(*cit)->partner())) {
	(*cit)->partner(partners[position].second);
	}
	}
	// must have a partner
	if(position<0) throw Exception() << "Failed to find partner in "
	<< "PartnerFinder::setEWInitialEvolutionScales"
	<< (**cit) << Exception::eventerror;
	// Calculate the evolution scales for all possible pairs of of particles
	vector<pair<Energy,Energy> > scales;
	for(unsigned int ix=0;ix< partners.size();++ix) {
	scales.push_back(calculateInitialEvolutionScales(ShowerPPair(*cit,partners[ix].second),
	isDecayCase, 0));
	}
	// store all the possible partners
	for(unsigned int ix=0;ix<partners.size();++ix) {
	(**cit).addPartner(ShowerParticle::EvolutionPartner(partners[ix].second,
	partners[ix].first,
	ShowerPartnerType::EW,
	scales[ix].first));
	}
	// set scales
	(**cit).scales().EW = scales[position].first;
	}
	}

	vector< pair<double, tShowerParticlePtr> >
	PartnerFinder::findEWPartners(tShowerParticlePtr particle,
	const ShowerParticleVector &particles,
	const bool isDecayCase ) {
	vector< pair<double, tShowerParticlePtr> > partners;
	for(ShowerParticlePtr partner : particles) {
	if(particle==partner \|\| !weaklyInteracting(partner->dataPtr()))
	continue;
	partners.push_back(make_pair(1.,partner));
	}
	return partners;
	}

	pair<Energy,Energy>
	PartnerFinder::calculateFinalFinalScales(
	const Lorentz5Momentum & p1,
	const Lorentz5Momentum & p2, int key)
	{
	static const double eps=1e-7;
	// Using JHEP 12(2003)045 we find that we need ktilde = 1/2(1+b-c+lambda)
	// ktilde = qtilde^2/Q^2 therefore qtilde = sqrt(ktilde*Q^2)
	// find momenta in rest frame of system
	// calculate quantities for the scales
	Energy2 Q2 = (p1+p2).m2();
	double b = p1.mass2()/Q2;
	double c = p2.mass2()/Q2;
	if(b<0.) {
	if(b<-eps) {
	throw Exception() << "Negative Mass squared b = " << b
	<< "in PartnerFinder::calculateFinalFinalScales()"
	<< Exception::eventerror;
	}
	b = 0.;
	}
	if(c<0.) {
	if(c<-eps) {
	throw Exception() << "Negative Mass squared c = " << c
	<< "in PartnerFinder::calculateFinalFinalScales()"
	<< Exception::eventerror;
	}
	c = 0.;
	}
	// KMH & PR - 16 May 2008 - swapped lambda calculation from
	// double lam=2.*p1.vect().mag()/Q; to sqrt(kallen(1,b,c)),
	// which should be identical for p1 & p2 onshell in their COM
	// but in the inverse construction for the Nason method, this
	// was not the case, leading to misuse.
	const double lam=sqrt((1.+sqrt(b)+sqrt(c))*(1.-sqrt(b)-sqrt(c))
	(sqrt(b)-1.-sqrt(c))(sqrt(c)-1.-sqrt(b)));

	Energy firstQ, secondQ;
	double kappab(0.), kappac(0.);
	//key = 0; // symmetric case pre-selection
	switch(key) {
	case 0: // symmetric case
	firstQ = sqrt(0.5Q2(1.+b-c+lam));
	secondQ = sqrt(0.5Q2(1.-b+c+lam));
	break;
	case 1: // maximum emission from both legs
	kappab=4.(1.-2.sqrt(c)-b+c);
	kappac=4.(1.-2.sqrt(b)-c+b);
	firstQ = sqrt(Q2*kappab);
	secondQ = sqrt(Q2*kappac);
	break;
	default:
	assert(false);
	}
	// calculate the scales
	return pair<Energy,Energy>(firstQ, secondQ);
	}

	pair<Energy,Energy>
	PartnerFinder::calculateInitialFinalScales(const Lorentz5Momentum& pb, const Lorentz5Momentum& pc,
	const bool isDecayCase) {
	if(!isDecayCase) {
	// In this case from JHEP 12(2003)045 we find the conditions
	// ktilde_b = (1+c) and ktilde_c = (1+2c)
	// We also find that c = m_c^2/Q^2. The process is a+b->c where
	// particle a is not colour connected (considered as a colour singlet).
	// Therefore we simply find that q_b = sqrt(Q^2+m_c^2) and
	// q_c = sqrt(Q^2+2 m_c^2)
	// We also assume that the first particle in the pair is the initial
	// state particle and the second is the final state one c
	const Energy2 mc2 = sqr(pc.mass());
	const Energy2 Q2 = -(pb-pc).m2();
	return { sqrt(Q2+mc2), sqrt(Q2+2*mc2) };
	}
	else {
	// In this case from JHEP 12(2003)045 we find, for the decay
	// process b->c+a(neutral), the condition
	// (ktilde_b-1)(ktilde_c-c)=(1/4)sqr(1-a+c+lambda).
	// We also assume that the first particle in the pair is the initial
	// state particle (b) and the second is the final state one (c).
	// - We find maximal phase space coverage through emissions from
	// c if we set ktilde_c = 4.*(sqr(1.-sqrt(a))-c)
	// - We find the most 'symmetric' way to populate the phase space
	// occurs for (ktilde_b-1)=(ktilde_c-c)=(1/2)*(1-a+c+lambda)
	// - We find the most 'smooth' way to populate the phase space
	// occurs for...
	Energy2 mb2(sqr(pb.mass()));
	double a=(pb-pc).m2()/mb2;
	double c=sqr(pc.mass())/mb2;
	double lambda = 1. + aa + cc - 2.a - 2.c - 2.ac;
	lambda = sqrt(max(lambda,0.));
	const double PROD = 0.25*sqr(1. - a + c + lambda);
	int key = 0;
	double ktilde_b, ktilde_c, cosi = 0.;
	switch (key) {
	case 0: // the 'symmetric' choice
	ktilde_c = 0.5*(1-a+c+lambda) + c ;
	ktilde_b = 1.+PROD/(ktilde_c-c) ;
	break;
	case 1: // the 'maximal' choice
	ktilde_c = 4.0*(sqr(1.-sqrt(a))-c);
	ktilde_b = 1.+PROD/(ktilde_c-c) ;
	break;
	case 2: // the 'smooth' choice
	c = max(c,1.*GeV2/mb2);
	cosi = (sqr(1-sqrt(c))-a)/lambda;
	ktilde_b = 2.0/(1.0-cosi);
	ktilde_c = (1.0-a+c+lambda)(1.0+c-a-lambdacosi)/(2.0*(1.0+cosi));
	break;
	}
	return { sqrt(mb2ktilde_b), sqrt(mb2ktilde_c) };
	}
	}

	pair<Energy,Energy>
	PartnerFinder::calculateInitialInitialScales(const Lorentz5Momentum& p1, const Lorentz5Momentum& p2) {
	// This case is quite simple. From JHEP 12(2003)045 we find the condition
	// that ktilde_b = ktilde_c = 1. In this case we have the process
	// b+c->a so we need merely boost to the CM frame of the two incoming
	// particles and then qtilde is equal to the energy in that frame
	const Energy Q = sqrt((p1+p2).m2());
	return {Q,Q};
	}
	diff --git a/Shower/QTilde/SplittingFunctions/SplittingFunction.cc b/Shower/QTilde/SplittingFunctions/SplittingFunction.cc
	--- a/Shower/QTilde/SplittingFunctions/SplittingFunction.cc
	+++ b/Shower/QTilde/SplittingFunctions/SplittingFunction.cc
	@@ -1,1048 +1,1049 @@
	// -- C++ --
	//
	// SplittingFunction.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2019 The Herwig Collaboration
	//
	// Herwig is licenced under version 3 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 SplittingFunction class.
	//

	#include "SplittingFunction.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Utilities/EnumIO.h"
	#include "Herwig/Shower/QTilde/Base/ShowerParticle.h"
	#include "ThePEG/Utilities/DescribeClass.h"

	using namespace Herwig;

	DescribeAbstractClass<SplittingFunction,Interfaced>
	describeSplittingFunction ("Herwig::SplittingFunction","");

	void SplittingFunction::Init() {

	static ClassDocumentation<SplittingFunction> documentation
	("The SplittingFunction class is the based class for 1->2 splitting functions"
	" in Herwig");

	static Switch<SplittingFunction,ColourStructure> interfaceColourStructure
	("ColourStructure",
	"The colour structure for the splitting function",
	&SplittingFunction::_colourStructure, Undefined, false, false);
	static SwitchOption interfaceColourStructureTripletTripletOctet
	(interfaceColourStructure,
	"TripletTripletOctet",
	"3 -> 3 8",
	TripletTripletOctet);
	static SwitchOption interfaceColourStructureOctetOctetOctet
	(interfaceColourStructure,
	"OctetOctetOctet",
	"8 -> 8 8",
	OctetOctetOctet);
	static SwitchOption interfaceColourStructureOctetTripletTriplet
	(interfaceColourStructure,
	"OctetTripletTriplet",
	"8 -> 3 3bar",
	OctetTripletTriplet);
	static SwitchOption interfaceColourStructureTripletOctetTriplet
	(interfaceColourStructure,
	"TripletOctetTriplet",
	"3 -> 8 3",
	- TripletOctetTriplet);
	+ TripletOctetTriplet);
	static SwitchOption interfaceColourStructureSextetSextetOctet
	(interfaceColourStructure,
	"SextetSextetOctet",
	"6 -> 6 8",
	SextetSextetOctet);
	-
	+
	static SwitchOption interfaceColourStructureChargedChargedNeutral
	(interfaceColourStructure,
	"ChargedChargedNeutral",
	"q -> q 0",
	ChargedChargedNeutral);
	static SwitchOption interfaceColourStructureNeutralChargedCharged
	(interfaceColourStructure,
	"NeutralChargedCharged",
	"0 -> q qbar",
	NeutralChargedCharged);
	static SwitchOption interfaceColourStructureChargedNeutralCharged
	(interfaceColourStructure,
	"ChargedNeutralCharged",
	"q -> 0 q",
	ChargedNeutralCharged);
	static SwitchOption interfaceColourStructureEW
	(interfaceColourStructure,
	"EW",
	"q -> q W/Z",
	EW);

	- static Switch<SplittingFunction,ShowerInteraction>
	+ static Switch<SplittingFunction,ShowerInteraction>
	interfaceInteractionType
	("InteractionType",
	"Type of the interaction",
	- &SplittingFunction::_interactionType,
	+ &SplittingFunction::_interactionType,
	ShowerInteraction::UNDEFINED, false, false);
	static SwitchOption interfaceInteractionTypeQCD
	(interfaceInteractionType,
	"QCD","QCD",ShowerInteraction::QCD);
	static SwitchOption interfaceInteractionTypeQED
	(interfaceInteractionType,
	"QED","QED",ShowerInteraction::QED);
	static SwitchOption interfaceInteractionTypeEW
	(interfaceInteractionType,
	"EW","EW",ShowerInteraction::EW);

	static Switch<SplittingFunction,bool> interfaceAngularOrdered
	("AngularOrdered",
	"Whether or not this interaction is angular ordered, "
	"normally only g->q qbar and gamma-> f fbar are the only ones which aren't.",
	&SplittingFunction::angularOrdered_, true, false, false);
	static SwitchOption interfaceAngularOrderedYes
	(interfaceAngularOrdered,
	"Yes",
	"Interaction is angular ordered",
	true);
	static SwitchOption interfaceAngularOrderedNo
	(interfaceAngularOrdered,
	"No",
	"Interaction isn't angular ordered",
	false);

	static Switch<SplittingFunction,unsigned int> interfaceScaleChoice
	("ScaleChoice",
	"The scale choice to be used",
	&SplittingFunction::scaleChoice_, 2, false, false);
	static SwitchOption interfaceScaleChoicepT
	(interfaceScaleChoice,
	"pT",
	"pT of the branching",
	0);
	static SwitchOption interfaceScaleChoiceQ2
	(interfaceScaleChoice,
	"Q2",
	"Q2 of the branching",
	1);
	static SwitchOption interfaceScaleChoiceFromAngularOrdering
	(interfaceScaleChoice,
	"FromAngularOrdering",
	"If angular order use pT, otherwise Q2",
	2);

	static Switch<SplittingFunction,bool> interfaceStrictAO
	("StrictAO",
	"Whether or not to apply strict angular-ordering,"
	" i.e. for QED even in QCD emission, and vice versa",
	&SplittingFunction::strictAO_, true, false, false);
	static SwitchOption interfaceStrictAOYes
	(interfaceStrictAO,
	"Yes",
	"Apply strict ordering",
	true);
	static SwitchOption interfaceStrictAONo
	(interfaceStrictAO,
	"No",
	"Don't apply strict ordering",
	false);
	}

	void SplittingFunction::persistentOutput(PersistentOStream & os) const {
	os << oenum(_interactionType)
	<< oenum(_colourStructure) << _colourFactor
	<< angularOrdered_ << scaleChoice_ << strictAO_;
	}

	void SplittingFunction::persistentInput(PersistentIStream & is, int) {
	is >> ienum(_interactionType)
	>> ienum(_colourStructure) >> _colourFactor
	>> angularOrdered_ >> scaleChoice_ >> strictAO_;
	}

	void SplittingFunction::colourConnection(tShowerParticlePtr parent,
	tShowerParticlePtr first,
	tShowerParticlePtr second,
	- ShowerPartnerType partnerType,
	+ ShowerPartnerType partnerType,
	const bool back) const {
	if(_colourStructure==TripletTripletOctet) {
	if(!back) {
	- ColinePair cparent = ColinePair(parent->colourLine(),
	+ ColinePair cparent = ColinePair(parent->colourLine(),
	parent->antiColourLine());
	// ensure input consistency
	- assert(( cparent.first && !cparent.second &&
	- partnerType==ShowerPartnerType::QCDColourLine) \|\|
	- ( !cparent.first && cparent.second &&
	+ assert(( cparent.first && !cparent.second &&
	+ partnerType==ShowerPartnerType::QCDColourLine) \|\|
	+ ( !cparent.first && cparent.second &&
	partnerType==ShowerPartnerType::QCDAntiColourLine));
	// q -> q g
	if(cparent.first) {
	ColinePtr newline=new_ptr(ColourLine());
	cparent.first->addColoured(second);
	newline->addColoured ( first);
	newline->addAntiColoured (second);
	}
	// qbar -> qbar g
	else {
	ColinePtr newline=new_ptr(ColourLine());
	cparent.second->addAntiColoured(second);
	newline->addColoured(second);
	newline->addAntiColoured(first);
	}
	// Set progenitor
	first->progenitor(parent->progenitor());
	second->progenitor(parent->progenitor());
	}
	else {
	- ColinePair cfirst = ColinePair(first->colourLine(),
	+ ColinePair cfirst = ColinePair(first->colourLine(),
	first->antiColourLine());
	// ensure input consistency
	- assert(( cfirst.first && !cfirst.second &&
	- partnerType==ShowerPartnerType::QCDColourLine) \|\|
	- ( !cfirst.first && cfirst.second &&
	+ assert(( cfirst.first && !cfirst.second &&
	+ partnerType==ShowerPartnerType::QCDColourLine) \|\|
	+ ( !cfirst.first && cfirst.second &&
	partnerType==ShowerPartnerType::QCDAntiColourLine));
	// q -> q g
	if(cfirst.first) {
	ColinePtr newline=new_ptr(ColourLine());
	cfirst.first->addAntiColoured(second);
	newline->addColoured(second);
	newline->addColoured(parent);
	}
	// qbar -> qbar g
	else {
	ColinePtr newline=new_ptr(ColourLine());
	cfirst.second->addColoured(second);
	newline->addAntiColoured(second);
	newline->addAntiColoured(parent);
	}
	// Set progenitor
	parent->progenitor(first->progenitor());
	second->progenitor(first->progenitor());
	}
	}
	else if(_colourStructure==OctetOctetOctet) {
	if(!back) {
	- ColinePair cparent = ColinePair(parent->colourLine(),
	+ ColinePair cparent = ColinePair(parent->colourLine(),
	parent->antiColourLine());
	// ensure input consistency
	assert(cparent.first&&cparent.second);
	// ensure first gluon is hardest
	- if( first->id()==second->id() && parent->showerKinematics()->z()<0.5 )
	+ if( first->id()==second->id() && parent->showerKinematics()->z()<0.5 )
	swap(first,second);
	// colour line radiates
	if(partnerType==ShowerPartnerType::QCDColourLine) {
	// The colour line is radiating
	ColinePtr newline=new_ptr(ColourLine());
	cparent.first->addColoured(second);
	cparent.second->addAntiColoured(first);
	newline->addColoured(first);
	newline->addAntiColoured(second);
	}
	// anti colour line radiates
	else if(partnerType==ShowerPartnerType::QCDAntiColourLine) {
	ColinePtr newline=new_ptr(ColourLine());
	cparent.first->addColoured(first);
	cparent.second->addAntiColoured(second);
	newline->addColoured(second);
	newline->addAntiColoured(first);
	}
	else
	assert(false);
	}
	else {
	- ColinePair cfirst = ColinePair(first->colourLine(),
	+ ColinePair cfirst = ColinePair(first->colourLine(),
	first->antiColourLine());
	// ensure input consistency
	assert(cfirst.first&&cfirst.second);
	// The colour line is radiating
	if(partnerType==ShowerPartnerType::QCDColourLine) {
	ColinePtr newline=new_ptr(ColourLine());
	cfirst.first->addAntiColoured(second);
	cfirst.second->addAntiColoured(parent);
	newline->addColoured(parent);
	newline->addColoured(second);
	}
	// anti colour line radiates
	else if(partnerType==ShowerPartnerType::QCDAntiColourLine) {
	ColinePtr newline=new_ptr(ColourLine());
	cfirst.first->addColoured(parent);
	cfirst.second->addColoured(second);
	newline->addAntiColoured(second);
	newline->addAntiColoured(parent);
	}
	else
	assert(false);
	- }
	+ }
	}
	else if(_colourStructure == OctetTripletTriplet) {
	if(!back) {
	- ColinePair cparent = ColinePair(parent->colourLine(),
	+ ColinePair cparent = ColinePair(parent->colourLine(),
	parent->antiColourLine());
	// ensure input consistency
	assert(cparent.first&&cparent.second);
	cparent.first ->addColoured ( first);
	cparent.second->addAntiColoured(second);
	// Set progenitor
	first->progenitor(parent->progenitor());
	second->progenitor(parent->progenitor());
	}
	else {
	- ColinePair cfirst = ColinePair(first->colourLine(),
	+ ColinePair cfirst = ColinePair(first->colourLine(),
	first->antiColourLine());
	// ensure input consistency
	assert(( cfirst.first && !cfirst.second) \|\|
	(!cfirst.first && cfirst.second));
	// g -> q qbar
	if(cfirst.first) {
	ColinePtr newline=new_ptr(ColourLine());
	cfirst.first->addColoured(parent);
	newline->addAntiColoured(second);
	- newline->addAntiColoured(parent);
	+ newline->addAntiColoured(parent);
	}
	// g -> qbar q
	else {
	ColinePtr newline=new_ptr(ColourLine());
	cfirst.second->addAntiColoured(parent);
	newline->addColoured(second);
	newline->addColoured(parent);
	}
	// Set progenitor
	parent->progenitor(first->progenitor());
	second->progenitor(first->progenitor());
	}
	}
	else if(_colourStructure == TripletOctetTriplet) {
	if(!back) {
	- ColinePair cparent = ColinePair(parent->colourLine(),
	+ ColinePair cparent = ColinePair(parent->colourLine(),
	parent->antiColourLine());
	// ensure input consistency
	- assert(( cparent.first && !cparent.second) \|\|
	+ assert(( cparent.first && !cparent.second) \|\|
	(!cparent.first && cparent.second));
	// q -> g q
	if(cparent.first) {
	ColinePtr newline=new_ptr(ColourLine());
	cparent.first->addColoured(first);
	newline->addColoured (second);
	newline->addAntiColoured( first);
	}
	// qbar -> g qbar
	else {
	ColinePtr newline=new_ptr(ColourLine());
	cparent.second->addAntiColoured(first);
	newline->addColoured ( first);
	- newline->addAntiColoured(second);
	+ newline->addAntiColoured(second);
	}
	// Set progenitor
	first->progenitor(parent->progenitor());
	second->progenitor(parent->progenitor());
	}
	else {
	- ColinePair cfirst = ColinePair(first->colourLine(),
	+ ColinePair cfirst = ColinePair(first->colourLine(),
	first->antiColourLine());
	// ensure input consistency
	assert(cfirst.first&&cfirst.second);
	// q -> g q
	if(parent->id()>0) {
	cfirst.first ->addColoured(parent);
	cfirst.second->addColoured(second);
	}
	else {
	cfirst.first ->addAntiColoured(second);
	cfirst.second->addAntiColoured(parent);
	}
	// Set progenitor
	parent->progenitor(first->progenitor());
	second->progenitor(first->progenitor());
	}
	}
	else if(_colourStructure==SextetSextetOctet) {
	//make sure we're not doing backward evolution
	assert(!back);

	//make sure something sensible
	assert(parent->colourLine() \|\| parent->antiColourLine());
	-
	+
	//get the colour lines or anti-colour lines
	bool isAntiColour=true;
	ColinePair cparent;
	if(parent->colourLine()) {
	- cparent = ColinePair(const_ptr_cast<tColinePtr>(parent->colourInfo()->colourLines()[0]),
	+ cparent = ColinePair(const_ptr_cast<tColinePtr>(parent->colourInfo()->colourLines()[0]),
	const_ptr_cast<tColinePtr>(parent->colourInfo()->colourLines()[1]));
	isAntiColour=false;
	}
	else {
	- cparent = ColinePair(const_ptr_cast<tColinePtr>(parent->colourInfo()->antiColourLines()[0]),
	+ cparent = ColinePair(const_ptr_cast<tColinePtr>(parent->colourInfo()->antiColourLines()[0]),
	const_ptr_cast<tColinePtr>(parent->colourInfo()->antiColourLines()[1]));
	}
	-
	+
	//check for sensible input
	// assert(cparent.first && cparent.second);

	// sextet has 2 colour lines
	if(!isAntiColour) {
	//pick at random which of the colour topolgies to take
	double topology = UseRandom::rnd();
	if(topology < 0.25) {
	ColinePtr newline=new_ptr(ColourLine());
	cparent.first->addColoured(second);
	cparent.second->addColoured(first);
	newline->addColoured(first);
	newline->addAntiColoured(second);
	}
	else if(topology >=0.25 && topology < 0.5) {
	ColinePtr newline=new_ptr(ColourLine());
	cparent.first->addColoured(first);
	cparent.second->addColoured(second);
	newline->addColoured(first);
	- newline->addAntiColoured(second);
	+ newline->addAntiColoured(second);
	}
	else if(topology >= 0.5 && topology < 0.75) {
	ColinePtr newline=new_ptr(ColourLine());
	cparent.first->addColoured(second);
	- cparent.second->addColoured(first);
	- newline->addColoured(first);
	- newline->addAntiColoured(second);
	+ cparent.second->addColoured(first);
	+ newline->addColoured(first);
	+ newline->addAntiColoured(second);
	}
	else {
	ColinePtr newline=new_ptr(ColourLine());
	cparent.first->addColoured(first);
	cparent.second->addColoured(second);
	newline->addColoured(first);
	newline->addAntiColoured(second);
	}
	}
	// sextet has 2 anti-colour lines
	else {
	double topology = UseRandom::rnd();
	if(topology < 0.25){
	ColinePtr newline=new_ptr(ColourLine());
	cparent.first->addAntiColoured(second);
	cparent.second->addAntiColoured(first);
	newline->addAntiColoured(first);
	newline->addColoured(second);
	}
	else if(topology >=0.25 && topology < 0.5) {
	ColinePtr newline=new_ptr(ColourLine());
	cparent.first->addAntiColoured(first);
	cparent.second->addAntiColoured(second);
	newline->addAntiColoured(first);
	- newline->addColoured(second);
	+ newline->addColoured(second);
	}
	else if(topology >= 0.5 && topology < 0.75) {
	ColinePtr newline=new_ptr(ColourLine());
	cparent.first->addAntiColoured(second);
	cparent.second->addAntiColoured(first);
	newline->addAntiColoured(first);
	- newline->addColoured(second);
	+ newline->addColoured(second);
	}
	else {
	ColinePtr newline=new_ptr(ColourLine());
	cparent.first->addAntiColoured(first);
	cparent.second->addAntiColoured(second);
	newline->addAntiColoured(first);
	newline->addColoured(second);
	}
	- }
	+ }
	}
	else if(_colourStructure == ChargedChargedNeutral) {
	if(!parent->data().coloured()) return;
	if(!back) {
	- ColinePair cparent = ColinePair(parent->colourLine(),
	+ ColinePair cparent = ColinePair(parent->colourLine(),
	parent->antiColourLine());
	// q -> q g
	if(cparent.first) {
	cparent.first->addColoured(first);
	}
	// qbar -> qbar g
	if(cparent.second) {
	cparent.second->addAntiColoured(first);
	}
	}
	else {
	- ColinePair cfirst = ColinePair(first->colourLine(),
	+ ColinePair cfirst = ColinePair(first->colourLine(),
	first->antiColourLine());
	// q -> q g
	if(cfirst.first) {
	cfirst.first->addColoured(parent);
	}
	// qbar -> qbar g
	if(cfirst.second) {
	cfirst.second->addAntiColoured(parent);
	}
	}
	}
	else if(_colourStructure == ChargedNeutralCharged) {
	if(!parent->data().coloured()) return;
	if(!back) {
	- ColinePair cparent = ColinePair(parent->colourLine(),
	+ ColinePair cparent = ColinePair(parent->colourLine(),
	parent->antiColourLine());
	// q -> q g
	if(cparent.first) {
	cparent.first->addColoured(second);
	}
	// qbar -> qbar g
	if(cparent.second) {
	cparent.second->addAntiColoured(second);
	}
	}
	else {
	if (second->dataPtr()->iColour()==PDT::Colour3 ) {
	ColinePtr newline=new_ptr(ColourLine());
	newline->addColoured(second);
	newline->addColoured(parent);
	}
	else if (second->dataPtr()->iColour()==PDT::Colour3bar ) {
	ColinePtr newline=new_ptr(ColourLine());
	newline->addAntiColoured(second);
	newline->addAntiColoured(parent);
	}
	}
	}
	else if(_colourStructure == NeutralChargedCharged ) {
	if(!back) {
	if(first->dataPtr()->coloured()) {
	ColinePtr newline=new_ptr(ColourLine());
	if(first->dataPtr()->iColour()==PDT::Colour3) {
	newline->addColoured (first );
	newline->addAntiColoured(second);
	}
	else if (first->dataPtr()->iColour()==PDT::Colour3bar) {
	newline->addColoured (second);
	newline->addAntiColoured(first );
	}
	else
	assert(false);
	}
	}
	- else {
	- ColinePair cfirst = ColinePair(first->colourLine(),
	+ else {
	+ ColinePair cfirst = ColinePair(first->colourLine(),
	first->antiColourLine());
	// gamma -> q qbar
	if(cfirst.first) {
	cfirst.first->addAntiColoured(second);
	}
	// gamma -> qbar q
	else if(cfirst.second) {
	cfirst.second->addColoured(second);
	}
	- else
	+ else
	assert(false);
	}
	}
	else if(_colourStructure == EW) {
	if(!parent->data().coloured()) return;
	if(!back) {
	- ColinePair cparent = ColinePair(parent->colourLine(),
	+ ColinePair cparent = ColinePair(parent->colourLine(),
	parent->antiColourLine());
	// q -> q g
	if(cparent.first) {
	cparent.first->addColoured(first);
	}
	// qbar -> qbar g
	if(cparent.second) {
	cparent.second->addAntiColoured(first);
	}
	}
	else {
	- ColinePair cfirst = ColinePair(first->colourLine(),
	+ ColinePair cfirst = ColinePair(first->colourLine(),
	first->antiColourLine());
	// q -> q g
	if(cfirst.first) {
	cfirst.first->addColoured(parent);
	}
	// qbar -> qbar g
	if(cfirst.second) {
	cfirst.second->addAntiColoured(parent);
	}
	}
	}
	else {
	assert(false);
	}
	}

	void SplittingFunction::doinit() {
	Interfaced::doinit();
	assert(_interactionType!=ShowerInteraction::UNDEFINED);
	assert((_colourStructure>0&&_interactionType==ShowerInteraction::QCD) \|\|
	(_colourStructure<0&&(_interactionType==ShowerInteraction::QED \|\|
	_interactionType==ShowerInteraction::EW)) );
	if(_colourFactor>0.) return;
	// compute the colour factors if need
	if(_colourStructure==TripletTripletOctet) {
	_colourFactor = 4./3.;
	}
	else if(_colourStructure==OctetOctetOctet) {
	_colourFactor = 3.;
	}
	else if(_colourStructure==OctetTripletTriplet) {
	_colourFactor = 0.5;
	}
	else if(_colourStructure==TripletOctetTriplet) {
	_colourFactor = 4./3.;
	}
	else if(_colourStructure==SextetSextetOctet) {
	_colourFactor = 10./3.;
	}
	else if(_colourStructure<0) {
	_colourFactor = 1.;
	}
	else {
	assert(false);
	}
	}

	bool SplittingFunction::checkColours(const IdList & ids) const {
	if(_colourStructure==TripletTripletOctet) {
	if(ids[0]!=ids[1]) return false;
	if((ids[0]->iColour()==PDT::Colour3\|\|ids[0]->iColour()==PDT::Colour3bar) &&
	ids[2]->iColour()==PDT::Colour8) return true;
	return false;
	}
	else if(_colourStructure==OctetOctetOctet) {
	for(unsigned int ix=0;ix<3;++ix) {
	if(ids[ix]->iColour()!=PDT::Colour8) return false;
	}
	return true;
	}
	else if(_colourStructure==OctetTripletTriplet) {
	if(ids[0]->iColour()!=PDT::Colour8) return false;
	if(ids[1]->iColour()==PDT::Colour3&&ids[2]->iColour()==PDT::Colour3bar)
	return true;
	if(ids[1]->iColour()==PDT::Colour3bar&&ids[2]->iColour()==PDT::Colour3)
	return true;
	return false;
	}
	else if(_colourStructure==TripletOctetTriplet) {
	if(ids[0]!=ids[2]) return false;
	if((ids[0]->iColour()==PDT::Colour3\|\|ids[0]->iColour()==PDT::Colour3bar) &&
	ids[1]->iColour()==PDT::Colour8) return true;
	return false;
	}
	else if(_colourStructure==SextetSextetOctet) {
	if(ids[0]!=ids[1]) return false;
	if((ids[0]->iColour()==PDT::Colour6 \|\| ids[0]->iColour()==PDT::Colour6bar) &&
	ids[2]->iColour()==PDT::Colour8) return true;
	return false;
	}
	else if(_colourStructure==ChargedChargedNeutral) {
	if(ids[0]!=ids[1]) return false;
	if(ids[2]->iCharge()!=0) return false;
	if(ids[0]->iCharge()==ids[1]->iCharge()) return true;
	return false;
	}
	else if(_colourStructure==ChargedNeutralCharged) {
	if(ids[0]!=ids[2]) return false;
	if(ids[1]->iCharge()!=0) return false;
	if(ids[0]->iCharge()==ids[2]->iCharge()) return true;
	return false;
	}
	else if(_colourStructure==NeutralChargedCharged) {
	if(ids[1]->id()!=-ids[2]->id()) return false;
	if(ids[0]->iCharge()!=0) return false;
	if(ids[1]->iCharge()==-ids[2]->iCharge()) return true;
	return false;
	}
	else {
	assert(false);
	}
	return false;
	}

	namespace {

	bool hasColour(tPPtr p) {
	PDT::Colour colour = p->dataPtr()->iColour();
	return colour==PDT::Colour3 \|\| colour==PDT::Colour8 \|\| colour == PDT::Colour6;
	}

	bool hasAntiColour(tPPtr p) {
	PDT::Colour colour = p->dataPtr()->iColour();
	return colour==PDT::Colour3bar \|\| colour==PDT::Colour8 \|\| colour == PDT::Colour6bar;
	}
	-
	+
	}

	void SplittingFunction::evaluateFinalStateScales(ShowerPartnerType partnerType,
	Energy scale, double z,
	tShowerParticlePtr parent,
	tShowerParticlePtr emitter,
	tShowerParticlePtr emitted) {
	// identify emitter and emitted
	double zEmitter = z, zEmitted = 1.-z;
	bool bosonSplitting(false);
	// special for g -> gg, particle highest z is emitter
	if(emitter->id() == emitted->id() && emitter->id() == parent->id() &&
	zEmitted > zEmitter) {
	swap(zEmitted,zEmitter);
	swap( emitted, emitter);
	}
	// otherwise if particle ID same
	else if(emitted->id()==parent->id()) {
	swap(zEmitted,zEmitter);
	swap( emitted, emitter);
	}
	// no real emitter/emitted
	else if(emitter->id()!=parent->id()) {
	bosonSplitting = true;
	}
	// may need to add angularOrder flag here
	// now the various scales
	// QED
	if(partnerType==ShowerPartnerType::QED) {
	assert(colourStructure()==ChargedChargedNeutral \|\|
	colourStructure()==ChargedNeutralCharged \|\|
	- colourStructure()==NeutralChargedCharged );
	+ colourStructure()==NeutralChargedCharged \|\|
	+ colourStructure()==EW);
	// normal case
	if(!bosonSplitting) {
	- assert(colourStructure()==ChargedChargedNeutral \|\|
	- colourStructure()==ChargedNeutralCharged );
	+ assert(colourStructure()==ChargedChargedNeutral);
	// set the scales
	// emitter
	emitter->scales().QED = zEmitter*scale;
	emitter->scales().QED_noAO = scale;
	if(strictAO_)
	emitter->scales().QCD_c = min(zEmitter*scale,parent->scales().QCD_c );
	else
	emitter->scales().QCD_c = min( scale,parent->scales().QCD_c );
	emitter->scales().QCD_c_noAO = min(scale,parent->scales().QCD_c_noAO );
	if(strictAO_)
	emitter->scales().QCD_ac = min(zEmitter*scale,parent->scales().QCD_ac );
	else
	emitter->scales().QCD_ac = min( scale,parent->scales().QCD_ac );
	emitter->scales().QCD_ac_noAO = min(scale,parent->scales().QCD_ac_noAO);
	emitter->scales().EW = min(scale,parent->scales().EW );
	- // emitted
	+ // emitted
	emitted->scales().QED = zEmitted*scale;
	emitted->scales().QED_noAO = scale;
	emitted->scales().QCD_c = ZERO;
	emitted->scales().QCD_c_noAO = ZERO;
	emitted->scales().QCD_ac = ZERO;
	emitted->scales().QCD_ac_noAO = ZERO;
	emitted->scales().EW = min(scale,parent->scales().EW );
	}
	// gamma -> f fbar
	else {
	- assert(colourStructure()==NeutralChargedCharged );
	+ assert(colourStructure()==NeutralChargedCharged \|\|
	+ colourStructure()==EW);
	// emitter
	emitter->scales().QED = zEmitter*scale;
	emitter->scales().QED_noAO = scale;
	if(hasColour(emitter)) {
	emitter->scales().QCD_c = zEmitter*scale;
	emitter->scales().QCD_c_noAO = scale;
	}
	if(hasAntiColour(emitter)) {
	emitter->scales().QCD_ac = zEmitter*scale;
	emitter->scales().QCD_ac_noAO = scale;
	}
	emitter->scales().EW = zEmitter*scale;
	- // emitted
	+ // emitted
	emitted->scales().QED = zEmitted*scale;
	emitted->scales().QED_noAO = scale;
	if(hasColour(emitted)) {
	emitted->scales().QCD_c = zEmitted*scale;
	emitted->scales().QCD_c_noAO = scale;
	}
	if(hasAntiColour(emitted)) {
	emitted->scales().QCD_ac = zEmitted*scale;
	emitted->scales().QCD_ac_noAO = scale;
	}
	emitted->scales().EW = zEmitted*scale;
	}
	}
	// QCD
	else if (partnerType==ShowerPartnerType::QCDColourLine \|\|
	partnerType==ShowerPartnerType::QCDAntiColourLine) {
	// normal case eg q -> q g and g -> g g
	if(!bosonSplitting) {
	if(strictAO_)
	emitter->scales().QED = min(zEmitter*scale,parent->scales().QED );
	else
	emitter->scales().QED = min( scale,parent->scales().QED );
	emitter->scales().QED_noAO = min(scale,parent->scales().QED_noAO);
	emitter->scales().EW = min(scale,parent->scales().EW );
	if(partnerType==ShowerPartnerType::QCDColourLine) {
	emitter->scales().QCD_c = zEmitter*scale;
	emitter->scales().QCD_c_noAO = scale;
	emitter->scales().QCD_ac = min(zEmitter*scale,parent->scales().QCD_ac );
	emitter->scales().QCD_ac_noAO = min( scale,parent->scales().QCD_ac_noAO);
	}
	else {
	emitter->scales().QCD_c = min(zEmitter*scale,parent->scales().QCD_c );
	emitter->scales().QCD_c_noAO = min( scale,parent->scales().QCD_c_noAO );
	emitter->scales().QCD_ac = zEmitter*scale;
	emitter->scales().QCD_ac_noAO = scale;
	}
	- // emitted
	+ // emitted
	emitted->scales().QED = ZERO;
	emitted->scales().QED_noAO = ZERO;
	emitted->scales().QCD_c = zEmitted*scale;
	emitted->scales().QCD_c_noAO = scale;
	emitted->scales().QCD_ac = zEmitted*scale;
	emitted->scales().QCD_ac_noAO = scale;
	emitted->scales().EW = min(scale,parent->scales().EW );
	}
	// g -> q qbar
	else {
	// emitter
	if(emitter->dataPtr()->charged()) {
	emitter->scales().QED = zEmitter*scale;
	emitter->scales().QED_noAO = scale;
	}
	emitter->scales().EW = zEmitter*scale;
	emitter->scales().QCD_c = zEmitter*scale;
	emitter->scales().QCD_c_noAO = scale;
	emitter->scales().QCD_ac = zEmitter*scale;
	emitter->scales().QCD_ac_noAO = scale;
	- // emitted
	+ // emitted
	if(emitted->dataPtr()->charged()) {
	emitted->scales().QED = zEmitted*scale;
	emitted->scales().QED_noAO = scale;
	}
	emitted->scales().EW = zEmitted*scale;
	emitted->scales().QCD_c = zEmitted*scale;
	emitted->scales().QCD_c_noAO = scale;
	emitted->scales().QCD_ac = zEmitted*scale;
	emitted->scales().QCD_ac_noAO = scale;
	}
	}
	else if(partnerType==ShowerPartnerType::EW) {
	// EW
	emitter->scales().EW = zEmitter*scale;
	emitted->scales().EW = zEmitted*scale;
	// QED
	// W radiation AO
	if(emitted->dataPtr()->charged()) {
	emitter->scales().QED = zEmitter*scale;
	emitter->scales().QED_noAO = scale;
	emitted->scales().QED = zEmitted*scale;
	emitted->scales().QED_noAO = scale;
	}
	// Z don't
	else {
	emitter->scales().QED = min(scale,parent->scales().QED );
	emitter->scales().QED_noAO = min(scale,parent->scales().QED_noAO);
	emitted->scales().QED = ZERO;
	emitted->scales().QED_noAO = ZERO;
	}
	// QCD
	emitter->scales().QCD_c = min(scale,parent->scales().QCD_c );
	emitter->scales().QCD_c_noAO = min(scale,parent->scales().QCD_c_noAO );
	emitter->scales().QCD_ac = min(scale,parent->scales().QCD_ac );
	emitter->scales().QCD_ac_noAO = min(scale,parent->scales().QCD_ac_noAO);
	emitted->scales().QCD_c = ZERO;
	emitted->scales().QCD_c_noAO = ZERO;
	emitted->scales().QCD_ac = ZERO;
	emitted->scales().QCD_ac_noAO = ZERO;
	}
	else
	assert(false);
	}

	void SplittingFunction::evaluateInitialStateScales(ShowerPartnerType partnerType,
	Energy scale, double z,
	tShowerParticlePtr parent,
	tShowerParticlePtr spacelike,
	tShowerParticlePtr timelike) {
	// scale for time-like child
	Energy AOScale = (1.-z)*scale;
	// QED
	if(partnerType==ShowerPartnerType::QED) {
	if(parent->id()==spacelike->id()) {
	// parent
	parent ->scales().QED = scale;
	parent ->scales().QED_noAO = scale;
	parent ->scales().QCD_c = min(scale,spacelike->scales().QCD_c );
	parent ->scales().QCD_c_noAO = min(scale,spacelike->scales().QCD_c_noAO );
	parent ->scales().QCD_ac = min(scale,spacelike->scales().QCD_ac );
	parent ->scales().QCD_ac_noAO = min(scale,spacelike->scales().QCD_ac_noAO);
	// timelike
	timelike->scales().QED = AOScale;
	timelike->scales().QED_noAO = scale;
	timelike->scales().QCD_c = ZERO;
	timelike->scales().QCD_c_noAO = ZERO;
	timelike->scales().QCD_ac = ZERO;
	timelike->scales().QCD_ac_noAO = ZERO;
	}
	else if(parent->id()==timelike->id()) {
	parent ->scales().QED = scale;
	parent ->scales().QED_noAO = scale;
	if(hasColour(parent)) {
	parent ->scales().QCD_c = scale;
	parent ->scales().QCD_c_noAO = scale;
	}
	if(hasAntiColour(parent)) {
	parent ->scales().QCD_ac = scale;
	parent ->scales().QCD_ac_noAO = scale;
	}
	- // timelike
	+ // timelike
	timelike->scales().QED = AOScale;
	timelike->scales().QED_noAO = scale;
	if(hasColour(timelike)) {
	timelike->scales().QCD_c = AOScale;
	timelike->scales().QCD_c_noAO = scale;
	}
	if(hasAntiColour(timelike)) {
	timelike->scales().QCD_ac = AOScale;
	timelike->scales().QCD_ac_noAO = scale;
	}
	}
	else {
	parent ->scales().QED = scale;
	parent ->scales().QED_noAO = scale;
	parent ->scales().QCD_c = ZERO ;
	parent ->scales().QCD_c_noAO = ZERO ;
	parent ->scales().QCD_ac = ZERO ;
	parent ->scales().QCD_ac_noAO = ZERO ;
	- // timelike
	+ // timelike
	timelike->scales().QED = AOScale;
	timelike->scales().QED_noAO = scale;
	if(hasColour(timelike)) {
	timelike->scales().QCD_c = min(AOScale,spacelike->scales().QCD_ac );
	timelike->scales().QCD_c_noAO = min( scale,spacelike->scales().QCD_ac_noAO);
	}
	if(hasAntiColour(timelike)) {
	timelike->scales().QCD_ac = min(AOScale,spacelike->scales().QCD_c );
	timelike->scales().QCD_ac_noAO = min( scale,spacelike->scales().QCD_c_noAO );
	}
	}
	}
	// QCD
	else if (partnerType==ShowerPartnerType::QCDColourLine \|\|
	partnerType==ShowerPartnerType::QCDAntiColourLine) {
	- // timelike
	+ // timelike
	if(timelike->dataPtr()->charged()) {
	timelike->scales().QED = AOScale;
	timelike->scales().QED_noAO = scale;
	}
	if(hasColour(timelike)) {
	timelike->scales().QCD_c = AOScale;
	timelike->scales().QCD_c_noAO = scale;
	}
	if(hasAntiColour(timelike)) {
	timelike->scales().QCD_ac = AOScale;
	timelike->scales().QCD_ac_noAO = scale;
	}
	if(parent->id()==spacelike->id()) {
	parent ->scales().QED = min(scale,spacelike->scales().QED );
	parent ->scales().QED_noAO = min(scale,spacelike->scales().QED_noAO );
	parent ->scales().QCD_c = min(scale,spacelike->scales().QCD_c );
	parent ->scales().QCD_c_noAO = min(scale,spacelike->scales().QCD_c_noAO );
	parent ->scales().QCD_ac = min(scale,spacelike->scales().QCD_ac );
	parent ->scales().QCD_ac_noAO = min(scale,spacelike->scales().QCD_ac_noAO);
	}
	else {
	if(parent->dataPtr()->charged()) {
	parent ->scales().QED = scale;
	parent ->scales().QED_noAO = scale;
	}
	if(hasColour(parent)) {
	parent ->scales().QCD_c = scale;
	parent ->scales().QCD_c_noAO = scale;
	}
	if(hasAntiColour(parent)) {
	parent ->scales().QCD_ac = scale;
	parent ->scales().QCD_ac_noAO = scale;
	}
	}
	}
	else if(partnerType==ShowerPartnerType::EW) {
	- if(abs(spacelike->id())!=ParticleID::Wplus &&
	+ if(abs(spacelike->id())!=ParticleID::Wplus &&
	spacelike->id() !=ParticleID::Z0 ) {
	// QCD scales
	parent ->scales().QCD_c = min(scale,spacelike->scales().QCD_c );
	parent ->scales().QCD_c_noAO = min(scale,spacelike->scales().QCD_c_noAO );
	parent ->scales().QCD_ac = min(scale,spacelike->scales().QCD_ac );
	parent ->scales().QCD_ac_noAO = min(scale,spacelike->scales().QCD_ac_noAO);
	timelike->scales().QCD_c = ZERO;
	timelike->scales().QCD_c_noAO = ZERO;
	timelike->scales().QCD_ac = ZERO;
	timelike->scales().QCD_ac_noAO = ZERO;
	// QED scales
	if(timelike->id()==ParticleID::Z0) {
	parent ->scales().QED = min(scale,spacelike->scales().QED );
	parent ->scales().QED_noAO = min(scale,spacelike->scales().QED_noAO );
	timelike->scales().QED = ZERO;
	timelike->scales().QED_noAO = ZERO;
	}
	else {
	parent ->scales().QED = scale;
	parent ->scales().QED_noAO = scale;
	timelike->scales().QED = AOScale;
	timelike->scales().QED_noAO = scale;
	}
	// EW scales
	parent ->scales().EW = scale;
	timelike->scales().EW = AOScale;
	}
	else assert(false);
	}
	else
	assert(false);
	}

	void SplittingFunction::evaluateDecayScales(ShowerPartnerType partnerType,
	Energy scale, double z,
	tShowerParticlePtr parent,
	tShowerParticlePtr spacelike,
	tShowerParticlePtr timelike) {
	assert(parent->id()==spacelike->id());
	// angular-ordered scale for 2nd child
	Energy AOScale = (1.-z)*scale;
	// QED
	if(partnerType==ShowerPartnerType::QED) {
	// timelike
	timelike->scales().QED = AOScale;
	timelike->scales().QED_noAO = scale;
	timelike->scales().QCD_c = ZERO;
	timelike->scales().QCD_c_noAO = ZERO;
	timelike->scales().QCD_ac = ZERO;
	timelike->scales().QCD_ac_noAO = ZERO;
	timelike->scales().EW = ZERO;
	// spacelike
	spacelike->scales().QED = scale;
	spacelike->scales().QED_noAO = scale;
	spacelike->scales().EW = max(scale,parent->scales().EW );
	}
	// QCD
	else if(partnerType==ShowerPartnerType::QCDColourLine \|\|
	partnerType==ShowerPartnerType::QCDAntiColourLine) {
	- // timelike
	+ // timelike
	timelike->scales().QED = ZERO;
	timelike->scales().QED_noAO = ZERO;
	timelike->scales().QCD_c = AOScale;
	timelike->scales().QCD_c_noAO = scale;
	timelike->scales().QCD_ac = AOScale;
	timelike->scales().QCD_ac_noAO = scale;
	timelike->scales().EW = ZERO;
	// spacelike
	spacelike->scales().QED = max(scale,parent->scales().QED );
	spacelike->scales().QED_noAO = max(scale,parent->scales().QED_noAO );
	spacelike->scales().EW = max(scale,parent->scales().EW );
	}
	else if(partnerType==ShowerPartnerType::EW) {
	// EW
	timelike->scales().EW = AOScale;
	spacelike->scales().EW = max(scale,parent->scales().EW );
	// QCD
	timelike->scales().QCD_c = ZERO;
	timelike->scales().QCD_c_noAO = ZERO;
	timelike->scales().QCD_ac = ZERO;
	timelike->scales().QCD_ac_noAO = ZERO;
	timelike->scales().EW = ZERO;
	// QED
	timelike->scales().QED = ZERO;
	timelike->scales().QED_noAO = ZERO;
	spacelike->scales().QED = max(scale,parent->scales().QED );
	spacelike->scales().QED_noAO = max(scale,parent->scales().QED_noAO );
	}
	- else
	+ else
	assert(false);
	spacelike->scales().QCD_c = max(scale,parent->scales().QCD_c );
	spacelike->scales().QCD_c_noAO = max(scale,parent->scales().QCD_c_noAO );
	spacelike->scales().QCD_ac = max(scale,parent->scales().QCD_ac );
	spacelike->scales().QCD_ac_noAO = max(scale,parent->scales().QCD_ac_noAO);
	}

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 9:10 PM (1 d, 41 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3806210
Default Alt Text: (67 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions