No OneTemporary
Actions

Size

72 KB

Subscribers

None

View Options

	diff --git a/Config/RopeUserHooks.h b/Config/RopeUserHooks.h
	--- a/Config/RopeUserHooks.h
	+++ b/Config/RopeUserHooks.h
	@@ -1,264 +1,346 @@
	#ifndef THEP8I_RopeUserHooks_H
	#define THEP8I_RopeUserHooks_H

	#include "Pythia.h"
	#include "ThePEG/EventRecord/Particle.h"
	#include "TheP8I/Hadronization/ParameterHandler.h"
	#include "TheP8I/Hadronization/Ropewalk.h"
	#include "ThePEG/Utilities/Throw.h"



	namespace TheP8I {

	class RopeUserHooks : public Pythia8::UserHooks {

	// Convenient typedefs
	typedef map<string,double> PytPars;
	typedef map<Energy2, Ropewalk::Dipole *> DipMass;
	typedef map<tcPPtr,DipMass> FlavourEnd;
	public:

	RopeUserHooks() : _ph(NULL), _h(-1.0), _m0(0GeV), _r0(1femtometer) {
	}

	~RopeUserHooks() {
	}

	virtual bool canChangeFragPar() {
	return true;
	}

	virtual bool doChangeFragPar(Pythia8::StringFlav* flavPtr, Pythia8::StringZ* zPtr, Pythia8::StringPT * pTPtr, int endFlavour, double m2Had, vector<int> iParton) {
	// Not in use here!
	return false;
	}

	virtual bool doChangeFragPar(Pythia8::StringFlav* flavPtr, Pythia8::StringZ* zPtr, Pythia8::StringPT * pTPtr, int endFlavour, double m2Had) {

	// Get new parameters
	PytPars newPar = fetchParameters(endFlavour,m2Had);
	if(newPar.find("null") != newPar.end()) {
	Throw<RopeException>()
	<< "Problem fetching parameters in RopeUserHook. "
	<< "Ropes switched off in this string."
	<< Exception::warning;
	_h = 1.0;
	newPar = fetchParameters(endFlavour,m2Had);
	_h = -1.0;
	}
	for(PytPars::iterator itr = newPar.begin(); itr!=newPar.end();++itr) settingsPtr->parm(itr->first,itr->second);

	// Re-initialize all three
	flavPtr->init(*settingsPtr,rndmPtr);
	zPtr->init(settingsPtr,particleDataPtr,rndmPtr);
	pTPtr->init(settingsPtr,particleDataPtr,rndmPtr);

	return true;
	}

	void setParameterHandler(ParameterHandler * ph){
	_ph = ph;
	}


	void setEnhancement(double h){
	// Will overrule others if set.
	_h = h;
	}

	+ bool setDipoles(pair<ColourSinglet , vector<Ropewalk::Dipole > > * dm, Energy m0, Length r0, bool throwHadr = false, double alpha = 1.0){
	+
	+ // If we set the dipoles, we should also use them.
	+ _h = -1.0;
	+ _m0 = m0;
	+ _r0 = r0;
	+ _throwHadr = throwHadr;
	+ _alpha = alpha;
	+
	+ vector<Ropewalk::Dipole*> dip = dm->second;
	+ // Closed gluon loop -- we will return false and let StringFragmentation figure out what to do
	+ if(dip.front()->pa->id() == 21 && dip.front()->pc->id() == 21 &&
	+ dip.back()->pa->id() == 21 && dip.back()->pc->id() == 21)
	+ return false;
	+
	+
	+ // Get flavour of quark ends -- get the one that is not a gluon.
	+ // We always assign left to the front of the
	+ // dipole vector and right to the back.
	+ DipMass leftC;
	+ DipMass rightC;
	+ leftC.clear();
	+ rightC.clear();
	+ if(dip.size()==1){
	+ _leftP = dip[0]->pc;
	+ _rightP = dip[0]->pa;
	+ }
	+ else{
	+ _leftP = (dip.front()->pa->id() != 21 ? dip.front()->pa : dip.front()->pc);
	+ _rightP = (dip.back()->pc->id() != 21 ? dip.back()->pc : dip.back()->pa);
	+ }
	+ if (_leftP->id() == _rightP->id() )
	+ Throw<RopeException>()
	+ << "Flavours in strings corrupted. "
	+ << "This is a serious error - please contact the authors."
	+ << Exception::abortnow;
	+
	+ // A dipole gets all of quark mass, half of gluon
	+ // The end quarks energy m0
	+
	+ // Start by adding half of the first quark
	+ // from both left and right
	+
	+ LorentzMomentum momL(0GeV,0GeV,0GeV,0GeV);
	+ LorentzMomentum momR(0GeV,0GeV,0GeV,0GeV);
	+
	+
	+ for(size_t i = 0, N = dip.size(); i < N; ++i ){
	+ Ropewalk::Dipole * dPtrL = dip[i];
	+ Ropewalk::Dipole * dPtrR = dip[N - i - 1];
	+
	+ if(!dPtrL \|\| dPtrL->broken \|\| !dPtrR \|\| dPtrR->broken){
	+ Throw<RopeException>()
	+ << "Broken dipole in RopeUserHooks. "
	+ << "This is a serious error - please contact the authors."
	+ << Exception::abortnow;
	+ return false;
	+ }
	+ momL += 0.5*(dPtrL->pc->momentum());
	+ momL += 0.5*(dPtrL->pa->momentum());
	+ momR += 0.5*(dPtrR->pc->momentum());
	+ momR += 0.5*(dPtrR->pa->momentum());
	+
	+ if(i == 0){
	+ momL += 0.5 * _leftP->momentum();
	+ momR += 0.5 * _rightP->momentum();
	+ }
	+ if(i == N - 1){
	+ momR += 0.5 * _leftP->momentum();
	+ momL += 0.5 * _rightP->momentum();
	+ }
	+ leftC.insert(DipMass::value_type(momL.m2(),dPtrL));
	+ rightC.insert(DipMass::value_type(momR.m2(),dPtrR));
	+ }
	+
	+ flavourDipoles.clear();
	+ flavourDipoles.insert( FlavourEnd::value_type(_leftP,leftC) );
	+ flavourDipoles.insert( FlavourEnd::value_type(_rightP,rightC) );
	+
	+ return true;
	+ }
	+
	bool setDipoles(pair<ColourSinglet * const, vector<Ropewalk::Dipole > > dm, Energy m0, Length r0, bool throwHadr = false, double alpha = 1.0){

	// If we set the dipoles, we should also use them.
	_h = -1.0;
	_m0 = m0;
	_r0 = r0;
	_throwHadr = throwHadr;
	_alpha = alpha;

	vector<Ropewalk::Dipole*> dip = dm->second;
	// Closed gluon loop -- we will return false and let StringFragmentation figure out what to do
	if(dip.front()->pa->id() == 21 && dip.front()->pc->id() == 21 &&
	dip.back()->pa->id() == 21 && dip.back()->pc->id() == 21)
	return false;


	// Get flavour of quark ends -- get the one that is not a gluon.
	// We always assign left to the front of the
	// dipole vector and right to the back.
	DipMass leftC;
	DipMass rightC;
	leftC.clear();
	rightC.clear();
	if(dip.size()==1){
	_leftP = dip[0]->pc;
	_rightP = dip[0]->pa;
	}
	else{
	_leftP = (dip.front()->pa->id() != 21 ? dip.front()->pa : dip.front()->pc);
	_rightP = (dip.back()->pc->id() != 21 ? dip.back()->pc : dip.back()->pa);
	}
	if (_leftP->id() == _rightP->id() )
	Throw<RopeException>()
	<< "Flavours in strings corrupted. "
	<< "This is a serious error - please contact the authors."
	<< Exception::abortnow;

	// A dipole gets all of quark mass, half of gluon
	// The end quarks energy m0

	// Start by adding half of the first quark
	// from both left and right

	LorentzMomentum momL(0GeV,0GeV,0GeV,0GeV);
	LorentzMomentum momR(0GeV,0GeV,0GeV,0GeV);


	for(size_t i = 0, N = dip.size(); i < N; ++i ){
	Ropewalk::Dipole * dPtrL = dip[i];
	Ropewalk::Dipole * dPtrR = dip[N - i - 1];

	if(!dPtrL \|\| dPtrL->broken \|\| !dPtrR \|\| dPtrR->broken){
	Throw<RopeException>()
	<< "Broken dipole in RopeUserHooks. "
	<< "This is a serious error - please contact the authors."
	<< Exception::abortnow;
	return false;
	}
	momL += 0.5*(dPtrL->pc->momentum());
	momL += 0.5*(dPtrL->pa->momentum());
	momR += 0.5*(dPtrR->pc->momentum());
	momR += 0.5*(dPtrR->pa->momentum());

	if(i == 0){
	momL += 0.5 * _leftP->momentum();
	momR += 0.5 * _rightP->momentum();
	}
	if(i == N - 1){
	momR += 0.5 * _leftP->momentum();
	momL += 0.5 * _rightP->momentum();
	}
	leftC.insert(DipMass::value_type(momL.m2(),dPtrL));
	rightC.insert(DipMass::value_type(momR.m2(),dPtrR));
	}

	flavourDipoles.clear();
	flavourDipoles.insert( FlavourEnd::value_type(_leftP,leftC) );
	flavourDipoles.insert( FlavourEnd::value_type(_rightP,rightC) );

	return true;
	}

	public:

	/**
	* Exception class.
	*/
	struct RopeException: public Exception {};

	private:

	PytPars fetchParameters(int endFlavour, double m2Had){

	// Did we remember to set the callback pointer?
	// Do not just construct a new one, as there will be no sensible initial values, resulting
	// in very subtle errors.
	if(!(_ph)){
	Throw<RopeException>()
	<< "Missing parameter handler in RopeUserHooks. "
	<< "This is a serious error - please contact the authors."
	<< Exception::abortnow;
	PytPars p;
	p.insert(make_pair<const string,double>("null",0.));
	return p;
	}

	// If we have manually set enhancement, we should just use that, and avoid complicated behavior.
	if(_h > 0) return _ph->GetEffectiveParameters(_h);

	// Test the string ends...
	if( (endFlavour == _leftP->id() && endFlavour == _rightP->id()) \|\|
	(endFlavour != _leftP->id() && endFlavour != _rightP->id()) ){
	Throw<RopeException>()
	<< "Flavours in strings corrupted. "
	<< "This is a serious error - please contact the authors."
	<< Exception::abortnow;
	PytPars p;
	p.insert(make_pair<const string,double>("null",0.));
	return p;
	}

	// Find the right end
	tcPPtr thisEnd = (_leftP->id() == endFlavour ? _leftP : _rightP);
	FlavourEnd::iterator feItr = flavourDipoles.find( thisEnd);

	if( feItr == flavourDipoles.end()){
	Throw<RopeException>()
	<< "Could not find dipole end. "
	<< "Ropes switched off in this string."
	<< Exception::warning;
	PytPars p;
	p.insert(make_pair<const string,double>("null",0.));
	return p;
	}
	DipMass dms = feItr->second;

	// Find the first dipole where the total invariant mass sq. from the string
	// (left or right) exceeds m2Had
	DipMass::iterator dmsItr = dms.lower_bound(m2Had*GeV2);
	if( dmsItr == dms.end()){
	Throw<RopeException>()
	<< "Could not get correct dipole. mpythia: " << m2Had
	<< " " << (--dmsItr)->first/GeV2 << ". "
	<< "Ropes switched off in this string."
	<< Exception::warning;
	PytPars p;
	p.insert(make_pair<const string,double>("null",0.));
	return p;
	}

	Ropewalk::Dipole * dPtr = dmsItr->second;
	if(!dPtr){
	Throw<RopeException>()
	<< "Missing dipole in RopeUserHooks. "
	<< "This is a serious error - please contact the authors."
	<< Exception::abortnow;
	PytPars p;
	p.insert(make_pair<const string,double>("null",0.));
	return p;
	}
	// Find out how long in we are on the dipole
	double dipFrac;
	double mBig = dmsItr->first/GeV2;

	if(m2Had == 0) dipFrac = 0;
	else if( dmsItr == dms.begin() ) dipFrac = sqrt(m2Had/mBig);
	else{
	double mSmall = double((--dmsItr)->first/GeV2);
	dipFrac = (sqrt(m2Had) - sqrt(mSmall)) / (sqrt(mBig) - sqrt(mSmall));
	}
	// We need dipFrag to be fraction from the pc parton in the dipole.
	// Check if the starting end is pa or pc. If it is already a pc, do noting,
	// otherwise take fraction from other side
	if( thisEnd == dms.lower_bound(0*GeV2)->second->pa) dipFrac = 1 - dipFrac;
	dPtr->reinit(dipFrac,_r0,_m0);
	if (_throwHadr){
	(*dPtr).hadr = true;
	return _ph->GetEffectiveParameters(dPtr->firstKappa(_alpha));
	}
	// cout << "p = " << dPtr->p << " q = " << dPtr->q << " kappa " << << endl;
	return _ph->GetEffectiveParameters(dPtr->kappaEnhancement());

	}

	tcPPtr _leftP, _rightP;
	FlavourEnd flavourDipoles;
	ParameterHandler * _ph;
	double _h, _alpha;
	Energy _m0;
	Length _r0;
	bool _throwHadr;
	};

	}
	#endif
	diff --git a/Hadronization/Ropewalk.cc b/Hadronization/Ropewalk.cc
	--- a/Hadronization/Ropewalk.cc
	+++ b/Hadronization/Ropewalk.cc
	@@ -1,560 +1,565 @@
	// -- C++ --
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the Ropewalk class.
	//

	#include "Ropewalk.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/EventRecord/Particle.h"
	#include "ThePEG/EventRecord/Step.h"
	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Repository/CurrentGenerator.h"
	#include "ThePEG/PDT/EnumParticles.h"
	#include "ThePEG/Utilities/Selector.h"
	#include "ThePEG/Utilities/SimplePhaseSpace.h"
	#include "ThePEG/Utilities/UtilityBase.h"
	#include "ThePEG/Utilities/Throw.h"
	#include "ThePEG/Handlers/StepHandler.h"



	using namespace TheP8I;

	Ropewalk::Ropewalk(const vector<ColourSinglet> & singlets, Length width,
	Energy scale, double jDP, bool throwaway, bool deb)
	: R0(width), m0(scale), junctionDiquarkProb(jDP), debug(deb),
	strl0(0.0), strl(0.0), avh(0.0), avp(0.0), avq(0.0) {
	for ( int is = 0, Ns = singlets.size(); is < Ns; ++is ){
	stringdipoles[cloneToFinal(singlets[is])];
	}
	getDipoles();

	lambdaBefore = lambdaSum();
	if ( debug ) CurrentGenerator::log() << singlets.size() << "s "
	<< dipoles.size() << "d ";
	calculateOverlaps();
	if(throwaway) doBreakups();

	}

	Ropewalk::~Ropewalk() {
	for ( DipoleMap::iterator it = stringdipoles.begin();
	it != stringdipoles.end(); ++it ) delete it->first;
	}

	ColourSinglet * Ropewalk::cloneToFinal(const ColourSinglet & cs) {
	tcParticleSet final;
	for ( int i = 0, N = cs.partons().size(); i < N; ++i )
	final.insert(cs.partons()[i]->final());
	tcPPtr p = *final.begin();
	if ( p->colourLine() ) return new ColourSinglet(p->colourLine(), final);
	if ( p->antiColourLine() ) return new ColourSinglet(p->antiColourLine(), final);
	Throw<ColourException>()
	<< "Cloning ColourSinglets failed in Ropewalk. "
	<< "This is a serious error - please contact the authors."
	<< Exception::abortnow;
	return 0;
	}

	LorentzPoint Ropewalk::
	propagate(const LorentzPoint & b, const LorentzMomentum & p) const {
	return b + p/(p.e()*m0/hbarc);
	// return b;
	}

	void Ropewalk::getDipoles() {

	// First extract all dipoles from all strings (pieces).
	for ( DipoleMap::iterator it = stringdipoles.begin();
	it != stringdipoles.end(); ++it ) {
	ColourSinglet & string = *it->first;
	for ( int isp = 1, Nsp = string.nPieces(); isp <= Nsp; ++isp ) {
	const ColourSinglet::StringPiece & sp = string.piece(isp);
	if ( sp.size() < 2 ) continue;
	bool forward = sp[0]->colourLine() &&
	sp[0]->colourLine() == sp[1]->antiColourLine();
	for ( int i = 0, N = sp.size(); i < (N - 1); ++i ) {
	if ( forward )
	dipoles.push_back(Dipole(sp[i], sp[i + 1], string));
	else
	dipoles.push_back(Dipole(sp[i + 1], sp[i], string));
	}
	if ( !forward && sp.front()->colourLine() &&
	sp.front()->colourLine() == sp.back()->antiColourLine() )
	dipoles.push_back(Dipole(sp.front(), sp.back(), string));
	else if ( forward && sp.front()->antiColourLine() &&
	sp.front()->antiColourLine() == sp.back()->colourLine() )
	dipoles.push_back(Dipole(sp.back(), sp.front(), string));
	}
	}

	for ( int i = 0, N = dipoles.size(); i < N; ++i ) {
	stringdipoles[dipoles[i].string].push_back(&dipoles[i]);
	strl0 += dipoles[i].yspan(1.0*GeV);
	}

	}

	double Ropewalk::limitrap(const LorentzMomentum & p) const {
	if ( p.z() == ZERO ) return 0.0;
	Energy mt = sqrt(max(p.perp2() + p.m2(), sqr(m0)));
	double rap = log((p.t() + abs(p.z()))/mt);
	return p.z() > ZERO? rap: -rap;
	}

	Ropewalk::OverlappingDipole::
	OverlappingDipole(const Dipole & d, const LorentzRotation R, const Ropewalk * rw)
	: dipole(&d), dir(1) {
	// Get the boost to other dipole's rest frame and calculate
	// rapidities.
	bc = R*rw->propagate(d.pc->vertex(), d.pc->momentum());
	ba = R*rw->propagate(d.pa->vertex(), d.pa->momentum());
	yc = rw->limitrap(R*d.pc->momentum());
	ya = rw->limitrap(R*d.pa->momentum());
	if ( yc < ya ) dir = -1;
	}

	void Ropewalk::calculateOverlaps() {

	// Then calculate all overlaps between pairs of dipoles.
	for ( int i1 = 0, Nd = dipoles.size(); i1 < Nd; ++i1 ) {
	Dipole & d1 = dipoles[i1];
	if ( d1.s() <= sqr(m0) ) continue;
	// Get the boost to dipoles rest frame and calculate rapidities.
	LorentzMomentum pc1 = d1.pc->momentum();
	LorentzMomentum pa1 = d1.pa->momentum();
	LorentzRotation R = Utilities::boostToCM(make_pair(&pc1, &pa1));
	d1.bc = R*propagate(d1.pc->vertex(), d1.pc->momentum());
	d1.ba = R*propagate(d1.pa->vertex(), d1.pa->momentum());
	double yc1 = limitrap(pc1);
	double ya1 = limitrap(pa1);
	double n = 0.0;
	double m = 0.0;
	if ( yc1 <= ya1 ) continue; // don't care about too small strings.
	for ( int i2 = 0; i2 < Nd; ++i2 ) {
	if ( i1 == i2 ) continue;
	// Boost to the rest frame of dipole 1.
	if ( dipoles[i2].s() <= sqr(m0) ) continue;
	OverlappingDipole od(dipoles[i2], R, this);
	// Ignore if not overlapping in rapidity.
	if ( min(od.yc, od.ya) > yc1 \|\| max(od.yc, od.ya) < ya1 \|\| od.yc == od.ya )
	continue;
	// Calculate rapidity overlap.
	double yomax = min(max(od.yc, od.ya), yc1);
	double yomin = max(min(od.yc, od.ya), ya1);
	// Sample a random point in the rapidity overlap region and get
	// positions in space and check overlap.
	double yfrac = (yomax - yomin)/(yc1 - ya1);
	double y = UseRandom::rnd(yomin, yomax);
	if ( !od.overlap(y, d1.ba + (d1.bc - d1.ba)*(y - ya1)/(yc1 - ya1), R0) )
	yfrac = 0.0;
	// Sum overlaps.
	if ( od.dir > 0 ) {
	n += yfrac;
	} else {
	m += yfrac;
	}
	od.yfrac = yfrac*od.dir;
	d1.overlaps.push_back(od);
	}
	d1.n = int(n + UseRandom::rnd());
	d1.m = int(m + UseRandom::rnd());
	d1.initWeightMultiplet();
	avp += d1.pd1.yspan(1.0GeV);
	avq += d1.qd1.yspan(1.0GeV);
	}
	}

	double Ropewalk::Dipole::reinit(double ry, Length R0, Energy m0) {
	// First calculate the rapidity in this restframe.
	double y = yspan(m0)*(ry - 0.5);
	// Then get the position in space of the string at this rapidity.
	LorentzPoint b = ba + (bc - ba)*ry;
	p = 1;
	q = 0;
	// Now go through all overlapping dipoles.
	for ( int i = 0, N = overlaps.size(); i < N; ++i ) {
	if ( overlaps[i].dipole->broken ) continue;
	if( overlaps[i].dipole->hadr ) continue;
	if ( overlaps[i].overlap(y, b, R0) ) {
	if ( overlaps[i].dir > 0 ) ++p;
	else ++q;
	}
	}
	return kappaEnhancement();
	}

	void Ropewalk::Dipole::initMultiplet() {
	typedef pair<int,int> Plet;

	int ns = 0;
	int ms = 0;
	while ( ns < n \|\| ms < m ) {
	Plet diff;
	double octetveto = double(ns + ms + 1 - p - q)/double(ns + ms + 1);
	if ( double(n - ns) > UseRandom::rnd()*double(n + m - ns - ms) ) {
	Selector<Plet> sel;
	sel.insert(multiplicity(p + 1, q), Plet(1, 0));
	sel.insert(multiplicity(p, q - 1)*octetveto, Plet(0, -1));
	sel.insert(multiplicity(p - 1, q + 1), Plet(-1, 1));
	diff = sel[UseRandom::rnd()];
	++ns;
	} else {
	Selector<Plet> sel;
	sel.insert(multiplicity(p, q + 1), Plet(0, 1));
	sel.insert(multiplicity(p - 1, q)*octetveto, Plet(-1, 0));
	sel.insert(multiplicity(p + 1, q - 1), Plet(1, -1));
	diff = sel[UseRandom::rnd()];
	++ms;
	}
	if ( diff.first == -diff.second ) nb++;
	p += diff.first;
	q += diff.second;
	}
	/* * ATTENTION * maybe better if Christian implements this */

	}

	void Ropewalk::Dipole::initNewMultiplet() {
	typedef pair<int,int> Plet;

	int ns = 0;
	int ms = 0;
	for ( int i = 0, N = overlaps.size(); i < N; ++i ) {
	if ( abs(overlaps[i].yfrac) < UseRandom::rnd() ) continue;
	Plet diff(0,0);
	double octetveto = double(ns + ms + 1 - p - q)/double(ns + ms + 1);
	if ( overlaps[i].dir > 0 ) {
	Selector<Plet> sel;
	sel.insert(multiplicity(p + 1, q), Plet(1, 0));
	sel.insert(multiplicity(p, q - 1)*octetveto, Plet(0, -1));
	sel.insert(multiplicity(p - 1, q + 1), Plet(-1, 1));
	diff = sel[UseRandom::rnd()];
	++ns;
	} else {
	Selector<Plet> sel;
	sel.insert(multiplicity(p, q + 1), Plet(0, 1));
	sel.insert(multiplicity(p - 1, q)*octetveto, Plet(-1, 0));
	sel.insert(multiplicity(p + 1, q - 1), Plet(1, -1));
	diff = sel[UseRandom::rnd()];
	++ms;
	}
	if ( diff.first == -diff.second ) nb++;
	p += diff.first;
	q += diff.second;
	}
	n = ns;
	m = ms;
	/* * ATTENTION * maybe better if Christian implements this */

	}

	void Ropewalk::Dipole::initWeightMultiplet() {
	typedef pair<int,int> Plet;

	int ns = 0;
	int ms = 0;
	double mab = sqrt(s())/GeV;
	for ( int i = 0, N = overlaps.size(); i < N; ++i ) {
	if ( abs(overlaps[i].yfrac) < UseRandom::rnd() ) continue;
	const Dipole * od = overlaps[i].dipole;
	double mcd = sqrt(od->s())/GeV;
	double w8 = 1.0/sqr(mab*mcd); //The weight for 3 + 3bar -> 8
	double w6 = w8; // The weight for 3 + 3 -> 6
	Plet diff(0,0);
	double octetveto = double(ns + ms + 1 - p - q)/double(ns + ms + 1);
	if ( overlaps[i].dir > 0 ) {
	double mac = (pc->momentum() + od->pc->momentum()).m()/GeV;
	double mbd = (pa->momentum() + od->pa->momentum()).m()/GeV;
	double w1 = octetveto/sqr(mac*mbd); // The weight for 3 + 3bar -> 1
	double w3 = 1.0/(mabmcdmac*mbd); //The weight for 3 + 3 -> 3bar
	Selector<Plet> sel;
	sel.insert(multiplicity(p + 1, q) +
	multiplicity(p, q - 1)*w8/(w8 + w1) +
	multiplicity(p - 1, q + 1)*w6/(w6 + w3), Plet(1, 0));
	sel.insert(multiplicity(p, q - 1)*w1/(w1 + w8), Plet(0, -1));
	sel.insert(multiplicity(p - 1, q + 1)*w3/(w3 + w6), Plet(-1, 1));
	diff = sel[UseRandom::rnd()];
	++ns;
	} else {
	double mac = (pa->momentum() + od->pc->momentum()).m()/GeV;
	double mbd = (pc->momentum() + od->pa->momentum()).m()/GeV;
	if ( mac*mbd <= 0.0 ) {
	q++;
	continue;
	}
	double w1 = octetveto/sqr(mac*mbd); // The weight for 3 + 3bar -> 1
	double w3 = 1.0/(mabmcdmac*mbd); //The weight for 3 + 3 -> 3bar
	Selector<Plet> sel;
	sel.insert(multiplicity(p, q + 1) +
	multiplicity(p - 1, q)*w8/(w8 + w1) +
	multiplicity(p + 1, q - 1)*w6/(w6 + w3), Plet(0, 1));
	sel.insert(multiplicity(p - 1, q)*w1/(w1 + w8), Plet(-1, 0));
	sel.insert(multiplicity(p + 1, q - 1)*w3/(w3 + w6), Plet(1, -1));
	diff = sel[UseRandom::rnd()];
	++ms;
	}
	if ( diff.first == -diff.second ) nb++;
	p += diff.first;
	q += diff.second;
	}
	n = ns;
	m = ms;
	/* * ATTENTION * maybe better if Christian implements this */

	}

	double Ropewalk::getNb(){
	double nba = 0;
	- for ( int i = 0, N = dipoles.size(); i < N; ++i )
	- // nba += double(dipoles[i].nb);
	- nba += double(dipoles[i].n + dipoles[i].m + 1 - dipoles[i].p - dipoles[i].q);
	+ for ( int i = 0, N = dipoles.size(); i < N; ++i ){
	+ nba += double(dipoles[i].nb);
	+ //cout << dipoles[i].n << " " << dipoles[i].m << " " << dipoles[i].p << " " << dipoles[i].q << endl;
	+ //cout << double(dipoles[i].n + dipoles[i].m + 1 - dipoles[i].p - dipoles[i].q) << endl;
	+ //cout << " --- " << endl;
	+ //nba += double(dipoles[i].n + dipoles[i].m + 1 - dipoles[i].p - dipoles[i].q);
	+ }
	+ // cout << nba << endl;
	return nba;
	}

	double Ropewalk::getkW(){
	double kw = 0;
	for ( int i = 0, N = dipoles.size(); i < N; ++i )
	kw += dipoles[i].kappaEnhancement()*log(dipoles[i].s()/sqr(m0));
	return kw;
	}

	pair<double,double> Ropewalk::getMN(){
	pair<double,double> ret = make_pair<double,double>(0,0);
	for (int i = 0, N = dipoles.size(); i < N; ++i){
	ret.first += dipoles[i].m;
	ret.second += dipoles[i].n;
	}
	return ret;
	}

	double Ropewalk::lambdaSum(){
	double ret = 0;
	for (int i = 0, N = dipoles.size(); i < N; ++i)
	ret += dipoles[i].s() > 0.1*GeV2 ? log(dipoles[i].s()/sqr(m0)) : 0;
	return ret;
	}


	double Ropewalk::Dipole::breakupProbability() const {
	if ( !breakable() \|\| n + m <= 0.0 \|\| n + m + 1 == p + q ) return 0.0;
	double sum = 0.0;
	for (int j = 0, N = overlaps.size(); j < N; ++j )
	if ( overlaps[j].dipole->breakable() )
	sum += abs(overlaps[j].yfrac);
	if ( sum <= 0.0 ) return 1.0;
	return double(n + m + 1 - p - q)/(sum + 1.0);
	}

	Step & Ropewalk::step() const {
	return *(CurrentGenerator()->currentStepHandler()->newStep());
	}

	void Ropewalk::doBreakups() {

	typedef multimap<double, const Dipole *> BMap;
	BMap breakups;

	// First order alldipoles in increasing probability of breaking
	for ( int i = 0, N = dipoles.size(); i < N; ++i )
	breakups.insert(make_pair(dipoles[i].breakupProbability(), &dipoles[i]));

	// Then start breaking in reverse order
	mspec.clear();
	for ( BMap::reverse_iterator it = breakups.rbegin();
	it != breakups.rend(); ++it ) {
	const Dipole & d = *it->second;
	if ( d.breakupProbability() < UseRandom::rnd() ) continue;
	mspec.push_back(d.s()/GeV2);
	double hinv = 1.0 / d.kappaEnhancement();
	// The default parameters should be set
	// properly from the outside. Now just Pythia8 defaults
	double rho = pow(0.19, hinv);
	double x = pow(1.0,hinv);
	double y = pow(0.027,hinv);

	Selector<int> qsel;
	qsel.clear();
	if(UseRandom::rnd() < junctionDiquarkProb){
	// Take just ordinary quarks
	qsel.insert(1.0, ParticleID::ubar);
	qsel.insert(1.0, ParticleID::dbar);
	qsel.insert(rho, ParticleID::sbar);
	}
	else {
	// Take diquarks
	// TODO (Pythia 8.2) Set status code of these diquarks to 74
	qsel.insert(1.0, ParticleID::ud_0);
	qsel.insert(3.0*y, ParticleID::dd_1);
	qsel.insert(3.0*y, ParticleID::ud_1);
	qsel.insert(3.0*y, ParticleID::uu_1);
	qsel.insert(rho*x, ParticleID::su_0);
	qsel.insert(3xrho*y, ParticleID::su_1);
	qsel.insert(3yxxrho*rho, ParticleID::ss_1);
	qsel.insert(rho*x, ParticleID::sd_0);
	qsel.insert(3xrho*y, ParticleID::sd_1);
	}
	int flav = qsel[UseRandom::rnd()];
	PPtr dic = CurrentGenerator()->getParticle(-flav);
	PPtr dia = CurrentGenerator()->getParticle(flav);
	//cout << (dic->momentum() + dia->momentum()).m2() / GeV2 << " " << d.s() / GeV2 << endl;
	if((dic->momentum() + dia->momentum()).m2() > d.s()) continue;
	// Get boost to dipole rest frame and place the di-quarks there.
	LorentzRotation R = Utilities::getBoostFromCM(make_pair(d.pc, d.pa));

	SimplePhaseSpace::CMS(dic, dia, d.s(), 1.0, 0.0);
	dia->transform(R);
	dic->transform(R);
	// Add the di-quarks as children to the decayed gluons and split
	// the resulting string.
	if ( !( step().addDecayProduct(d.pc, dic, true) &&
	step().addDecayProduct(d.pa, dia, true) &&
	step().addDecayProduct(d.pc, dia, false) &&
	step().addDecayProduct(d.pa, dic, false) ) )
	Throw<ColourException>()
	<< "Adding decay products failed in Ropewalk. "
	<< "This is a serious error - please contact the authors."
	<< Exception::abortnow;
	tcParticleSet pset(d.string->partons().begin(), d.string->partons().end());
	pset.erase(d.pc);
	pset.erase(d.pa);
	pset.insert(dic);
	pset.insert(dia);
	// Reset the neighboring dipoles to point to the new diquarks.
	d.broken = true;
	ColourSinglet * olds = d.string;
	const vector<Dipole *> & oldd = stringdipoles[olds];
	for ( int id = 0, Nd = oldd.size(); id < Nd; ++id ) {
	if ( oldd[id]->pc == d.pa ) oldd[id]->pc = dia;
	if ( oldd[id]->pa == d.pc ) oldd[id]->pa = dic;
	}
	// remove the old string and insert the new ones fixing the
	// pointers from the dipolmes.
	vector<ColourSinglet> newstrings =
	ColourSinglet::getSinglets(pset.begin(), pset.end());
	for ( int is = 0, Ns = newstrings.size(); is < Ns; ++is ) {
	ColourSinglet * news = new ColourSinglet(newstrings[is]);
	tcParticleSet pset(news->partons().begin(), news->partons().end());
	vector<Dipole *> & newd = stringdipoles[news];
	for ( int id = 0, Nd = oldd.size(); id < Nd; ++id ) {
	if ( !oldd[id]->broken && pset.find(oldd[id]->pc) != pset.end() ) {
	newd.push_back(oldd[id]);
	oldd[id]->string = news;
	}
	}
	sort(newd, news->partons()[0]);
	}
	stringdipoles.erase(olds);
	delete olds;
	}
	}

	vector< pair<ColourSinglet,double> > Ropewalk::getSinglets(double DeltaY) const {
	vector< pair<ColourSinglet,double> > ret;
	ret.reserve(stringdipoles.size());
	double toty = 0.0;
	double totyh = 0.0;
	double toth = 0.0;
	double toth2 = 0.0;
	double minh = 10.0;
	double maxh = 0.0;
	// Calculate the kappa-enhancement of the remaining strings.
	for ( DipoleMap::const_iterator it = stringdipoles.begin();
	it != stringdipoles.end(); ++it ) {
	double sumyh = 0.0;
	double sumy = 0.0;
	for ( int id = 0, Nd = it->second.size(); id < Nd; ++id ) {
	double y = it->second[id]->yspan(m0);
	double h = it->second[id]->kappaEnhancement();
	avh += hit->second[id]->yspan(1.0GeV);
	strl += it->second[id]->yspan(1.0*GeV);
	if ( DeltaY > 0.0 &&
	abs(it->second[id]->pc->eta()) > DeltaY &&
	abs(it->second[id]->pa->eta()) > DeltaY ) continue;
	sumy += y;
	sumyh += y*h;
	}
	toty += sumy;
	totyh += sumyh;

	if ( sumy > 0.0 ) {
	double h = 0.1int(10.0sumyh/sumy + UseRandom::rnd());
	ret.push_back(make_pair(*it->first, h));
	toth += sumyh/sumy;
	toth2 += sqr(sumyh/sumy);
	minh = min(minh, sumyh/sumy);
	maxh = max(maxh, sumyh/sumy);
	} else {
	ret.push_back(make_pair(*it->first, 1.0));
	toth += 1.0;
	toth2 += 1.0;
	minh = min(minh, 1.0);
	maxh = max(maxh, 1.0);
	}
	}
	if ( debug )
	CurrentGenerator::log() << ret.size() << "ns "
	<< totyh/max(toty, 0.00001) << "hwa "
	<< minh << "<h "
	<< toth/ret.size() << "ha "
	<< maxh << ">h "
	<< sqrt(toth2/ret.size() - sqr(toth/ret.size())) << "hsd " << endl;

	if ( avh > 0.0 ) avh /= strl;
	if ( avp > 0.0 ) avp /= strl0;
	if ( avq > 0.0 ) avq /= strl0;

	return ret;

	}


	void Ropewalk::sort(vector<Dipole *> & dips, tcPPtr p) {
	if ( p->colourLine() ) {
	map<tcPPtr, Dipole *> dmap;
	Dipole * current = 0;
	for ( int i = 0, N = dips.size(); i < N; ++i ) {
	if ( dips[i]->pa == p ) current = dips[i];
	else dmap[dips[i]->pa] = dips[i];
	}
	dips.clear();
	while ( current ) {
	dips.push_back(current);
	map<tcPPtr, Dipole *>::iterator it = dmap.find(current->pc);
	if ( it == dmap.end() ) current = 0;
	else {
	current = it->second;
	dmap.erase(it);
	}
	}
	if ( !dmap.empty() )
	Throw<ColourException>()
	<< "Failed to sort dipoles in Ropewalk. "
	<< "This is a serious error - please contact the authors."
	<< Exception::abortnow;
	} else {
	map<tcPPtr, Dipole *> dmap;
	Dipole * current = 0;
	for ( int i = 0, N = dips.size(); i < N; ++i ) {
	if ( dips[i]->pc == p ) current = dips[i];
	else dmap[dips[i]->pc] = dips[i];
	}
	dips.clear();
	while ( current ) {
	dips.push_back(current);
	map<tcPPtr, Dipole *>::iterator it = dmap.find(current->pa);
	if ( it == dmap.end() ) current = 0;
	else {
	current = it->second;
	dmap.erase(it);
	}
	}
	if ( !dmap.empty() )
	Throw<ColourException>()
	<< "Failed to sort dipoles in Ropewalk. "
	<< "This is a serious error - please contact the authors."
	<< Exception::abortnow;
	}
	}

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

	#include "StringFragmentation.h"
	#include "Ropewalk.h"
	#include "RandomAverageHandler.h"
	#include "RandomHandler.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/Utilities/EnumIO.h"
	#include "ThePEG/Utilities/HoldFlag.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	-
	+#include <algorithm>
	using namespace TheP8I;


	StringFragmentation::StringFragmentation()
	: pythia(), fScheme(0), stringR0(0.5femtometer), stringm0(1.0GeV), junctionDiquark(1.0), alpha(1.0), average(true),
	stringyCut(2.0), stringpTCut(6*GeV), fragmentationMass(0.134),
	baryonSuppression(0.5), window(false), throwaway(false), _eventShapes(0),
	useThrustAxis(0), opHandler(0), maxTries(2), doAnalysis(0), analysisPath(""),
	#include "StringFragmentation-init.h"
	{
	}

	StringFragmentation::~StringFragmentation() {

	if ( opHandler ) delete opHandler;

	}

	void StringFragmentation::handle(EventHandler & eh, const tPVector & tagged,
	const Hint & h) {
	++nev;
	static unsigned long lastEventId = 0;
	bool secondary = false;
	// Get the event
	tPVector all = RemnantDecayer::decayRemnants(tagged, *newStep());
	HoldFlag<int> oldFS(fScheme, fScheme);
	// Prevent funny behavior if hadronizing decays of eg. Upsilon.
	if ( newStep()->collision()->uniqueId == lastEventId ){
	secondary = true;
	if(fScheme == 4 \|\| fScheme == 5 \|\| fScheme == 1) fScheme = 99;
	else fScheme = 0;
	}
	else
	lastEventId = newStep()->collision()->uniqueId;

	if(_eventShapes) _eventShapes->reset(all);
	vector<ColourSinglet> singlets;
	singlets.clear();

	if ( theCollapser && !secondary )
	singlets = theCollapser->collapse(all, newStep());
	else
	singlets = ColourSinglet::getSinglets(all.begin(), all.end());
	// Goto correct hadronization scheme
	switch(fScheme){
	case 0: // Pythia
	{
	hadronizeSystems(*opHandler->GetPythiaPtr(1.0,nev%10000==0), singlets, all);
	break;
	}
	case 1: // For playing around and printing stuff
	{
	- Ropewalk ropewalk(singlets, stringR0, hbarc/stringR0, baryonSuppression, throwaway, false);
	+ if(throwaway){
	+ Ropewalk ropewalk_init(singlets, stringR0, stringm0, baryonSuppression, throwaway, false);
	+ vector< pair<ColourSinglet,double> > allStrings = ropewalk_init.getSinglets(stringyCut);
	+ tPVector allParticles;
	+ for(vector< pair<ColourSinglet,double> >::iterator sItr = allStrings.begin(); sItr != allStrings.end(); ++sItr)
	+ for(tcPVector::iterator pItr = sItr->first.partons().begin(); pItr != sItr->first.partons().end(); ++pItr)
	+ allParticles.push_back(const_ptr_cast<tPPtr>(*pItr));
	+ singlets = theCollapser->collapse(allParticles, newStep());
	+ }
	+ Ropewalk ropewalk(singlets, stringR0, stringm0, baryonSuppression, false, false);
	+

	vector< pair<ColourSinglet,double> > strings = ropewalk.getSinglets(stringyCut);

	static ofstream out(( generator()->runName() + ".txt").c_str());

	if(!isnan(ropewalk.lambdaSum()+ropewalk.getkW()+ropewalk.getNb()))
	out << generator()->currentEvent()->weight() << " " << ropewalk.lambdaSum() << " " << ropewalk.getkW() << " " << ropewalk.getNb() << endl;
	//vector<double> mspec = ropewalk.mspec;
	//for(int i = 0, N = mspec.size(); i < N; ++i) out << generator()->currentEvent()->weight() << " " << mspec[i] << endl;
	break;

	}
	case 2: // Ropewalk/Dipole
	{
	Ropewalk ropewalk(singlets, stringR0, stringm0, junctionDiquark, throwaway, false);
	vector< pair<ColourSinglet,double> > strings = ropewalk.getSinglets(stringyCut);
	vector<ColourSinglet> toHadronization;
	for(int i = 0, N = strings.size(); i < N; ++i ){
	Pythia8Interface * pytp = opHandler->GetPythiaPtr(strings[i].second,nev%10000==0);
	toHadronization.clear();
	toHadronization.push_back(strings[i].first);
	hadronizeSystems(*pytp,toHadronization,all);
	}
	break;
	}
	case 3: // Pipes with average
	{
	RandomAverageHandler avg_enhancer(throwaway);
	// TODO: Implement throwaway functionality in this
	avg_enhancer.clear();
	vector<StringPipe> pipes;
	// Make all the pipes
	for(vector<ColourSinglet>::iterator sItr = singlets.begin(); sItr!=singlets.end(); ++sItr)
	pipes.push_back(StringPipe(&(*sItr),stringR0,_eventShapes));

	avg_enhancer.SetEvent(pipes);
	for(vector<StringPipe>::iterator it = pipes.begin(); it!=pipes.end(); ++it){
	vector<ColourSinglet> toHadronization;
	double h = avg_enhancer.KappaEnhancement(*it);
	toHadronization.push_back((it).GetSingletPtr());
	forcerun = true;
	hadronizeSystems(*(opHandler->GetPythiaPtr(h,nev%10000==0)),toHadronization,all);
	}
	break;
	}
	case 4: // Average kappa over whole strings
	{
	Ropewalk ropewalk(singlets, stringR0, stringm0, baryonSuppression,throwaway, false);
	vector< pair<ColourSinglet,double> > strings = ropewalk.getSinglets(stringyCut);
	vector<ColourSinglet> toHadronization;
	double avge = 0;
	for(int i = 0, N = strings.size(); i < N; ++i ){
	avge += strings[i].second;
	pythia.getRopeUserHooksPtr()->setEnhancement(strings[i].second);
	toHadronization.clear();
	toHadronization.push_back(strings[i].first);
	hadronizeSystems(pythia,toHadronization,all);
	}
	// ofstream out(( "/scratch/galette/bierlich/work/dipsy/pprange/" + generator()->runName() + "/stringplot.txt").c_str(),ios::app);
	// out << generator()->currentEvent()->weight() << " " << strings.size() << " " << avge << endl;
	break;
	}
	case 5: // Altering kappa per dipole
	{
	- Ropewalk ropewalk(singlets, stringR0, stringm0, baryonSuppression, throwaway, false);
	+
	+ if(throwaway){
	+ Ropewalk ropewalk_init(singlets, stringR0, stringm0, baryonSuppression, throwaway, false);
	+ vector< pair<ColourSinglet,double> > allStrings = ropewalk_init.getSinglets(stringyCut);
	+ tPVector allParticles;
	+ for(vector< pair<ColourSinglet,double> >::iterator sItr = allStrings.begin(); sItr != allStrings.end(); ++sItr)
	+ for(tcPVector::iterator pItr = sItr->first.partons().begin(); pItr != sItr->first.partons().end(); ++pItr)
	+ allParticles.push_back(const_ptr_cast<tPPtr>(*pItr));
	+ singlets = theCollapser->collapse(allParticles, newStep());
	+ }
	+ Ropewalk ropewalk(singlets, stringR0, stringm0, baryonSuppression, false, false);
	vector<ColourSinglet> toHadronization;
	vector< pair<ColourSinglet,double> > gloops = ropewalk.getSinglets(stringyCut);
	- int i = 0;
	+ // Let's see if we can get the ratios to drop
	+ // at high pT by introducing an ordering of the strings
	+ // Next 8 lines is the old version
	+ /*int i = 0;
	for(map<ColourSinglet , vector<Ropewalk::Dipole > >::iterator itr = ropewalk.begin(); itr!=ropewalk.end(); ++itr, ++i){
	if(!pythia.getRopeUserHooksPtr()->setDipoles(&(*itr),stringm0, stringR0, !average, alpha))
	pythia.getRopeUserHooksPtr()->setEnhancement(gloops[i].second);
	toHadronization.clear();
	toHadronization.push_back(*(itr->first));
	hadronizeSystems(pythia,toHadronization,all);
	}
	-
	+ */
	+ vector<pair<ColourSinglet , vector<Ropewalk::Dipole > > > funMap;
	+ for(Ropewalk::DipoleMap::iterator itr = ropewalk.begin(); itr != ropewalk.end(); ++itr)
	+ funMap.push_back(make_pair<ColourSinglet , vector<Ropewalk::Dipole > >(itr->first,itr->second));
	+ std::sort(funMap.begin(),funMap.end(),TheP8I::sorting_principle);
	+ int i = 0;
	+ for(vector<pair<ColourSinglet , vector<Ropewalk::Dipole > > >::iterator itr = funMap.begin(); itr!=funMap.end(); ++itr, ++i){
	+ if(!pythia.getRopeUserHooksPtr()->setDipoles(&(*itr),stringm0, stringR0, !average, alpha))
	+ pythia.getRopeUserHooksPtr()->setEnhancement(gloops[i].second);
	+ toHadronization.clear();
	+ toHadronization.push_back(*(itr->first));
	+ hadronizeSystems(pythia,toHadronization,all);
	+ }
	+
	break;
	}
	case 99: // New Pythia
	{
	pythia.getRopeUserHooksPtr()->setEnhancement(1.0);
	hadronizeSystems(pythia, singlets, all);
	break;
	}
	default:
	cout << "We should really not be here. This is bad..." << endl;
	}

	// fScheme = oldFS;

	// Do short analysis here:
	if(doAnalysis){
	ofstream out((analysisPath + "analysis.txt").c_str(),ios::app);
	// out << "<event>" << endl;
	// out << PrintStringsToMatlab(singlets) << endl;
	// out << "</event>" << endl;
	out.close();
	}

	}

	void StringFragmentation::dofinish(){
	if(doAnalysis!=0){
	YODA::mkWriter("yoda");
	YODA::WriterYODA::write(analysisPath + generator()->runName() + "1D.yoda",_histograms.begin(),_histograms.end());
	YODA::WriterYODA::write(analysisPath + generator()->runName() + "2D.yoda",_histograms2D.begin(),_histograms2D.end());

	}
	}

	bool StringFragmentation::
	hadronizeSystems(Pythia8Interface & pyt, const vector<ColourSinglet> & singlets, const tPVector & all) {

	TheP8EventShapes * _es = NULL;
	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 Pythia8 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 Pythia8 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, (_es ? _es->yT(p.momentum()) : p.momentum().rapidity()));
	//ymax = max(ymax, p.momentum().rapidity());
	ymin = min(ymin,(_es ? _es->yT(p.momentum()) : 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((_es ? _es->yT(p.momentum()) : p.momentum().rapidity()) ) -(_es ? _es->yT(singlets[i].partons().back()->momentum())
	: singlets[i].partons().back()->momentum().rapidity());
	//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((_es ? _es->yT(p.momentum()) : p.momentum().rapidity()) ) -(_es ? _es->yT(singlets[i].partons()[j-i]->momentum())
	: singlets[i].partons()[j-i]->momentum().rapidity());

	//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();
	pyt.errorlist();
	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();
	}
	if(window){
	// cout << stringyCut << " " << stringpTCut/GeV << endl;
	if(forcerun) forcerun = false;
	else{
	for ( int i = 1, N = event.size(); i < N; ++i ) {
	tPPtr p = pyt.getParticle(i);
	// cout << p->momentum().perp()/GeV << " " << p->rapidity() << " " << p->id() << endl;
	if (p->momentum().perp() > stringpTCut && abs(p->rapidity()) < stringyCut ){
	forcerun = true;
	return false;
	}
	}
	}
	}
	// 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 ( isnan(p->momentum().perp()/GeV) ){
	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);
	}
	}
	return true;
	}

	string StringFragmentation::PrintStringsToMatlab(vector<ColourSinglet>& singlets) {
	stringstream ss;
	vector<vector<double> > drawit;
	vector<char> names;
	for(int i=0;i<26;++i) names.push_back(char(i+97));
	vector<char>::iterator nItr = names.begin();


	for(vector<ColourSinglet>::iterator sItr = singlets.begin(); sItr!=singlets.end(); ++sItr){
	drawit.clear();
	for (tcPVector::iterator pItr = sItr->partons().begin(); pItr!=sItr->partons().end();++pItr) {
	vector<double> tmp;
	tmp.clear();
	tmp.push_back((*pItr)->rapidity());
	tmp.push_back((*pItr)->vertex().x()/femtometer);
	tmp.push_back((*pItr)->vertex().y()/femtometer);
	drawit.push_back(tmp);
	}
	ss << *nItr << " = [";
	for(int i=0, N=drawit.size();i<N;++i){
	ss << drawit[i].at(0) << ", " << drawit[i].at(1) << ", " << drawit[i].at(2) << ";" << endl;
	}
	ss << "]';\n\n" << endl;
	++nItr;
	}
	return ss.str();
	}


	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");

	if(fScheme == 4 \|\| fScheme == 5 \|\| fScheme == 1){ // We don't need multiple Pythia objects anymore!
	pythia.enableHooks();
	pythia.init(*this,theAdditionalP8Settings);
	if(pythia.version() > 0 && pythia.version() - 1.234 > 1.0){
	cout << "Pythia version: " << pythia.version() << endl;
	cout << "The chosen fragmentation scheeme requires a tweaked "
	<< "Pythia v. 1.234 for option. I will default you to "
	"option 'pythia'." << endl;
	fScheme = 0;
	}
	else{
	PytPars p;

	p.insert(make_pair<const string,double>("StringPT:sigma",theStringPT_sigma));
	p.insert(make_pair<const string,double>("StringZ:aLund",theStringZ_aLund));
	p.insert(make_pair<const string,double>("StringZ:bLund",theStringZ_bLund));
	p.insert(make_pair<const string,double>("StringFlav:probStoUD",theStringFlav_probStoUD));
	p.insert(make_pair<const string,double>("StringFlav:probSQtoQQ",theStringFlav_probSQtoQQ));
	p.insert(make_pair<const string,double>("StringFlav:probQQ1toQQ0",theStringFlav_probQQ1toQQ0));
	p.insert(make_pair<const string,double>("StringFlav:probQQtoQ",theStringFlav_probQQtoQ));

	phandler.init(fragmentationMass*fragmentationMass,baryonSuppression,p);
	pythia.getRopeUserHooksPtr()->setParameterHandler(&phandler);
	}
	}
	if( !( fScheme == 4 \|\| fScheme == 5 \|\| fScheme == 1) ){ // Not 'else' as it can change above...

	vector<string> moresettings = theAdditionalP8Settings;
	moresettings.push_back("StringPT:sigma = " + convert(theStringPT_sigma));
	moresettings.push_back("StringZ:aLund = " + convert(theStringZ_aLund));
	moresettings.push_back("StringZ:bLund = " + convert(theStringZ_bLund));
	moresettings.push_back("StringFlav:probStoUD = " +
	convert(theStringFlav_probStoUD));
	moresettings.push_back("StringFlav:probSQtoQQ = " +
	convert(theStringFlav_probSQtoQQ));
	moresettings.push_back("StringFlav:probQQ1toQQ0 = " +
	convert(theStringFlav_probQQ1toQQ0));
	moresettings.push_back("StringFlav:probQQtoQ = " +
	convert(theStringFlav_probQQtoQ));
	moresettings.push_back("OverlapStrings:fragMass = " +
	convert(fragmentationMass));
	moresettings.push_back("OverlapStrings:baryonSuppression = " + convert(baryonSuppression));

	// Initialize the OverlapPythia Handler
	if ( opHandler ) delete opHandler;
	opHandler = new OverlapPythiaHandler(this,moresettings);

	}

	// Should we do event shapes?
	if(_eventShapes) delete _eventShapes;
	_eventShapes = NULL;
	if(useThrustAxis==1){
	_eventShapes = new TheP8EventShapes();
	}

	if( doAnalysis ) {
	_histograms.push_back(new YODA::Histo1D(100,1,15,"/h_lowpT","h_lowpT"));
	_histograms.push_back(new YODA::Histo1D(100,1,15,"/h_midpT","h_midpT"));
	_histograms.push_back(new YODA::Histo1D(100,1,15,"/h_highpT","h_highpT"));
	_histograms2D.push_back(new YODA::Histo2D(10,0.,15,10,0,15,"/m_vs_pT","m_vs_pT"));
	}
	nev = 0;

	}

	void StringFragmentation::persistentOutput(PersistentOStream & os) const {
	os
	#include "StringFragmentation-output.h"
	<< fScheme << ounit(stringR0,femtometer) << ounit(stringm0,GeV) << junctionDiquark << alpha << average << stringyCut
	<< ounit(stringpTCut,GeV) << fragmentationMass << baryonSuppression
	<< oenum(window) << oenum(throwaway) << useThrustAxis << maxTries
	<< theCollapser << doAnalysis << analysisPath;

	}

	void StringFragmentation::persistentInput(PersistentIStream & is, int) {
	is
	#include "StringFragmentation-input.h"
	>> fScheme >> iunit(stringR0,femtometer) >> iunit(stringm0,GeV) >> junctionDiquark >> alpha >> average >> stringyCut
	>> iunit(stringpTCut,GeV) >> fragmentationMass >> baryonSuppression
	>> ienum(window) >> ienum(throwaway) >> useThrustAxis >> maxTries
	>> theCollapser >> doAnalysis >> analysisPath;
	}

	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 Switch<StringFragmentation,int> interfaceFragmentationScheme
	("FragmentationScheme",
	"Different options for how to handle overlapping strings.",
	&StringFragmentation::fScheme, 0, true, false);
	static SwitchOption interfaceFragmentationSchemepythia
	(interfaceFragmentationScheme,
	"pythia",
	"Plain old Pythia fragmentation",
	0);
	static SwitchOption interfaceFragmentationSchemedep1
	(interfaceFragmentationScheme,
	"dep1",
	"Not sure about this.",
	1);
	static SwitchOption interfaceFragmentationSchemedep2
	(interfaceFragmentationScheme,
	"dep2",
	"Not sure about this.",
	2);
	static SwitchOption interfaceFragmentationSchemenone
	(interfaceFragmentationScheme,
	"none",
	"Plain old Pythia fragmentation",
	0);
	static SwitchOption interfaceFragmentationSchemepipe
	(interfaceFragmentationScheme,
	"pipe",
	"Each string is enclosed by a cylinder. Effective parameters are calculated from the overlap of cylinders.",
	3);
	static SwitchOption interfaceFragmentationSchemeaverage
	(interfaceFragmentationScheme,
	"average",
	"The overlap is calculated for each dipole in all strings, and for each string the effective parameters are obtained from the average overlap.",
	4);
	static SwitchOption interfaceFragmentationSchemedipole
	(interfaceFragmentationScheme,
	"dipole",
	"Effective parameters are calculated for each breakup by determining the overlap of the corresponding dipole with other dipoles.",
	5);

	static Parameter<StringFragmentation,int> interfaceUseThrustAxis
	("UseThrustAxis",
	"Whether or not to use rapidity wrt. Thrust Axis when counting strings "
	"(put 1 or 0).",
	&StringFragmentation::useThrustAxis, 0, 0, 0,
	true, false, Interface::lowerlim);

	static Parameter<StringFragmentation,int> interfaceDoAnalysis
	("DoAnalysis",
	"Can do a short analysis where histograms of the effective parameters are "
	"plotted, and the strings in (bx,by,y)-space are printed for a couple of events. "
	"(put 1 or 0).",
	&StringFragmentation::doAnalysis, 0, 0, 0,
	true, false, Interface::lowerlim);

	static Parameter<StringFragmentation,string> interfaceAnalysisPath
	("AnalysisPath",
	"Set the system path where the (optional) analysis output should appear "
	" (directory must exist!)" ,
	&StringFragmentation::analysisPath, "");

	static Switch<StringFragmentation,bool> interfaceThrowAway
	("ThrowAway",
	"Throw away strings with weight:"
	"Use 1 - (p + q)/(m + n) for (0,-1) configuration." ,
	&StringFragmentation::throwaway, false, true, false);
	static SwitchOption interfaceThrowAwayTrue
	(interfaceThrowAway,"True","enabled.",true);
	static SwitchOption interfaceThrowAwayFalse
	(interfaceThrowAway,"False","disabled.",false);

	static Switch<StringFragmentation,bool> interfaceWindow
	("StringWindow",
	"Enable the 'window'-cut procedure, off by default."
	"Parameters will only affect if this is switched on. " ,
	&StringFragmentation::window, false, true, false);
	static SwitchOption interfaceWindowTrue
	(interfaceWindow,"True","enabled.",true);
	static SwitchOption interfaceWindowFalse
	(interfaceWindow,"False","disabled.",false);

	static Parameter<StringFragmentation,Energy> interfaceStringpTCut
	("StringpTCut",
	"No enhancement of strings with a constituent pT higher than StringpTCut (GeV), and within "
	" StringyCut (in GeV).",
	&StringFragmentation::stringpTCut, GeV, 6.0*GeV,
	0.0GeV, 0GeV,
	true, false, Interface::lowerlim);

	static Parameter<StringFragmentation,Energy> interfaceStringm0
	("Stringm0",
	".",
	&StringFragmentation::stringm0, GeV, 1.0*GeV,
	0.0GeV, 0GeV,
	true, false, Interface::lowerlim);

	static Parameter<StringFragmentation,double> interfaceJunctionDiquark
	("JunctionDiquark",
	"Suppress diquark production in breaking of junctions with this amount",
	&StringFragmentation::junctionDiquark, 0, 1.0,
	0.0, 0,
	true, false, Interface::lowerlim);

	static Parameter<StringFragmentation,double> interfaceAlpha
	("Alpha",
	"Boost the enhanced string tension with additional factor",
	&StringFragmentation::alpha, 0, 1.0,
	1.0, 0,
	true, false, Interface::lowerlim);

	static Switch<StringFragmentation,bool> interfaceAverage
	("Average",
	"Disable to try and hadronise strings with individual tensions"
	"instead of average value." ,
	&StringFragmentation::average, false, true, false);
	static SwitchOption interfaceAverageTrue
	(interfaceAverage,"True","enabled.",true);
	static SwitchOption interfaceAverageFalse
	(interfaceAverage,"False","disabled.",false);

	static Parameter<StringFragmentation,double> interfaceStringyCut
	("StringyCut",
	"No enhancement of strings with a constituent pT higher than StringpTCut (GeV), and within "
	" StringyCut.",
	&StringFragmentation::stringyCut, 0, 2.0,
	0.0, 0,
	true, false, Interface::lowerlim);

	static Parameter<StringFragmentation,Length> interfaceStringR0
	("StringR0",
	"In the string overlap model the string R0 dictates the minimum radius "
	" of a string (in fm).",
	&StringFragmentation::stringR0, femtometer, 0.5*femtometer,
	0.0femtometer, 0femtometer,
	true, false, Interface::lowerlim);

	static Parameter<StringFragmentation,double> interfaceFragmentationMass
	("FragmentationMass",
	"Set the mass used for the f(z) integral in the overlap string model "
	" default is pion mass (units GeV).",
	&StringFragmentation::fragmentationMass, 0, 0.135,
	0., 1.,
	true, false, Interface::lowerlim);

	static Parameter<StringFragmentation,double> interfaceBaryonSuppression
	("BaryonSuppression",
	"Set the fudge factor used for baryon suppression in the overlap string model. "
	" One day this should be calculated properly. This day is not today. Probably around 2...",
	&StringFragmentation::baryonSuppression, 0, 0.5,
	0.0, 1.0,
	true, false, Interface::limited);

	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);

	}

	string StringFragmentation::convert(double d) {
	ostringstream os;
	os << d;
	return os.str();
	}
	diff --git a/Hadronization/StringFragmentation.h b/Hadronization/StringFragmentation.h
	--- a/Hadronization/StringFragmentation.h
	+++ b/Hadronization/StringFragmentation.h
	@@ -1,307 +1,313 @@
	// -- 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"
	#include "OverlapPythiaHandler.h"
	#include "TheP8EventShapes.h"
	#include <sstream>
	#include "YODA/Histo1D.h"
	#include "YODA/Histo2D.h"
	#include "YODA/Writer.h"
	#include "YODA/WriterYODA.h"

	namespace TheP8I {

	+ bool sorting_principle(pair<ThePEG::ColourSinglet, vector<TheP8I::Ropewalk::Dipole> > c1, pair<ThePEG::ColourSinglet, vector<TheP8I::Ropewalk::Dipole> > c2){
	+ return (c1.first->momentum().perp() < c2.first->momentum().perp());
	+ }
	+
	+
	using namespace ThePEG;

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

	typedef map<string, double> PytPars;

	friend class Ropewalk::Dipole;

	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);

	/** @name Functions used by the persistent I/O system. */
	//@{

	/**
	* Let the given Pythia8Interface hadronize the ColourSinglet
	* systems provided.
	*/
	bool hadronizeSystems(Pythia8Interface & pyt, const vector<ColourSinglet> & singlets,
	const tPVector & all);


	void hadronizeTweak(Pythia8Interface & pyt, const ColourSinglet & singlet,
	const tPVector & all);
	/**
	* 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;

	/**
	* Internal switch for the fragmentation scheme
	*/
	int fScheme;

	/**
	* Object that handles calculation of new Pythia parameters
	*/
	ParameterHandler phandler;

	/**
	* The intrinsic string radius
	*/
	Length stringR0;

	/**
	* The assumed mass for calculating rapidities in the reeperbahn scheme
	*/
	Energy stringm0;

	/**
	* Supression of diquarks emerging from breaking junctions
	*/
	double junctionDiquark;

	double alpha;

	bool average;


	/**
	* If window is set, this determines the abs(y) window for which no enhancement will be made
	*/
	double stringyCut;

	/**
	* If window is set, this determines the pT window for which no enhancement will be made
	*/
	Energy stringpTCut;

	/**
	* Assumed m0 value in calculation of normalization of f(z) the Lund splitting function.
	*/
	double fragmentationMass;

	/**
	* Parameter surpressing baryonic content further in calculation of effective parameters
	*/
	double baryonSuppression;

	// Force hadronize systems to run to the end, regardless of windows
	bool forcerun;

	/**
	* Can provide a window in y and pT where no enhancement is made, in order to not include
	* very hard gluons in central rapidity region in a rope
	*/
	bool window;

	/**
	* Throw away some string entirely while making ropes
	*/
	bool throwaway;

	TheP8EventShapes * _eventShapes;
	/**
	* Use thrust axis to calculate rapidities; useful for LEP validation
	*/
	int useThrustAxis;

	/**
	* Additional interfaces to Pythia8 objects in case the string
	* overlap model is to be used.
	*/
	OverlapPythiaHandler * opHandler;

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

	// Keep track of how many events we already hadronized using this object, in order to clean up
	int nev;

	// Can do a small analysis
	int doAnalysis;
	vector<YODA::Histo1D*> _histograms;
	vector<YODA::Histo2D*> _histograms2D;
	string analysisPath;

	void dofinish();

	// Can print the strings of the event in a matlab friendly format
	string PrintStringsToMatlab(vector<ColourSinglet>& singlets);

	string convert(double d);
	#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: Mon, Jan 20, 8:50 PM (23 h, 46 s)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4241820
Default Alt Text: (72 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions