No OneTemporary
Actions

Size

82 KB

Subscribers

None

View Options

	diff --git a/DIPSY/Emitter.cc b/DIPSY/Emitter.cc
	--- a/DIPSY/Emitter.cc
	+++ b/DIPSY/Emitter.cc
	@@ -1,942 +1,942 @@
	// -- C++ --
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the Emitter class.
	//

	//Commented

	#include "Emitter.h"
	#include "Parton.h"
	#include "ShadowParton.h"
	#include "Dipole.h"
	#include "DipoleState.h"
	#include "EffectiveParton.h"
	#include "ThePEG/Utilities/Current.h"
	#include "DipoleEventHandler.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Utilities/Throw.h"
	#include "ThePEG/Utilities/DebugItem.h"

	#include <iostream>

	#ifdef ThePEG_TEMPLATES_IN_CC_FILE
	// #include "Emitter.tcc"
	#endif

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

	#include "gsl/gsl_sf_bessel.h"

	using namespace DIPSY;

	Emitter::~Emitter() {}

	/*
	* Just a shortcut to access the input parameter from the Current handler.
	*/
	InvEnergy Emitter::rMax() const {
	return theRMax > 0.0*InvGeV? theRMax: Current<DipoleEventHandler>()->rMax();
	}

	/*
	* Shortcut to access running alpha from the Current handler.
	*/
	double Emitter::alphaS(InvEnergy r) const {
	return Current<DipoleEventHandler>()->alphaS(r);
	}

	/*
	* The veto check for the overestimates done previously.
	* This calculates the overestimated emission probability, and checks the
	* actual probability it should have been in this point, and then returns
	* true or false randomly depending on the ratio.
	* Note that this is not the final distribution, as there are still more
	* vetos done later.
	*/
	bool Emitter::OEVeto(DipolePtr dip, double y0, Parton::Point p) const {
	//set up dipole sizes.
	InvEnergy R13 = (dip->effectivePartons().first->position() - p).pt();
	InvEnergy R23 = (dip->effectivePartons().second->position() - p).pt();
	InvEnergy R12 = dip->size();

	//Too large dipoles messes up the bessel functions, and will not make it
	//through the p- veto anyways, so can be vetoed already now, to avoid
	//errors from the bessel function.
	if(R13*GeV > 100.0) return true;
	if(R23*GeV > 100.0) return true;
	if(R12*GeV > 100.0) return false;

	//some more set up.
	double bess13 = gsl_sf_bessel_K1(R13/rMax());
	double bess23 = gsl_sf_bessel_K1(R23/rMax());

	//this is the correct probability distribution used in the paper.
	Energy2 correctDist = alphaBar(min(min(R13,R23),R12))/(M_PI2.0sqr(rMax()))*
	(sqr(bess13) + sqr(bess23) - (sqr(R13) + sqr(R23) - sqr(R12))/
	(R13R23)bess13*bess23);

	// If we create dipoles which are larger than the original, the
	// original partons will be distributed as dpt2/pt4. Here we may
	// include an extra fudge factor to emulate a matrix element
	// correction.
	if ( R13 > R12 && R23 > R12 && Current<DipoleEventHandler>()->fudgeME() )
	correctDist *= 1.0 - 1.0/(1.0 + cosh(dip->partons().first->y() -
	dip->partons().second->y()));

	//Now calculate the overestimated distribution in the same way as it is done in gerateY() and generateXT().
	double Coe = alphaBar(R12/2.0)/M_PI*
	sqr(gsl_sf_bessel_K1(R12/(2.0rMax())))
	sqr(R12/(2.0rMax()))
	(R12/(2.0*theRScale)+2.0);
	if(R12>2.0*rMax()) //double check that this is ok!!!!
	Coe = sqrt(R12/(2.0rMax()))*
	exp(R12/(rMax())-2.0);

	//this is the overestime
	Energy2 overestimate = Coe(1.0/(sqr(R13)(R13/theRScale+2.0)) +
	1.0/(sqr(R23)*(R23/theRScale+2.0)));

	//if the overestimate is ok, this should never happen.
	if( correctDist/overestimate > 1.0 )
	Throw<EmitterException>()
	<< "In DIPSY::Emitter " << name() << ": XT keep prob is " << correctDist/overestimate
	<< " (r12 = " << R12GeV << ", Rr13 = " << R13GeV << ", r23 = " << R23*GeV << ")"
	<< Exception::warning;

	//generate a random number and check if veto or not.
	if(correctDist/overestimate>UseRandom::rnd())
	return false;
	else
	return true;
	}

	/*
	* The veto check for the overestimates done previously using shadows.
	* This calculates the overestimated emission probability, and checks the
	* actual probability it should have been in this point, and then returns
	* true or false randomly depending on the ratio.
	* Note that this is not the final distribution, as there are still more
	* vetos done later.
	*/
	bool Emitter::OEVeto(DipolePtr dip, tSPartonPtr sp1, tSPartonPtr sp2,
	double y0, Parton::Point p) const {
	//set up dipole sizes.
	InvEnergy R13 = (dip->partons().first->position() - p).pt();
	InvEnergy R23 = (dip->partons().second->position() - p).pt();
	InvEnergy R12 = dip->size();

	//Too large dipoles messes up the bessel functions, and will not make it
	//through the p- veto anyways, so can be vetoed already now, to avoid
	//errors from the bessel function.
	if ( R13GeV > 100.0 \|\| R23GeV > 100.0 ) return true;
	if ( R12*GeV > 100.0 ) return false;

	//some more set up.
	double bess13 = gsl_sf_bessel_K1(R13/rMax());
	double bess23 = gsl_sf_bessel_K1(R23/rMax());

	//this is the correct probability distribution used in the paper.
	Energy2 correctDist = alphaBar(min(min(R13,R23),R12))/(M_PI2.0sqr(rMax()))*
	(sqr(bess13) + sqr(bess23) - (sqr(R13) + sqr(R23) - sqr(R12))/
	(R13R23)bess13*bess23);

	// If we create dipoles which are larger than the original, the
	// original partons will be distributed as dpt2/pt4. Here we may
	// include an extra fudge factor to emulate a matrix element
	// correction.
	if ( R13 > R12 && R23 > R12 && Current<DipoleEventHandler>()->fudgeME() )
	correctDist *= 1.0 - 1.0/(1.0 + cosh(dip->partons().first->y() -
	dip->partons().second->y()));

	//Now calculate the overestimated distribution in the same way as it is done in gerateY() and generateXT().
	double Coe = alphaBar(R12/2.0)/M_PI*
	sqr(gsl_sf_bessel_K1(R12/(2.0rMax())))
	sqr(R12/(2.0rMax()))
	(R12/(2.0*theRScale)+2.0);
	if(R12>2.0*rMax()) //double check that this is ok!!!!
	Coe = sqrt(R12/(2.0rMax()))*
	exp(R12/(rMax())-2.0);

	//this is the overestime
	Energy2 overestimate = Coe(1.0/(sqr(R13)(R13/theRScale+2.0)) +
	1.0/(sqr(R23)*(R23/theRScale+2.0)));

	//if the overestimate is ok, this should never happen.
	//TODO: write a proper error message.
	if ( correctDist/overestimate > 1.0 )
	Throw<EmitterException>()
	<< "In DIPSY::Emitter " << name() << ": XT keep prob is " << correctDist/overestimate
	<< " (r12 = " << R12GeV << ", Rr13 = " << R13GeV << ", r23 = " << R23*GeV << ")"
	<< Exception::warning;

	//generate a random number and check if veto or not.
	return !UseRandom::rndbool(correctDist/overestimate);

	}

	/*
	* This corrects the overestimated emission probability in Y used in generateY().
	* The correct emission rate is calculated for the rapidity in question and
	* compared to the overestimate rate used in generateY().
	* Note that this is not the final correct y-distribution, as there are still several
	* vetos to be checked, it is just a step closer to the correct distribution.
	*/
	bool Emitter::YVeto(double y,DipolePtr dip, double Coe, double rateOE) const {
	tEffectivePartonPtr p1 = dip->effectivePartons().first;
	tEffectivePartonPtr p2 = dip->effectivePartons().second;

	//calculate the correct rate at the rapidity y in question.
	//TODO: member calculating cutoff, interface before/after recoil.
	//Before recoil corresponds to order wrt propagator, after wrt final gluon.
	InvEnergy cutoff = 0.990.5
	min(pTScale()/(p1->plus()*exp(y) + p1->pT().pt()/2.0),
	min(pTScale()/(p2->plus()*exp(y) + p2->pT().pt()/2.0),
	dip->size()/4.0)); //if ordered AFTER recoil
	// InvEnergy cutoff = 0.99*
	// min(2.0/3.0pTScale()exp(-y)/p1->plus(),
	// min(2.0/3.0pTScale()exp(-y)/p2->plus(),
	// dip->size()/4.0)); //if ordered BEFORE recoil

	//calculate overestimated rate according the what is used in generateY().
	double rate = 2.0M_PICoelog(1.0 + 2.0theRScale/cutoff);

	//generate rnd() and veto.
	if(rate/rateOE > UseRandom::rnd()) //if exp(-...) big enough
	return false; //keep
	else
	return true; //else veto
	}

	/*
	* This corrects the overestimated emission probability in Y used in generateY().
	* The correct emission rate is calculated for the rapidity in question and
	* compared to the overestimate rate used in generateY().
	* Note that this is not the final correct y-distribution, as there are still several
	* vetos to be checked, it is just a step closer to the correct distribution.
	*/
	bool Emitter::YVeto(double y, DipolePtr dip, tSPartonPtr sp1, tSPartonPtr sp2,
	double Coe, double rateOE) const {
	//calculate the correct rate at the rapidity y in question.
	//TODO: member calculating cutoff, interface before/after recoil.
	//Before recoil corresponds to order wrt propagator, after wrt final gluon.
	InvEnergy cutoff = 0.990.5
	- min(pTScale()/(sp1->plus()*exp(y) + sp1->pT().pt()/2.0),
	- min(pTScale()/(sp2->plus()*exp(y) + sp2->pT().pt()/2.0),
	+ min(pTScale()/(sp1->plus0()*exp(y) + sp1->pT0().pt()/2.0),
	+ min(pTScale()/(sp2->plus0()*exp(y) + sp2->pT0().pt()/2.0),
	dip->size()/4.0));

	//calculate overestimated rate according the what is used in generateY().
	double rate = 2.0M_PICoelog(1.0 + 2.0theRScale/cutoff);

	return !UseRandom::rndbool(rate/rateOE);

	}

	/*
	* Generates a rapidity for the next emission.
	* The distribution f(y) in this member is using is an overestimate, and will be
	* vetoed down to correct distribution in the checks of the x_T distribution.
	*
	* Already the overestimate f(y) is too messy to generate directly though,
	* so it has to be done through a further overestimate g(y) > f(y) which is
	* then vetoed (only depending on y) down to f(y) with YVeto().
	* Look at the writeup for more details on f() and g().
	*/
	double Emitter::
	generateY(DipolePtr dip, double ymin, double ymax) const {
	//set up shorter names.
	tEffectivePartonPtr p1 = dip->effectivePartons().first;
	tEffectivePartonPtr p2 = dip->effectivePartons().second;
	InvEnergy r = dip->size();

	//To large dipoles will mess up the bessel function, so just shut them
	// off by emitting after the maximum y.
	if ( r*GeV > 100.0 ) return ymax + 1.0;

	//Calculate the cutoff in transverse distance. See writeup for details.
	InvEnergy cutoff = 0.990.5
	min(pTScale()/(p1->plus()*exp(ymin) + p1->pT().pt()/2.0),
	min(pTScale()/(p2->plus()*exp(ymin) + p2->pT().pt()/2.0),
	dip->size()/4.0)); //if ordered AFTER recoil
	// InvEnergy cutoff = 0.99*
	// min(2.0/3.0pTScale()exp(-ymin)/p1->plus(),
	// min(2.0/3.0pTScale()exp(-ymin)/p2->plus(),
	// dip->size()/4.0)); //if ordered BEFORE recoil

	//Calculate the overestimated coefficient. See writeup for details.
	double Coe = alphaBar(r/2.0)/M_PI*
	sqr(gsl_sf_bessel_K1(r/(2.0rMax())))
	sqr(r/(2.0rMax()))
	(r/(2.0*theRScale)+2.0);

	//For assymptotically large dipoles, the overestimated emission rate grows
	//faster than Coe, so add an extra factor for large dipoles.
	//this will have to be removed in the veto.
	if(r>2.0*rMax()) //double check that this is ok!!!
	Coe = sqrt(r/(2.0rMax()))*
	exp(r/(rMax())-2.0);

	//calculate the overestimated rate of emission (per unit rapidity).
	double rateOE = 2.0M_PICoelog(1.0 + 2.0theRScale/cutoff);

	//generate the rapidity. (see writeup for details.)
	double y = ymin + log(1.0-log(UseRandom::rnd())/rateOE);

	//veto to get to correct distribution, recur if vetoed.
	if(YVeto(y,dip,Coe,rateOE*exp(y-ymin)) && y < ymax)
	return generateY(dip, y, ymax);
	else{
	return y;
	}
	}

	/*
	* Generates a rapidity for the next emission using shadows.
	* The distribution f(y) in this member is using is an overestimate, and will be
	* vetoed down to correct distribution in the checks of the x_T distribution.
	*
	* Already the overestimate f(y) is too messy to generate directly though,
	* so it has to be done through a further overestimate g(y) > f(y) which is
	* then vetoed (only depending on y) down to f(y) with YVeto().
	* Look at the writeup for more details on f() and g().
	*/
	double Emitter::generateY(DipolePtr dip, tSPartonPtr sp1, tSPartonPtr sp2,
	double ymin, double ymax) const {
	//set up shorter names.
	InvEnergy r = dip->size();

	//To large dipoles will mess up the bessel function, so just shut them
	// off by emitting after the maximum y.
	if ( r*GeV > 100.0 ) return ymax + 1.0;

	//Calculate the cutoff in transverse distance. See writeup for details.
	InvEnergy cutoff = 0.990.5
	- min(pTScale()/(sp1->plus()*exp(ymin) + sp1->pT().pt()/2.0),
	- min(pTScale()/(sp2->plus()*exp(ymin) + sp2->pT().pt()/2.0), r/4.0));
	+ min(pTScale()/(sp1->plus0()*exp(ymin) + sp1->pT0().pt()/2.0),
	+ min(pTScale()/(sp2->plus0()*exp(ymin) + sp2->pT0().pt()/2.0), r/4.0));

	//Calculate the overestimated coefficient. See writeup for details.
	double Coe = alphaBar(r/2.0)/M_PI*
	sqr(gsl_sf_bessel_K1(r/(2.0rMax())))
	sqr(r/(2.0rMax()))
	(r/(2.0*theRScale)+2.0);

	//For assymptotically large dipoles, the overestimated emission rate grows
	//faster than Coe, so add an extra factor for large dipoles.
	//this will have to be removed in the veto.
	if(r>2.0*rMax()) //double check that this is ok!!!
	Coe = sqrt(r/(2.0rMax()))*
	exp(r/(rMax())-2.0);

	//calculate the overestimated rate of emission (per unit rapidity).
	double rateOE = 2.0M_PICoelog(1.0 + 2.0theRScale/cutoff);

	//generate the rapidity. (see writeup for details.)
	double y = ymin + log(1.0-log(UseRandom::rnd())/rateOE);

	//veto to get to correct distribution, recur if vetoed.
	if( YVeto(y, dip, sp1, sp2, Coe, rateOE*exp(y-ymin)) && y < ymax )
	return generateY(dip, sp1, sp2, y, ymax);

	return y;

	}

	/*
	* Generate the transverse position. First generate the distance r
	* larger than the cutoff used in generateY(), then decide around which
	* parton, and then the angle. Distributions found in writeup.
	*/
	Parton::Point Emitter::generateXT(DipolePtr dip, double y0) const {
	//set up
	tEffectivePartonPtr p1 = dip->effectivePartons().first;
	tEffectivePartonPtr p2 = dip->effectivePartons().second;

	//Calculate the cutoff.
	//TODO: make cutoff member, taking dipole and y as arguments.
	InvEnergy cutoff = 0.990.5
	min(pTScale()/(p1->plus()*exp(y0) + p1->pT().pt()/2.0),
	min(pTScale()/(p2->plus()*exp(y0) + p2->pT().pt()/2.0),
	dip->size()/4.0)); //if ordered AFTER recoil
	// InvEnergy cutoff = 0.99*
	// min(2.0/3.0pTScale()exp(-y0)/p1->plus(),
	// min(2.0/3.0pTScale()exp(-y0)/p2->plus(),
	// dip->size()/4.0)); //if ordered BEFORE recoil

	//The normalisation. See writeup for details.
	double C2 = 2.0/log(1.0+2.0*theRScale/cutoff); //normalises g(z)

	//Generate the distance.
	InvEnergy r = 2.0theRScale/(exp(2.0UseRandom::rnd()/C2)-1.0);
	//tested to give correct distribution C2/(r^3/theRScale+2r^2)

	//generate the angle.
	double phi = 2.0M_PIUseRandom::rnd();

	//decide which parton, and create the Point to be returned.
	if(UseRandom::rnd()>0.5)
	return Parton::Point(p1->position().first + cos(phi)*r,
	p1->position().second + sin(phi)*r);
	else
	return Parton::Point(p2->position().first + cos(phi)*r,
	p2->position().second + sin(phi)*r);
	}

	/*
	* Generate the transverse position. First generate the distance r
	* larger than the cutoff used in generateY(), then decide around which
	* parton, and then the angle. Distributions found in writeup.
	*/
	Parton::Point
	Emitter::generateXT(DipolePtr dip, tSPartonPtr sp1, tSPartonPtr sp2, double y0) const {

	//Calculate the cutoff.
	//TODO: make cutoff member, taking dipole and y as arguments.
	InvEnergy cutoff = 0.990.5
	- min(pTScale()/(sp1->plus()*exp(y0) + sp1->pT().pt()/2.0),
	- min(pTScale()/(sp2->plus()*exp(y0) + sp2->pT().pt()/2.0),
	+ min(pTScale()/(sp1->plus0()*exp(y0) + sp1->pT0().pt()/2.0),
	+ min(pTScale()/(sp2->plus0()*exp(y0) + sp2->pT0().pt()/2.0),
	dip->size()/4.0)); //if ordered AFTER recoil

	//The normalisation. See writeup for details.
	double C2 = 2.0/log(1.0 + 2.0*theRScale/cutoff); //normalises g(z)

	//Generate the distance.
	InvEnergy r = 2.0theRScale/(exp(2.0UseRandom::rnd()/C2)-1.0);
	//tested to give correct distribution C2/(r^3/theRScale+2r^2)

	//generate the angle.
	double phi = 2.0M_PIUseRandom::rnd();

	//decide which parton, and create the Point to be returned.
	tPartonPtr pp =
	UseRandom::rndbool()? dip->partons().first: dip->partons().second;
	return Point(pp->position().first + cos(phi)*r,
	pp->position().second + sin(phi)*r);

	}

	/*
	* the main function to be called from outside. Other functions are mainly to help
	* this one. To get the emission right, several layers of overestimates and vetoes
	* are used.
	*/
	void Emitter::generate(Dipole & dip, double ymin, double ymax) const {

	if ( dip.partons().first->shadow() ) {
	generateWithShadows(dip, ymin, ymax);
	return;
	}

	//Lowest allowed emission rapidity. Not before any of the partons
	//in the dipole, and not before the supplied ymin.
	double y0 = max(ymin,max(dip.partons().first->y(),dip.partons().second->y()));

	//define the resolved effective partons, as function of the dipole size.
	//Later, when the transverse position is decided, this will have to
	//be recalculated with the actual resolution ranges.
	//This is an overestimate, as later lowering the range can only
	//decrease the allowed phase space through less p+, and increased p_T.
	dip.effectivePartons(EffectiveParton::create(*(dip.partons().first),
	dip.size()/2.0),
	EffectiveParton::create(*(dip.partons().second),
	dip.size()/2.0));

	//Some renaming to save space later.
	tEffectivePartonPtr p1 = dip.effectivePartons().first;
	tEffectivePartonPtr p2 = dip.effectivePartons().second;

	//We will go through the while loop until something passes all the vetoes
	//or pass the maximum allowed rapidity ymax. However, put a limit
	//on the number of trials not to get stuck. So define the count to keep
	//track. Possibly a for loop would be more appropriate at this point...
	while ( true ) {
	//reset the range to the overestimate, as we do not know
	//what x_T will eb chosen this trial.
	p1->setRange( dip.size()/2.0 );
	p2->setRange( dip.size()/2.0 );

	//generate a rapidity, according to an overestimate.
	y0 = generateY(& dip, y0, ymax);


	//if the generated rapidity is above the limit, return an empty pointer
	if ( y0 > ymax ) {
	dip.generatedGluon(new_ptr(Parton()));
	dip.generatedY(y0);
	break;
	}

	//generate a transverse position, using an overestimated
	//phase space.
	Parton::Point p = generateXT(& dip, y0);

	//Bring down the overestimate in the distribution
	//to get the right amplitude. After this, only phase space limitations
	//left.
	if(OEVeto(& dip,y0,p)) {
	continue;
	}

	//Set up notation.
	InvEnergy r13 = (p1->position() - p).pt();
	InvEnergy r23 = (p2->position() - p).pt();
	InvEnergy r12 = dip.size();

	//this is to what extent the first parent emitted, and to what extent the
	//second one. Here, it is some kind of superposition where they are emitting
	//together, so both parents recoil a bit, etc.
	double P1 = sqr(r23)/(sqr(r23) + sqr(r13));
	double P2 = sqr(r13)/(sqr(r13) + sqr(r23));

	//check how much is actually resolved, small emissions resolve more
	if ( Current<DipoleEventHandler>()->emitter().rangeMode() != 1 ) {
	if ( r13 < r12/2.0 )
	p1->setRange( r13 );
	if ( r23 < r12/2.0 )
	p2->setRange( r23 );
	}

	//The transverse momentum of the new parton from the parents.
	TransverseMomentum v1 = pTScale()*(p - p1->position())/sqr(r13);
	TransverseMomentum v2 = pTScale()*(p - p2->position())/sqr(r23);

	//calculate the 4-momentum of the emitted gluon.
	TransverseMomentum pt = v1 + v2;
	Energy pplus = pt.pt()*exp(-y0);
	Energy pminus = pt.pt2()/pplus;

	//check that there's enough p+
	if ( p1->plus() - P1pplus <= 0.0GeV ){
	continue;
	}
	if ( p2->plus() - P2pplus <= 0.0GeV ){
	continue;
	}

	//calculate the new 4-momenta of the recoiled parents.
	TransverseMomentum pt1 = p1->pT() - v1;
	TransverseMomentum pt2 = p2->pT() - v2;
	Energy plus1 = p1->plus() - P1*pplus;
	Energy plus2 = p2->plus() - P2*pplus;

	double y1 = log(pt1.pt()/plus1);
	double y2 = log(pt2.pt()/plus2);
	Energy minus1 = pt1.pt2()/plus1;
	Energy minus2 = pt2.pt2()/plus2;

	//this option use the p- ordering of the non-effective parton,
	//so p- should be calculated from the non-effective parton.
	if ( theMinusOrderingMode == 1 ) {
	minus1 = pt1.pt()*exp(dip.partons().first->y());
	minus2 = pt2.pt()*exp(dip.partons().second->y());
	}

	//in the first option, p- is required to be fully ordered with both parents.
	//in the second option, the parton has to be ordered mainly to the parton that
	//emitted it most, so it is allowed to be unordered by a factor P1 and P2.
	//PSInflation is a plain factor that the emission is allowed to be unordered by.
	if ( bothOrderedEvo() ) {
	if ( pminusthePSInflation < max(minus1, minus2)thePMinusOrdering ) continue;
	}
	else {
	if ( pminusthePSInflation < max(P1minus1, P2minus2)thePMinusOrdering ) continue;
	}

	//check rapidity ordered with the two recoiled parents.
	if ( y1 > y0 ) {
	continue;
	}
	if ( y2 > y0 ) {
	continue;
	}

	//check that none of the internal partons recoiled over ymax.
	if ( p1->recoilsOverYMax( v1, P1*pplus, ymax ) ) {
	continue;
	}
	if ( p2->recoilsOverYMax( v2, P2*pplus, ymax ) ) {
	continue;
	}

	//This is a self-consistency check.
	//If things are donw correctly, then emissions at the limit of the
	//cutoff in x_T used in generateXT() and generateY() should all be vetoed.
	//if emissions at the limit go through, it means that we are missing part of
	//the phase space just within the cutoff, which is bad.
	//
	//to check this, the cutoff is decreased by a factor 0.99 in generateY() etc above,
	//to create the occasional emission in a region that should always be vetoed if the
	//overestimates are fair. This is a check at the end so that none of these emission
	//between 1 and 0.99 of the cutoff went through.
	//TODO: match this with cutoff above
	//TODO: proper error message.
	if ( min(r13, r23) < min(2.0/3.0pTScale()/(p1->plus()exp(y0)),
	min(2.0/3.0pTScale()/(p2->plus()exp(y0)),
	dip.size()/4.0)) )
	Throw<EmitterException>()
	<< "In DIPSY::Emitter " << name() << ": Emission below r-cutoff passed vetos!"
	<< " (r12 = " << ounit(r12, InvGeV) << ", r13 = " << ounit(r13, InvGeV)
	<< ", r23 = " << ounit(r23, InvGeV)
	<< "v1/pt = " << double(min(v1.pt(), v2.pt())/pt.pt()) << ")"
	<< Exception::warning;

	// passed. Set up a new parton for the dipole with the information of the emission.
	PartonPtr gluon = new_ptr(Parton());
	gluon->position(p);
	dip.generatedGluon(gluon);
	dip.generatedY(y0);

	//leave the while loop.
	break;
	}
	}

	/*
	* the main function to be called from outside when using
	* shadows. Other functions are mainly to help this one. To get the
	* emission right, several layers of overestimates and vetoes are used.
	*/
	void Emitter::generateWithShadows(Dipole & dip, double ymin, double ymax) const {

	static DebugItem trace("DIPSY::Trace", 9);

	if ( trace ) cerr << "Gen " << dip.tag() << endl;

	// Handy pointers to the partons.
	tPartonPtr p1 = dip.partons().first;
	tPartonPtr p2 = dip.partons().second;

	//First check how far down the history we must consider previous shadows.
	tSPartonPtr sp1 =
	dip.partons().first->shadow()->resolve(dip.size2()/4.0, dip.partons().second);
	tSPartonPtr sp2 =
	dip.partons().second->shadow()->resolve(dip.size2()/4.0, dip.partons().first);

	//Lowest allowed emission rapidity. Not before any of the partons
	//in the dipole, and not before the supplied ymin.
	double y0 = max(ymin, max(p1->y(), p2->y()));

	//We will go through the while loop until something passes all the vetoes
	//or pass the maximum allowed rapidity ymax. However, put a limit
	//on the number of trials not to get stuck. So define the count to keep
	//track. Possibly a for loop would be more appropriate at this point...
	while ( true ) {
	//generate a rapidity, according to an overestimate.
	y0 = generateY(&dip, sp1, sp2, y0, ymax);


	//if the generated rapidity is above the limit, return an empty pointer
	if ( y0 > ymax ) {
	dip.generatedGluon(new_ptr(Parton()));
	dip.generatedY(y0);
	return;
	}

	//generate a transverse position, using an overestimated
	//phase space.
	Parton::Point p = generateXT(&dip, sp1, sp2, y0);

	//Bring down the overestimate in the distribution
	//to get the right amplitude. After this, only phase space limitations
	//left.
	if ( OEVeto(&dip, sp1, sp2, y0, p) ) {
	continue;
	}

	//Set up notation.
	InvEnergy2 r13 = (p1->position() - p).pt2();
	InvEnergy2 r23 = (p2->position() - p).pt2();
	InvEnergy2 r12 = dip.size2();

	// We need to already here decide which parton is emitting.
	bool firstsplit = UseRandom::rndbool(r23/(r23 + r13));
	tPartonPtr pe = firstsplit? p1: p2;
	tPartonPtr pr = firstsplit? p2: p1;
	tSPartonPtr spe = firstsplit? sp1: sp2;
	InvEnergy2 re3 = firstsplit? r13: r23;

	// The emitting parton may be further resolved by the emission.
	tSPartonPtr spre = re3 < r12/4.0? pe->shadow()->resolve(re3, pr): spe;

	//The transverse momentum of the new parton from the emitter.
	TransverseMomentum pt = pTScale()*(p - pe->position())/re3;
	Energy pplus = pt.pt()*exp(-y0);
	Energy pminus = pt.pt2()/pplus;

	//check that there's enough p+
	- if ( pplus > spre->plus() ) continue;
	+ if ( pplus > spre->plus0() ) continue;

	//calculate the new 4-momenta of the recoiled parents.
	TransverseMomentum pte = spre->pT() - pt;
	- Energy pluse = spre->plus() - pplus;
	+ Energy pluse = spre->plus0() - pplus;

	double ye = log(pte.pt()/pluse);
	Energy minuse = pte.pt2()/pluse;
	if ( theMinusOrderingMode == 2 ) {
	ShadowParton::Propagator prop =
	spre->propagator(min(re3, r12/4.0), pr, -1);
	if ( prop.fail ) continue;
	pluse = prop.p.plus() - pplus;
	if ( pluse <= ZERO ) continue;
	pte = TransverseMomentum(prop.p) - pt;
	ye = log(pte.pt()/pluse);
	minuse = pte.pt2()/pluse;

	// * TODO * fix masses!
	}
	else if ( theMinusOrderingMode == 1 )
	minuse = pte.pt()*exp(pe->y());

	// Ensure approximate minus ordering.
	if ( pminusthePSInflation < minusethePMinusOrdering ) continue;

	//check rapidity ordered with the two recoiled parents.
	if ( ye > y0 ) continue;

	// passed. Set up a new parton for the dipole with the information of the emission.
	PartonPtr gluon = new_ptr(Parton());
	gluon->position(p);
	gluon->mainParent(pe);
	dip.generatedGluon(gluon);
	dip.generatedY(y0);

	// And we're done!
	return;
	}
	}

	/*
	* Maked the dipole actually emit the dipole prepared in generate(), setting up
	* parents, children, changing 4-momenta etc.
	*/
	void Emitter::emit(Dipole & d) const {
	//some notation.
	if ( d.partons().first->shadow() ) {
	emitWithShadows(d);
	return;
	}
	PartonPtr p = d.generatedGluon();
	InvEnergy r13 = (d.partons().first->position() - p->position()).pt();
	InvEnergy r23 = (d.partons().second->position() - p->position()).pt();

	//The recoils.
	//TODO: use recoil().
	TransverseMomentum v1 = pTScale()*(p->position() - d.partons().first->position())/sqr(r13);
	TransverseMomentum v2 = pTScale()*(p->position() - d.partons().second->position())/sqr(r23);

	//the ratio each parent emitted the gluon with.
	double P1 = sqr(r23)/(sqr(r13) + sqr(r23));
	double P2 = 1.0 - P1;

	//set 4-momenta for emission.
	p->pT(v1 + v2);
	p->y(d.generatedY());
	p->plus(p->pT().pt()*exp(-d.generatedY()));
	p->minus(p->pT().pt()*exp(d.generatedY()));
	p->oY(p->y());

	//change the dipole structure, generate new colour etc.
	d.splitDipole(P2);

	//perform the recoil on the effective partons.
	//the resolution ranges should be set since the generation of the emission.
	d.effectivePartons().first->recoil( v1, P1*p->plus() );
	d.effectivePartons().second->recoil( v2, P2*p->plus() );
	}

	void Emitter::emitWithShadows(Dipole & d) const {
	static DebugItem trace("DIPSY::Trace", 9);

	if ( trace ) cerr << "Emit " << d.tag();

	//some notation.
	PartonPtr p = d.generatedGluon();

	// Decide which parton is the emitter.
	bool em1 = ( p->mainParent() == d.partons().first );
	tPartonPtr emitter = em1? d.partons().first: d.partons().second;
	tPartonPtr recoiler = em1? d.partons().second: d.partons().first;
	InvEnergy2 res = emitter->dist2(*p);
	tSPartonPtr sem = emitter->shadow()->resolve(min(res, d.size2()/4.0), recoiler);

	TransverseMomentum rec =
	pTScale()*(p->position() - emitter->position())/res;

	//change the dipole structure, generate new colour etc.
	d.splitDipole(em1? 0.0: 1.0);

	//set 4-momenta for emission.
	p->pT(rec);
	p->y(d.generatedY());
	p->plus(p->pT().pt()*exp(-d.generatedY()));
	p->minus(p->pT().pt()*exp(d.generatedY()));
	p->oY(p->y());

	if ( theMinusOrderingMode == 2 ) {
	ShadowParton::Propagator prop =
	sem->propagator(min(res, d.size2()/4.0), recoiler, -1);
	emitter->pT(TransverseMomentum(prop.p) - rec);
	emitter->plus(prop.p.plus() - p->plus());
	emitter->y(log(emitter->mt()/emitter->plus()));
	emitter->minus(emitter->mt2()/emitter->plus());
	} else {
	emitter->pT(sem->pT() - rec);
	emitter->plus(sem->plus() - p->plus());
	emitter->y(log(emitter->pT().pt()/emitter->plus()));
	emitter->minus(emitter->mt2()/emitter->plus());
	}

	emitter->shadow()->setupEmission(emitter, p, *recoiler);

	if ( trace ) cerr << " -> " << d.children().first->tag()
	<< " + " << d.children().second->tag() << endl;

	}

	void Emitter::persistentOutput(PersistentOStream & os) const {
	os << ounit(thePSInflation, 1.0)
	<< ounit(thePMinusOrdering, 1.0)
	<< ounit(theRScale, InvGeV)
	<< ounit(theRMax, InvGeV)
	<< ounit(thePTScale, 1.0)
	<< ounit(theBothOrderedEvo, true)
	<< ounit(theBothOrderedInt, true)
	<< ounit(theBothOrderedFS, true)
	<< ounit(theRangeMode, 1)
	<< ounit(theMinusOrderingMode, 1);
	}

	void Emitter::persistentInput(PersistentIStream & is, int) {
	is >> iunit(thePSInflation, 1.0)
	>> iunit(thePMinusOrdering, 1.0)
	>> iunit(theRScale, InvGeV)
	>> iunit(theRMax, InvGeV)
	>> iunit(thePTScale, 1.0)
	>> iunit(theBothOrderedEvo, true)
	>> iunit(theBothOrderedInt, true)
	>> iunit(theBothOrderedFS, true)
	>> iunit(theRangeMode, 1)
	>> iunit(theMinusOrderingMode, 1);
	}


	// Static variable needed for the type description system in ThePEG.
	#include "ThePEG/Utilities/DescribeClass.h"
	DescribeClass<Emitter,HandlerBase>
	describeDIPSYEmitter("DIPSY::Emitter", "libAriadne5.so libDIPSY.so");

	void Emitter::Init() {

	static ClassDocumentation<Emitter> documentation
	("The Emitter class is responsible for generating and performing "
	"emissions from dipoles. This base class does the default emission "
	"strategy.");


	static Parameter<Emitter,InvEnergy> interfaceRScale
	("RScale",
	"Constant used in overestimate of emission probability.",
	&Emitter::theRScale, InvGeV, 1.0InvGeV, 0.0InvGeV, 0.0*InvGeV,
	true, false, Interface::lowerlim);

	static Parameter<Emitter,double> interfacePSInflation
	("PSInflation",
	"How much p+- ordering should be enforced. 1 is normal ordering, high numbers "
	"are no ordering(energy must still be conserved though), "
	"low number demands stronger ordering",
	&Emitter::thePSInflation, 1.0, 1.0, 0.0, 0.0,
	true, false, Interface::lowerlim);

	static Parameter<Emitter,double> interfacePMinusOrdering
	("PMinusOrdering",
	"An extra factor strengthening the p- ordering on top of the PSInflation."
	"A large value will suppress large dipoles more.",
	&Emitter::thePMinusOrdering, 1.0, 1.0, 0.0, 0.0,
	true, false, Interface::lowerlim);

	static Parameter<Emitter,double> interfacePTScale
	("PTScale",
	"If pT is 1/r or 2/r.",
	&Emitter::thePTScale, 1.0, 1.0, 0.0, 0.0,
	true, false, Interface::lowerlim);

	static Parameter<Emitter,InvEnergy> interfaceRMax
	("RMax",
	"The confinement scale (in iverse GeV). If set to zero, "
	"the value of <interface>DipoleEventHandler::RMax</interface> of the "
	"controlling event handler will be used.",
	&Emitter::theRMax, InvGeV, 0.0InvGeV, 0.0InvGeV, 0*InvGeV,
	true, false, Interface::lowerlim);

	static Switch<Emitter,bool> interfaceBothOrderedEvo
	("BothOrderedEvo",
	"If an emission should be fully ordered with both its parents."
	"otherwise it will be weighted according to how close they are.",
	&Emitter::theBothOrderedEvo, false, true, false);
	static SwitchOption interfaceBothOrderedEvoTrue
	(interfaceBothOrderedEvo,"True","both parents fully ordered.",true);
	static SwitchOption interfaceBothOrderedEvoFalse
	(interfaceBothOrderedEvo,"False","weighted by distance.",false);

	static Switch<Emitter,bool> interfaceBothOrderedInt
	("BothOrderedInt",
	"If an emission should be fully ordered with both its parents."
	"otherwise it will be weighted according to how close they are.",
	&Emitter::theBothOrderedInt, false, true, false);
	static SwitchOption interfaceBothOrderedIntTrue
	(interfaceBothOrderedInt,"True","both parents fully ordered.",true);
	static SwitchOption interfaceBothOrderedIntFalse
	(interfaceBothOrderedInt,"False","weighted by distance.",false);

	static Switch<Emitter,bool> interfaceBothOrderedFS
	("BothOrderedFS",
	"If an emission should be fully ordered with both its parents."
	"otherwise it will be weighted according to how close they are.",
	&Emitter::theBothOrderedFS, false, true, false);
	static SwitchOption interfaceBothOrderedFSTrue
	(interfaceBothOrderedFS,"True","both parents fully ordered.",true);
	static SwitchOption interfaceBothOrderedFSFalse
	(interfaceBothOrderedFS,"False","weighted by distance.",false);

	static Switch<Emitter,int> interfaceRangeMode
	("RangeMode",
	"How the range of coherenent emissions is determined.",
	&Emitter::theRangeMode, 0, true, false);
	static SwitchOption interfaceRangeModeMin
	(interfaceRangeMode,
	"Min",
	"minimum of half mother dipole and distance to emission. Default.",
	0);
	static SwitchOption interfaceRangeModeMax
	(interfaceRangeMode,
	"Mother",
	"Half distance of mother dipole, no matter how close the emission is.",
	1);

	static Switch<Emitter,int> interfaceMinusOrderingMode
	("MinusOrderingMode",
	"Sets how the ordering in p- is done in the virtual cascade, in relation to effective partons mainly.",
	&Emitter::theMinusOrderingMode, 0, true, false);
	static SwitchOption interfaceMinusOrderingModeEffectiveParton
	(interfaceMinusOrderingMode,
	"EffectiveParton",
	"Uses the momentum of the effective parton to order, both plus and pt.",
	0);
	static SwitchOption interfaceMinusOrderingModeEffectivePT
	(interfaceMinusOrderingMode,
	"EffectivePT",
	"Uses the pt of the effective parton, but the rapidity of the single emitting parton.",
	1);

	static SwitchOption interfaceMinusOrderingModeTrueShadow
	(interfaceMinusOrderingMode,
	"TrueShadow",
	"Use the full shadow mechanism for resolved emissions to get the "
	"incoming propagator.",
	2);

	}

	diff --git a/DIPSY/ShadowParton.cc b/DIPSY/ShadowParton.cc
	--- a/DIPSY/ShadowParton.cc
	+++ b/DIPSY/ShadowParton.cc
	@@ -1,676 +1,669 @@
	// -- C++ --
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the ShadowParton class.
	//

	#include "ShadowParton.h"
	#include "Parton.h"
	#include "DipoleState.h"
	#include "DipoleEventHandler.h"
	#include "Dipole.h"
	#include "ImpactParameters.h"
	#include "ThePEG/EventRecord/Particle.h"
	#include "ThePEG/Repository/CurrentGenerator.h"
	#include "ThePEG/Utilities/Debug.h"
	#include "ThePEG/Utilities/Throw.h"
	#include "ThePEG/Utilities/EnumIO.h"
	#include "ThePEG/Utilities/Current.h"

	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "ThePEG/Utilities/DebugItem.h"

	using namespace DIPSY;

	ShadowParton::ShadowParton(Parton & p)
	- : theOriginal(&p), thePosition(p.position()), thePlus(p.plus()),
	+ : theOriginal(&p), thePlus(p.plus()),
	thePT(p.pT()), thePT0(p.pT()), theMass(p.mass()), theMinus(p.minus()),
	theY(p.y()), theY0(p.y()), theFlavour(p.flavour()),
	- theDifference(first), theDipoles(p.dipoles()), theEmissionFactor(ZERO),
	+ theDipoles(p.dipoles()), theEmissionFactor(ZERO),
	theRes(ZERO), isLocked(false), hasInteracted(false), isOnShell(false),
	isValence(p.valence()), memememe(false) {}

	ShadowParton::~ShadowParton() {}

	void ShadowParton::setupParton() {
	if ( tSPartonPtr sp = original()->shadow() ) {
	sp->theNext = this;
	thePrevious = sp;
	}
	original()->shadow(this);
	}

	SPartonPtr ShadowParton::createValence(Parton & p) {
	SPartonPtr sp = new_ptr(ShadowParton(p));
	sp->setupParton();
	sp->isValence = true;
	return sp;
	}

	void ShadowParton::setupSwing(Parton & pa1, Parton & pa2, Parton & pb1, Parton & pb2) {
	// SPartonPtr sa1 = new_ptr(ShadowParton(pa1));
	// sa1->setupParton();
	// SPartonPtr sa2 = new_ptr(ShadowParton(pa2));
	// sa2->setupParton();
	// SPartonPtr sb1 = new_ptr(ShadowParton(pb1));
	// sb1->setupParton();
	// SPartonPtr sb2 = new_ptr(ShadowParton(pb2));
	// sb2->setupParton();

	- // sa1->theDifference = sa2->theDifference = sb1->theDifference = sb2->theDifference = swing;
	// sa1->setupNeighbors();
	// sa2->setupNeighbors();
	// sb1->setupNeighbors();
	// sb2->setupNeighbors();

	}

	void ShadowParton::setupEmission(Parton & emitter, Parton & produced, Parton & recoiler) {
	SPartonPtr se = new_ptr(ShadowParton(emitter));
	se->setupParton();
	SPartonPtr sp = new_ptr(ShadowParton(produced));
	sp->setupParton();
	// SPartonPtr sr = new_ptr(ShadowParton(recoiler));
	// sr->setupParton();

	- se->theDifference = emission;
	se->theSibling = sp;
	se->setupNeighbors();

	- sp->theDifference = emission;
	sp->theParent = this;
	sp->theSibling = se;
	theChild = sp;
	sp->setupNeighbors();

	- // sr->theDifference = recoil;
	// sr->setupNeighbors();

	InvEnergy2 d2 = se->theRes = sp->theRes = se->dist2(*sp);
	se->theEmissionFactor = sp->theEmissionFactor = alphaS(d2)*d2/UseRandom::rnd();
	- se->pT0(sp->pT0());
	- se->y0(sp->y0());
	+ // se->pT0(sp->pT0());
	+ // se->y0(sp->y0());
	}

	void ShadowParton::setupNeighbors() {
	if ( dipoles().first )
	theNeighbors.first = dipoles().first->partons().first->shadow();
	if ( dipoles().second )
	theNeighbors.second = dipoles().second->partons().second->shadow();
	}

	void ShadowParton::lock() {
	if ( mother() ) mother()->lock();
	isLocked = true;
	}

	void ShadowParton::unlock() {
	if ( mother() && locked() ) mother()->unlock();
	isLocked = false;
	}

	tSPartonPtr ShadowParton::last() {
	if ( next() ) return next();
	return this;
	}

	tSPartonPtr ShadowParton::initial() {
	if ( previous() ) return previous();
	return this;
	}

	tSPartonPtr ShadowParton::valenceMother() {
	if ( mother() ) return mother()->valenceMother();
	return this;
	}

	tcSPartonPtr ShadowParton::last() const {
	if ( next() ) return next();
	return this;
	}

	double ShadowParton::alphaS(InvEnergy2 r2) {
	return Current<DipoleEventHandler>()->alphaS(sqrt(r2));
	}

	// tcSPartonPtr ShadowParton::resolve(InvEnergy2 r2, tPartonPtr stopp) const {
	// if ( resolved(r2, stopp) \|\| !mother() ) return this;
	// return mother()->resolve(r2, stopp);
	// }

	tSPartonPtr ShadowParton::resolve(InvEnergy2 r2) {
	if ( resolved(r2) \|\| !mother() ) return this;
	return mother()->resolve(r2);
	}

	tSPartonPtr ShadowParton::resolveInteraction(InvEnergy2 r2) {
	if ( resolved(r2)
	\|\| !mother()
	\|\| forceEmission()
	\|\| onShell()
	\|\| ( sibling() && sibling()->interacted() ) ) return this;
	return mother()->resolveInteraction(r2);
	}

	// void ShadowParton::flagOnShell(InvEnergy2 r2, tPartonPtr stopp) {
	// tSPartonPtr resolved = resolve(r2, stopp);
	// if ( resolved->onShell() ) return;
	// resolved->onShell(true);
	// if ( resolved->sibling() ) {
	// r2 = resolved->res();
	// resolved->sibling()->onShell(true);
	// }
	// if ( resolved->mother() ) resolved->mother()->flagOnShell(r2, stopp);
	// // * TODO * probably the whole function should be deleted
	// // else dipoleState().flagShadowsOnShell();
	// }

	void ShadowParton::pTplus(const TransverseMomentum & qt, Energy qp) {
	pT(qt);
	plus(qp);
	minus(mt2()/plus());
	y(log(mt()/plus()));
	}

	void ShadowParton::setOnShell(int mode) {
	if ( mode < 0 ) return;
	onShell(true);
	if ( mode <= 0 ) return;
	original()->plus(plus());
	original()->pT(pT());
	original()->minus(original()->mt2()/plus());
	original()->y(log(original()->mt()/original()->plus()));
	original()->onShell(true);
	}

	void ShadowParton::unsetOnShell() {
	if ( onShell() ) {
	original()->onShell(false);
	onShell(false);
	}
	}

	void ShadowParton::interact() {
	if ( interacted() ) return;
	hasInteracted = true;
	if ( mother() ) mother()->interact();
	}

	void ShadowParton::resetInteracted() {
	if ( !mother() ) {
	pT(pT0());
	y(y0());
	plus(mt()*exp(-y()));
	minus(mt()*exp(y()));
	}
	hasInteracted = false;
	theInteractions.clear();
	onShell(false);
	if ( next() ) next()->resetInteracted();
	if ( child() ) child()->resetInteracted();
	}

	void ShadowParton::reset() {
	if ( !mother() ) {
	pT(pT0());
	y(y0());
	plus(mt()*exp(-y()));
	minus(mt()*exp(y()));
	}
	onShell(false);
	theInteractions.clear();
	hasInteracted = false;
	theInteractionRoots.clear();
	if ( next() ) next()->reset();
	if ( child() ) child()->reset();
	}

	void ShadowParton::insertInteraction(int i) {
	// If this parton has been set on-shell by a previous interaction we
	// do not have to look further. Similarly if the sibling has
	// interacted.
	if ( onShell()
	\|\| ( sibling() && sibling()->interacted() )
	\|\| !mother() ) {
	interactionRoot(i);
	}
	// Otherwise we try with the mother instead.
	else mother()->insertInteraction(i);
	}

	void ShadowParton::rejectInteraction(int i) {
	if ( hasInteractionRoot(i) ) theInteractionRoots.erase(i);
	else mother()->rejectInteraction(i);
	}

	void ShadowParton::acceptInteraction(int i) {
	hasInteracted = true;
	if ( !hasInteractionRoot(i) ) mother()->acceptInteraction(i);
	}

	void ShadowParton::makeIncoming() {
	if ( interacted() && onShell() ) onShell(false);
	else if ( mother() ) mother()->makeIncoming();
	}

	bool ShadowParton::setEmissionMomenta(const LorentzMomentum & p,
	bool forced) {
	// * TODO * Think this through! If the emission unresolved but
	// anyway should be done because a valence need to be put on shell,
	// or a subsequent emission needs it, the generated pt is ignored
	// and the emitted parton and the continuing propagator will share
	// the incoming pt, while the generated rapidity of the emission
	// remains untouched.
	- TransverseMomentum qT = pT0();
	+ TransverseMomentum qT = parent()? pT0(): sibling()->pT0();
	if ( forced ) qT = TransverseMomentum(-p.x()/2.0, -p.y()/2.0);

	if ( parent() ) {
	// If this parton was emitted, pretend it was actually he emitter
	// and continues as a propagator. It will retain its original
	// rapidity, but will get the transverse momentum of the incoming
	// propagator.
	pT(TransverseMomentum(p.x(), p.y()) + qT);
	plus(mt0()*exp(y0()));
	y(y0());
	// The sibling will get minus the original transverse momentum and
	// what ever is left of the positive light-cone momenta. (return
	// false if not enough to put on shell).
	sibling()->pT(-qT);
	sibling()->plus(p.plus() - mt()*exp(-y()));
	if ( sibling()->plus() <= ZERO ) return false;
	sibling()->minus(mt2()/sibling()->plus());
	sibling()->y(log(sibling()->mt()/sibling()->plus()));
	// The propagators negative light-cone momentum is trivial,
	minus(p.minus() - sibling()->minus());
	} else {
	// If this parton was assumed in the evolution, life is simpler.
	sibling()->pT(-qT);
	- sibling()->plus(sibling()->mt()*exp(-y0()));
	+ sibling()->plus(sibling()->mt()*exp(-sibling()->y0()));
	sibling()->minus(sibling()->mt2()/sibling()->plus());
	- sibling()->y(y0());
	+ sibling()->y(sibling()->y0());
	pT(TransverseMomentum(p.x(), p.y()) + qT);
	plus(p.plus() - sibling()->plus());
	if ( plus() <= ZERO ) return false;
	minus(mt2()/plus());
	y(log(mt()/plus()));
	}
	return true;
	}

	// ShadowParton::Propagator ShadowParton::
	// oldpropagator(InvEnergy2 r2, tPartonPtr stopp, int mode) {
	// // If we reach a parton that has been set on shell, we will simply
	// // return its momentum and set it off-shell.
	// if ( onShell() ) {
	// if ( mode >= 0 ) unsetOnShell();
	// return momentum();
	// }
	// // If this was a valence, ask the DipoleState about its momentum.
	// if ( !mother() ) return dipoleState().incomingMomentum(this, mode);

	// // Now, if this emission was not resolved then simply return the
	// // propagator of the mother. However, if the emission is needed for
	// // subsequent interactions or if the emitted parton was a valence,
	// // we need to force it. with a special procedure.
	// bool forced = forceEmission();
	// bool unresolved = !resolved(r2, stopp);

	// if ( unresolved && !forced )
	// return setMomentum(mother()->propagator(r2, stopp, mode));

	// // Get the propagator of the mother with the new resolution
	// // scale. But don't put anything on-shell yet in case the emission
	// // is rejected later.
	// Propagator prop = mother()->propagator(res(), stopp, -1);

	// // If it is not possible to set the emitted parton on shell, flag
	// // the propagator as failed.
	// if ( !setEmissionMomenta(prop.p, false) ) {
	// prop.fail = true;
	// // If the emitted parton is necessary for a subsequent
	// // interaction, return the failed propagator.
	// if ( forced ) return prop;
	// // Otherwise just ignore the emission.
	// return setMomentum(mother()->propagator(r2, stopp, mode));

	// }

	// // Now depending on the colour line of the emitted parton, it must
	// // be ordered in light-cone momenta with previous partons on the
	// // same side except if it partakes in a subsequent interaction.
	// if ( colourSibling() ) {
	// if ( sibling()->interactionRoot() ) {
	// prop.colpos = Constants::MaxEnergy;
	// prop.colneg = ZERO;
	// } else {
	// // Check if a resolved emission was not ordered, in that case
	// // flag it unresolved.
	// if ( sibling()->plus() > prop.colpos \|\|
	// sibling()->minus() < prop.colneg )
	// unresolved = true;
	// else {
	// prop.colpos = sibling()->plus();
	// prop.colneg = sibling()->minus();
	// }
	// }
	// } else {
	// if ( sibling()->interactionRoot() ) {
	// prop.acopos = Constants::MaxEnergy;
	// prop.aconeg = ZERO;
	// } else {
	// if ( sibling()->plus() > prop.acopos \|\|
	// sibling()->minus() < prop.aconeg )
	// unresolved = true;
	// else {
	// prop.acopos = sibling()->plus();
	// prop.aconeg = sibling()->minus();
	// }
	// }
	// }

	// // If the emission was found to be unresolved we need to recalculate
	// // the incoming momentum.
	// if ( unresolved ) {
	// prop = mother()->propagator(r2, stopp, mode);
	// // If the unresolved emission was a valence we need to emit it
	// // anyway with a special treatment.
	// if ( forced ) {
	// if ( !setEmissionMomenta(prop.p, true) ) {
	// prop.fail = true;
	// return prop;
	// }
	// } else
	// return setMomentum(prop);
	// }
	// // In any case we need to redo the incoming momentum to mark it final.
	// if ( mode >= 0 ) prop = mother()->propagator(res(), stopp, mode);

	// // Now we are done. Depending on the mode we will put the emitted
	// // parton on shell.
	// sibling()->setOnShell(mode);
	// prop.p -= sibling()->momentum();

	// return prop;
	// }

	ShadowParton::Propagator ShadowParton::
	propagator(InvEnergy2 r2, int mode) {

	static DebugItem unorderlock("DIPSY::UnorderLock", 6);

	Propagator prop;

	// If we reach a parton that has been set on shell, we will simply
	// return its momentum and set it off-shell.
	if ( onShell() ) {
	if ( mode >= 0 ) unsetOnShell();
	return momentum();
	}
	// If this was a valence, ask the DipoleState about its momentum.
	if ( !mother() ) return dipoleState().incomingMomentum(this, mode);

	// Now, if this emission was not resolved then simply return the
	// propagator of the mother. However, if the emission is needed for
	// subsequent interactions or if the emitted parton was a valence,
	// we need to force it. with a special procedure.
	bool forced = forceEmission();
	bool unresolved = !resolved(r2);


	// First we check if a resolved emission is possible.
	if ( !unresolved ) {

	// Get the propagator of the mother with the new resolution
	// scale. But don't put anything on-shell yet in case the emission
	// is rejected later.
	prop = mother()->propagator(res(), -1);

	// First check that it was at all kinematically possible perform
	// the emission.
	if ( setEmissionMomenta(prop.p, false) ) {

	// OK. That seemed to work, but we must also check the ordering.
	// Now depending on the colour line of the emitted parton, it
	// must be ordered in light-cone momenta with previous partons
	// on the same side. If that is not the case, te emissionis
	// marked unresolved.
	if ( ! ( unorderlock && sibling() && sibling()->locked() ) ) {
	if ( colourSibling() ) {
	if ( sibling()->plus() > prop.colpos \|\|
	sibling()->minus() < prop.colneg )
	unresolved = true;
	} else{
	if ( sibling()->plus() > prop.acopos \|\|
	sibling()->minus() < prop.aconeg )
	unresolved = true;
	}
	}
	} else {
	// Since the emission could not be performed, we flaf it
	// unresolved.
	unresolved = true;
	}
	}

	// If the emission was really resolved, we get the incoming
	// propagator again, but this time put stuff on-shell. Then we
	// return the outgoing propagator, after setting the ordering for
	// the subsequent emissions, and we're done.
	if ( !unresolved ) {
	if ( mode >= 0 ) prop = mother()->propagator(res(), mode);
	if ( colourSibling() ) {
	prop.colpos = sibling()->plus();
	prop.colneg = sibling()->minus();
	} else {
	prop.acopos = sibling()->plus();
	prop.aconeg = sibling()->minus();
	}
	sibling()->setOnShell(mode);
	prop.p -= sibling()->momentum();
	return prop;
	}

	// OK, so it was unresolved. Get the incoming propagator with the
	// previous resolution scale. If it was not forced, then let it go.
	prop = mother()->propagator(r2, mode);
	if ( !forced ) return prop;

	// However, if we have to force it we need to check that it really
	if ( !setEmissionMomenta(prop.p, true) ) {
	// If it doesn't work, just give up the whole thing and fail the
	// whole propagator chain.
	prop.fail = true;
	return prop;
	}

	// If it does work we set things on-shell, return the propagator and
	// we're done.
	sibling()->setOnShell(mode);
	prop.p -= sibling()->momentum();

	return prop;
	}

	tSPartonPtr ShadowParton::findFirstOnShell() {
	if ( original()->onShell() ) return this;
	if ( !next() ) return tSPartonPtr();
	if ( next()->colourSibling() ) {
	tSPartonPtr ch = child()->findFirstOnShell();
	if ( ch ) return ch;
	}
	return next()->findFirstOnShell();
	}

	tSPartonPtr ShadowParton::findSecondOnShell() {
	if ( original()->onShell() ) return this;
	if ( !next() ) return tSPartonPtr();
	if ( !next()->colourSibling() ) {
	tSPartonPtr ch = child()->findSecondOnShell();
	if ( ch ) return ch;
	}
	return next()->findSecondOnShell();
	}

	Ariadne5::ClonePtr ShadowParton::clone() const {
	return new_ptr(*this);
	}

	-void ShadowParton::mirror(double y0) {
	- y(2.0*y0 - y());
	+void ShadowParton::mirror(double yf) {
	+ y(2.0*yf - y());
	swap(thePlus, theMinus);
	- if ( next() ) next()->mirror(y0);
	- if ( child() ) child()->mirror(y0);
	+ if ( next() ) next()->mirror(yf);
	+ if ( child() ) child()->mirror(yf);
	}

	void ShadowParton::translate(const ImpactParameters & b) {
	- thePosition = b.translate(thePosition);
	thePT = b.rotatePT(thePT);
	if ( next() ) next()->translate(b);
	if ( child() ) child()->translate(b);
	}

	void ShadowParton::propagateShadowMomenta() {
	if ( onShell() ) {
	original()->plus(plus());
	original()->pT(pT());
	original()->y(log(original()->mt()/original()->plus()));
	original()->onShell(true);
	}
	if ( next() ) next()->propagateShadowMomenta();
	if ( child() ) child()->propagateShadowMomenta();
	}

	void ShadowParton::rebind(const TranslationMap & trans) {
	thePrevious = trans.translate(thePrevious);
	theParent = trans.translate(theParent);
	theSibling = trans.translate(theSibling);
	theNext = trans.translate(theNext);
	theChild = trans.translate(theChild);
	theDipoles.first = trans.translate(theDipoles.first);
	theDipoles.second = trans.translate(theDipoles.second);
	theNeighbors.first = trans.translate(theNeighbors.first);
	theNeighbors.second = trans.translate(theNeighbors.second);
	}


	DipoleState& ShadowParton::dipoleState() const {
	if ( dipoles().first ) return dipoles().first->dipoleState();
	else return dipoles().second->dipoleState();
	}


	double ShadowParton::pTScale() const {
	- return dipoleState().handler().emitter().pTScale();
	+ return Current<DipoleEventHandler>()->emitter().pTScale();
	}


	Energy ShadowParton::mass() const {
	if ( theMass < ZERO )
	theMass = CurrentGenerator::current().getParticleData(flavour())->mass();
	return theMass;
	}

	bool ShadowParton::partonOnShell() const {
	return onShell() \|\| ( previous() && previous()->partonOnShell() );
	}

	void ShadowParton::debugme() {
	// cout << "data for ShadowParton " << this << endl;
	// cout << "thePosition1: " << thePosition.first*GeV
	// << ", thePosition2: " << thePosition.second*GeV
	// << ", thePlus: " << thePlus/GeV
	// << ", thePT: " << thePT.pt()/GeV
	// << ", theMinus: " << theMinus/GeV
	// << ", theY: " << theY
	// << ", theFlavour: " << theFlavour
	// << ", theDipoles1: " << theDipoles.first
	// << ", theDipoles2: " << theDipoles.second
	// << ", hasInteracted: " << hasInteracted
	// << ", isOnShell: " << isOnShell
	// << ", isValence: " << isValence
	// << ", theMass: " << theMass/GeV << endl;
	memememe = true;
	valenceMother()->debugTree("");
	memememe = false;

	}

	void ShadowParton::debugTree(string indent) {
	- // if ( mother() && diff() != emission && next() ) {
	- // next()->debugTree(indent);
	- // return;
	- // }
	cerr << indent << (memememe? "==": "--")
	<< (valence()? "v": "-") << (locked()? "%": "-");
	if ( interactionRoot() ) {
	set<int>::iterator i = theInteractionRoots.begin();
	cerr << "{" << *i;
	while ( ++i != theInteractionRoots.end() )
	cerr << "," << *i;
	cerr << "}";
	}
	cerr << (interacted()? "+": "-");
	if ( interacting() ) {
	set<int>::iterator i = theInteractions.begin();
	cerr << "[" << *i;
	while ( ++i != theInteractions.end() )
	cerr << "," << *i;
	cerr << "]";
	}
	cerr << (onShell()? "*": "-") << (memememe? "=>": "->");
	if ( onShell() )
	cerr << "[" << pT().x()/GeV << ", "
	<< pT().y()/GeV << ", "
	<< plus()/GeV << "]";
	else if ( interacted() )
	cerr << "(" << pT().x()/GeV << ", "
	<< pT().y()/GeV << ", "
	<< plus()/GeV << ")";
	cerr << endl;
	if ( indent.size() > 1 && indent[indent.size() - 1] == '\\' )
	indent[indent.size() - 1] = ' ';
	if ( next() && child() ) {
	if ( child()->interacted() ) next()->debugTree(indent + " +");
	else next()->debugTree(indent + " \|");
	}
	if ( next() && !child() ) next()->debugTree(indent + " ");
	if ( child() ) child()->debugTree(indent + " \\");
	}
	/*

	----->
	\|---->
	\| \|---->
	\| \---->
	\| \|---->
	\| \---->
	\---->
	*/


	void ShadowParton::persistentOutput(PersistentOStream & os) const {
	- os << theOriginal << ounit(thePosition, InvGeV) << ounit(thePlus, GeV)
	- << ounit(thePT, GeV) << ounit(theMass, GeV) << ounit(theMinus, GeV)
	- << theY << theFlavour << thePrevious << theParent << theSibling
	+ os << theOriginal << ounit(thePlus, GeV)
	+ << ounit(thePT, GeV) << ounit(thePT0, GeV) << ounit(theMass, GeV)
	+ << ounit(theMinus, GeV) << theY << theY0 << theFlavour
	+ << thePrevious << theParent << theSibling
	<< theNext << theChild << theDipoles << theNeighbors
	<< ounit(theEmissionFactor, 1.0/GeV2) << ounit(theRes, 1.0/GeV2)
	- << hasInteracted << isOnShell << isValence << oenum(theDifference);
	+ << hasInteracted << isOnShell << isValence;
	}

	void ShadowParton::persistentInput(PersistentIStream & is, int) {
	- is >> theOriginal >> iunit(thePosition, InvGeV) >> iunit(thePlus, GeV)
	- >> iunit(thePT, GeV) >> iunit(theMass, GeV)>> iunit(theMinus, GeV)
	- >> theY >> theFlavour >> thePrevious >> theParent >> theSibling
	+ is >> theOriginal >> iunit(thePlus, GeV)
	+ >> iunit(thePT, GeV) >> iunit(thePT0, GeV) >> iunit(theMass, GeV)
	+ >> iunit(theMinus, GeV) >> theY >> theY0 >> theFlavour
	+ >> thePrevious >> theParent >> theSibling
	>> theNext >> theChild >> theDipoles >> theNeighbors
	>> iunit(theEmissionFactor, 1.0/GeV2) >> iunit(theRes, 1.0/GeV2)
	- >> hasInteracted >> isOnShell >> isValence >> ienum(theDifference);
	+ >> hasInteracted >> isOnShell >> isValence;
	}

	DescribeClass<ShadowParton,Ariadne5::CloneBase>
	describeDIPSYShadowParton("DIPSY::ShadowParton", "libAriadne5.so libDIPSY.so");
	// Definition of the static class description member.

	void ShadowParton::Init() {}

	diff --git a/DIPSY/ShadowParton.h b/DIPSY/ShadowParton.h
	--- a/DIPSY/ShadowParton.h
	+++ b/DIPSY/ShadowParton.h
	@@ -1,1038 +1,1003 @@
	// -- C++ --
	#ifndef DIPSY_ShadowParton_H
	#define DIPSY_ShadowParton_H
	//
	// This is the declaration of the ShadowParton class.
	//

	#include "ThePEG/Config/ThePEG.h"
	#include "ThePEG/Vectors/Transverse.h"
	#include "ThePEG/PDT/EnumParticles.h"
	#include "ThePEG/Repository/UseRandom.h"
	#include "Ariadne/Config/CloneBase.h"
	#include "ShadowParton.fh"
	#include "Parton.h"
	#include "Dipole.fh"
	#include "DipoleXSec.fh"

	namespace DIPSY {

	using namespace ThePEG;

	class ImpactParameters;

	/**
	* Here is the documentation of the ShadowParton class.
	*/
	class ShadowParton: public Ariadne5::CloneBase {

	public:

	/**
	* Typedef for position in transverse coordinate space.
	*/
	typedef Transverse<InvEnergy> Point;

	/**
	* A pair of partons.
	*/
	typedef pair<tSPartonPtr,tSPartonPtr> tSPartonPair;

	/**
	* A pair of dipoles.
	*/
	typedef pair<tDipolePtr,tDipolePtr> tDipolePair;

	/**
	* The Dipole is a friend.
	*/
	friend class Dipole;

	/**
	- * Enum difference from the previous version of the same parton.
	- */
	- enum Difference {
	- first = 0, /*< No difference - this is the first instance. /
	- emission = 1, /*< The previous instance emitted. /
	- recoil = 2, /*< The previous received a recoil from a neighboring emission /
	- swing = 3, /*< The previous underwent a swing. /
	- interactn = 4 /*< The previous underwent an interaction /
	- };
	-
	- /**
	* Helper class for traversing along a propagator to an interaction
	*/
	struct Propagator {

	/**
	* Constructor.
	*/
	Propagator(): colpos(Constants::MaxEnergy), colneg(ZERO),
	acopos(Constants::MaxEnergy), aconeg(ZERO), fail(false) {}

	/**
	* Construct from momentum.
	*/
	Propagator(const LorentzMomentum & mom)
	: p(mom), colpos(Constants::MaxEnergy), colneg(ZERO),
	acopos(Constants::MaxEnergy), aconeg(ZERO), fail(false) {}

	/**
	* The Momentum of the propagator.
	*/
	LorentzMomentum p;

	/**
	* The positive light-cone momentum of the previous emission on the colour line.
	*/
	Energy colpos;

	/**
	* Thenegative light-cone momentum of the previous emission on the colour line.
	*/
	Energy colneg;

	/**
	* The positive light-cone momentum of the previous emission on
	* the anti-colour line.
	*/
	Energy acopos;

	/**
	* Thenegative light-cone momentum of the previous emission on the
	* anti-colour line.
	*/
	Energy aconeg;

	/**
	* Indicate that the propagator has failed.
	*/
	bool fail;

	};

	/**
	* Helper class to temporarily protect a propagator.
	*/
	struct Lock {

	/**
	* The only constructor. Locking the propagator to a given parton.
	*/
	Lock(tPartonPtr p): sp(p->shadow()) {
	sp->lock();
	}

	/**
	* The destructor, unlocking the propagator.
	*/
	~Lock() {
	sp->unlock();
	}

	/**
	* The parton to which the propagator should be locked.
	*/
	tSPartonPtr sp;

	};


	public:

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* The default constructor.
	*/
	inline ShadowParton()
	- : thePosition(Point()), thePlus(ZERO), theMass(-1GeV), theMinus(0GeV),
	- theY(0), theY0(0), theFlavour(ParticleID::g), theDifference(first),
	+ : thePlus(ZERO), theMass(-1GeV), theMinus(0GeV),
	+ theY(0), theY0(0), theFlavour(ParticleID::g),
	theEmissionFactor(ZERO), theRes(ZERO), isLocked(false),
	hasInteracted(false), isOnShell(false), isValence(false), memememe(false) {}


	/**
	* Construct from ordinary parton.
	*/
	ShadowParton(Parton & p);

	/**
	* The destructor.
	*/
	virtual ~ShadowParton();

	/**
	* Create a valence shadow parton from a given ordinary parton.
	*/
	static SPartonPtr createValence(Parton & p);

	/**
	* Easy access to alphaS.
	*/
	static double alphaS(InvEnergy2 r2);

	/**
	* Setup shadow partons in an emission.
	*/
	void setupEmission(Parton & emitter, Parton & produced, Parton & recoiler);

	/**
	* Setup shadow partons in a swing.
	*/
	static void setupSwing(Parton & pa1, Parton & pa2, Parton & pb1, Parton & pb2);

	/**
	* Setup pointers to and from neighboring partons.
	*/
	void setupParton();

	/**
	* Setup neighboring shadow partons.
	*/
	void setupNeighbors();

	//@}

	protected:

	/** @name The virtual functions to be overridden in sub-classes. */
	//@{
	/**
	* Return a simple clone of this object. Should be implemented as
	* <code>return new_ptr(*this);</code> by a derived class.
	*/
	virtual Ariadne5::ClonePtr clone() const;

	/**
	* Rebind pointers to other CloneBase objects. Called after a number
	* of interconnected CloneBase objects have been cloned, so that
	* the cloned objects will refer to the cloned copies afterwards.
	*
	* @param trans a TranslationMap relating the original objects to
	* their respective clones.
	*/
	virtual void rebind(const TranslationMap & trans);
	//@}

	public:

	/**
	* Return the DipoleState to which this parton belongs. If it was
	* originally belonged to one state which was subsequently merged
	* into another, the parton will belong to the latter state.
	*/
	DipoleState & dipoleState() const;

	/**
	* finds the pTscale through the dipolestates emitter.
	*/
	double pTScale() const;

	/**
	* Calculate the squared transverse distance to the given parton.
	*/
	inline InvEnergy2 dist2(const ShadowParton & p) const {
	return (position() - p.position()).pt2();
	}

	/**
	* Calculate the transverse distance to the given parton.
	*/
	inline InvEnergy dist(const ShadowParton & p) const {
	return sqrt(dist2(p));
	}

	/**
	* Produce a ThePEG::Particle corresponding to this parton.
	*/
	PPtr produceParticle() const;

	/**
	* The mass of this parton.
	*/
	Energy mass() const;

	/**
	* The transverse momentum squared of this parton.
	*/
	Energy2 pt2() const {
	return pT().pt2();
	}

	/**
	* The transverse mass squared of this parton.
	*/
	Energy2 mt2() const {
	return pt2() + sqr(mass());
	}

	/**
	* The transverse mass of this parton in the emission.
	*/
	Energy2 mt02() const {
	return pT0().pt2() + sqr(mass());
	}

	/**
	* The transverse mass of this parton in the emission.
	*/
	Energy mt0() const {
	return sqrt(mt02());
	}

	/**
	* The transverse momentum of this parton.
	*/
	Energy pt() const {
	return sqrt(pt2());
	}

	/**
	* The transverse momentum of this parton.
	*/
	Energy mt() const {
	return sqrt(mt2());
	}

	/**
	* The final-state momentum of this particle.
	*/
	LorentzMomentum momentum() const {
	return lightCone(plus(), mt2()/plus(), pT());
	}

	/** @name Simple access functions. */
	//@{
	/**
	* Get the position in impact parameter space.
	*/
	inline const Point & position() const {
	- return thePosition;
	+ return original()->position();
	}

	/**
	* Get the positive light-cone momentum.
	*/
	inline Energy plus() const {
	return thePlus;
	}

	/**
	* Get the original positive light-cone momentum.
	*/
	inline Energy plus0() const {
	return pT0().pt()*exp(-y0());
	}

	/**
	* Get the transverse momentum.
	*/
	inline TransverseMomentum pT() const {
	return thePT;
	}

	/**
	* Get the transverse momentum.
	*/
	inline TransverseMomentum pT0() const {
	return thePT0;
	}

	/**
	* Get the accumulated negative light-cone momentum deficit.
	*/
	inline Energy minus() const {
	return theMinus;
	}

	/**
	* Get the rapidity.
	*/
	inline double y() const {
	return theY;
	}

	/**
	* Get the rapidity generated in the emission.
	*/
	inline double y0() const {
	return theY0;
	}

	/**
	* Get the flavour of this parton.
	*/
	inline long flavour() const {
	return theFlavour;
	}

	/**
	- * Indicate what happened to the previous instance of this parton.
	- */
	- Difference diff() const {
	- return theDifference;
	- }
	-
	- /**
	* Get the original DIPSY parton.
	*/
	tPartonPtr original() const {
	return theOriginal;
	}

	/**
	* Get the previous instance of this parton befor it was emitted or
	* swinged.
	*/
	inline tSPartonPtr previous() const {
	return thePrevious;
	}

	/**
	* Get the parent parton if this was emitted.
	*/
	inline tSPartonPtr parent() const {
	return theParent;
	}

	/**
	* Get the sibling of this parton if any.
	*/
	inline tSPartonPtr sibling() const {
	return theSibling;
	}

	/**
	* Get the next version of this parton if emitted or swinged.
	*/
	inline tSPartonPtr next() const {
	return theNext;
	}

	/**
	* Get the first version of this parton if emitted or swinged.
	*/
	tSPartonPtr initial();

	/**
	* Get the first (grand)mother of this parton.
	*/
	tSPartonPtr valenceMother();

	/**
	* Get the last version of this parton if emitted or swinged.
	*/
	tSPartonPtr last();

	/**
	* Get the last version of this parton if emitted or swinged.
	*/
	tcSPartonPtr last() const;

	/**
	* Get the parton this has emitted if any.
	*/
	inline tSPartonPtr child() const {
	return theChild;
	}

	/**
	* Get the mother parton, either previous() or parent().
	*/
	inline tSPartonPtr mother() const {
	return parent()? parent(): previous();
	}

	/**
	* Get the connecting dipoles.
	*/
	inline tDipolePair dipoles() const {
	return theDipoles;
	}

	/**
	* Indicate that the corresponding emission is on a protected
	* propagator. Any other propagator must assume that this emission
	* must be resolved.
	*/
	inline bool locked() const {
	return isLocked;
	}

	/**
	* Indicate that the propagator leading up to this parton is protected.
	*/
	void lock();

	/**
	* Indicate that the propagator leading up to thi parton is no
	* longer protected.
	*/
	void unlock();

	/**
	* Indicate if this parton has interacted.
	*/
	inline bool interacted() const {
	return hasInteracted;
	}

	/**
	* Return true if the sibling is on the colour side of this parton.
	*/
	inline bool colourSibling() const {
	return sibling() == theNeighbors.first;
	}

	/**
	* Return if this shadow was the root of an interacting propagator.
	*/
	inline bool interactionRoot() const {
	return theInteractionRoots.size() > 0;
	}

	/**
	* Return if this shadow was the root of an interacting propagator.
	*/
	inline bool hasInteractionRoot(int i) const {
	return theInteractionRoots.find(i) != theInteractionRoots.end();
	}

	/**
	* Return if this shadow was directly involved in the interaction.
	*/
	inline int interacting() const {
	return theInteractions.size() > 0;
	}

	/**
	* Return if this parton is on shell.
	*/
	inline bool onShell() const {
	return isOnShell;
	}

	/**
	* Return true if the actual parton should be put on-shell;
	*/
	bool partonOnShell() const;

	/**
	* Set if this parton is on shell.
	*/
	inline void onShell(bool b) {
	isOnShell = b;
	}

	/**
	* Indicate that this parton is the root of an interaction.
	*/
	inline void interactionRoot(int i) {
	theInteractionRoots.insert(i);
	}

	/**
	* Set if this shadow is directly involved in an interaction.
	*/
	inline void interacting(int i) {
	theInteractions.insert(i);
	}

	/**
	* Return if this parton is a valence parton.
	*/
	inline bool valence() const {
	return isValence;
	}

	/**
	* Set if this parton is a valence parton.
	*/
	inline void valence(bool b) {
	isValence = b;
	}

	/**
	- * Set the position in impact parameter space.
	- */
	- inline void position(const Point & x) {
	- thePosition = x;
	- }
	-
	- /**
	* Set the positive light-cone momentum.
	*/
	inline void plus(Energy x) {
	thePlus = x;
	}

	/**
	* Set the transverse momentum.
	*/
	inline void pT(const TransverseMomentum & x) {
	thePT = x;
	}

	/**
	* Set the momentum using transverse momentum and positive
	* light-cone momentum. Set also other components assuming real.
	*/
	void pTplus(const TransverseMomentum & qt, Energy qp);

	/**
	* Set the transverse momentum.
	*/
	inline void pT0(const TransverseMomentum & x) {
	thePT0 = x;
	}

	/**
	* Set the accumulated negative light-cone momentum deficit.
	*/
	inline void minus(Energy x) {
	theMinus = x;
	}

	/**
	* Set the rapidity.
	*/
	inline void y(double x) {
	theY = x;
	}

	/**
	* Set the rapidity.
	*/
	inline void y0(double x) {
	theY0 = x;
	}

	/**
	* Set the flavour of this parton.
	*/
	inline void flavour(long x) {
	theFlavour = (x == 0? ParticleID::g: x);
	}

	/**
	* Set the parent parton.
	*/
	inline void parent(tSPartonPtr p) {
	theParent = p;
	}

	/**
	* Set the connecting dipoles.
	*/
	inline void dipoles(tDipolePair x) {
	theDipoles = x;
	}

	/**
	* Indicate that this parton has interacted and mark all parent and
	* original partons if needed.
	*/
	void interact();

	/**
	* The emission factor for this parton if emitted or emitter.
	*/
	InvEnergy2 emissionFactor() const {
	return theEmissionFactor;
	}

	/**
	* The resolusion scale factor for this parton if emitted or emitter.
	*/
	InvEnergy2 res() const {
	return theRes;
	}

	/**
	* Return true if this parton and its sibling are resolved.
	*/
	bool resolved(InvEnergy2 r2) const {
	return ( ( sibling() && ( sibling()->locked()
	\|\| emissionFactor() > alphaS(r2)*r2 ) ) );
	}

	/**
	* Return true if this parton and its sibling are resolved. An
	* emission involving \a stopp is always resolved.
	*/
	// bool resolved(InvEnergy2 r2, tPartonPtr stopp) const {
	// return ( ( sibling() && ( sibling()->original() == stopp
	// \|\| sibling()->locked()
	// \|\| emissionFactor() > alphaS(r2)*r2 ) ) );
	// }

	/**
	* Return true if the sibling must be put on-shell even if it was not
	* resolved, because it will partake in a subsequent interaction or is
	* a valence.
	*/
	bool forceEmission() const {
	return sibling() && ( sibling()->valence()
	\|\| sibling()->interactionRoot() );
	}


	/**
	* Go back in the history of this shadow parton and find the
	* original parton which would be resolved at a given size or would
	* be the interaction root. Always resolve emissions involving \a
	* stopp.
	*/
	tSPartonPtr resolveInteraction(InvEnergy2 r2, tPartonPtr stopp) {
	Lock safeguard(stopp);
	return resolveInteraction(r2);
	}

	/**
	* Go back in the history of this shadow parton and find the
	* original parton which would be resolved at a given size or would
	* be the interaction root.
	*/
	tSPartonPtr resolveInteraction(InvEnergy2 r2);

	/**
	* Go back in the history of this shadow parton and find the
	* original parton which would be resolved at a given size. Always
	* resolve emissions involving \a stopp.
	*/
	tSPartonPtr resolve(InvEnergy2 r2, tPartonPtr stopp) {
	Lock safeguard(stopp);
	return resolve(r2);
	}

	/**
	* Go back in the history of this shadow parton and find the
	* original parton which would be resolved at a given size. Always
	* resolve emissions involving \a stopp. (const version).
	*/
	// tcSPartonPtr resolve(InvEnergy2 r2, tPartonPtr stopp) const;

	/**
	* Go back in the history and flag all partons which are resoved by
	* the given size to be set on-shell. Always
	* resolve emissions involving \a stopp.
	*/
	// void flagOnShell(InvEnergy2 r2, tPartonPtr stopp);

	/**
	* If this parton is not on shell, find an emitted parton
	* with the same colour lines (modulo swing) that is.
	*/
	tSPartonPtr findFirstOnShell();

	/**
	* If this parton is not on shell, find an emitted parton
	* with the same colour lines (modulo swing) that is.
	*/
	tSPartonPtr findSecondOnShell();

	/**
	* Flag on-shell if \a mode >= 0 and copy to original parton if \a
	* mode > 0.
	*/
	void setOnShell(int mode);

	/**
	* If prevoíously set on-shell, clear the flag.
	*/
	void unsetOnShell();

	/**
	* Go forward in the evolution and reset all interaction flags.
	*/
	void reset();

	/**
	* Go forward in the evolution and reset all interacted flags.
	*/
	void resetInteracted();

	/**
	* Set the momenta in an emission, assuming the incoming momentum \a p.
	*/
	bool setEmissionMomenta(const LorentzMomentum & p, bool valencefix);

	protected:

	/**
	* Go back in the history of this shadow parton and find the
	* original parton which would be resolved at a given size. Always
	* resolve emissions involving \a stopp.
	*/
	tSPartonPtr resolve(InvEnergy2 r2);

	public:

	/**
	* Go back in the history and find the momentum of this incoming
	* parton. Always resolve emissions involving \a stopp. If \a mode
	* >= 0 flag partons which would be put on shell. if \a mode > 0
	* also propagate the momentum to the original() parton.
	*/
	Propagator oldpropagator(InvEnergy2 r2, tPartonPtr stopp, int mode);
	Propagator propagator(InvEnergy2 r2, tPartonPtr stopp, int mode) {
	Lock safeguard(stopp);
	return propagator(r2, mode);
	}
	Propagator propagator(InvEnergy2 r2, int mode);

	Propagator setMomentum(const Propagator & prop) {
	plus(prop.p.plus());
	pT(TransverseMomentum(prop.p.x(), prop.p.y()));
	return prop;
	}

	/**
	* Recursively tell all this and any offspring to propagate their
	* momentum to the original partons if set on shell.
	*/
	void propagateShadowMomenta();

	/**
	* Indicate that this parton should be tested for an
	* interaction. Insert the corresponding interaction object at the
	* correct place in the chain. The primary interaction is always at
	* one of the incoming shadows. Secondary interactions are inserted
	* at a parton which has been or should be made real by a previous
	* emission. If this parton already had an interaction, simply leave
	* it as it is.
	*/
	void insertInteraction(int i);

	/**
	* Indicate that a previously tested interaction is accepted. This
	* is done by marking all propagators interacted() in the chain down
	* to the one where the interaction was found.
	*/
	void acceptInteraction(int i);

	/**
	* Make this and all its anceseters incoming.
	*/
	void makeIncoming();

	/**
	* Indicate that the interaction that has been tested was rejected
	* by simply removing the interaction root inserted by
	* insertInteraction().
	*/
	void rejectInteraction(int i);


	/**
	* Change direction before producing collision.
	*/
	void mirror(double yframe);

	/**
	* Translate according to the impagt parameter.
	*/
	void translate(const ImpactParameters & b);

	/**
	* Prints all the member variables to cerr.
	*/
	void debugme();

	/**
	* Debug the structure of the shadow tree.
	*/
	void debugTree(string indent);

	//@}

	public:

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

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

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

	private:

	/**
	* The original DIPSY parton
	*/
	tPartonPtr theOriginal;

	/**
	- * The position in impact parameter space.
	- */
	- Point thePosition;
	-
	- /**
	* The positive light-cone momentum.
	*/
	Energy thePlus;

	/**
	* The transverse momentum.
	*/
	TransverseMomentum thePT;

	/**
	* The transverse momentum generated in the emission.
	*/
	TransverseMomentum thePT0;

	/**
	* The mass of this parton.
	*/
	mutable Energy theMass;

	/**
	* The accumulated negative light-cone momentum deficit.
	*/
	Energy theMinus;

	/**
	* The rapidity.
	*/
	double theY;

	/**
	* The rapidity generated in the emission.
	*/
	double theY0;

	/**
	* The flavour of this parton.
	*/
	long theFlavour;

	/**
	- * Indicate what happened to the previous instance of this parton.
	- */
	- Difference theDifference;
	-
	- /**
	* The state of this parton, if any, before the swing or the emission of this
	* parton's sibling. Is aways null if theParent exists.
	*/
	tSPartonPtr thePrevious;

	/**
	* The parent parton if this has emitted or swinged. Is aways null if
	* thePrevious exists.
	*/
	tSPartonPtr theParent;

	/**
	* The sibling of this parton if any.
	*/
	tSPartonPtr theSibling;

	/**
	* The state of this parton after the emission of the child. Is null
	* if this parton has not emitted.
	*/
	SPartonPtr theNext;

	/**
	* The child this parton has emitted, if any.
	*/
	SPartonPtr theChild;

	/**
	* The dipoles connected to this parton.
	*/
	tDipolePair theDipoles;

	/**
	* The partons connected to this.
	*/
	tSPartonPair theNeighbors;

	/**
	* The emission factor for this parton if emitted or emitter.
	*/
	InvEnergy2 theEmissionFactor;

	/**
	* The resolution factor for this parton if emitted or emitter.
	*/
	InvEnergy2 theRes;

	/**
	* Indicate that the corresponding emission is on a protected
	* propagator. Any other propagator must assume that this emission
	* must be resolved.
	*/
	bool isLocked;

	/**
	* Indicate if this parton has interacted.
	*/
	bool hasInteracted;

	/**
	* Indicate if this parton is on shell, that is if it has
	* been supplied with the neccesary p+ or p- for left or rightmoving
	* particles respectively.
	*/
	bool isOnShell;

	/**
	* This is root of an interaction chain.
	*/
	set<int> theInteractionRoots;

	/**
	* The interactions this was directly involved with.
	*/
	set<int> theInteractions;

	/**
	* Indicate if this parton is a valence parton.
	*/
	bool isValence;

	protected:

	/**
	* Exception class for bad kinematics.
	*/
	struct ShadowPartonKinematicsException: public Exception {};

	public:

	bool memememe;

	private:

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

	};

	}

	#endif /* DIPSY_ShadowParton_H */

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 7:35 PM (1 d, 10 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3805850
Default Alt Text: (82 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions