No OneTemporary
Actions

Size

167 KB

Subscribers

None

View Options

	diff --git a/Hadronization/ClusterFissioner.cc b/Hadronization/ClusterFissioner.cc
	--- a/Hadronization/ClusterFissioner.cc
	+++ b/Hadronization/ClusterFissioner.cc
	@@ -1,3616 +1,3683 @@
	// -- C++ --
	//
	// ClusterFissioner.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.
	//
	//
	// Thisk is the implementation of the non-inlined, non-templated member
	// functions of the ClusterFissioner class.
	//

	#include "ClusterFissioner.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 "Herwig/Utilities/Kinematics.h"
	#include "Cluster.h"
	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include <ThePEG/Utilities/DescribeClass.h>
	#include "ThePEG/Interface/ParMap.h"

	#include "Herwig/Utilities/AlphaS.h"

	#include <boost/numeric/ublas/matrix.hpp>
	#include <boost/numeric/ublas/io.hpp>
	#include <boost/numeric/ublas/lu.hpp>
	#include <cassert>
	#include <vector>
	using namespace Herwig;

	DescribeClass<ClusterFissioner,Interfaced>
	describeClusterFissioner("Herwig::ClusterFissioner","Herwig.so");

	// Initialize counters
	unsigned int ClusterFissioner::_counterMaxLoopViolations=0;
	unsigned int ClusterFissioner::_counterPaccGreater1=0;
	unsigned int ClusterFissioner::_counterFissionMatrixElement=0;
	ClusterFissioner::ClusterFissioner() :
	_clMaxLight(3.35*GeV),
	_clMaxDiquark(3.35*GeV),
	_clMaxExotic(3.35*GeV),
	_clPowLight(2.0),
	_clPowDiquark(2.0),
	_clPowExotic(2.0),
	_pSplitLight(1.0),
	_pSplitExotic(1.0),
	_phaseSpaceWeights(0),
	_dim(4),
	_fissionCluster(0),
	_kinematicThresholdChoice(0),
	_pwtDIquark(0.0),
	_diquarkClusterFission(0),
	_btClM(1.0*GeV),
	_iopRem(1),
	_kappa(1.0e15*GeV/meter),
	_enhanceSProb(0),
	_m0Fission(2.*GeV),
	_massMeasure(0),
	_probPowFactor(4.0),
	_probShift(0.0),
	_kinThresholdShift(1.0*sqr(GeV)),
	_strictDiquarkKinematics(0),
	_covariantBoost(false),
	_allowHadronFinalStates(2),
	_massSampler(0),
	_phaseSpaceSamplerCluster(0),
	_phaseSpaceSamplerConstituent(0),
	_matrixElement(0),
	_fissionApproach(0),
	- _powerLawPower(-2.0),
	+ _powerLawPower(0.0),
	_maxLoopFissionMatrixElement(5000000),
	_safetyFactorMatrixElement(10.0),
	_epsilonIR(1.0),
	_failModeMaxLoopMatrixElement(0),
	_writeOut(0),
	_hadronizingStrangeDiquarks(2)
	{
	}
	ClusterFissioner::~ClusterFissioner(){
	if (_counterMaxLoopViolations){
	std::cout << "WARNING: There have been "
	<< _counterMaxLoopViolations
	<< "/" << _counterFissionMatrixElement
	<< " Maximum tries violations during the Simulation! ("
	<< 100.0*double(_counterMaxLoopViolations)/double(_counterFissionMatrixElement)
	<< " %)\n"
	<< "You may want to reduce ClusterFissioner:SafetyFactorOverEst (=>potentially Pacc>=1) "
	<< "and/or increase ClusterFissioner:MaxLoopMatrixElement (=>slower performance)\n";
	}
	if (_counterPaccGreater1){
	std::cout << "WARNING: There have been "
	<< _counterPaccGreater1
	<< "/" << _counterFissionMatrixElement
	<< " Pacc>=1 exceptions during the Simulation! ("
	<< 100.0*double(_counterPaccGreater1)/double(_counterFissionMatrixElement)
	<< " %)\n"
	<< "You may want to increase ClusterFissioner:SafetyFactorOverEst to reduce this "
	<< "\nNOTE: look at the log file and look for the larges Pacc accepted and multiply "
	<< "\n ClusterFissioner:SafetyFactorOverEst by this factor";
	}
	}
	IBPtr ClusterFissioner::clone() const {
	return new_ptr(*this);
	}

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

	void ClusterFissioner::persistentOutput(PersistentOStream & os) const {
	os << ounit(_clMaxLight,GeV) << ounit(_clMaxHeavy,GeV) << ounit(_clMaxDiquark,GeV) << ounit(_clMaxExotic,GeV)
	<< _clPowLight << _clPowHeavy << _clPowDiquark << _clPowExotic
	<< _pSplitLight << _pSplitHeavy << _pSplitExotic
	<< _fissionCluster << _fissionPwt
	<< _pwtDIquark
	<< _diquarkClusterFission
	<< ounit(_btClM,GeV)
	<< _iopRem << ounit(_kappa, GeV/meter)
	<< _enhanceSProb << ounit(_m0Fission,GeV) << _massMeasure
	<< _dim << _phaseSpaceWeights
	<< _hadronSpectrum << _kinematicThresholdChoice
	<< _probPowFactor << _probShift << ounit(_kinThresholdShift,sqr(GeV))
	<< _strictDiquarkKinematics
	<< _covariantBoost
	<< _allowHadronFinalStates
	<< _massSampler
	<< _phaseSpaceSamplerCluster
	<< _phaseSpaceSamplerConstituent
	<< _matrixElement
	<< _fissionApproach
	<< _powerLawPower
	<< _maxLoopFissionMatrixElement
	<< _safetyFactorMatrixElement
	<< _epsilonIR
	<< _failModeMaxLoopMatrixElement
	<< _writeOut
	<< _hadronizingStrangeDiquarks
	;
	}

	void ClusterFissioner::persistentInput(PersistentIStream & is, int) {
	is >> iunit(_clMaxLight,GeV) >> iunit(_clMaxHeavy,GeV) >> iunit(_clMaxDiquark,GeV) >> iunit(_clMaxExotic,GeV)
	>> _clPowLight >> _clPowHeavy >> _clPowDiquark >> _clPowExotic
	>> _pSplitLight >> _pSplitHeavy >> _pSplitExotic
	>> _fissionCluster >> _fissionPwt
	>> _pwtDIquark
	>> _diquarkClusterFission
	>> iunit(_btClM,GeV)
	>> _iopRem >> iunit(_kappa, GeV/meter)
	>> _enhanceSProb >> iunit(_m0Fission,GeV) >> _massMeasure
	>> _dim >> _phaseSpaceWeights
	>> _hadronSpectrum >> _kinematicThresholdChoice
	>> _probPowFactor >> _probShift >> iunit(_kinThresholdShift,sqr(GeV))
	>> _strictDiquarkKinematics
	>> _covariantBoost
	>> _allowHadronFinalStates
	>> _massSampler
	>> _phaseSpaceSamplerCluster
	>> _phaseSpaceSamplerConstituent
	>> _matrixElement
	>> _fissionApproach
	>> _powerLawPower
	>> _maxLoopFissionMatrixElement
	>> _safetyFactorMatrixElement
	>> _epsilonIR
	>> _failModeMaxLoopMatrixElement
	>> _writeOut
	>> _hadronizingStrangeDiquarks
	;
	}
	/*
	namespace{
	void printV(Lorentz5Momentum p) {
	std::cout << "("<<p.e()/GeV<<"\|"<<p.vect().x()/GeV<<","<<p.vect().y()/GeV<<","<<p.vect().z()/GeV<<") Mass = "<<p.mass()/GeV<<" m = "<<p.m()/GeV<<"\n";
	}
	}
	*/
	void ClusterFissioner::doinit() {
	Interfaced::doinit();
	if (_writeOut){
	std::ofstream out("data_CluFis.dat", std::ios::out);
	out.close();
	}
	// Energy m1=1.2*GeV;
	// Energy m2=2.3*GeV;
	// std::cout <<std::setprecision(18)<< sqr(m1+m2)/GeV2 << std::endl;
	// std::cout <<std::setprecision(18)<< ((m1+m2)*(m1+m2))/GeV2 << std::endl;
	// std::cout <<std::setprecision(18)<< ((sqr(m1+m2))-((m1+m2)*(m1+m2)))/GeV2 << std::endl;
	/* TESTs for Two particle boost
	* TODO
	* */
	/*
	Energy M,M1,M2;
	Energy m,m1,m2;

	M=20.1*GeV;
	M1=10.1*GeV;
	M2=5.1*GeV;
	m2=0.325*GeV;
	m1=0.625*GeV;
	m=1.1*GeV;

	// double cosTheta1=0.3;
	// double phi1=0.7*3.14;
	// double cosTheta2=-0.4;
	// double phi2=-0.2*3.14;

	// Axis LABz(0,0,1);
	// Energy Plab=80.1*GeV;

	// Lorentz5Momentum P(M,Plab*LABz);
	// Energy PCOM=P.boost(-P.boostVector()).vect().mag();


	double cosTheta=0.0;
	double phi=0.7*3.14;

	Energy P1=1.2*GeV;
	Energy P2=10.2*GeV;
	Energy Q=1.2*GeV;
	Lorentz5Momentum p1(m1,P1Axis(cosThetacos(phi),cosTheta*sin(phi),sin(acos(cosTheta))));
	Lorentz5Momentum p2(m1*2,Momentum3(-p1.vect().x(),-p1.vect().y(),P2));
	Lorentz5Momentum qOld(m,Q*sampleDirectionIsotropic());
	Lorentz5Momentum qNew=qOld;
	Lorentz5Momentum pC=p1+p2;
	printV(p1);
	printV(p2);
	printV(p1+p2);
	printV(pC);
	printV(qOld);
	printV(qNew);
	// p1.boost(-pC.boostVector());
	// p2.boost(-pC.boostVector());
	// q.boost(-pC.boostVector());
	// printV(q);
	std::vector<Lorentz5Momentum *> v;
	v.push_back(&qNew);
	Kinematics::BoostIntoTwoParticleFrame(pC.m(),p1,p2,v);
	// Kinematics::twoBodDecay(pC,);
	printV(qOld);
	printV(qNew);
	// std::cout << p1 << std::endl;
	// std::cout << p2 << std::endl;
	// std::cout << p1+p2 << std::endl;
	// std::cout << Q << std::endl;




	*/








	// TODO: Some User warnings/errors but not complete list
	if (_matrixElement!=0 && _fissionApproach==0) generator()->logWarning( Exception(
	"For non-trivial MatrixElement you need to enable FissionApproach=New or Hybrid\n",
	Exception::warning));
	for ( const long& id : spectrum()->heavyHadronizingQuarks() ) {
	if ( _pSplitHeavy.find(id) == _pSplitHeavy.end() \|\|
	_clPowHeavy.find(id) == _clPowHeavy.end() \|\|
	_clMaxHeavy.find(id) == _clMaxHeavy.end() )
	throw InitException() << "not all parameters have been set for heavy quark cluster fission";
	}
	// for ( const long& id : spectrum()->lightHadronizingDiquarks() ) {
	// if (
	// _clPowDiquark.find(id) == _clPowDiquark.end() \|\|
	// _clMaxDiquark.find(id) == _clMaxDiquark.end() )
	// throw InitException() << "not all parameters have been set for diquarks quark cluster fission";
	// }
	// for default Pwts not needed to initialize
	if (_fissionCluster==0) return;
	for ( const long& id : spectrum()->lightHadronizingQuarks() ) {
	if ( _fissionPwt.find(id) == _fissionPwt.end() )
	// check that all relevant weights are set
	throw InitException() << "fission weights for light quarks have not been set";
	}
	/*
	// Not needed since we set Diquark weights from quark weights
	for ( const long& id : spectrum()->lightHadronizingDiquarks() ) {
	if ( _fissionPwt.find(id) == _fissionPwt.end() )
	throw InitException() << "fission weights for light diquarks have not been set";
	}*/
	double pwtDquark=_fissionPwt.find(ParticleID::d)->second;
	double pwtUquark=_fissionPwt.find(ParticleID::u)->second;
	double pwtSquark=_fissionPwt.find(ParticleID::s)->second;
	// ERROR: TODO makeDiquarkID is protected function ?
	// _fissionPwt[spectrum()->makeDiquarkID(ParticleID::d,ParticleID::d,3)] = _pwtDIquark * pwtDquark * pwtDquark;
	// _fissionPwt[spectrum()->makeDiquarkID(ParticleID::u,ParticleID::d,1)] = 0.5 * _pwtDIquark * pwtUquark * pwtDquark;
	// _fissionPwt[spectrum()->makeDiquarkID(ParticleID::u,ParticleID::u,3)] = _pwtDIquark * pwtUquark * pwtUquark;
	// _fissionPwt[spectrum()->makeDiquarkID(ParticleID::s,ParticleID::d,1)] = 0.5 * _pwtDIquark * pwtSquark * pwtDquark;
	// _fissionPwt[spectrum()->makeDiquarkID(ParticleID::s,ParticleID::u,1)] = 0.5 * _pwtDIquark * pwtSquark * pwtUquark;
	// _fissionPwt[spectrum()->makeDiquarkID(ParticleID::s,ParticleID::s,3)] = _pwtDIquark * pwtSquark * pwtSquark;
	// TODO better solution for this magic number alternative
	_fissionPwt[1103] = _pwtDIquark * pwtDquark * pwtDquark;
	_fissionPwt[2101] = 0.5 * _pwtDIquark * pwtUquark * pwtDquark;
	_fissionPwt[2203] = _pwtDIquark * pwtUquark * pwtUquark;
	if (_hadronizingStrangeDiquarks>0) {
	_fissionPwt[3101] = 0.5 * _pwtDIquark * pwtSquark * pwtDquark;
	_fissionPwt[3201] = 0.5 * _pwtDIquark * pwtSquark * pwtUquark;
	if (_hadronizingStrangeDiquarks==2) {
	_fissionPwt[3303] = _pwtDIquark* pwtSquark * pwtSquark;
	}
	}
	}

	void ClusterFissioner::Init() {

	static ClassDocumentation<ClusterFissioner> documentation
	("Class responsibles for chopping up the clusters");

	static Reference<ClusterFissioner,HadronSpectrum> interfaceHadronSpectrum
	("HadronSpectrum",
	"Set the Hadron spectrum for this cluster fissioner.",
	&ClusterFissioner::_hadronSpectrum, false, false, true, false);

	// ClMax for light, Bottom, Charm and exotic (e.g. Susy) quarks
	static Parameter<ClusterFissioner,Energy>
	interfaceClMaxLight ("ClMaxLight","cluster max mass for light quarks (unit [GeV])",
	&ClusterFissioner::_clMaxLight, GeV, 3.35GeV, ZERO, 100.0GeV,
	false,false,false);
	static Parameter<ClusterFissioner,Energy>
	interfaceClMaxDiquark ("ClMaxDiquark","cluster max mass for light hadronizing diquarks (unit [GeV])",
	&ClusterFissioner::_clMaxDiquark, GeV, 3.35GeV, ZERO, 100.0GeV,
	false,false,false);
	static ParMap<ClusterFissioner,Energy> interfaceClMaxHeavy
	("ClMaxHeavy",
	"ClMax for heavy quarks",
	&ClusterFissioner::_clMaxHeavy, GeV, -1, 3.35GeV, ZERO, 100.0GeV,
	false, false, Interface::upperlim);
	// static ParMap<ClusterFissioner,Energy> interfaceClMaxDiquark
	// ("ClMaxDiquark",
	// "ClMax for light hadronizing diquarks",
	// &ClusterFissioner::_clMaxDiquark, GeV, -1, 3.35GeV, ZERO, 100.0GeV,
	// false, false, Interface::upperlim);
	static Parameter<ClusterFissioner,Energy>
	interfaceClMaxExotic ("ClMaxExotic","cluster max mass for exotic quarks (unit [GeV])",
	&ClusterFissioner::_clMaxExotic, GeV, 3.35GeV, ZERO, 100.0GeV,
	false,false,false);

	// ClPow for light, Bottom, Charm and exotic (e.g. Susy) quarks
	static Parameter<ClusterFissioner,double>
	interfaceClPowLight ("ClPowLight","cluster mass exponent for light quarks",
	&ClusterFissioner::_clPowLight, 0, 2.0, 0.0, 10.0,false,false,false);
	static ParMap<ClusterFissioner,double> interfaceClPowHeavy
	("ClPowHeavy",
	"ClPow for heavy quarks",
	&ClusterFissioner::_clPowHeavy, -1, 1.0, 0.0, 10.0,
	false, false, Interface::upperlim);
	// static ParMap<ClusterFissioner,double> interfaceClPowDiquark
	// ("ClPowDiquark",
	// "ClPow for light hadronizing diquarks",
	// &ClusterFissioner::_clPowDiquark, -1, 1.0, 0.0, 10.0,
	// false, false, Interface::upperlim);
	static Parameter<ClusterFissioner,double>
	interfaceClPowDiquark ("ClPowDiquark","cluster mass exponent for light hadronizing diquarks",
	&ClusterFissioner::_clPowDiquark, 0, 2.0, 0.0, 10.0,false,false,false);
	static Parameter<ClusterFissioner,double>
	interfaceClPowExotic ("ClPowExotic","cluster mass exponent for exotic quarks",
	&ClusterFissioner::_clPowExotic, 0, 2.0, 0.0, 10.0,false,false,false);

	// PSplit for light, Bottom, Charm and exotic (e.g. Susy) quarks
	static Parameter<ClusterFissioner,double>
	interfacePSplitLight ("PSplitLight","cluster mass splitting param for light quarks",
	&ClusterFissioner::_pSplitLight, 0, 1.0, 0.0, 10.0,false,false,false);
	static ParMap<ClusterFissioner,double> interfacePSplitHeavy
	("PSplitHeavy",
	"PSplit for heavy quarks",
	&ClusterFissioner::_pSplitHeavy, -1, 1.0, 0.0, 10.0,
	false, false, Interface::upperlim);
	static Parameter<ClusterFissioner,double>
	interfacePSplitExotic ("PSplitExotic","cluster mass splitting param for exotic quarks",
	&ClusterFissioner::_pSplitExotic, 0, 1.0, 0.0, 10.0,false,false,false);


	static Switch<ClusterFissioner,int> interfaceFission
	("Fission",
	"Option for different Fission options",
	&ClusterFissioner::_fissionCluster, 1, false, false);
	static SwitchOption interfaceFissionDefault
	(interfaceFission,
	"Default",
	"Normal cluster fission which depends on the hadron spectrum class.",
	0);
	static SwitchOption interfaceFissionNew
	(interfaceFission,
	"New",
	"Alternative cluster fission which does not depend on the hadron spectrum class",
	1);
	static SwitchOption interfaceFissionNewDiquarkSuppression
	(interfaceFission,
	"NewDiquarkSuppression",
	"Alternative cluster fission which does not depend on the hadron spectrum class"
	" and includes a suppression of AlphaS^2(Mc) for Diquark Production during "
	"Cluster Fission",
	-1);


	static Switch<ClusterFissioner,int> interfaceFissionApproach
	("FissionApproach",
	"Option for different Cluster Fission approaches",
	&ClusterFissioner::_fissionApproach, 0, false, false);
	static SwitchOption interfaceFissionApproachDefault
	(interfaceFissionApproach,
	"Default",
	"Default Herwig-7.3.0 cluster fission without restructuring",
	0);
	static SwitchOption interfaceFissionApproachNew
	(interfaceFissionApproach,
	"New",
	"New cluster fission which allows to choose MassSampler"
	", PhaseSpaceSampler and MatrixElement",
	1);
	static SwitchOption interfaceFissionApproachHybrid
	(interfaceFissionApproach,
	"Hybrid",
	"New cluster fission which allows to choose MassSampler"
	", PhaseSpaceSampler and MatrixElement, but uses Default"
	" Approach for BeamClusters",
	2);
	static SwitchOption interfaceFissionApproachPreservePert
	(interfaceFissionApproach,
	"PreservePert",
	"New cluster fission which allows to choose MassSampler"
	", PhaseSpaceSampler and MatrixElement, but uses Default"
	" Approach for BeamClusters",
	3);
	static SwitchOption interfaceFissionApproachPreserveNonPert
	(interfaceFissionApproach,
	"PreserveNonPert",
	"New cluster fission which allows to choose MassSampler"
	", PhaseSpaceSampler and MatrixElement, but uses Default"
	" Approach for BeamClusters",
	4);
	static SwitchOption interfaceFissionApproachPreserveFirstPert
	(interfaceFissionApproach,
	"PreserveFirstPert",
	"New cluster fission which allows to choose MassSampler"
	", PhaseSpaceSampler and MatrixElement, but uses Default"
	" Approach for BeamClusters",
	5);

	// Switch C->H1,C2 C->H1,H2 on or off
	static Switch<ClusterFissioner,int> interfaceAllowHadronFinalStates
	("AllowHadronFinalStates",
	"Option for allowing hadron final states of cluster fission",
	&ClusterFissioner::_allowHadronFinalStates, 0, false, false);
	static SwitchOption interfaceAllowHadronFinalStatesNone
	(interfaceAllowHadronFinalStates,
	"None",
	"Option for disabling hadron final states of cluster fission",
	0);
	static SwitchOption interfaceAllowHadronFinalStatesSemiHadronicOnly
	(interfaceAllowHadronFinalStates,
	"SemiHadronicOnly",
	"Option for allowing hadron final states of cluster fission of type C->H1,C2 "
	"(NOT YET USABLE WITH MatrixElement only use option None)",
	1);
	static SwitchOption interfaceAllowHadronFinalStatesAll
	(interfaceAllowHadronFinalStates,
	"All",
	"Option for allowing hadron final states of cluster fission "
	"of type C->H1,C2 or C->H1,H2 "
	"(NOT YET USABLE WITH MatrixElement only use option None)",
	2);

	// Mass Sampler Switch
	static Switch<ClusterFissioner,int> interfaceMassSampler
	("MassSampler",
	"Option for different Mass sampling options",
	&ClusterFissioner::_massSampler, 0, false, false);
	static SwitchOption interfaceMassSamplerDefault
	(interfaceMassSampler,
	"Default",
	"Choose H7.2.3 default mass sampling using PSplitX",
	0);
	static SwitchOption interfaceMassSamplerUniform
	(interfaceMassSampler,
	"Uniform",
	"Choose Uniform Mass sampling in M1,M2 space",
	1);
	static SwitchOption interfaceMassSamplerFlatPhaseSpace
	(interfaceMassSampler,
	"FlatPhaseSpace",
	"Choose Flat Phase Space sampling of Mass in M1,M2 space (Recommended)",
	2);
	+ static SwitchOption interfaceMassSamplerPhaseSpacePowerLaw
	+ (interfaceMassSampler,
	+ "PhaseSpacePowerLaw",
	+ "Choose Phase Space times a power law sampling of Mass in M1,M2 space.",
	+ 3);
	+
	+ static Parameter<ClusterFissioner,double>
	+ interfaceMassPowerLaw("MassPowerLaw",
	+ "power of mass power law modulation of FlatPhaseSpace",
	+ &ClusterFissioner::_powerLawPower, 0.0, -10.0, 0.0, 10.0,
	+ false,false,false);

	// Phase Space Sampler Switch
	static Switch<ClusterFissioner,int> interfacePhaseSpaceSamplerCluster
	("PhaseSpaceSamplerCluster",
	"Option for different phase space sampling options",
	&ClusterFissioner::_phaseSpaceSamplerCluster, 0, false, false);
	static SwitchOption interfacePhaseSpaceSamplerClusterAligned
	(interfacePhaseSpaceSamplerCluster,
	"Aligned",
	"Draw the momentum of child clusters to be aligned"
	" the relative momentum of the mother cluster",
	0);
	static SwitchOption interfacePhaseSpaceSamplerClusterIsotropic
	(interfacePhaseSpaceSamplerCluster,
	"Isotropic",
	"Draw the momentum of child clusters to be isotropic"
	" the relative momentum of the mother cluster",
	1);
	static SwitchOption interfacePhaseSpaceSamplerClusterTchannel
	(interfacePhaseSpaceSamplerCluster,
	"Tchannel",
	"Draw the momentum of child clusters to be t-channel like"
	" the relative momentum of the mother cluster"
	" i.e. cosTheta ~ 1/(1+A-cosTheta)^2",
	2);

	static Switch<ClusterFissioner,int> interfacePhaseSpaceSamplerConstituents
	("PhaseSpaceSamplerConstituents",
	"Option for different phase space sampling options",
	&ClusterFissioner::_phaseSpaceSamplerConstituent, 0, false, false);
	static SwitchOption interfacePhaseSpaceSamplerConstituentsAligned
	(interfacePhaseSpaceSamplerConstituents,
	"Aligned",
	"Draw the momentum of the constituents of the child clusters "
	"to be aligned relative to the direction of the child cluster",
	0);
	static SwitchOption interfacePhaseSpaceSamplerConstituentsIsotropic
	(interfacePhaseSpaceSamplerConstituents,
	"Isotropic",
	"Draw the momentum of the constituents of the child clusters "
	"to be isotropic relative to the direction of the child cluster",
	1);
	static SwitchOption interfacePhaseSpaceSamplerConstituentsTchannel
	(interfacePhaseSpaceSamplerConstituents,
	"Exponential",
	"Draw the momentum of the constituents of the child clusters "
	"to be exponential relative to the direction of the child cluster"
	" i.e. cosTheta ~ exp(lam*(cosTheta-1)",
	2);



	// Matrix Element Choice Switch
	static Switch<ClusterFissioner,int> interfaceMatrixElement
	("MatrixElement",
	"Option for different Matrix Element options",
	&ClusterFissioner::_matrixElement, 0, false, false);
	static SwitchOption interfaceMatrixElementDefault
	(interfaceMatrixElement,
	"Default",
	"Choose trivial matrix element i.e. whatever comes from the mass and "
	"phase space sampling",
	0);
	static SwitchOption interfaceMatrixElementSoftQQbarFinalFinal
	(interfaceMatrixElement,
	"SoftQQbarFinalFinal",
	"Choose Soft p1,p2->q1,q2,g*->q1,q2,q,qbar matrix element"
	"NOTE: Here the matrix element depends on qi.q(bar)",
	1);
	static SwitchOption interfaceMatrixElementSoftQQbarInitialFinal
	(interfaceMatrixElement,
	"SoftQQbarInitialFinal",
	"Choose Soft p1,p2->q1,q2,g*->q1,q2,q,qbar matrix element "
	"NOTE: Here the matrix element depends on pi.q(bar)",
	2);
	static SwitchOption interfaceMatrixElementTest
	(interfaceMatrixElement,
	"Test",
	"Choose Soft p1,p2->q1,q2,g*->q1,q2,q,qbar matrix element "
	"NOTE: Here the matrix element depends on pi.q(bar)",
	4);
	static SwitchOption interfaceMatrixElementSoftQQbarInitialFinalTchannel
	(interfaceMatrixElement,
	"SoftQQbarInitialFinalTchannel",
	"Choose Soft p1,p2->q1,q2,g*->q1,q2,q,qbar matrix element "
	"NOTE: Here the matrix element depends on pi.q(bar)",
	5);
	// Technical Max loop parameter for New ClusterFission approach and Matrix Element
	// Rejection sampling
	static Parameter<ClusterFissioner,int> interfaceMaxLoopMatrixElement
	("MaxLoopMatrixElement",
	"Technical Parameter for how many tries are allowed to sample the "
	"Cluster Fission matrix element before reverting to fissioning "
	"using the default Fission Aproach",
	&ClusterFissioner::_maxLoopFissionMatrixElement, 5000000, 100, 1e8,
	false, false, Interface::limited);
	static Switch<ClusterFissioner,int> interfaceFailModeMaxLoopMatrixElement
	("FailModeMaxLoopMatrixElement",
	"How to fail if we try more often than MaxLoopMatrixElement",
	&ClusterFissioner::_failModeMaxLoopMatrixElement, 0, false, false);
	static SwitchOption interfaceDiquarkClusterFissionOldFission
	(interfaceFailModeMaxLoopMatrixElement,
	"OldFission",
	"Use old Cluster Fission Aproach if we try more often than MaxLoopMatrixElement",
	0);
	static SwitchOption interfaceDiquarkClusterFissionRejectEvent
	(interfaceFailModeMaxLoopMatrixElement,
	"RejectEvent",
	"Reject Event if we try more often than MaxLoopMatrixElement",
	1);

	static Switch<ClusterFissioner,int> interfaceDiquarkClusterFission
	("DiquarkClusterFission",
	"Allow clusters to fission to 1 or 2 diquark Clusters or Turn off diquark fission completely",
	&ClusterFissioner::_diquarkClusterFission, 0, false, false);
	static SwitchOption interfaceDiquarkClusterFissionAll
	(interfaceDiquarkClusterFission,
	"All",
	"Allow diquark clusters and baryon clusters to fission to new diquark Clusters",
	2);
	static SwitchOption interfaceDiquarkClusterFissionOnlyBaryonClusters
	(interfaceDiquarkClusterFission,
	"OnlyBaryonClusters",
	"Allow only baryon clusters to fission to new diquark Clusters",
	1);
	static SwitchOption interfaceDiquarkClusterFissionNo
	(interfaceDiquarkClusterFission,
	"No",
	"Don't allow clusters to fission to new diquark Clusters",
	0);
	static SwitchOption interfaceDiquarkClusterFissionOff
	(interfaceDiquarkClusterFission,
	"Off",
	"Don't allow clusters fission to draw diquarks ",
	-1);

	static ParMap<ClusterFissioner,double> interfaceFissionPwt
	("FissionPwt",
	"The weights for quarks in the fission process.",
	&ClusterFissioner::_fissionPwt, -1, 1.0, 0.0, 10.0,
	false, false, Interface::upperlim);

	static Switch<ClusterFissioner,int> interfaceRemnantOption
	("RemnantOption",
	"Option for the treatment of remnant clusters",
	&ClusterFissioner::_iopRem, 1, false, false);
	static SwitchOption interfaceRemnantOptionSoft
	(interfaceRemnantOption,
	"Soft",
	"Both clusters produced in the fission of the beam cluster"
	" are treated as soft clusters.",
	0);
	static SwitchOption interfaceRemnantOptionHard
	(interfaceRemnantOption,
	"Hard",
	"Only the cluster containing the remnant is treated as a soft cluster.",
	1);
	static SwitchOption interfaceRemnantOptionVeryHard
	(interfaceRemnantOption,
	"VeryHard",
	"Even remnant clusters are treated as hard, i.e. all clusters the same",
	2);

	static Parameter<ClusterFissioner,Energy> interfaceBTCLM
	("SoftClusterFactor",
	"Parameter for the mass spectrum of remnant clusters",
	&ClusterFissioner::_btClM, GeV, 1.GeV, 0.1GeV, 10.0*GeV,
	false, false, Interface::limited);


	static Parameter<ClusterFissioner,Tension> interfaceStringTension
	("StringTension",
	"String tension used in vertex displacement calculation",
	&ClusterFissioner::_kappa, GeV/meter,
	1.0e15*GeV/meter, ZERO, ZERO,
	false, false, Interface::lowerlim);

	static Switch<ClusterFissioner,int> interfaceEnhanceSProb
	("EnhanceSProb",
	"Option for enhancing strangeness",
	&ClusterFissioner::_enhanceSProb, 0, false, false);
	static SwitchOption interfaceEnhanceSProbNo
	(interfaceEnhanceSProb,
	"No",
	"No strangeness enhancement.",
	0);
	static SwitchOption interfaceEnhanceSProbScaled
	(interfaceEnhanceSProb,
	"Scaled",
	"Scaled strangeness enhancement",
	1);
	static SwitchOption interfaceEnhanceSProbExponential
	(interfaceEnhanceSProb,
	"Exponential",
	"Exponential strangeness enhancement",
	2);

	static Switch<ClusterFissioner,int> interfaceMassMeasure
	("MassMeasure",
	"Option to use different mass measures",
	&ClusterFissioner::_massMeasure,0,false,false);
	static SwitchOption interfaceMassMeasureMass
	(interfaceMassMeasure,
	"Mass",
	"Mass Measure",
	0);
	static SwitchOption interfaceMassMeasureLambda
	(interfaceMassMeasure,
	"Lambda",
	"Lambda Measure",
	1);

	static Parameter<ClusterFissioner,Energy> interfaceFissionMassScale
	("FissionMassScale",
	"Cluster fission mass scale",
	&ClusterFissioner::_m0Fission, GeV, 2.0GeV, 0.1GeV, 50.*GeV,
	false, false, Interface::limited);

	static Parameter<ClusterFissioner,double> interfaceProbPowFactor
	("ProbablityPowerFactor",
	"Power factor in ClausterFissioner bell probablity function",
	&ClusterFissioner::_probPowFactor, 2.0, 0.001, 20.0,
	false, false, Interface::limited);

	static Parameter<ClusterFissioner,double> interfaceProbShift
	("ProbablityShift",
	"Shifts from the center in ClausterFissioner bell probablity function",
	&ClusterFissioner::_probShift, 0.0, -10.0, 10.0,
	false, false, Interface::limited);

	static Parameter<ClusterFissioner,Energy2> interfaceKineticThresholdShift
	("KineticThresholdShift",
	"Shifts from the kinetic threshold in ClusterFissioner",
	&ClusterFissioner::_kinThresholdShift, sqr(GeV), 0.sqr(GeV), -10.0sqr(GeV), 10.0*sqr(GeV),
	false, false, Interface::limited);

	static Switch<ClusterFissioner,int> interfaceKinematicThreshold
	("KinematicThreshold",
	"Option for using static or dynamic kinematic thresholds in cluster splittings",
	&ClusterFissioner::_kinematicThresholdChoice, 0, false, false);
	static SwitchOption interfaceKinematicThresholdStatic
	(interfaceKinematicThreshold,
	"Static",
	"Set static kinematic thresholds for cluster splittings.",
	0);
	static SwitchOption interfaceKinematicThresholdDynamic
	(interfaceKinematicThreshold,
	"Dynamic",
	"Set dynamic kinematic thresholds for cluster splittings.",
	1);
	static SwitchOption interfaceKinematicThresholdProbabilistic
	(interfaceKinematicThreshold,
	"Probabilistic",
	"Set Probabilistic kinematic threshold for cluster fission.",
	2);



	static Switch<ClusterFissioner,bool> interfaceCovariantBoost
	("CovariantBoost",
	"Use single Covariant Boost for Cluster Fission",
	&ClusterFissioner::_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<ClusterFissioner,int> interfaceStrictDiquarkKinematics
	("StrictDiquarkKinematics",
	"Option for selecting different selection criterions of diquarks for ClusterFission",
	&ClusterFissioner::_strictDiquarkKinematics, 0, false, false);
	static SwitchOption interfaceStrictDiquarkKinematicsLoose
	(interfaceStrictDiquarkKinematics,
	"Loose",
	"No kinematic threshold for diquark selection except for Mass bigger than 2 baryons",
	0);
	static SwitchOption interfaceStrictDiquarkKinematicsStrict
	(interfaceStrictDiquarkKinematics,
	"Strict",
	"Resulting clusters are at least as heavy as 2 lightest baryons",
	1);

	static Parameter<ClusterFissioner,double> interfacePwtDIquark
	("PwtDIquark",
	"specific probability for choosing a d diquark",
	&ClusterFissioner::_pwtDIquark, 0.0, 0.0, 10.0,
	false, false, Interface::limited);

	static Switch<ClusterFissioner,int> interfacePhaseSpaceWeights
	("PhaseSpaceWeights",
	"Include phase space weights.",
	&ClusterFissioner::_phaseSpaceWeights, 0, false, false);
	static SwitchOption interfacePhaseSpaceWeightsNo
	(interfacePhaseSpaceWeights,
	"No",
	"Do not include the effect of cluster phase space",
	0);
	static SwitchOption interfacePhaseSpaceWeightsYes
	(interfacePhaseSpaceWeights,
	"Yes",
	"Do include the effect of cluster fission phase space "
	"related to constituent masses."
	"Note: Need static Threshold choice",
	1);
	static SwitchOption interfacePhaseSpaceWeightsUseHadronMasses
	(interfacePhaseSpaceWeights,
	"UseHadronMasses",
	"Do include the effect of cluster fission phase space "
	"related to hadron masses."
	"Note: Need static Threshold choice",
	2);
	static SwitchOption interfacePhaseSpaceWeightsNoConstituentMasses
	(interfacePhaseSpaceWeights,
	"NoConstituentMasses",
	"Do not include the effect of cluster fission phase space "
	"related to constituent masses."
	"Note: Need static Threshold choice",
	3);

	static Parameter<ClusterFissioner,double>
	interfaceDim ("Dimension","Dimension in which phase space weights are calculated",
	&ClusterFissioner::_dim, 0, 4.0, 0.0, 10.0,false,false,false);

	// Matrix Element Choice Switch
	static Parameter<ClusterFissioner,double>
	interfaceSafetyFactorOverEst("SafetyFactorOverEst","Safety factor with which the numerical overestimate is calculated",
	&ClusterFissioner::_safetyFactorMatrixElement, 0, 10.0, 1.0, 1000.0,false,false,false);
	// Matrix Element IR cutoff
	static Parameter<ClusterFissioner,double>
	interfaceMatrixElementIRCutOff("MatrixElementIRCutOff","IR CutOff for MatrixElement with IR divergences. "
	"for MatrixElement SoftQQbarInitialFinalTchannel e.g. such that 1/(t^2)->1/((t-_epsilonIR*mGluonConst^2)^2)",
	&ClusterFissioner::_epsilonIR, 0, 0.01, 1e-6, 1000.0,false,false,false);

	// Matrix Element Choice Switch
	static Switch<ClusterFissioner,int> interfaceWriteOut
	("WriteOut",
	"Option for different Matrix Element options",
	&ClusterFissioner::_writeOut, 0, false, false);
	static SwitchOption interfaceWriteOutNo
	(interfaceWriteOut,
	"No",
	"Choose trivial matrix element i.e. whatever comes from the mass and "
	"phase space sampling",
	0);
	static SwitchOption interfaceWriteOutYes
	(interfaceWriteOut,
	"Yes",
	"Choose Soft q1,q2->q1,q2,g*->q1,q2,q,qbar matrix element",
	1);

	// Allowing for strange diquarks in the ClusterFission
	static Switch<ClusterFissioner,unsigned int> interfaceHadronizingStrangeDiquarks
	("HadronizingStrangeDiquarks",
	"Option for adding strange diquarks to Cluster Fission (if Fission = New or Hybrid is enabled)",
	&ClusterFissioner::_hadronizingStrangeDiquarks, 0, false, false);
	static SwitchOption interfaceHadronizingStrangeDiquarksNo
	(interfaceHadronizingStrangeDiquarks,
	"No",
	"No strangeness containing diquarks during Cluster Fission",
	0);
	static SwitchOption interfaceHadronizingStrangeDiquarksOnlySingleStrange
	(interfaceHadronizingStrangeDiquarks,
	"OnlySingleStrange",
	"Only one strangeness containing diquarks during Cluster Fission i.e. su,sd",
	1);
	static SwitchOption interfaceHadronizingStrangeDiquarksAll
	(interfaceHadronizingStrangeDiquarks,
	"All",
	"All strangeness containing diquarks during Cluster Fission i.e. su,sd,ss",
	2);

	}

	tPVector ClusterFissioner::fission(ClusterVector & clusters, bool softUEisOn) {
	// return if no clusters
	if (clusters.empty()) return tPVector();

	/*****************
	* Loop over the (input) collection of cluster pointers, and store in
	* the vector splitClusters all the clusters that need to be split
	* (these are beam clusters, if soft underlying event is off, and
	* heavy non-beam clusters).
	********************/

	stack<ClusterPtr> splitClusters;
	for(ClusterVector::iterator it = clusters.begin() ;
	it != clusters.end() ; ++it) {
	/**************
	* Skip 3-component clusters that have been redefined (as 2-component
	* clusters) or not available clusters. The latter check is indeed
	* redundant now, but it is used for possible future extensions in which,
	* for some reasons, some of the clusters found by ClusterFinder are tagged
	* straight away as not available.
	**************/
	if((it)->isRedefined() \|\| !(it)->isAvailable()) continue;
	// if the cluster is a beam cluster add it to the vector of clusters
	// to be split or if it is heavy
	if((it)->isBeamCluster() \|\| isHeavy(it)) splitClusters.push(*it);
	}
	tPVector finalhadrons;
	cut(splitClusters, clusters, finalhadrons, softUEisOn);
	return finalhadrons;
	}

	void ClusterFissioner::cut(stack<ClusterPtr> & clusterStack,
	ClusterVector &clusters, tPVector & finalhadrons,
	bool softUEisOn) {
	/**************************************************
	* This method does the splitting of the cluster pointed by cluPtr
	* and "recursively" by all of its cluster children, if heavy. All of these
	* new children clusters are added (indeed the pointers to them) to the
	* collection of cluster pointers collecCluPtr. The method works as follows.
	* Initially the vector vecCluPtr contains just the input pointer to the
	* cluster to be split. Then it will be filled "recursively" by all
	* of the cluster's children that are heavy enough to require, in their turn,
	* to be split. In each loop, the last element of the vector vecCluPtr is
	* considered (only once because it is then removed from the vector).
	* This approach is conceptually recursive, but avoid the overhead of
	* a concrete recursive function. Furthermore it requires minimal changes
	* in the case that the fission of an heavy cluster could produce more
	* than two cluster children as assumed now.
	*
	* Draw the masses: for normal, non-beam clusters a power-like mass dist
	* is used, whereas for beam clusters a fast-decreasing exponential mass
	* dist is used instead (to avoid many iterative splitting which could
	* produce an unphysical large transverse energy from a supposed soft beam
	* remnant process).
	****************************************/
	// Here we recursively loop over clusters in the stack and cut them
	while (!clusterStack.empty()) {
	// take the last element of the vector
	ClusterPtr iCluster = clusterStack.top(); clusterStack.pop();
	// split it
	cutType ct = iCluster->numComponents() == 2 ?
	cutTwo(iCluster, finalhadrons, softUEisOn) :
	cutThree(iCluster, finalhadrons, softUEisOn);

	// There are cases when we don't want to split, even if it fails mass test
	if(!ct.first.first \|\| !ct.second.first) {
	// if an unsplit beam cluster leave if for the underlying event
	if(iCluster->isBeamCluster() && softUEisOn)
	iCluster->isAvailable(false);
	continue;
	}
	// check if clusters
	ClusterPtr one = dynamic_ptr_cast<ClusterPtr>(ct.first.first);
	ClusterPtr two = dynamic_ptr_cast<ClusterPtr>(ct.second.first);
	// is a beam cluster must be split into two clusters
	if(iCluster->isBeamCluster() && (!one\|\|!two) && softUEisOn) {
	iCluster->isAvailable(false);
	continue;
	}

	// There should always be a intermediate quark(s) from the splitting
	assert(ct.first.second && ct.second.second);
	/// \todo sort out motherless quark pairs here. Watch out for 'quark in final state' errors
	iCluster->addChild(ct.first.first);
	// iCluster->addChild(ct.first.second);
	// ct.first.second->addChild(ct.first.first);

	iCluster->addChild(ct.second.first);
	// iCluster->addChild(ct.second.second);
	// ct.second.second->addChild(ct.second.first);

	// Sometimes the clusters decay C -> H + C' or C -> H + H' rather then C -> C' + C''
	if(one) {
	clusters.push_back(one);
	if(one->isBeamCluster() && softUEisOn)
	one->isAvailable(false);
	if(isHeavy(one) && one->isAvailable())
	clusterStack.push(one);
	}
	if(two) {
	clusters.push_back(two);
	if(two->isBeamCluster() && softUEisOn)
	two->isAvailable(false);
	if(isHeavy(two) && two->isAvailable())
	clusterStack.push(two);
	}
	}
	}
	ClusterFissioner::cutType
	ClusterFissioner::cutTwo(ClusterPtr & cluster, tPVector & finalhadrons,
	bool softUEisOn) {
	switch (_fissionApproach)
	{
	case 0:
	return cutTwoDefault(cluster, finalhadrons, softUEisOn);
	break;
	case 1:
	return cutTwoNew(cluster, finalhadrons, softUEisOn);
	break;
	case 2:
	if (cluster->isBeamCluster()) {
	return cutTwoDefault(cluster, finalhadrons, softUEisOn);
	}
	else {
	return cutTwoNew(cluster, finalhadrons, softUEisOn);
	}
	break;
	case 3:
	{
	bool isPert=false;
	for (unsigned int i = 0; i < cluster->numComponents(); i++)
	{
	if (cluster->isPerturbative(i))
	isPert=true;
	}
	if (isPert)
	return cutTwoDefault(cluster, finalhadrons, softUEisOn);
	else
	return cutTwoNew(cluster, finalhadrons, softUEisOn);
	break;
	}
	case 4:
	{
	bool isPert=false;
	for (unsigned int i = 0; i < cluster->numComponents(); i++)
	{
	if (cluster->isPerturbative(i))
	isPert=true;
	}
	if (!isPert)
	return cutTwoDefault(cluster, finalhadrons, softUEisOn);
	else
	return cutTwoNew(cluster, finalhadrons, softUEisOn);
	break;
	}
	case 5:
	{
	bool isPert=true;
	for (unsigned int i = 0; i < cluster->numComponents(); i++)
	{
	if (!cluster->isPerturbative(i)) {
	isPert=false;
	break;
	}
	}
	if (isPert)
	return cutTwoDefault(cluster, finalhadrons, softUEisOn);
	else
	return cutTwoNew(cluster, finalhadrons, softUEisOn);
	break;
	}
	default:
	assert(false);
	}
	}

	ClusterFissioner::cutType
	ClusterFissioner::cutTwoDefault(ClusterPtr & cluster, tPVector & finalhadrons,
	bool softUEisOn) {
	// need to make sure only 2-cpt clusters get here
	assert(cluster->numComponents() == 2);
	tPPtr ptrQ1 = cluster->particle(0);
	tPPtr ptrQ2 = cluster->particle(1);
	Energy Mc = cluster->mass();
	assert(ptrQ1);
	assert(ptrQ2);

	// And check if those particles are from a beam remnant
	bool rem1 = cluster->isBeamRemnant(0);
	bool rem2 = cluster->isBeamRemnant(1);
	// workout which distribution to use
	bool soft1(false),soft2(false);
	switch (_iopRem) {
	case 0:
	soft1 = rem1 \|\| rem2;
	soft2 = rem2 \|\| rem1;
	break;
	case 1:
	soft1 = rem1;
	soft2 = rem2;
	break;
	}
	// Initialization for the exponential ("soft") mass distribution.
	static const int max_loop = 1000;
	int counter = 0;
	Energy Mc1 = ZERO, Mc2 = ZERO,m1=ZERO,m2=ZERO,m=ZERO;
	tcPDPtr toHadron1, toHadron2;
	PPtr newPtr1 = PPtr ();
	PPtr newPtr2 = PPtr ();
	bool succeeded = false;
	Lorentz5Momentum pClu1, pClu2, pQ1, pQone, pQtwo, pQ2;
	do
	{
	succeeded = false;
	++counter;
	// get a flavour for the qqbar pair
	drawNewFlavour(newPtr1,newPtr2,cluster);
	// check for right ordering
	assert (ptrQ2);
	assert (newPtr2);
	assert (ptrQ2->dataPtr());
	assert (newPtr2->dataPtr());
	if(cantMakeHadron(ptrQ1, newPtr1) \|\| cantMakeHadron(ptrQ2, newPtr2)) {
	swap(newPtr1, newPtr2);
	// check again
	if(cantMakeHadron(ptrQ1, newPtr1) \|\| cantMakeHadron(ptrQ2, newPtr2)) {
	throw Exception()
	<< "ClusterFissioner cannot split the cluster ("
	<< ptrQ1->PDGName() << ' ' << ptrQ2->PDGName()
	<< ") into hadrons.\n" << Exception::runerror;
	}
	}
	// Check that new clusters can produce particles and there is enough
	// phase space to choose the drawn flavour
	m1 = ptrQ1->data().constituentMass();
	m2 = ptrQ2->data().constituentMass();
	m = newPtr1->data().constituentMass();
	// Do not split in the case there is no phase space available
	if(Mc < m1+m + m2+m) continue;

	pQ1.setMass(m1);
	pQone.setMass(m);
	pQtwo.setMass(m);
	pQ2.setMass(m2);

	// pair<Energy,Energy> res = drawNewMasses(Mc, soft1, soft2, pClu1, pClu2,
	// ptrQ1, pQ1, newPtr1, pQone,
	// newPtr2, pQtwo, ptrQ2, pQ2);
	double weightMasses = drawNewMassesDefault(Mc, soft1, soft2, pClu1, pClu2,
	ptrQ1, pQ1, newPtr1, pQone,
	newPtr2, pQtwo, ptrQ2, pQ2);
	if (weightMasses==0.0)
	continue;

	// derive the masses of the children
	Mc1 = pClu1.mass();
	Mc2 = pClu2.mass();
	// static kinematic threshold
	if(_kinematicThresholdChoice == 0) {
	if (Mc1 < m1+m \|\| Mc2 < m+m2 \|\| Mc1+Mc2 > Mc) continue;
	if (_phaseSpaceWeights==2 &&
	( Mc1 < spectrum()->massLightestHadronPair(ptrQ1->dataPtr(),newPtr1->dataPtr())
	\|\| Mc2 < spectrum()->massLightestHadronPair(ptrQ2->dataPtr(),newPtr2->dataPtr()) ))
	continue;
	// dynamic kinematic threshold
	}
	else if(_kinematicThresholdChoice == 1) {
	bool C1 = ( sqr(Mc1) )/( sqr(m1) + sqr(m) + _kinThresholdShift ) < 1.0 ? true : false;
	bool C2 = ( sqr(Mc2) )/( sqr(m2) + sqr(m) + _kinThresholdShift ) < 1.0 ? true : false;
	bool C3 = ( sqr(Mc1) + sqr(Mc2) )/( sqr(Mc) ) > 1.0 ? true : false;

	if( C1 \|\| C2 \|\| C3 ) continue;
	}
	- if ( _phaseSpaceWeights && phaseSpaceVeto(Mc,Mc1,Mc2,m,m1,m2, ptrQ1, ptrQ2, newPtr1) ) {
	+ if ( _phaseSpaceWeights && phaseSpaceVeto(Mc,Mc1,Mc2,m,m1,m2, ptrQ1, ptrQ2, newPtr1, _powerLawPower) ) {
	// reduce counter as it regards only the mass sampling
	counter--;
	continue;
	}

	/**************************
	* New (not present in Fortran Herwig):
	* check whether the fragment masses Mc1 and Mc2 are above the
	* threshold for the production of the lightest pair of hadrons with the
	* right flavours. If not, then set by hand the mass to the lightest
	* single hadron with the right flavours, in order to solve correctly
	* the kinematics, and (later in this method) create directly such hadron
	* and add it to the children hadrons of the cluster that undergoes the
	* fission (i.e. the one pointed by iCluPtr). Notice that in this special
	* case, the heavy cluster that undergoes the fission has one single
	* cluster child and one single hadron child. We prefer this approach,
	* rather than to create a light cluster, with the mass set equal to
	* the lightest hadron, and let then the class LightClusterDecayer to do
	* the job to decay it to that single hadron, for two reasons:
	* First, because the sum of the masses of the two constituents can be,
	* in this case, greater than the mass of that hadron, hence it would
	* be impossible to solve the kinematics for such two components, and
	* therefore we would have a cluster whose components are undefined.
	* Second, the algorithm is faster, because it avoids the reshuffling
	* procedure that would be necessary if we used LightClusterDecayer
	* to decay the light cluster to the lightest hadron.
	****************************/
	// override chosen masses if needed
	toHadron1 = _hadronSpectrum->chooseSingleHadron(ptrQ1->dataPtr(), newPtr1->dataPtr(),Mc1);
	if(toHadron1) { Mc1 = toHadron1->mass(); pClu1.setMass(Mc1); }
	toHadron2 = _hadronSpectrum->chooseSingleHadron(ptrQ2->dataPtr(), newPtr2->dataPtr(),Mc2);
	if(toHadron2) { Mc2 = toHadron2->mass(); pClu2.setMass(Mc2); }
	// if a beam cluster not allowed to decay to hadrons
	if(cluster->isBeamCluster() && (toHadron1\|\|toHadron2) && softUEisOn)
	continue;
	// Check if the decay kinematics is still possible: if not then
	// force the one-hadron decay for the other cluster as well.
	if(Mc1 + Mc2 > Mc) {
	if(!toHadron1) {
	toHadron1 = _hadronSpectrum->chooseSingleHadron(ptrQ1->dataPtr(), newPtr1->dataPtr(),Mc-Mc2);
	if(toHadron1) { Mc1 = toHadron1->mass(); pClu1.setMass(Mc1); }
	}
	else if(!toHadron2) {
	toHadron2 = _hadronSpectrum->chooseSingleHadron(ptrQ2->dataPtr(), newPtr2->dataPtr(),Mc-Mc1);
	if(toHadron2) { Mc2 = toHadron2->mass(); pClu2.setMass(Mc2); }
	}
	}
	succeeded = (Mc >= Mc1+Mc2);
	}
	while (!succeeded && counter < max_loop);

	if(counter >= max_loop) {
	static const PPtr null = PPtr();
	return cutType(PPair(null,null),PPair(null,null));
	}

	// Determined the (5-components) momenta (all in the LAB frame)
	Lorentz5Momentum pClu = cluster->momentum(); // known
	Lorentz5Momentum p0Q1 = ptrQ1->momentum(); // known (mom Q1 before fission)
	calculateKinematics(pClu,p0Q1,toHadron1,toHadron2,
	pClu1,pClu2,pQ1,pQone,pQtwo,pQ2);

	/******************
	* The previous methods have determined the kinematics and positions
	* of C -> C1 + C2.
	* In the case that one of the two product is light, that means either
	* decayOneHadronClu1 or decayOneHadronClu2 is true, then the momenta
	* of the components of that light product have not been determined,
	* and a (light) cluster will not be created: the heavy father cluster
	* decays, in this case, into a single (not-light) cluster and a
	* single hadron. In the other, "normal", cases the father cluster
	* decays into two clusters, each of which has well defined components.
	* Notice that, in the case of components which point to particles, the
	* momenta of the components is properly set to the new values, whereas
	* we do not change the momenta of the pointed particles, because we
	* want to keep all of the information (that is the new momentum of a
	* component after the splitting, which is contained in the _momentum
	* member of the Component class, and the (old) momentum of that component
	* before the splitting, which is contained in the momentum of the
	* pointed particle). Please not make confusion of this only apparent
	* inconsistency!
	********************/
	LorentzPoint posC,pos1,pos2;
	posC = cluster->vertex();
	calculatePositions(pClu, posC, pClu1, pClu2, pos1, pos2);
	cutType rval;
	if(toHadron1) {
	rval.first = produceHadron(toHadron1, newPtr1, pClu1, pos1);
	finalhadrons.push_back(rval.first.first);
	}
	else {
	rval.first = produceCluster(ptrQ1, newPtr1, pClu1, pos1, pQ1, pQone, rem1);
	}
	if(toHadron2) {
	rval.second = produceHadron(toHadron2, newPtr2, pClu2, pos2);
	finalhadrons.push_back(rval.second.first);
	}
	else {
	rval.second = produceCluster(ptrQ2, newPtr2, pClu2, pos2, pQ2, pQtwo, rem2);
	}
	return rval;
	}

	namespace {
	int areNotSame(const Lorentz5Momentum & p1,const Lorentz5Momentum & p2){
	double tol=1e-5;
	if (abs(p1.vect().x()-p2.vect().x())>tol*GeV){
	std::cout << "disagreeing x " <<std::setprecision(-log10(tol)+2) << abs(p1.vect().x()-p2.vect().x())/GeV<< std::endl;
	return 1;
	}
	if (abs(p1.vect().y()-p2.vect().y())>tol*GeV){
	std::cout << "disagreeing y " <<std::setprecision(-log10(tol)+2)<< abs(p1.vect().y()-p2.vect().y())/GeV<< std::endl;
	return 2;
	}
	if (abs(p1.vect().z()-p2.vect().z())>tol*GeV){
	std::cout << "disagreeing z " <<std::setprecision(-log10(tol)+2)<< abs(p1.vect().z()-p2.vect().z())/GeV<< std::endl;
	return 3;
	}
	if (abs(p1.e()-p2.e())>tol*GeV){
	std::cout << "disagreeing e " <<std::setprecision(-log10(tol)+2)<< abs(p1.e()-p2.e())/GeV<< std::endl;
	return 4;
	}
	if (abs(p1.m()-p2.m())>tol*GeV){
	std::cout << "disagreeing m " <<std::setprecision(-log10(tol)+2)<< abs(p1.m()-p2.m())/GeV<< std::endl;
	return 4;
	}
	if (abs(p1.m()-p1.mass())>tol*GeV){
	std::cout << "disagreeing p1 m - mass " <<std::setprecision(-log10(tol)+2)<< abs(p1.m()-p1.mass())/GeV<< std::endl;
	return 4;
	}
	if (abs(p2.m()-p2.mass())>tol*GeV){
	std::cout << "disagreeing p2 m - mass " <<std::setprecision(-log10(tol)+2)<< abs(p2.m()-p2.mass())/GeV<< std::endl;
	return 4;
	}
	return 0;
	}
	}
	ClusterFissioner::cutType
	ClusterFissioner::cutTwoNew(ClusterPtr & cluster, tPVector & finalhadrons,
	bool softUEisOn) {
	// need to make sure only 2-cpt clusters get here
	assert(cluster->numComponents() == 2);
	tPPtr ptrQ1 = cluster->particle(0);
	tPPtr ptrQ2 = cluster->particle(1);
	Energy Mc = cluster->mass();
	// TODO BEGIN Changed comp to default
	if ( Mc < spectrum()->massLightestHadronPair(ptrQ1->dataPtr(),ptrQ2->dataPtr())) {
	static const PPtr null = PPtr();
	return cutType(PPair(null,null),PPair(null,null));
	}
	// TODO END Changed comp to default
	assert(ptrQ1);
	assert(ptrQ2);

	// And check if those particles are from a beam remnant
	bool rem1 = cluster->isBeamRemnant(0);
	bool rem2 = cluster->isBeamRemnant(1);
	// workout which distribution to use
	bool soft1(false),soft2(false);
	switch (_iopRem) {
	case 0:
	soft1 = rem1 \|\| rem2;
	soft2 = rem2 \|\| rem1;
	break;
	case 1:
	soft1 = rem1;
	soft2 = rem2;
	break;
	}
	// Initialization for the exponential ("soft") mass distribution.
	int counter_MEtry = 0;
	Energy Mc1 = ZERO;
	Energy Mc2 = ZERO;
	Energy m = ZERO;
	Energy m1 = ptrQ1->data().constituentMass();
	Energy m2 = ptrQ2->data().constituentMass();
	Energy mMin = getParticleData(ParticleID::d)->constituentMass();
	// Minimal threshold for non-zero Mass PhaseSpace
	if ( Mc < (m1 + m2 + 2*mMin )) {
	static const PPtr null = PPtr();
	return cutType(PPair(null,null),PPair(null,null));
	}
	tcPDPtr toHadron1, toHadron2;
	PPtr newPtr1 = PPtr ();
	PPtr newPtr2 = PPtr ();
	Lorentz5Momentum pClu1, pClu2, pQ1, pQone, pQtwo, pQ2;
	Lorentz5Momentum pClu = cluster->momentum(); // known
	const Lorentz5Momentum p0Q1 = ptrQ1->momentum(); // known (mom Q1 before fission)
	const Lorentz5Momentum p0Q2 = ptrQ2->momentum(); // known (mom Q2 before fission)
	// where to return to in case of rejected sample
	enum returnTo {
	FlavourSampling,
	MassSampling,
	PhaseSpaceSampling,
	MatrixElementSampling,
	Done
	};
	// start with FlavourSampling
	returnTo retTo=FlavourSampling;
	int letFissionToXHadrons = _allowHadronFinalStates;
	// if a beam cluster not allowed to decay to hadrons
	if (cluster->isBeamCluster() && softUEisOn) letFissionToXHadrons = 0;
	// ### Flavour, Mass, PhaseSpace and MatrixElement Sampling loop until accepted: ###
	bool escape = false;
	double weightMasses,weightPhaseSpaceAngles;
	do
	{
	switch (retTo)
	{
	case FlavourSampling:
	{
	// ## Flavour sampling and kinematic constraints ##
	drawNewFlavour(newPtr1,newPtr2,cluster);
	// get a flavour for the qqbar pair
	// check for right ordering
	assert (ptrQ2);
	assert (newPtr2);
	assert (ptrQ2->dataPtr());
	assert (newPtr2->dataPtr());
	// careful if DiquarkClusters can exist
	bool Q1diq = DiquarkMatcher::Check(ptrQ1->id());
	bool Q2diq = DiquarkMatcher::Check(ptrQ2->id());
	bool newQ1diq = DiquarkMatcher::Check(newPtr1->id());
	bool newQ2diq = DiquarkMatcher::Check(newPtr2->id());
	bool diqClu1 = Q1diq && newQ1diq;
	bool diqClu2 = Q2diq && newQ2diq;
	// DEBUG only:
	// std::cout << "Make Clusters: ( " << ptrQ1->PDGName() << " " << newPtr1->PDGName() << " ), ( "
	// << ptrQ2->PDGName() << " " << newPtr2->PDGName() << " )\n";
	// check if Hadron formation is possible
	if (!( diqClu1 \|\| diqClu2 )
	&& (cantMakeHadron(ptrQ1, newPtr1) \|\| cantMakeHadron(ptrQ2, newPtr2))) {
	swap(newPtr1, newPtr2);
	// check again
	if(cantMakeHadron(ptrQ1, newPtr1) \|\| cantMakeHadron(ptrQ2, newPtr2)) {
	throw Exception()
	<< "ClusterFissioner cannot split the cluster ("
	<< ptrQ1->PDGName() << ' ' << ptrQ2->PDGName()
	<< ") into hadrons.\n"
	<< "drawn Flavour: "<< newPtr1->PDGName()<<"\n"<< Exception::runerror;
	}
	}
	else if ( diqClu1 \|\| diqClu2 ){
	bool swapped=false;
	if ( !diqClu1 && cantMakeHadron(ptrQ1,newPtr1) ) {
	swap(newPtr1, newPtr2);
	swapped=true;
	}
	if ( !diqClu2 && cantMakeHadron(ptrQ2,newPtr2) ) {
	assert(!swapped);
	swap(newPtr1, newPtr2);
	}
	if ( diqClu2 && diqClu1 && ptrQ1->id()*newPtr1->id()>0 ) {
	assert(!swapped);
	swap(newPtr1, newPtr2);
	}
	if (!diqClu1) assert(!cantMakeHadron(ptrQ1,newPtr1));
	if (!diqClu2) assert(!cantMakeHadron(ptrQ2,newPtr2));
	}
	// Check that new clusters can produce particles and there is enough
	// phase space to choose the drawn flavour
	m = newPtr1->data().constituentMass();
	// Do not split in the case there is no phase space available + permille security
	double eps=0.1;
	if(Mc < (1+eps)*(m1 + m + m2 + m)) {
	retTo = FlavourSampling;
	// escape if no flavour phase space possibile without fission
	if (fabs((m - mMin)/GeV) < 1e-3) {
	escape = true;
	retTo = Done;
	}
	continue;
	}
	pQ1.setMass(m1);
	pQone.setMass(m);
	pQtwo.setMass(m);
	pQ2.setMass(m2);
	// Determined the (5-components) momenta (all in the LAB frame)
	// p0Q1.setMass(ptrQ1->mass()); // known (mom Q1 before fission)
	// p0Q1.rescaleEnergy(); // TODO check if needed
	// p0Q2.setMass(ptrQ2->mass()); // known (mom Q2 before fission)
	// p0Q2.rescaleEnergy();// TODO check if needed
	pClu.rescaleMass();

	Energy MLHP1 = spectrum()->hadronPairThreshold(ptrQ1->dataPtr(),newPtr1->dataPtr());
	Energy MLHP2 = spectrum()->hadronPairThreshold(ptrQ2->dataPtr(),newPtr2->dataPtr());

	Energy MLH1 = _hadronSpectrum->lightestHadron(ptrQ1->dataPtr(),newPtr1->dataPtr())->mass();
	Energy MLH2 = _hadronSpectrum->lightestHadron(ptrQ2->dataPtr(),newPtr2->dataPtr())->mass();

	bool canBeSingleHadron1 = (m1 + m) < MLHP1;
	bool canBeSingleHadron2 = (m2 + m) < MLHP2;

	Energy mThresh1 = (m1 + m);
	Energy mThresh2 = (m2 + m);

	if (canBeSingleHadron1) mThresh1 = MLHP1;
	if (canBeSingleHadron2) mThresh2 = MLHP2;

	switch (letFissionToXHadrons)
	{
	case 0:
	{
	// Option None: only C->C1,C2 allowed
	// check if mass phase space is non-zero
	// resample or escape if only allowed mass phase space is for C->H1,H2 or C->H1,C2
	if ( Mc < (mThresh1 + mThresh2)) {
	// escape if not even the lightest flavour phase space is possibile
	// TODO make this independent of explicit u/d quark
	if (
	fabs((m - getParticleData(ThePEG::ParticleID::u)->constituentMass())/GeV) < 1e-3
	\|\|
	fabs((m - getParticleData(ThePEG::ParticleID::d)->constituentMass())/GeV) < 1e-3
	) {
	escape = true;
	retTo = Done;
	continue;
	}
	else {
	retTo = FlavourSampling;
	continue;
	}
	}
	break;
	}
	case 1:
	{
	// Option SemiHadronicOnly: C->H,C allowed
	// NOTE: TODO implement matrix element for this
	// resample or escape if only allowed mass phase space is for C->H1,H2
	// First case is for ensuring the enough mass to be available and second one rejects disjoint mass regions
	if ( ( (canBeSingleHadron1 && canBeSingleHadron2)
	&& Mc < (mThresh1 + mThresh2) )
	\|\|
	( (canBeSingleHadron1 \|\| canBeSingleHadron2)
	&& (canBeSingleHadron1 ? Mc-(m2+m) < MLH1:false \|\| canBeSingleHadron2 ? Mc-(m1+m) < MLH2:false) )
	){
	// escape if not even the lightest flavour phase space is possibile
	// TODO make this independent of explicit u/d quark
	if (
	fabs((m - getParticleData(ThePEG::ParticleID::u)->constituentMass())/GeV) < 1e-3
	\|\|
	fabs((m - getParticleData(ThePEG::ParticleID::d)->constituentMass())/GeV) < 1e-3
	) {
	escape = true;
	retTo = Done;
	continue;
	}
	else {
	retTo = FlavourSampling;
	continue;
	}
	}
	break;
	}
	case 2:
	{
	// Option All: C->H,C and C->H,H allowed
	// NOTE: TODO implement matrix element for this
	// Mass Phase space for all option can always be found if cluster massive enough to go
	// to the lightest 2 hadrons
	break;
	}
	default:
	assert(false);
	}
	// Note: want to fallthrough (in C++17 one could uncomment
	// the line below to show that this is intentional)
	[[fallthrough]];
	}
	/*
	* MassSampling choices:
	* - Default (default)
	* - Uniform
	* - FlatPhaseSpace
	* - SoftMEPheno
	* */
	case MassSampling:
	{
	weightMasses = drawNewMasses(Mc, soft1, soft2, pClu1, pClu2,
	ptrQ1, pQ1, newPtr1, pQone,
	newPtr2, pQtwo, ptrQ2, pQ2);

	// TODO IF C->C1,C2 (and not C->C,H or H1,H2) masses sampled and is in PhaseSpace must push through
	// because otherwise no matrix element
	if(weightMasses==0.0) {
	// TODO check which option is better
	retTo = FlavourSampling;
	// retTo = MassSampling;
	continue;
	}

	// derive the masses of the children
	Mc1 = pClu1.mass();
	Mc2 = pClu2.mass();
	// static kinematic threshold
	if(_kinematicThresholdChoice != 1 )
	{
	// NOTE: that the squared thresholds below are
	// numerically more stringent which is needed
	// kaellen function calculation
	if(sqr(Mc1) < sqr(m1+m) \|\| sqr(Mc2) < sqr(m+m2) \|\| sqr(Mc1+Mc2) > sqr(Mc)){
	// TODO check which option is better
	// retTo = FlavourSampling;
	retTo = MassSampling;
	continue;
	}
	// dynamic kinematic threshold
	}
	else if(_kinematicThresholdChoice == 1) {
	bool C1 = ( sqr(Mc1) )/( sqr(m1) + sqr(m) + _kinThresholdShift ) < 1.0 ? true : false;
	bool C2 = ( sqr(Mc2) )/( sqr(m2) + sqr(m) + _kinThresholdShift ) < 1.0 ? true : false;
	bool C3 = ( sqr(Mc1) + sqr(Mc2) )/( sqr(Mc) ) > 1.0 ? true : false;

	if( C1 \|\| C2 \|\| C3 ) {
	// TODO check which option is better
	// retTo = FlavourSampling;
	retTo = MassSampling;
	continue;
	}
	}
	else
	assert(false);
	/**************************
	* New (not present in Fortran Herwig):
	* check whether the fragment masses Mc1 and Mc2 are above the
	* threshold for the production of the lightest pair of hadrons with the
	* right flavours. If not, then set by hand the mass to the lightest
	* single hadron with the right flavours, in order to solve correctly
	* the kinematics, and (later in this method) create directly such hadron
	* and add it to the children hadrons of the cluster that undergoes the
	* fission (i.e. the one pointed by iCluPtr). Notice that in this special
	* case, the heavy cluster that undergoes the fission has one single
	* cluster child and one single hadron child. We prefer this approach,
	* rather than to create a light cluster, with the mass set equal to
	* the lightest hadron, and let then the class LightClusterDecayer to do
	* the job to decay it to that single hadron, for two reasons:
	* First, because the sum of the masses of the two constituents can be,
	* in this case, greater than the mass of that hadron, hence it would
	* be impossible to solve the kinematics for such two components, and
	* therefore we would have a cluster whose components are undefined.
	* Second, the algorithm is faster, because it avoids the reshuffling
	* procedure that would be necessary if we used LightClusterDecayer
	* to decay the light cluster to the lightest hadron.
	****************************/
	// override chosen masses if needed
	toHadron1 = _hadronSpectrum->chooseSingleHadron(ptrQ1->dataPtr(), newPtr1->dataPtr(),Mc1);
	if ( letFissionToXHadrons == 0 && toHadron1 ) {
	// reject mass samples which would force C->C,H or C->H1,H2 fission if desired
	//
	// Check if Mc1max < MLHP1, in which case we might need to choose a different flavour
	// else resampling the masses should be sufficient
	Energy MLHP1 = _hadronSpectrum->massLightestHadronPair(ptrQ1->dataPtr(), newPtr1->dataPtr());
	Energy MLHP2 = _hadronSpectrum->massLightestHadronPair(ptrQ2->dataPtr(), newPtr2->dataPtr());
	Energy Mc1max = (m2+m) > MLHP2 ? (Mc-(m2+m)):(Mc-MLHP2);
	// for avoiding inf loops set min threshold
	if ( Mc1max - MLHP1 < 1e-2*GeV ) {
	retTo = FlavourSampling;
	}
	else {
	retTo = MassSampling;
	}
	continue;
	}
	if(toHadron1) {
	Mc1 = toHadron1->mass();
	pClu1.setMass(Mc1);
	}
	toHadron2 = _hadronSpectrum->chooseSingleHadron(ptrQ2->dataPtr(), newPtr2->dataPtr(),Mc2);
	if ( letFissionToXHadrons == 0 && toHadron2 ) {
	// reject mass samples which would force C->C,H or C->H1,H2 fission if desired
	//
	// Check if Mc2max < MLHP2, in which case we might need to choose a different flavour
	// else resampling the masses should be sufficient
	Energy MLHP1 = _hadronSpectrum->massLightestHadronPair(ptrQ1->dataPtr(), newPtr1->dataPtr());
	Energy MLHP2 = _hadronSpectrum->massLightestHadronPair(ptrQ2->dataPtr(), newPtr2->dataPtr());
	Energy Mc2max = (m1+m) > MLHP1 ? (Mc-(m1+m)):(Mc-MLHP1);
	// for avoiding inf loops set min threshold
	if ( Mc2max - MLHP2 < 1e-2*GeV ) {
	retTo = FlavourSampling;
	}
	else {
	retTo = MassSampling;
	}
	continue;
	}
	if(toHadron2) {
	Mc2 = toHadron2->mass();
	pClu2.setMass(Mc2);
	}
	if (letFissionToXHadrons == 1 && (toHadron1 && toHadron2) ) {
	// reject mass samples which would force C->H1,H2 fission if desired
	// TODO check which option is better
	// retTo = FlavourSampling;
	retTo = MassSampling;
	continue;
	}
	// Check if the decay kinematics is still possible: if not then
	// force the one-hadron decay for the other cluster as well.
	// NOTE: that the squared thresholds below are
	// numerically more stringent which is needed
	// kaellen function calculation
	if(sqr(Mc1 + Mc2) > sqr(Mc)) {
	// reject if we would need to create a C->H,H process if letFissionToXHadrons<2
	if (letFissionToXHadrons < 2) {
	// TODO check which option is better
	// retTo = FlavourSampling;
	retTo = MassSampling;
	continue;
	}

	// TODO forbid other cluster!!!! to be also at hadron
	if(!toHadron1) {
	toHadron1 = _hadronSpectrum->chooseSingleHadron(ptrQ1->dataPtr(), newPtr1->dataPtr(),Mc-Mc2);
	// toHadron1 = _hadronSpectrum->chooseSingleHadron(ptrQ1->dataPtr(), newPtr1->dataPtr(),ZERO);
	if(toHadron1) {
	Mc1 = toHadron1->mass();
	pClu1.setMass(Mc1);
	}
	}
	else if(!toHadron2) {
	toHadron2 = _hadronSpectrum->chooseSingleHadron(ptrQ2->dataPtr(), newPtr2->dataPtr(),Mc-Mc1);
	// toHadron2 = _hadronSpectrum->chooseSingleHadron(ptrQ2->dataPtr(), newPtr2->dataPtr(),ZERO);
	if(toHadron2) {
	Mc2 = toHadron2->mass();
	pClu2.setMass(Mc2);
	}
	}
	}
	// NOTE: that the squared thresholds below are
	// numerically more stringent which is needed
	// kaellen function calculation
	if (sqr(Mc) <= sqr(Mc1+Mc2)){
	// escape if already lightest quark drawn
	if (
	fabs((m - getParticleData(ThePEG::ParticleID::u)->constituentMass())/GeV) < 1e-3
	\|\|
	fabs((m - getParticleData(ThePEG::ParticleID::d)->constituentMass())/GeV) < 1e-3
	) {
	// escape = true;
	retTo = Done;
	}
	else {
	// Try again with lighter quark
	retTo = FlavourSampling;
	}
	continue;
	}
	// Note: want to fallthrough (in C++17 one could uncomment
	// the line below to show that this is intentional)
	[[fallthrough]];
	}
	/*
	* PhaseSpaceSampling choices:
	* - FullyAligned (default)
	* - AlignedIsotropic
	* - FullyIsotropic
	* */
	case PhaseSpaceSampling:
	{
	// ### Sample the Phase Space with respect to Matrix Element: ###
	// TODO insert here PhaseSpace sampler
	weightPhaseSpaceAngles = drawKinematics(pClu,p0Q1,p0Q2,toHadron1,toHadron2,
	pClu1,pClu2,pQ1,pQone,pQtwo,pQ2);
	{ //sanity
	if (areNotSame(pClu1+pClu2,pClu)) {
	std::cout << "PCLU" << std::endl;
	}
	if (areNotSame(pClu1,pQ1+pQone)) {
	std::cout << "PCLU1" << std::endl;
	}
	if (areNotSame(pClu2,pQ2+pQtwo)) {
	std::cout << "PCLU2" << std::endl;
	}
	if (areNotSame(pClu,pQ1+pQone+pQ2+pQtwo)) {
	std::cout << "PTOT" << std::endl;
	}
	}

	if(weightPhaseSpaceAngles==0.0) {
	// TODO check which option is better
	retTo = MassSampling;
	continue;
	}
	// Should be precise i.e. no rejection expected
	// Note: want to fallthrough (in C++17 one could uncomment
	// the line below to show that this is intentional)
	[[fallthrough]];
	}
	/*
	* MatrixElementSampling choices:
	* - Default (default)
	* - SoftQQbar
	* */
	case MatrixElementSampling:
	{
	counter_MEtry++;
	// TODO maybe bridge this to work more neatly
	// Ignore matrix element for C->C,H or C->H1,H2 fission
	if (toHadron1 \|\| toHadron2) {
	retTo = Done;
	break;
	}
	// Actual MatrixElement evaluated at sampled PhaseSpace point
	double SQME = calculateSQME(p0Q1,p0Q2,pQ1,pQone,pQ2,pQtwo);
	// Total overestimate of MatrixElement independent
	// of the PhaseSpace point and independent of M1,M2
	double SQMEoverEstimate = calculateSQME_OverEstimate(Mc,m1,m2,m);
	// weight for MatrixElement*PhaseSpace must be in [0:1]
	double weightSQME = SQME/SQMEoverEstimate;
	assert(weightSQME>0.0);
	// weight(M1,M2) for M1M2(Two body PhaseSpace)**3 should be in [0,1]
	double weightFlatPS = weightFlatPhaseSpace(Mc, Mc1, Mc2, m, m1, m2, ptrQ1, ptrQ2, newPtr1);
	if (weightFlatPS<0 \|\| weightFlatPS>1.0){
	throw Exception() << "weightFlatPS = "<< weightFlatPS << " > 1 or negative in ClusterFissioner::cutTwo"
	<< "Mc = " << Mc/GeV
	<< "Mc1 = " << Mc1/GeV
	<< "Mc2 = " << Mc2/GeV
	<< "m1 = " << m1/GeV
	<< "m2 = " << m2/GeV
	<< "m = " << m/GeV
	<< "SQME = " << SQME
	<< "SQME_OE = " << SQMEoverEstimate
	<< "PS = " << weightFlatPS
	<< Exception::runerror;
	}
	// current phase space point is distributed according to weightSamp
	double Pacc = weightFlatPS * weightSQME;
	double SamplingWeight = weightMasses * weightPhaseSpaceAngles;
	if (!(SamplingWeight >= 0 ) \|\| std::isnan(SamplingWeight) \|\| std::isinf(SamplingWeight)){
	throw Exception() << "SamplingWeight = "<< SamplingWeight << " in ClusterFissioner::cutTwo"
	<< "Mc = " << Mc/GeV
	<< "Mc1 = " << Mc1/GeV
	<< "Mc2 = " << Mc2/GeV
	<< "m1 = " << m1/GeV
	<< "m2 = " << m2/GeV
	<< "m = " << m/GeV
	<< "weightMasses = " << weightMasses
	<< "weightPhaseSpaceAngles = " << weightPhaseSpaceAngles
	<< "PS = " << weightFlatPS
	<< Exception::runerror;
	}
	Pacc/=SamplingWeight;
	if (!(Pacc >= 0 ) \|\| std::isnan(Pacc) \|\| std::isinf(Pacc)){
	throw Exception() << "Pacc = "<< Pacc << " < 0 in ClusterFissioner::cutTwo"
	<< "Mc = " << Mc/GeV
	<< "Mc1 = " << Mc1/GeV
	<< "Mc2 = " << Mc2/GeV
	<< "m1 = " << m1/GeV
	<< "m2 = " << m2/GeV
	<< "m = " << m/GeV
	<< "SQME = " << SQME
	<< "SQME_OE = " << SQMEoverEstimate
	<< "PS = " << weightFlatPS
	<< Exception::runerror;
	}
	assert(Pacc >= 0.0);
	if (_matrixElement!=0) {
	if ( Pacc > 1.0){
	countPaccGreater1();
	throw Exception() << "Pacc = "<< Pacc << " > 1 in ClusterFissioner::cutTwo"
	<< "Mc = " << Mc/GeV
	<< "Mc1 = " << Mc1/GeV
	<< "Mc2 = " << Mc2/GeV
	<< "m1 = " << m1/GeV
	<< "m2 = " << m2/GeV
	<< "m = " << m/GeV
	<< Exception::warning;
	}
	}
	static int first=_writeOut;
	if (_matrixElement==0 \|\| UseRandom::rnd()<Pacc) { //Always accept a sample if trivial matrix element
	if (_writeOut) {
	// std::cout << "\nAccept Pacc = "<<Pacc<<"\n";
	std::ofstream out("data_CluFis.dat", std::ios::app \| std::ios::out);
	Lorentz5Momentum p0Q1tmp(p0Q1);
	Lorentz5Momentum pClu1tmp(pClu1);
	Lorentz5Momentum pClu2tmp(pClu2);
	p0Q1tmp.boost(-pClu.boostVector());
	pClu1tmp.boost(-pClu.boostVector());
	double cosThetaP1C1=p0Q1tmp.vect().cosTheta(pClu1tmp.vect());
	out << Pacc << "\t"
	<< Mc/GeV << "\t"
	<< pClu1.mass()/GeV << "\t"
	<< pClu2.mass()/GeV << "\t"
	<< pQonepQtwo/(pQone.mass()pQtwo.mass()) << "\t"
	<< pQ1pQ2/(pQ1.mass()pQ2.mass()) << "\t"
	<< pQ1pQtwo/(pQ1.mass()pQtwo.mass()) << "\t"
	<< pQ2pQone/(pQ2.mass()pQone.mass()) << "\t"
	<< p0Q1(pQone+pQtwo)/(p0Q1.mass()(pQone.mass()+pQtwo.mass())) << "\t"
	<< p0Q2(pQone+pQtwo)/(p0Q2.mass()(pQone.mass()+pQtwo.mass())) << "\t"
	<< m/GeV << "\t"
	<< cosThetaP1C1
	<< "\n";
	out.close();
	Energy MbinStart=2.0*GeV;
	Energy MbinEnd=91.0*GeV;
	Energy dMbin=1.0*GeV;
	Energy MbinLow;
	Energy MbinHig;
	if ( fabs((m-getParticleData(ParticleID::d)->constituentMass())/GeV)<1e-5
	&& fabs((m1-getParticleData(ParticleID::d)->constituentMass())/GeV)<1e-5
	&& fabs((m2-getParticleData(ParticleID::d)->constituentMass())/GeV)<1e-5) {
	int cnt = 0;
	for (Energy MbinIt=MbinStart; MbinIt < MbinEnd; MbinIt+=dMbin)
	{
	if (first){
	first=0;
	std::cout << "\nFirst\n" << std::endl;
	int ctr=0;
	for (Energy MbinIt=MbinStart; MbinIt < MbinEnd; MbinIt+=dMbin)
	{
	std::string name = "WriteOut/data_CluFis_BinM_"+std::to_string(ctr)+".dat";
	std::ofstream out2(name,std::ios::out);
	out2 << "# Binned from "<<MbinIt/GeV << "\tto\t" << (MbinIt+dMbin)/GeV<<"\n";
	out2 << "# m=m1=m2=0.325\n";
	out2 << "# (1) Pacc\t"
	<< "(2) Mc/GeV\t"
	<< "(3) M1/GeV\t"
	<< "(4) M2/GeV\t"
	<< "(5) q.qbar/(mq^2)\t"
	<< "(6) q1.q2/(m1*m2)\t"
	<< "(7) q1.qbar/(m1*mq)\t"
	<< "(8) q2.q/(m2*mq)\t"
	<< "(9) p1.(q+qbar)/(m12mq)\t"
	<< "(10) p2.(q+qbar)/(m22mq)\t"
	<< "(11) m/GeV\t"
	<< "(12) cos(theta_p1_Q1)\n";
	out2.close();
	ctr++;
	}
	// first=0;
	}
	MbinLow = MbinIt;
	MbinHig = MbinLow+dMbin;
	if (Mc>MbinLow && Mc<MbinHig) {
	std::string name = "WriteOut/data_CluFis_BinM_"+std::to_string(cnt)+".dat";
	std::ofstream out2(name, std::ios::app );
	Lorentz5Momentum p0Q1tmp(p0Q1);
	Lorentz5Momentum pClu1tmp(pClu1);
	Lorentz5Momentum pClu2tmp(pClu2);
	p0Q1tmp.boost(-pClu.boostVector());
	pClu1tmp.boost(-pClu.boostVector());
	double cosThetaP1C1=p0Q1tmp.vect().cosTheta(pClu1tmp.vect());
	out2 << Pacc << "\t"
	<< Mc/GeV << "\t"
	<< pClu1.mass()/GeV << "\t"
	<< pClu2.mass()/GeV << "\t"
	<< pQonepQtwo/(pQone.mass()pQtwo.mass()) << "\t"
	<< pQ1pQ2/(pQ1.mass()pQ2.mass()) << "\t"
	<< pQ1pQtwo/(pQ1.mass()pQtwo.mass()) << "\t"
	<< pQ2pQone/(pQ2.mass()pQone.mass()) << "\t"
	<< p0Q1(pQone+pQtwo)/(p0Q1.mass()(pQone.mass()+pQtwo.mass())) << "\t"
	<< p0Q2(pQone+pQtwo)/(p0Q2.mass()(pQone.mass()+pQtwo.mass())) << "\t"
	<< m/GeV << "\t"
	<< cosThetaP1C1
	<< "\n";
	out2.close();
	break;
	}
	// if (Mc<MbinIt) break;
	cnt++;
	}
	}
	}
	retTo=Done;
	break;
	}
	retTo = MassSampling;
	continue;
	}
	default:
	{
	assert(false);
	}
	}
	}
	while (retTo!=Done && !escape && counter_MEtry < _maxLoopFissionMatrixElement);

	if(escape) {
	// happens if we get at too light cluster to begin with
	static const PPtr null = PPtr();
	return cutType(PPair(null,null),PPair(null,null));
	}

	if(counter_MEtry >= _maxLoopFissionMatrixElement) {
	countMaxLoopViolations();
	// happens if we get too massive clusters where the matrix element
	// is very inefficiently sampled
	std::stringstream warning;
	warning
	<< "Matrix Element rejection sampling tried more than "
	<< counter_MEtry
	<< " times.\nMcluster = " << Mc/GeV
	<< "GeV\nm1 = " << m1/GeV
	<< "GeV\nm2 = " << m2/GeV
	<< "GeV\nm = " << m/GeV
	<< "GeV\n"
	<< "isBeamCluster = " << cluster->isBeamCluster()
	<<"Using default as ClusterFissioner::cutTwoDefault as a fallback.\n"
	<< "*This Exception should not happen too often!* ";
	generator()->logWarning( Exception(warning.str(),Exception::warning));
	switch (_failModeMaxLoopMatrixElement)
	{
	case 0:
	// OldFission
	return cutTwoDefault(cluster,finalhadrons,softUEisOn);
	break;
	case 1:
	// RejectEvent
	throw Exception() << warning.str()
	<< Exception::eventerror;
	break;
	default:
	assert(false);
	}
	}
	assert(abs(Mc1-pClu1.m())<1e-2*GeV);
	assert(abs(Mc2-pClu2.m())<1e-2*GeV);
	assert(abs(Mc1-pClu1.mass())<1e-2*GeV);
	assert(abs(Mc2-pClu2.mass())<1e-2*GeV);
	{
	Lorentz5Momentum pp1(p0Q1);
	Lorentz5Momentum pclu1(pClu1);
	pclu1.boost(-pClu.boostVector());
	pp1.boost(-pClu.boostVector());
	// std::ofstream out("testCF.dat", std::ios::app);
	// out << pclu1.vect().cosTheta(pp1.vect())<<"\n";
	// out.close();
	}
	// Count successfull ClusterFission
	countFissionMatrixElement();
	// ==> full sample generated
	/******************
	* The previous methods have determined the kinematics and positions
	* of C -> C1 + C2.
	* In the case that one of the two product is light, that means either
	* decayOneHadronClu1 or decayOneHadronClu2 is true, then the momenta
	* of the components of that light product have not been determined,
	* and a (light) cluster will not be created: the heavy father cluster
	* decays, in this case, into a single (not-light) cluster and a
	* single hadron. In the other, "normal", cases the father cluster
	* decays into two clusters, each of which has well defined components.
	* Notice that, in the case of components which point to particles, the
	* momenta of the components is properly set to the new values, whereas
	* we do not change the momenta of the pointed particles, because we
	* want to keep all of the information (that is the new momentum of a
	* component after the splitting, which is contained in the _momentum
	* member of the Component class, and the (old) momentum of that component
	* before the splitting, which is contained in the momentum of the
	* pointed particle). Please not make confusion of this only apparent
	* inconsistency!
	********************/
	LorentzPoint posC,pos1,pos2;
	posC = cluster->vertex();
	calculatePositions(pClu, posC, pClu1, pClu2, pos1, pos2);
	cutType rval;
	if(toHadron1) {
	rval.first = produceHadron(toHadron1, newPtr1, pClu1, pos1);
	finalhadrons.push_back(rval.first.first);
	}
	else {
	rval.first = produceCluster(ptrQ1, newPtr1, pClu1, pos1, pQ1, pQone, rem1);
	}
	if(toHadron2) {
	rval.second = produceHadron(toHadron2, newPtr2, pClu2, pos2);
	finalhadrons.push_back(rval.second.first);
	}
	else {
	rval.second = produceCluster(ptrQ2, newPtr2, pClu2, pos2, pQ2, pQtwo, rem2);
	}
	return rval;
	}

	ClusterFissioner::cutType
	ClusterFissioner::cutThree(ClusterPtr & cluster, tPVector & finalhadrons,
	bool softUEisOn) {
	// need to make sure only 3-cpt clusters get here
	assert(cluster->numComponents() == 3);
	// extract quarks
	tPPtr ptrQ[3] = {cluster->particle(0),cluster->particle(1),cluster->particle(2)};
	assert( ptrQ[0] && ptrQ[1] && ptrQ[2] );
	// find maximum mass pair
	Energy mmax(ZERO);
	Lorentz5Momentum pDiQuark;
	int iq1(-1),iq2(-1);
	Lorentz5Momentum psum;
	for(int q1=0;q1<3;++q1) {
	psum+= ptrQ[q1]->momentum();
	for(int q2=q1+1;q2<3;++q2) {
	Lorentz5Momentum ptest = ptrQ[q1]->momentum()+ptrQ[q2]->momentum();
	ptest.rescaleMass();
	Energy mass = ptest.m();
	if(mass>mmax) {
	mmax = mass;
	pDiQuark = ptest;
	iq1 = q1;
	iq2 = q2;
	}
	}
	}
	// and the spectators
	int iother(-1);
	for(int ix=0;ix<3;++ix) if(ix!=iq1&&ix!=iq2) iother=ix;
	assert(iq1>=0&&iq2>=0&&iother>=0);

	// And check if those particles are from a beam remnant
	bool rem1 = cluster->isBeamRemnant(iq1);
	bool rem2 = cluster->isBeamRemnant(iq2);
	// workout which distribution to use
	bool soft1(false),soft2(false);
	switch (_iopRem) {
	case 0:
	soft1 = rem1 \|\| rem2;
	soft2 = rem2 \|\| rem1;
	break;
	case 1:
	soft1 = rem1;
	soft2 = rem2;
	break;
	}
	// Initialization for the exponential ("soft") mass distribution.
	static const int max_loop = 1000;
	int counter = 0;
	Energy Mc1 = ZERO, Mc2 = ZERO, m1=ZERO, m2=ZERO, m=ZERO;
	tcPDPtr toHadron;
	bool toDiQuark(false);
	PPtr newPtr1 = PPtr(),newPtr2 = PPtr();
	PDPtr diquark;
	bool succeeded = false;
	Lorentz5Momentum pClu1, pClu2, pQ1, pQone, pQtwo, pQ2;
	do {
	succeeded = false;
	++counter;

	// get a flavour for the qqbar pair
	drawNewFlavour(newPtr1,newPtr2,cluster);

	// randomly pick which will be (anti)diquark and which a mesonic cluster
	if(UseRandom::rndbool()) {
	swap(iq1,iq2);
	swap(rem1,rem2);
	}
	// check first order
	if(cantMakeHadron(ptrQ[iq1], newPtr1) \|\| !spectrum()->canMakeDiQuark(ptrQ[iq2], newPtr2)) {
	swap(newPtr1,newPtr2);
	}
	// check again
	if(cantMakeHadron(ptrQ[iq1], newPtr1) \|\| !spectrum()->canMakeDiQuark(ptrQ[iq2], newPtr2)) {
	throw Exception()
	<< "ClusterFissioner cannot split the cluster ("
	<< ptrQ[iq1]->PDGName() << ' ' << ptrQ[iq2]->PDGName()
	<< ") into a hadron and diquark.\n" << Exception::runerror;
	}
	// Check that new clusters can produce particles and there is enough
	// phase space to choose the drawn flavour
	m1 = ptrQ[iq1]->data().constituentMass();
	m2 = ptrQ[iq2]->data().constituentMass();
	m = newPtr1->data().constituentMass();
	// Do not split in the case there is no phase space available
	if(mmax < m1+m + m2+m) continue;

	pQ1.setMass(m1);
	pQone.setMass(m);
	pQtwo.setMass(m);
	pQ2.setMass(m2);

	double weightMasses = drawNewMassesDefault(mmax, soft1, soft2, pClu1, pClu2,
	ptrQ[iq1], pQ1, newPtr1, pQone,
	newPtr2, pQtwo, ptrQ[iq1], pQ2);

	if (weightMasses == 0.0) continue;

	Mc1 = pClu1.mass();
	Mc2 = pClu2.mass();

	if(Mc1 < m1+m \|\| Mc2 < m+m2 \|\| Mc1+Mc2 > mmax) continue;

	if ( _phaseSpaceWeights && phaseSpaceVeto(mmax,Mc1,Mc2,m,m1,m2) ) {
	// reduce counter as it regards only the mass sampling
	counter--;
	continue;
	}

	// check if need to force meson clster to hadron
	toHadron = _hadronSpectrum->chooseSingleHadron(ptrQ[iq1]->dataPtr(), newPtr1->dataPtr(),Mc1);
	if(toHadron) { Mc1 = toHadron->mass(); pClu1.setMass(Mc1); }
	// check if need to force diquark cluster to be on-shell
	toDiQuark = false;
	diquark = spectrum()->makeDiquark(ptrQ[iq2]->dataPtr(), newPtr2->dataPtr());
	if(Mc2 < diquark->constituentMass()) {
	Mc2 = diquark->constituentMass(); pClu2.setMass(Mc2);
	toDiQuark = true;
	}
	// if a beam cluster not allowed to decay to hadrons
	if(cluster->isBeamCluster() && toHadron && softUEisOn)
	continue;
	// Check if the decay kinematics is still possible: if not then
	// force the one-hadron decay for the other cluster as well.
	if(Mc1 + Mc2 > mmax) {
	if(!toHadron) {
	toHadron = _hadronSpectrum->chooseSingleHadron(ptrQ[iq1]->dataPtr(), newPtr1->dataPtr(),mmax-Mc2);
	if(toHadron) { Mc1 = toHadron->mass(); pClu1.setMass(Mc1); }
	}
	else if(!toDiQuark) {
	Mc2 = _hadronSpectrum->massLightestHadron(ptrQ[iq2]->dataPtr(), newPtr2->dataPtr()); pClu2.setMass(Mc2);
	toDiQuark = true;
	}
	}
	succeeded = (mmax >= Mc1+Mc2);
	}
	while (!succeeded && counter < max_loop);
	// check no of tries
	if(counter >= max_loop) return cutType();

	// Determine the (5-components) momenta (all in the LAB frame)
	Lorentz5Momentum p0Q1 = ptrQ[iq1]->momentum();
	calculateKinematics(pDiQuark,p0Q1,toHadron,toDiQuark,
	pClu1,pClu2,pQ1,pQone,pQtwo,pQ2);
	// positions of the new clusters
	LorentzPoint pos1,pos2;
	Lorentz5Momentum pBaryon = pClu2+ptrQ[iother]->momentum();
	calculatePositions(cluster->momentum(), cluster->vertex(), pClu1, pBaryon, pos1, pos2);
	// first the mesonic cluster/meson
	cutType rval;
	if(toHadron) {
	rval.first = produceHadron(toHadron, newPtr1, pClu1, pos1);
	finalhadrons.push_back(rval.first.first);
	}
	else {
	rval.first = produceCluster(ptrQ[iq1], newPtr1, pClu1, pos1, pQ1, pQone, rem1);
	}
	if(toDiQuark) {
	rem2 \|= cluster->isBeamRemnant(iother);
	PPtr newDiQuark = diquark->produceParticle(pClu2);
	rval.second = produceCluster(newDiQuark, ptrQ[iother], pBaryon, pos2, pClu2,
	ptrQ[iother]->momentum(), rem2);
	}
	else {
	rval.second = produceCluster(ptrQ[iq2], newPtr2, pBaryon, pos2, pQ2, pQtwo, rem2,
	ptrQ[iother],cluster->isBeamRemnant(iother));
	}
	cluster->isAvailable(false);
	return rval;
	}

	ClusterFissioner::PPair
	ClusterFissioner::produceHadron(tcPDPtr hadron, tPPtr newPtr, const Lorentz5Momentum &a,
	const LorentzPoint &b) const {
	PPair rval;
	if(hadron->coloured()) {
	rval.first = (_hadronSpectrum->lightestHadron(hadron,newPtr->dataPtr()))->produceParticle();
	}
	else
	rval.first = hadron->produceParticle();
	rval.second = newPtr;
	rval.first->set5Momentum(a);
	rval.first->setVertex(b);
	return rval;
	}

	ClusterFissioner::PPair ClusterFissioner::produceCluster(tPPtr ptrQ, tPPtr newPtr,
	const Lorentz5Momentum & a,
	const LorentzPoint & b,
	const Lorentz5Momentum & c,
	const Lorentz5Momentum & d,
	bool isRem,
	tPPtr spect, bool remSpect) const {
	PPair rval;
	rval.second = newPtr;
	ClusterPtr cluster = !spect ? new_ptr(Cluster(ptrQ,rval.second)) : new_ptr(Cluster(ptrQ,rval.second,spect));
	rval.first = cluster;
	cluster->set5Momentum(a);
	cluster->setVertex(b);
	assert(cluster->particle(0)->id() == ptrQ->id());
	cluster->particle(0)->set5Momentum(c);
	cluster->particle(1)->set5Momentum(d);
	cluster->setBeamRemnant(0,isRem);
	if(remSpect) cluster->setBeamRemnant(2,remSpect);
	return rval;
	}

	/**
	* Calculate the phase space weight for M1M2(2 body PhaseSpace) ignore constituent masses
	*/
	double ClusterFissioner::weightFlatPhaseSpaceNoConstituentMasses(const Energy Mc, const Energy Mc1, const Energy Mc2) const {
	double M_temp = Mc/GeV;
	double M1_temp = Mc1/GeV;
	double M2_temp = Mc2/GeV;
	if (sqr(M_temp)<sqr(M1_temp+M2_temp)) {
	// This should be checked before
	throw Exception()
	<< "ClusterFissioner has not checked Masses properly\n"
	<< "Mc = " << M_temp << "\n"
	<< "Mc1 = " << M1_temp << "\n"
	<< "Mc2 = " << M2_temp << "\n"
	<< Exception::warning;
	return 0.0;
	}
	double lam = Kinematics::kaellen(M_temp, M1_temp, M2_temp);
	double ratio;
	// old weight of Jan without the Jacobi factor M1*M2 of the Mass integration
	// double PSweight = pow(lam,_dim-3.)pow(M1_tempM2_temp,2.-_dim);
	// new weight with the Jacobi factor M1*M2 of the Mass integration
	double PSweight = M1_tempM2_temppow(sqrt(lam),_dim-3.);
	// overestimate only possible for dim>=3.0
	assert(_dim>=3.0);
	// old overestimate of Jan without the Jacobi factor M1*M2 of the Mass integration
	// double overEstimate = _dim>=4.0 ? pow(M_temp,4._dim-14.):pow(M_temp,2(_dim-3.0))/pow((m1_temp+m_temp)*(m2_temp+m_temp),4.0-_dim);
	// new improved overestimate with the Jacobi factor M1*M2 of the Mass integration
	double overEstimate = pow(6.0sqrt(3.0), 3.0 - _dim)pow(M_temp, 2.*(_dim-2.));
	ratio = PSweight/overEstimate;
	if (!(ratio>=0)) std::cout << "ratio = " <<ratio<<" M "<<M_temp<<" M1 "<<M1_temp<<" M2 "<<M2_temp<<"\t"<<_dim<<"\t" << lam <<"\t"<< overEstimate<<"\n\n";
	if (!(ratio<=1)) std::cout << "ratio = " <<ratio<<" M "<<M_temp<<" M1 "<<M1_temp<<" M2 "<<M2_temp<<"\t"<<_dim<<"\t" << lam <<"\t"<< overEstimate<<"\n\n";
	// if (ratio > 0.9) std::cout << "ratio = " << ratio <<"\n";
	assert (ratio >= 0);
	assert (ratio <= 1);
	return ratio;
	}
	/**
	* Calculate the phase space weight for M1M2(2 body PhaseSpace)^3
	*/
	-double ClusterFissioner::weightFlatPhaseSpaceConstituentMasses(const Energy Mc, const Energy Mc1, const Energy Mc2,
	- const Energy m, const Energy m1, const Energy m2) const {
	+double ClusterFissioner::weightPhaseSpaceConstituentMasses(const Energy Mc, const Energy Mc1, const Energy Mc2,
	+ const Energy m, const Energy m1, const Energy m2, const double power) const {
	double M_temp = Mc/GeV;
	double M1_temp = Mc1/GeV;
	double M2_temp = Mc2/GeV;
	double m_temp = m/GeV;
	double m1_temp = m1/GeV;
	double m2_temp = m2/GeV;
	if (sqr(M_temp)<sqr(M1_temp+M2_temp)
	\|\| sqr(M1_temp)<sqr(m1_temp+m_temp)
	\|\| sqr(M2_temp)<sqr(m2_temp+m_temp)
	) {
	// This should be checked before
	throw Exception()
	<< "ClusterFissioner has not checked Masses properly\n"
	<< "Mc = " << M_temp << "\n"
	<< "Mc1 = " << M1_temp << "\n"
	<< "Mc2 = " << M2_temp << "\n"
	<< "m1 = " << m1_temp << "\n"
	<< "m2 = " << m2_temp << "\n"
	<< "m = " << m_temp << "\n"
	<< Exception::warning;
	return 0.0;
	}
	double lam1 = Kinematics::kaellen(M1_temp, m1_temp, m_temp);
	double lam2 = Kinematics::kaellen(M2_temp, m2_temp, m_temp);
	double lam3 = Kinematics::kaellen(M_temp, M1_temp, M2_temp);
	double ratio;
	// old weight of Jan without the Jacobi factor M1*M2 of the Mass integration
	// double PSweight = pow(lam1lam2lam3,_dim-3.)pow(M1_tempM2_temp,2.-_dim);
	// new weight with the Jacobi factor M1*M2 of the Mass integration
	double PSweight = pow(lam1lam2lam3,(_dim-3.)/2.0)pow(M1_tempM2_temp,3.-_dim);
	// overestimate only possible for dim>=3.0
	assert(_dim>=3.0);
	// old overestimate of Jan without the Jacobi factor M1*M2 of the Mass integration
	// double overEstimate = _dim>=4.0 ? pow(M_temp,4._dim-14.):pow(M_temp,2(_dim-3.0))/pow((m1_temp+m_temp)*(m2_temp+m_temp),4.0-_dim);
	// new improved overestimate with the Jacobi factor M1*M2 of the Mass integration
	double overEstimate = pow(6.0sqrt(3.0), 3.0 - _dim)pow(M_temp, 4.*_dim-12.);
	ratio = PSweight/overEstimate;
	if (!(ratio>=0)) std::cout << "ratio = " <<ratio<<" M "<<M_temp<<" M1 "<<M1_temp<<" M2 "<<M2_temp<<" m1 "<<m1_temp<<" m2 "<<m2_temp<<" m "<<m_temp<<"\t"<<_dim<<"\t" << lam1<<"\t"<< lam2<<"\t" << lam3 <<"\t"<< overEstimate<<"\n\n";
	+ if (power) {
	+ double powerLawOver = power<0 ? pow(Mc1Mc2/((m1+m)(m2+m)),power):pow(Mc1Mc2/((Mc-(m1+m))(Mc-(m2+m))),power);
	+ ratio*=powerLawOver;
	+ }
	// if (ratio > 0.9) std::cout << "ratio = " << ratio <<"\n";
	assert (ratio >= 0);
	assert (ratio <= 1);
	return ratio;
	}
	/**
	* Calculate the phase space weight for M1M2(2 body PhaseSpace)^3
	* using Hadron Masses
	*/
	double ClusterFissioner::weightFlatPhaseSpaceHadronMasses(const Energy Mc, const Energy Mc1, const Energy Mc2, tcPPtr pQ, tcPPtr pQ1, tcPPtr pQ2) const {
	auto LHP1 = spectrum()->lightestHadronPair(pQ1->dataPtr(),pQ->dataPtr());
	auto LHP2 = spectrum()->lightestHadronPair(pQ2->dataPtr(),pQ->dataPtr());
	if (sqr(Mc1)<sqr(LHP1.first->mass()+LHP1.second->mass()))
	return true;
	if (sqr(Mc2)<sqr(LHP2.first->mass()+LHP2.second->mass()))
	return true;
	// double weigthHadrons = Kinematics::pstarTwoBodyDecay(Mc1,m1,m)Kinematics::pstarTwoBodyDecay(Mc2,m2,m)/(Kinematics::pstarTwoBodyDecay(Mc1,LHP1.first->mass(),LHP1.second->mass())Kinematics::pstarTwoBodyDecay(Mc2,LHP2.first->mass(),LHP2.second->mass()));
	// double tot= weightFlatPhaseSpace(Mc, Mc1, Mc2, m, m1, m2)/weigthHadrons;
	// if (std::isinf(tot) \|\| std::isnan(tot) \|\| tot<0.0 \|\| tot>1.0)
	// std::cout << "tot = " <<
	double lam1 = sqrt(Kinematics::kaellen(Mc1/GeV, LHP1.first->mass()/GeV, LHP1.second->mass()/GeV));
	double lam2 = sqrt(Kinematics::kaellen(Mc2/GeV, LHP2.first->mass()/GeV, LHP2.second->mass()/GeV));
	double lam3 = sqrt(Kinematics::kaellen(Mc/GeV, Mc1/GeV, Mc2/GeV));
	double ratio;
	// old weight of Jan without the Jacobi factor M1*M2 of the Mass integration
	// double PSweight = pow(lam1lam2lam3,_dim-3.)pow(M1_tempM2_temp,2.-_dim);
	// new weight with the Jacobi factor M1*M2 of the Mass integration
	double PSweight = pow(lam1lam2lam3,_dim-3.)pow(Mc1Mc2/GeV2,3.-_dim);
	// overestimate only possible for dim>=3.0
	assert(_dim>=3.0);
	// old overestimate of Jan without the Jacobi factor M1*M2 of the Mass integration
	// double overEstimate = _dim>=4.0 ? pow(M_temp,4._dim-14.):pow(M_temp,2(_dim-3.0))/pow((m1_temp+m_temp)*(m2_temp+m_temp),4.0-_dim);
	// new improved overestimate with the Jacobi factor M1*M2 of the Mass integration
	double overEstimate = pow(6.0sqrt(3.0), 3.0 - _dim)pow(Mc/GeV, 4.*_dim-12.);
	ratio = PSweight/overEstimate;
	// if (!(ratio>=0)) std::cout << "M "<<Mc/GeV<<" M1 "<<Mc1/GeV<<" M2 "<<Mc2/GeV<<" m1 "<<m1_temp<<" m2 "<<m2_temp<<" m "<<m_temp<<"\n\n";
	// if (ratio > 0.9) std::cout << "ratio = " << ratio <<"\n";
	if (std::isinf(ratio) \|\| std::isnan(ratio) \|\| ratio<0.0 \|\| ratio>1.0)
	std::cout << "ratio = " << ratio<<std::endl;
	assert (ratio >= 0);
	assert (ratio <= 1);
	return ratio;
	}

	/**
	* Veto for the phase space weight
	* returns true if proposed Masses are rejected
	* else returns false
	*/
	bool ClusterFissioner::phaseSpaceVeto(const Energy Mc, const Energy Mc1, const Energy Mc2,
	- const Energy m, const Energy m1, const Energy m2, tcPPtr pQ1, tcPPtr pQ2, tcPPtr pQ) const {
	+ const Energy m, const Energy m1, const Energy m2, tcPPtr pQ1, tcPPtr pQ2, tcPPtr pQ, const double power) const {
	switch (_phaseSpaceWeights)
	{
	case 1:
	- return phaseSpaceVetoConstituentMasses(Mc, Mc1, Mc2, m, m1, m2);
	+ return phaseSpaceVetoConstituentMasses(Mc, Mc1, Mc2, m, m1, m2, power);
	case 2:
	return phaseSpaceVetoHadronPairs(Mc, Mc1, Mc2, pQ, pQ1, pQ2);
	case 3:
	return phaseSpaceVetoNoConstituentMasses(Mc, Mc1, Mc2);
	default:
	assert(false);
	}
	}

	/**
	* Veto for the phase space weight
	* returns true if proposed Masses are rejected
	* else returns false
	*/
	bool ClusterFissioner::phaseSpaceVetoConstituentMasses(const Energy Mc, const Energy Mc1, const Energy Mc2,
	- const Energy m, const Energy m1, const Energy m2) const {
	- return (UseRandom::rnd()>weightFlatPhaseSpaceConstituentMasses(Mc, Mc1, Mc2, m, m1, m2));
	+ const Energy m, const Energy m1, const Energy m2, const double power) const {
	+ return (UseRandom::rnd()>weightPhaseSpaceConstituentMasses(Mc, Mc1, Mc2, m, m1, m2, power));
	}
	bool ClusterFissioner::phaseSpaceVetoNoConstituentMasses(const Energy Mc, const Energy Mc1, const Energy Mc2) const {
	return (UseRandom::rnd()>weightFlatPhaseSpaceNoConstituentMasses(Mc, Mc1, Mc2));
	}

	bool ClusterFissioner::phaseSpaceVetoHadronPairs(const Energy Mc, const Energy Mc1, const Energy Mc2, tcPPtr pQ, tcPPtr pQ1, tcPPtr pQ2) const {
	return (UseRandom::rnd()>weightFlatPhaseSpaceHadronMasses(Mc, Mc1, Mc2, pQ, pQ1, pQ2));
	}
	/**
	* Calculate the masses and possibly kinematics of the cluster
	* fission at hand; if claculateKineamtics is perfomring non-trivial
	* steps kinematics calulcated here will be overriden. Currently resorts to the default
	* @return the potentially non-trivial distribution weight=f(M1,M2)
	* On Failure we return 0
	*/
	double ClusterFissioner::drawNewMasses(const Energy Mc, const bool soft1, const bool soft2,
	Lorentz5Momentum& pClu1, Lorentz5Momentum& pClu2,
	tcPPtr ptrQ1, const Lorentz5Momentum& pQ1,
	tcPPtr newPtr1, const Lorentz5Momentum& pQone,
	tcPPtr, const Lorentz5Momentum& pQtwo,
	tcPPtr ptrQ2, const Lorentz5Momentum& pQ2) const {
	// TODO add precise weightMS that could be used used for improving the rejection Sampling
	switch (_massSampler)
	{
	case 0:
	return drawNewMassesDefault(Mc, soft1, soft2, pClu1, pClu2, ptrQ1, pQ1, tcPPtr(), pQone, tcPPtr(),pQtwo, ptrQ2, pQ2);
	break;
	case 1:
	return drawNewMassesUniform(Mc, pClu1, pClu2, pQ1, pQone, pQ2);
	break;
	case 2:
	return drawNewMassesPhaseSpace(Mc, pClu1, pClu2, pQ1, pQone, pQ2, ptrQ1, newPtr1, ptrQ2);
	break;
	+ case 3:
	+ return drawNewMassesPhaseSpacePowerLaw(Mc, pClu1, pClu2, pQ1, pQone, pQ2, ptrQ1, newPtr1, ptrQ2);
	+ break;
	default:
	assert(false);
	}
	return 0;// failure
	}

	/**
	* Calculate the masses and possibly kinematics of the cluster
	* fission at hand; if claculateKineamtics is perfomring non-trivial
	* steps kinematics claulcated here will be overriden. Currentl;y resorts to the default
	*/
	double ClusterFissioner::drawNewMassesDefault(const Energy Mc, const bool soft1, const bool soft2,
	Lorentz5Momentum& pClu1, Lorentz5Momentum& pClu2,
	tcPPtr ptrQ1, const Lorentz5Momentum& pQ1,
	tcPPtr, const Lorentz5Momentum& pQone,
	tcPPtr, const Lorentz5Momentum& pQtwo,
	tcPPtr ptrQ2, const Lorentz5Momentum& pQ2) const {
	// power for splitting
	double exp1 = !spectrum()->isExotic(ptrQ1->dataPtr()) ? _pSplitLight : _pSplitExotic;
	double exp2 = !spectrum()->isExotic(ptrQ2->dataPtr()) ? _pSplitLight : _pSplitExotic;
	for ( const long& id : spectrum()->heavyHadronizingQuarks() ) {
	assert(_pSplitHeavy.find(id) != _pSplitHeavy.end());
	if ( spectrum()->hasHeavy(id,ptrQ1->dataPtr()) ) exp1 = _pSplitHeavy.find(id)->second;
	if ( spectrum()->hasHeavy(id,ptrQ2->dataPtr()) ) exp2 = _pSplitHeavy.find(id)->second;
	}

	Energy M1 = drawChildMass(Mc,pQ1.mass(),pQ2.mass(),pQone.mass(),exp1,soft1);
	Energy M2 = drawChildMass(Mc,pQ2.mass(),pQ1.mass(),pQtwo.mass(),exp2,soft2);

	pClu1.setMass(M1);
	pClu2.setMass(M2);

	return 1.0; // succeeds
	}

	/**
	* Sample the masses for flat phase space
	* */
	double ClusterFissioner::drawNewMassesUniform(const Energy Mc, Lorentz5Momentum& pClu1, Lorentz5Momentum& pClu2,
	const Lorentz5Momentum& pQ1,
	const Lorentz5Momentum& pQ,
	const Lorentz5Momentum& pQ2) const {

	Energy M1,M2;
	const Energy m1 = pQ1.mass();
	const Energy m2 = pQ2.mass();
	const Energy m = pQ.mass();
	const Energy M1min = m1 + m;
	const Energy M2min = m2 + m;
	const Energy M1max = Mc - M2min;
	const Energy M2max = Mc - M1min;

	assert(M1max-M1min>ZERO);
	assert(M2max-M2min>ZERO);

	double r1;
	double r2;

	int counter = 0;
	const int max_counter = 100;

	while (counter < max_counter) {
	r1 = UseRandom::rnd();
	r2 = UseRandom::rnd();

	M1 = (M1max-M1min)*r1 + M1min;
	M2 = (M2max-M2min)*r2 + M2min;

	counter++;
	if ( Mc > M1 + M2) break;
	}

	if (counter==max_counter
	\|\| Mc < M1 + M2
	\|\| M1 <= M1min
	\|\| M2 <= M2min ) return 0.0; // failure

	pClu1.setMass(M1);
	pClu2.setMass(M2);

	return 1.0; // succeeds
	}
	+/**
	+ * Sample the masses for flat phase space
	+ * */
	+double ClusterFissioner::drawNewMassesPhaseSpacePowerLaw(const Energy Mc,
	+ Lorentz5Momentum& pClu1, Lorentz5Momentum& pClu2,
	+ const Lorentz5Momentum& pQ1,
	+ const Lorentz5Momentum& pQ,
	+ const Lorentz5Momentum& pQ2,
	+ tcPPtr ptrQ1, tcPPtr ptrQ2, tcPPtr ptrQ) const {
	+ Energy M1,M2,MuS;
	+ const Energy m1 = pQ1.mass();
	+ const Energy m2 = pQ2.mass();
	+ const Energy m = pQ.mass();
	+ const Energy M1min = m1 + m;
	+ const Energy M2min = m2 + m;
	+ const Energy M1max = Mc - M2min;
	+ const Energy M2max = Mc - M1min;
	+
	+ assert(M1max-M1min>ZERO);
	+ assert(M2max-M2min>ZERO);
	+
	+ double r1;
	+ double r2;
	+
	+ int counter = 0;
	+ const int max_counter = 200;
	+
	+ while (counter < max_counter) {
	+ r1 = UseRandom::rnd();
	+ r2 = UseRandom::rnd();
	+
	+ M1 = (M1max-M1min)*r1 + M1min;
	+ M2 = (M2max-M2min)*r2 + M2min;
	+
	+ counter++;
	+ if (sqr(M1+M2)>sqr(Mc))
	+ continue;
	+
	+ // if (!phaseSpaceVeto(Mc,M1,M2,m,m1,m2) ) break; // For FlatPhaseSpace sampling vetoing
	+ if (!phaseSpaceVeto(Mc,M1,M2,m,m1,m2, ptrQ1, ptrQ2, ptrQ,_powerLawPower) ) break; // For FlatPhaseSpace sampling vetoing
	+ }
	+ if (counter==max_counter) return 0.0; // failure
	+
	+ pClu1.setMass(M1);
	+ pClu2.setMass(M2);
	+
	+ return weightPhaseSpaceConstituentMasses(Mc, M1, M2, m, m1, m2, _powerLawPower); // succeeds return weight
	+}
	+

	/**
	* Sample the masses for flat phase space
	* */
	double ClusterFissioner::drawNewMassesPhaseSpace(const Energy Mc,
	Lorentz5Momentum& pClu1, Lorentz5Momentum& pClu2,
	const Lorentz5Momentum& pQ1,
	const Lorentz5Momentum& pQ,
	const Lorentz5Momentum& pQ2,
	tcPPtr ptrQ1, tcPPtr ptrQ2, tcPPtr ptrQ ) const {
	Energy M1,M2,MuS;
	const Energy m1 = pQ1.mass();
	const Energy m2 = pQ2.mass();
	const Energy m = pQ.mass();
	const Energy M1min = m1 + m;
	const Energy M2min = m2 + m;
	const Energy M1max = Mc - M2min;
	const Energy M2max = Mc - M1min;

	assert(M1max-M1min>ZERO);
	assert(M2max-M2min>ZERO);

	double r1;
	double r2;

	int counter = 0;
	const int max_counter = 200;

	while (counter < max_counter) {
	r1 = UseRandom::rnd();
	r2 = UseRandom::rnd();

	M1 = (M1max-M1min)*r1 + M1min;
	M2 = (M2max-M2min)*r2 + M2min;

	counter++;
	if (sqr(M1+M2)>sqr(Mc))
	continue;

	// if (!phaseSpaceVeto(Mc,M1,M2,m,m1,m2) ) break; // For FlatPhaseSpace sampling vetoing
	if (!phaseSpaceVeto(Mc,M1,M2,m,m1,m2, ptrQ1, ptrQ2, ptrQ) ) break; // For FlatPhaseSpace sampling vetoing
	}
	if (counter==max_counter) return 0.0; // failure

	pClu1.setMass(M1);
	pClu2.setMass(M2);

	return weightFlatPhaseSpace(Mc, M1, M2, m, m1, m2, ptrQ1, ptrQ2, ptrQ); // succeeds return weight
	}

	void ClusterFissioner::drawNewFlavourDiquarks(PPtr& newPtrPos,PPtr& newPtrNeg,
	const ClusterPtr & clu) const {

	// Flavour is assumed to be only u, d, s, with weights
	// (which are not normalized probabilities) given
	// by the same weights as used in HadronsSelector for
	// the decay of clusters into two hadrons.

	unsigned hasDiquarks=0;
	assert(clu->numComponents()==2);
	tcPDPtr pD1=clu->particle(0)->dataPtr();
	tcPDPtr pD2=clu->particle(1)->dataPtr();
	bool isDiq1=DiquarkMatcher::Check(pD1->id());
	if (isDiq1) hasDiquarks++;
	bool isDiq2=DiquarkMatcher::Check(pD2->id());
	if (isDiq2) hasDiquarks++;
	assert(hasDiquarks<=2);
	Energy Mc=(clu->momentum().mass());
	// if (fabs(clu->momentum().massError() )>1e-14) std::cout << "Mass inconsistency CF : " << std::scientific << clu->momentum().massError() <<"\n";
	// Not allow yet Diquark Clusters
	// if ( hasDiquarks>=1 \|\| Mc < spectrum()->massLightestBaryonPair(pD1,pD2) )
	// return drawNewFlavour(newPtrPos,newPtrNeg);

	Energy minMass;
	double weight;
	// double factorPS;
	Selector<long> choice;
	// int countQ=0;
	// int countDiQ=0;
	// adding quark-antiquark pairs to the selection list
	for ( const long& id : spectrum()->lightHadronizingQuarks() ) {
	// TODO uncommenting below gives sometimes 0 selection possibility,
	// maybe need to be checked in the LightClusterDecayer and ColourReconnector
	// if (Mc < spectrum()->massLightestHadronPair(pD1,pD2)) continue;
	// countQ++;
	minMass=spectrum()->massLightestHadronPair(pD1,pD2);
	if (_fissionCluster==0) choice.insert(_hadronSpectrum->pwtQuark(id),id);
	else if (abs(_fissionCluster)==1) choice.insert(_fissionPwt.find(id)->second,id);
	else assert(false);
	}
	// adding diquark-antidiquark pairs to the selection list
	switch (hasDiquarks)
	{
	case 0:
	for ( const long& id : spectrum()->lightHadronizingDiquarks() ) {
	if (_strictDiquarkKinematics) {
	tPDPtr cand = getParticleData(id);
	Energy mH1=spectrum()->massLightestHadron(pD2,cand);
	Energy mH2=spectrum()->massLightestHadron(cand,pD1);
	// factorPS = Kinematics::pstarTwoBodyDecay(Mc,mH1,mH2)/(Mc/2.0);
	// factorPS = 1.0;
	minMass = mH1 + mH2;
	}
	else {
	minMass = spectrum()->massLightestBaryonPair(pD1,pD2);
	// factorPS = 1.0;
	}
	if (Mc < minMass) continue;
	// countDiQ++;
	if (_fissionCluster==0) weight = _hadronSpectrum->pwtQuark(id);
	else if (abs(_fissionCluster)==1) weight = _fissionPwt.find(id)->second;
	else assert(false);
	if (_fissionCluster==-1)
	weight=sqr(Herwig::Math::alphaS(Mc, 0.25GeV,3, 2));
	// weight/=factorPS;
	choice.insert(weight,id);
	}
	break;
	case 1:
	if (_diquarkClusterFission<1) break;
	for ( const long& id : spectrum()->lightHadronizingDiquarks() ) {
	tPDPtr diq = getParticleData(id);
	if (isDiq1)
	minMass = spectrum()->massLightestHadron(pD2,diq)
	+ spectrum()->massLightestBaryonPair(diq,pD1);
	else
	minMass = spectrum()->massLightestHadron(pD1,diq)
	+ spectrum()->massLightestBaryonPair(diq,pD2);
	if (Mc < minMass) continue;
	// countDiQ++;
	if (_fissionCluster==0) weight = _hadronSpectrum->pwtQuark(id);
	else if (abs(_fissionCluster)==1) weight = _fissionPwt.find(id)->second;
	else assert(false);
	if (_fissionCluster==-1)
	weight=sqr(Herwig::Math::alphaS(Mc, 0.25GeV,3, 2));
	choice.insert(weight,id);
	}
	break;
	case 2:
	if (_diquarkClusterFission<2) break;
	for ( const long& id : spectrum()->lightHadronizingDiquarks() ) {
	tPDPtr diq = getParticleData(id);
	if (Mc < spectrum()->massLightestBaryonPair(pD1,pD2)) {
	throw Exception() << "Found Diquark Cluster:\n" << *clu << "\nwith MassCluster = "
	<< ounit(Mc,GeV) <<" GeV MassLightestBaryonPair = "
	<< ounit(spectrum()->massLightestBaryonPair(pD1,pD2) ,GeV)
	<< " GeV cannot decay" << Exception::eventerror;
	}
	minMass = spectrum()->massLightestBaryonPair(pD1,diq)
	+ spectrum()->massLightestBaryonPair(diq,pD2);
	if (Mc < minMass) continue;
	// countDiQ++;
	if (_fissionCluster==0) weight = _hadronSpectrum->pwtQuark(id);
	else if (abs(_fissionCluster)==1) weight = _fissionPwt.find(id)->second;
	else assert(false);
	if (_fissionCluster==-1)
	weight=sqr(Herwig::Math::alphaS(Mc, 0.25GeV,3, 2));
	choice.insert(weight,id);
	}
	break;
	default:
	assert(false);
	}
	assert(choice.size()>0);
	long idNew = choice.select(UseRandom::rnd());
	newPtrPos = getParticle(idNew);
	newPtrNeg = getParticle(-idNew);
	assert(newPtrPos);
	assert(newPtrNeg);
	assert(newPtrPos->dataPtr());
	assert(newPtrNeg->dataPtr());

	}

	void ClusterFissioner::drawNewFlavourQuarks(PPtr& newPtrPos,PPtr& newPtrNeg) const {

	// Flavour is assumed to be only u, d, s, with weights
	// (which are not normalized probabilities) given
	// by the same weights as used in HadronsSelector for
	// the decay of clusters into two hadrons.


	Selector<long> choice;
	switch(abs(_fissionCluster)){
	case 0:
	for ( const long& id : spectrum()->lightHadronizingQuarks() )
	choice.insert(_hadronSpectrum->pwtQuark(id),id);
	break;
	case 1:
	for ( const long& id : spectrum()->lightHadronizingQuarks() )
	choice.insert(_fissionPwt.find(id)->second,id);
	break;
	default :
	assert(false);
	}
	long idNew = choice.select(UseRandom::rnd());
	newPtrPos = getParticle(idNew);
	newPtrNeg = getParticle(-idNew);
	assert (newPtrPos);
	assert(newPtrNeg);
	assert (newPtrPos->dataPtr());
	assert(newPtrNeg->dataPtr());

	}

	void ClusterFissioner::drawNewFlavourEnhanced(PPtr& newPtrPos,PPtr& newPtrNeg,
	Energy2 mass2) const {

	if ( spectrum()->gluonId() != ParticleID::g )
	throw Exception() << "strange enhancement only working with Standard Model hadronization"
	<< Exception::runerror;

	// Flavour is assumed to be only u, d, s, with weights
	// (which are not normalized probabilities) given
	// by the same weights as used in HadronsSelector for
	// the decay of clusters into two hadrons.

	double prob_d = 0.;
	double prob_u = 0.;
	double prob_s = 0.;
	double scale = abs(double(sqr(_m0Fission)/mass2));
	// Choose which splitting weights you wish to use
	switch(abs(_fissionCluster)){
	// 0: ClusterFissioner and ClusterDecayer use the same weights
	case 0:
	prob_d = _hadronSpectrum->pwtQuark(ParticleID::d);
	prob_u = _hadronSpectrum->pwtQuark(ParticleID::u);
	/* Strangeness enhancement:
	Case 1: probability scaling
	Case 2: Exponential scaling
	*/
	if (_enhanceSProb == 1)
	prob_s = (_maxScale < scale) ? 0. : pow(_hadronSpectrum->pwtQuark(ParticleID::s),scale);
	else if (_enhanceSProb == 2)
	prob_s = (_maxScale < scale) ? 0. : exp(-scale);
	break;
	/* 1: ClusterFissioner uses its own unique set of weights,
	i.e. decoupled from ClusterDecayer */
	case 1:
	prob_d = _fissionPwt.find(ParticleID::d)->second;
	prob_u = _fissionPwt.find(ParticleID::u)->second;
	if (_enhanceSProb == 1)
	prob_s = (_maxScale < scale) ? 0. : pow(_fissionPwt.find(ParticleID::s)->second,scale);
	else if (_enhanceSProb == 2)
	prob_s = (_maxScale < scale) ? 0. : exp(-scale);
	break;
	default:
	assert(false);
	}

	int choice = UseRandom::rnd3(prob_u, prob_d, prob_s);
	long idNew = 0;
	switch (choice) {
	case 0: idNew = ThePEG::ParticleID::u; break;
	case 1: idNew = ThePEG::ParticleID::d; break;
	case 2: idNew = ThePEG::ParticleID::s; break;
	}
	newPtrPos = getParticle(idNew);
	newPtrNeg = getParticle(-idNew);
	assert (newPtrPos);
	assert(newPtrNeg);
	assert (newPtrPos->dataPtr());
	assert(newPtrNeg->dataPtr());

	}

	Energy2 ClusterFissioner::clustermass(const ClusterPtr & cluster) const {
	Lorentz5Momentum pIn = cluster->momentum();
	Energy2 endpointmass2 = sqr(cluster->particle(0)->mass() +
	cluster->particle(1)->mass());
	Energy2 singletm2 = pIn.m2();
	// Return either the cluster mass, or the lambda measure
	return (_massMeasure == 0) ? singletm2 : singletm2 - endpointmass2;
	}

	Energy ClusterFissioner::drawChildMass(const Energy M, const Energy m1,
	const Energy m2, const Energy m,
	const double expt, const bool soft) const {

	/***************************
	* This method, given in input the cluster mass Mclu of an heavy cluster C,
	* made of consituents of masses m1 and m2, draws the masses Mclu1 and Mclu2
	* of, respectively, the children cluster C1, made of constituent masses m1
	* and m, and cluster C2, of mass Mclu2 and made of constituent masses m2
	* and m. The mass is extracted from one of the two following mass
	* distributions:
	* --- power-like ("normal" distribution)
	* d(Prob) / d(M^exponent) = const
	* where the exponent can be different from the two children C1 (exp1)
	* and C2 (exponent2).
	* --- exponential ("soft" distribution)
	* d(Prob) / d(M^2) = exp(-b*M)
	* where b = 2.0 / average.
	* Such distributions are limited below by the masses of
	* the constituents quarks, and above from the mass of decaying cluster C.
	* The choice of which of the two mass distributions to use for each of the
	* two cluster children is dictated by iRemnant (see below).
	* If the number of attempts to extract a pair of mass values that are
	* kinematically acceptable is above some fixed number (max_loop, see below)
	* the method gives up and returns false; otherwise, when it succeeds, it
	* returns true.
	*
	* These distributions have been modified from HERWIG:
	* Before these were:
	* Mclu1 = m1 + (Mclu - m1 - m2)*pow( rnd(), 1.0/exponent1 );
	* The new one coded here is a more efficient version, same density
	* but taking into account 'in phase space from' beforehand
	***************************/
	// hard cluster
	if(!soft) {
	return pow(UseRandom::rnd(pow((M-m1-m2-m)*UnitRemoval::InvE, expt),
	pow(m*UnitRemoval::InvE, expt)), 1./expt
	)*UnitRemoval::E + m1;
	}
	// Otherwise it uses a soft mass distribution
	else {
	static const InvEnergy b = 2.0 / _btClM;
	Energy max = M-m1-m2-2.0*m;
	double rmin = b*max;
	rmin = ( rmin < 50 ) ? exp(-rmin) : 0.;
	double r1;
	do {
	r1 = UseRandom::rnd(rmin, 1.0) * UseRandom::rnd(rmin, 1.0);
	}
	while (r1 < rmin);
	return m1 + m - log(r1)/b;
	}
	}
	std::pair<Axis,double> ClusterFissioner::sampleDirectionConstituents(
	const Lorentz5Momentum & pClu, const Energy Mcluster) const {
	switch (_phaseSpaceSamplerConstituent)
	{
	case 0:
	{
	// Aligned
	return std::make_pair(pClu.vect().unit(),1.0);
	}
	case 1:
	{
	// Isotropic
	return std::make_pair(sampleDirectionIsotropic(),1.0);
	}
	case 2:
	{
	// Exponential
	// New
	// double C_est_Exponential=1.18;
	// double power_est_Exponential=0.65;
	// Old
	double C_est_Exponential=0.63;
	double power_est_Exponential=0.89;
	double lambda=C_est_Exponential*pow(Mcluster/GeV,power_est_Exponential);
	return sampleDirectionExponential(pClu.vect().unit(), lambda);
	}
	default:
	assert(false);
	break;
	}
	return std::make_pair(Axis(),0.0); // Failure
	}

	std::pair<Axis,double> ClusterFissioner::sampleDirectionCluster(
	const Lorentz5Momentum & pQ,
	const Lorentz5Momentum & pClu) const {
	switch (_phaseSpaceSamplerCluster)
	{
	case 0:
	{
	// Aligned
	Axis dir;
	if (_covariantBoost)
	// in Covariant Boost the positive z-Axis is defined as the direction of
	// the pQ vector in the Cluster rest frame
	dir = Axis(0,0,1);
	else
	dir = sampleDirectionAligned(pQ, pClu);
	return std::make_pair(dir,1.0);
	}
	case 1:
	{
	// Isotropic
	return std::make_pair(sampleDirectionIsotropic(),1.0);
	}
	case 2:
	{
	// Tchannel
	Energy M=pClu.mass();
	// New
	// double C_est_Tchannel=3.66;
	// double power_est_Tchannel=-1.95;
	// Old
	double C_est_Tchannel=9.78;
	double power_est_Tchannel=-2.5;
	double A = C_est_Tchannel*pow(M/GeV,power_est_Tchannel);
	double Ainv = 1.0/(1.0+A);
	return sampleDirectionTchannel(sampleDirectionAligned(pQ,pClu),Ainv);
	}
	default:
	assert(false);
	}
	return std::make_pair(Axis(),0.0); // Failure
	}

	std::pair<Axis,double> ClusterFissioner::sampleDirectionExponential(const Axis dirQ, const double lambda) const {
	Axis FinalDir = dirQ;
	double cosTheta = Kinematics::sampleCosExp(lambda);
	double phi;
	// std::cout << "cosThetaSamp = "<< cosTheta <<"\tlambda = " <<lambda << std::endl;
	// If no change in angle keep the direction fixed
	if (fabs(cosTheta-1.0)>1e-14) {
	// rotate to sampled angles
	FinalDir.rotate(acos(cosTheta),dirQ.orthogonal());
	phi = UseRandom::rnd(-Constants::pi,Constants::pi);
	FinalDir.rotate(phi,dirQ);
	}
	// std::cout << "cosThetaTrue = "<< FinalDir.cosTheta(dirQ) <<"\tlambda = " <<lambda << std::endl;
	double weight = exp(lambda*(cosTheta-1.0));
	if (std::isnan(weight) \|\|std::isinf(weight) )
	std::cout << "lambda = " << lambda << "\t cos " << cosTheta<< std::endl;
	assert(!std::isnan(weight));
	assert(!std::isinf(weight));
	if (fabs(FinalDir.cosTheta(dirQ)-cosTheta)>1e-8) std::cout << "cosThetaTrue = "<< FinalDir.cosTheta(dirQ) <<"\tlambda = " <<lambda << std::endl;
	return std::make_pair(FinalDir,weight);
	}
	std::pair<Axis,double> ClusterFissioner::sampleDirectionTchannel(const Axis dirQ, const double Ainv) const {
	double cosTheta = Kinematics::sampleCosTchannel(Ainv);
	double phi;
	Axis FinalDir = dirQ;
	// If no change in angle keep the direction fixed
	if (fabs(cosTheta-1.0)>1e-14) {
	// rotate to sampled angles
	FinalDir.rotate(acos(cosTheta),dirQ.orthogonal());
	phi = UseRandom::rnd(-Constants::pi,Constants::pi);
	FinalDir.rotate(phi,dirQ);
	}
	double weight = 0.0;
	if (1.0-Ainv>1e-10)
	weight=1.0/pow(1.0/Ainv-cosTheta,2.0);
	else
	weight=1.0/pow((1.0-cosTheta)+(1.0-Ainv),2.0);
	if (std::isnan(weight) \|\|std::isinf(weight) )
	std::cout << "Ainv = " << Ainv << "\t cos " << cosTheta<< std::endl;
	assert(!std::isnan(weight));
	assert(!std::isinf(weight));
	if (fabs(FinalDir.cosTheta(dirQ)-cosTheta)>1e-8) std::cout << "cosThetaTrue = "<< FinalDir.cosTheta(dirQ) <<"\tAinv = " <<Ainv << std::endl;
	return std::make_pair(FinalDir,weight);
	}

	Axis ClusterFissioner::sampleDirectionAligned(const Lorentz5Momentum & pQ, const Lorentz5Momentum & pClu) const {
	Lorentz5Momentum pQinCOM(pQ);
	pQinCOM.setMass(pQ.m());
	pQinCOM.boost( -pClu.boostVector() ); // boost from LAB to C
	return pQinCOM.vect().unit();
	}

	Axis ClusterFissioner::sampleDirectionAlignedSmeared(const Lorentz5Momentum & pQ, const Lorentz5Momentum & pClu) const {
	Lorentz5Momentum pQinCOM(pQ);
	pQinCOM.setMass(pQ.m());
	pQinCOM.boost( -pClu.boostVector() ); // boost from LAB to C
	if (pQinCOM.vect().mag2()<=ZERO){
	return sampleDirectionAlignedSmeared(pClu-pQ,pClu);
	}
	const Axis dir = pQinCOM.vect().unit();
	double cluSmear = 0.5;
	double cosTheta;
	do {
	cosTheta = 1.0 + cluSmear*log( UseRandom::rnd() );
	}
	while (fabs(cosTheta)>1.0 \|\| std::isnan(cosTheta) \|\| std::isinf(cosTheta));

	if (!(dir.mag2()>ZERO)){
	std::cout << "\nDRI = 0"<< dir.mag2() << std::endl;
	}
	Axis dirSmeared = dir;
	double phi;
	if (fabs(cosTheta-1.0)>1e-10) {
	// rotate to sampled angles
	dirSmeared.rotate(acos(cosTheta),dir.orthogonal());
	phi = UseRandom::rnd(-M_PI,M_PI);
	dirSmeared.rotate(phi,dir);
	}
	if (!(dirSmeared.mag2()>ZERO)){
	std::cout << "\nDRI SMR = 0" <<cosTheta<< "\t" <<acos(cosTheta)<< "\t"<< phi<< "\t"<< dirSmeared.mag2() << std::endl;
	std::cout << dir.orthogonal() << std::endl;
	}
	return dirSmeared;
	}

	Axis ClusterFissioner::sampleDirectionIsotropic() const {
	double cosTheta = -1 + 2.0 * UseRandom::rnd();
	double sinTheta = sqrt(1.0-cosTheta*cosTheta);
	double Phi = 2.0 * Constants::pi * UseRandom::rnd();
	Axis uClusterUniform(cos(Phi)sinTheta, sin(Phi)sinTheta, cosTheta);
	return uClusterUniform.unit();
	}

	Axis ClusterFissioner::sampleDirectionSemiUniform(const Lorentz5Momentum & pQ, const Lorentz5Momentum & pClu) const {
	Axis dir = sampleDirectionAligned(pQ,pClu);
	Axis res;
	do {
	res=sampleDirectionIsotropic();
	}
	while (dir*res<0);
	return res;
	}
	namespace {
	double SoftFunction(
	const Lorentz5Momentum & pi,
	const Lorentz5Momentum & pj,
	const Lorentz5Momentum & q,
	const Lorentz5Momentum & qbar) {
	Energy2 mq2 = q*q;
	double numerator = -(pipj)(q*qbar+mq2)/sqr(GeV2);
	double denominator = sqr(qqbar+mq2)(pi(q+qbar))(pj(q+qbar))/sqr(GeV2GeV2);
	int cataniScheme = 1;
	// Different expressions which should yield similar results
	switch (cataniScheme) {
	case 0:
	numerator += ((piq)(pjqbar) + (pjq)(piqbar))/sqr(GeV2);
	break;
	case 1:
	numerator += - (0.5(pi(q-qbar))(pj(q-qbar)))/sqr(GeV2);
	break;
	default:
	assert(false);
	}
	double Iij = numerator/denominator;
	return Iij;
	}
	}
	/* SQME for p1,p2->C1(q1,q),C2(q2,qbar)
	* Note that:
	* p0Q1 -> p1
	* p0Q2 -> p2
	* pQ1 -> q1
	* pQone -> q
	* pQ2 -> q2
	* pQone -> qbar
	* */
	double ClusterFissioner::calculateSQME(
	const Lorentz5Momentum & p1,
	const Lorentz5Momentum & p2,
	const Lorentz5Momentum & q1,
	const Lorentz5Momentum & q,
	const Lorentz5Momentum & q2,
	const Lorentz5Momentum & qbar) const {
	double SQME;
	switch (_matrixElement)
	{
	case 0:
	SQME = 1.0;
	break;
	case 1:
	{
	/*
	// Energy2 p1p2 = p1*p2;
	Energy2 q1q2 = q1 * q2;
	Energy2 q1q = q1 * q ;
	Energy2 q2qbar = q2 * qbar;
	Energy2 q2q = q2 * q ;
	Energy2 q1qbar = q1 * qbar;
	Energy2 qqbar = q * qbar;
	Energy2 mq2 = q.mass2();
	double Numerator = q1q2 * (qqbar + mq2)/sqr(GeV2);
	Numerator += 0.5 * (q1q - q1qbar)*(q2q - q2qbar)/sqr(GeV2);
	double Denominator = sqr(qqbar + mq2)(q1q + q1qbar)(q2q + q2qbar)/sqr(sqr(GeV2));
	SQME = Numerator/Denominator;
	*/
	double I11 = SoftFunction(q1,q1,q,qbar);
	double I22 = SoftFunction(q2,q2,q,qbar);
	double I12 = SoftFunction(q1,q2,q,qbar);
	SQME = (I11+I22-2.0*I12);
	if (SQME<0.0) {
	std::cout << "I11 = " << I11<< std::endl;
	std::cout << "I22 = " << I22<< std::endl;
	std::cout << "2*I12 = " << I12<< std::endl;
	std::cout << "soft = " << I11+I22-2.0*I12<< std::endl;
	std::cout << "M = " << (p1+p2).m()/GeV<< std::endl;
	std::cout << "M1 = " << (q1+q).m()/GeV<< std::endl;
	std::cout << "M2 = " << (q2+qbar).m()/GeV<< std::endl;
	std::cout << "m1 = " << (q1).m()/GeV<< std::endl;
	std::cout << "m2 = " << (q2).m()/GeV<< std::endl;
	std::cout << "m = " << (q).m()/GeV<< std::endl;
	}
	break;
	}
	case 2:
	{
	/*
	Energy2 p1p2 = p1 * p2;
	Energy2 p1q = p1 * q ;
	Energy2 p2qbar = p2 * qbar;
	Energy2 p2q = p2 * q ;
	Energy2 p1qbar = p1 * qbar;
	Energy2 qqbar = q * qbar;
	Energy2 mq2 = q.mass2();
	double Numerator = p1p2 * (qqbar + mq2)/sqr(GeV2);
	Numerator += 0.5 * (p1q - p1qbar)*(p2q - p2qbar)/sqr(GeV2);
	double Denominator = sqr(qqbar + mq2)(p1q + p1qbar)(p2q + p2qbar)/sqr(sqr(GeV2));
	SQME = Numerator/Denominator;
	*/
	double I11 = SoftFunction(p1,p1,q,qbar);
	double I22 = SoftFunction(p2,p2,q,qbar);
	double I12 = SoftFunction(p1,p2,q,qbar);
	SQME = (I11+I22-2.0*I12);
	if (SQME<0.0) {
	std::cout << "I11 = " << I11<< std::endl;
	std::cout << "I22 = " << I22<< std::endl;
	std::cout << "2*I12 = " << I12<< std::endl;
	std::cout << "soft = " << I11+I22-2.0*I12<< std::endl;
	std::cout << "M = " << (p1+p2).m()/GeV<< std::endl;
	std::cout << "M1 = " << (q1+q).m()/GeV<< std::endl;
	std::cout << "M2 = " << (q2+qbar).m()/GeV<< std::endl;
	std::cout << "m1 = " << (q1).m()/GeV<< std::endl;
	std::cout << "m2 = " << (q2).m()/GeV<< std::endl;
	std::cout << "m = " << (q).m()/GeV<< std::endl;
	}
	break;
	}
	case 4:
	{
	SQME = sqr(sqr(GeV2))/(sqr(qqbar-qq)(p1(q1+q)-p1p1-sqrt((p1p1)(qq)))(p2(q2+qbar)-p2p2-sqrt((p2p2)(qq))));
	break;
	}
	case 5:
	{
	/*
	Energy2 p1p2 = p1 * p2;
	Energy2 p1q = p1 * q ;
	Energy2 p2qbar = p2 * qbar;
	Energy2 p2q = p2 * q ;
	Energy2 p1qbar = p1 * qbar;
	Energy2 qqbar = q * qbar;
	Energy2 mq2 = q.mass2();
	double Numerator = p1p2 * (qqbar + mq2)/sqr(GeV2);
	Numerator += 0.5 * (p1q - p1qbar)*(p2q - p2qbar)/sqr(GeV2);
	double Denominator = sqr(qqbar + mq2)(p1q + p1qbar)(p2q + p2qbar)/sqr(sqr(GeV2));
	// add test Tchannel Matrix Element interference with gluon mass regulated
	// Denominator = sqr(sqr(p1-q1)-_epsilonIRsqrt((p1p1)(p2*p2)))/sqr(GeV2);
	*/
	Energy mg=getParticleData(spectrum()->gluonId())->constituentMass();
	double Denominator = (sqr(p1-q1)-_epsilonIRmgmg)/(GeV2);
	Denominator = (sqr(p2-q2)-_epsilonIRmg*mg)/(GeV2);
	double Numerator = ((q1q2)(p1p2)+(p2q1)(p1q2))/sqr(GeV2);
	// double I11 = SoftFunction(p1,p1,q,qbar);
	// double I22 = SoftFunction(p2,p2,q,qbar);
	// double I12 = SoftFunction(p1,p2,q,qbar);
	double I11 = SoftFunction(q1,q1,q,qbar);
	double I22 = SoftFunction(q2,q2,q,qbar);
	double I12 = SoftFunction(q1,q2,q,qbar);
	SQME = Numerator(I11+I22-2.0I12)/Denominator;
	if (SQME<0.0) {
	std::cout << "I11 = " << I11<< std::endl;
	std::cout << "I22 = " << I22<< std::endl;
	std::cout << "2*I12 = " << I12<< std::endl;
	std::cout << "soft = " << I11+I22-2.0*I12<< std::endl;
	std::cout << "Numerator = " << Numerator<< std::endl;
	std::cout << "Denominator = " << Denominator<< std::endl;
	std::cout << "M = " << (p1+p2).m()/GeV<< std::endl;
	std::cout << "M1 = " << (q1+q).m()/GeV<< std::endl;
	std::cout << "M2 = " << (q2+qbar).m()/GeV<< std::endl;
	std::cout << "m1 = " << (q1).m()/GeV<< std::endl;
	std::cout << "m2 = " << (q2).m()/GeV<< std::endl;
	std::cout << "m = " << (q).m()/GeV<< std::endl;
	}
	break;
	}
	default:
	assert(false);
	}
	if (SQME < 0) throw Exception()
	<< "Squared Matrix Element = "<< SQME <<" < 0 in ClusterFissioner::calculateSQME() "
	<< Exception::runerror;
	return SQME;
	}
	/* Overestimate for SQME for p1,p2->C1(q1,q),C2(q2,qbar)
	* Note that:
	* p0Q1 -> p1 where Mass -> m1
	* p0Q2 -> p2 where Mass -> m2
	* pQ1 -> q1 where Mass -> m1
	* pQone -> q where Mass -> mq
	* pQ2 -> q2 where Mass -> m2
	* pQone -> qbar where Mass -> mq
	* */
	double ClusterFissioner::calculateSQME_OverEstimate(
	const Energy& Mc,
	const Energy& m1,
	const Energy& m2,
	const Energy& mq
	) const {
	double SQME_OverEstimate;
	switch (_matrixElement)
	{
	case 0:
	SQME_OverEstimate = 1.0;
	break;
	case 1:
	{
	// Fit factor for guess of best overestimate
	double A = 0.25;
	SQME_OverEstimate = _safetyFactorMatrixElementApow(mq/GeV,-4);
	break;
	}
	case 2:
	{
	// Fit factor for guess of best overestimate
	double A = 0.25;
	SQME_OverEstimate = _safetyFactorMatrixElementApow(mq/GeV,-4);
	break;
	}
	case 4:
	{
	// Fit factor for guess of best overestimate
	SQME_OverEstimate = _safetyFactorMatrixElementsqr(sqr(GeV2))/(sqr(mqmq)m1m2(m1+mq)(m2+mq));
	break;
	}
	case 5:
	{
	// Fit factor for guess of best overestimate
	// Energy MassMax = 8.550986578849006037e+00*GeV; // mass max of python
	// double wTotMax = 5.280507222727160297e+03; // at MassMax
	// double argExp = 55.811; // from fit of python package
	// double powLog = 0.08788; // from fit of python package

	Energy MassMax = 5.498037126562503119e+01*GeV; // mass max of python
	double wTotMax = 1.570045806367627112e+06; // at MassMax
	double argExp = 16.12; // from fit of python package
	double powLog = 0.3122; // from fit of python package
	// Energy2 p1p2=0.5(McMc-m1m1-m2m2);
	// SQME_OverEstimate = _safetyFactorMatrixElementApow(mq/GeV,-4)(2sqr(p1p2))/sqr(_epsilonIR(m1m2));
	Energy Mmin = (m1+m2+2*mq);
	// Energy MdblPrec = Mmin*exp(pow(pow(log(MassMax/Mmin),powLog)+600.0/argExp,1.0/powLog));
	// if (Mc >MdblPrec)
	// std::cout << "errroRRRRR R MC = " << Mc/GeV << "\t Mcmax "<< MdblPrec/GeV<< std::endl;
	if (MassMax<Mmin) MassMax=Mmin;
	double Overestimate = wTotMaxexp(argExp(pow(log(Mc/Mmin),powLog)-pow(log(MassMax/Mmin),powLog)));
	SQME_OverEstimate = _safetyFactorMatrixElementOverestimate;//Apow(mq/GeV,-4)(2sqr(p1p2))/sqr(_epsilonIR(m1*m2));
	if (SQME_OverEstimate==0 \|\| std::isinf(SQME_OverEstimate)\|\| std::isnan(SQME_OverEstimate)) {
	std::cout << "SQME_OverEstimate is " << SQME_OverEstimate << std::endl;
	std::cout << "Overestimate is " << Overestimate << std::endl;
	std::cout << "MC " << Mc/GeV << std::endl;
	std::cout << "m " << mq/GeV << std::endl;
	std::cout << "m1 " << m1/GeV << std::endl;
	std::cout << "m2 " << m2/GeV << std::endl;
	}
	break;
	}
	default:
	assert(false);
	}
	return SQME_OverEstimate;
	}

	double ClusterFissioner::drawKinematics(
	const Lorentz5Momentum & pClu,
	const Lorentz5Momentum & p0Q1,
	const Lorentz5Momentum & p0Q2,
	const bool toHadron1,
	const bool toHadron2,
	Lorentz5Momentum & pClu1,
	Lorentz5Momentum & pClu2,
	Lorentz5Momentum & pQ1,
	Lorentz5Momentum & pQbar,
	Lorentz5Momentum & pQ,
	Lorentz5Momentum & pQ2bar) const {
	double weightTotal=0.0;
	if (pClu.mass() < pClu1.mass() + pClu2.mass()
	\|\| pClu1.mass()<ZERO
	\|\| pClu2.mass()<ZERO ) {
	throw Exception() << "Impossible Kinematics in ClusterFissioner::drawKinematics() (A)\n"
	<< "Mc = "<< pClu.mass()/GeV <<" GeV\n"
	<< "Mc1 = "<< pClu1.mass()/GeV <<" GeV\n"
	<< "Mc2 = "<< pClu2.mass()/GeV <<" GeV\n"
	<< Exception::eventerror;
	}

	// Sample direction of the daughter clusters
	std::pair<Axis,double> directionCluster = sampleDirectionCluster(p0Q1, pClu);
	Axis DirToClu = directionCluster.first;
	weightTotal = directionCluster.second;
	if (0 && _writeOut) {
	ofstream out("WriteOut/test_Cluster.dat",std::ios::app);
	Lorentz5Momentum pQinCOM(p0Q1);
	pQinCOM.setMass(p0Q1.m());
	pQinCOM.boost( -pClu.boostVector() ); // boost from LAB to C
	out << DirToClu.cosTheta(pQinCOM.vect()) << "\n";
	}
	if (_covariantBoost) {
	const Energy M = pClu.mass();
	const Energy M1 = pClu1.mass();
	const Energy M2 = pClu2.mass();
	const Energy PcomClu=Kinematics::pstarTwoBodyDecay(M,M1,M2);
	Momentum3 pClu1sampVect( PcomClu*DirToClu);
	Momentum3 pClu2sampVect(-PcomClu*DirToClu);
	pClu1.setMass(M1);
	pClu1.setVect(pClu1sampVect);
	pClu1.rescaleEnergy();
	pClu2.setMass(M2);
	pClu2.setVect(pClu2sampVect);
	pClu2.rescaleEnergy();
	}
	else {
	Lorentz5Momentum pClutmp(pClu);
	pClutmp.boost(-pClu.boostVector());
	Kinematics::twoBodyDecay(pClutmp, pClu1.mass(), pClu2.mass(),DirToClu, pClu1, pClu2);
	}
	// In the case that cluster1 does not decay immediately into a single hadron,
	// calculate the momenta of Q1 (as constituent of C1) and Qbar in the
	// (parent) C1 frame first, where the direction of Q1 is u.vect().unit(),
	// and then boost back in the LAB frame.
	if(!toHadron1) {
	if (pClu1.m() < pQ1.mass() + pQbar.mass() ) {
	throw Exception() << "Impossible Kinematics in ClusterFissioner::drawKinematics() (B)"
	<< Exception::eventerror;
	}
	std::pair<Axis,double> direction1 = sampleDirectionConstituents(pClu1,pClu.mass());
	// Need to boost constituents first into the pClu rest frame
	// boost from Cluster1 rest frame to Cluster COM Frame
	// Kinematics::twoBodyDecay(pClu1, pQ1.mass(), pQbar.mass(), DirClu1, pQ1, pQbar);
	Kinematics::twoBodyDecay(pClu1, pQ1.mass(), pQbar.mass(), direction1.first, pQ1, pQbar);
	weightTotal*=direction1.second;
	if (0 && _writeOut) {
	ofstream out("WriteOut/test_Constituents1.dat",std::ios::app);
	Lorentz5Momentum pQ1inCOM(pQ1);
	pQ1inCOM.setMass(pQ1.m());
	pQ1inCOM.boost( -pClu1.boostVector() ); // boost from LAB to C
	out << pQ1inCOM.vect().cosTheta(pClu1.vect()) << "\n";
	}
	}

	// In the case that cluster2 does not decay immediately into a single hadron,
	// Calculate the momenta of Q and Q2bar (as constituent of C2) in the
	// (parent) C2 frame first, where the direction of Q is u.vect().unit(),
	// and then boost back in the LAB frame.
	if(!toHadron2) {
	if (pClu2.m() < pQ.mass() + pQ2bar.mass() ) {
	throw Exception() << "Impossible Kinematics in ClusterFissioner::drawKinematics() (C)"
	<< Exception::eventerror;
	}
	std::pair<Axis,double> direction2 = sampleDirectionConstituents(pClu2,pClu.mass());
	// boost from Cluster2 rest frame to Cluster COM Frame
	Kinematics::twoBodyDecay(pClu2, pQ2bar.mass(), pQ.mass(), direction2.first, pQ2bar, pQ);
	weightTotal*=direction2.second;
	if (0 && _writeOut) {
	ofstream out("WriteOut/test_Constituents2.dat",std::ios::app);
	Lorentz5Momentum pQ2inCOM(pQ2bar);
	pQ2inCOM.setMass(pQ2bar.m());
	pQ2inCOM.boost( -pClu2.boostVector() ); // boost from LAB to C
	out << pQ2inCOM.vect().cosTheta(pClu2.vect()) << "\n";
	}
	}
	// Boost all momenta from the Cluster COM frame to the Lab frame
	if (_covariantBoost) {
	std::vector<Lorentz5Momentum *> momenta;
	momenta.push_back(&pClu1);
	momenta.push_back(&pClu2);
	if (!toHadron1) {
	momenta.push_back(&pQ1);
	momenta.push_back(&pQbar);
	}
	if (!toHadron2) {
	momenta.push_back(&pQ);
	momenta.push_back(&pQ2bar);
	}
	Kinematics::BoostIntoTwoParticleFrame(pClu.mass(),p0Q1, p0Q2, momenta);
	}
	else {
	pClu1.boost(pClu.boostVector());
	pClu2.boost(pClu.boostVector());
	pQ1.boost(pClu.boostVector());
	pQ2bar.boost(pClu.boostVector());
	pQ.boost(pClu.boostVector());
	pQbar.boost(pClu.boostVector());
	}
	return weightTotal; // success
	}

	void ClusterFissioner::calculateKinematics(const Lorentz5Momentum & pClu,
	const Lorentz5Momentum & p0Q1,
	const bool toHadron1,
	const bool toHadron2,
	Lorentz5Momentum & pClu1,
	Lorentz5Momentum & pClu2,
	Lorentz5Momentum & pQ1,
	Lorentz5Momentum & pQbar,
	Lorentz5Momentum & pQ,
	Lorentz5Momentum & pQ2bar) const {

	/******************
	* This method solves the kinematics of the two body cluster decay:
	* C (Q1 Q2bar) ---> C1 (Q1 Qbar) + C2 (Q Q2bar)
	* In input we receive the momentum of C, pClu, and the momentum
	* of the quark Q1 (constituent of C), p0Q1, both in the LAB frame.
	* Furthermore, two boolean variables inform whether the two fission
	* products (C1, C2) decay immediately into a single hadron (in which
	* case the cluster itself is identify with that hadron) and we do
	* not have to solve the kinematics of the components (Q1,Qbar) for
	* C1 and (Q,Q2bar) for C2.
	* The output is given by the following momenta (all 5-components,
	* and all in the LAB frame):
	* pClu1 , pClu2 respectively of C1 , C2
	* pQ1 , pQbar respectively of Q1 , Qbar in C1
	* pQ , pQ2bar respectively of Q , Q2 in C2
	* The assumption, suggested from the string model, is that, in C frame,
	* C1 and its constituents Q1 and Qbar are collinear, and collinear to
	* the direction of Q1 in C (that is before cluster decay); similarly,
	* (always in the C frame) C2 and its constituents Q and Q2bar are
	* collinear (and therefore anti-collinear with C1,Q1,Qbar).
	* The solution is then obtained by using Lorentz boosts, as follows.
	* The kinematics of C1 and C2 is solved in their parent C frame,
	* and then boosted back in the LAB. The kinematics of Q1 and Qbar
	* is solved in their parent C1 frame and then boosted back in the LAB;
	* similarly, the kinematics of Q and Q2bar is solved in their parent
	* C2 frame and then boosted back in the LAB. In each of the three
	* "two-body decay"-like cases, we use the fact that the direction
	* of the motion of the decay products is known in the rest frame of
	* their parent. This is obvious for the first case in which the
	* parent rest frame is C; but it is also true in the other two cases
	* where the rest frames are C1 and C2. This is because C1 and C2
	* are boosted w.r.t. C in the same direction where their components,
	* respectively (Q1,Qbar) and (Q,Q2bar) move in C1 and C2 rest frame
	* respectively.
	* Of course, although the notation used assumed that C = (Q1 Q2bar)
	* where Q1 is a quark and Q2bar an antiquark, indeed everything remain
	* unchanged also in all following cases:
	* Q1 quark, Q2bar antiquark; --> Q quark;
	* Q1 antiquark , Q2bar quark; --> Q antiquark;
	* Q1 quark, Q2bar diquark; --> Q quark
	* Q1 antiquark, Q2bar anti-diquark; --> Q antiquark
	* Q1 diquark, Q2bar quark --> Q antiquark
	* Q1 anti-diquark, Q2bar antiquark; --> Q quark
	**************************/

	// Calculate the unit three-vector, in the C frame, along which
	// all of the constituents and children clusters move.
	Lorentz5Momentum u(p0Q1);
	u.boost( -pClu.boostVector() ); // boost from LAB to C
	// the unit three-vector is then u.vect().unit()

	// Calculate the momenta of C1 and C2 in the (parent) C frame first,
	// where the direction of C1 is u.vect().unit(), and then boost back in the
	// LAB frame.

	if (pClu.m() < pClu1.mass() + pClu2.mass() ) {
	throw Exception() << "Impossible Kinematics in ClusterFissioner::calculateKinematics() (A)"
	<< Exception::eventerror;
	}
	Kinematics::twoBodyDecay(pClu, pClu1.mass(), pClu2.mass(),
	u.vect().unit(), pClu1, pClu2);

	// In the case that cluster1 does not decay immediately into a single hadron,
	// calculate the momenta of Q1 (as constituent of C1) and Qbar in the
	// (parent) C1 frame first, where the direction of Q1 is u.vect().unit(),
	// and then boost back in the LAB frame.
	if(!toHadron1) {
	if (pClu1.m() < pQ1.mass() + pQbar.mass() ) {
	throw Exception() << "Impossible Kinematics in ClusterFissioner::calculateKinematics() (B)"
	<< Exception::eventerror;
	}
	Kinematics::twoBodyDecay(pClu1, pQ1.mass(), pQbar.mass(),
	u.vect().unit(), pQ1, pQbar);
	}

	// In the case that cluster2 does not decay immediately into a single hadron,
	// Calculate the momenta of Q and Q2bar (as constituent of C2) in the
	// (parent) C2 frame first, where the direction of Q is u.vect().unit(),
	// and then boost back in the LAB frame.
	if(!toHadron2) {
	if (pClu2.m() < pQ.mass() + pQ2bar.mass() ) {
	throw Exception() << "Impossible Kinematics in ClusterFissioner::calculateKinematics() (C)"
	<< Exception::eventerror;
	}
	Kinematics::twoBodyDecay(pClu2, pQ.mass(), pQ2bar.mass(),
	u.vect().unit(), pQ, pQ2bar);
	}
	}

	void ClusterFissioner::calculatePositions(const Lorentz5Momentum & pClu,
	const LorentzPoint & positionClu,
	const Lorentz5Momentum & pClu1,
	const Lorentz5Momentum & pClu2,
	LorentzPoint & positionClu1,
	LorentzPoint & positionClu2) const {
	// Determine positions of cluster children.
	// See Marc Smith's thesis, page 127, formulas (4.122) and (4.123).
	Energy Mclu = pClu.m();
	Energy Mclu1 = pClu1.m();
	Energy Mclu2 = pClu2.m();

	// Calculate the unit three-vector, in the C frame, along which
	// children clusters move.
	Lorentz5Momentum u(pClu1);
	u.boost( -pClu.boostVector() ); // boost from LAB to C frame

	// the unit three-vector is then u.vect().unit()

	Energy pstarChild = Kinematics::pstarTwoBodyDecay(Mclu,Mclu1,Mclu2);

	// First, determine the relative positions of the children clusters
	// in the parent cluster reference frame.

	Energy2 mag2 = u.vect().mag2();
	InvEnergy fact = mag2>ZERO ? 1./sqrt(mag2) : 1./GeV;

	Length x1 = ( 0.25Mclu + 0.5( pstarChild + (sqr(Mclu2) - sqr(Mclu1))/(2.0*Mclu)))/_kappa;
	Length t1 = Mclu/_kappa - x1;
	LorentzDistance distanceClu1( x1 * fact * u.vect(), t1 );

	Length x2 = (-0.25Mclu + 0.5(-pstarChild + (sqr(Mclu2) - sqr(Mclu1))/(2.0*Mclu)))/_kappa;
	Length t2 = Mclu/_kappa + x2;
	LorentzDistance distanceClu2( x2 * fact * u.vect(), t2 );

	// Then, transform such relative positions from the parent cluster
	// reference frame to the Lab frame.
	distanceClu1.boost( pClu.boostVector() );
	distanceClu2.boost( pClu.boostVector() );

	// Finally, determine the absolute positions in the Lab frame.
	positionClu1 = positionClu + distanceClu1;
	positionClu2 = positionClu + distanceClu2;

	}

	bool ClusterFissioner::ProbablityFunction(double scale, double threshold) {
	double cut = UseRandom::rnd(0.0,1.0);
	return 1./(1.+pow(abs((threshold-_probShift)/scale),_probPowFactor)) > cut ? true : false;
	}
	bool ClusterFissioner::ProbablityFunctionTest(double Mass, double threshold) {
	double cut = UseRandom::rnd(0.0,1.0);
	// double arg = scale/threshold;
	// double epsilon = 1e-5;
	// if (arg<(1.0+epsilon) ){
	// // std::cout << "arg = "<<arg << std::endl;
	// arg=1.0+epsilon;
	// }
	// return 1./(1.+_probPowFactor*pow(log(arg),-2.0)) > cut ? true : false;
	if ((Mass-threshold)<=0)
	return false;
	return 1.0/(1.0 + _probPowFactor*pow(1.0/(Mass-threshold),_clPowLight)) > cut ? true : false;
	}


	bool ClusterFissioner::isHeavy(tcClusterPtr clu) {
	// particle data for constituents
	tcPDPtr cptr[3]={tcPDPtr(),tcPDPtr(),tcPDPtr()};
	bool hasDiquark=0;
	for(size_t ix=0;ix<min(clu->numComponents(),3);++ix) {
	cptr[ix]=clu->particle(ix)->dataPtr();
	// Assuming diquark masses are ordered with larger id corresponding to larger masses
	if (DiquarkMatcher::Check(*(cptr[ix]))) {
	hasDiquark=true;
	break;
	}
	}
	// different parameters for exotic, bottom and charm clusters
	double clpow = !spectrum()->isExotic(cptr[0],cptr[1],cptr[1]) ? _clPowLight : _clPowExotic;
	Energy clmax = !spectrum()->isExotic(cptr[0],cptr[1],cptr[1]) ? _clMaxLight : _clMaxExotic;
	// if no heavy quark is found in the cluster, but diquarks are present use
	// different ClMax and ClPow
	if ( hasDiquark) {
	clpow = _clPowDiquark;
	clmax = _clMaxDiquark;
	}

	for ( const long& id : spectrum()->heavyHadronizingQuarks() ) {
	if ( spectrum()->hasHeavy(id,cptr[0],cptr[1],cptr[1])) {
	clpow = _clPowHeavy[id];
	clmax = _clMaxHeavy[id];
	}
	}
	// required test for SUSY clusters, since aboveCutoff alone
	// cannot guarantee (Mc > m1 + m2 + 2*m) in cut()
	static const Energy minmass
	= getParticleData(ParticleID::d)->constituentMass();
	bool aboveCutoff = false, canSplitMinimally = false;
	// static kinematic threshold
	if(_kinematicThresholdChoice == 0) {
	aboveCutoff = (
	pow(clu->mass()*UnitRemoval::InvE , clpow)
	>
	pow(clmax*UnitRemoval::InvE, clpow)
	+ pow(clu->sumConstituentMasses()*UnitRemoval::InvE, clpow)
	);

	canSplitMinimally = clu->mass() > clu->sumConstituentMasses() + 2.0 * minmass;
	}
	// dynamic kinematic threshold
	else if(_kinematicThresholdChoice == 1) {
	//some smooth probablity function to create a dynamic thershold
	double scale = pow(clu->mass()/GeV , clpow);
	double threshold = pow(clmax/GeV, clpow)
	+ pow(clu->sumConstituentMasses()/GeV, clpow);
	aboveCutoff = ProbablityFunction(scale,threshold);

	scale = clu->mass()/GeV;
	threshold = clu->sumConstituentMasses()/GeV + 2.0 * minmass/GeV;

	canSplitMinimally = ProbablityFunction(scale,threshold);
	}
	// probablistic kinematic threshold
	else if(_kinematicThresholdChoice == 2) {
	//some smooth probablity function to create a dynamic thershold
	// double scale = pow(clu->mass()/GeV , clpow);
	// double threshold = pow(clmax/GeV, clpow)
	// + pow(clu->sumConstituentMasses()/GeV, clpow);
	// aboveCutoff = ProbablityFunction(scale,threshold);

	double Mass = clu->mass()/GeV;
	double threshold = clu->sumConstituentMasses()/GeV + 2.0 * minmass/GeV;
	aboveCutoff = ProbablityFunctionTest(Mass,threshold + clmax/GeV);

	canSplitMinimally = Mass - threshold>ZERO;
	}

	return aboveCutoff && canSplitMinimally;
	}
	diff --git a/Hadronization/ClusterFissioner.h b/Hadronization/ClusterFissioner.h
	--- a/Hadronization/ClusterFissioner.h
	+++ b/Hadronization/ClusterFissioner.h
	@@ -1,812 +1,821 @@
	// -- C++ --
	//
	// ClusterFissioner.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_ClusterFissioner_H
	#define HERWIG_ClusterFissioner_H

	#include <ThePEG/Interface/Interfaced.h>
	#include "CluHadConfig.h"
	#include "ClusterFissioner.fh"
	#include "HadronSpectrum.h"

	namespace Herwig {
	using namespace ThePEG;

	//class Cluster; // forward declaration

	/** \ingroup Hadronization
	* \class ClusterFissioner
	* \brief This class handles clusters which are too heavy.
	* \author Philip Stephens
	* \author Alberto Ribon
	* \author Stefan Gieseke
	*
	* This class does the job of chopping up either heavy clusters or beam
	* clusters in two lighter ones. The procedure is repeated recursively until
	* all of the cluster children have masses below some threshold values.
	*
	* For the beam remnant clusters, at the moment what is done is the following.
	* In the case that the soft underlying event is switched on, the
	* beam remnant clusters are tagged as not available,
	* therefore they will not be treated at all during the hadronization.
	* In the case instead that the soft underlying event is switched off,
	* then the beam remnant clusters are treated exactly as "normal" clusters,
	* with the only exception of the mass spectrum used to generate the
	* cluster children masses. For non-beam clusters, the masses of the cluster
	* children are draw from a power-like mass distribution; for beam clusters,
	* according to the value of the flag _IOpRem, either both
	* children masses are draw from a fast-decreasing exponential mass
	* distribution (case _IOpRem == 0, or, indendently by
	* _IOpRem, in the special case that the beam cluster contains two
	* beam remnants), or one mass from the exponential distribution (corresponding
	* of the cluster child with the beam remnant) and the other with the usual
	* power-like distribution (case _IOpRem == 1, which is the
	* default one, as in Herwig 6.3).
	*
	* The reason behind the use of a fast-decreasing exponential distribution
	* is that to avoid a large transverse energy from the many sequential
	* fissions that would otherwise occur due to the typical large cluster
	* mass of beam clusters. Using instead an exponential distribution
	* the masses of the two cluster children will be very small (order of
	* GeV).
	*
	* The rationale behind the implementation of the splitting of clusters
	* has been to preserve all of the information about such splitting
	* process. More explicitly a ThePEG::Step class is passed in and the
	* new clusters are added to the step as the decay products of the
	* heavy cluster. This approach has the twofold
	* advantage to provide all of the information that could be needed
	* (expecially in future developments), without any information loss,
	* and furthermore it allows a better debugging.
	*
	* @see \ref ClusterFissionerInterfaces "The interfaces"
	* defined for ClusterFissioner.
	*/
	class ClusterFissioner: public Interfaced {

	public:

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* Default constructor.
	*/
	ClusterFissioner();
	/**
	* Default destructor.
	*/
	virtual ~ClusterFissioner();

	//@}

	/** Splits the clusters which are too heavy.
	*
	* Split either heavy clusters or beam clusters recursively until all
	* children have mass below some threshold. Heavy clusters are those that
	* satisfy the condition
	* \f[ M^P > C^P + S^P \f]
	* where \f$ M \f$ is the clusters mass, \f$ P \f$ is the parameter
	* ClPow, \f$ C \f$ is the parameter ClMax and \f$ S \f$ is the
	* sum of the clusters constituent partons.
	* For beam clusters, they are split only if the soft underlying event
	* is switched off, otherwise these clusters will be tagged as unavailable
	* and they will not be treated by the hadronization altogether.
	* In the case beam clusters will be split, the procedure is exactly
	* the same as for normal non-beam clusters, with the only exception
	* of the mass spectrum from which to draw the masses of the two
	* cluster children (see method drawChildrenMasses for details).
	*/
	tPVector fission(ClusterVector & clusters, bool softUEisOn);

	/**
	* Return the hadron spectrum
	*/
	Ptr<HadronSpectrum>::tptr spectrum() const {
	return _hadronSpectrum;
	}

	public:

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

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

	/**
	* Standard Init function used to initialize the interfaces.
	*/
	static void Init();

	protected:

	/** @name Clone Methods. */
	//@{
	/**
	* Make a simple clone of this object.
	* @return a pointer to the new object.
	*/
	virtual IBPtr clone() const;

	/** Make a clone of this object, possibly modifying the cloned object
	* to make it sane.
	* @return a pointer to the new object.
	*/
	virtual IBPtr fullclone() const;
	//@}

	private:

	/**
	* Private and non-existent assignment operator.
	*/
	ClusterFissioner & operator=(const ClusterFissioner &) = delete;

	/**
	* This method directs the splitting of the heavy clusters
	*
	* This method does the splitting of the clusters and all of its cluster
	* children, if heavy. All of these new children clusters are added to the
	* collection of clusters. The method works as follows.
	* Initially the vector contains just the stack of input pointers to the
	* clusters to be split. Then it will be filled recursively by all
	* of the cluster's children that are heavy enough to require
	* to be split. In each loop, the last element of the vector is
	* considered (only once because it is then removed from the vector).
	*
	* \todo is the following still true?
	* For normal, non-beam clusters, a power-like mass distribution
	* is used, whereas for beam clusters a fast-decreasing exponential mass
	* distribution is used instead. This avoids many iterative splitting which
	* could produce an unphysical large transverse energy from a supposed
	* soft beam remnant process.
	*/
	void cut(stack<ClusterPtr> &,
	ClusterVector&, tPVector & finalhadrons, bool softUEisOn);

	public:

	/**
	* Definition for easy passing of two particles.
	*/
	typedef pair<PPtr,PPtr> PPair;

	/**
	* Definition for use in the cut function.
	*/
	typedef pair<PPair,PPair> cutType;

	/**
	* Splits the input cluster.
	*
	* Split the input cluster (which can be either an heavy non-beam
	* cluster or a beam cluster). The result is two pairs of particles. The
	* first element of each pair is new cluster/hadron, while the second
	* element of each pair is the particle drawn from the vacuum to create
	* the new cluster/hadron.
	* Notice that this method treats also beam clusters by using a different
	* mass spectrum used to generate the cluster child masses (see method
	* drawChildMass).
	*/
	//@{
	/**
	* Split two-component cluster
	*/
	virtual cutType cutTwo(ClusterPtr &, tPVector & finalhadrons, bool softUEisOn);

	/**
	* Split two-component cluster using Default FissionApproach
	*/
	virtual cutType cutTwoDefault(ClusterPtr &, tPVector & finalhadrons, bool softUEisOn);

	/**
	* Split two-component cluster using New FissionApproach
	*/
	virtual cutType cutTwoNew(ClusterPtr &, tPVector & finalhadrons, bool softUEisOn);


	/**
	* Split three-component cluster
	*/
	virtual cutType cutThree(ClusterPtr &, tPVector & finalhadrons, bool softUEisOn);
	//@}
	public:

	/**
	* Produces a hadron and returns the flavour drawn from the vacuum.
	*
	* This routine produces a new hadron. It
	* also sets the momentum and vertex to the values given.
	*/
	PPair produceHadron(tcPDPtr hadron, tPPtr newPtr, const Lorentz5Momentum &a,
	const LorentzPoint &b) const;
	protected:

	/**
	* Produces a cluster from the flavours passed in.
	*
	* This routine produces a new cluster with the flavours given by ptrQ and newPtr.
	* The new 5 momentum is a and the parent momentum are c and d. C is for the
	* ptrQ and d is for the new particle newPtr. rem specifies whether the existing
	* particle is a beam remnant or not.
	*/
	PPair produceCluster(tPPtr ptrQ, tPPtr newPtr, const Lorentz5Momentum &a,
	const LorentzPoint &b, const Lorentz5Momentum &c,
	const Lorentz5Momentum &d, const bool rem,
	tPPtr spect=tPPtr(), bool remSpect=false) const;

	/**
	* Returns the new quark-antiquark pair
	* needed for fission of a heavy cluster. Equal probabilities
	* are assumed for producing u, d, or s pairs.
	*/
	void drawNewFlavourQuarks(PPtr& newPtrPos,PPtr& newPtrNeg) const;

	/**
	* Returns the new quark-antiquark pair or diquark -
	* antidiquark pair needed for fission of a heavy cluster.
	*/
	void drawNewFlavourDiquarks(PPtr& newPtrPos,PPtr& newPtrNeg,
	const ClusterPtr & clu) const;

	/**
	* Returns the new quark-antiquark pair
	* needed for fission of a heavy cluster. Equal probabilities
	* are assumed for producing u, d, or s pairs.
	* Extra argument is used when performing strangeness enhancement
	*/
	void drawNewFlavourEnhanced(PPtr& newPtrPos,PPtr& newPtrNeg, Energy2 mass2) const;


	/**
	* Produces the mass of a child cluster.
	*
	* Draw the masses \f$M'\f$ of the the cluster child produced
	* by the fission of an heavy cluster (of mass M). m1, m2 are the masses
	* of the constituents of the cluster; m is the mass of the quark extract
	* from the vacuum (together with its antiparticle). The algorithm produces
	* the mass of the cluster formed with consituent m1.
	* Two mass distributions can be used for the child cluster mass:
	* -# power-like mass distribution ("normal" mass) with power exp
	* \f[ M' = {\rm rnd}((M-m_1-m_2-m)^P, m^p)^{1/P} + m_1 \f]
	* where \f$ P \f$ is a parameter of the model and \f$ \rm{rnd} \f$ is
	* the function:
	* \f[ \rm{rnd}(a,b) = (1-r)a + r b \f]
	* and here \f$ r \f$ is a random number [0,1].
	* -# fast-decreasing exponential mass distribution ("soft" mass) with
	* rmin. rmin is given by
	* \f[ r_{\rm min} = \exp(-b (M - m_1 - m_2 - 2 m)) \f]
	* where \f$ b \f$ is a parameter of the model. The generated mass is
	* given by
	* \f[ M' = m_1 + m - \frac{\log\left(
	* {\rm rnd}(r_{\rm min}, 1-r_{\rm min})\right)}{b} \f].
	*
	* The choice of which mass distribution should be used for each of the two
	* cluster children is dictated by the parameter soft.
	*/
	Energy drawChildMass(const Energy M, const Energy m1, const Energy m2,
	const Energy m, const double exp, const bool soft) const;

	/**
	* Determine the positions of the two children clusters.
	*
	* This routine generates the momentum of the decay products. It also
	* generates the momentum in the lab frame of the partons drawn out of
	* the vacuum.
	*/
	void calculatePositions(const Lorentz5Momentum &pClu,
	const LorentzPoint & positionClu,
	const Lorentz5Momentum & pClu1,
	const Lorentz5Momentum & pClu2,
	LorentzPoint & positionClu1,
	LorentzPoint & positionClu2 ) const;

	protected:

	/**
	* Dimension used to calculate phase space weights
	*/
	double dim() const {return _dim;}

	/**
	* Access to soft-cluster parameter
	*/
	Energy btClM() const {return _btClM;}

	/**
	* Function that returns either the cluster mass or the lambda measure
	*/
	Energy2 clustermass(const ClusterPtr & cluster) const;

	/**
	* old calculateKinematics function for Default FissionApproach
	*/
	void calculateKinematics(const Lorentz5Momentum & pClu,
	const Lorentz5Momentum & p0Q1,
	const bool toHadron1,
	const bool toHadron2,
	Lorentz5Momentum & pClu1,
	Lorentz5Momentum & pClu2,
	Lorentz5Momentum & pQ1,
	Lorentz5Momentum & pQbar,
	Lorentz5Momentum & pQ,
	Lorentz5Momentum & pQ2bar) const;

	/**
	* Draw a new flavour for the given cluster; currently defaults to
	* the default model
	*/
	virtual void drawNewFlavour(PPtr& newPtr1, PPtr& newPtr2, const ClusterPtr & cluster) const {
	if (_enhanceSProb == 0){
	if (_diquarkClusterFission>=0) drawNewFlavourDiquarks(newPtr1,newPtr2,cluster);
	else drawNewFlavourQuarks(newPtr1,newPtr2);
	}
	else {
	drawNewFlavourEnhanced(newPtr1,newPtr2,clustermass(cluster));
	}
	}

	/**
	* Calculate the masses and possibly kinematics of the cluster
	* fission at hand; if claculateKineamtics is perfomring non-trivial
	* steps kinematics claulcated here will be overriden. Currentl;y resorts to the default
	* @return the potentially non-trivial distribution weight=f(M1,M2)
	* On Failure we return 0
	*/
	double drawNewMasses(const Energy Mc, const bool soft1, const bool soft2,
	Lorentz5Momentum& pClu1, Lorentz5Momentum& pClu2,
	tcPPtr ptrQ1, const Lorentz5Momentum& pQ1,
	tcPPtr, const Lorentz5Momentum& pQone,
	tcPPtr, const Lorentz5Momentum& pQtwo,
	tcPPtr ptrQ2, const Lorentz5Momentum& pQ2) const;
	/**
	* Calculate the masses and possibly kinematics of the cluster
	* fission at hand; if claculateKineamtics is perfomring non-trivial
	* steps kinematics claulcated here will be overriden. Currentl;y resorts to the default
	*/
	double drawNewMassesDefault(const Energy Mc, const bool soft1, const bool soft2,
	Lorentz5Momentum& pClu1, Lorentz5Momentum& pClu2,
	tcPPtr ptrQ1, const Lorentz5Momentum& pQ1,
	tcPPtr, const Lorentz5Momentum& pQone,
	tcPPtr, const Lorentz5Momentum& pQtwo,
	tcPPtr ptrQ2, const Lorentz5Momentum& pQ2) const;

	/**
	* Sample the masses for flat phase space
	* */
	double drawNewMassesUniform(const Energy Mc,
	Lorentz5Momentum & pClu1, Lorentz5Momentum & pClu2,
	const Lorentz5Momentum & pQ1,
	const Lorentz5Momentum & pQ,
	const Lorentz5Momentum & pQ2) const;

	/**
	* Sample the masses for flat phase space
	* */
	double drawNewMassesPhaseSpace(const Energy Mc,
	Lorentz5Momentum & pClu1, Lorentz5Momentum & pClu2,
	const Lorentz5Momentum & pQ1,
	const Lorentz5Momentum & pQ,
	const Lorentz5Momentum & pQ2,
	tcPPtr ptrQ1, tcPPtr ptrQ2, tcPPtr ptrQ ) const;
	+ /**
	+ * Sample the masses for flat phase space modulated by a power law
	+ * */
	+ double drawNewMassesPhaseSpacePowerLaw(const Energy Mc,
	+ Lorentz5Momentum& pClu1, Lorentz5Momentum& pClu2,
	+ const Lorentz5Momentum& pQ1,
	+ const Lorentz5Momentum& pQ,
	+ const Lorentz5Momentum& pQ2,
	+ tcPPtr ptrQ1, tcPPtr ptrQ2, tcPPtr ptrQ) const;

	/**
	* Calculate the final kinematics of a heavy cluster decay C->C1 +
	* C2, if not already performed by drawNewMasses
	* @return returns false if failes
	*/
	double drawKinematics(const Lorentz5Momentum &pClu,
	const Lorentz5Momentum &p0Q1,
	const Lorentz5Momentum &p0Q2,
	const bool toHadron1, const bool toHadron2,
	Lorentz5Momentum &pClu1, Lorentz5Momentum &pClu2,
	Lorentz5Momentum &pQ1, Lorentz5Momentum &pQb,
	Lorentz5Momentum &pQ, Lorentz5Momentum &pQ2b) const;

	/**
	* Calculation of the squared matrix element from the drawn
	* momenta
	* @return value of the squared matrix element in units of
	* 1/GeV4
	* */
	double calculateSQME(
	const Lorentz5Momentum & p1,
	const Lorentz5Momentum & p2,
	const Lorentz5Momentum & q1,
	const Lorentz5Momentum & q,
	const Lorentz5Momentum & q2,
	const Lorentz5Momentum & qbar) const;

	/**
	* Calculation of the overestimate for the squared matrix
	* element independent on M1 and M2
	* @return value of the overestimate squared matrix element
	* in units of 1/GeV4
	* */
	double calculateSQME_OverEstimate(
	const Energy& Mc,
	const Energy& m1,
	const Energy& m2,
	const Energy& mq) const;
	protected:

	/** @name Access members for child classes. */
	//@{
	/**
	* Access to the hadron selector
	*/
	HadronSpectrumPtr hadronSpectrum() const {return _hadronSpectrum;}
	//@}

	protected:

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

	//@}

	private:

	std::pair<Axis,double> sampleDirectionCluster(const Lorentz5Momentum & pQ, const Lorentz5Momentum & pClu) const;
	std::pair<Axis,double> sampleDirectionConstituents(const Lorentz5Momentum & pClu, const Energy Mcluster) const;
	/**
	* Samples the direction of Cluster Fission products uniformly
	**/
	Axis sampleDirectionIsotropic() const;

	/**
	* Samples the direction of Cluster Fission products uniformly
	* but only accepts those flying in the direction of pRelCOM
	**/
	Axis sampleDirectionSemiUniform(const Lorentz5Momentum & pQ, const Lorentz5Momentum & pClu) const;

	/**
	* Samples the direction of Cluster Fission products according to Default
	* fully aligned D = 1+1 Fission
	* */
	Axis sampleDirectionAligned(const Lorentz5Momentum & pQ, const Lorentz5Momentum & pClu) const;

	/**
	* Samples the direction of Cluster Fission products according to Default
	* Tchannel like distribution
	* */
	std::pair<Axis,double> sampleDirectionTchannel(const Axis dirQ, const double Ainv) const;

	/**
	* Samples the direction of direction from an exponential distribution
	* */
	std::pair<Axis,double> sampleDirectionExponential(const Axis dirQ, const double lambda) const;

	/**
	* Samples the direction of Cluster Fission products according to smeared direction along the cluster
	* */
	Axis sampleDirectionAlignedSmeared(const Lorentz5Momentum & pQ, const Lorentz5Momentum & pClu) const;

	/**
	* Smooth probability for dynamic threshold cuts:
	* @scale the current scale, e.g. the mass of the cluster,
	* @threshold the physical threshold,
	*/
	bool ProbablityFunction(double scale, double threshold);
	bool ProbablityFunctionTest(double Mass, double threshold);

	/**
	* Check if a cluster is heavy enough to split again
	*/
	bool isHeavy(tcClusterPtr );

	/**
	* Check if a cluster is heavy enough to be at least kinematically able to split
	*/
	bool canSplitMinimally(tcClusterPtr, Energy);

	/**
	* Check if can't make a hadron from the partons
	*/
	inline bool cantMakeHadron(tcPPtr p1, tcPPtr p2) {
	return ! spectrum()->canBeHadron(p1->dataPtr(), p2->dataPtr());
	}
	/**
	* Flat PhaseSpace weight for ClusterFission
	*/
	double weightFlatPhaseSpace(const Energy Mc, const Energy Mc1, const Energy Mc2,
	const Energy m, const Energy m1, const Energy m2,
	tcPPtr pQ, tcPPtr pQ1, tcPPtr pQ2) const {
	switch (_phaseSpaceWeights)
	{
	case 1:
	- return weightFlatPhaseSpaceConstituentMasses(Mc, Mc1, Mc2, m, m1, m2);
	+ return weightPhaseSpaceConstituentMasses(Mc, Mc1, Mc2, m, m1, m2, 0.0);
	case 2:
	return weightFlatPhaseSpaceHadronMasses(Mc, Mc1, Mc2, pQ, pQ1, pQ2);
	case 3:
	return weightFlatPhaseSpaceNoConstituentMasses(Mc, Mc1, Mc2);
	default:
	assert(false);
	}
	};
	/**
	- * Flat PhaseSpace weight for ClusterFission using constituent masses
	+ * PhaseSpace weight for ClusterFission using constituent masses
	*/
	- double weightFlatPhaseSpaceConstituentMasses(const Energy Mc, const Energy Mc1, const Energy Mc2,
	- const Energy m, const Energy m1, const Energy m2) const;
	+ double weightPhaseSpaceConstituentMasses(const Energy Mc, const Energy Mc1, const Energy Mc2,
	+ const Energy m, const Energy m1, const Energy m2, const double power=0.0) const;
	/**
	* Flat PhaseSpace weight for ClusterFission using lightest hadron masses
	*/
	double weightFlatPhaseSpaceHadronMasses(const Energy Mc, const Energy Mc1, const Energy Mc2,
	tcPPtr pQ, tcPPtr pQ1, tcPPtr pQ2) const;
	double weightFlatPhaseSpaceNoConstituentMasses(const Energy Mc, const Energy Mc1, const Energy Mc2) const;


	/**
	* Calculate a veto for the phase space weight
	*/
	bool phaseSpaceVeto(const Energy Mc, const Energy Mc1, const Energy Mc2,
	- const Energy m, const Energy m1, const Energy m2, tcPPtr pQ1=tcPPtr(), tcPPtr pQ2=tcPPtr(), tcPPtr pQ=tcPPtr()) const;
	+ const Energy m, const Energy m1, const Energy m2, tcPPtr pQ1=tcPPtr(), tcPPtr pQ2=tcPPtr(), tcPPtr pQ=tcPPtr(), const double power = 0.0) const;
	bool phaseSpaceVetoConstituentMasses(const Energy Mc, const Energy Mc1, const Energy Mc2,
	- const Energy m, const Energy m1, const Energy m2) const;
	+ const Energy m, const Energy m1, const Energy m2, const double power = 0.0) const;
	bool phaseSpaceVetoNoConstituentMasses(const Energy Mc, const Energy Mc1, const Energy Mc2) const;

	bool phaseSpaceVetoHadronPairs(const Energy Mc, const Energy Mc1, const Energy Mc2,
	tcPPtr pQ1, tcPPtr pQ2, tcPPtr pQconst) const;


	void countPaccGreater1(){
	_counterPaccGreater1++;
	};
	void countMaxLoopViolations(){
	_counterMaxLoopViolations++;
	};
	void countFissionMatrixElement(){
	_counterFissionMatrixElement++;
	};


	/**
	* A pointer to a Herwig::HadronSpectrum object for generating hadrons.
	*/
	HadronSpectrumPtr _hadronSpectrum;

	/**
	* @name The Cluster max mass,dependant on which quarks are involved, used to determine when
	* fission will occur.
	*/
	//@{
	Energy _clMaxLight;
	Energy _clMaxDiquark;
	map<long,Energy> _clMaxHeavy;
	Energy _clMaxExotic;
	//@}
	/**
	* @name The power used to determine when cluster fission will occur.
	*/
	//@{
	double _clPowLight;
	double _clPowDiquark;
	map<long,double> _clPowHeavy;
	double _clPowExotic;
	//@}
	/**
	* @name The power, dependant on whic quarks are involved, used in the cluster mass generation.
	*/
	//@{
	double _pSplitLight;
	map<long,double> _pSplitHeavy;
	double _pSplitExotic;

	/**
	* Weights for alternative cluster fission
	*/
	map<long,double> _fissionPwt;

	/**
	* Include phase space weights
	*/
	int _phaseSpaceWeights;

	/**
	* Dimensionality of phase space weight
	*/
	double _dim;

	/**
	* Flag used to determine between normal cluster fission and alternative cluster fission
	*/
	int _fissionCluster;

	/**
	* Flag to choose static or dynamic kinematic thresholds in cluster splittings
	*/
	int _kinematicThresholdChoice;

	/**
	* Pwt weight for drawing diquark
	*/
	double _pwtDIquark;

	/**
	* allow clusters to fission to 1 (or 2) diquark clusters or not
	*/
	int _diquarkClusterFission;

	//@}
	/**
	* Parameter used (2/b) for the beam cluster mass generation.
	* Currently hard coded value.
	*/
	Energy _btClM;

	/**
	* Flag used to determine what distributions to use for the cluster masses.
	*/
	int _iopRem;

	/**
	* The string constant
	*/
	Tension _kappa;

	/**
	* Width of the gaussian sampling for the FluxTube Kinematics
	*/
	double _fluxTubeWidth;

	/**
	* Flag that switches between no strangeness enhancement, scaling enhancement,
	* and exponential enhancement (in numerical order)
	*/
	int _enhanceSProb;

	/**
	* Parameter that governs the strangeness enhancement scaling
	*/
	Energy _m0Fission;

	/**
	* Flag that switches between mass measures used in strangeness enhancement:
	* cluster mass, or the lambda measure - ( m_{clu}^2 - (m_q + m_{qbar})^2 )
	*/
	int _massMeasure;

	/**
	* Constant variable which stops the scale from being to large, and not worth
	* calculating
	*/
	const double _maxScale = 20.;

	/**
	* Power factor in ClausterFissioner bell probablity function
	*/
	double _probPowFactor;

	/**
	* Shifts from the center in ClausterFissioner bell probablity function
	*/
	double _probShift;

	/**
	* Shifts from the kinetic threshold in ClausterFissioner
	*/
	Energy2 _kinThresholdShift;

	/**
	* Flag for strict diquark selection according to kinematics
	*/
	int _strictDiquarkKinematics;

	/**
	* Use Covariant boost in ClusterFissioner
	*/
	bool _covariantBoost;

	/**
	* Flag for allowing Hadron Final states in Cluster Fission
	*/
	int _allowHadronFinalStates;

	/**
	* Choice of Mass sampling for ClusterFissioner
	* i.e. rejection sampling starting from MassPresampling
	* Chooses the to be sampled Mass distribution
	* Note: This ideal distribution would be ideally
	* exactly the integral over the angles as a
	* function of M1,M2
	*/
	int _massSampler;

	/**
	* Choice of Phase Space sampling for ClusterFissioner
	* for child cluster directions and constituent directions
	*/
	int _phaseSpaceSamplerCluster;
	int _phaseSpaceSamplerConstituent;

	/**
	* Choice of Matrix Element for ClusterFissioner
	* Note: The choice of the matrix element requires
	* to provide one overestimate as a function
	* of M1,M2 and another overestimate of the
	* previous overestimate independent of
	* M1 and M2 i.e. only dependent on M,m1,m2,m
	*/
	int _matrixElement;

	/**
	* Choice of ClusterFissioner Approach
	*/
	int _fissionApproach;

	/**
	* Power for MassPreSampler = PowerLaw
	*/
	double _powerLawPower;

	/**
	* Choice of ClusterFissioner Approach
	* Technical Parameter for how many tries are allowed to sample the
	* Cluster Fission matrix element before reverting to fissioning
	* using the default Fission Aproach
	*/
	int _maxLoopFissionMatrixElement;

	/*
	* Safety factor for a better overestimate of the matrix Element
	*
	* */
	double _safetyFactorMatrixElement;

	/*
	* IR cutOff for IR divergent diagramms
	* */
	double _epsilonIR;

	/*
	* Failing mode for MaxLoopMatrixElement
	* */
	int _failModeMaxLoopMatrixElement;

	int _writeOut;

	/*
	* flag for allowing strange Diquarks to be produced during
	* Cluster Fission
	* */
	unsigned int _hadronizingStrangeDiquarks;


	private:
	// For MatrixElement only:

	// Counter for how many times we accept events with Pacc>1
	static unsigned int _counterPaccGreater1;
	// Counter for how many times we escape the sampling by
	// tries > _maxLoopFissionMatrixElement
	// NOTE: This maxing out either rejects event or uses old
	// ClusterFission
	static unsigned int _counterMaxLoopViolations;
	static unsigned int _counterFissionMatrixElement;
	};

	}

	#endif /* HERWIG_ClusterFissioner_H */

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 9:06 PM (22 h, 37 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3803401
Default Alt Text: (167 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions