No OneTemporary
Actions

Size

68 KB

Subscribers

None

View Options

	diff --git a/Config/Makefile.am b/Config/Makefile.am
	--- a/Config/Makefile.am
	+++ b/Config/Makefile.am
	@@ -1,13 +1,13 @@
	mySOURCES = RndmEngine.cc Pythia8Interface.cc

	-DOCFILES = TheP8I.h RndmEngine.h Pythia8Interface.h
	+DOCFILES = TheP8I.h RndmEngine.h Pythia8Interface.h RopeUserHooks.h

	INCLUDEFILES = $(DOCFILES)

	noinst_LTLIBRARIES = libTheP8IConfig.la
	libTheP8IConfig_la_SOURCES = $(mySOURCES) $(INCLUDEFILES)

	Pythia8Interface.cc: interfaces.pl

	include $(top_srcdir)/Config/Makefile.aminclude

	diff --git a/Config/Pythia8Interface.cc b/Config/Pythia8Interface.cc
	--- a/Config/Pythia8Interface.cc
	+++ b/Config/Pythia8Interface.cc
	@@ -1,211 +1,217 @@
	// -- 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(): pythia(NULL), hooks(NULL), doHooks(false) {}

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

	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;

	}

	void Pythia8Interface::errorlist() {
	pythia->info.errorStatistics();
	}

	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;
	+ if( hooks ) delete hooks;

	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);
	+
	+ if(doHooks){
	+ hooks = new RopeUserHooks();
	+ pythia->setUserHooksPtr(hooks);
	+ }

	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/Config/Pythia8Interface.h b/Config/Pythia8Interface.h
	--- a/Config/Pythia8Interface.h
	+++ b/Config/Pythia8Interface.h
	@@ -1,196 +1,215 @@
	// -- C++ --
	+
	+
	#ifndef THEP8I_Pythia8Interface_H
	#define THEP8I_Pythia8Interface_H
	//
	// This is the declaration of the Pythia8Interface class.
	//

	#include "ThePEG/Utilities/Throw.h"
	#include "TheP8I/Config/TheP8I.h"
	#include "TheP8I/Config/RndmEngine.h"
	#include "ThePEG/Utilities/ObjectIndexer.h"
	+#include "RopeUserHooks.h"
	#include "Pythia.h"
	-
	namespace TheP8I {

	/**
	* The Pythia8Interface class is a wrapper around a static Pythia8
	* object. All communication between ThePEG and Pythia8 is handled by
	* static functions of the Pythia8Interface class.
	*/
	class Pythia8Interface: public Base {

	public:

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

	/**
	* The destructor.
	*/
	~Pythia8Interface();
	//@}

	public:

	/**
	* Initialize the main Pythia8 object setting all parameters which
	* are not default by copying them from the corresponding
	* interfaces..
	*/
	void init(const Interfaced &, const vector<string> &);

	/**
	* Set parameters which are not default in the main Pythia8 object
	* by copying them from the corresponding interfaces..
	*/
	void setParameters(const Interfaced &, const vector<string> &);

	/**
	* Access the local Pythia8 object.
	*/
	Pythia8::Pythia & operator()() {
	return *pythia;
	}

	/**
	* Access the event object of the local Pythia8 object.
	*/
	Pythia8::Event & event() {
	return pythia->event;
	}

	/**
	* Check if the local Pythia8 object has been created.
	*/
	bool created() const {
	return pythia != 0;
	}

	/**
	* Prepare for introducing a new event.
	*/
	void clearEvent() {
	event().reset();
	colourIndex.clear();
	colourIndex(0, tColinePtr());
	particleIndex.clear();
	}

	/**
	* Add a ThePEG::Particle to the Pythia8 event (if it hasn't been
	* added already) and return its index.
	*/
	int addParticle(tPPtr p, int status, int mother1, int mother2);

	/**
	* Add a ThePEG::Particle to the Pythia8 event (if it hasn't been
	* added already) and return its index.
	*/
	int addParticle(tcPPtr p, int status, int mother1, int mother2) {
	return addParticle(const_ptr_cast<tPPtr>(p), status, mother1, mother2);
	}

	/**
	* Get the colour Pythia8 colour index for a given colour
	* line. Generate a new index if the line has not been seen before.
	*/
	int addColourLine(tColinePtr c);

	/**
	* Return index of the given particle \a p.
	* @return -1 if not found.
	*/
	int indexOf(tPPtr p) {
	return particleIndex.included(p)? particleIndex(p): -1;
	}

	/**
	* Return a particle corresponding to the given entry in the Pythia8 event.
	*/
	PPtr getParticle(int idx);

	/**
	* After the event has been filled, let Pythia8 do its thing. Return
	* true if everything was OK.
	*/
	bool go();

	+ /*
	+ * Get a pointer to the UserHooks object (for ropes), pointer is NULL if not initialized
	+ */
	+ RopeUserHooks * getRopeUserHooksPtr(){
	+ return hooks;
	+ }
	+
	+ void enableHooks(){
	+ doHooks = true;
	+ }
	+
	void errorlist();
	/**
	* Retrieve the version number of Pythia.
	*/
	double version() const {
	return pythia? pythia->settings.parm("Pythia:versionNumber"): -1.0;
	}

	/**
	* Print out debugging information.
	*/
	void debug();

	/**
	* Perform all static initializations.
	*/
	static void Init();

	private:

	/**
	* The main Pythia8 object.
	*/
	Pythia8::Pythia * pythia;

	/**
	+ * Pythia 8 user hooks object for fooling around with ropes.
	+ */
	+ RopeUserHooks * hooks;
	+
	+ bool doHooks;
	+ /**
	* Association between ColourLines and colour indices.
	*/
	ObjectIndexer<int,ColourLine> colourIndex;

	/**
	* Association between Particles and indices.
	*/
	ObjectIndexer<int,Particle> particleIndex;

	/**
	* True if the static members have been initialized.
	*/
	static bool initialized;

	/**
	* The place where Pythia 8 is installed.
	*/
	static string installDir;

	/**
	* The interface to the random generatorof ThePEG.
	*/
	static RndmEngine rnd;

	/**
	* Internal Exception class.
	*/
	struct Pythia8InitException: public Exception {};

	/**
	* Internal Exception class.
	*/
	struct Pythia8ExecException: public Exception {};

	private:

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

	};

	}


	#endif /* THEP8I_Pythia8Interface_H */
	diff --git a/Config/RopeUserHooks.h b/Config/RopeUserHooks.h
	new file mode 100644
	--- /dev/null
	+++ b/Config/RopeUserHooks.h
	@@ -0,0 +1,94 @@
	+
	+#ifndef THEP8I_RopeUserHooks_H
	+#define THEP8I_RopeUserHooks_H
	+
	+#include "Pythia.h"
	+#include "TheP8I/Hadronization/ParameterHandler.h"
	+#include "TheP8I/Hadronization/Ropewalk.h"
	+
	+
	+namespace TheP8I {
	+
	+class RopeUserHooks : public Pythia8::UserHooks {
	+
	+
	+// Convenient typedef
	+typedef map<string,double> PytPars;
	+
	+ public:
	+
	+ RopeUserHooks() : _ph(NULL), _rw(NULL), _h(-1.0) {
	+ }
	+
	+ ~RopeUserHooks() {
	+ }
	+
	+ virtual bool canChangeFragPar() {return true;}
	+
	+
	+ virtual bool doChangeFragPar(Pythia8::StringFlav* flavPtr, Pythia8::StringZ* zPtr, Pythia8::StringPT * pTPtr, pair<int,int> steps, double rap) {
	+
	+ // Get new parameters
	+ PytPars newPar = fetchParameters(steps,rap);
	+ for(PytPars::iterator itr = newPar.begin(); itr!=newPar.end();++itr) settingsPtr->parm(itr->first,itr->second);
	+
	+ // debugging purposes, TODO: will go later
	+ //settingsPtr->listAll();
	+
	+ // 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){ // This is for testing, TODO: Change over a singlet
	+ _h = h;
	+ }
	+
	+ void setRopewalk(Ropewalk * rw){
	+ _rw = rw;
	+ // If we set the ropewalk, we should also use it.
	+ _h = -1.0;
	+ }
	+
	+ private:
	+
	+ PytPars fetchParameters(pair<int,int> steps, double rap){
	+
	+ // Did we remember to set the callback pointer? If we don't, we're gonna have a bad time.
	+ // Do not just construct a new one, as there will be no sensible initial values, resulting
	+ // in very subtle errors.
	+ if(!(_ph)){
	+ cout << "You failed to set the ParameterHandler callback pointer in your user-hook. Shame on you." << endl;
	+ 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);
	+
	+ if(!(_rw)){
	+ cout << "You failed to set the Ropewalk callback pointer in your user-hook. Shame on you." << endl;
	+ PytPars p;
	+ p.insert(make_pair<const string,double>("null",0.));
	+ return p;
	+ }
	+ // Make ropewalk stuff here!
	+
	+ return _ph->GetEffectiveParameters(1.0);
	+
	+ }
	+
	+ ParameterHandler * _ph;
	+ Ropewalk * _rw;
	+ double _h;
	+};
	+
	+}
	+#endif
	diff --git a/Hadronization/Makefile.am b/Hadronization/Makefile.am
	--- a/Hadronization/Makefile.am
	+++ b/Hadronization/Makefile.am
	@@ -1,21 +1,21 @@
	if PYTHIA8WORKS

	-mySOURCES = StringFragmentation.cc TheP8EventShapes.cc StringPipe.cc OverlapPythiaHandler.cc RandomAverageHandler.cc RandomHandler.cc Ropewalk.cc
	+mySOURCES = StringFragmentation.cc TheP8EventShapes.cc StringPipe.cc OverlapPythiaHandler.cc ParameterHandler.cc RandomAverageHandler.cc RandomHandler.cc Ropewalk.cc

	-DOCFILES = StringFragmentation.h TheP8EventShapes.h StringPipe.h OverlapPythiaHandler.h RandomAverageHandler.h RandomHandler.h Ropewalk.h
	+DOCFILES = StringFragmentation.h TheP8EventShapes.h StringPipe.h OverlapPythiaHandler.h ParameterHandler.h RandomAverageHandler.h RandomHandler.h Ropewalk.h

	AUTOGENERATED = StringFragmentation-var.h StringFragmentation-init.h \
	StringFragmentation-input.h StringFragmentation-output.h \
	StringFragmentation-interfaces.h

	INCLUDEFILES = $(DOCFILES) $(AUTOGENERATED)

	AM_CPPFLAGS += -I/usr/local/new/include/

	noinst_LTLIBRARIES = libTheP8IString.la

	libTheP8IString_la_SOURCES = $(mySOURCES) $(INCLUDEFILES)

	endif

	include $(top_srcdir)/Config/Makefile.aminclude
	diff --git a/Hadronization/ParameterHandler.cc b/Hadronization/ParameterHandler.cc
	new file mode 100644
	--- /dev/null
	+++ b/Hadronization/ParameterHandler.cc
	@@ -0,0 +1,134 @@
	+// -- C++ --
	+//
	+// This is the implementation of the non-inlined, non-templated member
	+// functions of the ParameterHandler class.
	+//
	+
	+#include "ParameterHandler.h"
	+
	+using namespace TheP8I;
	+
	+typedef map<string,double> PytPars;
	+
	+ParameterHandler::ParameterHandler(){
	+
	+}
	+
	+ParameterHandler::~ParameterHandler(){
	+
	+}
	+
	+void ParameterHandler::init(double m2, double bsparameter, PytPars pars){
	+ _m2 = m2;
	+ _bsparameter = bsparameter;
	+ _parameters.insert(make_pair<double,PytPars>(1.0,pars));
	+
	+ for (PytPars::iterator itr = pars.begin(); itr!=pars.end(); ++itr){
	+
	+ if(itr->first.find("StringPT:sigma")!=string::npos) sigma = itr->second;
	+
	+ else if(itr->first.find("StringZ:aLund")!=string::npos) a = itr->second;
	+
	+ else if(itr->first.find("StringZ:bLund")!=string::npos) b = itr->second;
	+
	+ else if(itr->first.find("StringFlav:probStoUD")!=string::npos) rho = itr->second;
	+
	+ else if(itr->first.find("StringFlav:probSQtoQQ")!=string::npos) x = itr->second;
	+
	+ else if(itr->first.find("StringFlav:probQQ1toQQ0")!=string::npos) y = itr->second;
	+
	+ else if(itr->first.find("StringFlav:probQQtoQ")!=string::npos) xi = itr->second;
	+
	+ else{
	+ cout << "Broken arrow!" << endl;
	+ }
	+ }
	+ // Sanity check of parameters
	+ if (rho < 0 \|\| rho > 1 \|\| x < 0 \|\| x > 1 \|\| y < 0 \|\| y > 1 \|\| xi < 0 \|\|
	+ xi > 1 \|\| sigma < 0 \|\| sigma > 1 \|\| a < 0.0 \|\| a > 2.0 \|\| b < 0.2 \|\| b > 2.0){
	+ cout << "Did you set up sensible initial Pythia 8 values? Remember:" << endl;
	+ cout << "0 < a < 2; 0.2 < b < 2; 0 < sigma < 1; 0 < xi < 1; 0 < y < 1; 0 < x < 1; 0 < rho < 1" << endl;
	+ }
	+}
	+
	+PytPars ParameterHandler::GetEffectiveParameters(double h){
	+ map<double,PytPars>::iterator it = _parameters.find(h);
	+ if(it!=_parameters.end()) return it->second;
	+
	+ if(!CalculateEffectiveParameters(h)){
	+ cout << "Something went wrong calculating effective Pythia parameters!" << endl;
	+ return _parameters.find(1.0)->second;
	+ }
	+ if(!InsertEffectiveParameters(h)){
	+ cout << "Something went wrong inserting effective Pythia parameters!" << endl;
	+ return _parameters.find(1.0)->second;
	+ }
	+ return GetEffectiveParameters(h);
	+ }
	+
	+bool ParameterHandler::InsertEffectiveParameters(double h){
	+ PytPars p;
	+ p.insert(make_pair<const string,double>("StringPT:sigma",sigma_eff));
	+ p.insert(make_pair<const string,double>("StringZ:aLund",sigma_eff));
	+ p.insert(make_pair<const string,double>("StringZ:bLund",b_eff));
	+ p.insert(make_pair<const string,double>("StringFlav:probStoUD",rho_eff));
	+ p.insert(make_pair<const string,double>("StringFlav:probSQtoQQ",x_eff));
	+ p.insert(make_pair<const string,double>("StringFlav:probQQ1toQQ0",y_eff));
	+ p.insert(make_pair<const string,double>("StringFlav:probQQtoQ",xi_eff));
	+ return (_parameters.insert(make_pair<double,PytPars>(h,p)).second);
	+}
	+
	+bool ParameterHandler::CalculateEffectiveParameters(double h){
	+ if (h <= 0) return false;
	+
	+ double hinv = 1.0 / h;
	+
	+ rho_eff = pow(rho, hinv);
	+ x_eff = pow(x, hinv);
	+ y_eff = pow(3.0 * y, hinv) / 3.0;
	+ sigma_eff = sigma * sqrt(h);
	+
	+ double C1 = _bsparameter * xi * (2 + rho) / (3 + 2xrho + 9y + 6xrhoy + xxrho*rho +
	+ 3yxxrho*rho);
	+ xi_eff = (pow(C1, hinv)/_bsparameter) * (3 + 2x_effrho_eff + 9*y_eff +
	+ 6x_effrho_effy_eff + x_effx_effrho_effrho_eff +
	+ 3y_effx_effx_effrho_eff*rho_eff) / (2 + rho_eff);
	+ if (xi_eff > 1.0) xi_eff = 1.0;
	+
	+ C1 = b / (2 + rho);
	+ b_eff = C1 * (2 + rho_eff);
	+ if (b_eff < 0.2) b_eff = 0.2;
	+ if (b_eff > 2.0) b_eff = 2.0;
	+
	+ double N = IFragmentationF(a, b);
	+ double N_eff = IFragmentationF(a, b_eff);
	+ int s = (N - N_eff) < 0 ? -1 : 1;
	+ double da = 0.1;
	+ a_eff = a - s * da;
	+ do {
	+ N_eff = IFragmentationF(a_eff, b_eff);
	+ if (((N - N_eff) < 0 ? -1 : 1) != s) {
	+ s = (N - N_eff) < 0 ? -1 : 1;
	+ da /= 10.0;
	+ }
	+ a_eff -= s * da;
	+ if (a_eff < 0.0) {a_eff = 0.0; break;}
	+ if (a_eff > 2.0) {a_eff = 2.0; break;}
	+ } while (da > 0.00001);
	+
	+
	+
	+ return true;
	+}
	+
	+double ParameterHandler::IFragmentationF(double a, double b){
	+ const int N = 100000;
	+ const double dz = 1.0 / N; //increment of z variable
	+ double Integral = 0.0;
	+ for (double z = dz; z < 1; z += dz) {
	+ //numerical integral calculation
	+ Integral += pow(1 - z, a) * exp(- b * m2 / z) / z;
	+ }
	+ return Integral * dz;
	+ }
	+
	diff --git a/Hadronization/ParameterHandler.h b/Hadronization/ParameterHandler.h
	new file mode 100644
	--- /dev/null
	+++ b/Hadronization/ParameterHandler.h
	@@ -0,0 +1,57 @@
	+// -- C++ --
	+
	+#ifndef THEP8I_ParameterHandler_H
	+#define THEP8I_ParameterHandler_H
	+
	+#include "TheP8I/Config/TheP8I.h"
	+
	+
	+namespace TheP8I {
	+
	+using namespace ThePEG;
	+
	+class ParameterHandler {
	+
	+
	+public:
	+ typedef map<string,double> PytPars;
	+
	+ /** @name Standard constructors and destructors. */
	+ //@{
	+ /**
	+ * The default constructor.
	+ */
	+ ParameterHandler();
	+
	+ /**
	+ * The destructor.
	+ */
	+ virtual ~ParameterHandler();
	+ //@}
	+
	+ void init(double m2, double bsparameter, PytPars pars);
	+
	+ PytPars GetEffectiveParameters(double h);
	+
	+protected:
	+
	+private:
	+
	+ bool InsertEffectiveParameters(double h);
	+
	+ bool CalculateEffectiveParameters(double h);
	+
	+ double IFragmentationF(double a, double b);
	+
	+ // sets of parameters, ordered in h
	+ map<double, PytPars> _parameters;
	+
	+ double a, b, rho, x, y, xi, sigma, m2;
	+ double a_eff, b_eff, rho_eff, x_eff, y_eff, xi_eff, sigma_eff;
	+ double _m2, _bsparameter;
	+};
	+
	+}
	+
	+
	+#endif
	\ No newline at end of file
	diff --git a/Hadronization/Ropewalk.h b/Hadronization/Ropewalk.h
	--- a/Hadronization/Ropewalk.h
	+++ b/Hadronization/Ropewalk.h
	@@ -1,235 +1,235 @@
	// -- C++ --
	#ifndef TheP8I_Ropewalk_H
	#define TheP8I_Ropewalk_H
	//
	// This is the declaration of the Ropewalk class.
	//

	-#include "TheP8I/Config/Pythia8Interface.h"
	+#include "TheP8I/Config/TheP8I.h"
	#include "ThePEG/EventRecord/Particle.h"
	#include "ThePEG/EventRecord/ColourSinglet.h"

	namespace TheP8I {

	using namespace ThePEG;

	/**
	* Here is the documentation of the Ropewalk class.
	*/
	class Ropewalk {

	public:

	/**
	* Container class containing information about individual dipole overlaps.
	*/
	struct Dipole {

	/**
	* Constructor.
	*/
	Dipole(tcPPtr p1, tcPPtr p2, ColourSinglet & str)
	: pc(p1), pa(p2), string(&str), n(0), m(0), p(1), q(0) {}

	/**
	* Return true if this dipole is breakable. Only dipoles between
	* gluons (not belonging to other dipoles which have broken) and
	* which has a mass above 2 GeV can break.
	*/
	bool breakable() const {
	return pc->id() == ParticleID::g && pa->id() == ParticleID::g &&
	!pc->decayed() && !pa->decayed() && s() > 4.0*GeV2;
	}

	/**
	* Calculate the probability that this dipole should break. Taking
	* into account the number of negative steps in the colour space
	* and the fraction of overlapping dipoles which are able to
	* break.
	*/
	double breakupProbability() const;

	/**
	* Set and return the effective multiplet.
	*/
	void initMultiplet();

	/**
	* Calculate the multiplicity given pp and qq;
	*/
	static double multiplicity(int pp, int qq) {
	return ( pp < 0 \|\| qq < 0 \|\| pp + qq == 0 )? 0.0:
	0.5(pp + 1)(qq + 1)*(pp + qq + 2);
	}

	/**
	* Calculate the kappa-enhancement.
	*/
	double kappaEnhancement() const {
	return 0.25(p + q + 3 - pq/double(p + q));
	}

	/**
	* Return the squared invariant mass of this dipole.
	*/
	Energy2 s() const {
	return (pc->momentum() + pa->momentum()).m2();
	}

	/**
	* Return the effective rapidityspan of this dipole.
	*/
	double yspan(Energy m0) const {
	return s() > sqr(m0)? 0.5*log(s()/sqr(m0)): 0.0;
	}

	/**
	* The coloured particle.
	*/
	tcPPtr pc;

	/**
	* The anti-coloured particle.
	*/
	tcPPtr pa;

	/**
	* The overlapping dipoles with the amount of overlap (negative
	* means anti overlap).
	*/
	vector< pair<const Dipole *, double> > overlaps;

	/**
	* The string to which this Dipole belongs.
	*/
	ColourSinglet * string;

	/**
	* The summed parallel (n) and anti-parallel (m) overlaps from
	* other dipoles.
	*/
	int n, m;

	/**
	* The multiplet.
	*/
	int p, q;

	};


	/**
	* Convenient typedef
	*/
	typedef map<ColourSinglet , vector<Dipole > > DipoleMap;

	public:

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* The default constructor.
	*/
	Ropewalk(const vector<ColourSinglet> & singlets, Length width,
	Energy scale, bool throwaway = true, bool deb = false);

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


	/**
	* Return all ColourSinglet objects with their kappa enhancement.
	*/
	vector< pair<ColourSinglet,double> > getSinglets(double DeltaY = 0.0) const;

	protected:

	/**
	* Extract all dipoles.
	*/
	void getDipoles();

	/**
	* Calculate all overlaps and initialize multiplets in all dipoles.
	*/
	void calculateOverlaps();

	/**
	* Break dipoles (and strings)
	*/
	void doBreakups();

	/**
	* Propagate a parton a finite time inthe lab-system.
	*/
	LorentzPoint propagate(const LorentzPoint & b,
	const LorentzMomentum & p) const;

	/**
	* Calculate the rapidity of a Lorentz vector but limit the
	* transverse mass to be above m0.
	*/
	double limitrap(const LorentzMomentum & p) const;

	/**
	* Return the current step.
	*/
	Step & step() const;

	/**
	* Make a copy of a ColourSinglet making sure all partons are final.
	*/
	static ColourSinglet * cloneToFinal(const ColourSinglet & cs);

	private:

	/**
	* The assumed radius of a string.
	*/
	Length R0;

	/**
	* The assumed mass for calculating rapidities
	*/
	Energy m0;

	/**
	* The vector of all Dipoles
	*/
	vector<Dipole> dipoles;

	/**
	* All Dipoles in all string
	*/
	DipoleMap stringdipoles;

	/**
	* Debug flag.
	*/
	bool debug;

	public:

	mutable double strl0;
	mutable double strl;
	mutable double avh;
	mutable double avp;
	mutable double avq;


	private:

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

	};

	}

	#endif /* TheP8I_Ropewalk_H */
	diff --git a/Hadronization/StringFragmentation.cc b/Hadronization/StringFragmentation.cc
	--- a/Hadronization/StringFragmentation.cc
	+++ b/Hadronization/StringFragmentation.cc
	@@ -1,740 +1,761 @@
	// -- 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/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	-#include <algorithm> // std::sort
	+
	using namespace TheP8I;


	StringFragmentation::StringFragmentation()
	: pythia(), stringR0(0.5femtometer), stringm0(1.0GeV), stringyCut(2.0), stringpTCut(6*GeV), fragmentationMass(0.134), baryonSuppression(2.), window(false), throwaway(false),
	useThrustAxis(0), fragmentationScheme("none"), 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());

	// Prevent funny behavior if hadronizing decays of eg. Upsilon.
	if ( newStep()->collision()->uniqueId == lastEventId ){
	secondary = true;
	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: // Pipes with full random walk, and without any further ado
	{
	RandomHandler enhancer(throwaway);
	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));

	enhancer.SetEvent(pipes);

	forcerun = false;
	// Hadronize all pipes
	for(vector<StringPipe>::iterator it = pipes.begin(); it!=pipes.end(); ++it){
	vector<ColourSinglet> toHadronization;
	double h = enhancer.KappaEnhancement(*it);
	if(h > -998.0){
	toHadronization.push_back((it).GetSingletPtr());
	if(!hadronizeSystems(*(opHandler->GetPythiaPtr(h,nev%10000==0)),toHadronization,all))
	hadronizeSystems(*(opHandler->GetPythiaPtr(1.0,nev%10000==0)),toHadronization,all);
	}
	toHadronization.clear();
	}
	break;
	}
	case 2: // Ropewalk/Dipole
	{
	Ropewalk ropewalk(singlets, stringR0, stringm0, 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: // Testing hadronization with tweaked version of Pythia.
	{
	Ropewalk ropewalk(singlets, stringR0, stringm0, throwaway, false);
	- vector< pair<ColourSinglet,double> > strings = ropewalk.getSinglets();
	+ 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);
	- hadronizeTweak(*pytp,strings[i].first,all);
	+ pythia.getRopeUserHooksPtr()->setEnhancement(strings[i].second);
	+ toHadronization.clear();
	+ toHadronization.push_back(strings[i].first);
	+ hadronizeSystems(pythia,toHadronization,all);
	}
	-
	break;
	}
	+
	default:
	cout << "We should really not be here. This is bad..." << endl;
	}


	// 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;
	}


	void StringFragmentation::
	hadronizeTweak(Pythia8Interface & pyt, const ColourSinglet & singlet, const tPVector & all) {

	Pythia8::Event & event = pyt.event();
	pyt.clearEvent();
	-
	if (singlet.nPieces() > 1) Throw<StringFragError>() << name() << " Junction string. Abort. " << Exception::runerror;
	else for(int j = 0, M = singlet.partons().size(); j < M; ++j) pyt.addParticle(singlet.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;
	CurrentGenerator::Redirect stdout2(cout, true);

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

	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");
	- 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));
	-
	- nev = 0;
	- // Initialize the OverlapPythia Handler
	- if ( opHandler ) delete opHandler;
	- opHandler = new OverlapPythiaHandler(this,moresettings);
	-
	- // Should we do event shapes?
	- _eventShapes = NULL;
	- if(useThrustAxis==1){
	- _eventShapes = new TheP8EventShapes();
	- }

	// Select hadronization scheme
	if(fragmentationScheme.find("pythia")!=string::npos \|\|
	fragmentationScheme.find("none")!=string::npos) fScheme = 0;

	else if(fragmentationScheme.find("rope")!=string::npos \|\|
	fragmentationScheme.find("test")!=string::npos) fScheme = 1;

	else if(fragmentationScheme.find("dipole")!=string::npos \|\|
	fragmentationScheme.find("ropewalk")!=string::npos) fScheme = 2;

	else if(fragmentationScheme.find("average")!=string::npos) fScheme = 3;

	- else if(fragmentationScheme.find("stringy")!=string::npos){
	- fScheme = 4;
	- cout << "You have selected hadronization with the tweaked version of Pythia!\n ThePEG needs to be configured with the tweaked Pythia 8 for this to work. Please ensure that this is the case." << endl;
	- }
	+ else if(fragmentationScheme.find("stringy")!=string::npos) fScheme = 4;
	+
	else{
	cout << "You failed to provide a hadronization scheme. I shall use Vanilla Pythia 8." << endl;
	fScheme = 0;
	}

	+ 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){ // 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 << "Reconfigure with tweaked Pythia v. 1.234 for option 'stringy'. 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){ // 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(20,1,10, "/enhancement" ,"enhancement"));
	- _histograms.push_back(new YODA::Histo1D(100,0,1., "/parameterA", "parameterA"));
	- _histograms.push_back(new YODA::Histo1D(100,0,1., "/parameterB", "parameterB"));
	- _histograms.push_back(new YODA::Histo1D(100,0,1., "/parameterRHO", "parameterRHO"));
	- _histograms.push_back(new YODA::Histo1D(100,0,1., "/parameterXI", "parameterXI"));
	- _histograms.push_back(new YODA::Histo1D(100,0,1., "/parameterX", "parameterX"));
	- _histograms.push_back(new YODA::Histo1D(100,0,1., "/parameterY", "parameterY"));
	-*/
	_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"
	<< ounit(stringR0,femtometer) << ounit(stringm0,GeV) << stringyCut << ounit(stringpTCut,GeV) << fragmentationMass << baryonSuppression << ounit(window,true) << ounit(throwaway,true) << useThrustAxis << fragmentationScheme << maxTries << theCollapser << doAnalysis << analysisPath;

	}

	void StringFragmentation::persistentInput(PersistentIStream & is, int) {
	is
	#include "StringFragmentation-input.h"
	>> iunit(stringR0,femtometer) >> iunit(stringm0,GeV) >> stringyCut >> iunit(stringpTCut,GeV) >> fragmentationMass >> baryonSuppression >> iunit(window,true) >> iunit(throwaway,true) >> useThrustAxis >> fragmentationScheme >> 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 Parameter<StringFragmentation, string> interfaceFragmentationScheme
	("FragmentationScheme",
	" Choose between hadronization schemes." ,
	&StringFragmentation::fragmentationScheme, "");

	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> 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, 2.0,
	1., 10.,
	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);

	}

	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,292 +1,299 @@
	// -- 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 {

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

	/**
	* 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;
	/**
	* Choose which fragmentation scheme to use
	*/
	string fragmentationScheme;

	/**
	* 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: Sat, May 3, 6:08 AM (1 d, 3 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4976495
Default Alt Text: (68 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions