No OneTemporary
Actions

Size

31 KB

Subscribers

None

View Options

	diff --git a/BE/BoseEinsteinHandler.cc b/BE/BoseEinsteinHandler.cc
	--- a/BE/BoseEinsteinHandler.cc
	+++ b/BE/BoseEinsteinHandler.cc
	@@ -1,168 +1,182 @@
	// -- C++ --
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the BoseEinsteinHandler class.
	//

	#include "BoseEinsteinHandler.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/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"

	using namespace TheP8I;

	BoseEinsteinHandler::BoseEinsteinHandler()
	: pythia(), doContinueDecays(true),
	#include "BoseEinsteinHandler-init.h"
	{}

	BoseEinsteinHandler::~BoseEinsteinHandler() {}

	void BoseEinsteinHandler::
	performDecay(tPPtr parent, Step & s) const {
	if ( parent->data().width() < theBoseEinstein_widthSep*GeV &&
	parent->id() != ParticleID::K0 &&
	parent->id() != ParticleID::Kbar0 ) return;
	ParticleVector children = Decayer::DecayParticle(parent, s, maxLoop());
	for ( int i = 0, N = children.size(); i < N; ++i )
	if ( !children[i]->data().stable() ) performDecay(children[i], s);
	}

	void BoseEinsteinHandler::handle(EventHandler & eh, const tPVector & tagged,
	const Hint & h) {
	if ( !pythia.created() ) {
	pythia.init(*this, theAdditionalP8Settings);
	theBE.init(&(pythia().info), pythia().settings, pythia().particleData);
	}

	// First go through to see which of the tagged particles should be
	// decayed and/or shifted: Exit if none are found.
	tPVector parents;
	for ( int i = 0, N = tagged.size(); i < N; ++i )
	if ( tagged[i] ) parents.push_back(tagged[i]);

	if ( parents.empty() ) return;

	// Create a new step, decay all particles which have a width larger
	// than the given limit
	// (<interface>BoseEinstein_widthSep</interface>) and add their
	// children in the new step.
	for ( int i = 0, N = parents.size(); i < N; ++i )
	if ( !parents[i]->data().stable() )
	performDecay(newStep()->find(parents[i]->final()), *newStep());

	// Put all particles (or their children fi they have decayed) in the
	// Pythia8::Event.
	pythia.clearEvent();
	particles.clear();
	for ( int i = 0, N = parents.size(); i < N; ++i ) addParticle(parents[i]);

	// Perform the BE-shift
	theBE.shiftEvent(pythia.event());

	// Insert all shifted particles in the Step.
	for ( int i = 0, N = pythia.event().size(); i < N; ++i ) {
	if ( !pythia.event()[i].isFinal() ) continue;
	int im = pythia.event()[i].mother1();
	if ( im <= 0 ) continue;
	Lorentz5Momentum p(pythia.event()[i].px()*GeV,
	pythia.event()[i].py()*GeV,
	pythia.event()[i].pz()*GeV,
	pythia.event()[i].e()*GeV,
	pythia.event()[i].m()*GeV);
	if ( particles[im]->momentum().x() != p.x() \|\|
	particles[im]->momentum().y() != p.y() \|\|
	particles[im]->momentum().z() != p.z() \|\|
	particles[im]->momentum().e() != p.e() ) {
	tPPtr newp = newStep()->copyParticle(particles[im]);
	newp->set5Momentum(p);
	particles[im] = newp;
	}
	}

	// Now go through all particles and decay them further if requested.
	if ( !doContinueDecays ) return;
	for ( int i = 0, N = particles.size(); i < N; ++i )
	if ( particles[i] && !particles[i]->data().stable() )
	DecayHandler::performDecay(particles[i], *newStep());

	}

	void BoseEinsteinHandler::addParticle(tPPtr p) {
	if ( p->next() ) addParticle(p->next());
	else if ( p->decayed() )
	for ( int i = 0, N = p->children().size(); i < N; ++i )
	addParticle(p->children()[i]);
	else {
	int idx = pythia.addParticle(p, 1, 0, 0);
	particles.resize(idx + 1);
	particles[idx] = p;
	}
	}

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

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


	void BoseEinsteinHandler::doinitrun() {
	DecayHandler::doinitrun();
	theAdditionalP8Settings.push_back("ProcessLevel:all = off");
	theAdditionalP8Settings.push_back("HadronLevel:Decay = off");
	+ theAdditionalP8Settings.push_back("Check:event = off");
	+ theAdditionalP8Settings.push_back("Next:numberCount = 0");
	+ theAdditionalP8Settings.push_back("Next:numberShowLHA = 0");
	+ theAdditionalP8Settings.push_back("Next:numberShowInfo = 0");
	+ theAdditionalP8Settings.push_back("Next:numberShowProcess = 0");
	+ theAdditionalP8Settings.push_back("Next:numberShowEvent = 0");
	+ theAdditionalP8Settings.push_back("Init:showChangedSettings = 0");
	+ theAdditionalP8Settings.push_back("Init:showAllSettings = 0");
	+ theAdditionalP8Settings.push_back("Init:showChangedParticleData = 0");
	+ theAdditionalP8Settings.push_back("Init:showChangedResonanceData = 0");
	+ theAdditionalP8Settings.push_back("Init:showAllParticleData = 0");
	+ theAdditionalP8Settings.push_back("Init:showOneParticleData = 0");
	+ theAdditionalP8Settings.push_back("Init:showProcesses = 0");
	+
	pythia.init(*this, theAdditionalP8Settings);
	theBE.init(&(pythia().info), pythia().settings, pythia().particleData);
	}

	void BoseEinsteinHandler::persistentOutput(PersistentOStream & os) const {
	os
	#include "BoseEinsteinHandler-output.h"
	<< doContinueDecays;
	}

	void BoseEinsteinHandler::persistentInput(PersistentIStream & is, int) {
	is
	#include "BoseEinsteinHandler-input.h"
	>> doContinueDecays;
	}

	ClassDescription<BoseEinsteinHandler> BoseEinsteinHandler::initBoseEinsteinHandler;
	// Definition of the static class description member.

	void BoseEinsteinHandler::Init() {

	#include "BoseEinsteinHandler-interfaces.h"



	static Switch<BoseEinsteinHandler,bool> interfaceContinueDecays
	("ContinueDecays",
	"Continue the decay chains after the Bose-Einstein shifting.",
	&BoseEinsteinHandler::doContinueDecays, true, true, false);
	static SwitchOption interfaceContinueDecaysYes
	(interfaceContinueDecays,
	"Yes",
	"Continue decays.",
	true);
	static SwitchOption interfaceContinueDecaysNo
	(interfaceContinueDecays,
	"No",
	"Do not continue decays.",
	false);


	}

	diff --git a/Config/Pythia8Interface.cc b/Config/Pythia8Interface.cc
	--- a/Config/Pythia8Interface.cc
	+++ b/Config/Pythia8Interface.cc
	@@ -1,206 +1,207 @@
	// -- C++ --
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the Pythia8Interface class.
	//

	#include "Pythia8Interface.h"
	#include "ThePEG/Repository/CurrentGenerator.h"
	#include "ThePEG/Handlers/StepHandler.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Interface/Parameter.h"

	using namespace TheP8I;

	Pythia8Interface::Pythia8Interface(): pythia(0) {}

	Pythia8Interface::~Pythia8Interface() {
	if ( pythia ) delete pythia;
	}

	bool Pythia8Interface::initialized = false;

	string Pythia8Interface::installDir = PYTHIA8_DIR;

	RndmEngine Pythia8Interface::rnd;

	void Pythia8Interface::Init() {
	if ( initialized ) return;

	if ( installDir[installDir.length() - 1] != '/' ) installDir += '/';

	initialized = true;

	}

	bool Pythia8Interface::go() {
	try {
	return pythia->next();
	} catch ( ... ) {
	Throw<Pythia8ExecException>()
	<< "Pythia8 threw an exception while executing 'Pythia::next()'."
	<< Exception::runerror;
	}
	return false;
	}

	void Pythia8Interface::
	setParameters(const Interfaced & handler, const vector<string> & additional) {
	if ( !pythia ) return;
	InterfaceMap ifs = BaseRepository::getInterfaces(typeid(handler));
	for ( InterfaceMap::iterator it = ifs.begin(); it != ifs.end(); ++it ) {
	string name = it->first;
	string::size_type i = name.find('_');
	ostringstream cmd;
	if ( i == string::npos ) continue;
	while ( ( i = name.find('_') ) != string::npos ) name[i] = ':';
	if ( const SwitchBase * si = dynamic_cast<const SwitchBase *>(it->second) ) {
	if ( si->get(handler) == si->def(handler) ) continue;
	cmd << name << " = " << si->get(handler);
	}
	else if ( const ParameterTBase<double> * pi =
	dynamic_cast<const ParameterTBase<double> *>(it->second) ) {
	if ( pi->tget(handler) == pi->tdef(handler) ) continue;
	cmd << name << " = " << pi->tget(handler);
	}
	else if ( const ParameterTBase<int> * pi =
	dynamic_cast<const ParameterTBase<int> *>(it->second) ) {
	if ( pi->tget(handler) == pi->tdef(handler) ) continue;
	cmd << name << " = " << pi->tget(handler);
	}
	else
	continue;
	pythia->readString(cmd.str());
	}
	for ( int i = 0, N = additional.size(); i < N; ++i )
	pythia->readString(additional[i]);
	}

	void Pythia8Interface::debug() {
	event().list();
	}

	void Pythia8Interface::
	init(const Interfaced & handler, const vector<string> & additional) {
	if ( pythia ) delete pythia;

	Init();

	CurrentGenerator::Redirect stdout(cout);

	pythia = new Pythia8::Pythia(installDir + "xmldoc");

	pythia->setRndmEnginePtr(&rnd);

	CurrentGenerator eg;
	for ( ParticleMap::const_iterator it = eg.current().particles().begin();
	it != eg.current().particles().end(); ++it ) {
	ParticleData & pd = *it->second;
	ostringstream oss;
	if ( pd.id() == 0 ) Throw<Pythia8InitException>()
	<< "The particle " << pd.name()
	<< " was found to have an PDG id of zero." << Exception::warning;
	oss << pd.id();
	if ( pythia->particleData.isParticle(pd.id()) )
	oss << " all ";
	else
	oss << " new ";
	oss << pd.PDGName() << " ";
	if ( pd.CC() ) {
	if ( pd.id() < 0 ) continue;
	ParticleData & apd = *pd.CC();
	oss << apd.PDGName() << " ";
	if ( pd.iSpin() != apd.iSpin() \|\|
	pd.iColour() != -apd.iColour() \|\|
	pd.iCharge() != -apd.iCharge() \|\|
	pd.mass() != apd.mass() \|\|
	pd.width() != apd.width() \|\|
	pd.massMax() != apd.massMax() \|\|
	pd.massMin() != apd.massMin() \|\|
	pd.cTau() != apd.cTau() )
	Throw<Pythia8InitException>()
	<< "The particle " << pd.PDGName()
	<< " did not have the the same properties as its anti-particle. "
	<< "Pythia8 will be initialized with the anti-particle having "
	<< "the same properties as the particle." << Exception::warning;
	}
	else {
	oss << "void ";
	}
	oss << pd.iSpin() << " "
	<< pd.iCharge() << " "
	<< ( pd.iColour() == PDT::Colour3? 1:
	( pd.iColour() == PDT::Colour3bar? -1:
	( pd.iColour() == PDT::Colour8? 2: 0) ) ) << " "
	<< pd.mass()/GeV << " "
	<< pd.width()/GeV << " "
	<< pd.massMin()/GeV << " "
	<< pd.massMax()/GeV << " "
	<< pd.cTau()/millimeter;
	if ( !pythia->particleData.readString(oss.str(), false) )
	Throw<Pythia8InitException>()
	<< "Something went wrong when initializing the particle data table "
	<< "of Pythia8. The settings line was:\n" << oss.str()
	<< "\n" << Exception::runerror;
	+ pythia->particleData.hasChanged(pd.id(), false);
	}

	setParameters(handler, additional);

	pythia->init();

	}


	int Pythia8Interface::addParticle(tPPtr p, int status, int mother1, int mother2) {
	if ( particleIndex.included(p) ) return particleIndex(p);
	int idx = event().size();
	particleIndex(idx, p);
	int pos = event().append(p->id(), status, mother1, mother2, 0, 0,
	addColourLine(p->colourLine()),
	addColourLine(p->antiColourLine()),
	p->momentum().x()/GeV, p->momentum().y()/GeV,
	p->momentum().z()/GeV, p->momentum().e()/GeV,
	p->momentum().mass()/GeV);
	if ( !(p->vertex() == LorentzPoint()) ) {
	event()[pos].vProd(p->vertex().x()/millimeter, p->vertex().y()/millimeter,
	p->vertex().z()/millimeter, p->vertex().t()/millimeter);
	event()[pos].tau(p->lifeLength().tau()/millimeter);
	}

	return idx;
	}

	int Pythia8Interface::addColourLine(tColinePtr c) {
	if ( colourIndex.included(c) ) return colourIndex(c);
	int ctag = event().nextColTag();
	colourIndex(ctag, c);
	return ctag;
	}

	PPtr Pythia8Interface::getParticle(int idx) {
	if ( particleIndex.included(idx) ) return particleIndex.find(idx);
	tcPDPtr pd = CurrentGenerator::current().getParticleData(event()[idx].id());
	if ( !pd ) return PPtr();
	PPtr p = pd->produceParticle(Lorentz5Momentum(event()[idx].px()*GeV,
	event()[idx].py()*GeV,
	event()[idx].pz()*GeV,
	event()[idx].e()*GeV,
	event()[idx].m()*GeV));
	if ( event()[idx].hasVertex() )
	p->setVertex(LorentzPoint(event()[idx].xProd()*millimeter,
	event()[idx].yProd()*millimeter,
	event()[idx].zProd()*millimeter,
	event()[idx].tProd()*millimeter));
	if ( event()[idx].tau() > 0.0 && event()[idx].m() > 0.0 )
	p->setLifeLength(event()[idx].tau()*millimeter
	*p->momentum()/p->momentum().mass());


	tColinePtr c = colourIndex(event()[idx].col());
	if ( c ) c->addColoured(p);
	c = colourIndex(event()[idx].acol());
	if ( c ) c->addAntiColoured(p);
	particleIndex(idx, p);
	return p;

	}
	diff --git a/Hadronization/StringFragmentation.cc b/Hadronization/StringFragmentation.cc
	--- a/Hadronization/StringFragmentation.cc
	+++ b/Hadronization/StringFragmentation.cc
	@@ -1,458 +1,470 @@
	// -- 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/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];
	}

	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");
	+ theAdditionalP8Settings.push_back("Next:numberCount = 0");
	+ theAdditionalP8Settings.push_back("Next:numberShowLHA = 0");
	+ theAdditionalP8Settings.push_back("Next:numberShowInfo = 0");
	+ theAdditionalP8Settings.push_back("Next:numberShowProcess = 0");
	+ theAdditionalP8Settings.push_back("Next:numberShowEvent = 0");
	+ theAdditionalP8Settings.push_back("Init:showChangedSettings = 0");
	+ theAdditionalP8Settings.push_back("Init:showAllSettings = 0");
	+ theAdditionalP8Settings.push_back("Init:showChangedParticleData = 0");
	+ theAdditionalP8Settings.push_back("Init:showChangedResonanceData = 0");
	+ theAdditionalP8Settings.push_back("Init:showAllParticleData = 0");
	+ theAdditionalP8Settings.push_back("Init:showOneParticleData = 0");
	+ theAdditionalP8Settings.push_back("Init:showProcesses = 0");
	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);

	}

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 2:38 PM (1 d, 8 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3804774
Default Alt Text: (31 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions