No OneTemporary
Actions

Size

33 KB

Subscribers

None

View Options

	diff --git a/Hadronization/StringFragmentation.cc b/Hadronization/StringFragmentation.cc
	--- a/Hadronization/StringFragmentation.cc
	+++ b/Hadronization/StringFragmentation.cc
	@@ -1,458 +1,673 @@
	// -- C++ --
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the StringFragmentation class.
	//

	#include "StringFragmentation.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Interface/ParVector.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/PDT/RemnantDecayer.h"
	#include "ThePEG/EventRecord/Particle.h"
	#include "ThePEG/EventRecord/Step.h"
	#include "ThePEG/Handlers/EventHandler.h"
	+#include "ThePEG/Repository/EventGenerator.h"
	+#include "ThePEG/Utilities/DebugItem.h"


	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"

	using namespace TheP8I;

	StringFragmentation::StringFragmentation()
	: pythia(), overlapN(0), overlapY(1.0), maxTries(2),
	#include "StringFragmentation-init.h"
	{}

	StringFragmentation::~StringFragmentation() {}

	void StringFragmentation::handle(EventHandler & eh, const tPVector & tagged,
	const Hint & h) {
	if ( !pythia.created() ) pythia.init(*this, theAdditionalP8Settings);
	if ( overlapN > 0 &&
	( overlapPythias.empty() \|\| !overlapPythias.back()->created() ) )
	createOverlapPythias();

	tPVector all = RemnantDecayer::decayRemnants(tagged, *newStep());
	vector<ColourSinglet> singlets;
	if ( theCollapser ) singlets = theCollapser->collapse(all, newStep());
	else singlets = ColourSinglet::getSinglets(all.begin(), all.end());


	Pythia8Interface * pytp = &pythia;

	if ( overlapN > 0 ) {
	int ncross = 0;
	for ( int i = 0, N = singlets.size(); i < N; ++i ) {
	if ( singlets[i].nPieces() > 1 ) {
	for ( int ip = 1, Np = singlets[i].nPieces(); ip <= Np; ++ip ) {
	double y = singlets[i].piece(ip)[0]->rapidity();
	int pos = y > 0.0? 1: -1;
	for ( int j = 1, M = singlets[i].piece(ip).size(); j < M; ++j ) {
	y = singlets[i].piece(ip)[j]->rapidity();
	if ( y > overlapY ) {
	if ( pos < 0 ) ++ncross;
	pos = 1;
	}
	if ( y < -overlapY ) {
	if ( pos > 0 ) ++ncross;
	pos = -1;
	}
	}
	}
	} else {
	double y = singlets[i].partons()[0]->rapidity();
	if ( singlets[i].partons()[0]->data().iColour() == PDT::Colour8 )
	y = singlets[i].partons().back()->rapidity();
	int pos = y > 0.0? 1: -1;
	for ( int j = 0, M = singlets[i].partons().size(); j < M; ++j ) {
	y = singlets[i].partons()[j]->rapidity();
	if ( y > overlapY ) {
	if ( pos < 0 ) ++ncross;
	pos = 1;
	}
	if ( y < -overlapY ) {
	if ( pos > 0 ) ++ncross;
	pos = -1;
	}
	}
	}
	}
	int ipyt = max(min(ncross, overlapN), 1) - 1;
	pytp = overlapPythias[ipyt];
	}

	+ for ( int i = 0, N = singlets.size(); i < N; ++i )
	+ hadronizeSystems(*pytp, vector<ColourSinglet>(1, singlets[i]), all);
	+ return;
	+
	Pythia8Interface & pyt = *pytp;
	Pythia8::Event & event = pyt.event();

	pyt.clearEvent();

	for ( int i = 0, N = singlets.size(); i < N; ++i ) {

	if ( singlets[i].nPieces() == 3 ) {
	// Simple junction.
	// Save place where we will store dummy particle.
	int nsave = event.size();
	event.append(22, -21, 0, 0, 0, 0, 0, 0, 0.0, 0.0, 0.0, 0.0);
	int npart = 0;
	for ( int ip = 1; ip <= 3; ++ip )
	for ( int j = 0, M = singlets[i].piece(ip).size(); j < M; ++j ) {
	pyt.addParticle(singlets[i].piece(ip)[j], 23, nsave, 0);
	++npart;
	}
	event[nsave].daughters(nsave + 1, nsave + npart);
	}
	else if ( singlets[i].nPieces() == 5 ) {
	// Double, connected junction
	// Save place where we will store dummy beam particles.
	int nb1 = event.size();
	int nb2 = nb1 + 1;
	event.append(2212, -21, 0, 0, nb1 + 2, nb1 + 4, 0, 0,
	0.0, 0.0, 0.0, 0.0);
	event.append(-2212, -21, 0, 0, nb2 + 4, nb2 + 6, 0, 0,
	0.0, 0.0, 0.0, 0.0);

	// Find the string piece connecting the junctions, and the
	// other loose string pieces.
	int connector = 0;
	int q1 = 0;
	int q2 = 0;
	int aq1 = 0;
	int aq2 = 0;
	for ( int ip = 1; ip <= 5; ++ip ) {
	if ( singlets[i].sink(ip).first && singlets[i].source(ip).first )
	connector = ip;
	else if ( singlets[i].source(ip).first ) {
	if ( q1 ) q2 = ip;
	else q1 = ip;
	}
	else if ( singlets[i].sink(ip).first ) {
	if ( aq1 ) aq2 = ip;
	else aq1 = ip;
	}
	}
	if ( !connector \|\| !q1 \|\| !q2 \|\| !aq1 \|\| ! aq2 )
	Throw<StringFragError>()
	<< name() << " found complicated junction string. Although Pytha8 can "
	<< "hadronize junction strings, this one was too complicated."
	<< Exception::runerror;

	// Insert the partons of the loose triplet ends.
	int start = event.size();
	for ( int j = 0, M = singlets[i].piece(q1).size(); j < M; ++j )
	pyt.addParticle(singlets[i].piece(q1)[j], 23, nb1, 0);
	for ( int j = 0, M = singlets[i].piece(q2).size(); j < M; ++j )
	pyt.addParticle(singlets[i].piece(q2)[j], 23, nb1, 0);
	// Insert dummy triplet incoming parton with correct colour code.
	int col1 = singlets[i].piece(connector).empty()? event.nextColTag():
	pyt.addColourLine(singlets[i].piece(connector).front()->colourLine());
	int dum1 = event.size();
	event.append(2, -21, nb1, 0, 0, 0, col1, 0, 0.0, 0.0, 0.0, 0.0, 0.0);
	event[nb1].daughters(start, start + singlets[i].piece(q1).size() +
	singlets[i].piece(q2).size());

	// Insert the partons of the loose anti-triplet ends.
	start = event.size();
	for ( int j = 0, M = singlets[i].piece(aq1).size(); j < M; ++j )
	pyt.addParticle(singlets[i].piece(aq1)[j], 23, nb2, 0);
	for ( int j = 0, M = singlets[i].piece(aq2).size(); j < M; ++j )
	pyt.addParticle(singlets[i].piece(aq2)[j], 23, nb2, 0);
	// Insert dummy anti-triplet incoming parton with correct colour code.
	int col2 = singlets[i].piece(connector).empty()? col1:
	pyt.addColourLine(singlets[i].piece(connector).back()->antiColourLine());
	int dum2 = event.size();
	event.append(-2, -21, nb2, 0, 0, 0, 0, col2, 0.0, 0.0, 0.0, 0.0, 0.0);
	event[nb2].daughters(start, start + singlets[i].piece(aq1).size() +
	singlets[i].piece(aq2).size());

	// Add the partons from the connecting string piece.
	for ( int j = 0, M = singlets[i].piece(connector).size(); j < M; ++j )
	pyt.addParticle(singlets[i].piece(connector)[j], 23, dum1, dum2);
	}
	else if ( singlets[i].nPieces() > 1 ) {
	// We don't know how to handle other junctions yet.
	Throw<StringFragError>()
	<< name() << " found complicated junction string. Although Pytha8 can "
	<< "hadronize junction strings, that interface is not ready yet."
	<< Exception::runerror;
	} else {
	// Normal string
	for ( int j = 0, M = singlets[i].partons().size(); j < M; ++j )
	pyt.addParticle(singlets[i].partons()[j], 23, 0, 0);
	}

	}

	for ( int i = 0, N = all.size(); i < N; ++i )
	if ( !all[i]->coloured() )
	pyt.addParticle(all[i], 23, 0, 0);

	int oldsize = event.size();

	Pythia8::Event saveEvent = event;
	int ntry = maxTries;
	CurrentGenerator::Redirect stdout(cout, false);
	while ( !pyt.go() && --ntry ) event = saveEvent;
	if ( !ntry ) {
	CurrentGenerator::Redirect stdout2(cout, true);
	event.list();
	cout << "ThePEG event listing:\n" << *eh.currentEvent();
	Throw<StringFragError>()
	<< "Pythia8 failed to hadronize partonic state:\n"
	<< stdout2.str() << "This event will be discarded!\n"
	<< Exception::warning;
	throw Veto();
	}
	// event.list(cerr);
	map<tPPtr, set<tPPtr> > children;
	map<tPPtr, set<tPPtr> > parents;
	for ( int i = 1, N = event.size(); i < N; ++i ) {
	tPPtr p = pyt.getParticle(i);
	int d1 = event[i].daughter1();
	if ( d1 <= 0 ) continue;
	children[p].insert(pyt.getParticle(d1));
	parents[pyt.getParticle(d1)].insert(p);
	int d2 = event[i].daughter2();
	if ( d2 > 0 ) {
	children[p].insert(pyt.getParticle(d2));
	parents[pyt.getParticle(d2)].insert(p);
	}
	for ( int di = d1 + 1; di < d2; ++di ) {
	children[p].insert(pyt.getParticle(di));
	parents[pyt.getParticle(di)].insert(p);
	}
	}

	for ( int i = oldsize, N = event.size(); i < N; ++i ) {
	PPtr p = pyt.getParticle(i);
	set<tPPtr> & pars = parents[p];
	if ( !p ) {
	Throw<StringFragError>()
	<< "Failed to reconstruct hadronized state from Pythia8:\n"
	<< stdout.str() << "This event will be discarded!\n" << Exception::warning;
	throw Veto();
	}
	if ( pars.empty() ) {
	Pythia8::Particle & pyp = event[i];
	if ( pyp.mother1() > 0 ) pars.insert(pyt.getParticle(pyp.mother1()));
	if ( pyp.mother2() > 0 ) pars.insert(pyt.getParticle(pyp.mother2()));
	for ( int im = pyp.mother1() + 1; im < pyp.mother2(); ++im )
	pars.insert(pyt.getParticle(im));
	if ( pars.empty() ) {
	Throw<StringFragError>()
	<< "Failed to reconstruct hadronized state from Pythia8:\n"
	<< stdout.str() << "This event will be discarded!\n" << Exception::warning;
	throw Veto();
	}
	}
	if ( pars.size() == 1 ) {
	tPPtr par = *pars.begin();
	if ( children[par].size() == 1 &&
	*children[par].begin() == p && par->id() == p->id() )
	newStep()->setCopy(par, p);
	else
	newStep()->addDecayProduct(par, p);
	} else {
	newStep()->addDecayProduct(pars.begin(), pars.end(), p);
	}
	}

	// for ( int i = oldsize, N = event.size(); i < N; ++i ) {
	// PPtr p = pyt.getParticle(i);
	// if ( !p ) continue;
	// Pythia8::Particle & pyp = event[i];
	// if ( pyp.mother1() <= 0 ) continue;
	// if ( pyp.mother1() == pyp.mother2() ) {
	// if ( !newStep()->setCopy(pyt.getParticle(pyp.mother1()), p) )
	// cerr << i << " " << pyp.mother1() << " " << pyp.mother2() << endl;
	// }
	// else if ( pyp.mother2() <= 0 ) {
	// if ( !newStep()->addDecayProduct(pyt.getParticle(pyp.mother1()), p) )
	// cerr << i << " " << pyp.mother1() << " " << pyp.mother2() << endl;
	// }
	// else if ( pyp.mother2() > pyp.mother1() ) {
	// tPVector parents;
	// for ( int j = pyp.mother1(); j <= pyp.mother2(); ++j ) {
	// tPPtr par = pyt.getParticle(j);
	// if ( par ) parents.push_back(par);
	// }
	// if ( !newStep()->addDecayProduct(parents.begin(), parents.end(), p) )
	// cerr << i << " " << pyp.mother1() << " " << pyp.mother2() << endl;
	// }
	// }
	}

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

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


	void StringFragmentation::doinitrun() {
	HadronizationHandler::doinitrun();
	theAdditionalP8Settings.push_back("ProcessLevel:all = off");
	theAdditionalP8Settings.push_back("HadronLevel:Decay = off");
	theAdditionalP8Settings.push_back("Check:event = off");
	pythia.init(*this, theAdditionalP8Settings);
	if ( overlapN > 0 ) createOverlapPythias();
	}

	void StringFragmentation::createOverlapPythias() {
	overlapPythias.clear();
	overlapPythias.resize(overlapN);
	if ( overlapSigma.empty() ) overlapSigma.push_back(0.23723);
	overlapSigma.resize(overlapN, overlapSigma.back());
	if ( overlapA.empty() ) overlapA.push_back(0.3);
	overlapA.resize(overlapN, overlapA.back());
	if ( overlapB.empty() ) overlapB.push_back(0.3);
	overlapB.resize(overlapN, overlapB.back());
	if ( overlapProbStoUD.empty() ) overlapProbStoUD.push_back(0.19);
	overlapProbStoUD.resize(overlapN, overlapProbStoUD.back());
	if ( overlapProbSQtoQQ.empty() ) overlapProbSQtoQQ.push_back(1.0);
	overlapProbSQtoQQ.resize(overlapN, overlapProbSQtoQQ.back());
	if ( overlapProbQQ1toQQ0.empty() ) overlapProbQQ1toQQ0.push_back(0.027);
	overlapProbQQ1toQQ0.resize(overlapN, overlapProbQQ1toQQ0.back());
	if ( overlapProbQQtoQ.empty() ) overlapProbQQtoQ.push_back(0.09);
	overlapProbQQtoQ.resize(overlapN, overlapProbQQtoQ.back());

	for ( int ipyt = 0; ipyt < overlapN; ++ipyt ) {
	vector<string> moresettings = theAdditionalP8Settings;
	moresettings.push_back("StringPT:sigma = " + convert(overlapSigma[ipyt]));
	moresettings.push_back("StringZ:aLund = " + convert(overlapA[ipyt]));
	moresettings.push_back("StringZ:bLund = " + convert(overlapB[ipyt]));
	moresettings.push_back("StringFlav:probStoUD = " +
	convert(overlapProbStoUD[ipyt]));
	moresettings.push_back("StringFlav:probSQtoQQ = " +
	convert(overlapProbSQtoQQ[ipyt]));
	moresettings.push_back("StringFlav:probQQ1toQQ0 = " +
	convert(overlapProbQQ1toQQ0[ipyt]));
	moresettings.push_back("StringFlav:probQQtoQ = " +
	convert(overlapProbQQtoQ[ipyt]));
	overlapPythias[ipyt] = new_ptr(Pythia8Interface());
	overlapPythias[ipyt]->init(*this, moresettings);
	}
	}

	string StringFragmentation::convert(double d) {
	ostringstream os;
	os << d;
	return os.str();
	}

	void StringFragmentation::persistentOutput(PersistentOStream & os) const {
	os
	#include "StringFragmentation-output.h"
	<< overlapN << overlapY << overlapSigma << overlapA << overlapB
	<< overlapProbStoUD << overlapProbSQtoQQ << overlapProbQQ1toQQ0
	<< overlapProbQQtoQ << maxTries << theCollapser;
	}

	void StringFragmentation::persistentInput(PersistentIStream & is, int) {
	is
	#include "StringFragmentation-input.h"
	>> overlapN >> overlapY >> overlapSigma >> overlapA >> overlapB
	>> overlapProbStoUD >> overlapProbSQtoQQ >> overlapProbQQ1toQQ0
	>> overlapProbQQtoQ >> maxTries >> theCollapser;
	overlapPythias.clear();
	}

	ClassDescription<StringFragmentation> StringFragmentation::initStringFragmentation;
	// Definition of the static class description member.

	void StringFragmentation::Init() {

	#include "StringFragmentation-interfaces.h"


	static Reference<StringFragmentation,ClusterCollapser> interfaceCollapser
	("Collapser",
	"A ThePEG::ClusterCollapser object used to collapse colour singlet "
	"clusters which are too small to fragment. If no object is given the "
	"MinistringFragmentetion of Pythia8 is used instead.",
	&StringFragmentation::theCollapser, true, false, true, true, false);

	static Parameter<StringFragmentation,int> interfaceOverlapN
	("OverlapN",
	"If positive, use the string overlap model with different "
	"parameters for up to overlapN strings.",
	&StringFragmentation::overlapN, 0, 0, 0,
	true, false, Interface::lowerlim);

	static Parameter<StringFragmentation,double> interfaceOverlapY
	("OverlapY",
	"In the string overlap model the number of strings spanning "
	"the mid-rapidity bin with width overlapY selects which "
	"hadronization parameters to use.",
	&StringFragmentation::overlapY, 1.0, 0.0, 0,
	true, false, Interface::lowerlim);

	static ParVector<StringFragmentation,double> interfaceOverlapSigma
	("OverlapSigma",
	"Different sigmas used for the string overlap model. The first "
	"entry is used for the case of no overlapping strings, the second "
	"for two overlapping strings, and so on.",
	&StringFragmentation::overlapSigma, -1, 0.23723, 0.0, 0,
	true, false, Interface::lowerlim);

	static ParVector<StringFragmentation,double> interfaceOverlapA
	("OverlapA",
	"Different a-parameters used for the string overlap model. The first "
	"entry is used for the case of no overlapping strings, the second "
	"for two overlapping strings, and so on.",
	&StringFragmentation::overlapA, -1, 0.3, 0.0, 0,
	true, false, Interface::lowerlim);

	static ParVector<StringFragmentation,double> interfaceOverlapB
	("OverlapB",
	"Different b-parameters used for the string overlap model. The first "
	"entry is used for the case of no overlapping strings, the second "
	"for two overlapping strings, and so on.",
	&StringFragmentation::overlapB, -1, 0.8, 0.0, 0,
	true, false, Interface::lowerlim);

	static ParVector<StringFragmentation,double> interfaceOverlapProbStoUD
	("OverlapProbStoUD",
	"Different s/ud suppression parameters used for the string overlap "
	"model. The first entry is used for the case of no overlapping "
	"strings, the second for two overlapping strings, and so on.",
	&StringFragmentation::overlapProbStoUD, -1, 0.19, 0.0, 0,
	true, false, Interface::lowerlim);

	static ParVector<StringFragmentation,double> interfaceOverlapProbSQtoQQ
	("OverlapProbSQtoQQ",
	"Different sq/qq suppression parameters used for the string overlap "
	"model. The first entry is used for the case of no overlapping "
	"strings, the second for two overlapping strings, and so on.",
	&StringFragmentation::overlapProbSQtoQQ, -1, 1.0, 0.0, 0,
	true, false, Interface::lowerlim);

	static ParVector<StringFragmentation,double> interfaceOverlapProbQQ1toQQ0
	("OverlapProbQQ1toQQ0",
	"Different qq_1/qq_0 suppression parameters used for the string overlap "
	"model. The first entry is used for the case of no overlapping "
	"strings, the second for two overlapping strings, and so on.",
	&StringFragmentation::overlapProbQQ1toQQ0, -1, 0.027, 0.0, 0,
	true, false, Interface::lowerlim);

	static ParVector<StringFragmentation,double> interfaceOverlapProbQQtoQ
	("OverlapProbQQtoQ",
	"Different s/ud suppression parameters used for the string overlap "
	"model. The first entry is used for the case of no overlapping "
	"strings, the second for two overlapping strings, and so on.",
	&StringFragmentation::overlapProbQQtoQ, -1, 0.09, 0.0, 0,
	true, false, Interface::lowerlim);

	static Parameter<StringFragmentation,int> interfaceMaxTries
	("MaxTries",
	"Sometimes Pythia gives up on an event too easily. We therefore allow "
	"it to re-try a couple of times.",
	&StringFragmentation::maxTries, 2, 1, 0,
	true, false, Interface::lowerlim);

	}

	+
	+void StringFragmentation::
	+hadronizeSystems(Pythia8Interface & pyt, const vector<ColourSinglet> & singlets, const tPVector & all) {
	+
	+ Pythia8::Event & event = pyt.event();
	+ pyt.clearEvent();
	+
	+ for ( int i = 0, N = singlets.size(); i < N; ++i ) {
	+
	+ if ( singlets[i].nPieces() == 3 ) {
	+ // Simple junction.
	+ // Save place where we will store dummy particle.
	+ int nsave = event.size();
	+ event.append(22, -21, 0, 0, 0, 0, 0, 0, 0.0, 0.0, 0.0, 0.0);
	+ int npart = 0;
	+ for ( int ip = 1; ip <= 3; ++ip )
	+ for ( int j = 0, M = singlets[i].piece(ip).size(); j < M; ++j ) {
	+ pyt.addParticle(singlets[i].piece(ip)[j], 23, nsave, 0);
	+ ++npart;
	+ }
	+ event[nsave].daughters(nsave + 1, nsave + npart);
	+ }
	+ else if ( singlets[i].nPieces() == 5 ) {
	+ // Double, connected junction
	+ // Save place where we will store dummy beam particles.
	+ int nb1 = event.size();
	+ int nb2 = nb1 + 1;
	+ event.append(2212, -21, 0, 0, nb1 + 2, nb1 + 4, 0, 0,
	+ 0.0, 0.0, 0.0, 0.0);
	+ event.append(-2212, -21, 0, 0, nb2 + 4, nb2 + 6, 0, 0,
	+ 0.0, 0.0, 0.0, 0.0);
	+
	+ // Find the string piece connecting the junctions, and the
	+ // other loose string pieces.
	+ int connector = 0;
	+ int q1 = 0;
	+ int q2 = 0;
	+ int aq1 = 0;
	+ int aq2 = 0;
	+ for ( int ip = 1; ip <= 5; ++ip ) {
	+ if ( singlets[i].sink(ip).first && singlets[i].source(ip).first )
	+ connector = ip;
	+ else if ( singlets[i].source(ip).first ) {
	+ if ( q1 ) q2 = ip;
	+ else q1 = ip;
	+ }
	+ else if ( singlets[i].sink(ip).first ) {
	+ if ( aq1 ) aq2 = ip;
	+ else aq1 = ip;
	+ }
	+ }
	+ if ( !connector \|\| !q1 \|\| !q2 \|\| !aq1 \|\| ! aq2 )
	+ Throw<StringFragError>()
	+ << name() << " found complicated junction string. Although Pytha8 can "
	+ << "hadronize junction strings, this one was too complicated."
	+ << Exception::runerror;
	+
	+ // Insert the partons of the loose triplet ends.
	+ int start = event.size();
	+ for ( int j = 0, M = singlets[i].piece(q1).size(); j < M; ++j )
	+ pyt.addParticle(singlets[i].piece(q1)[j], 23, nb1, 0);
	+ for ( int j = 0, M = singlets[i].piece(q2).size(); j < M; ++j )
	+ pyt.addParticle(singlets[i].piece(q2)[j], 23, nb1, 0);
	+ // Insert dummy triplet incoming parton with correct colour code.
	+ int col1 = singlets[i].piece(connector).empty()? event.nextColTag():
	+ pyt.addColourLine(singlets[i].piece(connector).front()->colourLine());
	+ int dum1 = event.size();
	+ event.append(2, -21, nb1, 0, 0, 0, col1, 0, 0.0, 0.0, 0.0, 0.0, 0.0);
	+ event[nb1].daughters(start, start + singlets[i].piece(q1).size() +
	+ singlets[i].piece(q2).size());
	+
	+ // Insert the partons of the loose anti-triplet ends.
	+ start = event.size();
	+ for ( int j = 0, M = singlets[i].piece(aq1).size(); j < M; ++j )
	+ pyt.addParticle(singlets[i].piece(aq1)[j], 23, nb2, 0);
	+ for ( int j = 0, M = singlets[i].piece(aq2).size(); j < M; ++j )
	+ pyt.addParticle(singlets[i].piece(aq2)[j], 23, nb2, 0);
	+ // Insert dummy anti-triplet incoming parton with correct colour code.
	+ int col2 = singlets[i].piece(connector).empty()? col1:
	+ pyt.addColourLine(singlets[i].piece(connector).back()->antiColourLine());
	+ int dum2 = event.size();
	+ event.append(-2, -21, nb2, 0, 0, 0, 0, col2, 0.0, 0.0, 0.0, 0.0, 0.0);
	+ event[nb2].daughters(start, start + singlets[i].piece(aq1).size() +
	+ singlets[i].piece(aq2).size());
	+
	+ // Add the partons from the connecting string piece.
	+ for ( int j = 0, M = singlets[i].piece(connector).size(); j < M; ++j )
	+ pyt.addParticle(singlets[i].piece(connector)[j], 23, dum1, dum2);
	+ }
	+ else if ( singlets[i].nPieces() > 1 ) {
	+ // We don't know how to handle other junctions yet.
	+ Throw<StringFragError>()
	+ << name() << " found complicated junction string. Although Pytha8 can "
	+ << "hadronize junction strings, that interface is not ready yet."
	+ << Exception::runerror;
	+ } else {
	+ // Normal string
	+ for ( int j = 0, M = singlets[i].partons().size(); j < M; ++j )
	+ pyt.addParticle(singlets[i].partons()[j], 23, 0, 0);
	+ }
	+
	+ }
	+
	+ for ( int i = 0, N = all.size(); i < N; ++i )
	+ if ( !all[i]->coloured() )
	+ pyt.addParticle(all[i], 23, 0, 0);
	+
	+ int oldsize = event.size();
	+
	+ Pythia8::Event saveEvent = event;
	+ int ntry = maxTries;
	+ CurrentGenerator::Redirect stdout(cout, false);
	+ while ( !pyt.go() && --ntry ) event = saveEvent;
	+ if ( !ntry ) {
	+ static DebugItem printfailed("TheP8I::PrintFailed", 10);
	+ if ( printfailed ) {
	+ double ymax = -1000.0;
	+ double ymin = 1000.0;
	+ double sumdy = 0.0;
	+ for ( int i = 0, N = singlets.size(); i < N; ++i )
	+ for ( int j = 0, M = singlets[i].partons().size(); j < M; ++j ) {
	+ const Particle & p = *singlets[i].partons()[j];
	+ ymax = max(ymax, p.momentum().rapidity());
	+ ymin = min(ymin, p.momentum().rapidity());
	+ cerr << setw(5) << j << setw(14) << p.momentum().rapidity()
	+ << setw(14) << p.momentum().perp()/GeV
	+ << setw(14) << p.momentum().m()/GeV;
	+ if ( j == 0 && p.data().iColour() == PDT::Colour8 ) {
	+ cerr << setw(14) << (p.momentum() + singlets[i].partons().back()->momentum()).m()/GeV;
	+ sumdy += abs(p.momentum().rapidity() - singlets[i].partons().back()->momentum().rapidity());
	+ }
	+ else if ( j > 0 ) {
	+ cerr << setw(14) << (p.momentum() + singlets[i].partons()[j-1]->momentum()).m()/GeV;
	+ sumdy += abs(p.momentum().rapidity() - singlets[i].partons()[j-1]->momentum().rapidity());
	+ if ( j + 1 < M )
	+ cerr << setw(14) << (p.momentum() + singlets[i].partons()[j-1]->momentum()).m()*
	+ (p.momentum() + singlets[i].partons()[j+1]->momentum()).m()/
	+ (p.momentum() + singlets[i].partons()[j+1]->momentum() +
	+ singlets[i].partons()[j-1]->momentum()).m()/GeV;
	+ }
	+ cerr << endl;
	+ }
	+ cerr << setw(14) << ymax - ymin << setw(14) << sumdy/(ymax - ymin) << endl;
	+ }
	+ CurrentGenerator::Redirect stdout2(cout, true);
	+ event.list();
	+ cout << "ThePEG event listing:\n" << *(generator()->currentEvent());
	+ Throw<StringFragError>()
	+ << "Pythia8 failed to hadronize partonic state:\n"
	+ << stdout2.str() << "This event will be discarded!\n"
	+ << Exception::warning;
	+ throw Veto();
	+ }
	+ // event.list(cerr);
	+ map<tPPtr, set<tPPtr> > children;
	+ map<tPPtr, set<tPPtr> > parents;
	+ for ( int i = 1, N = event.size(); i < N; ++i ) {
	+ tPPtr p = pyt.getParticle(i);
	+ int d1 = event[i].daughter1();
	+ if ( d1 <= 0 ) continue;
	+ children[p].insert(pyt.getParticle(d1));
	+ parents[pyt.getParticle(d1)].insert(p);
	+ int d2 = event[i].daughter2();
	+ if ( d2 > 0 ) {
	+ children[p].insert(pyt.getParticle(d2));
	+ parents[pyt.getParticle(d2)].insert(p);
	+ }
	+ for ( int di = d1 + 1; di < d2; ++di ) {
	+ children[p].insert(pyt.getParticle(di));
	+ parents[pyt.getParticle(di)].insert(p);
	+ }
	+ }
	+
	+ for ( int i = oldsize, N = event.size(); i < N; ++i ) {
	+ PPtr p = pyt.getParticle(i);
	+ set<tPPtr> & pars = parents[p];
	+ if ( !p ) {
	+ Throw<StringFragError>()
	+ << "Failed to reconstruct hadronized state from Pythia8:\n"
	+ << stdout.str() << "This event will be discarded!\n" << Exception::warning;
	+ throw Veto();
	+ }
	+ if ( pars.empty() ) {
	+ Pythia8::Particle & pyp = event[i];
	+ if ( pyp.mother1() > 0 ) pars.insert(pyt.getParticle(pyp.mother1()));
	+ if ( pyp.mother2() > 0 ) pars.insert(pyt.getParticle(pyp.mother2()));
	+ for ( int im = pyp.mother1() + 1; im < pyp.mother2(); ++im )
	+ pars.insert(pyt.getParticle(im));
	+ if ( pars.empty() ) {
	+ Throw<StringFragError>()
	+ << "Failed to reconstruct hadronized state from Pythia8:\n"
	+ << stdout.str() << "This event will be discarded!\n" << Exception::warning;
	+ throw Veto();
	+ }
	+ }
	+ if ( pars.size() == 1 ) {
	+ tPPtr par = *pars.begin();
	+ if ( children[par].size() == 1 &&
	+ *children[par].begin() == p && par->id() == p->id() )
	+ newStep()->setCopy(par, p);
	+ else
	+ newStep()->addDecayProduct(par, p);
	+ } else {
	+ newStep()->addDecayProduct(pars.begin(), pars.end(), p);
	+ }
	+ }
	+
	+}
	+
	diff --git a/Hadronization/StringFragmentation.h b/Hadronization/StringFragmentation.h
	--- a/Hadronization/StringFragmentation.h
	+++ b/Hadronization/StringFragmentation.h
	@@ -1,261 +1,269 @@
	// -- C++ --
	#ifndef THEP8I_StringFragmentation_H
	#define THEP8I_StringFragmentation_H
	//
	// This is the declaration of the StringFragmentation class.
	//

	#include "ThePEG/Handlers/HadronizationHandler.h"
	#include "ThePEG/Handlers/ClusterCollapser.h"
	#include "TheP8I/Config/Pythia8Interface.h"

	namespace TheP8I {

	using namespace ThePEG;

	/**
	* Here is the documentation of the StringFragmentation class.
	*
	* @see \ref StringFragmentationInterfaces "The interfaces"
	* defined for StringFragmentation.
	*/
	class StringFragmentation: public HadronizationHandler {

	public:

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* The default constructor.
	*/
	StringFragmentation();

	/**
	* The destructor.
	*/
	virtual ~StringFragmentation();
	//@}

	public:

	/**
	* The main function called by the EventHandler class to
	* perform the Hadronization step.
	* @param eh the EventHandler in charge of the Event generation.
	* @param tagged if not empty these are the only particles which should
	* be considered by the StepHandler.
	* @param hint a Hint object with possible information from previously
	* performed steps.
	* @throws Veto if the StepHandler requires the current step to be
	* discarded.
	* @throws Stop if the generation of the current Event should be stopped
	* after this call.
	* @throws Exception if something goes wrong.
	*/
	virtual void handle(EventHandler & eh, const tPVector & tagged,
	const Hint & hint);

	+
	+ /**
	+ * Let the given Pythia8Interface hadronize the ColourSinglet
	+ * systems provided.
	+ */
	+ void hadronizeSystems(Pythia8Interface & pyt, const vector<ColourSinglet> & singlets,
	+ const tPVector & all);
	+
	/** @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. Called in the run phase just before
	* a run begins.
	*/
	virtual void doinitrun();
	//@}

	private:

	/**
	* The interface to the Pythia 8 object,
	*/
	Pythia8Interface pythia;

	/**
	* If positive, use the string overlap model with different
	* parameters for up to overlapN strings.
	*/
	int overlapN;

	/**
	* In the string overlap model the number of strings spanning the
	* mid-rapidity bin with width overlapY selects which hadronization
	* parameters to use.
	*/
	double overlapY;

	/**
	* Additional interfaces to Pythia8 objects in case the string
	* overlap model is to be used.
	*/
	vector<Ptr<Pythia8Interface>::ptr> overlapPythias;

	/**
	* Different sigmas used for the string overlap model.
	*/
	vector<double> overlapSigma;

	/**
	* Different a-parameters used for the string overlap model.
	*/
	vector<double> overlapA;

	/**
	* Different a-parameters used for the string overlap model.
	*/
	vector<double> overlapB;

	/**
	* Different s/ud suppression parameters used for the string overlap
	* model.
	*/
	vector<double> overlapProbStoUD;

	/**
	* Different sq/qq suppression parameters used for the string overlap
	* model.
	*/
	vector<double> overlapProbSQtoQQ;

	/**
	* Different qq_1/qq_0 suppression parameters used for the string overlap
	* model.
	*/
	vector<double> overlapProbQQ1toQQ0;

	/**
	* Different q/qq suppression parameters used for the string overlap
	* model.
	*/
	vector<double> overlapProbQQtoQ;

	/**
	* Sometimes Pythia gives up on an event too easily. We allow it to
	* re-try a couple of times.
	*/
	int maxTries;

	/**
	* Little helper function which probably should go into
	* ThePEG::StringUtils.
	*/
	static string convert(double d);

	/**
	* Create and initialize Pythia8 objects for overlapping strings.
	*/
	void createOverlapPythias();

	#include "StringFragmentation-var.h"

	/**
	* The object used to avoid too small strings in the hadronization.
	*/
	ClusterCollapserPtr theCollapser;

	public:

	/**
	* Exception class declaration.
	*/
	struct StringFragError: public Exception {};

	private:

	/**
	* The static object used to initialize the description of this class.
	* Indicates that this is a concrete class with persistent data.
	*/
	static ClassDescription<StringFragmentation> initStringFragmentation;

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

	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

	/** This template specialization informs ThePEG about the
	* base classes of StringFragmentation. */
	template <>
	struct BaseClassTrait<TheP8I::StringFragmentation,1> {
	/** Typedef of the first base class of StringFragmentation. */
	typedef HadronizationHandler NthBase;
	};

	/** This template specialization informs ThePEG about the name of
	* the StringFragmentation class and the shared object where it is defined. */
	template <>
	struct ClassTraits<TheP8I::StringFragmentation>
	: public ClassTraitsBase<TheP8I::StringFragmentation> {
	/** Return a platform-independent class name */
	static string className() { return "TheP8I::StringFragmentation"; }
	/**
	* The name of a file containing the dynamic library where the class
	* StringFragmentation is implemented. It may also include several, space-separated,
	* libraries if the class StringFragmentation depends on other classes (base classes
	* excepted). In this case the listed libraries will be dynamically
	* linked in the order they are specified.
	*/
	static string library() { return "libTheP8I.so"; }
	};

	/** @endcond */

	}

	#endif /* THEP8I_StringFragmentation_H */

File Metadata

Mime Type: text/x-diff
Expires: Thu, Apr 24, 6:38 AM (1 d, 21 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4880421
Default Alt Text: (33 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions