No OneTemporary
Actions

Size

29 KB

Subscribers

None

View Options

	diff --git a/Hadronization/PartonSplitter.cc b/Hadronization/PartonSplitter.cc
	--- a/Hadronization/PartonSplitter.cc
	+++ b/Hadronization/PartonSplitter.cc
	@@ -1,542 +1,596 @@
	// -- C++ --
	//
	// PartonSplitter.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2019 The Herwig Collaboration
	//
	// Herwig is licenced under version 3 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the PartonSplitter class.
	//

	#include "PartonSplitter.h"
	#include <ThePEG/Interface/ClassDocumentation.h>
	#include <ThePEG/Interface/Reference.h>
	#include <ThePEG/Interface/Switch.h>
	#include <ThePEG/Persistency/PersistentOStream.h>
	#include <ThePEG/Persistency/PersistentIStream.h>
	#include <ThePEG/PDT/EnumParticles.h>
	#include <ThePEG/EventRecord/Step.h>
	#include "ThePEG/Interface/Parameter.h"
	#include <ThePEG/Repository/EventGenerator.h>
	#include <ThePEG/Repository/CurrentGenerator.h>
	#include "ThePEG/Repository/UseRandom.h"
	#include "Herwig/Utilities/Kinematics.h"
	#include <ThePEG/Utilities/DescribeClass.h>
	#include "ClusterHadronizationHandler.h"
	#include <ThePEG/EventRecord/Particle.h>
	#include <ThePEG/PDT/PDT.h>
	#include "ThePEG/Interface/ParMap.h"

	using namespace Herwig;

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

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

	void PartonSplitter::persistentOutput(PersistentOStream & os) const {
	os << _quarkSelector << ounit(_gluonDistance,femtometer)
	- << _splitGluon << _splitPwt
	+ << _splitGluon << _splitPwt << _kinematics
	<< _enhanceSProb << ounit(_m0,GeV) << _massMeasure
	<< _hadronSpectrum;
	}

	void PartonSplitter::persistentInput(PersistentIStream & is, int) {
	is >> _quarkSelector >> iunit(_gluonDistance,femtometer)
	- >> _splitGluon >> _splitPwt
	+ >> _splitGluon >> _splitPwt >> _kinematics
	>>_enhanceSProb >> iunit(_m0,GeV) >> _massMeasure
	>> _hadronSpectrum;
	}

	DescribeClass<PartonSplitter,Interfaced>
	describePartonSplitter("Herwig::PartonSplitter","Herwig.so");


	void PartonSplitter::Init() {

	static ClassDocumentation<PartonSplitter> documentation
	("This class is reponsible of the nonperturbative splitting of partons");

	static Switch<PartonSplitter,int> interfaceSplit
	("Split",
	"Option for different splitting options",
	&PartonSplitter::_splitGluon, 0, false, false);
	static SwitchOption interfaceSplitDefault
	(interfaceSplit,
	"Kinematic",
	"All kineamtically allowed quarks will contribute to gluon splitting.",
	0);
	static SwitchOption interfaceSplitAll
	(interfaceSplit,
	"Light",
	"Alternative cluster splitting where only lighht quarks can be drawn.",
	1);

	+ static Switch<PartonSplitter,int> interfaceKinematics
	+ ("Kinematics",
	+ "Option for different splitting Kinematic options",
	+ &PartonSplitter::_kinematics, 0, false, false);
	+ static SwitchOption interfaceKinematicsIsotropic
	+ (interfaceKinematics,
	+ "Isotropic",
	+ "All kinematics are done isotropically (default)",
	+ 0);
	+ static SwitchOption interfaceKinematicsAligned
	+ (interfaceKinematics,
	+ "Aligned",
	+ "All kinematics are done aligned to the gluon"
	+ "NOTE: This option is experimental and just for debugging!!",
	+ 1);
	+
	+
	static ParMap<PartonSplitter,double> interfaceSplitPwt
	("SplitPwt",
	"The splitting weights for individual quarks.",
	&PartonSplitter::_splitPwt, -1, 1.0, 0.0, 10.0,
	false, false, Interface::upperlim);

	static Switch<PartonSplitter,int> interfaceEnhanceSProb
	("EnhanceSProb",
	"Option to enhance the strangeness weight (MassSplit Switch needs to be on)",
	&PartonSplitter::_enhanceSProb,0,false,false);
	static SwitchOption interfaceEnhanceSProbNo
	(interfaceEnhanceSProb,
	"No",
	"No scaling for strangeness",
	0);
	static SwitchOption interfaceEnhanceSProbScale
	(interfaceEnhanceSProb,
	"Scaled",
	"Power-law scaling for strangeness",
	1);
	static SwitchOption interfaceEnhanceSProbExp
	(interfaceEnhanceSProb,
	"Exponential",
	"Exponential suppression for strangeness",
	2);

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

	static Parameter<PartonSplitter,Energy> interfaceMassScale
	("MassScale",
	"Mass scale for g->qqb strangeness enhancement",
	&PartonSplitter::_m0, GeV, 20.GeV, 1.GeV, 1000.*GeV,
	false, false, Interface::limited);

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

	}

	void PartonSplitter::split(PVector & tagged) {
	// set the gluon c tau once and for all
	static bool first = true;

	if(first) {
	_gluonDistance = hbarc*getParticleData(spectrum()->gluonId())->constituentMass()/
	ClusterHadronizationHandler::currentHandler()->minVirtuality2();
	first = false;
	}
	// Copy incoming for the (possible sorting and) splitting
	PVector particles = tagged;
	// Switch to enhance strangeness
	if (_enhanceSProb >= 1){
	colourSorted(particles);
	}

	PVector newtag;
	// Loop over all of the particles in the event.
	// Loop counter for colourSorted
	for(PVector::iterator pit = particles.begin(); pit!=particles.end(); ++pit) {
	// only considering gluons so add other particles to list of particles
	if( (**pit).data().id() != spectrum()->gluonId() ) {
	newtag.push_back(*pit);
	continue;
	}
	Energy2 Q02 = 0.99(*pit).momentum().m2();
	// should not have been called for massless or space-like gluons
	if((*pit).momentum().m2() <= 0.0sqr(MeV) ) {
	throw Exception()
	<< "Spacelike or massless gluon m2= " << (**pit).momentum().m2()/GeV2
	<< "GeV2 in PartonSplitter::split()"
	<< Exception::eventerror;
	}
	// time like gluon gets split
	PPtr ptrQ = PPtr();
	PPtr ptrQbar = PPtr();
	if (_enhanceSProb == 0){
	splitTimeLikeGluon(*pit,ptrQ,ptrQbar);
	}
	else {
	size_t i = pit - particles.begin();
	massSplitTimeLikeGluon(*pit, ptrQ, ptrQbar, i);
	}

	ptrQ->scale(Q02);
	ptrQbar->scale(Q02);

	(*pit)->colourLine()->addColoured(ptrQ);
	(*pit)->addChild(ptrQ);
	newtag.push_back(ptrQ);

	(*pit)->antiColourLine()->addAntiColoured(ptrQbar);
	(*pit)->addChild(ptrQbar);
	newtag.push_back(ptrQbar);

	// set the life length of gluon
	Length distance = UseRandom::rndExp(_gluonDistance);
	(pit).setLifeLength((distance/(pit).mass())(*pit).momentum());
	// assume quarks same position as gluon
	ptrQ ->setVertex((**pit).decayVertex());
	ptrQ ->setLifeLength(Lorentz5Distance());
	ptrQbar->setVertex((**pit).decayVertex());
	ptrQbar->setLifeLength(Lorentz5Distance());
	}
	swap(tagged,newtag);
	}

	void PartonSplitter::splitTimeLikeGluon(tcPPtr ptrGluon,
	PPtr & ptrQ,
	PPtr & ptrQbar){
	// select the quark flavour
	tPDPtr quark = _quarkSelector.select(UseRandom::rnd());
	// Solve the kinematics of the two body decay G --> Q + Qbar
	Lorentz5Momentum momentumQ;
	Lorentz5Momentum momentumQbar;
	- double cosThetaStar = UseRandom::rnd( -1.0 , 1.0 );
	- using Constants::pi;
	- double phiStar = UseRandom::rnd( -pi , pi );

	Energy constituentQmass = quark->constituentMass();

	if (ptrGluon->momentum().m() < 2.0*constituentQmass) {
	throw Exception() << "Impossible Kinematics in PartonSplitter::splitTimeLikeGluon()"
	<< Exception::eventerror;
	}
	- Kinematics::twoBodyDecay(ptrGluon->momentum(), constituentQmass,
	- constituentQmass, cosThetaStar, phiStar, momentumQ,
	- momentumQbar );
	+ switch (_kinematics)
	+ {
	+ case 0:
	+ {
	+ // Isotropic default
	+ double cosThetaStar = UseRandom::rnd( -1.0 , 1.0 );
	+ using Constants::pi;
	+ double phiStar = UseRandom::rnd( -pi , pi );
	+ Kinematics::twoBodyDecay(ptrGluon->momentum(), constituentQmass,
	+ constituentQmass, cosThetaStar, phiStar, momentumQ,
	+ momentumQbar );
	+ break;
	+ }
	+ case 1:
	+ {
	+ // Aligned
	+ Kinematics::twoBodyDecay(ptrGluon->momentum(), constituentQmass,
	+ constituentQmass, ptrGluon->momentum().vect().unit(), momentumQ,
	+ momentumQbar );
	+ break;
	+ }
	+ default:
	+ assert(false);
	+ }
	+

	// Create quark and anti-quark particles of the chosen flavour
	// and set they 5-momentum (the mass is the constituent one).
	ptrQ = new_ptr(Particle(quark ));
	ptrQbar = new_ptr(Particle(quark->CC()));

	ptrQ ->set5Momentum( momentumQ );
	ptrQbar ->set5Momentum( momentumQbar );

	}

	void PartonSplitter::doinit() {
	Interfaced::doinit();
	switch(_splitGluon){
	case 0: {
	// calculate the probabilties for the gluon to branch into each quark type
	// based on the available phase-space, as in fortran.
	Energy mg=getParticleData(spectrum()->gluonId())->constituentMass();
	for( const long& ix : spectrum()->hadronizingQuarks() ) {
	PDPtr quark = getParticleData(ix);
	Energy pcm = Kinematics::pstarTwoBodyDecay(mg,quark->constituentMass(),
	quark->constituentMass());
	if(pcm>ZERO) _quarkSelector.insert(pcm/GeV,quark);
	}
	break;
	}
	case 1:
	for ( const long& ix : spectrum()->lightHadronizingQuarks() ) {
	PDPtr quark = getParticleData(ix);
	if (_splitPwt.find(ix) == _splitPwt.end() )
	throw InitException() << "no split weight found for quark flavour "
	<< ix;
	_quarkSelector.insert(_splitPwt[ix],quark);
	}
	break;
	}
	if(_quarkSelector.empty())
	throw InitException() << "At least one quark must have constituent mass less "
	<< "then the constituent mass of the gluon in "
	<< "PartonSplitter::doinit()" << Exception::runerror;
	}

	// Method to colour sort the event and calculate the masses of the
	// pre-clusters
	// Convention is to have the triplet of the colour singlet first,
	// then all gluons, then the antitriplet (and repeat for all the
	// colour-singlets in the event)
	void PartonSplitter::colourSorted(PVector& tagged) {
	// Set up the output
	PVector sorted;
	// Reset the storage variables for doing the mass-based strangeness
	// enhancement
	_colSingletSize.resize(0);
	_colSingletm2.resize(0);
	// Variable to exit out of gluon loops
	bool gluonLoop = false;
	// Partons left to consider
	PVector left = tagged;
	// Loop over all the tagged particles
	while (int(left.size()) > 0){
	// Pick the first particle available
	PPtr p = left[0];
	// Temporary holding pen for momenta
	Lorentz5Momentum tempMom(ZERO, ZERO, ZERO, ZERO);
	// Don't do anything if the particle has no colour
	// Simply add it back into the list of particles
	if ( !p->coloured() ) {
	sorted.push_back(p);
	// nparts is the index of the particle after adding it to the sorted list
	int nparts = 1;
	// Add on the last entry of colSingletSize if the vector is not empty
	// This is essentially the index the particle will have once it
	// Get placed into the new colour sorted event
	nparts += (_colSingletSize.empty()) ? 0 : _colSingletSize.back();
	tempMom += p->momentum();
	Energy2 singletm2 = tempMom.m2();
	// Store the number of particles and mass of the colour-singlet
	_colSingletSize.push_back(nparts);
	_colSingletm2.push_back(singletm2);
	// Remove the particle from the list of particles left
	left.erase(remove(left.begin(),left.end(), p), left.end());
	}
	// Temporary holding pen for partons
	PVector temp;
	// Variable to sum end-point masses i.e. triplets and anti-triplets
	Energy endPointMass = ZERO;
	// If the particle in question is a gluon, search up it's antiColourLine
	// until we get to the triplet.
	// Note there are situations where we have a gluon loop
	if ( p->hasColour() && p->hasAntiColour() ){
	// Search up its anticolour line until we get to the start particle
	tPPtr partner = p->antiColourLine()->endParticle();
	// Dummy index used to loop over the anticolour line trace
	tPPtr dummy = partner;
	// Store the partner particle
	temp.push_back(partner);
	tempMom += partner->momentum();
	left.erase(remove(left.begin(),left.end(), partner), left.end());
	// While loop continues while we still reach a particle with with
	// anti-colour, i.e. a gluon
	while ( dummy->hasAntiColour() ){
	dummy = dummy->antiColourLine()->endParticle();
	// Check that we haven't already included it via colour indices
	// If we have, it is a gluon loop.
	if ( find(left.begin(), left.end(), dummy) == left.end() ) {
	gluonLoop = true;
	break;
	}
	// Store the dummy partons in reverse
	temp.push_back(dummy);
	tempMom += dummy->momentum();
	// Remove counted ones from the remaining list
	left.erase(remove(left.begin(),left.end(), dummy), left.end());
	}
	// Number of particles in this colour singlets so far
	int nparts = int(temp.size());
	// Insert the new particles in the reverse order
	sorted.insert(sorted.end(), temp.rbegin(), temp.rend());
	endPointMass += ((temp.back())->mass());
	// If it is a gluon loop we've already looped over the colour-singlet
	// in its entirety, so we need to end early
	if (gluonLoop){
	// Store the index of the final entry
	nparts += (_colSingletSize.empty()) ? 0 : _colSingletSize.back();
	// Insert the new particles in the correct order
	// i.e. triplet first, then the gluons we have seen so far
	Energy2 singletm2 = tempMom.m2();
	_colSingletSize.push_back(nparts);
	_colSingletm2.push_back(singletm2);
	continue;
	}
	// If it is not a gluon loop, we now need to trace the colourLine
	// down, until we reach the triplet.
	// Works similarly to the first half
	// Reset the temp PVector
	temp.resize(0);
	// Push back the particle in question
	temp.push_back(p);
	// tempMom hasn't been reset, add the particle we started with
	tempMom += p->momentum();
	left.erase(remove(left.begin(),left.end(), p), left.end());
	// Search down its colour line until we get to the end particle
	dummy = p->colourLine()->startParticle();
	temp.push_back(dummy);
	tempMom += dummy->momentum();
	left.erase(remove(left.begin(),left.end(), dummy), left.end());
	while ( dummy->hasColour() ){
	dummy = dummy->colourLine()->startParticle();
	temp.push_back(dummy);
	tempMom += dummy->momentum();
	left.erase(remove(left.begin(),left.end(), dummy), left.end());
	}
	endPointMass += ((temp.back())->mass());
	// Update size of colour singlet
	nparts += int(temp.size());
	nparts += (_colSingletSize.empty()) ? 0 : _colSingletSize.back();
	// Insert the new particles in the correct order
	Energy2 singletm2 = tempMom.m2();
	sorted.insert(sorted.end(), temp.begin(), temp.end());
	endPointMass += ((sorted.back())->mass());
	_colSingletSize.push_back(nparts);
	// Chooses whether to use invariant mass of singlet
	// or to use the lambda measure i.e. m^2 - (\sum m_i)^2
	Energy2 m2 = (_massMeasure == 0) ? singletm2 : singletm2 - sqr(endPointMass);
	_colSingletm2.push_back(m2);
	}
	// Else if it's a quark
	else if ( p->hasColour() ) {
	// Search up its colour line until we get to the start particle
	// Works the same way as the second half of the gluon handling
	tPPtr partner = p->colourLine()->startParticle();
	tPPtr dummy = partner;
	temp.push_back(p);
	endPointMass += ((temp.back())->mass());
	temp.push_back(partner);
	tempMom += p->momentum();
	tempMom += partner->momentum();
	left.erase(remove(left.begin(),left.end(), p), left.end());
	left.erase(remove(left.begin(),left.end(), partner), left.end());
	while ( dummy->hasColour() ){
	dummy = dummy->colourLine()->startParticle();
	temp.push_back(dummy);
	tempMom += dummy->momentum();
	left.erase(remove(left.begin(),left.end(), dummy), left.end());
	}
	// Number of particles in this colour singlets
	int nparts = int(temp.size());
	nparts += (_colSingletSize.empty()) ? 0 : _colSingletSize.back();
	// Insert the new particles in the correct order
	Energy2 singletm2 = tempMom.m2();
	sorted.insert(sorted.end(), temp.begin(), temp.end());
	endPointMass += ((sorted.back())->mass());
	_colSingletSize.push_back(nparts);
	Energy2 m2 = (_massMeasure == 0) ? singletm2 : singletm2 - sqr(endPointMass);
	_colSingletm2.push_back(m2);
	}
	// Else it's an antiquark
	else if ( p->hasAntiColour() ) {
	// Search along anti-colour line, storing particles, and reversing the order
	// at the end
	// Works in the same way as the first half of the gluon handling
	tPPtr partner = p->antiColourLine()->endParticle();
	tPPtr dummy = partner;
	temp.push_back(p);
	endPointMass += ((temp.back())->mass());
	temp.push_back(partner);
	tempMom += p->momentum();
	tempMom += partner->momentum();
	left.erase(remove(left.begin(),left.end(), p), left.end());
	left.erase(remove(left.begin(),left.end(), partner), left.end());
	while ( dummy->hasAntiColour() ){
	dummy = dummy->antiColourLine()->endParticle();
	temp.push_back(dummy);
	tempMom += dummy->momentum();
	left.erase(remove(left.begin(),left.end(), dummy), left.end());
	}
	// Number of particles in this colour singlets
	int nparts = int(temp.size());
	nparts += (_colSingletSize.empty()) ? 0 : _colSingletSize.back();
	// Insert the particles in the reverse order
	Energy2 singletm2 = tempMom.m2();
	sorted.insert(sorted.end(), temp.rbegin(), temp.rend());
	endPointMass += ((temp.back())->mass());
	_colSingletSize.push_back(nparts);
	Energy2 m2 = (_massMeasure == 0) ? singletm2 : singletm2 - sqr(endPointMass);
	_colSingletm2.push_back(m2);
	}
	}

	// Check that the algorithm hasn't missed any particles.
	assert( sorted.size() == tagged.size() );

	swap(sorted,tagged);

	}

	void PartonSplitter::massSplitTimeLikeGluon(tcPPtr ptrGluon,
	PPtr & ptrQ,
	PPtr & ptrQbar, size_t i){

	// only works in the Standard Model
	if ( ptrGluon->data().id() != ParticleID::g )
	throw Exception() << "strange enhancement only working with Standard Model hadronization"
	<< Exception::runerror;

	// select the quark flavour
	tPDPtr quark;
	long idNew=0;

	if ( ptrGluon->momentum().m() <
	2.0 *getParticle(ThePEG::ParticleID::s)->data().constituentMass() ) {
	throw Exception() << "Impossible Kinematics in PartonSplitter::massSplitTimeLikeGluon()"
	<< Exception::runerror;
	}
	// Only allow light quarks u,d,s with the probabilities
	double prob_d = _splitPwt[ParticleID::d];
	double prob_u = _splitPwt[ParticleID::u];
	double prob_s = enhanceStrange(i);
	int choice = UseRandom::rnd3(prob_u, prob_d, prob_s);
	switch(choice) {
	case 0: idNew = ParticleID::u; break;
	case 1: idNew = ParticleID::d; break;
	case 2: idNew = ParticleID::s; break;
	}
	ptrQ = getParticle(idNew);
	ptrQbar = getParticle(-idNew);

	// Solve the kinematics of the two body decay G --> Q + Qbar
	Lorentz5Momentum momentumQ;
	Lorentz5Momentum momentumQbar;
	- double cosThetaStar = UseRandom::rnd( -1.0 , 1.0 );
	- using Constants::pi;
	- double phiStar = UseRandom::rnd( -pi , pi );
	-
	- Energy constituentQmass;
	- constituentQmass = ptrQ->data().constituentMass();
	+ Energy constituentQmass = ptrQ->data().constituentMass();

	if (ptrGluon->momentum().m() < 2.0*constituentQmass) {
	throw Exception() << "Impossible Kinematics in PartonSplitter::massSplitTimeLikeGluon()"
	<< Exception::eventerror;
	}
	- Kinematics::twoBodyDecay(ptrGluon->momentum(), constituentQmass,
	- constituentQmass, cosThetaStar, phiStar, momentumQ,
	- momentumQbar );
	+
	+ switch (_kinematics)
	+ {
	+ case 0:
	+ {
	+ // Isotropic default
	+ double cosThetaStar = UseRandom::rnd( -1.0 , 1.0 );
	+ using Constants::pi;
	+ double phiStar = UseRandom::rnd( -pi , pi );
	+ Kinematics::twoBodyDecay(ptrGluon->momentum(), constituentQmass,
	+ constituentQmass, cosThetaStar, phiStar, momentumQ,
	+ momentumQbar );
	+ break;
	+ }
	+ case 1:
	+ {
	+ // Aligned
	+ Kinematics::twoBodyDecay(ptrGluon->momentum(), constituentQmass,
	+ constituentQmass, ptrGluon->momentum().vect().unit(), momentumQ,
	+ momentumQbar );
	+ break;
	+ }
	+ default:
	+ assert(false);
	+ }
	+
	// Create quark and anti-quark particles of the chosen flavour
	// and set they 5-momentum (the mass is the constituent one).
	ptrQ ->set5Momentum( momentumQ );
	ptrQbar ->set5Momentum( momentumQbar );

	}

	double PartonSplitter::enhanceStrange(size_t i){

	// Get the m2 of the relevant colour singlet
	// First we need to get the index of the colour-singlet the indexed i
	// parton has
	auto const it = lower_bound(_colSingletSize.begin(), _colSingletSize.end(), i);

	// Get the index of the colourSinglet mass
	int indx = distance(_colSingletSize.begin(), it);

	Energy2 mass2 = _colSingletm2[indx];
	Energy2 m2 = _m0*_m0;

	// Scaling strangeness enhancement
	if (_enhanceSProb == 1){
	// If the mass is significantly smaller than the characteristic mass,
	// just set the prob to 0
	double scale = double(m2/mass2);
	return (_maxScale < scale \|\| scale < 0.) ? 0. : pow(_splitPwt[ParticleID::s],scale);
	}
	// Exponential strangeness enhancement
	else if (_enhanceSProb == 2){
	double scale = double(m2/mass2);
	return (_maxScale < scale \|\| scale < 0.) ? 0. : exp(-scale);
	}
	else
	return _splitPwt[ParticleID::s];
	}

	diff --git a/Hadronization/PartonSplitter.h b/Hadronization/PartonSplitter.h
	--- a/Hadronization/PartonSplitter.h
	+++ b/Hadronization/PartonSplitter.h
	@@ -1,240 +1,245 @@
	// -- C++ --
	//
	// PartonSplitter.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_PartonSplitter_H
	#define HERWIG_PartonSplitter_H

	#include "CluHadConfig.h"
	#include <ThePEG/Interface/Interfaced.h>
	#include <ThePEG/Utilities/Selector.h>
	#include "PartonSplitter.fh"
	#include "HadronSpectrum.h"

	namespace Herwig {


	using namespace ThePEG;


	/** \ingroup Hadronization
	* \class PartonSplitter
	* \brief This class splits the gluons from the end of the shower.
	* \author Philip Stephens
	* \author Alberto Ribon
	*
	* This class does all of the nonperturbative parton splittings needed
	* immediately after the end of the showering (both initial and final),
	* as very first step of the cluster hadronization.
	*
	* the quarks are attributed with different weights for the splitting
	* by default only the splitting in u and d quarks is allowed
	* the option "set /Herwig/Hadronization/PartonSplitter:Split 1"
	* allows for additional splitting into s quarks based on some weight
	* in order for that to work the mass of the strange quark has to be changed
	* from the default value s.t. m_g > 2m_s
	*
	*
	* * @see \ref PartonSplitterInterfaces "The interfaces"
	* defined for PartonSplitter.
	*/
	class PartonSplitter: public Interfaced {

	public:

	/**
	* Default constructor
	*/
	PartonSplitter() :
	Interfaced(),
	_gluonDistance(ZERO),
	_splitGluon(0),
	+ _kinematics(0),
	_enhanceSProb(0),
	_m0(10.*GeV),
	_massMeasure(0) {}

	/**
	* This method does the nonperturbative splitting of:
	* time-like gluons. At the end of the shower the gluons should be
	* on a "physical" mass shell and should therefore be time-like.
	* @param tagged The tagged particles to be split
	* @return The particles which were not split and the products of splitting.
	*/
	void split(PVector & tagged);

	public:

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

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

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

	protected:

	/** @name 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;
	//@}

	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:

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

	/**
	* Non-perturbatively split a time-like gluon,
	* if something goes wrong null pointers are returned.
	* @param gluon The gluon to be split
	* @param quark The quark produced in the splitting
	* @param anti The antiquark produced in the splitting
	*/
	virtual void splitTimeLikeGluon(tcPPtr gluon, PPtr & quark, PPtr & anti);


	/**
	* Non-perturbatively split a time-like gluon using a mass-based
	* strangeness enhancement,
	* if something goes wrong null pointers are returned.
	* @param gluon The gluon to be split
	* @param quark The quark produced in the splitting
	* @param anti The antiquark produced in the splitting
	* @param gluon's index in the colour-sorted Particle vector
	*/
	void massSplitTimeLikeGluon(tcPPtr gluon, PPtr & quark, PPtr & anti,
	size_t i);

	/**
	* Colour-sort the event into grouped colour-singlets.
	* Convention: triplet - gluons - antitriplet, repeat for
	* all other colour-singlets or colourless particles.
	*/
	void colourSorted(PVector& tagged);

	/**
	* Method to calculate the probability for producing strange quarks
	*/
	double enhanceStrange(size_t i);

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

	private:

	/**
	* Different variables for the weight for splitting into various
	* light quark species.
	*/
	map<long,double> _splitPwt;

	/**
	* The selector to pick the type of quark
	*/
	Selector<PDPtr,double> _quarkSelector;

	/**
	* c tau for gluon decays
	*/
	Length _gluonDistance;

	/**
	* Flag used to determine between normal gluon splitting and alternative gluon splitting
	*/
	int _splitGluon;

	/**
	+ * Flag for choosing the Kinematics of the splitting
	+ */
	+ int _kinematics;
	+ /**
	* Vector that stores the index in the event record of the anti-triplet
	* of the colour-singlets, or of a colourless particle
	*/
	vector<int> _colSingletSize;

	/**
	* Vector that stores the masses of the colour-singlets or of a colourless
	* particle
	*/
	vector<Energy2> _colSingletm2;

	/**
	* Option to choose which functional form to use for strangeness
	* production
	*/
	int _enhanceSProb;

	/**
	* Characteristic mass scale for strangeness enhancement
	*/
	Energy _m0;

	/**
	* 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 that stops the scale in strangeness enhancement from
	* becoming too large
	*/
	const double _maxScale = 20.;

	/**
	* The hadron spectrum to consider
	*/
	Ptr<HadronSpectrum>::ptr _hadronSpectrum;

	};

	}

	#endif /* HERWIG_PartonSplitter_H */

File Metadata

Mime Type: text/x-diff
Expires: Sat, May 3, 5:52 AM (17 h, 17 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4982842
Default Alt Text: (29 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions