No OneTemporary
Actions

Size

80 KB

Subscribers

None

View Options

	diff --git a/Shower/Powheg/Matching/PowhegHandler.cc b/Shower/Powheg/Matching/PowhegHandler.cc
	--- a/Shower/Powheg/Matching/PowhegHandler.cc
	+++ b/Shower/Powheg/Matching/PowhegHandler.cc
	@@ -1,1788 +1,1807 @@
	// -- C++ --
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the PowhegHandler class.
	//

	#include "PowhegHandler.h"
	#include "ThePEG/Utilities/CFileLineReader.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "Herwig++/Shower/Base/KinematicsReconstructor.h"
	#include "Herwig++/Shower/Base/PartnerFinder.h"
	#include "Herwig++/Shower/Base/MECorrectionBase.h"
	#include "Herwig++/Utilities/Histogram.h"
	#include "QTildeSudakovIntegrator.h"
	#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
	#include "ThePEG/PDT/EnumParticles.h"
	#include "Herwig++/PDF/HwRemDecayer.h"
	#include <queue>
	#include "ThePEG/Repository/UseRandom.h"


	using namespace Herwig;

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

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

	void PowhegHandler::persistentOutput(PersistentOStream & os) const {
	os << _alphaS << _sudopt << _sudname << _jetMeasureMode << _allowedInitial
	<< _allowedFinal << _matrixElement << _lepton << _highestMult << _reweightOff << _testSudakovs <<_yini
	<< _alphaSMG << _npoint << ounit( _max_qtilde, GeV ) << ounit( _max_pt_cut, GeV )
	<< ounit( _min_pt_cut, GeV ) << _clusterOption;
	}

	void PowhegHandler::persistentInput(PersistentIStream & is, int) {
	is >> _alphaS >> _sudopt >> _sudname >> _jetMeasureMode >> _allowedInitial
	>> _allowedFinal >> _matrixElement >> _lepton >> _highestMult >> _reweightOff >> _testSudakovs >> _yini
	>> _alphaSMG >> _npoint >> iunit( _max_qtilde, GeV ) >> iunit( _max_pt_cut, GeV )
	>> iunit( _min_pt_cut, GeV ) >> _clusterOption;
	}

	ClassDescription<PowhegHandler> PowhegHandler::initPowhegHandler;
	// Definition of the static class description member.

	void PowhegHandler::Init() {

	static ClassDocumentation<PowhegHandler> documentation
	("The PowhegHandler class manages the implementation of the CKKW approach using"
	"the truncated shower.");

	static Reference<PowhegHandler,ShowerAlpha> interfaceShowerAlpha
	("ShowerAlpha",
	"The object calculating the strong coupling constant",
	&PowhegHandler::_alphaS, false, false, true, false, false);

	static Switch<PowhegHandler,unsigned int> interfaceSudakovOption
	("SudakovOption",
	"Option for the initialisation of the Sudakov tables",
	&PowhegHandler::_sudopt, 0, false, false);
	static SwitchOption interfaceSudakovOptionWrite
	(interfaceSudakovOption,
	"Write",
	"Calculate the Sudakov and write the table to a file",
	1);
	static SwitchOption interfaceSudakovOptionRead
	(interfaceSudakovOption,
	"Read",
	"Read the Sudakov table from a file",
	2);
	static SwitchOption interfaceSudakovOptionCompute
	(interfaceSudakovOption,
	"Compute",
	"Calculate the Sudakov but don't write the table",
	0);

	static Parameter<PowhegHandler,string> interfaceSudakovName
	("SudakovName",
	"Name for the file containing the Sudakov form factors",
	&PowhegHandler::_sudname, "sudakov.data",
	false, false);

	static Switch<PowhegHandler, unsigned int> ifaceJetMeasureMode
	("JetMeasure",
	"Choice of the jet measure algorithm",
	&PowhegHandler::_jetMeasureMode, 1, false, false);

	static SwitchOption Durham
	(ifaceJetMeasureMode,"Durham","Durham jet algorithm", 0);

	static SwitchOption ShowerPt
	(ifaceJetMeasureMode,"ShowerPt","ShowerPt", 1);

	static SwitchOption LUCLUS
	(ifaceJetMeasureMode,"LUCLUS","LUCLUS jet algorithm", 2);

	static Parameter<PowhegHandler,double> interfaceMergeScale
	("MergeScale",
	"The CKKW merging scale, yini",
	&PowhegHandler::_yini, 0.001, 0.0, 1.0,
	false, false, Interface::limited );

	static Parameter<PowhegHandler,double> interfaceAlphaSMG
	("alphaSMG",
	"The fixed alphas used in MG event generation",
	&PowhegHandler::_alphaSMG, 0.118, 0.0, 1.0,
	false, false, Interface::limited );

	static Switch<PowhegHandler,bool> interfaceLepton
	("Lepton",
	"Whether is a hadron-hadron or lepton-lepton collision",
	&PowhegHandler::_lepton, true, false, false);
	static SwitchOption interfaceLeptonLeptonic
	(interfaceLepton,
	"Leptonic",
	"Leptonic collision",
	true);
	static SwitchOption interfaceLeptonHadronic
	(interfaceLepton,
	"Hadronic",
	"Hadronic collision",
	false);

	static Switch<PowhegHandler,bool> interfaceHighestMultiplicity
	("HighestMult",
	"Whether we are treating the highest mutliplicity treatment",
	&PowhegHandler::_highestMult, false, false, false);

	static SwitchOption interfaceMultHighest
	(interfaceHighestMultiplicity,
	"Highest",
	"highest multiplicity",
	true);

	static SwitchOption interfaceMultNotHighest
	(interfaceHighestMultiplicity,
	"NotHighest",
	"Not the highest multiplicity",
	false);

	static Switch<PowhegHandler,bool> interfaceReweight
	("ReweightOff",
	"Whether to switch off the sudakov reweighting",
	&PowhegHandler::_reweightOff, false, false, false);
	static SwitchOption interfaceReweightOff
	(interfaceReweight,
	"Off",
	"No Sudakov reweighting",
	true);
	static SwitchOption interfaceReweightOn
	(interfaceReweight,
	"On",
	"Do Sudakov reweighting",
	false);

	static Switch<PowhegHandler,bool> interfaceTestSudakov
	("testSudakov",
	"Whether to output Sudakov test histograms",
	&PowhegHandler::_testSudakovs, false, false, false);

	static SwitchOption interfaceTestSudakovOff
	(interfaceTestSudakov,
	"Off",
	"No Sudakov testing",
	false);
	static SwitchOption interfaceTestSudakovOn
	(interfaceTestSudakov,
	"On",
	"Do Sudakov testing",
	true);


	static Reference<PowhegHandler,MEBase> interfaceMatrixElement
	("MatrixElement",
	"The matrix element class for the core 2->2 process",
	&PowhegHandler::_matrixElement, false, false, true, true, false);

	static Parameter<PowhegHandler,unsigned int> interfaceInterpPoints
	("InterpolatorPoints",
	"The number of points used for sudakov interpolation tables",
	&PowhegHandler::_npoint, 10, 0, 1000000,
	false, false, Interface::limited );

	static Parameter<PowhegHandler, Energy> interfaceMaxQTilde
	("maxQTilde",
	"The maximum QTilde scale for sudakov interpolation tables",
	&PowhegHandler::_max_qtilde, GeV, 91.2GeV, 1.GeV, 1000000.*GeV,
	false, false, Interface::limited);

	static Parameter<PowhegHandler, Energy> interfaceMaxPtCut
	("maxPtCut",
	"The maximum pt cut for sudakov interpolation tables",
	&PowhegHandler::_max_pt_cut, GeV, 45.6GeV, 1.GeV, 1000000.*GeV,
	false, false, Interface::limited);

	static Parameter<PowhegHandler, Energy> interfaceMinPtCut
	("minPtCut",
	"The minimum pt cut for sudakov interpolation tables",
	&PowhegHandler::_min_pt_cut, GeV, 0.GeV, 0.GeV, 1000000.*GeV,
	false, false, Interface::limited);

	static Switch<PowhegHandler, unsigned int> ifaceClusterOption
	("ClusterOption",
	"Choice of the clustering scheme",
	&PowhegHandler::_clusterOption, 0, false, false);

	static SwitchOption angularOrdered
	(ifaceClusterOption,"angularOrdered","make all histories and require angular ordering", 0);

	static SwitchOption jetClustered
	(ifaceClusterOption,"jetClustered", "cluster according to jet algorithm", 1);

	}

	double PowhegHandler::reweightCKKW(int minMult, int maxMult) {
	// cluster the event

	if( _clusterOption == 0 )
	_theHardTree = doClusteringOrdered();
	else
	_theHardTree = doClustering();

	// return if fails
	if(!_theHardTree)
	return 0.;

	// compute the Sudakov weight
	double SudWgt = _lepton ? sudakovWeight( _theHardTree ) : 1.;

	//update the sub process
	if(_lepton) {
	ParticleVector outgoing = lastXCombPtr()->subProcess()->outgoing();
	for(unsigned int ix=0;ix<outgoing.size();++ix) {
	lastXCombPtr()->subProcess()->removeEntry(outgoing[ix]);
	tParticleVector parents=outgoing[ix]->parents();
	for(unsigned int iy=0;iy<parents.size();++iy)
	parents[iy]->abandonChild(outgoing[ix]);
	}
	// add new ones based on the HardTree
	map<ColinePtr,ColinePtr> colourMap;
	for(set<HardBranchingPtr>::const_iterator it=_theHardTree->branchings().begin();
	it!=_theHardTree->branchings().end();++it) {
	if((**it).incoming()) continue;
	PPtr newParticle = new_ptr(Particle((**it).branchingParticle()->dataPtr()));
	newParticle->set5Momentum((**it).showerMomentum());
	//do colour connections
	if((**it).branchingParticle()->colourLine()) {
	map<ColinePtr,ColinePtr>::iterator loc
	= colourMap.find((**it).branchingParticle()->colourLine());
	if(loc!=colourMap.end()) loc->second->addColoured(newParticle);
	else {
	ColinePtr newLine=new_ptr(ColourLine());
	colourMap[(**it).branchingParticle()->colourLine()]=newLine;
	newLine->addColoured(newParticle);
	}
	}
	if((**it).branchingParticle()->antiColourLine()) {
	map<ColinePtr,ColinePtr>::iterator loc
	= colourMap.find((**it).branchingParticle()->antiColourLine());
	if(loc!=colourMap.end()) loc->second->addAntiColoured(newParticle);
	else {
	ColinePtr newLine=new_ptr(ColourLine());
	colourMap[(**it).branchingParticle()->antiColourLine()]=newLine;
	newLine->addAntiColoured(newParticle);
	}
	}
	lastXCombPtr()->subProcess()->addOutgoing(newParticle);
	}
	}
	else {
	set<HardBranchingPtr>::const_iterator it;
	map<ColinePtr,ColinePtr> colourMap;
	ParticleVector outgoing;
	PPair incoming;
	for(it=_theHardTree->branchings().begin();
	it!=_theHardTree->branchings().end();++it) {
	PPtr newParticle = new_ptr(Particle((**it).branchingParticle()->dataPtr()));
	newParticle->set5Momentum((**it).showerMomentum());
	if((**it).branchingParticle()->colourLine()) {
	map<ColinePtr,ColinePtr>::iterator loc
	= colourMap.find((**it).branchingParticle()->colourLine());
	if(loc!=colourMap.end()) loc->second->addColoured(newParticle);
	else {
	ColinePtr newLine=new_ptr(ColourLine());
	colourMap[(**it).branchingParticle()->colourLine()]=newLine;
	newLine->addColoured(newParticle);
	}
	}
	if((**it).branchingParticle()->antiColourLine()) {
	map<ColinePtr,ColinePtr>::iterator loc
	= colourMap.find((**it).branchingParticle()->antiColourLine());
	if(loc!=colourMap.end()) loc->second->addAntiColoured(newParticle);
	else {
	ColinePtr newLine=new_ptr(ColourLine());
	colourMap[(**it).branchingParticle()->antiColourLine()]=newLine;
	newLine->addAntiColoured(newParticle);
	}
	}
	if((**it).incoming()) {
	if(lastXCombPtr()->subProcess()->incoming().first->momentum().z()/
	newParticle->momentum().z()>0.)
	incoming.first = newParticle;
	else
	incoming.second = newParticle;
	}
	else
	outgoing.push_back(newParticle);
	}
	SubProPtr newSubProcess=
	new_ptr(SubProcess(incoming,
	lastXCombPtr()->subProcess()->collision(),
	lastXCombPtr()->subProcess()->handler()));
	for(unsigned int ix=0;ix<outgoing.size();++ix)
	newSubProcess->addOutgoing(outgoing[ix]);
	lastXCombPtr()->subProcess(newSubProcess);
	}
	if( ! _reweightOff ){
	if ( SudWgt / 5. > 1 ) cerr << SudWgt <<"\n";
	return SudWgt / 5.;
	}

	else{
	return 1.;
	}
	}
	void PowhegHandler::dofinish() {
	ShowerHandler::dofinish();
	+
	+ if( _clusterOption == 1 ){
	+ cout<<"\n---------\n"
	+ <<"proportion of unordered trees created = "
	+ << ( 1. - double( _ordered_trees_created )
	+ / double( _trees_created ) ) * 100.
	+ <<" %\n--------\n";
	+ }
	+
	string fname = generator()->filename() + string("-")
	+ string("wgts.top");
	ofstream output(fname.c_str());

	using namespace HistogramOptions;

	_hSud->topdrawOutput(output,Frame,
	"RED",
	"Sudakov wgts",
	"",
	"freq",
	"",
	"wgt",
	"");


	_halphaS->topdrawOutput(output,Frame,
	"RED",
	"AlphaS wgts",
	"",
	"freq",
	"",
	"wgt",
	"");
	}


	void PowhegHandler::doinitrun() {
	+ _trees_created = 0;
	+ _ordered_trees_created = 0;
	+
	ShowerHandler::doinitrun();
	_s = sqr( generator()->maximumCMEnergy() );

	_hSud = new_ptr(Histogram(0.,2.,100));
	_halphaS = new_ptr(Histogram(0.,2.,100));

	ofstream sudFileOutput;
	if(_sudopt==1) sudFileOutput.open(_sudname.c_str());

	// integrator for the outer integral
	GaussianIntegrator outer;
	// get the final-state branchings from the evolver
	if(_sudopt!=2) {

	for(BranchingList::const_iterator
	it = evolver()->splittingGenerator()->finalStateBranchings().begin();
	it != evolver()->splittingGenerator()->finalStateBranchings().end(); ++it) {

	//skip sudakovs involving tops
	if( abs( it->second.second[0] ) == 6 \|\|
	abs( it->second.second[1] ) == 6 \|\|
	abs( it->second.second[2] ) == 6 ) continue;

	Ptr<QTildeSudakovIntegrator>::pointer integrator =
	new_ptr( QTildeSudakovIntegrator(it->second, _jetMeasureMode, _s ) );

	Energy qtildemax = _max_qtilde;
	Energy qtildemin = integrator->minimumScale();

	//initialise sudakov values on grid ij ( pt_i, qtilde_j )
	vector< double > dummy( _npoint, 0. );
	vector< vector< double > > sud( _npoint, dummy );
	vector< Energy > ptCut;
	vector< Energy > scale;

	//fill scales at start
	for( unsigned int ix = 0; ix < _npoint; ++ix ){
	ptCut.push_back( _min_pt_cut + double( ix ) * ( _max_pt_cut - _min_pt_cut ) / double( _npoint - 1 ) );
	scale.push_back( qtildemin + double( ix ) * ( qtildemax - qtildemin ) / double( _npoint - 1 ) );
	}
	//fill sud integrals
	for( unsigned int ix = 0; ix < _npoint; ++ix ){
	sud[ix][0] = 0.;
	for( unsigned int jx = 1; jx < _npoint; ++jx ) {
	double currentSud = integrator->value( scale[ jx ], scale[ jx - 1 ], ptCut[ix] );
	sud[ix][jx] = ( sud[ix][ jx - 1 ] + currentSud );
	}
	}
	//exponentiate to the Sudakov
	for( unsigned int ix = 0; ix < _npoint; ++ix ) {
	for( unsigned int jx = 0; jx < _npoint; ++jx ) {
	sud[ix][jx] = exp( - sud[ix][jx] );
	}
	}

	Interpolator2d< double, Energy, Energy >::Ptr theInterpolator =
	new_ptr( Interpolator2d<double,Energy,Energy>( sud, ptCut, scale ) );

	_fbranchings.insert( make_pair( it->first, make_pair( theInterpolator, qtildemin ) ) );

	//write current sud grid to selected output file
	if(_sudopt==1) {

	sudFileOutput << it->second.second[0] << "\t"
	<< it->second.second[1] << "\t"
	<< it->second.second[2] << "\n";

	//output the grid size i * j
	sudFileOutput << ptCut.size() << "\t" << scale.size() << "\n";

	for( unsigned int jx = 0;jx < scale.size(); ++jx ){
	//write a row
	for( unsigned int ix = 0; ix < ptCut.size(); ++ix ){
	sudFileOutput << sud[ix][jx] << "\t";
	}
	sudFileOutput << "\n";
	}
	//output pts in a line
	for( unsigned int ix = 0; ix < ptCut.size(); ++ix ){
	sudFileOutput << ptCut[ix] / GeV << "\t";
	}
	sudFileOutput << "\n";
	//output scales in a line
	for( unsigned int ix = 0; ix < scale.size(); ++ix ){
	sudFileOutput << scale[ix] / GeV << "\t";
	}
	sudFileOutput << "\n";
	}

	}
	sudFileOutput.close();
	}

	else {
	CFileLineReader file(_sudname);
	while(file.readline()) {
	string line = file.getline();
	istringstream is;
	is.str(line);
	IdList ids(3);
	//GET NAMES AND IDS FROM FIRST LINE
	is >> ids[0] >> ids[1] >> ids[2];

	file.readline();
	unsigned int isize;
	unsigned int jsize;

	//GET THE NUMBER OF POINTS IN THIS
	is.str(file.getline());
	is >> isize >> jsize;

	//initialise vectors and matrix to the correct size
	vector< double > dummy( jsize, 0. );
	vector< vector< double > > sud( isize, dummy );
	vector< Energy > pt( isize );
	vector< Energy > scale( jsize );

	//read in matrix of sud values sud( pt_i, scale_i )

	//read a column- different line is different qtilde
	for( unsigned int jx = 0; jx < jsize; ++jx ) {
	//read a horizontal line diferent entry ids different pt
	file.readline();
	is.str(file.getline());
	for( unsigned int ix = 0; ix < isize; ++ix )
	is >> sud[ix][jx];
	}
	//read in pt values
	file.readline();
	is.str(file.getline());
	for( unsigned int ix = 0; ix < isize; ++ix ){
	double val;
	is >> val;
	pt[ix] = val * GeV;
	}
	file.readline();
	is.str(file.getline());
	for( unsigned int jx = 0; jx < jsize; ++jx ){
	double val;
	is >> val;
	scale[jx] = val * GeV;
	}

	//find branching list with matching ids
	BranchingList::const_iterator it,
	start = evolver()->splittingGenerator()->finalStateBranchings().lower_bound(ids[0]),
	end = evolver()->splittingGenerator()->finalStateBranchings().upper_bound(ids[0]);
	for( it = start; it != end; ++it ) {
	if( it->second.second[0] == ids[0] &&
	it->second.second[1] == ids[1] &&
	it->second.second[2] == ids[2] ) {
	// construct the Interpolators
	Interpolator2d< double, Energy, Energy >::Ptr theInterpolator =
	new_ptr( Interpolator2d< double, Energy, Energy >( sud, pt, scale ) );
	Energy qtildemin = scale[0];
	_fbranchings.insert( make_pair( it->first, make_pair( theInterpolator, qtildemin ) ) );
	break;
	}

	if( it == end ) {
	cerr << "sud read error: could not find correct branching list\n";
	}
	}
	}
	}
	if( _testSudakovs ) testSudakovs();
	}

	void PowhegHandler::doinit() throw(InitException) {
	ShowerHandler::doinit();
	// extract the allowed branchings
	// final-state
	for(BranchingList::const_iterator
	it = evolver()->splittingGenerator()->finalStateBranchings().begin();
	it != evolver()->splittingGenerator()->finalStateBranchings().end(); ++it) {
	pair<long,long> prod(make_pair(it->second.second[1],it->second.second[2]));
	_allowedFinal.insert(make_pair(prod,it->second));
	swap(prod.first,prod.second);
	_allowedFinal.insert(make_pair(prod,it->second));
	}
	// initial-state
	for(BranchingList::const_iterator
	it = evolver()->splittingGenerator()->initialStateBranchings().begin();
	it != evolver()->splittingGenerator()->initialStateBranchings().end(); ++it) {
	_allowedInitial.insert(make_pair(it->second.second[0],it->second));
	}

	}

	void PowhegHandler::cascade() {
	ShowerHandler::cascade();
	}

	double PowhegHandler::getJetMeasure(ShowerParticlePtr part_i,
	ShowerParticlePtr part_j){
	double yij;
	double costheta = part_i->momentum().vect().dot( part_j->momentum().vect() )
	/ part_i->momentum().vect().mag() / part_j->momentum().vect().mag();
	switch( _jetMeasureMode ){
	case 0:
	if( sqr( part_i->momentum().e() ) > sqr( part_j->momentum().e() ) )
	yij = 2. * sqr( part_j->momentum().e() ) * ( 1. - costheta ) / _s ;
	else
	yij = 2. * sqr( part_i->momentum().e() ) * ( 1. - costheta ) / _s ;
	break;
	case 2:
	yij = 2. * sqr( part_i->momentum().e() * part_j->momentum().e() /
	( part_i->momentum().e() + part_j->momentum().e() ) )/_s
	* ( 1. - costheta );
	break;
	default:
	yij = 1.;
	break;
	}
	return yij;
	}

	//given two particles returns value of durham jet algorithm
	bool PowhegHandler::splittingAllowed( ShowerParticlePtr part_i,
	ShowerParticlePtr part_j,
	int qq_pairs ) {
	// g 2 q qbar or an incorrect qq type
	if ( abs ( part_i->id() ) < 7 && abs ( part_j->id() ) < 7 ) {
	if ( abs ( part_i->id() ) != abs ( part_j->id() ) ) return false;
	if ( ( part_i->id() < 0 && part_j->id() < 0 ) \|\|
	( part_i->id() > 0 && part_j->id() > 0 ) ) return false;
	if ( qq_pairs < 2 ) return false;
	}
	return true;
	}

	// finds the id of the emitting particle and sudakov for the desired clustering
	// also swaps order of children pointers as required (not pointers passed by reference)
	SudakovPtr PowhegHandler::getSud( long & emmitter_id,
	ShowerParticlePtr & part_i,
	ShowerParticlePtr & part_j ) {
	// g 2 q qbar or an incorrect qq type
	if ( abs ( part_i->id() ) < 7 && abs ( part_j->id() ) < 7 ) {
	if ( abs ( part_i->id() ) != abs ( part_j->id() ) ) return SudakovPtr();
	if ( ( part_i->id() < 0 && part_j->id() < 0 ) \|\|
	( part_i->id() > 0 && part_j->id() > 0 ) ) return SudakovPtr();
	//if the q and qbar are the wrong way round then switch order
	if ( part_j->id() > part_i->id() ) swap( part_i, part_j );
	emmitter_id = 21;
	}
	// q/qbar 2 q/qbar g
	else if ( abs ( part_i->id() ) < 7 \|\| abs ( part_j->id() ) < 7 ) {
	if( abs ( part_i->id() ) < 7 ){
	emmitter_id = part_i->id();
	}
	else {
	emmitter_id = part_j->id();
	}
	}
	// g 2 g g
	else {
	emmitter_id = 21;
	}
	BranchingList branchings =
	evolver()->splittingGenerator()->finalStateBranchings();

	//cycle through list of all branchings with the correct abs ( emmitter_id )
	for(BranchingList::const_iterator cit = branchings.lower_bound( abs(emmitter_id) );
	cit != branchings.upper_bound( abs(emmitter_id) ); ++cit ) {
	IdList ids = cit->second.second;
	if( abs( ids[0] ) == abs( emmitter_id ) ) {
	if( abs(ids[1]) == abs(part_i->id()) &&
	abs(ids[2]) == abs(part_j->id()) ) {
	return cit->second.first;
	}
	if( abs( ids[1] ) == abs( part_j->id() ) &&
	abs( ids[2] ) == abs( part_i->id() ) ) {
	swap( part_i, part_j );
	return cit->second.first;
	}
	}
	}
	return SudakovPtr();
	}
	//have the & in here so that can remove pointer - is this right??
	bool PowhegHandler::fillProtoTrees( map< ShowerParticlePtr, HardBranchingPtr > theParticles,
	ProtoTreePtr currentProtoTree ){
	if( theParticles.size() < 3 ) return true;
	//find number of qqpairs - count quarks - used
	int no_qqbar = 0;
	for( map<ShowerParticlePtr, HardBranchingPtr>::iterator itc = theParticles.begin();
	itc != theParticles.end(); itc++ )
	if( itc->first->id() > 0 && itc->first->id() < 7 ) no_qqbar ++;

	HardBranchingPtr currentBranching;
	ProtoTreePtr newProtoTree;
	map< ShowerParticlePtr, HardBranchingPtr > newParticles;
	for( map< ShowerParticlePtr, HardBranchingPtr >::iterator ita = theParticles.begin();
	ita != theParticles.end(); ita++ ) {
	for( map< ShowerParticlePtr, HardBranchingPtr >::iterator itb = theParticles.begin();
	itb != ita; itb++) {
	if( ! splittingAllowed( ita->first, itb->first, no_qqbar ) ) continue;

	currentBranching = getCluster( make_pair( ita->first, itb->first ),
	theParticles );

	set< HardBranchingPtr > newTreeBranchings = currentProtoTree->getBranchings();

	newProtoTree = new_ptr( ProtoTree( newTreeBranchings ) );
	//now remove the hard branchings of clustered particles from newTree

	//currentProtoTree should contain which ever hardbranchings are in theParticles
	if( ! newProtoTree->removeBranching( ita->second ) )
	cerr<<"fill proto tree problem!!! can't find clustered in newProtoTree 1 \n"
	<< "couldn't find: " << ita->second <<"\n";
	if( ! newProtoTree->removeBranching( itb->second ) )
	cerr<<"fill proto tree problem!!! can't find clustered in newProtoTree 1 \n"
	<< "couldn't find: " << itb->second <<"\n";

	newProtoTree->addBranching( currentBranching );

	newParticles = theParticles;
	if( newParticles.find( ita->first ) != newParticles.end() )
	newParticles.erase( ita->first );
	else
	cerr<<"fill newParticles problem!!! can't find clustered in newParticles \n";
	if( newParticles.find( itb->first ) != newParticles.end() )
	newParticles.erase( itb->first );
	else
	cerr<<"fill newParticles problem!!! can't find clustered in newParticles \n";
	newParticles.insert( make_pair( currentBranching->branchingParticle(),
	currentBranching ) );

	if( ! repeatProtoTree( newProtoTree ) ) _proto_trees.insert( newProtoTree );

	//remove the current tree if it hasn't already been removed
	if( _proto_trees.find( currentProtoTree ) != _proto_trees.end() )
	_proto_trees.erase( currentProtoTree );
	//do recursion
	fillProtoTrees( newParticles, newProtoTree );
	}
	}
	return true;
	}

	HardBranchingPtr PowhegHandler::getCluster( pair< ShowerParticlePtr, ShowerParticlePtr > clusterPair,
	map< ShowerParticlePtr, HardBranchingPtr > theParticles ){
	for( map< HardBranchingPtr , pair< ShowerParticlePtr, ShowerParticlePtr > >::const_iterator
	cit = _all_clusters.begin(); cit != _all_clusters.end(); ++cit ){
	if( ( cit->second.first == clusterPair.first && cit->second.second == clusterPair.second ) \|\|
	( cit->second.first == clusterPair.second && cit->second.second == clusterPair.first ) ){
	return cit->first;
	}
	}
	//branching not found create -- with sudakov and everything
	long thePartId;
	tcPDPtr particle_data;
	SudakovPtr theSudakov = getSud( thePartId, clusterPair.first, clusterPair.second );
	if( !theSudakov ){
	cerr << "can't find the sudakov in: \n"
	<< *clusterPair.first<<"\n"
	<< *clusterPair.second<<"\n";
	}
	Lorentz5Momentum pairMomentum = clusterPair.first->momentum() +
	clusterPair.second->momentum();
	pairMomentum.setMass( 0. * MeV );
	particle_data = getParticleData( thePartId );

	//creates emitter particle
	ShowerParticlePtr clustered = new_ptr( ShowerParticle( particle_data, true ) );
	clustered->set5Momentum( pairMomentum );

	HardBranchingPtr clusteredBranch( new_ptr( HardBranching( clustered, theSudakov,
	HardBranchingPtr(), false ) ) );
	_all_clusters.insert( make_pair( clusteredBranch, clusterPair ) );

	//join children
	clusteredBranch->addChild( theParticles.find( clusterPair.first )->second );
	clusteredBranch->addChild( theParticles.find( clusterPair.second )->second );

	return clusteredBranch;
	}

	bool PowhegHandler::repeatProtoTree( ProtoTreePtr currentProtoTree ){
	//loop over all prototrees and see how many hardbranchings of curentProtoTree are found in each
	for( set< ProtoTreePtr >::const_iterator cit = _proto_trees.begin();
	cit != _proto_trees.end(); ++cit ){
	unsigned int no_matches = 0;
	for( set< HardBranchingPtr >::const_iterator ckt
	= currentProtoTree->getBranchings().begin(); ckt != currentProtoTree->getBranchings().end(); ckt++ ){
	if( (cit)->getBranchings().find( ckt ) != (*cit)->getBranchings().end() )
	no_matches++;
	}
	if( no_matches == currentProtoTree->getBranchings().size() ) return true;
	}
	return false;
	}

	bool PowhegHandler::simpleColConnections( ProtoTreePtr theProtoTree ){
	set< HardBranchingPtr > currentProtoTree = theProtoTree->getBranchings();
	if( currentProtoTree.size() != 2 ) {
	cerr<<"\n\nwrong size of proto tree: "
	<< currentProtoTree.size() <<"\n\n\n";
	return false;
	}
	//colourline to join up q and qbar
	ColinePtr newline = new_ptr( ColourLine() );
	for( set< HardBranchingPtr >::iterator it = currentProtoTree.begin();
	it != currentProtoTree.end(); ++it ){
	(*it)->branchingParticle()->resetColour();
	if( (*it)->branchingParticle()->dataPtr()->iColour()
	== PDT::Colour3 )
	newline->addColoured( (*it)->branchingParticle() );

	else if( (*it)->branchingParticle()->dataPtr()->iColour()
	== PDT::Colour3bar )
	newline->addAntiColoured( (*it)->branchingParticle() );

	else cerr<< "\n\n\nClustered back to gluon\n\n";
	}
	return true;
	}

	bool PowhegHandler::simpleColConnections( HardTreePtr theHardTree ){
	set< HardBranchingPtr > particles = theHardTree->branchings();
	//colourline to join up q and qbar
	ColinePtr newline = new_ptr( ColourLine() );
	for( set< HardBranchingPtr >::iterator it = particles.begin();
	it != particles.end(); ++it ){
	if( (*it)->incoming() ) continue;
	(*it)->branchingParticle()->resetColour();
	if( (*it)->branchingParticle()->dataPtr()->iColour()
	== PDT::Colour3 )
	newline->addColoured( (*it)->branchingParticle() );

	else if( (*it)->branchingParticle()->dataPtr()->iColour()
	== PDT::Colour3bar )
	newline->addAntiColoured( (*it)->branchingParticle() );
	else cerr<< "\n\n\nClustered back to gluon\n\n";
	}
	return true;
	}

	HardTreePtr PowhegHandler::doClusteringOrdered() {
	if(!_lepton) {
	return generalClustering();
	}
	ParticleVector theParts = lastXCombPtr()->subProcess()->outgoing();

	//initialise global variables
	_all_clusters.clear();
	_proto_trees.clear();
	_hardTrees.clear();

	//make an intermediate and add to subprocess if not read in
	if(lastXCombPtr()->subProcess()->intermediates().empty()) {
	return HardTreePtr();
	PPair theIncomings = lastXCombPtr()->subProcess()->incoming();
	//set intermediate to Z
	long intermediate_id = 23;
	PPtr theIntermediate = new_ptr( Particle( getParticleData( intermediate_id ) ) );
	theIntermediate->set5Momentum( theIncomings.first->momentum() +
	theIncomings.second->momentum() );
	//add the intermediate - parent/child relations should be updated
	lastXCombPtr()->subProcess()->addIntermediate( theIntermediate );
	cerr<<"added intermediate\n"
	<< *theIntermediate<<"\n";
	}

	PPtr vb = lastXCombPtr()->subProcess()->intermediates()[0];

	map <ShowerParticlePtr,HardBranchingPtr> theParticles;
	tcPDPtr particle_data;
	ShowerParticlePtr vBoson = new_ptr( ShowerParticle( *vb, 1, false, false ) );

	//loops through the FS particles and create hardBranchings
	for( unsigned int i = 0; i < theParts.size(); i++){
	ShowerParticlePtr currentParticle =
	new_ptr( ShowerParticle( *theParts[i], 1, true, false ) );
	HardBranchingPtr currentBranching = new_ptr( HardBranching( currentParticle, SudakovPtr(),
	HardBranchingPtr(),false ) );
	theParticles.insert( make_pair( currentParticle, currentBranching ) );
	}
	//create and initialise the first tree
	ProtoTreePtr initialProtoTree = new_ptr( ProtoTree() );
	for( map<ShowerParticlePtr, HardBranchingPtr>::iterator ita = theParticles.begin();
	ita != theParticles.end(); ita++ ){
	initialProtoTree->addBranching( ita->second );
	}
	_proto_trees.insert( initialProtoTree );

	fillProtoTrees( theParticles, initialProtoTree );

	double totalWeight = 0.;

	//create a hardtree from each proto tree and fill _hardTrees with angular ordered configs
	for( set< ProtoTreePtr >::const_iterator cit = _proto_trees.begin();
	cit != _proto_trees.end(); ++cit ){
	simpleColConnections( *cit );
	//vector boson branching
	vector<HardBranchingPtr> spaceBranchings;
	//all branchings
	vector<HardBranchingPtr> theBranchings;
	//fill theBranchings
	for( set< HardBranchingPtr >::const_iterator cjt = (*cit)->getBranchings().begin();
	cjt != (*cit)->getBranchings().end(); ++cjt )
	theBranchings.push_back( *cjt );
	spaceBranchings.push_back( new_ptr( HardBranching( vBoson, SudakovPtr(),
	HardBranchingPtr(),
	true ) ) );
	theBranchings.push_back( spaceBranchings.back() );
	HardTreePtr powhegtree = new_ptr( HardTree( theBranchings,
	spaceBranchings ) );
	// Calculate the shower variables
	// if momentum deconstruction fails then continue and ignore
	if( ! evolver()->showerModel()->kinematicsReconstructor()
	->deconstructDecayJets( powhegtree, evolver() ) ) continue;

	//only insert angular ordered hardTrees
	if( powhegtree->checkHardOrdering() ) {
	//find the wgt and fill _hardTrees map
	powhegtree->findNodes();
	double treeWeight = sudakovWeight( powhegtree );
	treeWeight *= splittingFnWeight( powhegtree );
	_hardTrees.push_back( make_pair( powhegtree, treeWeight ) );
	totalWeight += treeWeight;
	}
	}

	if( _hardTrees.empty() )
	return HardTreePtr();

	//now just choose a hardTree from
	long treeIndex;
	do{
	treeIndex = UseRandom::irnd( _hardTrees.size() );
	} while ( _hardTrees[ treeIndex ].second / totalWeight < UseRandom::rnd() );

	//re-do momentum deconstruction (gets overridden by other trees otherwise)
	simpleColConnections( _hardTrees[ treeIndex ].first );
	if( ! evolver()->showerModel()->kinematicsReconstructor()
	->deconstructDecayJets( _hardTrees[ treeIndex ].first, evolver() ) )
	cerr<<"\n\nproblem doing momentum decon in selected tree \n\n";

	// return HardTreePtr();
	return _hardTrees[ treeIndex ].first;
	}

	void PowhegHandler::fixColours(tPPtr parent, tPPtr child1, tPPtr child2) {
	// the different possible cases
	if(parent->dataPtr()->iColour()==PDT::Colour3&&
	child1->dataPtr()->iColour()==PDT::Colour3&&
	child2->dataPtr()->iColour()==PDT::Colour8) {
	child2->colourLine()->addColoured(parent);
	ColinePtr temp = child2->antiColourLine();
	temp->addColoured(child1);
	child1->colourLine()->join(temp);
	}
	else if(parent->dataPtr()->iColour()==PDT::Colour3&&
	child2->dataPtr()->iColour()==PDT::Colour3&&
	child1->dataPtr()->iColour()==PDT::Colour8) {
	child1->colourLine()->addColoured(parent);
	ColinePtr temp = child1->antiColourLine();
	temp->addColoured(child2);
	child2->colourLine()->join(temp);
	}
	else if(parent->dataPtr()->iColour()==PDT::Colour3bar&&
	child1->dataPtr()->iColour()==PDT::Colour3bar&&
	child2->dataPtr()->iColour()==PDT::Colour8) {
	child2->antiColourLine()->addAntiColoured(parent);
	ColinePtr temp = child1->antiColourLine();
	temp->addColoured(child2);
	child2->colourLine()->join(temp);
	}
	else if(parent->dataPtr()->iColour()==PDT::Colour3bar&&
	child2->dataPtr()->iColour()==PDT::Colour3bar&&
	child1->dataPtr()->iColour()==PDT::Colour8) {
	child1->antiColourLine()->addAntiColoured(parent);
	ColinePtr temp = child2->antiColourLine();
	temp->addColoured(child1);
	child1->colourLine()->join(temp);
	}
	else if(parent->dataPtr()->iColour()==PDT::Colour8&&
	child1->dataPtr()->iColour()==PDT::Colour8&&
	child2->dataPtr()->iColour()==PDT::Colour8) {
	if(UseRandom::rndbool(0.5)) {
	child1->colourLine()->addColoured(parent);
	child2->antiColourLine()->addAntiColoured(parent);
	ColinePtr temp = child1->antiColourLine();
	temp->addColoured(child2);
	child2->colourLine()->join(temp);
	}
	else {
	child2->colourLine()->addColoured(parent);
	child1->antiColourLine()->addAntiColoured(parent);
	ColinePtr temp = child2->antiColourLine();
	temp->addColoured(child1);
	child1->colourLine()->join(temp);
	}
	}
	else if(parent->dataPtr()->iColour()==PDT::Colour8&&
	child1->dataPtr()->iColour()==PDT::Colour3&&
	child2->dataPtr()->iColour()==PDT::Colour3bar) {
	child1->colourLine()->addColoured(parent);
	child2->antiColourLine()->addAntiColoured(parent);
	}
	else if(parent->dataPtr()->iColour()==PDT::Colour8&&
	child1->dataPtr()->iColour()==PDT::Colour3bar&&
	child2->dataPtr()->iColour()==PDT::Colour3) {
	child2->colourLine()->addColoured(parent);
	child1->antiColourLine()->addAntiColoured(parent);
	}
	else {
	throw Exception() << "Unknown colour in PowhegHandler::fixColours()"
	<< Exception::runerror;
	}
	}

	HardTreePtr PowhegHandler::generalClustering() {
	if(!_matrixElement)
	throw Exception() << "PowhegHandler::generalClustering()"
	<< " must have a MatrixElement object for the core "
	<< "2->2 process" << Exception::runerror;
	PPair incoming = lastXCombPtr()->subProcess()->incoming();
	ParticleVector outgoing = lastXCombPtr()->subProcess()->outgoing();
	_s = lastXCombPtr()->lastS();
	// queue with the prototype trees
	std::queue<PrototypeTree> potentialTrees;
	// the base tree we'll make the others from
	PrototypeTree root;
	ShowerParticlePtr newParticle =
	new_ptr(ShowerParticle(incoming.first->dataPtr(),false));
	newParticle->set5Momentum(incoming.first->momentum());
	root.incoming.insert(new_ptr(PrototypeBranching(newParticle)));
	newParticle = new_ptr(ShowerParticle(incoming.second->dataPtr(),false));
	newParticle->set5Momentum(incoming.second->momentum());
	root.incoming.insert(new_ptr(PrototypeBranching(newParticle)));
	for(set<PrototypeBranchingPtr>::const_iterator it=root.incoming.begin();
	it!=root.incoming.end();++it) {
	root.currentSpaceLike.insert(*it);
	}
	for(unsigned int ix=0;ix<outgoing.size();++ix) {
	newParticle = new_ptr(ShowerParticle(outgoing[ix]->dataPtr(),true));
	newParticle->set5Momentum(outgoing[ix]->momentum());
	root.outgoing.insert(new_ptr(PrototypeBranching(newParticle)));
	}
	potentialTrees.push(root);
	// store the final potential trees
	list<PrototypeTree> trees;
	while (!potentialTrees.empty()) {
	PrototypeTree current = potentialTrees.front();
	bool found(false);
	// potential final-final mergings
	set<PrototypeBranchingPtr>::iterator it,jt;
	for(it=current.outgoing.begin();it!=current.outgoing.end();++it) {
	jt = it;
	++jt;
	for( ; jt!=current.outgoing.end();++jt) {
	pair<PrototypeBranchingPtr,PrototypeBranchingPtr>
	branch = make_pair(it,jt);
	BranchingElement allowed = allowedFinalStateBranching(branch);
	if(!allowed.first) continue;
	// copy the tree
	PrototypeTree newTree = current;
	map<PrototypeBranchingPtr,PrototypeBranchingPtr> pmap = newTree.reset();
	branch.first = pmap[branch.first ];
	branch.second = pmap[branch.second];
	// make the new branching
	// new particle first
	tcPDPtr newData = getParticleData(allowed.second[0]);
	Lorentz5Momentum newMomentum(branch.first ->particle->momentum()+
	branch.second->particle->momentum());
	if(!newData->CC()\|\|
	(branch.first ->particle->id()==allowed.second[1]&&
	branch.second->particle->id()==allowed.second[2])) {
	newParticle = new_ptr(ShowerParticle(newData,true));
	}
	else {
	newParticle = new_ptr(ShowerParticle(newData->CC(),true));
	}
	newParticle->set5Momentum(newMomentum);
	// then the branching
	PrototypeBranchingPtr newBranching(new_ptr(PrototypeBranching(newParticle)));
	branch.first ->parent =newBranching;
	branch.second->parent =newBranching;
	newBranching->children.push_back(branch.first );
	newBranching->children.push_back(branch.second);
	newBranching->sudakov = allowed.first;
	newTree.outgoing.erase(branch.first );
	newTree.outgoing.erase(branch.second);
	newTree.outgoing.insert(newBranching);
	// jet measure
	newTree.scales.push_back(hadronJetMeasure(branch.first ->particle->momentum(),
	branch.second->particle->momentum()));
	// insert in the relevant list
	if(newTree.outgoing.size()==2) trees.push_back(newTree);
	else potentialTrees.push(newTree);
	found = true;
	}
	}
	// initial-final mergings
	for(it=current.outgoing.begin();it!=current.outgoing.end();++it) {
	for(jt=current.currentSpaceLike.begin();
	jt!=current.currentSpaceLike.end();++jt) {
	pair<PrototypeBranchingPtr,PrototypeBranchingPtr>
	branch = make_pair(jt,it);
	BranchingElement allowed = allowedInitialStateBranching(branch);
	if(!allowed.first) continue;
	// copy the tree
	PrototypeTree newTree = current;
	map<PrototypeBranchingPtr,PrototypeBranchingPtr> pmap = newTree.reset();
	branch.first = pmap[branch.first ];
	branch.second = pmap[branch.second];
	// make the new branching
	// new particle first
	tcPDPtr newData = getParticleData(allowed.second[1]);
	Lorentz5Momentum newMomentum(branch.first ->particle->momentum()-
	branch.second->particle->momentum());
	if(!newData->CC()\|\|
	(branch.first ->particle->id()==allowed.second[0]&&
	branch.second->particle->id()==allowed.second[2])) {
	newParticle = new_ptr(ShowerParticle(newData,false));
	}
	else {
	newParticle = new_ptr(ShowerParticle(newData->CC(),false));
	}
	newParticle->set5Momentum(newMomentum);
	// then the branching
	PrototypeBranchingPtr newBranching(new_ptr(PrototypeBranching(newParticle)));
	newBranching->parent = branch.first;
	branch.second->parent = branch.first;
	branch.first->children.push_back(newBranching);
	branch.first->children.push_back(branch.second);
	newBranching->parent->sudakov = allowed.first;
	newTree.currentSpaceLike.erase(branch.first );
	newTree.outgoing .erase(branch.second);
	newTree.currentSpaceLike.insert(newBranching);
	// jet measure
	newTree.scales.push_back(hadronJetMeasure(branch.second->particle->momentum(),
	branch.second->particle->momentum(),false));
	if(branch.first ->particle->momentum().z()/
	branch.second->particle->momentum().z()>0.) newTree.scales.back()-=0.001*MeV2;
	// insert in the relevant list
	if(newTree.outgoing.size()==2) trees.push_back(newTree);
	else potentialTrees.push(newTree);
	found = true;
	}
	}
	// treated one branching so pop from the queue
	if(!found) trees.push_back(current);
	potentialTrees.pop();
	}
	// check the core process is allowed using the matrix element
	// and remove ones which aren't allowed
	list<PrototypeTree>::iterator it=trees.begin(),jt;
	while(it!=trees.end()) {
	DiagPtr diagram = getDiagram(*it);
	if(!diagram) it = trees.erase(it);
	else {
	it->diagram = diagram;
	++it;
	}
	}
	// finally for the moment select the one with the smallest pt for the first branching
	// now find the one with the minimum pt
	HardTreePtr newTree;
	while(!trees.empty()) {
	jt=trees.end();
	Energy2 minkT =1e30*GeV2;
	for(it=trees.begin();it!=trees.end();++it) {
	if(it->scales.back()<minkT) {
	minkT = it->scales.back();
	jt=it;
	}
	}
	// construct the hard tree
	newTree = (*jt).convert();
	// assign the beam particles
	setBeams(newTree);
	// construct the colour flow
	createColourFlow(newTree,jt->diagram);
	// Calculate the shower variables
	evolver()->showerModel()->kinematicsReconstructor()->
	deconstructDecayJets(newTree,evolver());
	if(checkTree(newTree)) break;
	trees.erase(jt);
	}
	// if no tree return an empty one
	if(trees.empty()) return HardTreePtr();
	// return the tree
	return newTree;
	}

	BranchingElement PowhegHandler::
	allowedFinalStateBranching(pair<PrototypeBranchingPtr,PrototypeBranchingPtr> & br) {
	// check with normal ID's
	pair<long,long> ptest = make_pair(br.first->particle->id(),br.second->particle->id());
	map<pair<long,long>,pair<SudakovPtr,IdList> >::const_iterator
	split = _allowedFinal.find(ptest);
	if(split!=_allowedFinal.end()) {
	if(split->second.second[1]!=ptest.first) swap(br.first,br.second);
	return split->second;
	}
	// check with CC
	if(br.first ->particle->dataPtr()->CC()) ptest.first *= -1;
	if(br.second->particle->dataPtr()->CC()) ptest.second *= -1;
	_allowedFinal.find(ptest);
	if(split!=_allowedFinal.end()) {
	if(split->second.second[1]!=ptest.first) swap(br.first,br.second);
	return split->second;
	}
	// not found found null pointer
	return make_pair(SudakovPtr(),IdList());
	}

	BranchingElement PowhegHandler::
	allowedInitialStateBranching(pair<PrototypeBranchingPtr,PrototypeBranchingPtr> & br) {
	// veto top
	if(abs(br.first ->particle->id())==ParticleID::t\|\|
	abs(br.second->particle->id())==ParticleID::t)
	return make_pair(SudakovPtr(),IdList());
	bool cc = br.first->particle->id()<0;
	pair<multimap<long, pair<SudakovPtr,IdList> >::const_iterator,
	multimap<long, pair<SudakovPtr,IdList> >::const_iterator>
	location = _allowedInitial.equal_range(abs(br.first->particle->id()));
	for(multimap<long, pair<SudakovPtr,IdList> >::const_iterator it=location.first;
	it!=location.second;++it) {
	long idtest = it->second.second[2];
	if(cc&&getParticleData(idtest)->CC()) idtest *= -1;
	if(idtest==br.second->particle->id()) return it->second;
	if(idtest==-br.second->particle->id()&&
	!br.first->particle->dataPtr()->CC()) return it->second;
	}
	// not found found null pointer
	return make_pair(SudakovPtr(),IdList());
	}

	DiagPtr PowhegHandler::getDiagram(const PrototypeTree & tree) {
	// extract the incoming particles
	set<PrototypeBranchingPtr>::const_iterator it=tree.currentSpaceLike.begin();
	tcPDPair incoming;
	incoming.first = (**it).particle->dataPtr();
	++it;
	incoming.second = (**it).particle->dataPtr();
	// and the outgoing particles
	multiset<tcPDPtr> outgoing;
	for(it=tree.outgoing.begin();it!=tree.outgoing.end();++it)
	outgoing.insert((**it).particle->dataPtr());
	// see if the process is allowed
	for(MEBase::DiagramVector::const_iterator dt = _matrixElement->diagrams().begin();
	dt!=_matrixElement->diagrams().end();++dt) {
	const cPDVector partons=(**dt).partons();
	// check incoming particles
	if(!((incoming.first==partons[0]&&incoming.second==partons[1])\|\|
	(incoming.first==partons[1]&&incoming.second==partons[0]))) continue;
	// check the number of outgoing
	if(partons.size()!=tree.outgoing.size()+2) return DiagPtr();
	// check the outgoing
	multiset<tcPDPtr> otemp(outgoing);
	multiset<tcPDPtr>::iterator it;
	for(unsigned int ix=2;ix<partons.size();++ix) {
	it=otemp.find(partons[ix]);
	if(it!=otemp.end()) otemp.erase(it);
	}
	if(!otemp.empty()) continue;
	return *dt;
	}
	return DiagPtr();
	}

	Energy2 PowhegHandler::hadronJetMeasure(const Lorentz5Momentum & p1,
	const Lorentz5Momentum & p2,
	bool final) {
	Energy2 output;
	if(final) {
	double deltay = p1.rapidity()-p2.rapidity();
	double deltaphi = p1.phi()-p2.phi();
	if(deltaphi<-Constants::pi) deltaphi += Constants::twopi;
	if(deltaphi> Constants::pi) deltaphi -= Constants::twopi;
	double deltaR = sqr(deltay)+sqr(deltaphi);
	output = min(p1.perp2(),p2.perp2())*deltaR;
	}
	else {
	output = p1.perp2();
	}
	return output;
	}

	HardBranchingPtr PrototypeBranching::convert() {
	if(!particle) {
	cerr << "testing don't have particle for the branching shit" << "\n";
	exit(0);
	}
	// create the new particle
	HardBranchingPtr hard=new_ptr(HardBranching(particle,sudakov,
	tHardBranchingPtr(),
	!particle->isFinalState()));
	// and the children
	for(unsigned int ix=0;ix<children.size();++ix) {
	hard->addChild(children[ix]->convert());
	hard->children().back()->parent(hard);
	}
	return hard;
	}

	HardTreePtr PrototypeTree::convert() {
	vector<HardBranchingPtr> branchings,spacelike;
	set<PrototypeBranchingPtr>::const_iterator it,jt;
	// incoming lines and spacelike inot the hard process
	for(it=incoming.begin();it!=incoming.end();++it) {
	spacelike.push_back((**it).convert());
	HardBranchingPtr br(spacelike.back());
	while (!br->children().empty()) {
	for(unsigned int ix=0;ix<br->children().size();++ix) {
	if(br->children()[ix]->incoming()) {
	br = br->children()[ix];
	break;
	}
	}
	}
	branchings.push_back(br);
	}
	// outgoing particles
	for(it=outgoing.begin();it!=outgoing.end();++it) {
	branchings.push_back((**it).convert());
	}
	HardTreePtr newTree = new_ptr(HardTree(branchings,spacelike));
	return newTree;
	}

	map<PrototypeBranchingPtr,PrototypeBranchingPtr> PrototypeTree::reset() {
	map<PrototypeBranchingPtr,PrototypeBranchingPtr> output;
	set<PrototypeBranchingPtr> newOutgoing;
	set<PrototypeBranchingPtr> newIncoming;
	set<PrototypeBranchingPtr> newSpaceLike;
	set<PrototypeBranchingPtr>::iterator it,jt;
	for(it=incoming.begin();it!=incoming.end();++it) {
	PrototypeBranchingPtr newBr = (**it).reset(PrototypeBranchingPtr(),output);
	newIncoming.insert(newBr);
	PrototypeBranchingPtr br=newBr;
	while(!br->children.empty()) {
	for(unsigned int ix=0;ix<br->children.size();++ix) {
	if(!br->children[ix]->particle->isFinalState()) {
	br = br->children[ix];
	break;
	}
	}
	}
	newSpaceLike.insert(br);
	}
	for(it=outgoing.begin();it!=outgoing.end();++it) {
	newOutgoing.insert((**it).reset(PrototypeBranchingPtr(),output));
	}
	outgoing = newOutgoing;
	incoming = newIncoming;
	currentSpaceLike = newSpaceLike;
	return output;
	}

	PrototypeBranchingPtr PrototypeBranching::
	reset(PrototypeBranchingPtr newParent,
	map<PrototypeBranchingPtr,PrototypeBranchingPtr> & pmap) {
	PrototypeBranchingPtr output(new_ptr(PrototypeBranching(particle)));
	pmap[this] = output;
	output->sudakov = sudakov;
	output->parent = newParent;
	for(unsigned int ix=0;ix<children.size();++ix) {
	output->children.push_back(children[ix]->reset(output,pmap));
	}
	return output;
	}

	void PowhegHandler::createColourFlow(HardTreePtr tree,
	DiagPtr diagram) {
	// first construct a set of on-shell momenta for the hard collison
	vector<Lorentz5Momentum> meMomenta;
	vector<tcPDPtr> mePartonData;
	PVector particles;
	set<HardBranchingPtr>::const_iterator it;
	for(it=tree->branchings().begin();it!=tree->branchings().end();++it) {
	if((**it).incoming()) {
	meMomenta.push_back((**it).branchingParticle()->momentum());
	mePartonData.push_back((**it).branchingParticle()->dataPtr());
	particles.push_back((**it).branchingParticle());
	}
	}
	for(it=tree->branchings().begin();it!=tree->branchings().end();++it) {
	if(!(**it).incoming()) {
	meMomenta.push_back((**it).branchingParticle()->momentum());
	mePartonData.push_back((**it).branchingParticle()->dataPtr());
	particles.push_back((**it).branchingParticle());
	}
	}
	// cerr << "testing number of partons\n";
	// for(unsigned int ix=0;ix<meMomenta.size();++ix) {
	// cerr << *particles[ix] << "\n";
	// }
	// boost the momenta to the CMF frame
	// compte boost to reset frame
	Lorentz5Momentum prest(meMomenta[0]+meMomenta[1]);
	LorentzRotation R(-prest.boostVector());
	// and then to put beams along the axis
	Lorentz5Momentum ptest = R*meMomenta[0];
	Axis axis(ptest.vect().unit());
	if(axis.perp2()>0.) {
	R.rotateZ(-axis.phi());
	R.rotateY(-acos(axis.z()));
	}
	const cPDVector partons=diagram->partons();
	// order of the incoming partons
	if(mePartonData[0]!=partons[0]) {
	swap(mePartonData[0],mePartonData[1]);
	swap(meMomenta[0],meMomenta[1]);
	swap(particles[0],particles[1]);
	}
	// order of the outgoing partons
	for(unsigned int ix=2;ix<partons.size();++ix) {
	for(unsigned int iy=ix;iy<meMomenta.size();++iy) {
	if(partons[ix]==mePartonData[iy]) {
	if(ix!=iy) {
	swap(mePartonData[ix],mePartonData[iy]);
	swap(meMomenta[ix],meMomenta[iy]);
	swap(particles[ix],particles[iy]);
	}
	break;
	}
	}
	}
	for(unsigned int ix=0;ix<meMomenta.size();++ix)
	meMomenta[ix].transform(R);
	PPair in(mePartonData[0]->produceParticle(meMomenta[0]),
	mePartonData[1]->produceParticle(meMomenta[1]));
	PVector out;
	for(unsigned int ix=2;ix<meMomenta.size();++ix) {
	out.push_back(mePartonData[ix]->produceParticle(meMomenta[ix]));
	}
	_matrixElement->setKinematics(in,out);
	_matrixElement->dSigHatDR();
	const ColourLines & cl = _matrixElement->selectColourGeometry(diagram);
	PVector slike;
	tPVector ret;
	slike.push_back(particles[0]);
	Ptr<Tree2toNDiagram>::pointer diagram2 =
	dynamic_ptr_cast<Ptr<Tree2toNDiagram>::pointer>(diagram);
	for ( int i = 1; i < diagram2->nSpace() - 1; ++i )
	slike.push_back(diagram2->allPartons()[i]->produceParticle());
	slike.push_back(particles[1]);
	ret = tPVector(slike.begin(), slike.end());
	int io = particles.size();
	PVector tlike(diagram2->allPartons().size() - diagram2->nSpace());
	for ( int i = diagram2->allPartons().size() - 1; i >= diagram2->nSpace(); --i ) {
	int it = i - diagram2->nSpace();
	pair<int,int> ch = diagram2->children(i);
	bool iso = ch.first < 0;
	if ( iso ) {
	tlike[it] = particles[--io];
	}
	else {
	Lorentz5Momentum p = tlike[ch.first - diagram2->nSpace()]->momentum() +
	tlike[ch.second - diagram2->nSpace()]->momentum();
	tlike[it] = diagram2->allPartons()[i]->produceParticle(p);
	}
	}
	ret.insert(ret.end(), tlike.begin(), tlike.end());
	cl.connect(ret);
	for(unsigned int ix=0;ix<ret.size();++ix) {
	PVector::iterator it = find(particles.begin(),particles.end(),ret[ix]);
	if(it==particles.end()) {
	ColinePtr line = ret[ix]->colourLine();
	if(line) line->removeColoured(ret[ix]);
	line = ret[ix]->antiColourLine();
	if(line) line->removeAntiColoured(ret[ix]);
	}
	}
	// now the colours of the rest of the particles
	for(set<HardBranchingPtr>::const_iterator it=tree->branchings().begin();
	it!=tree->branchings().end();++it) (**it).fixColours();
	}

	double PowhegHandler::Sud( Energy scale, long id, Energy pt_cut ){
	//upper limit on scale
	double sudwgt = 1.;
	//scale cut set to just below max_qtilde to avoid seg fault on boundary
	Energy scale_cut = sqrt( _s );
	multimap< long, pair< Interpolator2d< double, Energy, Energy >::Ptr, Energy > >::const_iterator cjt;
	for( cjt = _fbranchings.lower_bound( abs( id ) );
	cjt != _fbranchings.upper_bound( abs( id ) );
	++cjt ) {
	//check to see if we are within qtilde limits before calling interpolator
	if( scale < scale_cut && scale > cjt->second.second )
	sudwgt = ( cjt->second.first )( pt_cut, scale );
	else if ( scale < cjt->second.second )
	sudwgt = ( cjt->second.first )( pt_cut, cjt->second.second);
	else
	sudwgt = ( cjt->second.first )( pt_cut, scale_cut);

	//check for errors in sud value
	if( sudwgt == 0 ){ cerr<<"zero sud at scale = "<< scale / GeV
	<<", id = "<< id
	<<", pt_cut = " << pt_cut / GeV
	<<", scale_cut = " << scale_cut / GeV
	<<"\n";
	}
	if( sudwgt > 1.2 ){ cerr<<"large sud = "<< sudwgt<<" at scale = "<< scale / GeV
	<<", id = "<< id
	<<", pt_cut = " << pt_cut / GeV
	<<", scale_cut = " << scale_cut / GeV
	<<"\n";
	}
	if( isnan( sudwgt ) ){ cerr<<"nan sud at scale = "<< scale / GeV
	<<", id = "<< id
	<< ", pt_cut = " << pt_cut / GeV <<"\n";
	}
	if( isinf( sudwgt ) ){ cerr<<"inf sud at scale = "<< scale / GeV
	<<", id = "<< id
	<< ", pt_cut = " << pt_cut / GeV <<"\n";
	}
	}
	return sudwgt;
	}

	double PowhegHandler::splittingFnWeight( HardTreePtr theTree ){
	double splitFnWgt = 1.;
	for( map< HardBranchingPtr, Energy>::const_iterator cit = theTree->getNodes().begin();
	cit != theTree->getNodes().end(); ++cit ) {
	vector< long > ids;
	ids.push_back( cit->first->branchingParticle()->id() );
	if( ! cit->first->children().empty() ) {
	ids.push_back( cit->first->children()[0]->branchingParticle()->id() );
	ids.push_back( cit->first->children()[1]->branchingParticle()->id() );
	}
	else cerr<< "splittingFnWeight(): node with no children found\n";
	double z = cit->first->children()[0]->z();
	Energy2 t = z * ( 1. - z ) * sqr( cit->second );
	// cerr<<"z = "<< z<< ", qt = "<< cit->second / GeV <<"\n";
	splitFnWgt *= cit->first->sudakov()->splittingFn()->P( z, t, ids, true );
	}
	return splitFnWgt;
	}

	double PowhegHandler::sudakovWeight( HardTreePtr theTree ) {
	double SudWgt = 1.;
	Energy kt_cut;
	if( ! _highestMult ) kt_cut = sqrt( _yini * _s );
	else {
	kt_cut = theTree->lowestPt( _jetMeasureMode );
	if( kt_cut > _max_pt_cut )
	kt_cut = _max_pt_cut;
	}
	//external line weight
	for( map< ShowerParticlePtr, HardBranchingPtr >::const_iterator cit =
	theTree->getExternals().begin();
	cit != theTree->getExternals().end(); ++cit ) {

	Energy scale= sqrt( _s );
	if( cit->second ){
	if( cit->first == cit->second->children()[0]->branchingParticle() )
	scale = cit->second->scale() * cit->second->children()[0]->z();
	else if( cit->first == cit->second->children()[1]->branchingParticle() )
	scale = cit->second->scale() * cit->second->children()[1]->z();
	else cerr<<"could not find child in external HardBranching \n";
	}
	else{
	scale = cit->first->evolutionScale();
	}
	SudWgt *= Sud( scale, cit->first->id(), kt_cut );
	}
	if(SudWgt > 1.1) cerr<<"\n\n\nsudakov from externals > 1!!\n\n\n";
	//intermediate line wgts
	for( map< long, pair< Energy, Energy > >::const_iterator cit
	= theTree->getInternals().begin();
	cit != theTree->getInternals().end(); ++cit ) {

	Energy scale = cit->second.first;

	double internal_wgt = Sud( scale, cit->first, kt_cut );
	scale = cit->second.second;
	internal_wgt /= Sud( scale, cit->first, kt_cut );
	if(internal_wgt > 1.1 \|\| internal_wgt < 0.)cerr<<"\n\nbig internal weight of "<< internal_wgt
	<<"\nnum scale = "
	<<cit->second.first / GeV
	<<"\nden scale = "
	<<cit->second.second /GeV
	<<"\n\n";
	SudWgt *= internal_wgt;
	}
	if(SudWgt > 1.1 ) cerr<<"sud wgt is "<<SudWgt<<"\n";

	double alphaWgt = 1.;
	//alphaS weight
	for( map< HardBranchingPtr, Energy >::const_iterator cit = theTree->getNodes().begin();
	cit != theTree->getNodes().end(); ++cit ) {
	if( ! cit->first->children().empty() )
	alphaWgt *= _alphaS->value( sqr( cit->first->children()[0]->pT() ) ) / _alphaSMG;
	else cerr << "sudakovWeight(): node with no children \n";
	}

	if( SudWgt > 1.1 ) {
	cerr<<"\n\nweight exceeded 1 in PowhegHandler::reweight() !!! \n";
	cerr<<" alpha wgt = "<< alphaWgt
	<<"\n sudWgt = "<< SudWgt<<"\n\n";
	}
	return SudWgt*alphaWgt;
	}

	void PowhegHandler::setBeams(HardTreePtr tree) {
	PPair beams=lastXCombPtr()->lastParticles();
	if((**tree->incoming().begin()).branchingParticle()->momentum().z()/
	beams.first->momentum().z()<0.)
	swap(beams.first,beams.second);
	set<HardBranchingPtr>::iterator it = tree->incoming().begin();
	HardBranchingPtr br=*it;
	br->beam(beams.first);
	while (!br->children().empty()) {
	for(unsigned int ix=0;ix<br->children().size();++ix) {
	if(br->children()[ix]->incoming()) {
	br = br->children()[ix];
	break;
	}
	}
	br->beam(beams.first);
	}
	++it;
	br=*it;
	br->beam(beams.second);
	while (!br->children().empty()) {
	for(unsigned int ix=0;ix<br->children().size();++ix) {
	if(br->children()[ix]->incoming()) {
	br = br->children()[ix];
	break;
	}
	}
	br->beam(beams.second);
	}
	}

	bool PowhegHandler::checkTree(HardTreePtr tree) {
	set<HardBranchingPtr>::const_iterator it;
	bool reject = false;
	for(it=tree->incoming().begin();it!=tree->incoming().end();++it) {
	reject \|=checkBranching(*it);
	}
	for(it=tree->branchings().begin();it!=tree->branchings().end();++it) {
	if((**it).incoming()) continue;
	reject \|=checkBranching(*it);
	}
	return !reject;
	}

	bool PowhegHandler::checkBranching(HardBranchingPtr br) {
	static const double eps(1e-5);
	bool reject(false);
	for(vector<HardBranchingPtr>::const_iterator it=br->children().begin();
	it!=br->children().end();++it) {
	reject \|=checkBranching(*it);
	}
	reject \|= br->z()<-eps \|\| br->z()>1.+eps;
	return reject;
	}

	void PowhegHandler::testSudakovs(){
	ofstream sudTestOut;
	Energy deltaQt = 0.5 * GeV;
	//vector of pts to evaluate at and the colours they should be on plot
	vector< pair<Energy, string> > thePts;
	thePts.push_back( make_pair( 0.*GeV, string("BLACK" ) ) );
	thePts.push_back( make_pair( 2.*GeV, string("RED" ) ) );
	thePts.push_back( make_pair( 5.*GeV, string("BLUE" ) ) );
	thePts.push_back( make_pair( 10*GeV, string("GREEN" ) ) );

	sudTestOut.open( "sudTest.top" );
	multimap< long, pair< Interpolator2d< double, Energy, Energy >::Ptr, Energy > >::const_iterator cjt;
	for( cjt = _fbranchings.begin();
	cjt != _fbranchings.end();
	++cjt ) {
	//loop over pts
	sudTestOut <<"NEW FRAME \nSET WINDOW X 1.6 8 Y 3.5 9\nSET FONT DUPLEX\n"
	<<"TITLE TOP \"Sud Test "<< cjt->first
	<<": BLACK-pt=0GeV, RED-pt=2GeV, BLUE-pt=5GeV, GREEN-pt=10GeV\" \n"
	<<"CASE \"\" \nTITLE LEFT \"Sud(qtilde)\" \nCASE \" \" \n"
	<<"SET ORDER X Y DX \nTITLE BOTTOM \"qtilde / GeV \" \n";
	for( vector< pair< Energy, string > >::const_iterator cit = thePts.begin();
	cit != thePts.end(); ++cit ){
	Energy qtilde = cjt->second.second;
	while( qtilde < _max_qtilde ){
	double sud_val = (*cjt->second.first)( cit->first, qtilde );
	sudTestOut << qtilde / GeV <<"\t"<< sud_val <<"\t"<< deltaQt / 2. / GeV<<"\n";
	qtilde += deltaQt;
	}
	sudTestOut <<"HIST "<< cit->second <<"\n";
	}
	}
	sudTestOut.close();
	}


	HardTreePtr PowhegHandler::doClustering() {
	if(!_lepton) {
	return generalClustering();
	}

	ParticleVector theParts = lastXCombPtr()->subProcess()->outgoing();

	//make an intermediate and add to subprocess if not read in
	if(lastXCombPtr()->subProcess()->intermediates().empty()) {
	return HardTreePtr();
	PPair theIncomings = lastXCombPtr()->subProcess()->incoming();
	//set intermediate to Z
	long intermediate_id = 23;
	PPtr theIntermediate = new_ptr( Particle( getParticleData( intermediate_id ) ) );
	theIntermediate->set5Momentum( theIncomings.first->momentum() +
	theIncomings.second->momentum() );
	//add the intermediate - parent/child relations should be updated
	lastXCombPtr()->subProcess()->addIntermediate( theIntermediate );
	cerr<<"added intermediate\n"
	<< *theIntermediate<<"\n";

	}

	PPtr vb = lastXCombPtr()->subProcess()->intermediates()[0];

	map <ShowerParticlePtr,HardBranchingPtr> theParticles;
	tcPDPtr particle_data;
	ShowerParticlePtr vBoson = new_ptr( ShowerParticle( *vb, 1, false, false ) );
	//loops through the FS particles and create naon branchings
	for( unsigned int i = 0; i < theParts.size(); i++){
	ShowerParticlePtr currentParticle =
	new_ptr( ShowerParticle( *theParts[i], 1, true, false ) );
	theParticles.insert(make_pair(currentParticle,
	new_ptr( HardBranching( currentParticle, SudakovPtr(),
	HardBranchingPtr(),false ) ) ) );
	if(currentParticle->dataPtr()->iColour()==PDT::Colour3\|\|
	currentParticle->dataPtr()->iColour()==PDT::Colour8) {
	ColinePtr newline = new_ptr(ColourLine());
	newline->addColoured(currentParticle);
	}
	if(currentParticle->dataPtr()->iColour()==PDT::Colour3bar\|\|
	currentParticle->dataPtr()->iColour()==PDT::Colour8) {
	ColinePtr newline = new_ptr(ColourLine());
	newline->addAntiColoured(currentParticle);
	}
	}
	// loops clustering until we get down to qqbar
	while( theParticles.size() > 2 ){

	//get number of qqbar pairs
	int qq_pairs = 0;
	for( map<ShowerParticlePtr, HardBranchingPtr>::iterator ita = theParticles.begin();
	ita != theParticles.end() ; ita++ ) {
	if( ita->first->id() > 0 && ita->first->id() < 7 ) qq_pairs++;
	}

	double yij_min = 1.;
	pair< ShowerParticlePtr, ShowerParticlePtr > clusterPair;
	//loops over all pairs of particles in theParticles
	for( map<ShowerParticlePtr, HardBranchingPtr>::iterator ita = theParticles.begin();
	ita != theParticles.end() ; ita++ ) {
	for( map<ShowerParticlePtr, HardBranchingPtr>::iterator itb = theParticles.begin();
	itb != ita; itb++) {
	double yij = getJetMeasure( ita->first, itb->first );
	if( yij < yij_min && splittingAllowed( ita->first, itb->first, qq_pairs ) ) {
	clusterPair.first = ita->first;
	clusterPair.second = itb->first;
	yij_min = yij;
	}
	}
	}
	long thePartId;
	SudakovPtr theSudakov = getSud( thePartId,
	clusterPair.first, clusterPair.second );
	if( !theSudakov ){
	cerr << "can't find the sudakov in: \n";
	cerr << *clusterPair.first<<"\n"
	<< *clusterPair.second<<"\n";
	}
	Lorentz5Momentum pairMomentum = clusterPair.first->momentum() +
	clusterPair.second->momentum();
	pairMomentum.setMass(0.*MeV);
	particle_data = getParticleData( thePartId );

	//creates emitter particle
	ShowerParticlePtr clustered = new_ptr( ShowerParticle( particle_data, true ) );
	clustered->set5Momentum( pairMomentum );

	HardBranchingPtr clusteredBranch( new_ptr( HardBranching( clustered, theSudakov,
	HardBranchingPtr(), false ) ) );
	fixColours( clustered, clusterPair.first, clusterPair.second );
	theParticles.insert( make_pair( clustered, clusteredBranch ) );

	//add children
	clusteredBranch->addChild(theParticles.find(clusterPair.first )->second);
	clusteredBranch->addChild(theParticles.find(clusterPair.second)->second);
	theParticles.erase( clusterPair.first );
	theParticles.erase( clusterPair.second );
	}
	vector<HardBranchingPtr> theBranchings;
	for( map<ShowerParticlePtr, HardBranchingPtr>::iterator it =
	theParticles.begin();
	it != theParticles.end(); ++it )
	theBranchings.push_back( it->second );

	// fix for e+e- to match up the colours of the q qbar pair
	if(theBranchings[0]->branchingParticle()->dataPtr()->iColour()==PDT::Colour3) {
	ColinePtr temp = theBranchings[1]->branchingParticle()->antiColourLine();
	temp->addColoured(theBranchings[0]->branchingParticle());
	theBranchings[0]->branchingParticle()->colourLine()->join(temp);
	}
	else {
	ColinePtr temp = theBranchings[0]->branchingParticle()->antiColourLine();
	temp->addColoured(theBranchings[1]->branchingParticle());
	theBranchings[1]->branchingParticle()->colourLine()->join(temp);
	}
	theBranchings[0]->colourPartner(theBranchings[1]);
	theBranchings[1]->colourPartner(theBranchings[0]);
	vector<HardBranchingPtr> spaceBranchings;
	spaceBranchings.push_back( new_ptr( HardBranching( vBoson, SudakovPtr(),
	HardBranchingPtr(),
	true ) ) );
	theBranchings.push_back( spaceBranchings.back() );
	HardTreePtr powhegtree = new_ptr( HardTree( theBranchings,
	spaceBranchings ) );

	// Calculate the shower variables
	evolver()->showerModel()->kinematicsReconstructor()->
	deconstructDecayJets( powhegtree, evolver() );

	+ //keep track of proportion of trees that are ordered
	+ _trees_created++;
	+ if( powhegtree->checkHardOrdering() )
	+ _ordered_trees_created++;
	+
	+ powhegtree->findNodes();
	+
	return powhegtree;
	}

	diff --git a/Shower/Powheg/Matching/PowhegHandler.h b/Shower/Powheg/Matching/PowhegHandler.h
	--- a/Shower/Powheg/Matching/PowhegHandler.h
	+++ b/Shower/Powheg/Matching/PowhegHandler.h
	@@ -1,483 +1,495 @@
	// -- C++ --
	#ifndef HERWIG_PowhegHandler_H
	#define HERWIG_PowhegHandler_H
	//
	// This is the declaration of the PowhegHandler class.
	//

	#include "Herwig++/Utilities/Histogram.h"
	#include "Herwig++/Shower/ShowerHandler.h"
	#include "Herwig++/Shower/Base/Evolver.h"
	#include "Herwig++/Utilities/Interpolator.h"
	#include "Herwig++/Utilities/Interpolator2d.h"
	#include "Herwig++/Utilities/GaussianIntegrator.h"
	#include "Herwig++/Shower/Base/ShowerParticle.h"
	#include "ThePEG/Config/Pointers.h"
	#include "Herwig++/Shower/Base/HardTree.h"
	#include "Herwig++/Shower/Couplings/ShowerAlpha.h"
	#include "ThePEG/MatrixElement/MEBase.h"
	#include "ThePEG/MatrixElement/DiagramBase.fh"

	namespace Herwig {
	class PowhegHandler;
	class PrototypeBranching;
	class ProtoTree;
	}

	//declaration of thepeg ptr
	namespace ThePEG {

	ThePEG_DECLARE_POINTERS(Herwig::PowhegHandler,PowhegHandlerPtr);
	ThePEG_DECLARE_POINTERS(Herwig::PrototypeBranching,PrototypeBranchingPtr);
	ThePEG_DECLARE_POINTERS(Herwig::ProtoTree,ProtoTreePtr);

	}

	namespace Herwig {

	using namespace ThePEG;

	class PrototypeBranching : public Base {
	public:
	PrototypeBranching(ShowerParticlePtr in) : particle(in) {}
	ShowerParticlePtr particle;
	tSudakovPtr sudakov;
	tPrototypeBranchingPtr parent;
	vector<PrototypeBranchingPtr> children;
	HardBranchingPtr convert();
	PrototypeBranchingPtr reset(PrototypeBranchingPtr,
	map<PrototypeBranchingPtr,PrototypeBranchingPtr> &);
	};

	struct PrototypeTree {
	set<PrototypeBranchingPtr> outgoing;
	set<PrototypeBranchingPtr> incoming;
	set<PrototypeBranchingPtr> currentSpaceLike;
	vector<Energy2> scales;
	HardTreePtr convert();
	map<PrototypeBranchingPtr,PrototypeBranchingPtr> reset();
	DiagPtr diagram;
	};

	class ProtoTree:public Base {
	public:
	ProtoTree() {}
	ProtoTree(const set< HardBranchingPtr > & newBranchings) : theBranchings( newBranchings ) {}

	void addBranching( HardBranchingPtr Branching ){
	theBranchings.insert( Branching );
	}
	bool removeBranching( HardBranchingPtr Branching ){
	if( theBranchings.find( Branching ) != theBranchings.end() ){
	theBranchings.erase( Branching );
	return true;
	}
	else return false;
	}
	const set< HardBranchingPtr > & getBranchings() const {
	return theBranchings;
	}
	private:
	set< HardBranchingPtr > theBranchings;
	};

	/**
	* Here is the documentation of the PowhegHandler class.
	*
	* @see \ref PowhegHandlerInterfaces "The interfaces"
	* defined for PowhegHandler.
	*/
	class PowhegHandler: public ShowerHandler {

	public:

	/**
	* The default constructor.
	*/
	PowhegHandler() : _npoint(10), _sudopt(0), _sudname("sudakov.data"), _jetMeasureMode(1), _lepton(true), _reweightOff(false),
	_highestMult(false), _testSudakovs(false),
	_yini(0.001), _alphaSMG(0.118), _max_qtilde( 91.2GeV ), _max_pt_cut( 45.6GeV ), _min_pt_cut( 0.*GeV ), _clusterOption( 0 ) {}

	/**
	* Perform CKKW reweighting
	*/
	virtual double reweightCKKW(int minMult, int maxMult);

	/**
	* Main method for the cascade
	*/
	virtual void cascade();

	public:

	/**
	* access to the hard tree object
	*/
	inline HardTreePtr getHardTree() {
	return _theHardTree;
	}

	/**
	* access to the merging scale
	*/
	inline Energy getMergeScale() {
	return sqrt( _yini * _s );
	}

	inline bool highestMult(){
	return _highestMult;
	}

	/**
	* access to the jet measure definition being used
	*/
	inline unsigned int jetMeasureMode() {
	return _jetMeasureMode;
	}

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

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

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

	protected:

	/** @name Standard Interfaced functions. */
	//@{

	virtual void dofinish();



	/**
	* Initialize this object. Called in the run phase just before
	* a run begins.
	*/
	virtual void doinitrun();
	//@}

	/** @name Standard Interfaced functions. */
	//@{
	/**
	* Initialize this object after the setup phase before saving an
	* EventGenerator to disk.
	* @throws InitException if object could not be initialized properly.
	*/
	virtual void doinit() throw(InitException);

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

	private:

	//a list of all clusters used in finding all shower histories
	map< HardBranchingPtr , pair< ShowerParticlePtr, ShowerParticlePtr > > _all_clusters;

	//a set of prototrees (which are a set of HardBranchings)
	//can do find on set to see if it contains a certain hardbranching
	set< ProtoTreePtr > _proto_trees;

	//all possible shower configurations that are angular ordered
	vector< pair< HardTreePtr, double > > _hardTrees;

	//just connect up the progenitors in currentProtoTree
	bool simpleColConnections( ProtoTreePtr currentProtoTree );

	bool simpleColConnections( HardTreePtr currentHardTree );


	//recursive method to find all trees
	//does a single clustering on the current tree adding more prototrees if more than one possible branching
	bool fillProtoTrees( map< ShowerParticlePtr, HardBranchingPtr >,
	ProtoTreePtr );

	//looks to see if a cluster of the cluster_particles already exists in _all_clusterings
	//if so returns the pointer to that hardBranching if not creates the hardBranchings, adds it
	//to _all_branchings and returns the pointer
	HardBranchingPtr getCluster( pair< ShowerParticlePtr, ShowerParticlePtr >, map< ShowerParticlePtr, HardBranchingPtr > );

	//checks whether a prototree containing the same branchings exists already in _proto_trees in
	//which case the current tree is a repeat and should be removed (and not recursed)
	bool repeatProtoTree( ProtoTreePtr currentProtoTree );

	/**
	* Clusters the partons and creates a branching history
	* by combining the 2 particles with smallest
	* jet measure out of all allowed pairings until we are left
	* with \f$q\bar{q}\f$.
	*/
	HardTreePtr doClustering( );

	/**
	* Creates all cluster histories and selects one
	* according to its shower probability
	*/
	HardTreePtr doClusteringOrdered( );

	double Sud( Energy scale, long id, Energy pt_cut );


	HardTreePtr generalClustering();
	BranchingElement allowedFinalStateBranching
	(pair<PrototypeBranchingPtr,PrototypeBranchingPtr> &);
	BranchingElement allowedInitialStateBranching
	(pair<PrototypeBranchingPtr,PrototypeBranchingPtr> &);
	DiagPtr getDiagram(const PrototypeTree &);
	Energy2 hadronJetMeasure(const Lorentz5Momentum & p1,
	const Lorentz5Momentum & p2,
	bool final=true);
	void createColourFlow(HardTreePtr,DiagPtr);
	void setBeams(HardTreePtr);
	bool checkTree(HardTreePtr);
	bool checkBranching(HardBranchingPtr);

	/**
	* Calculate the Sudakov weight
	*/
	double sudakovWeight( HardTreePtr );

	/**
	* Calculate the splitting function weight
	*/
	double splittingFnWeight( HardTreePtr );


	/**
	* Checks to see that the splitting is allowed.
	*/
	bool splittingAllowed( ShowerParticlePtr part_i,
	ShowerParticlePtr part_j,
	int qq_pairs);

	/**
	* Sort's out the colour lines
	*/
	void fixColours(tPPtr parent, tPPtr child1, tPPtr child2);

	/**
	* Checks to see that the splitting is allowed and finds the
	* Sudakov form factor for the splitting.
	*/
	SudakovPtr getSud(long & emmitter_id,
	ShowerParticlePtr & part_i,
	ShowerParticlePtr & part_j ) ;

	/**
	* Returns the durham jet measure, yij, for the two particles.
	*/
	double getJetMeasure(ShowerParticlePtr part_i, ShowerParticlePtr part_j);


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

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

	/**
	* Outputs some sudakov test histograms
	*/
	void testSudakovs();

	private:

	/**
	* The chosen hard tree
	*/
	HardTreePtr _theHardTree;

	/**
	* Number of points for the interpolation tables
	*/
	unsigned int _npoint;

	/**
	* Centre of mass energy
	*/
	Energy2 _s;

	/**
	* Map containing the sudakovs for the final-state particles
	*/
	multimap< long, pair< Interpolator2d< double, Energy, Energy >::Ptr, Energy > > _fbranchings;

	/**
	* Pointer to the object calculating the strong coupling
	*/
	ShowerAlphaPtr _alphaS;

	/**
	* Option for the Sudakov table
	*/
	unsigned int _sudopt;

	/**
	* Filename for the Sudakov table
	*/
	string _sudname;

	/**
	* The jet measure algorithm we are using.
	*/
	unsigned int _jetMeasureMode;

	/**
	* Whether \f$e^+e^-\f$ or hadron-hadron
	*/
	bool _lepton;

	/**
	* Whether sudakov reweighting is switched off
	*/
	bool _reweightOff;

	/**
	* Whether we are treating the highest multiplicity contribution
	*/
	bool _highestMult;

	/**
	* Whether the sudakovs should be tested
	*/
	bool _testSudakovs;

	/**
	* The allowed final-state branchings
	*/
	map<pair<long,long>,pair<SudakovPtr,IdList> > _allowedFinal;

	/**
	* The allowed initial-state branchings
	*/
	multimap<long, pair<SudakovPtr,IdList> > _allowedInitial;

	/**
	* The matrix element for the core process
	*/
	MEPtr _matrixElement;

	/**
	* The ckkw merging scale
	*/
	double _yini;

	/**
	* The fixed alphaS value that was used to generate mad graph events
	*/
	double _alphaSMG;

	/**`
	* Histogram of sudakov weights
	*/
	HistogramPtr _hSud;

	/**
	* Histogram of alphaS weights
	*/
	HistogramPtr _halphaS;

	/**
	* maximum qtilde scale for sudakov interpolation tables
	*/
	Energy _max_qtilde;

	/**
	* minimum qtilde scale for sudakov interpolation tables
	*/
	Energy _min_qtilde;

	/**
	* maximum pt cut for sudakov interpolation tables
	*/
	Energy _max_pt_cut;

	/**
	* minimum pt cut for sudakov interpolation tables
	*/
	Energy _min_pt_cut;

	/**
	* which clustering scheme to use
	*/
	unsigned int _clusterOption;

	+ /**
	+ * total number of hard trees created in this run
	+ */
	+ unsigned int _trees_created;
	+
	+ /**
	+ * the number of ordered hard trees
	+ */
	+ unsigned int _ordered_trees_created;
	+
	+
	+
	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

	/** This template specialization informs ThePEG about the name of
	* the PowhegHandler class and the shared object where it is defined. */
	template <>
	struct ClassTraits<Herwig::PowhegHandler>
	: public ClassTraitsBase<Herwig::PowhegHandler> {
	/** Return a platform-independent class name */
	static string className() { return "Herwig::PowhegHandler"; }
	/**
	* The name of a file containing the dynamic library where the class
	* PowhegHandler is implemented. It may also include several, space-separated,
	* libraries if the class PowhegHandler 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 "HwPowhegShower.so"; }
	};

	/** @endcond */

	}

	#endif /* HERWIG_PowhegHandler_H */

File Metadata

Mime Type: text/x-diff
Expires: Sun, Feb 23, 2:37 PM (15 h, 50 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4486662
Default Alt Text: (80 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions