No OneTemporary
Actions

Size

66 KB

Subscribers

None

View Options

	diff --git a/Hadronization/HadronSelector.cc b/Hadronization/HadronSelector.cc
	--- a/Hadronization/HadronSelector.cc
	+++ b/Hadronization/HadronSelector.cc
	@@ -1,823 +1,823 @@
	// -- C++ --
	//
	// HadronSelector.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the HadronSelector class.
	//

	#include "HadronSelector.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Interface/ParVector.h"
	#include "ThePEG/Interface/RefVector.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include <ThePEG/PDT/EnumParticles.h>
	#include <ThePEG/Repository/EventGenerator.h>
	#include <ThePEG/Repository/CurrentGenerator.h>
	#include <ThePEG/Repository/Repository.h>
	#include "CheckId.h"
	#include <ThePEG/Utilities/DescribeClass.h>

	using namespace Herwig;

	DescribeAbstractClass<HadronSelector,Interfaced>
	describeHadronSelector("Herwig::HadronSelector","");

	namespace {
	// // debug helper
	// void dumpTable(const HadronSelector::HadronTable & tbl) {
	// typedef HadronSelector::HadronTable::const_iterator TableIter;
	// for (TableIter it = tbl.begin(); it != tbl.end(); ++it) {
	// cerr << it->first.first << ' '
	// << it->first.second << '\n';
	// for (HadronSelector::KupcoData::const_iterator jt = it->second.begin();
	// jt != it->second.end(); ++jt) {
	// cerr << '\t' << *jt << '\n';
	// }
	// }
	// }

	- bool weightIsLess (pair<tcPDPtr,double> a, pair<tcPDPtr,double> b) {
	+ bool weightIsLess (pair<long,double> a, pair<long,double> b) {
	return a.second < b.second;
	}


	}

	ostream & Herwig::operator<< (ostream & os,
	const HadronSelector::HadronInfo & hi ) {
	os << std::scientific << std::showpoint
	<< std::setprecision(4)
	<< setw(2)
	<< hi.id << '\t'
	// << hi.ptrData << ' '
	<< hi.swtef << '\t'
	<< hi.wt << '\t'
	<< hi.overallWeight << '\t'
	<< ounit(hi.mass,GeV);
	return os;
	}

	HadronSelector::HadronSelector(unsigned int opt)
	: _pwtDquark( 1.0 ),_pwtUquark( 1.0 ),_pwtSquark( 1.0 ),_pwtCquark( 1.0 ),
	_pwtBquark( 1.0 ),_pwtDIquark( 1.0 ),
	_weight1S0(Nmax,1.),_weight3S1(Nmax,1.),_weight1P1(Nmax,1.),_weight3P0(Nmax,1.),
	_weight3P1(Nmax,1.),_weight3P2(Nmax,1.),_weight1D2(Nmax,1.),_weight3D1(Nmax,1.),
	_weight3D2(Nmax,1.),_weight3D3(Nmax,1.),
	_repwt(Lmax,vector<vector<double> >(Jmax,vector<double>(Nmax))),
	_sngWt( 1.0 ),_decWt( 1.0 ),
	_topt(opt),_trial(0)
	{
	// The mixing angles
	// the ideal mixing angle
	const double idealAngleMix = atan( sqrt(0.5) ) * 180.0 / Constants::pi;
	// \eta-\eta' mixing angle
	_etamix = -23.0;
	// phi-omega mixing angle
	_phimix = +36.0;
	// h_1'-h_1 mixing angle
	_h1mix = idealAngleMix;
	// f_0(1710)-f_0(1370) mixing angle
	_f0mix = idealAngleMix;
	// f_1(1420)-f_1(1285)\f$ mixing angle
	_f1mix = idealAngleMix;
	// f'_2-f_2\f$ mixing angle
	_f2mix = +26.0;
	// eta_2(1870)-eta_2(1645) mixing angle
	_eta2mix = idealAngleMix;
	// phi(???)-omega(1650) mixing angle
	_omhmix = idealAngleMix;
	// phi_3-omega_3 mixing angle
	_ph3mix = +28.0;
	// eta(1475)-eta(1295) mixing angle
	_eta2Smix = idealAngleMix;
	// phi(1680)-omega(1420) mixing angle
	_phi2Smix = idealAngleMix;
	}

	void HadronSelector::persistentOutput(PersistentOStream & os) const {
	os << _partons << _pwtDquark << _pwtUquark << _pwtSquark
	<< _pwtCquark << _pwtBquark << _pwtDIquark
	<< _etamix << _phimix << _h1mix << _f0mix << _f1mix << _f2mix
	<< _eta2mix << _omhmix << _ph3mix << _eta2Smix << _phi2Smix
	<< _weight1S0 << _weight3S1 << _weight1P1 << _weight3P0 << _weight3P1
	<< _weight3P2 << _weight1D2 << _weight3D1 << _weight3D2 << _weight3D3
	<< _forbidden << _sngWt << _decWt << _repwt << _pwt
	<< _table;
	}

	void HadronSelector::persistentInput(PersistentIStream & is, int) {
	is >> _partons >> _pwtDquark >> _pwtUquark >> _pwtSquark
	>> _pwtCquark >> _pwtBquark
	>> _pwtDIquark>> _etamix >> _phimix >> _h1mix >> _f0mix >> _f1mix >> _f2mix
	>> _eta2mix >> _omhmix >> _ph3mix >> _eta2Smix >> _phi2Smix
	>> _weight1S0 >> _weight3S1 >> _weight1P1 >> _weight3P0 >> _weight3P1
	>> _weight3P2 >> _weight1D2 >> _weight3D1 >> _weight3D2 >> _weight3D3
	>> _forbidden >> _sngWt >> _decWt >> _repwt >> _pwt
	>> _table;
	}

	void HadronSelector::Init() {

	static ClassDocumentation<HadronSelector> documentation
	("There is no documentation for the HadronSelector class");

	static Parameter<HadronSelector,double>
	interfacePwtDquark("PwtDquark","Weight for choosing a quark D",
	&HadronSelector::_pwtDquark, 0, 1.0, 0.0, 10.0,
	false,false,false);

	static Parameter<HadronSelector,double>
	interfacePwtUquark("PwtUquark","Weight for choosing a quark U",
	&HadronSelector::_pwtUquark, 0, 1.0, 0.0, 10.0,
	false,false,false);

	static Parameter<HadronSelector,double>
	interfacePwtSquark("PwtSquark","Weight for choosing a quark S",
	&HadronSelector::_pwtSquark, 0, 1.0, 0.0, 10.0,
	false,false,false);

	static Parameter<HadronSelector,double>
	interfacePwtCquark("PwtCquark","Weight for choosing a quark C",
	&HadronSelector::_pwtCquark, 0, 1.0, 0.0, 10.0,
	false,false,false);

	static Parameter<HadronSelector,double>
	interfacePwtBquark("PwtBquark","Weight for choosing a quark B",
	&HadronSelector::_pwtBquark, 0, 1.0, 0.0, 10.0,
	false,false,false);

	static Parameter<HadronSelector,double>
	interfacePwtDIquark("PwtDIquark","Weight for choosing a DIquark",
	&HadronSelector::_pwtDIquark, 0, 1.0, 0.0, 100.0,
	false,false,false);

	static Parameter<HadronSelector,double>
	interfaceSngWt("SngWt","Weight for singlet baryons",
	&HadronSelector::_sngWt, 0, 1.0, 0.0, 10.0,
	false,false,false);

	static Parameter<HadronSelector,double>
	interfaceDecWt("DecWt","Weight for decuplet baryons",
	&HadronSelector::_decWt, 0, 1.0, 0.0, 10.0,
	false,false,false);

	static RefVector<HadronSelector,ParticleData> interfacePartons
	("Partons",
	"The partons which are to be considered as the consistuents of the hadrons.",
	&HadronSelector::_partons, -1, false, false, true, false, false);

	static RefVector<HadronSelector,ParticleData> interfaceForbidden
	("Forbidden",
	"The PDG codes of the particles which cannot be produced in the hadronization.",
	&HadronSelector::_forbidden, -1, false, false, true, false, false);

	//
	// mixing angles
	//
	// the ideal mixing angle
	const double idealAngleMix = atan( sqrt(0.5) ) * 180.0 / Constants::pi;

	static Parameter<HadronSelector,double> interface11S0Mixing
	("11S0Mixing",
	"The mixing angle for the I=0 mesons from the 1 1S0 multiplet,"
	" i.e. eta and etaprime.",
	&HadronSelector::_etamix, -23., -180., 180.,
	false, false, Interface::limited);

	static Parameter<HadronSelector,double> interface13S1Mixing
	("13S1Mixing",
	"The mixing angle for the I=0 mesons from the 1 3S1 multiplet,"
	" i.e. phi and omega.",
	&HadronSelector::_phimix, +36., -180., 180.,
	false, false, Interface::limited);

	static Parameter<HadronSelector,double> interface11P1Mixing
	("11P1Mixing",
	"The mixing angle for the I=0 mesons from the 1 1P1 multiplet,"
	" i.e. h_1' and h_1.",
	&HadronSelector::_h1mix, idealAngleMix, -180., 180.,
	false, false, Interface::limited);

	static Parameter<HadronSelector,double> interface13P0Mixing
	("13P0Mixing",
	"The mixing angle for the I=0 mesons from the 1 3P0 multiplet,"
	" i.e. f_0(1710) and f_0(1370).",
	&HadronSelector::_f0mix, idealAngleMix, -180., 180.,
	false, false, Interface::limited);

	static Parameter<HadronSelector,double> interface13P1Mixing
	("13P1Mixing",
	"The mixing angle for the I=0 mesons from the 1 3P1 multiplet,"
	" i.e. f_1(1420) and f_1(1285).",
	&HadronSelector::_f1mix, idealAngleMix, -180., 180.,
	false, false, Interface::limited);

	static Parameter<HadronSelector,double> interface13P2Mixing
	("13P2Mixing",
	"The mixing angle for the I=0 mesons from the 1 3P2 multiplet,"
	" i.e. f'_2 and f_2.",
	&HadronSelector::_f2mix, 26.0, -180., 180.,
	false, false, Interface::limited);

	static Parameter<HadronSelector,double> interface11D2Mixing
	("11D2Mixing",
	"The mixing angle for the I=0 mesons from the 1 1D2 multiplet,"
	" i.e. eta_2(1870) and eta_2(1645).",
	&HadronSelector::_eta2mix, idealAngleMix, -180., 180.,
	false, false, Interface::limited);

	static Parameter<HadronSelector,double> interface13D0Mixing
	("13D0Mixing",
	"The mixing angle for the I=0 mesons from the 1 3D0 multiplet,"
	" i.e. eta_2(1870) phi(?) and omega(1650).",
	&HadronSelector::_omhmix, idealAngleMix, -180., 180.,
	false, false, Interface::limited);

	static Parameter<HadronSelector,double> interface13D1Mixing
	("13D1Mixing",
	"The mixing angle for the I=0 mesons from the 1 3D1 multiplet,"
	" i.e. phi_3 and omega_3.",
	&HadronSelector::_ph3mix, 28.0, -180., 180.,
	false, false, Interface::limited);

	static Parameter<HadronSelector,double> interface21S0Mixing
	("21S0Mixing",
	"The mixing angle for the I=0 mesons from the 2 1S0 multiplet,"
	" i.e. eta(1475) and eta(1295).",
	&HadronSelector::_eta2Smix, idealAngleMix, -180., 180.,
	false, false, Interface::limited);

	static Parameter<HadronSelector,double> interface23S1Mixing
	("23S1Mixing",
	"The mixing angle for the I=0 mesons from the 1 3S1 multiplet,"
	" i.e. phi(1680) and omega(1420).",
	&HadronSelector::_phi2Smix, idealAngleMix, -180., 180.,
	false, false, Interface::limited);
	//
	// the meson weights
	//
	static ParVector<HadronSelector,double> interface1S0Weights
	("1S0Weights",
	"The weights for the 1S0 multiplets start with n=1.",
	&HadronSelector::_weight1S0, Nmax, 1.0, 0.0, 100.0,
	false, false, Interface::limited);

	static ParVector<HadronSelector,double> interface3S1Weights
	("3S1Weights",
	"The weights for the 3S1 multiplets start with n=1.",
	&HadronSelector::_weight3S1, Nmax, 1.0, 0.0, 100.0,
	false, false, Interface::limited);

	static ParVector<HadronSelector,double> interface1P1Weights
	("1P1Weights",
	"The weights for the 1P1 multiplets start with n=1.",
	&HadronSelector::_weight1P1, Nmax, 1.0, 0.0, 100.0,
	false, false, Interface::limited);

	static ParVector<HadronSelector,double> interface3P0Weights
	("3P0Weights",
	"The weights for the 3P0 multiplets start with n=1.",
	&HadronSelector::_weight3P0, Nmax, 1.0, 0.0, 100.0,
	false, false, Interface::limited);

	static ParVector<HadronSelector,double> interface3P1Weights
	("3P1Weights",
	"The weights for the 3P1 multiplets start with n=1.",
	&HadronSelector::_weight3P1, Nmax, 1.0, 0.0, 100.0,
	false, false, Interface::limited);

	static ParVector<HadronSelector,double> interface3P2Weights
	("3P2Weights",
	"The weights for the 3P2 multiplets start with n=1.",
	&HadronSelector::_weight3P2, Nmax, 1.0, 0.0, 100.0,
	false, false, Interface::limited);

	static ParVector<HadronSelector,double> interface1D2Weights
	("1D2Weights",
	"The weights for the 1D2 multiplets start with n=1.",
	&HadronSelector::_weight1D2, Nmax, 1.0, 0.0, 100.0,
	false, false, Interface::limited);

	static ParVector<HadronSelector,double> interface3D1Weights
	("3D1Weights",
	"The weights for the 3D1 multiplets start with n=1.",
	&HadronSelector::_weight3D1, Nmax, 1.0, 0.0, 100.0,
	false, false, Interface::limited);

	static ParVector<HadronSelector,double> interface3D2Weights
	("3D2Weights",
	"The weights for the 3D2 multiplets start with n=1.",
	&HadronSelector::_weight3D2, Nmax, 1.0, 0.0, 100.0,
	false, false, Interface::limited);

	static ParVector<HadronSelector,double> interface3D3Weights
	("3D3Weights",
	"The weights for the 3D3 multiplets start with n=1.",
	&HadronSelector::_weight3D3, Nmax, 1.0, 0.0, 100.0,
	false, false, Interface::limited);

	static Switch<HadronSelector,unsigned int> interfaceTrial
	("Trial",
	"A Debugging option to only produce certain types of hadrons",
	&HadronSelector::_trial, 0, false, false);
	static SwitchOption interfaceTrialAll
	(interfaceTrial,
	"All",
	"Produce all the hadrons",
	0);
	static SwitchOption interfaceTrialPions
	(interfaceTrial,
	"Pions",
	"Only produce pions",
	1);
	static SwitchOption interfaceTrialSpin2
	(interfaceTrial,
	"Spin2",
	"Only mesons with spin less than or equal to two are produced",
	2);
	static SwitchOption interfaceTrialSpin3
	(interfaceTrial,
	"Spin3",
	"Only hadrons with spin less tan or equal to three are produced",
	3);

	}

	double HadronSelector::mixingStateWeight(long id) {
	switch(id) {
	case ParticleID::eta: return 0.5*probabilityMixing(_etamix ,1);
	case ParticleID::etaprime: return 0.5*probabilityMixing(_etamix ,2);
	case ParticleID::phi: return 0.5*probabilityMixing(_phimix ,1);
	case ParticleID::omega: return 0.5*probabilityMixing(_phimix ,2);
	case ParticleID::hprime_1: return 0.5*probabilityMixing(_h1mix ,1);
	case ParticleID::h_1: return 0.5*probabilityMixing(_h1mix ,2);
	case 10331: return 0.5*probabilityMixing(_f0mix ,1);
	case 10221: return 0.5*probabilityMixing(_f0mix ,2);
	case ParticleID::fprime_1: return 0.5*probabilityMixing(_f1mix ,1);
	case ParticleID::f_1: return 0.5*probabilityMixing(_f1mix ,2);
	case ParticleID::fprime_2: return 0.5*probabilityMixing(_f2mix ,1);
	case ParticleID::f_2: return 0.5*probabilityMixing(_f2mix ,2);
	case 10335: return 0.5*probabilityMixing(_eta2mix ,1);
	case 10225: return 0.5*probabilityMixing(_eta2mix ,2);
	// missing phi member of 13D1 should be here
	case 30223: return 0.5*probabilityMixing(_omhmix ,2);
	case 337: return 0.5*probabilityMixing(_ph3mix ,1);
	case 227: return 0.5*probabilityMixing(_ph3mix ,2);
	case 100331: return 0.5*probabilityMixing(_eta2mix ,1);
	case 100221: return 0.5*probabilityMixing(_eta2mix ,2);
	case 100333: return 0.5*probabilityMixing(_phi2Smix,1);
	case 100223: return 0.5*probabilityMixing(_phi2Smix,2);
	default: return 1./3.;
	}
	}

	void HadronSelector::doinit() {
	Interfaced::doinit();
	// the default partons allowed
	// the quarks
	for ( int ix=1; ix<=5; ++ix ) {
	_partons.push_back(getParticleData(ix));
	}
	// the diquarks
	for(unsigned int ix=1;ix<=5;++ix) {
	for(unsigned int iy=1; iy<=ix;++iy) {
	_partons.push_back(getParticleData(CheckId::makeDiquarkID(ix,iy)));
	}
	}
	// set the weights for the various excited mesons
	// set all to one to start with
	for (int l = 0; l < Lmax; ++l ) {
	for (int j = 0; j < Jmax; ++j) {
	for (int n = 0; n < Nmax; ++n) {
	_repwt[l][j][n] = 1.0;
	}
	}
	}
	// set the others from the relevant vectors
	for( int ix=0;ix<max(int(_weight1S0.size()),int(Nmax));++ix)
	_repwt[0][0][ix]=_weight1S0[ix];
	for( int ix=0;ix<max(int(_weight3S1.size()),int(Nmax));++ix)
	_repwt[0][1][ix]=_weight3S1[ix];
	for( int ix=0;ix<max(int(_weight1P1.size()),int(Nmax));++ix)
	_repwt[1][1][ix]=_weight1P1[ix];
	for( int ix=0;ix<max(int(_weight3P0.size()),int(Nmax));++ix)
	_repwt[1][0][ix]=_weight3P0[ix];
	for( int ix=0;ix<max(int(_weight3P1.size()),int(Nmax));++ix)
	_repwt[1][1][ix]=_weight3P1[ix];
	for( int ix=0;ix<max(int(_weight3P2.size()),int(Nmax));++ix)
	_repwt[1][2][ix]=_weight3P2[ix];
	for( int ix=0;ix<max(int(_weight1D2.size()),int(Nmax));++ix)
	_repwt[2][2][ix]=_weight1D2[ix];
	for( int ix=0;ix<max(int(_weight3D1.size()),int(Nmax));++ix)
	_repwt[2][1][ix]=_weight3D1[ix];
	for( int ix=0;ix<max(int(_weight3D2.size()),int(Nmax));++ix)
	_repwt[2][2][ix]=_weight3D2[ix];
	for( int ix=0;ix<max(int(_weight3D3.size()),int(Nmax));++ix)
	_repwt[2][3][ix]=_weight3D3[ix];
	// weights for the different quarks etc
	for(unsigned int ix=0; ix<_partons.size(); ++ix) {
	- _pwt[_partons[ix]]=1.;
	+ _pwt[_partons[ix]->id()]=1.;
	}
	- _pwt[getParticleData(1)] = _pwtDquark;
	- _pwt[getParticleData(2)] = _pwtUquark;
	- _pwt[getParticleData(3)] = _pwtSquark;
	- _pwt[getParticleData(4)] = _pwtCquark;
	- _pwt[getParticleData(5)] = _pwtBquark;
	- _pwt[getParticleData(1103)] = _pwtDIquark * _pwtDquark * _pwtDquark;
	- _pwt[getParticleData(2101)] = 0.5 * _pwtDIquark * _pwtUquark * _pwtDquark;
	- _pwt[getParticleData(2203)] = _pwtDIquark * _pwtUquark * _pwtUquark;
	- _pwt[getParticleData(3101)] = 0.5 * _pwtDIquark * _pwtSquark * _pwtDquark;
	- _pwt[getParticleData(3201)] = 0.5 * _pwtDIquark * _pwtSquark * _pwtUquark;
	- _pwt[getParticleData(3303)] = _pwtDIquark * _pwtSquark * _pwtSquark;
	+ _pwt[1] = _pwtDquark;
	+ _pwt[2] = _pwtUquark;
	+ _pwt[3] = _pwtSquark;
	+ _pwt[4] = _pwtCquark;
	+ _pwt[5] = _pwtBquark;
	+ _pwt[1103] = _pwtDIquark * _pwtDquark * _pwtDquark;
	+ _pwt[2101] = 0.5 * _pwtDIquark * _pwtUquark * _pwtDquark;
	+ _pwt[2203] = _pwtDIquark * _pwtUquark * _pwtUquark;
	+ _pwt[3101] = 0.5 * _pwtDIquark * _pwtSquark * _pwtDquark;
	+ _pwt[3201] = 0.5 * _pwtDIquark * _pwtSquark * _pwtUquark;
	+ _pwt[3303] = _pwtDIquark * _pwtSquark * _pwtSquark;
	// Commenting out heavy di-quark weights
	- _pwt[getParticleData(4101)] = 0.0;
	- _pwt[getParticleData(4201)] = 0.0;
	- _pwt[getParticleData(4301)] = 0.0;
	- _pwt[getParticleData(4403)] = 0.0;
	- _pwt[getParticleData(5101)] = 0.0;
	- _pwt[getParticleData(5201)] = 0.0;
	- _pwt[getParticleData(5301)] = 0.0;
	- _pwt[getParticleData(5401)] = 0.0;
	- _pwt[getParticleData(5503)] = 0.0;
	+ _pwt[4101] = 0.0;
	+ _pwt[4201] = 0.0;
	+ _pwt[4301] = 0.0;
	+ _pwt[4403] = 0.0;
	+ _pwt[5101] = 0.0;
	+ _pwt[5201] = 0.0;
	+ _pwt[5301] = 0.0;
	+ _pwt[5401] = 0.0;
	+ _pwt[5503] = 0.0;
	// find the maximum
	- map<tcPDPtr,double>::iterator pit =
	+ map<long,double>::iterator pit =
	max_element(_pwt.begin(),_pwt.end(),weightIsLess);
	- double pmax = pit->second;
	+ const double pmax = pit->second;
	for(pit=_pwt.begin(); pit!=_pwt.end(); ++pit) {
	pit->second/=pmax;
	}
	// construct the hadron tables
	constructHadronTable();

	// for debugging
	// dumpTable(table());
	}

	void HadronSelector::constructHadronTable() {
	// initialise the table
	_table.clear();
	for(unsigned int ix=0; ix<_partons.size(); ++ix) {
	for(unsigned int iy=0; iy<_partons.size(); ++iy) {
	if (!(DiquarkMatcher::Check(_partons[ix]->id())
	&& DiquarkMatcher::Check(_partons[iy]->id())))
	_table[make_pair(_partons[ix]->id(),_partons[iy]->id())] = KupcoData();
	}
	}
	// get the particles from the event generator
	ParticleMap particles = generator()->particles();
	// loop over the particles
	double maxdd(0.),maxss(0.),maxrest(0.);
	for(ParticleMap::iterator it=particles.begin();
	it!=particles.end(); ++it) {
	long pid = it->first;
	tPDPtr particle = it->second;
	int pspin = particle->iSpin();
	// Don't include hadrons which are explicitly forbidden
	if(find(_forbidden.begin(),_forbidden.end(),particle)!=_forbidden.end())
	continue;
	// Don't include non-hadrons or antiparticles
	if(pid < 100) continue;
	// remove diffractive particles
	if(pspin == 0) continue;
	// K_0S and K_0L not made make K0 and Kbar0
	if(pid==ParticleID::K_S0\|\|pid==ParticleID::K_L0) continue;
	// Debugging options
	// Only include those with 2J+1 less than...5
	if(_trial==2 && pspin >= 5) continue;
	// Only include those with 2J+1 less than...7
	if(_trial==3 && pspin >= 7) continue;
	// Only include pions
	if(_trial==1 && pid!=111 && pid!=211) continue;
	// shouldn't be coloured
	if(particle->coloured()) continue;
	// Get the flavours
	const int x4 = (pid/1000)%10;
	const int x3 = (pid/100 )%10;
	const int x2 = (pid/10 )%10;
	const int x7 = (pid/1000000)%10;
	const bool wantSusy = x7 == 1 \|\| x7 == 2;
	int flav1;
	int flav2;
	// Skip non-hadrons (susy particles, etc...)
	if(x3 == 0 \|\| x2 == 0) continue;
	else if(x4 == 0) { // meson
	flav1 = x2;
	flav2 = x3;
	}
	else { // baryon
	flav1 = CheckId::makeDiquarkID(x2,x3);
	flav2 = x4;
	}
	if (wantSusy) flav2 += 1000000 * x7;
	HadronInfo a(pid,
	particle,
	specialWeight(pid),
	particle->mass());
	// set the weight to the number of spin states
	a.overallWeight = pspin;
	// identical light flavours
	if(flav1 == flav2 && flav1<=3) {
	// ddbar> uubar> admixture states
	if(flav1==1) {
	if(_topt != 0) a.overallWeight = 0.5a.swtef;
	_table[make_pair(1,1)].insert(a);
	_table[make_pair(2,2)].insert(a);
	if(_topt == 0 && a.overallWeight > maxdd) maxdd = a.overallWeight;
	}
	// load up ssbar> uubar> ddbar> admixture states
	else {
	a.wt = mixingStateWeight(pid);
	a.overallWeight *= a.wt;
	if(_topt != 0) a.overallWeight *= a.swtef;
	_table[make_pair(1,1)].insert(a);
	_table[make_pair(2,2)].insert(a);
	if(_topt == 0 && a.overallWeight > maxdd) maxdd = a.overallWeight;
	a.wt = (_topt != 0) ? 1.- 2.*a.wt : 1 - a.wt;
	if(a.wt > 0) {
	a.overallWeight = a.wt * a.swtef * pspin;
	_table[make_pair(3,3)].insert(a);
	if(_topt == 0 && a.overallWeight > maxss) maxss = a.overallWeight;
	}
	}
	}
	// light baryons with all quarks identical
	else if((flav1 == 1 && flav2 == 1103) \|\| (flav1 == 1103 && flav2 == 1) \|\|
	(flav1 == 2 && flav2 == 2203) \|\| (flav1 == 2203 && flav2 == 2) \|\|
	(flav1 == 3 && flav2 == 3303) \|\| (flav1 == 3303 && flav2 == 3)) {
	if(_topt != 0) a.overallWeight = 1.5a.swtef;
	_table[make_pair(flav1,flav2)].insert(a);
	_table[make_pair(flav2,flav1)].insert(a);
	if(_topt == 0 && a.overallWeight > maxrest) maxrest = a.overallWeight;
	}
	// all other cases
	else {
	if(_topt != 0) a.overallWeight *=a.swtef;
	_table[make_pair(flav1,flav2)].insert(a);
	if(flav1 != flav2) _table[make_pair(flav2,flav1)].insert(a);
	if(_topt == 0 && a.overallWeight > maxrest) maxrest = a.overallWeight;
	}
	}

	// Account for identical combos of diquark/quarks and symmetrical elements
	// e.g. U UD = D UU
	HadronTable::iterator tit;
	for(tit=_table.begin();tit!=_table.end();++tit) {
	if(tit->first.first>ParticleID::c) continue;
	if(!DiquarkMatcher::Check(tit->first.second)) continue;
	long k, l, sub;
	if(tit->first.second>=ParticleID::bd_0) {
	k = ParticleID::b;
	sub = ParticleID::bd_0/100;
	}
	else if(tit->first.second>=ParticleID::cd_0) {
	k = ParticleID::c;
	sub = ParticleID::cd_0/100;
	}
	else if(tit->first.second>=ParticleID::sd_0) {
	k = ParticleID::s;
	sub = ParticleID::sd_0/100;
	}
	else if(tit->first.second>=ParticleID::ud_0) {
	k = ParticleID::u;
	sub = ParticleID::ud_0/100;
	}
	else if(tit->first.second==ParticleID::dd_1) {
	k = ParticleID::d;
	sub = ParticleID::dd_1/100;
	}
	else continue;
	sub=tit->first.second/100-sub+1;
	if(sub > tit->first.first) {
	l = 1000sub+100tit->first.first+1;
	}
	else if(sub==tit->first.first) {
	l = 1000sub+ 100tit->first.first+3;
	}
	else {
	l = 100sub +1000tit->first.first+1;
	}
	if(tit->second.empty()) {
	pair<long,long> newpair(k,l);
	tit->second=_table[newpair];
	newpair=make_pair(tit->first.second,tit->first.first);
	_table[newpair]=tit->second;
	};
	}

	// normalise weights to one for first option
	if(_topt == 0) {
	HadronTable::const_iterator tit;
	KupcoData::iterator it;
	for(tit=_table.begin();tit!=_table.end();++tit) {
	double weight;
	if(tit->first.first==tit->first.second) {
	if(tit->first.first==1\|\|tit->first.first==2) weight=1./maxdd;
	else if (tit->first.first==3) weight=1./maxss;
	else weight=1./maxrest;
	}
	else weight=1./maxrest;
	for(it = tit->second.begin(); it!=tit->second.end(); ++it) {
	it->rescale(weight);
	}
	}
	}
	}

	double HadronSelector::specialWeight(long id) {
	const int pspin = id % 10;
	// Only K0L and K0S have pspin == 0, should
	// not get them until Decay step
	assert( pspin != 0 );
	// Baryon : J = 1/2 or 3/2
	if(pspin == 2) {
	// Singlet (Lambda-like) baryon
	if( (id/100)%10 < (id/10 )%10 ) return sqr(_sngWt);
	// octet
	else return 1.;
	}
	// Decuplet baryon
	else if (pspin == 4) {
	return sqr(_decWt);
	}
	// Meson
	else if(pspin % 2 == 1) {
	// Total angular momentum
	int j = (pspin - 1) / 2;
	// related to Orbital angular momentum l
	int nl = (id/10000 )%10;
	int l = -999;
	int n = (id/100000)%10; // Radial excitation
	if(j == 0) l = nl;
	else if(nl == 0) l = j - 1;
	else if(nl == 1 \|\| nl == 2) l = j;
	else if(nl == 3) l = j + 1;
	// Angular or Radial excited meson
	if((l\|\|j\|\|n) && l>=0 && l<Lmax && j<Jmax && n<Nmax) {
	return sqr(_repwt[l][j][n]);
	}
	}
	// rest is not excited or
	// has spin >= 5/2 (ispin >= 6), haven't got those
	return 1.0;
	}

	int HadronSelector::signHadron(tcPDPtr idQ1, tcPDPtr idQ2,
	tcPDPtr hadron) const {
	// This method receives in input three PDG ids, whose the
	// first two have proper signs (corresponding to particles, id > 0,
	// or antiparticles, id < 0 ), whereas the third one must
	// be always positive (particle not antiparticle),
	// corresponding to:
	// --- quark-antiquark, or antiquark-quark, or
	// quark-diquark, or diquark-quark, or
	// antiquark-antidiquark, or antidiquark-antiquark
	// for the first two input (idQ1, idQ2);
	// --- meson or baryon for the third input (idHad):
	// The method returns:
	// --- + 1 if the two partons (idQ1, idQ2) are exactly
	// the constituents for the hadron idHad;
	// --- - 1 if the two partons (idQ1, idQ2) are exactly
	// the constituents for the anti-hadron -idHad;
	// --- + 0 otherwise.
	// The method it is therefore useful to decide the
	// sign of the id of the produced hadron as appeared
	// in the vector _vecHad (where only hadron idHad > 0 are present)
	// given the two constituent partons.
	int sign = 0;
	long idHad = hadron->id();
	assert(idHad > 0);
	int chargeIn = idQ1->iCharge() + idQ2->iCharge();
	int chargeOut = hadron->iCharge();
	// same charge
	if( chargeIn == chargeOut && chargeIn !=0 ) sign = +1;
	else if(chargeIn == -chargeOut && chargeIn !=0 ) sign = -1;
	else if(chargeIn == 0 && chargeOut == 0 ) {
	// In the case of same null charge, there are four cases:
	// i) K0-like mesons, B0-like mesons, Bs-like mesons
	// the PDG convention is to consider them "antiparticle" (idHad < 0)
	// if the "dominant" (heavier) flavour (respectively, s, b)
	// is a quark (idQ > 0): for instance, B0s = (b, sbar) has id < 0
	// Remember that there is an important exception for K0L (id=130) and
	// K0S (id=310): they don't have antiparticles, therefore idHad > 0
	// always. We use below the fact that K0L and K0S are the unique
	// hadrons having 0 the first (less significant) digit of their id.
	// 2) D0-like mesons: the PDG convention is to consider them "particle"
	// (idHad > 0) if the charm flavour is carried by a c: (c,ubar) has id>0
	// 3) the remaining mesons should not have antiparticle, therefore their
	// sign is always positive.
	// 4) for baryons, that is when one of idQ1 and idQ2 is a (anti-) quark and
	// the other one is a (anti-) diquark the sign is negative when both
	// constituents are "anti", that is both with id < 0; positive otherwise.
	// meson
	if(abs(int(idQ1->iColour()))== 3 && abs(int(idQ2->iColour())) == 3 &&
	!DiquarkMatcher::Check(idQ1->id()) && !DiquarkMatcher::Check(idQ2->id()))
	{
	int idQa = abs(idQ1->id());
	int idQb = abs(idQ2->id());
	int dominant = idQ2->id();

	if(idQa > idQb) {
	swap(idQa,idQb);
	dominant = idQ1->id();
	}

	if((idQa==ParticleID::d && idQb==ParticleID::s) \|\|
	(idQa==ParticleID::d && idQb==ParticleID::b) \|\|
	(idQa==ParticleID::s && idQb==ParticleID::b)) {
	// idHad%10 is zero for K0L,K0S
	if (dominant < 0 \|\| idHad%10 == 0) sign = +1;
	else if(dominant > 0) sign = -1;
	}
	else if((idQa==ParticleID::u && idQb==ParticleID::c) \|\|
	(idQa==ParticleID::u && idQb==ParticleID::t) \|\|
	(idQa==ParticleID::c && idQb==ParticleID::t)) {
	if (dominant > 0) sign = +1;
	else if(dominant < 0) sign = -1;
	}
	else if(idQa==idQb) sign = +1;
	// sets sign for Susy particles
	else sign = (dominant > 0) ? +1 : -1;
	}
	// baryon
	else if(DiquarkMatcher::Check(idQ1->id()) \|\| DiquarkMatcher::Check(idQ2->id())) {
	if (idQ1->id() > 0 && idQ2->id() > 0) sign = +1;
	else if(idQ1->id() < 0 && idQ2->id() < 0) sign = -1;
	}
	}
	if (sign == 0) {
	cerr << "Could not work out sign for "
	<< idQ1->PDGName() << ' '
	<< idQ2->PDGName() << " => "
	<< hadron->PDGName() << '\n';
	assert(false);
	}
	return sign;
	}

	pair<tcPDPtr,tcPDPtr> HadronSelector::lightestHadronPair(tcPDPtr ptr1, tcPDPtr ptr2,
	tcPDPtr ptr3) const {
	// throw exception of id3!=0 as doesn't work
	if ( ptr3 ) throw Exception()
	<< "ptr3!=0 not yet implemented in HadronSelector::lightestHadronPair"
	<< Exception::abortnow;

	// charge
	int totalcharge = ptr1->iCharge() + ptr2->iCharge();
	if ( ptr3 ) totalcharge += ptr3->iCharge();

	tcPDPtr vIdHad1[2]={tcPDPtr(),tcPDPtr()},vIdHad2[2]={tcPDPtr(),tcPDPtr()};
	bool vOk[2] = {false, false};
	Energy vMassPair[2] = { ZERO, ZERO };
	for (int i = 0; i < 2; i++) {
	tcPDPtr idPartner = i==0 ? getParticleData(ParticleID::d) : getParticleData(ParticleID::u);
	// Change sign to idPartner (transform it into a anti-quark) if it is not
	// possible to form a meson or a baryon.
	assert (ptr1 && idPartner);
	if (!CheckId::canBeHadron(ptr1, idPartner)) idPartner = idPartner->CC();

	vIdHad1[i] = lightestHadron(ptr1, idPartner);
	vIdHad2[i] = lightestHadron(ptr2, idPartner->CC());
	if ( vIdHad1[i] && vIdHad2[i] &&
	vIdHad1[i]->iCharge() + vIdHad2[i]->iCharge() == totalcharge ) {
	vOk[i] = true;
	vMassPair[i] = vIdHad1[i]->mass() + vIdHad2[i]->mass();
	}
	}
	// Take the lightest pair compatible with charge conservation.
	if ( vOk[0] && ( ! vOk[1] \|\| vMassPair[0] <= vMassPair[1] ) ) {
	return make_pair(vIdHad1[0],vIdHad2[0]);
	}
	else if ( vOk[1] && ( ! vOk[0] \|\| vMassPair[1] < vMassPair[0] ) ) {
	return make_pair(vIdHad1[1],vIdHad2[1]);
	}
	else {
	return make_pair(tcPDPtr(),tcPDPtr());
	}
	}

	Energy HadronSelector::massLightestBaryonPair(tcPDPtr ptr1, tcPDPtr ptr2) const {
	// Make sure that we don't have any diquarks as input, return arbitrarily
	// large value if we do
	Energy currentSum = Constants::MaxEnergy;
	for(unsigned int ix=0; ix<_partons.size(); ++ix) {
	if(!DiquarkMatcher::Check(_partons[ix]->id())) continue;
	HadronTable::const_iterator
	tit1=_table.find(make_pair(abs(ptr1->id()),_partons[ix]->id())),
	tit2=_table.find(make_pair(_partons[ix]->id(),abs(ptr2->id())));
	if( tit1==_table.end() \|\| tit2==_table.end()) continue;
	if(tit1->second.empty()\|\|tit2->second.empty()) continue;
	Energy s = tit1->second.begin()->mass + tit2->second.begin()->mass;
	if(currentSum > s) currentSum = s;
	}
	return currentSum;
	}

	diff --git a/Hadronization/HadronSelector.h b/Hadronization/HadronSelector.h
	--- a/Hadronization/HadronSelector.h
	+++ b/Hadronization/HadronSelector.h
	@@ -1,813 +1,811 @@
	// -- C++ --
	//
	// HadronSelector.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef HERWIG_HadronSelector_H
	#define HERWIG_HadronSelector_H
	//
	// This is the declaration of the HadronSelector class.
	//

	#include "ThePEG/Interface/Interfaced.h"
	#include <ThePEG/Persistency/PersistentOStream.h>
	#include <ThePEG/Persistency/PersistentIStream.h>
	#include <ThePEG/PDT/ParticleData.h>
	#include <ThePEG/PDT/StandardMatchers.h>
	#include <ThePEG/Repository/EventGenerator.h>
	#include "HadronSelector.fh"

	namespace Herwig {

	using namespace ThePEG;

	/**\ingroup Hadronization
	* \class HadronSelector
	* \brief This class selects the hadron flavours of a cluster decay.
	* \author Philip Stephens
	* \author Alberto Ribon
	* \author Peter Richardson
	*
	* This is the base class for the selection of either a pair of hadrons, or
	* in some cases a single hadron. The different approaches which were
	* previously implemented in this class are now implemented in the
	* HwppSelector and Hw64Selector which inherit from this class.
	*
	* This class implements a number of methods which are needed by all models
	* and in addition contains the weights for the different meson multiplets and
	* mixing of the light \f$I=0\f$ mesons.
	*
	* @see \ref HadronSelectorInterfaces "The interfaces"
	* defined for HadronSelector.
	* @see HwppSelector
	* @see Hw64Selector
	*/
	class HadronSelector: public Interfaced {

	public:

	/** \ingroup Hadronization
	* \class HadronInfo
	* \brief Class used to store all the hadron information for easy access.
	* \author Philip Stephens
	*
	* Note that:
	* - the hadrons in _table can be filled in any ordered
	* w.r.t. the mass value, and flavours for different
	* groups (for instance, (u,s) hadrons don't need to
	* be placed after (d,s) or any other flavour), but
	* all hadrons with the same flavours must be consecutive
	* ( for instance you cannot alternate hadrons of type
	* (d,s) with those of flavour (u,s) ).
	* Furthermore, it is assumed that particle and antiparticle
	* have the same weights, and therefore only one of them
	* must be entered in the table: we have chosen to refer
	* to the particle, defined as PDG id > 0, although if
	* an anti-particle is provided in input it is automatically
	* transform to its particle, simply by taking the modulus
	* of its id.
	*/
	class HadronInfo {

	public:

	/**
	* Constructor
	* @param idin The PDG code of the hadron
	* @param datain The pointer to the ParticleData object
	* @param swtin The singlet/decuplet/orbital factor
	* @param massin The mass of the hadron
	*/
	HadronInfo(long idin=0, tPDPtr datain=tPDPtr(),
	double swtin=1., Energy massin=ZERO)
	: id(idin), ptrData(datain), swtef(swtin), wt(1.0), overallWeight(0.0),
	mass(massin)
	{}

	/**
	* Comparision operator on mass
	*/
	bool operator<(const HadronInfo &x) const {
	if(mass!=x.mass) return mass < x.mass;
	else return id < x.id;
	}

	/**
	* The hadrons id.
	*/
	long id;

	/**
	* pointer to ParticleData, to get the spin, etc...
	*/
	tPDPtr ptrData;

	/**
	* singlet/decuplet/orbital factor
	*/
	double swtef;

	/**
	* mixing factor
	*/
	double wt;

	/**
	* (2J+1)wt*swtef
	*/
	double overallWeight;

	/**
	* The hadrons mass
	*/
	Energy mass;

	/**
	* Rescale the weight for a given hadron
	*/
	void rescale(double x) const {
	const_cast<HadronInfo>(this)->overallWeight = x;
	}

	/**
	* Friend method used to print the value of a table element.
	*/
	friend PersistentOStream & operator<< (PersistentOStream & os,
	const HadronInfo & hi ) {
	os << hi.id << hi.ptrData << hi.swtef << hi.wt
	<< hi.overallWeight << ounit(hi.mass,GeV);
	return os;
	}

	/**
	* debug output
	*/
	friend ostream & operator<< (ostream & os, const HadronInfo & hi );

	/**
	* Friend method used to read in the value of a table element.
	*/
	friend PersistentIStream & operator>> (PersistentIStream & is,
	HadronInfo & hi ) {
	is >> hi.id >> hi.ptrData >> hi.swtef >> hi.wt
	>> hi.overallWeight >> iunit(hi.mass,GeV);
	return is;
	}
	};

	/** \ingroup Hadronization
	* \class Kupco
	* \brief Class designed to make STL routines easy to use.
	* \author Philip Stephens
	*
	* This class is used to generate a list of the hadron pairs which can
	* be produced that allows easy traversal and quick access.
	*/
	class Kupco {

	public:

	/**
	* Constructor
	* @param inidQ PDG code of the quark drawn from the vacuum.
	* @param inhad1 ParticleData for the first hadron produced.
	* @param inhad2 ParticleData for the second hadron produced.
	* @param inwgt The weight for the hadron pair
	*/
	Kupco(tcPDPtr inidQ,tcPDPtr inhad1,tcPDPtr inhad2, Energy inwgt)
	: idQ(inidQ),hadron1(inhad1),hadron2(inhad2),weight(inwgt)
	{}

	/**
	* id of the quark drawn from the vacuum.
	*/
	tcPDPtr idQ;

	/**
	* The ParticleData object for the first hadron produced.
	*/
	tcPDPtr hadron1;

	/**
	* The ParticleData object for the second hadron produced.
	*/
	tcPDPtr hadron2;

	/**
	* Weight factor of this componation.
	*/
	Energy weight;
	};

	public:

	/**
	* The helper classes
	*/
	//@{
	/**
	* The type is used to contain all the hadrons info of a given flavour.
	*/
	typedef set<HadronInfo> KupcoData;
	//@}

	/**
	* The hadron table type.
	*/
	typedef map<pair<long,long>,KupcoData> HadronTable;

	public:
	/**
	* The default constructor.
	*/
	HadronSelector(unsigned int);

	/**
	* Method to return a pair of hadrons given the PDG codes of
	* two or three constituents
	* @param cluMass The mass of the cluster
	* @param par1 The first constituent
	* @param par2 The second constituent
	* @param par3 The third constituent
	*/
	virtual pair<tcPDPtr,tcPDPtr> chooseHadronPair(const Energy cluMass, tcPDPtr par1,
	tcPDPtr par2,tcPDPtr par3 = PDPtr())
	= 0;

	/**
	* This returns the lightest pair of hadrons given by the flavours.
	*
	* Given the two (or three) constituents of a cluster, it returns
	* the two lightest hadrons with proper flavour numbers.
	* Furthermore, the first of the two hadrons must have the constituent with
	* par1, and the second must have the constituent with par2.
	* \todo At the moment it does nothing in the case that also par3 is present.
	*
	* The method is implemented by calling twice lightestHadron,
	* once with (par1,quarktopick->CC()) ,and once with (par2,quarktopick)
	* where quarktopick is either the pointer to
	* d or u quarks . In fact, the idea is that whatever the flavour of par1
	* and par2, no matter if (anti-)quark or (anti-)diquark, the lightest
	* pair of hadrons containing flavour par1 and par2 will have either
	* flavour d or u, being the lightest quarks.
	* The method returns the pair (PDPtr(),PDPtr()) if anything goes wrong.
	*
	* \todo The method assumes par3 == PDPtr() (otherwise we don't know how to proceed: a
	* possible, trivial way would be to randomly select two of the three
	* (anti-)quarks and treat them as a (anti-)diquark, reducing the problem
	* to two components as treated below.
	* In the normal (two components) situation, the strategy is the following:
	* treat in the same way the two possibilities: (d dbar) (i=0) and
	* (u ubar) (i=1) as the pair quark-antiquark necessary to form a
	* pair of hadrons containing the input flavour par1 and par2; finally,
	* select the one that produces the lightest pair of hadrons, compatible
	* with the charge conservation constraint.
	*/
	pair<tcPDPtr,tcPDPtr> lightestHadronPair(tcPDPtr ptr1, tcPDPtr ptr2,
	tcPDPtr ptr3 = PDPtr ()) const;

	/**
	* Returns the mass of the lightest pair of hadrons with the given particles
	* @param ptr1 is the first constituent
	* @param ptr2 is the second constituent
	* @param ptr3 is the third constituent
	*/
	Energy massLightestHadronPair(tcPDPtr ptr1, tcPDPtr ptr2,
	tcPDPtr ptr3 = PDPtr ()) const {
	pair<tcPDPtr,tcPDPtr> pairData = lightestHadronPair(ptr1, ptr2, ptr3);
	if ( ! pairData.first \|\| ! pairData.second ) return ZERO;
	return ( pairData.first->mass() + pairData.second->mass() );
	}

	/**
	* Returns the lightest hadron formed by the given particles.
	*
	* Given the id of two (or three) constituents of a cluster, it returns
	* the lightest hadron with proper flavour numbers.
	* At the moment it does nothing in the case that also 'ptr3' present.
	* @param ptr1 is the first constituent
	* @param ptr2 is the second constituent
	* @param ptr3 is the third constituent
	*/
	tcPDPtr lightestHadron(tcPDPtr ptr1, tcPDPtr ptr2,
	#ifndef NDEBUG
	tcPDPtr ptr3 = PDPtr ()) const {
	#else
	tcPDPtr = PDPtr ()) const {
	#endif
	// The method assumes ptr3 == 0
	// rest not implemented
	assert(ptr1 && ptr2 && !ptr3);
	// find entry in the table
	pair<long,long> ids = make_pair(abs(ptr1->id()),abs(ptr2->id()));
	HadronTable::const_iterator tit=_table.find(ids);
	// throw exception if flavours wrong
	if (tit==_table.end())
	throw Exception() << "Could not find "
	<< ids.first << ' ' << ids.second
	<< " in _table. "
	<< "In HadronSelector::lightestHadron()"
	<< Exception::eventerror;
	if(tit->second.empty())
	throw Exception() << "HadronSelector::lightestHadron "
	<< "could not find any hadrons containing "
	<< ptr1->id() << ' ' << ptr2->id() << '\n'
	<< tit->first.first << ' '
	<< tit->first.second << Exception::eventerror;
	// find the lightest hadron
	int sign = signHadron(ptr1,ptr2,tit->second.begin()->ptrData);
	tcPDPtr candidate = sign > 0 ?
	tit->second.begin()->ptrData : tit->second.begin()->ptrData->CC();
	// \todo 20 GeV limit is temporary fudge to let SM particles go through.
	// \todo Use isExotic instead?
	if (candidate->mass() > 20*GeV
	&& candidate->mass() < ptr1->constituentMass() + ptr2->constituentMass()) {
	generator()->log() << "HadronSelector::lightestHadron: "
	<< "chosen candidate " << candidate->PDGName()
	<< " is lighter than its constituents "
	<< ptr1->PDGName() << ", " << ptr2->PDGName() << '\n'
	<< candidate->mass()/GeV << " < " << ptr1->constituentMass()/GeV
	<< " + " << ptr2->constituentMass()/GeV << '\n'
	<< "Check your particle data tables.\n";
	assert(false);
	}
	return candidate;
	}

	/**
	* Return the nominal mass of the hadron returned by lightestHadron()
	* @param ptr1 is the first constituent
	* @param ptr2 is the second constituent
	* @param ptr3 is the third constituent
	*/
	Energy massLightestHadron(tcPDPtr ptr1, tcPDPtr ptr2,
	#ifndef NDEBUG
	tcPDPtr ptr3 = PDPtr ()) const {
	#else
	tcPDPtr = PDPtr ()) const {
	#endif
	// The method assumes ptr3 == empty
	assert(!(ptr3));
	// find entry in the table
	pair<long,long> ids(abs(ptr1->id()),abs(ptr2->id()));
	HadronTable::const_iterator tit=_table.find(ids);
	// throw exception if flavours wrong
	if(tit==_table.end()\|\|tit->second.empty())
	throw Exception() << "HadronSelector::massLightestHadron "
	<< "failed for particle" << ptr1->id() << " "
	<< ptr2->id()
	<< Exception::eventerror;
	// return the mass
	return tit->second.begin()->mass;
	}

	/**
	* Returns the mass of the lightest pair of baryons.
	* @param ptr1 is the first constituent
	* @param ptr2 is the second constituent
	*/
	Energy massLightestBaryonPair(tcPDPtr ptr1, tcPDPtr ptr2) const;

	/**
	* Return the weights for the different quarks and diquarks
	*/
	//@{
	/**
	* The down quark weight.
	*/
	double pwtDquark() const {
	return _pwtDquark;
	}

	/**
	* The up quark weight.
	*/
	double pwtUquark() const {
	return _pwtUquark;
	}

	/**
	* The strange quark weight.
	*/
	double pwtSquark() const {
	return _pwtSquark;
	}

	/**
	* The charm quark weight.
	*/
	double pwtCquark() const {
	return _pwtCquark;
	}

	/**
	* The bottom quark weight.
	*/
	double pwtBquark() const {
	return _pwtBquark;
	}

	/**
	* The diquark weight.
	*/
	double pwtDIquark() const {
	return _pwtDIquark;
	}
	//@}

	public:

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

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

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

	protected:

	/** @name Standard Interfaced functions. */
	//@{
	/**
	* Initialize this object after the setup phase before saving an
	* EventGenerator to disk.
	*
	* The array _repwt is initialized using the interfaces to set different
	* weights for different meson multiplets and the constructHadronTable()
	* method called to complete the construction of the hadron tables.
	*
	* @throws InitException if object could not be initialized properly.
	*/
	virtual void doinit();
	//@}

	protected:

	/**
	* Construct the table of hadron data
	* This is the main method to initialize the hadron data (mainly the
	* weights associated to each hadron, taking into account its spin,
	* eventual isoscalar-octect mixing, singlet-decuplet factor). This is
	* the method that one should update when new or updated hadron data is
	* available.
	*
	* This class implements the construction of the basic table but can be
	* overridden if needed in inheriting classes.
	*
	* The rationale for factors used for diquarks involving different quarks can
	* be can be explained by taking a prototype example that in the exact SU(2) limit,
	* in which:
	* \f[m_u=m_d\f]
	* \f[M_p=M_n=M_\Delta\f]
	* and we will have equal numbers of u and d quarks produced.
	* Suppose that we weight 1 the diquarks made of the same
	* quark and 1/2 those made of different quarks, the fractions
	* of u and d baryons (p, n, Delta) we get are the following:
	* - \f$\Delta^{++}\f$: 1 possibility only u uu with weight 1
	* - \f$\Delta^- \f$: 1 possibility only d dd with weight 1
	* - \f$p,\Delta^+ \f$: 2 possibilities u ud with weight 1/2
	* d uu with weight 1
	* - \f$n,\Delta^0 \f$: 2 possibilities d ud with weight 1/2
	* u dd with weight 1
	* In the latter two cases, we have to take into account the
	* fact that p and n have spin 1/2 whereas Delta+ and Delta0
	* have spin 3/2 therefore from phase space we get a double weight
	* for Delta+ and Delta0 relative to p and n respectively.
	* Therefore the relative amount of these baryons that is
	* produced is the following:
	* # p = # n = ( 1/2 + 1 ) * 1/3 = 1/2
	* # Delta++ = # Delta- = 1 = ( 1/2 + 1) * 2/3 # Delta+ = # Delta0
	* which is correct, and therefore the weight 1/2 for the
	* diquarks of different types of quarks is justified (at least
	* in this limit of exact SU(2) ).
	*/
	virtual void constructHadronTable();

	/**
	* Access to the table of hadrons
	*/
	HadronTable & table() {
	return _table;
	}

	/**
	* Access to the list of partons
	*/
	vector<PDPtr> & partons() {
	return _partons;
	}

	/**
	* Access the parton weights
	*/
	- double pwt(tcPDPtr p) {
	- assert(p);
	- if ( p->id() < 0 ) p = p->CC();
	- map<tcPDPtr,double>::iterator it = _pwt.find(p);
	+ double pwt(long pid) {
	+ map<long,double>::iterator it = _pwt.find(abs(pid));
	assert( it != _pwt.end() );
	return it->second;
	}

	/**
	* Methods for the mixing of \f$I=0\f$ mesons
	*/
	//@{
	/**
	* Return the probability of mixing for Octet-Singlet isoscalar mixing,
	* the probability of the
	* \f$\frac1{\sqrt{2}}(\|u\bar{u}\rangle + \|d\bar{d}\rangle)\f$ component
	* is returned.
	* @param angleMix The mixing angle in degrees (not radians)
	* @param order is 0 for no mixing, 1 for the first resonance of a pair,
	* 2 for the second one.
	* The mixing is defined so that for example with \f$\eta-\eta'\f$ mixing where
	* the mixing angle is \f$\theta=-23^0$ with $\eta\f$ as the first particle
	* and \f$\eta'\f$ the second one.
	* The convention used is
	* \f[\eta = \cos\theta\|\eta_{\rm octet }\rangle
	* -\sin\theta\|\eta_{\rm singlet}\rangle\f]
	* \f[\eta' = \sin\theta\|\eta_{\rm octet }\rangle
	* -\cos\theta\|\eta_{\rm singlet}\rangle\f]
	* with
	* \f[\|\eta_{\rm singlet}\rangle = \frac1{\sqrt{3}}
	* \left[\|u\bar{u}\rangle + \|d\bar{d}\rangle + \|s\bar{s}\rangle\right]\f]
	* \f[\|\eta_{\rm octet }\rangle = \frac1{\sqrt{6}}
	* \left[\|u\bar{u}\rangle + \|d\bar{d}\rangle - 2\|s\bar{s}\rangle\right]\f]
	*/
	double probabilityMixing(const double angleMix,
	const int order) {
	static double convert=Constants::pi/180.0;
	if (order == 1)
	return sqr( cos( angleMix*convert + atan( sqrt(2.0) ) ) );
	else if (order == 2)
	return sqr( sin( angleMix*convert + atan( sqrt(2.0) ) ) );
	else
	return 1.;
	}

	/**
	* Returns the weight of given mixing state.
	* @param id The PDG code of the meson
	*/
	virtual double mixingStateWeight(long id);
	//@}

	/**
	* Calculates a special weight specific to a given hadron.
	* @param id The PDG code of the hadron
	*/
	double specialWeight(long id);

	/**
	* This method returns the proper sign ( > 0 hadron; < 0 anti-hadron )
	* for the input PDG id idHad > 0, suppose to be made by the
	* two constituent particle pointers: par1 and par2 (both with proper sign).
	*/
	int signHadron(tcPDPtr ptr1, tcPDPtr ptr2, tcPDPtr hadron) const;

	private:

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

	private:

	/**
	* The PDG codes of the constituent particles allowed
	*/
	vector<PDPtr> _partons;

	/**
	* The PDG codes of the hadrons which cannot be produced in the hadronization
	*/
	vector<PDPtr> _forbidden;

	/**
	* The weights for the different quarks and diquarks
	*/
	//@{
	/**
	* The probability of producting a down quark.
	*/
	double _pwtDquark;

	/**
	* The probability of producting an up quark.
	*/
	double _pwtUquark;

	/**
	* The probability of producting a strange quark.
	*/
	double _pwtSquark;

	/**
	* The probability of producting a charm quark.
	*/
	double _pwtCquark;

	/**
	* The probability of producting a bottom quark.
	*/
	double _pwtBquark;

	/**
	* The probability of producting a diquark.
	*/
	double _pwtDIquark;

	/**
	* Weights for quarks and diquarks.
	*/
	- map<tcPDPtr,double> _pwt;
	+ map<long,double> _pwt;
	//@}

	/**
	* The mixing angles for the \f$I=0\f$ mesons containing light quarks
	*/
	//@{
	/**
	* The \f$\eta-\eta'\f$ mixing angle
	*/
	double _etamix;

	/**
	* The \f$\phi-\omega\f$ mixing angle
	*/
	double _phimix;

	/**
	* The \f$h_1'-h_1\f$ mixing angle
	*/
	double _h1mix;

	/**
	* The \f$f_0(1710)-f_0(1370)\f$ mixing angle
	*/
	double _f0mix;

	/**
	* The \f$f_1(1420)-f_1(1285)\f$ mixing angle
	*/
	double _f1mix;

	/**
	* The \f$f'_2-f_2\f$ mixing angle
	*/
	double _f2mix;

	/**
	* The \f$\eta_2(1870)-\eta_2(1645)\f$ mixing angle
	*/
	double _eta2mix;

	/**
	* The \f$\phi(???)-\omega(1650)\f$ mixing angle
	*/
	double _omhmix;

	/**
	* The \f$\phi_3-\omega_3\f$ mixing angle
	*/
	double _ph3mix;

	/**
	* The \f$\eta(1475)-\eta(1295)\f$ mixing angle
	*/
	double _eta2Smix;

	/**
	* The \f$\phi(1680)-\omega(1420)\f$ mixing angle
	*/
	double _phi2Smix;
	//@}

	/**
	* The weights for the various meson multiplets to be used to supress the
	* production of particular states
	*/
	//@{
	/**
	* The weights for the \f$\phantom{1}^1S_0\f$ multiplets
	*/
	vector<double> _weight1S0;

	/**
	* The weights for the \f$\phantom{1}^3S_1\f$ multiplets
	*/
	vector<double> _weight3S1;

	/**
	* The weights for the \f$\phantom{1}^1P_1\f$ multiplets
	*/
	vector<double> _weight1P1;

	/**
	* The weights for the \f$\phantom{1}^3P_0\f$ multiplets
	*/
	vector<double> _weight3P0;

	/**
	* The weights for the \f$\phantom{1}^3P_1\f$ multiplets
	*/
	vector<double> _weight3P1;

	/**
	* The weights for the \f$\phantom{1}^3P_2\f$ multiplets
	*/
	vector<double> _weight3P2;

	/**
	* The weights for the \f$\phantom{1}^1D_2\f$ multiplets
	*/
	vector<double> _weight1D2;

	/**
	* The weights for the \f$\phantom{1}^3D_1\f$ multiplets
	*/
	vector<double> _weight3D1;

	/**
	* The weights for the \f$\phantom{1}^3D_2\f$ multiplets
	*/
	vector<double> _weight3D2;

	/**
	* The weights for the \f$\phantom{1}^3D_3\f$ multiplets
	*/
	vector<double> _weight3D3;
	//@}

	/**
	* The weights for the excited meson multiplets
	*/
	vector<vector<vector<double> > > _repwt;

	/**
	* Singlet and Decuplet weights
	*/
	//@{
	/**
	* The singlet weight
	*/
	double _sngWt;

	/**
	* The decuplet weight
	*/
	double _decWt;
	//@}

	/**
	* The table of hadron data
	*/
	HadronTable _table;

	/**
	* Enums so arrays can be statically allocated
	*/
	//@{
	/**
	* Defines values for array sizes. L,J,N max values for excited mesons.
	*/
	enum MesonMultiplets { Lmax = 3, Jmax = 4, Nmax = 4};
	//@}

	/**
	* Option for the construction of the tables
	*/
	unsigned int _topt;

	/**
	* Which particles to produce for debugging purposes
	*/
	unsigned int _trial;
	};


	}

	#endif /* HERWIG_HadronSelector_H */
	diff --git a/Hadronization/Hw64Selector.cc b/Hadronization/Hw64Selector.cc
	--- a/Hadronization/Hw64Selector.cc
	+++ b/Hadronization/Hw64Selector.cc
	@@ -1,109 +1,109 @@
	// -- C++ --
	//
	// Hw64Selector.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the Hw64Selector class.
	//

	#include "Hw64Selector.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Repository/UseRandom.h"
	#include "CheckId.h"
	#include "Herwig++/Utilities/Kinematics.h"
	#include <ThePEG/Utilities/DescribeClass.h>

	using namespace Herwig;

	DescribeNoPIOClass<Hw64Selector,HadronSelector>
	describeHw64Selector("Herwig::Hw64Selector","");

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

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

	void Hw64Selector::Init() {

	static ClassDocumentation<Hw64Selector> documentation
	("There is no documentation for the Hw64Selector class");

	}

	pair<tcPDPtr,tcPDPtr> Hw64Selector::chooseHadronPair(const Energy cluMass,tcPDPtr par1,
	tcPDPtr par2,tcPDPtr)
	{
	bool diquark = !(DiquarkMatcher::Check(par1->id()) \|\| DiquarkMatcher::Check(par2->id()));



	pair<tcPDPtr,tcPDPtr> lighthad = lightestHadronPair(par1, par2);
	if(!lighthad.first \|\| !lighthad.second)
	throw Exception() << "Hw64Selector::chooseHadronPair "
	<< "We have 0's! First id = " << par1->id() << " second = "
	<< par2->id() << ". This is probably a problem with either"
	<< " undecayed heavy particles or colour connections"
	<< Exception::eventerror;
	// calculate maximum momentum
	Energy PCMax = Kinematics::pstarTwoBodyDecay(cluMass,lighthad.first->mass(),
	lighthad.second->mass());
	tcPDPtr had1 = tcPDPtr();
	tcPDPtr had2 = tcPDPtr();
	int ntry = 0;
	tcPDPtr quark = tcPDPtr();
	const int nmax = 5000;
	Energy p;
	do {
	quark = partons()[UseRandom::irnd(partons().size())];
	if(diquark && DiquarkMatcher::Check(quark->id())) continue;
	KupcoData::const_iterator it1,it2;
	- if(pwt(quark) <= UseRandom::rnd()) continue;
	+ if(pwt(quark->id()) <= UseRandom::rnd()) continue;
	pair<long,long> pid(abs(par1->id()),quark->id());
	do {
	it1 = table()[pid].begin();
	advance(it1,int(double(table()[pid].size())*UseRandom::rnd()));
	}
	while(it1 != table()[pid].end() && it1->overallWeight < UseRandom::rnd());

	had1 = it1->ptrData;
	pid = make_pair(quark->id(),abs(par2->id()));
	do {
	it2 = table()[pid].begin();
	advance(it2,int(double(table()[pid].size())*UseRandom::rnd()));
	}
	while(it2 != table()[pid].end() && it2->overallWeight < UseRandom::rnd());
	had2 = it2->ptrData;
	if(had1 && had2) {
	p = Kinematics::pstarTwoBodyDecay(cluMass, it1->mass, it2->mass);
	if(p/PCMax < UseRandom::rnd()) {
	had1 = had2 = tcPDPtr();
	ntry++;
	}
	}
	}
	while((!had1\|\| !had2) && ntry < nmax);
	if(ntry >= nmax) return lighthad;
	int signHad1 = 0;
	int signHad2 = 0;
	if(CheckId::canBeHadron(par1,quark->CC()) && CheckId::canBeHadron(quark,par2)) {
	signHad1 = signHadron(par1, quark->CC(), had1);
	signHad2 = signHadron(par2, quark, had2);
	}
	else if(CheckId::canBeHadron(par1,quark) && CheckId::canBeHadron(quark->CC(),par2)) {
	signHad1 = signHadron(par1, quark, had1);
	signHad2 = signHadron(par2, quark->CC(), had2);
	}

	else throw Exception() << "Hw64Selector::chooseHadronPair()"
	<< Exception::abortnow;
	return make_pair( signHad1 > 0 ? had1 : tcPDPtr(had1->CC()),
	signHad2 > 0 ? had2 : tcPDPtr(had2->CC()));
	}
	diff --git a/Hadronization/HwppSelector.cc b/Hadronization/HwppSelector.cc
	--- a/Hadronization/HwppSelector.cc
	+++ b/Hadronization/HwppSelector.cc
	@@ -1,188 +1,184 @@
	// -- C++ --
	//
	// HwppSelector.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the HwppSelector class.
	//

	#include "HwppSelector.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "Herwig++/Utilities/Kinematics.h"
	#include "ThePEG/Utilities/Selector.h"
	#include "ThePEG/Repository/UseRandom.h"
	#include "CheckId.h"
	#include <cassert>
	#include <ThePEG/Utilities/DescribeClass.h>

	using namespace Herwig;

	DescribeClass<HwppSelector,HadronSelector>
	describeHwppSelector("Herwig::HwppSelector","");

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

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

	void HwppSelector::doinit() {
	HadronSelector::doinit();
	}

	void HwppSelector::persistentOutput(PersistentOStream & os) const {
	os << _mode;
	}

	void HwppSelector::persistentInput(PersistentIStream & is, int) {
	is >> _mode;
	}

	void HwppSelector::Init() {

	static ClassDocumentation<HwppSelector> documentation
	("The HwppSelector class implements the Herwig++ algorithm for selecting"
	" the hadrons",
	"The hadronization used the selection algorithm described in \\cite{Kupco:1998fx}.",
	"%\\cite{Kupco:1998fx}\n"
	"\\bibitem{Kupco:1998fx}\n"
	" A.~Kupco,\n"
	" ``Cluster hadronization in HERWIG 5.9,''\n"
	" arXiv:hep-ph/9906412.\n"
	" %%CITATION = HEP-PH/9906412;%%\n"
	);
	// put useMe() only in correct place!


	static Switch<HwppSelector,unsigned int> interfaceMode
	("Mode",
	"Which algorithm to use",
	&HwppSelector::_mode, 1, false, false);
	static SwitchOption interfaceModeKupco
	(interfaceMode,
	"Kupco",
	"Use the Kupco approach",
	0);
	static SwitchOption interfaceModeHwpp
	(interfaceMode,
	"Hwpp",
	"Use the Herwig++ approach",
	1);

	}

	pair<tcPDPtr,tcPDPtr> HwppSelector::chooseHadronPair(const Energy cluMass,tcPDPtr par1,
	tcPDPtr par2,tcPDPtr )
	{
	// if either of the input partons is a diquark don't allow diquarks to be
	// produced
	bool diquark = !(DiquarkMatcher::Check(par1->id()) \|\| DiquarkMatcher::Check(par2->id()));
	bool quark = true;
	// if the Herwig++ algorithm
	if(_mode ==1) {
	if(cluMass > massLightestBaryonPair(par1,par2) &&
	UseRandom::rnd() > 1./(1.+pwtDIquark())) {
	diquark = true;
	quark = false;
	}
	else {
	useMe();
	diquark = false;
	quark = true;
	}
	}
	// weights for the different possibilities
	Energy weight, wgtsum(ZERO);
	// loop over all hadron pairs with the allowed flavours
	vector<Kupco> hadrons;
	+ hadrons.reserve(partons().size());
	for(unsigned int ix=0;ix<partons().size();++ix) {
	tcPDPtr quarktopick = partons()[ix];
	if(!quark && abs(int(quarktopick->iColour())) == 3
	&& !DiquarkMatcher::Check(quarktopick->id())) continue;
	if(!diquark && abs(int(quarktopick->iColour())) == 3
	&& DiquarkMatcher::Check(quarktopick->id())) continue;
	HadronTable::const_iterator
	tit1 = table().find(make_pair(abs(par1->id()),quarktopick->id()));
	HadronTable::const_iterator
	tit2 = table().find(make_pair(quarktopick->id(),abs(par2->id())));
	// If not in table skip
	if(tit1 == table().end()\|\|tit2==table().end()) continue;
	// tables empty skip
	- if(tit1->second.empty()\|\|tit2->second.empty()) continue;
	+ const KupcoData & T1 = tit1->second;
	+ const KupcoData & T2 = tit2->second;
	+ if(T1.empty()\|\|T2.empty()) continue;
	// if too massive skip
	- if(cluMass <= tit1->second.begin()->mass +
	- tit2->second.begin()->mass) continue;
	+ if(cluMass <= T1.begin()->mass +
	+ T2.begin()->mass) continue;
	// loop over the hadrons
	- KupcoData::iterator H1,H2;
	- for(H1 = tit1->second.begin();H1 != tit1->second.end(); ++H1) {
	- for(H2 = tit2->second.begin();H2 != tit2->second.end(); ++H2) {
	+ KupcoData::const_iterator H1,H2;
	+ for(H1 = T1.begin();H1 != T1.end(); ++H1) {
	+ for(H2 = T2.begin();H2 != T2.end(); ++H2) {
	// break if cluster too light
	if(cluMass < H1->mass + H2->mass) break;
	// calculate the weight
	- weight = pwt(quarktopick) * H1->overallWeight * H2->overallWeight *
	+ weight = pwt(quarktopick->id()) * H1->overallWeight * H2->overallWeight *
	Kinematics::pstarTwoBodyDecay(cluMass, H1->mass, H2->mass );
	int signQ = 0;
	assert (par1 && quarktopick);
	assert (par2);

	assert(quarktopick->CC());

	if(CheckId::canBeHadron(par1, quarktopick->CC())
	&& CheckId::canBeHadron(quarktopick, par2))
	signQ = +1;
	else if(CheckId::canBeHadron(par1, quarktopick)
	&& CheckId::canBeHadron(quarktopick->CC(), par2))
	signQ = -1;
	else {
	cerr << "Could not make sign for" << par1->id()<< " " << quarktopick->id()
	<< " " << par2->id() << "\n";
	assert(false);
	}

	if (signQ == -1)
	quarktopick = quarktopick->CC();
	// construct the object with the info
	Kupco a(quarktopick, H1->ptrData, H2->ptrData, weight);
	hadrons.push_back(a);
	wgtsum += weight;
	}
	}
	}
	if (hadrons.empty())
	return make_pair(tcPDPtr(),tcPDPtr());
	// select the hadron
	wgtsum *= UseRandom::rnd();
	unsigned int ix=0;
	do {
	wgtsum-= hadrons[ix].weight;
	++ix;
	}
	while(wgtsum > ZERO && ix < hadrons.size());
	if(ix == hadrons.size() && wgtsum > ZERO)
	return make_pair(tcPDPtr(),tcPDPtr());
	--ix;
	assert(hadrons[ix].idQ);
	int signHad1 = signHadron(par1, hadrons[ix].idQ->CC(), hadrons[ix].hadron1);
	int signHad2 = signHadron(par2, hadrons[ix].idQ, hadrons[ix].hadron2);
	- assert(!( signHad1 == 0 \|\| signHad2 == 0));
	-// throw Exception() << "HwppSelector::selectPair "
	-// << "***Inconsistent Hadron "
	-// << hadrons[ix].idQ->id() << " "
	-// << hadrons[ix].hadron1->id() << " "
	-// << hadrons[ix].hadron2->id() << " "
	-// << signHad1 << " " << signHad2
	-// << Exception::runerror;
	+ assert( signHad1 != 0 && signHad2 != 0 );
	return make_pair
	( signHad1 > 0 ? hadrons[ix].hadron1 : tcPDPtr(hadrons[ix].hadron1->CC()),
	signHad2 > 0 ? hadrons[ix].hadron2 : tcPDPtr(hadrons[ix].hadron2->CC()));
	}

File Metadata

Mime Type: text/x-diff
Expires: Sat, May 3, 5:40 AM (5 h, 31 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4982793
Default Alt Text: (66 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions