No OneTemporary
Actions

Size

55 KB

Subscribers

None

View Options

	diff --git a/Hadronization/ClusterDecayer.cc b/Hadronization/ClusterDecayer.cc
	--- a/Hadronization/ClusterDecayer.cc
	+++ b/Hadronization/ClusterDecayer.cc
	@@ -1,807 +1,813 @@
	// -- C++ --
	//
	// ClusterDecayer.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 ClusterDecayer class.
	//

	#include "ClusterDecayer.h"
	#include <ThePEG/Interface/ClassDocumentation.h>
	#include <ThePEG/Interface/Reference.h>
	#include <ThePEG/Interface/Parameter.h>
	#include <ThePEG/Interface/Switch.h>
	#include <ThePEG/Persistency/PersistentOStream.h>
	#include <ThePEG/Persistency/PersistentIStream.h>
	#include <ThePEG/PDT/EnumParticles.h>
	#include <ThePEG/Repository/EventGenerator.h>
	#include "Herwig/Utilities/Kinematics.h"
	#include "Cluster.h"
	#include <ThePEG/Utilities/DescribeClass.h>
	#include <ThePEG/Repository/UseRandom.h>
	#include "ThePEG/Interface/ParMap.h"

	#include "Herwig/Utilities/Histogram.h"

	using namespace Herwig;

	DescribeClass<ClusterDecayer,Interfaced>
	describeClusterDecayer("Herwig::ClusterDecayer","Herwig.so");

	ClusterDecayer::ClusterDecayer() :
	_clDirLight(1),
	_clDirExotic(1),
	_clSmrLight(0.0),
	_clSmrExotic(0.0),
	_onshell(false),
	_masstry(20),
	_covariantBoost(false),
	_kinematics(0)
	{}

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

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

	void ClusterDecayer::persistentOutput(PersistentOStream & os) const
	{
	os << _clDirLight << _clDirHeavy
	<< _clDirExotic << _clSmrLight << _clSmrHeavy
	<< _clSmrExotic << _onshell << _masstry
	<< _covariantBoost
	<< _kinematics
	<< _hadronSpectrum;
	}

	void ClusterDecayer::persistentInput(PersistentIStream & is, int) {
	is >> _clDirLight >> _clDirHeavy
	>> _clDirExotic >> _clSmrLight >> _clSmrHeavy
	>> _clSmrExotic >> _onshell >> _masstry
	>> _covariantBoost
	>> _kinematics
	>> _hadronSpectrum;
	}

	void ClusterDecayer::doinit() {
	Interfaced::doinit();
	for ( const long& id : spectrum()->heavyHadronizingQuarks() ) {
	if ( _clSmrHeavy.find(id) == _clSmrHeavy.end() \|\|
	_clDirHeavy.find(id) == _clDirHeavy.end() )
	throw InitException() << "not all parameters have been set for heavy quark cluster decay";
	}
	}

	void ClusterDecayer::Init() {

	static ClassDocumentation<ClusterDecayer> documentation
	("This class is responsible for the two-body decays of normal clusters");

	static Reference<ClusterDecayer,HadronSpectrum> interfaceHadronSpectrum
	("HadronSpectrum",
	"Set the hadron spectrum for this parton splitter.",
	&ClusterDecayer::_hadronSpectrum, false, false, true, false);

	//ClDir for light, Bottom, Charm and exotic (e.g Susy) quarks

	static Switch<ClusterDecayer,bool> interfaceClDirLight
	("ClDirLight",
	"Cluster direction for light quarks",
	&ClusterDecayer::_clDirLight, true, false, false);
	static SwitchOption interfaceClDirLightPreserve
	(interfaceClDirLight,
	"Preserve",
	"Preserve the direction of the quark from the perturbative stage"
	" as the direction of the meson produced from it",
	true);
	static SwitchOption interfaceClDirLightIsotropic
	(interfaceClDirLight,
	"Isotropic",
	"Assign the direction of the meson randomly",
	false);

	static ParMap<ClusterDecayer,bool> interfaceClDirHeavy
	("ClDirHeavy",
	"Preserve the direction of the given heavy quark from the perturbative stage.",
	&ClusterDecayer::_clDirHeavy, -1, true, false, true,
	false, false, Interface::limited);

	static Switch<ClusterDecayer,bool> interfaceClDirExotic
	("ClDirExotic",
	"Cluster direction for exotic quarks",
	&ClusterDecayer::_clDirExotic, true, false, false);
	static SwitchOption interfaceClDirExoticPreserve
	(interfaceClDirExotic,
	"Preserve",
	"Preserve the direction of the quark from the perturbative stage"
	" as the direction of the meson produced from it",
	true);
	static SwitchOption interfaceClDirExoticIsotropic
	(interfaceClDirExotic,
	"Isotropic",
	"Assign the direction of the meson randomly",
	false);


	static Switch<ClusterDecayer,bool> interfaceCovariantBoost
	("CovariantBoost",
	"Use single Covariant Boost for Cluster Fission",
	&ClusterDecayer::_covariantBoost, false, false, false);
	static SwitchOption interfaceCovariantBoostYes
	(interfaceCovariantBoost,
	"Yes",
	"Use Covariant boost",
	true);
	static SwitchOption interfaceCovariantBoostNo
	(interfaceCovariantBoost,
	"No",
	"Do NOT use Covariant boost",
	false);
	static Switch<ClusterDecayer,int> interfaceKinematics
	("Kinematics",
	"Choose kinematics of Cluster Decay",
	&ClusterDecayer::_kinematics, 0, true, false);
	static SwitchOption interfaceKinematicsDefault
	(interfaceKinematics,
	"Default",
	"default kinematics",
	0);
	static SwitchOption interfaceKinematicsTchannel
	(interfaceKinematics,
	"Tchannel",
	"t channel like kinematics using the matrix element 1/(p1.p3)^2",
	1);
	static SwitchOption interfaceKinematicsAllIsotropic
	(interfaceKinematics,
	"AllIsotropic",
	"All clusters decay isotropically",
	2);
	static SwitchOption interfaceKinematicsAllAligned
	(interfaceKinematics,
	"AllAligned",
	"All clusters decay in the direction of their constituents",
	3);
	static SwitchOption interfaceKinematicsTchannelHemisphere
	(interfaceKinematics,
	"TchannelHemisphere",
	"t channel like kinematics using the matrix element 1/(p1.p3)^2 "
	"But limited the sampling to the hemisphere of the original constituents",
	-1);
	static SwitchOption interfaceKinematicsAllIsotropicHemisphere
	(interfaceKinematics,
	"IsotropicHemisphere",
	"All clusters decay semi-isotropically "
	"limited the sampling to the hemisphere of the original constituents",
	-2);

	// ClSmr for ligth, Bottom, Charm and exotic (e.g. Susy) quarks
	static Parameter<ClusterDecayer,double>
	interfaceClSmrLight ("ClSmrLight", "cluster direction Gaussian smearing for light quark",
	&ClusterDecayer::_clSmrLight, 0, 0.0 , 0.0 , 2.0,false,false,false);
	static ParMap<ClusterDecayer,double> interfaceClSmrHeavy
	("ClSmrHeavy",
	"cluster direction Gaussian smearing for heavy quarks",
	&ClusterDecayer::_clSmrHeavy, -1, 0.0, 0.0, 2.0,
	false, false, Interface::limited);
	static Parameter<ClusterDecayer,double>
	interfaceClSmrExotic ("ClSmrExotic", "cluster direction Gaussian smearing for exotic quark",
	&ClusterDecayer::_clSmrExotic, 0, 0.0 , 0.0 , 2.0,false,false,false);

	static Switch<ClusterDecayer,bool> interfaceOnShell
	("OnShell",
	"Whether or not the hadrons produced should by on shell or generated using the"
	" mass generator.",
	&ClusterDecayer::_onshell, false, false, false);
	static SwitchOption interfaceOnShellOnShell
	(interfaceOnShell,
	"Yes",
	"Produce the hadrons on shell",
	true);
	static SwitchOption interfaceOnShellOffShell
	(interfaceOnShell,
	"No",
	"Generate the masses using the mass generator.",
	false);

	static Parameter<ClusterDecayer,unsigned int> interfaceMassTry
	("MassTry",
	"The number attempts to generate the masses of the hadrons produced"
	" in the cluster decay.",
	&ClusterDecayer::_masstry, 20, 1, 50,
	false, false, Interface::limited);

	}


	void ClusterDecayer::decay(const ClusterVector & clusters, tPVector & finalhadrons) {
	// Loop over all clusters, and if they are not too heavy (that is
	// intermediate clusters that have undergone to fission) or not
	// too light (that is final clusters that have been already decayed
	// into single hadron) then decay them into two hadrons.

	for (ClusterVector::const_iterator it = clusters.begin();
	it != clusters.end(); ++it) {
	if ((it)->isAvailable() && !(it)->isStatusFinal()
	&& (*it)->isReadyToDecay()) {
	pair<PPtr,PPtr> prod = decayIntoTwoHadrons(*it);
	if(!prod.first)
	throw Exception() << "Can't perform decay of cluster ("
	<< (*it)->particle(0)->dataPtr()->PDGName() << " "
	<< (*it)->particle(1)->dataPtr()->PDGName() << ")\nThreshold = "
	<< ounit(spectrum()->massLightestHadronPair((*it)->particle(0)->dataPtr(),
	(*it)->particle(1)->dataPtr()),GeV)
	<< "\nLightest hadron pair = ( "
	<< spectrum()->lightestHadronPair((*it)->particle(0)->dataPtr(),
	(*it)->particle(1)->dataPtr()).first->PDGName() << ",\t"
	<< spectrum()->lightestHadronPair((*it)->particle(0)->dataPtr(),
	(*it)->particle(1)->dataPtr()).second->PDGName() << " )\nCluster =\n"

	<< **it
	<< "in ClusterDecayer::decay()"
	<< Exception::eventerror;
	(*it)->addChild(prod.first);
	(*it)->addChild(prod.second);
	finalhadrons.push_back(prod.first);
	finalhadrons.push_back(prod.second);
	}
	}
	}

	static double sampleCosTchannel(double Ainv);
	static double sampleCosTchannel(double Ainv){
	/*
	* samples exactly from distribution 1/(1/Ainv-cosTheta)^2
	* where cosTheta in [-1:1]
	* */
	if (fabs(Ainv)<1e-14) return UseRandom::rnd()*2.0-1.0;
	double A=1.0/Ainv;
	double r=UseRandom::rnd();
	double cosTheta = r*2.0-1.0;

	return (1.0+A*cosTheta)/(A+cosTheta);
	}
	/*
	static double sampleCosRegular(double A, double R);
	static double sampleCosRegular(double A, double R){
	double zNeg=(A-1)/R;
	double zPos=(A+1)/R;
	double r=atan(zPos)+UseRandom::rnd()*(atan(zNeg)-atan(zPos));
	double cosTheta=A-R*tan(r);
	if (fabs(cosTheta)>1.0 \|\| std::isnan(cosTheta)){
	std::cout << "cosTheta = "<< cosTheta << std::endl;
	std::cout << "R = "<< R << std::endl;
	std::cout << "r = "<< r << std::endl;
	std::cout << "zNeg = "<< zNeg << std::endl;
	std::cout << "zPos = "<< zPos << std::endl;
	std::cout << "A = "<< A << std::endl;
	}
	return cosTheta;
	}
	*/
	pair<PPtr,PPtr> ClusterDecayer::decayIntoTwoHadrons(tClusterPtr ptr) {
	switch (abs(_kinematics))
	{
	case 0:
	// Default
	return decayIntoTwoHadronsDefault(ptr);
	case 1:
	// Tchannel- like
	return decayIntoTwoHadronsTchannel(ptr);
	case 2:
	// AllIsotropic
	return decayIntoTwoHadronsTchannel(ptr);
	case 3:
	// AllAligned
	return decayIntoTwoHadronsTchannel(ptr);
	default:
	assert(false);
	}
	}

	pair<PPtr,PPtr> ClusterDecayer::decayIntoTwoHadronsDefault(tClusterPtr ptr) {
	using Constants::pi;
	using Constants::twopi;
	// To decay the cluster into two hadrons one must distinguish between
	// constituent quarks (or diquarks) that originate from perturbative
	// processes (hard process or parton shower) from those that are generated
	// by the non-perturbative gluon splitting or from fission of heavy clusters.
	// In the latter case the two body decay is assumed to be isotropic.
	// In the former case instead, if proper flag are activated, the two body
	// decay is assumed to "remember" the direction of the constituents inside
	// the cluster, in the cluster frame. The actual smearing of the hadron
	// directions around the direction of the constituents, in the cluster
	// frame, can be different between non-b hadrons and b-hadrons, but is given
	// by the same functional form:
	// cosThetaSmearing = 1 + smearFactor * log( rnd() )
	// (repeated until cosThetaSmearing > -1)
	// where the factor smearFactor is different between b- and non-b hadrons.
	//
	// We need to require (at least at the moment, maybe in the future we
	// could change it) that the cluster has exactly two components.
	// If this is not the case, then send a warning because it is not suppose
	// to happen, and then return.
	if ( ptr->numComponents() != 2 ) {
	generator()->logWarning( Exception("ClusterDecayer::decayIntoTwoHadrons "
	"*Still cluster with not exactly 2 components* ",
	Exception::warning) );
	return pair<PPtr,PPtr>();
	}

	// Extract the id and particle pointer of the two components of the cluster.
	tPPtr ptr1 = ptr->particle(0);
	tPPtr ptr2 = ptr->particle(1);
	tcPDPtr ptr1data = ptr1->dataPtr();
	tcPDPtr ptr2data = ptr2->dataPtr();

	bool isHad1FlavSpecial = false;
	bool cluDirHad1 = _clDirLight;
	double cluSmearHad1 = _clSmrLight;
	bool isHad2FlavSpecial = false;
	bool cluDirHad2 = _clDirLight;
	double cluSmearHad2 = _clSmrLight;

	if (spectrum()->isExotic(ptr1data)) {
	isHad1FlavSpecial = true;
	cluDirHad1 = _clDirExotic;
	cluSmearHad1 = _clSmrExotic;
	} else {
	for ( const long& id : spectrum()->heavyHadronizingQuarks() ) {
	if ( spectrum()->hasHeavy(id,ptr1data) ) {
	cluDirHad1 = _clDirHeavy[id];
	cluSmearHad1 = _clSmrHeavy[id];
	}
	}
	}


	if (spectrum()->isExotic(ptr2data)) {
	isHad2FlavSpecial = true;
	cluDirHad2 = _clDirExotic;
	cluSmearHad2 = _clSmrExotic;
	} else {
	for ( const long& id : spectrum()->heavyHadronizingQuarks() ) {
	if ( spectrum()->hasHeavy(id,ptr2data) ) {
	cluDirHad2 = _clDirHeavy[id];
	cluSmearHad2 = _clSmrHeavy[id];
	}
	}
	}

	bool isOrigin1Perturbative = ptr->isPerturbative(0);
	bool isOrigin2Perturbative = ptr->isPerturbative(1);

	// We have to decide which, if any, of the two hadrons will have
	// the momentum, in the cluster parent frame, smeared around the
	// direction of its constituent (for Had1 is the one pointed by
	// ptr1, and for Had2 is the one pointed by ptr2).
	// This happens only if the flag _ClDirX is 1 and the constituent is
	// perturbative (that is not coming from nonperturbative gluon splitting
	// or cluster fission). In the case that both the hadrons satisfy this
	// two requirements (of course only one must be treated, because the other
	// one will have the momentum automatically fixed by the momentum
	// conservation) then more priority is given in the case of a b-hadron.
	// Finally, in the case that the two hadrons have same priority, then
	// we choose randomly, with equal probability, one of the two.

	int priorityHad1 = 0;
	if ( cluDirHad1 && isOrigin1Perturbative ) {
	priorityHad1 = isHad1FlavSpecial ? 2 : 1;
	}
	int priorityHad2 = 0;
	if ( cluDirHad2 && isOrigin2Perturbative ) {
	priorityHad2 = isHad2FlavSpecial ? 2 : 1;
	}
	if ( priorityHad2 && priorityHad1 == priorityHad2 && UseRandom::rndbool() ) {
	priorityHad2 = 0;
	}

	Lorentz5Momentum pClu = ptr->momentum();
	bool secondHad = false;
	Axis uSmear_v3;
	if ( priorityHad1 \|\| priorityHad2 ) {

	double cluSmear;
	Lorentz5Momentum pQ;
	Lorentz5Momentum pQbar;
	if ( priorityHad1 > priorityHad2 ) {
	pQ = ptr1->momentum();
	pQbar = ptr2->momentum();
	cluSmear = cluSmearHad1;
	} else {
	pQ = ptr2->momentum();
	pQbar = ptr1->momentum();
	cluSmear = cluSmearHad2;
	secondHad = true;
	}

	// To find the momenta of the two hadrons in the parent cluster frame
	// we proceed as follows. First, we determine the unit vector parallel
	// to the direction of the constituent in the cluster frame. Then we
	// have to smear such direction using the following prescription:
	// --- in theta angle w.r.t. such direction (not the z-axis),
	// the drawing of the cosine of such angle is done via:
	// 1.0 + cluSmear*log( rnd() )
	// (repeated until it gives a value above -1.0)
	// --- in phi angle w.r.t. such direction, the drawing is simply flat.
	// Then, given the direction in the parent cluster frame of one of the
	// two hadrons, it is straighforward to get the momenta of both hadrons
	// (always in the same parent cluster frame).

	if (_covariantBoost) {
	// Use correct boost for a 2 particle constituent cluster
	// the direction of the constituents is aligned with the z-Axis
	// since the boost is constructed as such
	Energy Pcom = Kinematics::pstarTwoBodyDecay(pClu.m(),pQ.m(),pQbar.m());
	pQ = secondHad ? Lorentz5Momentum(pQbar.m(), -PcomAxis(0,0,1)):Lorentz5Momentum(pQ.m(), PcomAxis(0,0,1));
	// pQbar = Lorentz5Momentum(pQbar.m(),-Pcom*Axis(0,0,1));
	}
	else {
	pQ.boost( -pClu.boostVector() ); // naive boost from Lab to Cluster frame
	}
	uSmear_v3 = pQ.vect().unit(); // Direction of the priority quark in COM frame
	// skip if cluSmear is too small
	if ( cluSmear > 0.001 ) {
	// generate the smearing angle
	double cosSmear;
	do cosSmear = 1.0 + cluSmear*log( UseRandom::rnd() );
	while ( cosSmear < -1.0 );
	double sinSmear = sqrt( 1.0 - sqr(cosSmear) );
	// calculate rotation to z axis
	LorentzRotation rot;
	double sinth(sqrt(1.-sqr(uSmear_v3.z())));
	if(abs(uSmear_v3.perp2()/uSmear_v3.z())>1e-10)
	rot.setRotate(-acos(uSmear_v3.z()),
	Axis(-uSmear_v3.y()/sinth,uSmear_v3.x()/sinth,0.));
	// + random azimuthal rotation
	rot.rotateZ(UseRandom::rnd()*twopi);
	// set direction in rotated frame
	Lorentz5Vector<double> ltemp(0.,sinSmear,cosSmear,0.);
	// rotate back
	rot.invert();
	ltemp *= rot;
	uSmear_v3 = ltemp.vect();
	}
	}
	else {

	// Isotropic decay: flat in cosTheta and phi.
	uSmear_v3 = Axis(1.0, 0.0, 0.0); // just to set the rho to 1
	uSmear_v3.setTheta( acos( UseRandom::rnd( -1.0 , 1.0 ) ) );
	uSmear_v3.setPhi( UseRandom::rnd( -pi , pi ) );

	}

	pair<tcPDPtr,tcPDPtr> dataPair
	= _hadronSpectrum->chooseHadronPair(ptr->mass(),
	ptr1data,
	ptr2data);
	if(dataPair.first == tcPDPtr() \|\|
	dataPair.second == tcPDPtr()) return pair<PPtr,PPtr>();

	// Create the two hadron particle objects with the specified id.
	PPtr ptrHad1,ptrHad2;
	// produce the hadrons on mass shell
	if(_onshell) {
	ptrHad1 = dataPair.first ->produceParticle(dataPair.first ->mass());
	ptrHad2 = dataPair.second->produceParticle(dataPair.second->mass());
	}
	// produce the hadrons with mass given by the mass generator
	else {
	unsigned int ntry(0);
	do {
	ptrHad1 = dataPair.first ->produceParticle();
	ptrHad2 = dataPair.second->produceParticle();
	++ntry;
	}
	while(ntry<_masstry&&ptrHad1->mass()+ptrHad2->mass()>ptr->mass());
	// if fails produce on shell and issue warning (should never happen??)
	if( ptrHad1->mass() + ptrHad2->mass() > ptr->mass() ) {
	generator()->log() << "Failed to produce off-shell hadrons in "
	<< "ClusterDecayer::decayIntoTwoHadrons producing hadrons "
	<< "on-shell" << endl;
	ptrHad1 = dataPair.first ->produceParticle(dataPair.first ->mass());
	ptrHad2 = dataPair.second->produceParticle(dataPair.second->mass());
	}
	}

	if (!ptrHad1 \|\| !ptrHad2) {
	ostringstream s;
	s << "ClusterDecayer::decayIntoTwoHadrons *Cannot create the two hadrons*"
	<< dataPair.first ->PDGName() << " and "
	<< dataPair.second->PDGName();
	cerr << s.str() << endl;
	generator()->logWarning( Exception(s.str(), Exception::warning) );
	} else {

	Lorentz5Momentum pHad1, pHad2; // 5-momentum vectors that we have to determine
	if ( secondHad ) uSmear_v3 *= -1.0;

	if (pClu.m() < ptrHad1->mass()+ptrHad2->mass() ) {
	throw Exception() << "Impossible Kinematics in ClusterDecayer::decayIntoTwoHadrons()"
	<< Exception::eventerror;
	}
	if (_covariantBoost) {
	// Use correct boost for a 2 particle constituent cluster
	Energy Pcom=Kinematics::pstarTwoBodyDecay(pClu.m(),ptrHad1->mass(),ptrHad2->mass());
	pHad1=Lorentz5Momentum(ptrHad1->mass(), Pcom*uSmear_v3);
	pHad2=Lorentz5Momentum(ptrHad2->mass(),-Pcom*uSmear_v3);
	const Lorentz5Momentum pQ1(ptr1data->constituentMass(),ptr->particle(0)->momentum().vect());
	const Lorentz5Momentum pQ2(ptr2data->constituentMass(),ptr->particle(1)->momentum().vect());
	- Kinematics::BoostIntoTwoParticleFrame(pClu.m(),pQ1, pQ2, pHad1, pHad2);
	+ std::vector<Lorentz5Momentum *> momenta;
	+ momenta.push_back(&pHad1);
	+ momenta.push_back(&pHad2);
	+ Kinematics::BoostIntoTwoParticleFrame(pClu.m(),pQ1, pQ2, momenta);
	}
	else {
	// 5-momentum vectors that we have to determine
	Kinematics::twoBodyDecay(pClu,ptrHad1->mass(),ptrHad2->mass(),uSmear_v3,
	pHad1,pHad2);
	}
	ptrHad1->set5Momentum(pHad1);
	ptrHad2->set5Momentum(pHad2);

	// Determine the positions of the two children clusters.
	LorentzPoint positionHad1 = LorentzPoint();
	LorentzPoint positionHad2 = LorentzPoint();
	calculatePositions(pClu, ptr->vertex(), pHad1, pHad2, positionHad1, positionHad2);
	ptrHad1->setVertex(positionHad1);
	ptrHad2->setVertex(positionHad2);
	}
	return pair<PPtr,PPtr>(ptrHad1,ptrHad2);
	}

	pair<PPtr,PPtr> ClusterDecayer::decayIntoTwoHadronsTchannel(tClusterPtr ptr) {
	using Constants::pi;
	using Constants::twopi;
	/* New test for cluster decay modelling according to matrix element
	* C-> q,qbar->h1,h2 using \|M\|^2=1/((p1-p3)^2-(m1-m3)^2)^2 which is Tchannel-like
	* */
	if ( ptr->numComponents() != 2 ) {
	generator()->logWarning( Exception("ClusterDecayer::decayIntoTwoHadrons "
	"*Still cluster with not exactly 2 components* ",
	Exception::warning) );
	return pair<PPtr,PPtr>();
	}

	// Extract the id and particle pointer of the two components of the cluster.
	tPPtr ptr1 = ptr->particle(0);
	tPPtr ptr2 = ptr->particle(1);
	tcPDPtr ptr1data = ptr1->dataPtr();
	tcPDPtr ptr2data = ptr2->dataPtr();

	Lorentz5Momentum pClu = ptr->momentum();

	pair<tcPDPtr,tcPDPtr> dataPair
	= _hadronSpectrum->chooseHadronPair(ptr->mass(),
	ptr1data,
	ptr2data);
	if(dataPair.first == tcPDPtr() \|\|
	dataPair.second == tcPDPtr()) return pair<PPtr,PPtr>();

	// Create the two hadron particle objects with the specified id.
	PPtr ptrHad1,ptrHad2;
	// produce the hadrons on mass shell
	if(_onshell) {
	ptrHad1 = dataPair.first ->produceParticle(dataPair.first ->mass());
	ptrHad2 = dataPair.second->produceParticle(dataPair.second->mass());
	}
	// produce the hadrons with mass given by the mass generator
	else {
	unsigned int ntry(0);
	do {
	ptrHad1 = dataPair.first ->produceParticle();
	ptrHad2 = dataPair.second->produceParticle();
	++ntry;
	}
	while(ntry<_masstry&&ptrHad1->mass()+ptrHad2->mass()>ptr->mass());
	// if fails produce on shell and issue warning (should never happen??)
	if( ptrHad1->mass() + ptrHad2->mass() > ptr->mass() ) {
	generator()->log() << "Failed to produce off-shell hadrons in "
	<< "ClusterDecayer::decayIntoTwoHadrons producing hadrons "
	<< "on-shell" << endl;
	ptrHad1 = dataPair.first ->produceParticle(dataPair.first ->mass());
	ptrHad2 = dataPair.second->produceParticle(dataPair.second->mass());
	}
	}

	if (!ptrHad1 \|\| !ptrHad2) {
	ostringstream s;
	s << "ClusterDecayer::decayIntoTwoHadrons *Cannot create the two hadrons*"
	<< dataPair.first ->PDGName() << " and "
	<< dataPair.second->PDGName();
	cerr << s.str() << endl;
	generator()->logWarning( Exception(s.str(), Exception::warning) );
	} else {

	if (pClu.m() < ptrHad1->mass()+ptrHad2->mass() ) {
	throw Exception() << "Impossible Kinematics in ClusterDecayer::decayIntoTwoHadrons()"
	<< Exception::eventerror;
	}
	Lorentz5Momentum pHad1,pHad2;
	Axis uSmear_v3;
	Lorentz5Momentum pQ = ptr1->momentum();
	Lorentz5Momentum pQbar = ptr2->momentum();

	if (_covariantBoost) {
	// Use correct boost for a 2 particle constituent cluster
	// the direction of the constituents is aligned with the z-Axis
	// since the boost is constructed as such
	Energy Pcom = Kinematics::pstarTwoBodyDecay(pClu.m(),pQ.m(),pQbar.m());
	pQ = Lorentz5Momentum(pQ.m(), Pcom*Axis(0,0,1));
	pQbar = Lorentz5Momentum(pQbar.m(),-Pcom*Axis(0,0,1));
	}
	else {
	pQ.boost( -pClu.boostVector() ); // naive boost from Lab to Cluster frame
	}
	const Axis dirQ = pQ.vect().unit(); // Direction of the quark in COM frame

	switch (abs(_kinematics))
	{
	case 1:
	{
	Energy mHad1 = ptrHad1->mass();
	Energy mHad2 = ptrHad2->mass();
	Energy mConst1 = ptr1data->constituentMass();
	Energy mConst2 = ptr2data->constituentMass();
	Energy PcomIni = Kinematics::pstarTwoBodyDecay(pClu.m(),mConst1,mConst2);
	Energy PcomFin = Kinematics::pstarTwoBodyDecay(pClu.m(),mHad1,mHad2);
	Energy EcomIni1 = sqrt(sqr(mConst1) + sqr(PcomIni));
	Energy EcomFin1 = sqrt(sqr(mHad1) + sqr(PcomFin));
	// double Ainv = PcomIniPcomFin/(EcomIni1EcomFin1);
	// New sampling from distribution 1/((p1-p3)^2-(m1-mHad1)^2)
	double Ainv = PcomIniPcomFin/(EcomIni1EcomFin1-mHad1*mConst1);
	double cosTheta,phi;
	do {
	uSmear_v3 = dirQ;
	cosTheta = sampleCosTchannel(Ainv);
	if (fabs(cosTheta-1.0)>1e-14) {
	// rotate to sampled angles
	uSmear_v3.rotate(acos(cosTheta),dirQ.orthogonal());
	phi = UseRandom::rnd(-pi,pi);
	uSmear_v3.rotate(phi,dirQ);
	}
	}
	// filter out not corresponding hemisphere
	while ( _kinematics<0 && dirQ.dot(uSmear_v3)<=0.0);
	break;
	}
	case 2:
	{
	// uSmear_v3 changed to isotropic direction
	Axis iso;
	do {
	iso = Axis(1.0, 0.0, 0.0); // just to set the rho to 1
	iso.setTheta( acos( UseRandom::rnd( -1.0 , 1.0 ) ) );
	iso.setPhi( UseRandom::rnd( -pi , pi ) );
	}
	// filter out not corresponding hemisphere
	while ( _kinematics<0 && dirQ.dot(iso)<=0.0);
	uSmear_v3 = iso;
	break;
	}
	case 3:
	{
	// set to be aligned
	uSmear_v3 = dirQ;
	break;
	}
	default:
	assert(false);
	}
	// sanity check
	if (_kinematics<0 && !(uSmear_v3.dot(dirQ)>0.0)) {
	std::cout << "uSmear_v3.dirQ "<< uSmear_v3.dot(dirQ) << std::endl;
	}
	assert(_kinematics>0 \|\| (_kinematics<0 && uSmear_v3.dot(dirQ)>0.0));
	// TODO simple Tchannel 1/(pConst1-pHad1)**4 distribution
	/*
	//sanity checks
	static int firstAsmall=1;
	static int firstAlarge=1;
	if (firstAsmall && fabs(A)<1){
	// Histogram hcosTheta(-1,1,100);
	Histogram hcosThetaRegular(-1,1,100);
	int N=10000;
	double R=1;
	std::cout << "\nA = "<<A <<"\n";
	std::cout << "\nR = "<<R <<"\n";
	for (int i = 0; i < N; i++)
	{
	double c=sampleCosRegular(A,R);
	hcosThetaRegular+=c;
	}
	ofstream out1("hcosThetaRegular_Asmall.dat", std::ios::out);
	hcosThetaRegular.topdrawOutput(out1);
	out1.close();
	ofstream out2("hcosThetaRegular_Asmall_RIVET.dat", std::ios::out);
	hcosThetaRegular.rivetOutput(out2);
	out2.close();
	firstAsmall=0;
	}
	static int firstAlarge=1;
	Ainv=0.8;
	if (firstAlarge && fabs(1./Ainv)>1){
	Histogram hcosThetaTchannel(-1,1,100);
	double A=1.0/Ainv;
	std::cout << "\nA = "<<A <<"\n";
	int N=10000;
	for (int i = 0; i < N; i++)
	{
	double c=sampleCosTchannel(Ainv);
	hcosThetaTchannel+=c;
	}
	ofstream out1("hcosThetaTchannel_Alarge.dat", std::ios::out);
	hcosThetaTchannel.topdrawOutput(out1);
	out1.close();
	ofstream out2("hcosThetaTchannel_Alarge_RIVET.dat", std::ios::out);
	hcosThetaTchannel.rivetOutput(out2);
	out2.close();
	firstAlarge=0;
	}
	*/
	if (_covariantBoost) {
	// Use correct boost for a 2 particle constituent cluster
	Energy PcomFin=Kinematics::pstarTwoBodyDecay(pClu.m(),ptrHad1->mass(),ptrHad2->mass());
	pHad1=Lorentz5Momentum(ptrHad1->mass(), PcomFin*uSmear_v3);
	pHad2=Lorentz5Momentum(ptrHad2->mass(),-PcomFin*uSmear_v3);
	const Lorentz5Momentum pQ1(ptr1data->constituentMass(),ptr->particle(0)->momentum().vect());
	const Lorentz5Momentum pQ2(ptr2data->constituentMass(),ptr->particle(1)->momentum().vect());
	- Kinematics::BoostIntoTwoParticleFrame(pClu.m(),pQ1, pQ2, pHad1, pHad2);
	+ std::vector<Lorentz5Momentum *> momenta;
	+ momenta.push_back(&pHad1);
	+ momenta.push_back(&pHad2);
	+ Kinematics::BoostIntoTwoParticleFrame(pClu.m(),pQ1, pQ2, momenta);
	}
	else {
	// 5-momentum vectors that we have to determine
	Kinematics::twoBodyDecay(pClu,ptrHad1->mass(),ptrHad2->mass(),uSmear_v3,
	pHad1,pHad2);
	}
	ptrHad1->set5Momentum(pHad1);
	ptrHad2->set5Momentum(pHad2);

	// Determine the positions of the two children clusters.
	LorentzPoint positionHad1 = LorentzPoint();
	LorentzPoint positionHad2 = LorentzPoint();
	calculatePositions(pClu, ptr->vertex(), pHad1, pHad2, positionHad1, positionHad2);
	ptrHad1->setVertex(positionHad1);
	ptrHad2->setVertex(positionHad2);
	}
	return pair<PPtr,PPtr>(ptrHad1,ptrHad2);
	}


	void ClusterDecayer::
	calculatePositions(const Lorentz5Momentum &pClu,
	const LorentzPoint &positionClu,
	const Lorentz5Momentum &,
	const Lorentz5Momentum &,
	LorentzPoint &positionHad1,
	LorentzPoint &positionHad2 ) const {
	// First, determine the relative positions of the children hadrons
	// with respect to their parent cluster, in the cluster reference frame,
	// assuming gaussian smearing with width inversely proportional to the
	// parent cluster mass.
	Length smearingWidth = hbarc / pClu.m();
	LorentzDistance distanceHad[2];
	for ( int i = 0; i < 2; i++ ) { // i==0 is the first hadron; i==1 is the second one
	Length delta[4]={ZERO,ZERO,ZERO,ZERO};
	// smearing of the four components of the LorentzDistance, two at the same time to improve speed
	for ( int j = 0; j < 3; j += 2 ) {
	delta[j] = UseRandom::rndGauss(smearingWidth, Length(ZERO));
	delta[j+1] = UseRandom::rndGauss(smearingWidth, Length(ZERO));
	}
	// set the distance
	delta[0] = abs(delta[0]) +sqrt(sqr(delta[1])+sqr(delta[2])+sqr(delta[3]));
	distanceHad[i] = LorentzDistance(delta[1],delta[2],delta[3],delta[0]);
	// Boost such relative positions of the children hadrons,
	// with respect to their parent cluster,
	// from the cluster reference frame to the Lab frame.
	distanceHad[i].boost(pClu.boostVector());
	}
	// Finally, determine the absolute positions of the children hadrons
	// in the Lab frame.
	positionHad1 = distanceHad[0] + positionClu;
	positionHad2 = distanceHad[1] + positionClu;
	}
	diff --git a/Utilities/Kinematics.cc b/Utilities/Kinematics.cc
	--- a/Utilities/Kinematics.cc
	+++ b/Utilities/Kinematics.cc
	@@ -1,545 +1,514 @@
	// -- C++ --
	//
	// Kinematics.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 Kinematics class.
	//

	#include "Kinematics.h"
	#include <ThePEG/Vectors/Lorentz5Vector.h>
	#include <ThePEG/Vectors/LorentzVector.h>
	#include <ThePEG/Vectors/LorentzRotation.h>
	#include <ThePEG/Repository/EventGenerator.h>
	#include <ThePEG/Repository/CurrentGenerator.h>
	#include <ThePEG/EventRecord/Event.h>

	#include <boost/numeric/ublas/matrix.hpp>
	#include <boost/numeric/ublas/io.hpp>
	#include <boost/numeric/ublas/lu.hpp>


	using namespace Herwig;
	using namespace std;
	using namespace ThePEG;

	/**
	* Boost consistently the Lorentz5Momenta in momenta (given in the pi and pj COM frame)
	* into the p0Q1 p0Q2 frame.
	* */
	void Kinematics::BoostIntoTwoParticleFrame(const Energy M, const Lorentz5Momentum & piLab,
	const Lorentz5Momentum & pjLab,
	std::vector<Lorentz5Momentum * > momenta) {
	double cosPhi=piLab.vect().cosTheta(pjLab.vect());
	double Phi=acos(cosPhi);
	double sinPhi=sin(Phi);
	// If Phi==0 use regular 1+1D boost
	double epsilon=std::numeric_limits<double>::epsilon();
	if (fabs(cosPhi-1.0)<=epsilon \|\| fabs(cosPhi+1.0)<=epsilon) {
	Lorentz5Momentum pClu(M,(piLab+pjLab).vect());
	Boost bv = pClu.boostVector();
	for (unsigned int it = 0; it < momenta.size(); it++) {
	momenta[it]->boost(bv);
	}
	return;
	}
	Energy Ei=piLab.e();
	Energy Ej=pjLab.e();
	if (std::isnan(Phi) \|\| std::isinf(Phi)) throw Exception() << "NAN or INF in Phi in Kinematics::BoostIntoTwoParticleFrame\n"
	<< Exception::runerror;
	Energy mi=piLab.mass();
	Energy mj=pjLab.mass();
	Energy2 mi2=mi*mi;
	Energy2 mj2=mj*mj;
	Energy Pi=piLab.vect().mag();
	Energy Pj=pjLab.vect().mag();
	assert(Pi>ZERO);

	std::vector<boost::numeric::ublas::vector<double>> momentaHat;

	std::vector<Energy> Masses;
	for (unsigned int it = 0; it < momenta.size(); it++) {
	Masses.push_back(momenta[it]->mass());
	momentaHat.push_back(boost::numeric::ublas::vector<double>(3));
	momentaHat[it](0) = momenta[it]->e()/GeV;
	momentaHat[it](1) = momenta[it]->x()/GeV;
	momentaHat[it](2) = momenta[it]->z()/GeV;
	}


	// Lorentz Matrix Lambda maps:
	// piHat = (Ecomi,0,0, Pcom) to piLab = (Ei, 0,0,Pi )
	// pjHat = (Ecomj,0,0,-Pcom) to pjLab = (Ej,Pjsin(phi),0,Pjcos(phi))
	// and therefore maps also correctly momentaHat to momentaOut into the Lab frame
	boost::numeric::ublas::matrix<double> Lambda(3,3);

	Energy pstar=Kinematics::pstarTwoBodyDecay(M,mi,mj);
	Energy2 A2=pstar*M;
	Energy2 B2=PiPjsinPhi;
	Energy B2divPi=Pj*sinPhi;
	Energy2 Deltaij=PiEj-EiPj*cosPhi;
	double delta2=mi2PjPjsinPhisinPhi/(Deltaij*Deltaij);
	double Lambda11=0;

	// better numerics
	if (delta2<1e-13) Lambda11 = Deltaij>=ZERO ? (1.0-0.5delta2):-(1.0-0.5delta2);
	else if (Deltaij!=ZERO) Lambda11= Deltaij>=ZERO ? 1.0/sqrt(1.0+delta2):-1.0/sqrt(1.0+delta2);


	if (std::isnan(A2/GeV2) \|\| std::isinf(A2/GeV2)) throw Exception() << "NAN in A2/GeV2\n"
	<< Exception::runerror;

	Lambda(0,0) = (Ei+Ej)/M;
	Lambda(0,1) = B2/A2;
	Lambda(0,2) = (Ei-Ej)/(2.0pstar)-((mi2-mj2)(Ei+Ej))/(2.0MA2);

	Lambda(1,0) = B2divPi/M;
	- Lambda(1,1) = Lambda11;
	+ Lambda(1,1) = Lambda11; // This should be just Deltaij/(M*pstar)
	Lambda(1,2) = -(MM-(mj2-mi2))B2divPi/(2.0MA2);

	Lambda(2,0) = (Pi+Pj*cosPhi)/M;
	Lambda(2,1) = Ei*B2divPi/A2;
	Lambda(2,2) = (A2A2-0.5Ei(Ej(MM-(mj2-mi2))-Ei(MM-(mi2-mj2))))/(PiM*A2);

	Axis zAxis(0,0,1);
	Axis xAxis(1,0,0);
	Lorentz5Momentum piRes(mi,Pi*zAxis);
	Lorentz5Momentum pjRes(mj,Pj(xAxissinPhi+zAxis*cosPhi));

	std::vector<Lorentz5Momentum> momentaRes;
	Lorentz5Momentum pClu1,pClu2;
	- // TODO FIX THIS ERROR consistently
	- // Horrible fix below but works partially
	- // throw Exception() << "Phi = "<< std::setprecision(16)<<fabs(pClu1.vect()*Axis(0,0,1)/pClu1.vect().mag())-1.0<<" in Kinematics::BoostIntoTwoParticleFrame\n" << Exception::warning;
	- // TODO in this case just rescale piLab pjLab correspondingly, but the directions shall not change
	- // pClu1 aligned with piLab
	boost::numeric::ublas::vector<double> momentaOut(3);
	- /*
	- if (
	- fabs(momenta[0]->x()/GeV) < 1e-8
	- && fabs(momenta[0]->y()/GeV) < 1e-8
	- && fabs(momenta[1]->x()/GeV) < 1e-8
	- && fabs(momenta[1]->y()/GeV) < 1e-8
	- && fabs((momenta[0]->z()+momenta[1]->z())/GeV) < 1e-8
	- ) {
	- isAligned = 2;
	- double factor = fabs(momenta[0]->z()/pstar);
	- double otherFactor[2];
	- if (momenta[0]->z() > ZERO) {
	- otherFactor[0] = (momenta[0]->e()-factor*sqrt(mi2+sqr(pstar)))/M;
	- otherFactor[1] = (momenta[1]->e()-factor*sqrt(mj2+sqr(pstar)))/M;
	- // doing the Boost into the Lab frame analytically:
	- pClu1 = factorpiRes + otherFactor[0](piRes + pjRes);
	- pClu2 = factorpjRes + otherFactor[1](piRes + pjRes);
	- }
	- else {
	- otherFactor[0] = (momenta[1]->e()-factor*sqrt(mi2+sqr(pstar)))/M;
	- otherFactor[1] = (momenta[0]->e()-factor*sqrt(mj2+sqr(pstar)))/M;
	- // doing the Boost into the Lab frame analytically:
	- pClu2 = factorpiRes + otherFactor[0](piRes + pjRes);
	- pClu1 = factorpjRes + otherFactor[1](piRes + pjRes);
	- }
	- momentaRes.push_back(pClu1);
	- momentaRes.push_back(pClu2);
	- } */
	-
	unsigned int iter = 0;
	bool isAligned;
	Momentum3 piHat(ZERO, ZERO, pstar);
	Momentum3 pjHat(ZERO, ZERO, -pstar);
	Momentum3 pAligned(ZERO, ZERO, ZERO);
	+ // TODO FIX THIS ERROR consistently
	+ // Horrible fix below but works partially
	+ // TODO in this case just rescale piLab pjLab correspondingly, but the directions shall not change
	+ // pClu1 aligned with piLab
	for (auto & pHat : momentaHat)
	{
	isAligned = false;
	if (momenta[iter]->vect().mag()>ZERO) {
	if (
	fabs(1.0 - momenta[iter]->vect().cosTheta(piHat)) < 1.0e-14
	) {
	isAligned = true;
	double factor = momenta[iter]->z()/pstar;
	assert(momenta[iter]->z()/pstar > 0);
	double otherFactor = (momenta[iter]->e()-factor*sqrt(mi2+sqr(pstar)))/M;
	// doing the Boost into the Lab frame analytically:
	pAligned = (factorpiRes.vect() + otherFactor(piRes.vect() + pjRes.vect()));
	} else if (
	fabs(1.0 - momenta[iter]->vect().cosTheta(pjHat)) < 1.0e-14
	) {
	isAligned = true;
	double factor = fabs(momenta[iter]->z()/pstar);
	double otherFactor = (momenta[iter]->e()-factor*sqrt(mj2+sqr(pstar)))/M;
	// doing the Boost into the Lab frame analytically:
	pAligned = (factorpjRes.vect() + otherFactor(piRes.vect() + pjRes.vect()));
	}
	}
	// doing the Boost into the Lab frame:
	if (!isAligned) {
	momentaOut = boost::numeric::ublas::prod(Lambda,pHat);
	momentaRes.push_back(Lorentz5Momentum(Masses[iter],
	GeV*Axis(momentaOut(1), double(momenta[iter]->y()/GeV), momentaOut(2))));
	}
	else {
	momentaRes.push_back(Lorentz5Momentum(Masses[iter], pAligned));
	}
	iter++;
	}
	// Computing the correct rotation, which maps pi/jRes into pi/jLab
	Axis omega1=piRes.vect().unit().cross(piLab.vect().unit());
	double cosAngle1=piRes.vect().unit()*piLab.vect().unit();
	double angle1=acos(cosAngle1);

	// Rotate piRes into piLab
	piRes.rotate(angle1, omega1);
	pjRes.rotate(angle1, omega1);
	// Correspondingly do the actual rotation on all momenta
	for(auto & pRes : momentaRes)
	pRes.rotate(angle1, omega1);

	Axis omega2=piRes.vect().unit();
	Momentum3 r1dim=(pjLab.vect()-piRes.vect().unit()(pjLab.vect()piRes.vect().unit()));
	Momentum3 r2dim=(pjRes.vect()-piRes.vect().unit()(pjRes.vect()piRes.vect().unit()));

	if (r1dim.mag()==ZERO \|\| r2dim.mag()==ZERO \|\| fabs(sinPhi)<1e-14) //trivial rotation so we are done
	{
	for (unsigned int i = 0; i < momentaRes.size(); i++) {
	// copy the final momenta
	*(momenta[i]) = momentaRes[i];
	}
	return;
	}
	Axis r1=r1dim.unit();
	Axis r2=r2dim.unit();

	// signs for 2nd rotation
	int signToPi = (piRes.vect()*pjLab.vect())/GeV2 > 0 ? 1:-1;
	int signToR1R2 = signToPi(r2.unit().cross(r1.unit())piRes.vect())/GeV> 0 ? 1:-1;
	double angle2=acos(r1.unit()*r2.unit());
	if (signToR1R2<0) angle2=-angle2;

	// Rotate pjRes into pjLab
	pjRes.rotate(angle2, signToPi*omega2);
	// Correspondingly do the actual rotation on all momenta
	for (unsigned int i = 0; i < momentaRes.size(); i++) {
	momentaRes[i].rotate(angle2, signToPi*omega2);
	// copy the final momenta
	*(momenta[i]) = momentaRes[i];
	}
	}
	-
	/**
	* Boost consistently the Lorentz5Momenta pClu1 and pClu2 (given in their COM frame)
	* into the p0Q1 p0Q2 frame.
	* */
	+/*
	void Kinematics::BoostIntoTwoParticleFrame(const Energy M, const Lorentz5Momentum & piLab,
	const Lorentz5Momentum & pjLab,
	Lorentz5Momentum & pClu1,
	Lorentz5Momentum & pClu2) {
	double cosPhi=piLab.vect().cosTheta(pjLab.vect());
	double Phi=acos(cosPhi);
	double sinPhi=sin(Phi);
	// If Phi==0 use regular 1+1D boost
	double epsilon=std::numeric_limits<double>::epsilon();
	if (fabs(cosPhi-1.0)<=epsilon \|\| fabs(cosPhi+1.0)<=epsilon) {
	Lorentz5Momentum pClu(M,(piLab+pjLab).vect());
	Boost bv = pClu.boostVector();
	pClu1.boost(bv);
	pClu2.boost(bv);
	return;
	}
	Energy Ei=piLab.e();
	Energy Ej=pjLab.e();
	if (std::isnan(Phi) \|\| std::isinf(Phi)) throw Exception() << "NAN or INF in Phi in Kinematics::BoostIntoTwoParticleFrame\n"
	<< Exception::runerror;
	Energy mi=piLab.mass();
	Energy mj=pjLab.mass();
	Energy2 mi2=mi*mi;
	Energy2 mj2=mj*mj;
	Energy Pi=piLab.vect().mag();
	Energy Pj=pjLab.vect().mag();
	assert(Pi>ZERO);

	boost::numeric::ublas::vector<double> Pclu1Hat(3);
	boost::numeric::ublas::vector<double> Pclu2Hat(3);

	const Energy M1 = pClu1.mass();
	const Energy M2 = pClu2.mass();

	Pclu1Hat(0)=pClu1.e()/GeV;
	Pclu1Hat(1)=pClu1.x()/GeV;
	Pclu1Hat(2)=pClu1.z()/GeV;

	Pclu2Hat(0)=pClu2.e()/GeV;
	Pclu2Hat(1)=pClu2.x()/GeV;
	Pclu2Hat(2)=pClu2.z()/GeV;

	// Lorentz Matrix Lambda maps:
	// piHat = (Ecomi,0,0, Pcom) to piLab = (Ei, 0,0,Pi )
	// pjHat = (Ecomj,0,0,-Pcom) to pjLab = (Ej,Pjsin(phi),0,Pjcos(phi))
	// and therefore maps also correctly Pclu1/2Hat to Pclu1/2 into the Lab frame
	boost::numeric::ublas::matrix<double> Lambda(3,3);

	Energy pstar=Kinematics::pstarTwoBodyDecay(M,mi,mj);
	Energy2 A2=pstar*M;
	Energy2 B2=PiPjsinPhi;
	Energy B2divPi=Pj*sinPhi;
	Energy2 Deltaij=PiEj-EiPj*cosPhi;
	double delta2=mi2PjPjsinPhisinPhi/(Deltaij*Deltaij);
	double Lambda11=0;

	// better numerics
	if (delta2<1e-13) Lambda11 = Deltaij>=ZERO ? (1.0-0.5delta2):-(1.0-0.5delta2);
	else if (Deltaij!=ZERO) Lambda11= Deltaij>=ZERO ? 1.0/sqrt(1.0+delta2):-1.0/sqrt(1.0+delta2);


	if (std::isnan(A2/GeV2) \|\| std::isinf(A2/GeV2)) throw Exception() << "NAN in A2/GeV2\n"
	<< Exception::runerror;

	Lambda(0,0) = (Ei+Ej)/M;
	Lambda(0,1) = B2/A2;
	Lambda(0,2) = (Ei-Ej)/(2.0pstar)-((mi2-mj2)(Ei+Ej))/(2.0MA2);

	Lambda(1,0) = B2divPi/M;
	Lambda(1,1) = Lambda11;
	Lambda(1,2) = -(MM-(mj2-mi2))B2divPi/(2.0MA2);

	Lambda(2,0) = (Pi+Pj*cosPhi)/M;
	Lambda(2,1) = Ei*B2divPi/A2;
	Lambda(2,2) = (A2A2-0.5Ei(Ej(MM-(mj2-mi2))-Ei(MM-(mi2-mj2))))/(PiM*A2);

	- /* Determinant calculation, but this is often far from 1 due to rounding errors
	- * */
	+ // Determinant calculation, but this is often far from 1 due to rounding errors
	// double det=0;
	// det += Lambda(0,0)Lambda(1,1)Lambda(3,2);
	// det += Lambda(1,0)Lambda(2,1)Lambda(0,2);
	// det += Lambda(2,0)Lambda(0,1)Lambda(1,2);

	// det -= Lambda(0,2)Lambda(1,1)Lambda(2,0);
	// det -= Lambda(1,0)Lambda(0,1)Lambda(2,2);
	// det -= Lambda(1,2)Lambda(2,1)Lambda(0,0);
	// if (fabs(det-1.0)>1e-3 \|\| std::isnan(det)) {
	// std::cout << "A2 = "<<std::setprecision(18)<< A2/(GeV2) << std::endl;
	// std::cout << "B2 = "<< B2/(GeV2)<< std::endl;
	// std::cout << "DET-1 = "<< det-1.0 << std::setprecision(3)<< std::endl;
	// std::cout << "Lambda:" << std::endl;
	// for (int i = 0; i < 3; i++)
	// {
	// for (int j = 0; j < 3; j++)
	// {
	// std::cout<< std::setprecision(18) << Lambda(i,j)<< std::setprecision(3) << "\t";
	// }
	// std::cout << "\n";
	// }
	// std::cout << "pi,pj in Lab" << std::endl;
	// std::cout << "Phi = " << std::setprecision(18) << Phi << std::endl;
	// std::cout << "Phi - pi = " << std::setprecision(18) << Phi - M_PI << std::endl;
	// }
	// TODO FIX THIS ERROR consistently
	// throw Exception() << "Phi = "<< std::setprecision(16)<<fabs(pClu1.vect()*Axis(0,0,1)/pClu1.vect().mag())-1.0<<" in Kinematics::BoostIntoTwoParticleFrame\n" << Exception::warning;
	// TODO in this case just rescale piLab pjLab correspondingly, but the directions shall not change
	// if ( pClu1.vect()*Axis(0,0,1)/pClu1.vect().mag() > 0) {
	// pClu1 aligned with piLab
	// double factor = pClu1;
	// Lorentz5Momenta pClu1Lab(pClu1.mass(),piLab.vect()*factor);
	// Lorentz5Momentum pClu(M,piLab+pjLab);
	// Boost bv = pClu.boostVector();
	// pClu1.boost(bv);
	// pClu2.boost(bv);
	// std::cout << "Phi " << std::endl;
	// return;

	// doing the Boost into the Lab frame:
	boost::numeric::ublas::vector<double> Pclu1=boost::numeric::ublas::prod(Lambda,Pclu1Hat);
	boost::numeric::ublas::vector<double> Pclu2=boost::numeric::ublas::prod(Lambda,Pclu2Hat);

	// Computing the correct rotation, which maps pi/jRes into pi/jLab
	Axis zAxis(0,0,1);
	Axis xAxis(1,0,0);

	Lorentz5Momentum piRes(mi,Pi*zAxis);
	Lorentz5Momentum pjRes(mj,Pj(xAxissinPhi+zAxis*cosPhi));

	Lorentz5Momentum pClu1Res(M1,GeV*Axis(Pclu1[1],double(pClu1.y()/GeV),Pclu1[2]));
	Lorentz5Momentum pClu2Res(M2,GeV*Axis(Pclu2[1],double(pClu2.y()/GeV),Pclu2[2]));
	Axis omega1=piRes.vect().unit().cross(piLab.vect().unit());
	double cosAngle1=piRes.vect().unit()*piLab.vect().unit();
	double angle1=acos(cosAngle1);

	// Rotate piRes into piLab
	piRes.rotate(angle1, omega1);
	pjRes.rotate(angle1, omega1);
	// Correspondingly do the actual rotation on pClu1Res and pClu2Res
	pClu1Res.rotate(angle1, omega1);
	pClu2Res.rotate(angle1, omega1);

	Axis omega2=piRes.vect().unit();
	Momentum3 r1dim=(pjLab.vect()-piRes.vect().unit()(pjLab.vect()piRes.vect().unit()));
	Momentum3 r2dim=(pjRes.vect()-piRes.vect().unit()(pjRes.vect()piRes.vect().unit()));

	if (r1dim.mag()==ZERO \|\| r2dim.mag()==ZERO \|\| fabs(sinPhi)<1e-14) //trivial rotation so we are done
	{
	pClu1=pClu1Res;
	pClu2=pClu2Res;
	return;
	}
	Axis r1=r1dim.unit();
	Axis r2=r2dim.unit();

	// signs for 2nd rotation
	int signToPi = (piRes.vect()*pjLab.vect())/GeV2 > 0 ? 1:-1;
	int signToR1R2 = signToPi(r2.unit().cross(r1.unit())piRes.vect())/GeV> 0 ? 1:-1;
	double angle2=acos(r1.unit()*r2.unit());
	if (signToR1R2<0) angle2=-angle2;

	// Rotate pjRes into pjLab
	pjRes.rotate(angle2, signToPi*omega2);
	// Correspondingly do the actual rotation on pClu1Res and pClu2Res
	pClu1Res.rotate(angle2, signToPi*omega2);
	pClu2Res.rotate(angle2, signToPi*omega2);
	pClu1=pClu1Res;
	pClu2=pClu2Res;
	}
	-
	+*/
	bool Kinematics::twoBodyDecayNoBoost(const Lorentz5Momentum & p,
	const Energy m1, const Energy m2,
	const double cosThetaK, const double phiK,
	Lorentz5Momentum & p1, Lorentz5Momentum & p2) {
	Energy M = p.mass();
	Energy2 M2 = M*M;
	Energy2 m12 = m1*m1;
	Energy2 m22 = m2*m2;
	double cph = cos(phiK);
	double sph = sin(phiK);
	Energy Q = p.vect().mag();
	Energy EQ = sqrt(M2+Q*Q);
	double B = sqr(M/Q);
	double sinThetaK = sqrt(1.-cosThetaK*cosThetaK);
	Energy Pcom = sqrt(Kinematics::kaellenV(M,m1,m2))/(2.0*M);
	double eta1 = (M2+m12-m22)/(2.0*M2);
	double eta2 = (M2-m12+m22)/(2.0*M2);
	Energy k = EQPcom/(Qsqrt(1.0+B-cosThetaK*cosThetaK));
	Energy k0 = -EQ*eta1;
	Energy E1 = sqrt(kk+m12+2.0cosThetaKQketa1+sqr(Qeta1));
	k0 += E1;
	Momentum3 kvec3(kcphsinThetaK,ksphsinThetaK,k*cosThetaK);
	Lorentz5Momentum kvec5(kvec3,k0);

	Axis xAx(1.0,0.0,0.0);
	Axis yAx(0.0,1.0,0.0);
	Axis zAx(0.0,0.0,1.0);
	Lorentz5Momentum pz (Q*zAx,EQ);
	p1 = pz*eta1+kvec5;
	p2 = pz*eta2-kvec5;
	// double Norm = sqrt(B)*(1.0+B)/2.0;
	// double PhaseSpaceVol = (1.0+B)Pcom/(8.0Q*Norm);
	double theta = acos(zAx.cosTheta(p.vect()));
	double phi = theta==0 ? 0.0:atan2(p.vect().y()/GeV,p.vect().x()/GeV);
	p1.rotate(theta, yAx);
	p1.rotate(phi, zAx);
	p2.rotate(theta, yAx);
	p2.rotate(phi, zAx);
	if (fabs((p1.m()-m1)/GeV)>1e-5) std::cout << "p1.m = "<<p1.m()/GeV << " m1 = "<<m1/GeV << std::endl;
	if (fabs((p2.m()-m2)/GeV)>1e-5) std::cout << "p2.m = "<<p2.m()/GeV << " m2 = "<<m2/GeV << std::endl;
	if (fabs(((p1+p2).m()-M)/GeV)>1e-5) std::cout << "(p1+p2).m = "<<(p1+p2).m()/GeV << " M = "<<M/GeV << std::endl;
	return true;
	}


	bool Kinematics::twoBodyDecay(const Lorentz5Momentum & p,
	const Energy m1, const Energy m2,
	const Axis & unitDir1,
	Lorentz5Momentum & p1, Lorentz5Momentum & p2) {
	Energy min=p.mass();
	if ( min >= m1 + m2 && m1 >= ZERO && m2 >= ZERO ) {
	Momentum3 pstarVector = unitDir1 * Kinematics::pstarTwoBodyDecay(min,m1,m2);
	p1 = Lorentz5Momentum(m1, pstarVector);
	p2 = Lorentz5Momentum(m2,-pstarVector);
	// boost from CM to LAB
	Boost bv = p.boostVector();
	double gammarest = p.e()/p.mass();
	p1.boost( bv, gammarest );
	p2.boost( bv, gammarest );
	return true;
	}
	return false;
	}

	/*****
	* This function, as the name implies, performs a three body decay. The decay
	* products are distributed uniformly in all three directions.
	****/
	bool Kinematics::threeBodyDecay(Lorentz5Momentum p0, Lorentz5Momentum &p1,
	Lorentz5Momentum &p2, Lorentz5Momentum &p3,
	double (*fcn)(Energy2,Energy2,Energy2,InvEnergy4)) {
	// Variables needed in calculation...named same as fortran version
	Energy a = p0.mass() + p1.mass();
	Energy b = p0.mass() - p1.mass();
	Energy c = p2.mass() + p3.mass();

	if(b < c) {
	CurrentGenerator::log()
	<< "Kinematics::threeBodyDecay() phase space problem\n"
	<< p0.mass()/GeV << " -> "
	<< p1.mass()/GeV << ' '
	<< p2.mass()/GeV << ' '
	<< p3.mass()/GeV << '\n';
	return false;
	}

	Energy d = abs(p2.mass()-p3.mass());
	Energy2 aa = sqr(a);
	Energy2 bb = sqr(b);
	Energy2 cc = sqr(c);
	Energy2 dd = sqr(d);
	Energy2 ee = (b-c)*(a-d);

	Energy2 a1 = 0.5 * (aa+bb);
	Energy2 b1 = 0.5 * (cc+dd);
	InvEnergy4 c1 = 4./(sqr(a1-b1));

	Energy2 ff;
	double ww;
	Energy4 pp,qq,rr;
	// Choose mass of subsystem 23 with prescribed distribution
	const unsigned int MAXTRY = 100;
	unsigned int ntry=0;
	do {
	// ff is the mass squared of the 23 subsystem
	ff = UseRandom::rnd()*(cc-bb)+bb;

	// pp is ((m0+m1)^2 - m23^2)((m0-m1)^2-m23)
	pp = (aa-ff)*(bb-ff);

	// qq is ((m2+m3)^2 - m23^2)(\|m2-m3\|^2-m23^2)
	qq = (cc-ff)*(dd-ff);

	// weight
	ww = (fcn != NULL) ? (*fcn)(ff,a1,b1,c1) : 1.0;
	ww = sqr(ww);
	rr = eeffUseRandom::rnd();
	++ntry;
	}
	while(ppqqww < rr*rr && ntry < MAXTRY );
	if(ntry >= MAXTRY) {
	CurrentGenerator::log() << "Kinematics::threeBodyDecay can't generate momenta"
	<< " after " << MAXTRY << " attempts\n";
	return false;
	}

	// ff is the mass squared of subsystem 23
	// do 2 body decays 0->1+23, 23->2+3
	double CosAngle, AzmAngle;
	Lorentz5Momentum p23;

	p23.setMass(sqrt(ff));

	generateAngles(CosAngle,AzmAngle);
	bool status = twoBodyDecay(p0,p1.mass(),p23.mass(),CosAngle,AzmAngle,p1,p23);

	generateAngles(CosAngle,AzmAngle);
	status &= twoBodyDecay(p23,p2.mass(),p3.mass(),CosAngle,AzmAngle,p2,p3);
	return status;
	}
	diff --git a/Utilities/Kinematics.h b/Utilities/Kinematics.h
	--- a/Utilities/Kinematics.h
	+++ b/Utilities/Kinematics.h
	@@ -1,151 +1,152 @@
	// -- C++ --
	//
	// Kinematics.h 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.
	//
	//
	#ifndef HERWIG_Kinematics_H
	#define HERWIG_Kinematics_H

	// This is the declaration of the Kinematics class.

	#include <ThePEG/Config/ThePEG.h>
	#include "ThePEG/Vectors/ThreeVector.h"
	#include "ThePEG/Vectors/LorentzRotation.h"
	#include "ThePEG/Repository/UseRandom.h"
	#include <ThePEG/Vectors/Lorentz5Vector.h>

	namespace Herwig {

	using namespace ThePEG;

	/** \ingroup Utilities
	* This is a namespace which provides some useful methods
	* for kinematics computation, as the two body decays.
	*
	* NB) For other useful kinematical methods (and probably even those
	* implemented in Kinematics class!):
	* @see UtilityBase
	*/
	namespace Kinematics {


	/**
	* Calculate the momenta for a two body decay
	* The return value indicates success or failure.
	* @param p The momentum of the decaying particle
	* @param m1 The mass of the first decay product
	* @param m2 The mass of the second decay product
	* @param unitDir1 Direction for the products in the rest frame of
	* the decaying particle
	* @param p1 The momentum of the first decay product
	* @param p2 The momentum of the second decay product
	*/
	bool twoBodyDecay(const Lorentz5Momentum & p,
	const Energy m1, const Energy m2,
	const Axis & unitDir1,
	Lorentz5Momentum & p1, Lorentz5Momentum & p2);

	/**
	* Calculate the momenta for a two body decay without performing a boost
	* The return value indicates success or failure.
	* @param p The momentum of the decaying particle
	* @param m1 The mass of the first decay product
	* @param m2 The mass of the second decay product
	* @param unitDir1 Direction for the products in the rest frame of
	* the decaying particle
	* @param p1 The momentum of the first decay product
	* @param p2 The momentum of the second decay product
	*/
	bool twoBodyDecayNoBoost(const Lorentz5Momentum & p,
	const Energy m1, const Energy m2,
	const double cosThetaK, const double phiK,
	Lorentz5Momentum & p1, Lorentz5Momentum & p2);

	/*
	* Kaellen functions
	* */
	inline double kaellen(double x, double y, double z){
	return (xx-(y+z)(y+z))(xx-(y-z)*(y-z));
	}

	inline Energy4 kaellenV(Energy x, Energy y, Energy z){
	return (xx-(y+z)(y+z))(xx-(y-z)*(y-z));
	}

	/**
	* Boost consistently the Lorentz5Momenta pClu1 and pClu2 (given in their COM frame)
	* into the p0Q1 p0Q2 frame.
	* */
	+ /*
	void BoostIntoTwoParticleFrame(const Energy M, const Lorentz5Momentum & piLab,
	const Lorentz5Momentum & pjLab,
	Lorentz5Momentum & pClu1,
	- Lorentz5Momentum & pClu2);
	+ Lorentz5Momentum & pClu2);*/
	void BoostIntoTwoParticleFrame(const Energy M, const Lorentz5Momentum & piLab,
	const Lorentz5Momentum & pjLab,
	std::vector<Lorentz5Momentum * > momenta);
	/**
	* It returns the unit 3-vector with the given cosTheta and phi.
	*/
	inline Axis unitDirection(const double cosTheta, const double phi) {
	return ( fabs( cosTheta ) <= 1.0 ?
	Axis( cos(phi)sqrt(1.0-cosThetacosTheta) ,
	sin(phi)sqrt(1.0-cosThetacosTheta) , cosTheta) : Axis() );
	}

	/**
	* Calculate the momenta for a two body decay
	* The return value indicates success or failure.
	* @param p The momentum of the decaying particle
	* @param m1 The mass of the first decay product
	* @param m2 The mass of the second decay product
	* @param cosThetaStar1 Polar angle in rest frame
	* @param phiStar1 Azimuthal angle in rest frame
	* @param p1 The momentum of the first decay product
	* @param p2 The momentum of the second decay product
	*/
	inline bool twoBodyDecay(const Lorentz5Momentum & p,
	const Energy m1, const Energy m2,
	const double cosThetaStar1,
	const double phiStar1,
	Lorentz5Momentum & p1, Lorentz5Momentum & p2) {
	return twoBodyDecay(p,m1,m2,unitDirection(cosThetaStar1,phiStar1),p1,p2);
	}

	/**
	* As the name implies, this takes the momentum p0 and does a flat three
	* body decay into p1..p3. The argument fcn is used to add additional
	* weights. If it is not used, the default is just flat in phasespace.
	* The return value indicates success or failure.
	*/
	bool threeBodyDecay(Lorentz5Momentum p0, Lorentz5Momentum &p1,
	Lorentz5Momentum &p2, Lorentz5Momentum &p3,
	double (*fcn)(Energy2,Energy2,Energy2,InvEnergy4) = NULL);

	/**
	* For the two body decay M -> m1 + m2 it gives the module of the
	* 3-momentum of the decay product in the rest frame of M.
	*/
	inline Energy pstarTwoBodyDecay(const Energy M,
	const Energy m1, const Energy m2) {
	return ( M > ZERO && m1 >=ZERO && m2 >= ZERO && M > m1+m2 ?
	Energy(sqrt(( sqr(M) - sqr(m1+m2) )*( sqr(M) - sqr(m1-m2) ))
	/ (2.0*M) ) : ZERO);
	}

	/**
	* This just generates angles. First flat -1..1, second flat 0..2Pi
	*/
	inline void generateAngles(double & ct, double & az) {
	ct = UseRandom::rnd()*2.0 - 1.0; // Flat from -1..1
	az = UseRandom::rnd()*Constants::twopi;
	}
	}

	}

	#endif /* HERWIG_Kinematics_H */

File Metadata

Mime Type: text/x-diff
Expires: Sun, Feb 23, 2:01 PM (1 h, 51 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4479888
Default Alt Text: (55 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions