No OneTemporary
Actions

Size

76 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,2095 +1,2096 @@
	// -- 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 << _reweightOpt << _testSudakovs <<_yini
	<< _alphaSMG << _npoint << ounit( _max_qtilde, GeV ) << ounit( _max_pt_cut, GeV )
	<< ounit( _min_pt_cut, GeV ) << _clusterOption << _dalitzOn << _qtildeDist;
	}

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

	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, unsigned int> interfaceReweight
	("ReweightOption",
	"Whether to switch off the sudakov reweighting",
	&PowhegHandler::_reweightOpt, 0, false, false);
	static SwitchOption interfaceReweightOff
	(interfaceReweight,
	"Off",
	"No Sudakov reweighting",
	1);
	static SwitchOption interfaceReweightOn
	(interfaceReweight,
	"On",
	"Do Sudakov reweighting",
	0);
	static SwitchOption interfaceReweightNoSud
	(interfaceReweight,
	"NoSud",
	"alphaS reweighting only",
	2);

	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 allHistories
	(ifaceClusterOption,"allHistories","make all histories and require angular ordering", 0);

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

	static SwitchOption ptChoice
	(ifaceClusterOption,"ptChoice", "choose ordered history with lowest total pt", 2);

	static SwitchOption highProbChoice
	(ifaceClusterOption,"highProbChoice", "choose ordered history with highest probability", 3);

	static Switch<PowhegHandler,bool> interfaceDalitz
	("Dalitz",
	"Switch for dalitz analysis of hard tree clustering (3 jets only)",
	&PowhegHandler::_dalitzOn, false, false, false);
	static SwitchOption interfaceDalitzOff
	(interfaceDalitz,
	"Off",
	"Dalitz analysis off",
	false);
	static SwitchOption interfaceDalitzOn
	(interfaceDalitz,
	"On",
	"Dalitz analysis on",
	true);

	static Switch<PowhegHandler,bool> interfaceQtildeDist
	("QtildeDist",
	"Switch for qtilde distribution hists from sudakov tables",
	&PowhegHandler::_qtildeDist, false, false, false);
	static SwitchOption interfaceQtildeDistOff
	(interfaceQtildeDist,
	"Off",
	"Qtilde analysis off",
	false);
	static SwitchOption interfaceQtildeDistOn
	(interfaceQtildeDist,
	"On",
	"Qtilde analysis on",
	true);
	}

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

	if( _clusterOption == 0 \|\| _clusterOption == 2 \|\| _clusterOption == 3 )
	_theHardTree = doClusteringOrdered();
	else
	_theHardTree = doClustering();

	//if highest mult set veto def
	if( _highestMult == true ) evolver()->setHighest( true );
	else evolver()->setHighest( false );


	// return if fails
	if(!_theHardTree)
	return 0.;
	//call dalitz analysis if asked for
	if( _dalitzOn ) getDalitz();

	// compute the Sudakov weight
	double SudWgt;
	if( _reweightOpt != 1 )
	SudWgt = _lepton ? sudakovWeight( _theHardTree ) : 1.;
	else SudWgt = 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);
	}
	//divide by constant 5 here is just a factor becasue alphas = 0.2 isn't high enough to guarrantee
	//that the alphaS weight is < 1
	if( ! _reweightOpt ){
	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());
	//output ddalitz analysis if switched on
	if( _dalitzOn ){
	string dalitz_fname = generator()->filename() + string("-") + string("dalitz.top");
	ofstream dalitz( dalitz_fname.c_str() );
	dalitz<<"SET WINDOW X 2 9 Y 2 9\n";
	dalitz<<"SET FONT DUPLEX\n";
	dalitz<<"SET LIMITS X 0 1 Y 0 1\n";
	dalitz<<"TITLE BOTTOM \"X011\" \n";
	dalitz<<"CASE \" X X\" \n";
	dalitz<<"TITLE LEFT \"X021\" \n";
	dalitz<<"CASE \" X X\" \n";
	for(int ix = 0; ix < _dalitz_from_q1.size(); ix++ )
	dalitz<< _dalitz_from_q1[ix].first <<"\t"<< _dalitz_from_q1[ix].second <<"\n";
	dalitz << "PLOT RED\n";
	for(int ix = 0; ix < _dalitz_from_q2.size(); ix++ )
	dalitz<< _dalitz_from_q2[ix].first <<"\t"<< _dalitz_from_q2[ix].second <<"\n";
	dalitz << "PLOT BLUE\n";

	dalitz<<"NEW FRAME \n";
	dalitz<<"SET WINDOW X 2 9 Y 2 9\n";
	dalitz<<"SET FONT DUPLEX\n";
	dalitz<<"SET LIMITS X 0 1 Y 0 1\n";
	dalitz<<"TITLE BOTTOM \"X011\" \n";
	dalitz<<"CASE \" X X\" \n";
	dalitz<<"TITLE LEFT \"X021\" \n";
	dalitz<<"CASE \" X X\" \n";
	for(int ix = 0; ix < _dalitz_from_q1.size(); ix++ )
	dalitz<< _dalitz_from_q1[ix].first <<"\t"<< _dalitz_from_q1[ix].second <<"\n";
	dalitz << "PLOT RED\n";

	dalitz<<"NEW FRAME \n";
	dalitz<<"SET WINDOW X 2 9 Y 2 9\n";
	dalitz<<"SET FONT DUPLEX\n";
	dalitz<<"SET LIMITS X 0 1 Y 0 1\n";
	dalitz<<"TITLE BOTTOM \"X011\" \n";
	dalitz<<"CASE \" X X\" \n";
	dalitz<<"TITLE LEFT \"X021\" \n";
	dalitz<<"CASE \" X X\" \n";
	for(int ix = 0; ix < _dalitz_from_q2.size(); ix++ )
	dalitz<< _dalitz_from_q2[ix].first <<"\t"<< _dalitz_from_q2[ix].second <<"\n";
	dalitz << "PLOT BLUE\n";

	cerr<<"no q1 dalitz = "<<_dalitz_from_q1.size()<<"\n";
	cerr<<"no q2 dalitz = "<<_dalitz_from_q2.size()<<"\n";

	dalitz.close();
	}

	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;

	_dalitz_from_q1.clear();
	_dalitz_from_q2.clear();

	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 ) {
	//the pt_cut here is pt in the jet measure variable used
	double currentSud = integrator->value( scale[ jx ], scale[ jx - 1 ], ptCut[ix] );
	sud[ix][jx] = ( sud[ix][ jx - 1 ] + currentSud );
	cerr<<jx<<"\t";
	}
	cerr<<ix<<"\n";
	}
	//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();
	if( _qtildeDist ) makeQtildeDist();
	}

	void PowhegHandler::doinit() {
	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 ) );

	//is this the highest multiplcity channel
	if( theParts.size() == _max_mult ) {
	// cerr<<"highest mult channel \n";
	_highestMult = true;
	}
	else _highestMult = 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();

	//the hardTreePtr that is to be returned
	HardTreePtr chosen_hardTree;

	//choose a hardTree from shower probability
	if( _clusterOption == 0 ){
	long treeIndex;
	do{

	} while ( _hardTrees[ treeIndex ].second / totalWeight < UseRandom::rnd() );
	chosen_hardTree = _hardTrees[ treeIndex ].first;
	}

	//choose hardtree with lowest pt
	else if( _clusterOption == 2 ){
	//set min pt to be large
	Energy min_pt = 9999999999.*GeV;
	for(unsigned int ix = 0; ix < _hardTrees.size(); ix++ ){
	if( _hardTrees[ix].first->totalPt() < min_pt ){
	min_pt = _hardTrees[ix].first->totalPt();
	chosen_hardTree = _hardTrees[ix].first;
	}
	}
	}

	//choose hardtree with highest prob (cluster option = 3)
	else{
	//set min pt to be large
	double max_prob = 0.;
	for(unsigned int ix = 0; ix < _hardTrees.size(); ix++ ){
	if( _hardTrees[ ix ].second > max_prob ){
	max_prob = _hardTrees[ ix ].second;
	chosen_hardTree = _hardTrees[ix].first;
	}
	if( _hardTrees[ ix ].second == max_prob && UseRandom::rndbool() ){
	max_prob = _hardTrees[ ix ].second;
	chosen_hardTree = _hardTrees[ix].first;
	}
	}
	}
	//re-do momentum deconstruction (has been overridden by other trees otherwise)
	simpleColConnections( chosen_hardTree );
	if( ! evolver()->showerModel()->kinematicsReconstructor()
	->deconstructDecayJets( chosen_hardTree, evolver() ) )
	cerr<<"\n\nproblem doing momentum decon in selected tree \n\n";
	if( ! chosen_hardTree ) {
	cerr<<"PowhegHandler::problem in choosing hard tree\n";
	return HardTreePtr();
	}

	return chosen_hardTree;
	}

	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.;
	Energy scale_cut = _max_qtilde;
	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

	//pt_cut called here is in the durham or luclus jet measure

	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 );
	splitFnWgt *= cit->first->sudakov()->splittingFn()->P( z, t, ids, true );
	}
	return splitFnWgt;
	}

	double PowhegHandler::sudakovWeight( HardTreePtr theTree ) {
	double SudWgt = 1.;
	//ktcut for sudakovs
	Energy kt_cut;
	if( ! _highestMult ) kt_cut = sqrt( _yini * _s );
	else {
	//does this lowest pt method return something in the luc/dur jet measure?
	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";
	}
	if( _reweightOpt == 0 )
	return SudWgt*alphaWgt;
	else
	return SudWgt;
	}

	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( sqrt(_s*0.001), string("RED" ) ) );
	thePts.push_back( make_pair( sqrt(_s*0.005 ), string("BLUE" ) ) );
	thePts.push_back( make_pair( sqrt(_s*0.01 ), string("GREEN" ) ) );
	thePts.push_back( make_pair( sqrt(_s*0.05 ), string("CYAN" ) ) );

	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:y_ms=0, RED:y_ms=0.001, BLUE:y_ms=0.005, GREEN:y_ms=0.01, CYAN:y_ms=0.05\" \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];

	//is this the highest multiplcity channel
	if( theParts.size() == _max_mult ) _highestMult = true;
	-
	+ else _highestMult = false;
	+
	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;
	}

	void PowhegHandler::getDalitz(){
	if( _theHardTree->getExternals().size() == 3 ){
	// cerr<<"\n\ndoing dalitz:\n";
	Energy total_energy = 0. * GeV;
	double x1_dal=0.;
	double x2_dal=0.;
	long parent_id = 0;
	for( map< ShowerParticlePtr, HardBranchingPtr >::const_iterator cit =
	_theHardTree->getExternals().begin();
	cit != _theHardTree->getExternals().end(); ++cit ) {
	total_energy += cit->first->momentum().t();

	if( abs( cit->first->id() ) < 7 ){
	if( cit->first->id() > 0 )
	x1_dal = 2. * cit->first->momentum().t() / GeV;
	else x2_dal = 2. * cit->first->momentum().t() / GeV;
	}
	}
	x1_dal /= total_energy / GeV;
	x2_dal /= total_energy / GeV;

	for( set< HardBranchingPtr >::const_iterator cit =
	_theHardTree->branchings().begin();
	cit != _theHardTree->branchings().end(); ++cit ) {
	if( (*cit)->incoming() ) continue;
	if( (*cit)->branchingParticle()->id() > 0 ){
	if( (cit)->children().size() == 2 && (cit)->children()[0]
	&& (*cit)->children()[1] ) {
	parent_id = (*cit)->branchingParticle()->id();
	}
	}
	else {

	if( (cit)->children().size() == 2 && (cit)->children()[0]
	&& (*cit)->children()[0] ) {
	parent_id = (*cit)->branchingParticle()->id();
	}
	}
	}
	if( parent_id > 0 )
	_dalitz_from_q1.push_back( make_pair(x1_dal, x2_dal) );
	else if( parent_id < 0 )
	_dalitz_from_q2.push_back( make_pair(x1_dal, x2_dal) );
	}
	}



	void PowhegHandler::makeQtildeDist(){
	//currently this only makes plots for up quarks

	ofstream sudTestOut;
	Energy deltaQt = 0.5 * GeV;
	//vector of pts to evaluate at and the colours they should be on plot
	vector< Energy > thePts;
	thePts.push_back( sqrt( 0.05 * _s ) );
	thePts.push_back( sqrt( 0.01 * _s ) );
	thePts.push_back( sqrt( 0.005 * _s ) );
	thePts.push_back( sqrt( 0.001 * _s ) );

	//iniialise the histograms
	vector< HistogramPtr > theHists;
	for(int ix = 0; ix < thePts.size(); ix ++ )
	theHists.push_back( new_ptr(Histogram(0.,100.,100) ) );

	sudTestOut.open( "QtildeDist.top" );

	//for each pt value a vector of qtilde and sud_bar(qtilde) can be obtained from 2d interpolator
	//from this can calculate sud(qtilde) and make a 1d interpolator going both ways from the vectors

	//iterator to the 2d interpolator of the up quark
	multimap< long, pair< Interpolator2d< double, Energy, Energy >::Ptr, Energy > >::const_iterator cjt =
	_fbranchings.find( long( 1 ) );

	//fill the vector of qtilde values once and for all
	//how to know exactly what the minimum value is - this came from integrator
	vector< Energy > scale;
	Energy qtildemin = cjt->second.second;
	Energy qtildemax = sqrt( _s );
	for( unsigned int ix = 0; ix < _npoint; ix++ )
	scale.push_back( qtildemin + double( ix ) * ( qtildemax - qtildemin ) / double( _npoint - 1 ) );

	//create a vector of pairs of 1d interpolators
	vector < pair< Interpolator< double, Energy >::Ptr, Interpolator< Energy, double >::Ptr > > sud_interp;

	for( unsigned int jx = 0; jx < thePts.size(); jx++ ){
	double sud_min = (* cjt->second.first )( thePts[jx], qtildemax );

	vector<double> sud;
	for( unsigned int ix = 0; ix < scale.size(); ix ++)
	sud.push_back( sud_min / (* cjt->second.first )( thePts[jx], scale[ix] ) );
	// construct the Interpolators
	Interpolator<double,Energy>::Ptr
	intq = new_ptr(Interpolator<double,Energy>(sud,scale,3));
	Interpolator<Energy,double>::Ptr
	ints = new_ptr(Interpolator<Energy,double>(scale,sud,3));
	sud_interp.push_back( make_pair( intq, ints ) );
	//make the interpolators both ways and add to the vectors
	}
	//fill histograms with the qtilde distributions
	for( unsigned int ix = 0; ix < sud_interp.size(); ix ++ ){
	for(unsigned int jx = 0; jx < 100000000; jx ++ ){
	double sud_wgt = UseRandom::rnd();
	Energy solved_scale = (* sud_interp[ix].second )( sud_wgt );
	(* theHists[ix] ) += solved_scale / GeV;
	}
	}
	//output the histograms
	for( unsigned int ix = 0; ix < theHists.size(); ix ++ ){
	using namespace HistogramOptions;
	theHists[ix]->topdrawOutput(sudTestOut,Frame,
	"RED",
	"qtilde distribution",
	"",
	"",
	"",
	"qtilde / GeV",
	"");
	}
	sudTestOut.close();
	}

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 8:37 PM (1 d, 3 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3806092
Default Alt Text: (76 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions