No OneTemporary
Actions

Size

59 KB

Subscribers

None

View Options

	diff --git a/Decay/EvtGen/EvtGenInterface.cc b/Decay/EvtGen/EvtGenInterface.cc
	--- a/Decay/EvtGen/EvtGenInterface.cc
	+++ b/Decay/EvtGen/EvtGenInterface.cc
	@@ -1,1012 +1,1025 @@
	// -- C++ --
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the EvtGenInterface class.
	//

	#include "EvtGenInterface.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/ParVector.h"
	#include "ThePEG/EventRecord/Particle.h"
	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/PDT/DecayMode.h"
	#include "Herwig/Decay/DecayVertex.h"
	#include "Herwig/Decay/GeneralDecayMatrixElement.h"

	#include "EvtGenBase/EvtRandom.hh"
	#include "EvtGenBase/EvtAbsRadCorr.hh"
	#include "EvtGenExternal/EvtExternalGenList.hh"
	#include "EvtGenBase/EvtScalarParticle.hh"
	#include "EvtGenBase/EvtDiracParticle.hh"
	#include "EvtGenBase/EvtVectorParticle.hh"
	#include "EvtGenBase/EvtRaritaSchwingerParticle.hh"
	#include "EvtGenBase/EvtTensorParticle.hh"
	#include "EvtGenBase/EvtStringParticle.hh"
	#include "EvtGenBase/EvtHighSpinParticle.hh"
	#include "EvtGenBase/EvtDecayTable.hh"

	#ifndef EVTGEN_PREFIX
	#error Makefile.am needs to define EVTGEN_PREFIX
	#endif

	using namespace Herwig;

	namespace {

	const string prefix=EVTGEN_PREFIX "";
	const string p8data=PYTHIA8DATA "";

	}

	EvtGenInterface::EvtGenInterface() : decayName_(prefix+"/share/DECAY_2010.DEC"),
	pdtName_(prefix+"/share/evt.pdl"),
	reDirect_(true), checkConv_(false),
	p8Data_(p8data)
	{}

	+EvtGenInterface::EvtGenInterface(const EvtGenInterface & x)
	+ : Interfaced(x), decayName_(x.decayName_), pdtName_(x.pdtName_),
	+ userDecays_(x.userDecays_), reDirect_(x.reDirect_),
	+ checkConv_(x.checkConv_), convID_(x.convID_),
	+ p8Data_(x.p8Data_), evtrnd_(x.evtrnd_),evtgen_(x.evtgen_)
	+{}
	+
	EvtGenInterface::~EvtGenInterface() {}

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

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

	void EvtGenInterface::persistentOutput(PersistentOStream & os) const {
	os << decayName_ << pdtName_ << reDirect_
	<< userDecays_ << checkConv_ << convID_ << p8Data_;
	}

	void EvtGenInterface::persistentInput(PersistentIStream & is, int) {
	is >> decayName_ >> pdtName_ >> reDirect_
	>> userDecays_ >> checkConv_ >> convID_ >> p8Data_;
	}


	// The following static variable is needed for the type
	// description system in ThePEG.
	DescribeClass<EvtGenInterface,Interfaced>
	describeHerwigEvtGenInterface("Herwig::EvtGenInterface",
	"HwEvtGenInterface.so");

	void EvtGenInterface::Init() {

	static ClassDocumentation<EvtGenInterface> documentation
	("The EvtGenInterface class is the main class for the use of the EvtGen "
	"decay package with Herwig");

	static Parameter<EvtGenInterface,string> interfaceDecay_File
	("Decay_File",
	"The name of the file for the EvtGen decays.",
	&EvtGenInterface::decayName_, prefix+"/share/DECAY_2010.DEC",
	false, false);

	static Parameter<EvtGenInterface,string> interfacePDT_File
	("PDTFile",
	"The name of the file for the EvtGen particle data.",
	&EvtGenInterface::pdtName_, prefix+"/share/evt.pdl",
	false, false);

	static Switch<EvtGenInterface,bool> interfaceRedirect
	("Redirect",
	"By default cerr and cout are redirected when EvtGen is running to"
	" allow us to catch errors, due to EvtGen's poor internal error "
	"handling. This can be a problem for debugging and so can be switched off.",
	&EvtGenInterface::reDirect_, true, false, false);
	static SwitchOption interfaceRedirectYes
	(interfaceRedirect,
	"Yes",
	"Redirect the output",
	true);
	static SwitchOption interfaceRedirectNo
	(interfaceRedirect,
	"No",
	"Don't redirect the output",
	false);

	static Switch<EvtGenInterface,bool> interfaceCheckConversion
	("CheckConversion",
	"Check the conversion of particles to and from EvtGen",
	&EvtGenInterface::checkConv_, false, false, false);
	static SwitchOption interfaceCheckConversionfalse
	(interfaceCheckConversion,
	"No",
	"Don't check the conversion",
	false);
	static SwitchOption interfaceCheckConversionCheck
	(interfaceCheckConversion,
	"Yes",
	"Check the conversion",
	true);

	static ParVector<EvtGenInterface,long> interfaceOutputModes
	("OutputModes",
	"Particles for which to output the EvtGen decay modes"
	" so they can be read into Herwig++",
	&EvtGenInterface::convID_, -1, long(0), 0, 0,
	false, false, Interface::nolimits);

	static ParVector<EvtGenInterface,string> interfaceUserDecays
	("UserDecays",
	"List of user decay files to be loaded",
	&EvtGenInterface::userDecays_, -1, "", "", "",
	false, false, Interface::nolimits);

	static Parameter<EvtGenInterface,string> interfacePythia8Data
	("Pythia8Data",
	"Location of the Pythia8 data directory",
	&EvtGenInterface::p8Data_, p8data,
	false, false);

	}

	void EvtGenInterface::doinitrun() {
	Interfaced::doinitrun();
	// redirect cerr and cout if needed
	std::streambuf *temp[2]={cout.rdbuf(),cerr.rdbuf()};
	- if(reDirect_) {
	- std::streambuf *psbuf = generator()->log().rdbuf();
	+ if(reDirect_ && ! generator()->useStdOut() ) {
	+ logFile_.open("EvtGen.log");
	+ std::streambuf *psbuf = logFile_.rdbuf();
	cout.rdbuf(psbuf);
	cerr.rdbuf(psbuf);
	}
	// output the EvtGen initialization info to the log file
	- generator()->log() << "Initializing EvtGen \n";
	+ cout << "Initializing EvtGen \n";
	// set up the random number generator for EvtGen
	- generator()->log() << "Setting EvtGen random number generator"
	+ cout << "Setting EvtGen random number generator"
	<< " to the Herwig one.\n";
	evtrnd_ = new EvtGenRandom(const_ptr_cast<Ptr<RandomGenerator>::pointer>
	(&(UseRandom::current())));
	EvtRandom::setRandomEngine(evtrnd_);
	// PHOTOS STUFF
	EvtAbsRadCorr* radCorrEngine = 0;
	std::list<EvtDecayBase*> extraModels;
	EvtExternalGenList genList(true,p8Data_,"gamma",true);
	radCorrEngine = genList.getPhotosModel();
	extraModels = genList.getListOfModels();
	// Initialize EvtGen
	evtgen_ = new EvtGen(decayName_.c_str(),pdtName_.c_str(),
	evtrnd_, radCorrEngine, &extraModels);
	// additional user decays
	for(unsigned int ix=0;ix<userDecays_.size();++ix)
	evtgen_->readUDecay(userDecays_[ix].c_str());
	// check the conversion of particles if needed
	if(checkConv_) checkConversion();
	// convert any particle decay modes which are needed
	for(unsigned int ix=0;ix<convID_.size();++ix) {
	outputEvtGenDecays(convID_[ix]);
	}
	// that's the lot
	- generator()->log() << "Finished initialisation of EvtGen \n";
	- if(reDirect_) {
	+ cout << "Finished initialisation of EvtGen \n";
	+ if(reDirect_ && ! generator()->useStdOut() ) {
	cout.rdbuf(temp[0]);
	cerr.rdbuf(temp[1]);
	}
	}

	+void EvtGenInterface::dofinish() {
	+ Interfaced::dofinish();
	+ if(logFile_.is_open())
	+ logFile_.close();
	+}

	ParticleVector EvtGenInterface::decay(const Particle &parent,
	bool recursive, const DecayMode & dm) const {
	// redirect cout to the log file
	ostringstream stemp;
	std::streambuf *temp[2]={cout.rdbuf(),cerr.rdbuf()};
	- if(reDirect_) {
	- std::streambuf *psbuf[2] ={generator()->log().rdbuf(),stemp.rdbuf()};
	+ if(reDirect_ && ! generator()->useStdOut() ) {
	+ std::streambuf *psbuf[2] ={logFile_.rdbuf(),stemp.rdbuf()};
	cout.rdbuf(psbuf[0]);
	cerr.rdbuf(psbuf[1]);
	}
	// generate the decay
	ParticleVector output;
	EvtId parID(EvtGenID(parent.id()));
	try {
	EvtDecayBase *decayer = NULL;
	// create evtgen particle for the parent
	EvtParticle* particle = EvtGenParticle(parent);
	EvtParticle* evtParent = NULL;
	// special if parent is the same to avoid doing mixing twice
	if(parent.parents().size()==1 && abs(parent.parents()[0]->id())==abs(parent.id())) {
	evtParent = EvtGenParticle(*parent.parents()[0]);
	particle->addDaug(evtParent);
	}
	// simplest case, evtgen selects mode and recursively decays, use their standard mechanism
	if(dm.wildProductMatcher() && recursive) {
	evtgen_->generateDecay(particle);
	}
	// otherwise we're in control
	else {
	// select the EvtGen decayer to use
	if(dm.wildProductMatcher()) {
	decayer = EvtDecayTable::getInstance()->getDecayFunc(particle);
	}
	// otherwise we should pick one
	else {
	int imode(EvtGenChannel(dm));
	decayer = EvtDecayTable::getInstance()->getDecay(parID.getAlias(),imode);
	particle->setChannel(imode);
	}
	// must be a decayer
	if(!decayer) throw Exception() << "Could find EvtGen decayer in EvtGen::decay()"
	<< Exception::runerror;
	// masses of the children
	bool massTreeOK(true);
	if ( particle->getNDaug() == 0 )
	massTreeOK = particle->generateMassTree();
	if ( ! massTreeOK ) {
	// delete the EvtGen particle
	particle->deleteDaughters();
	delete particle;
	particle = 0;
	if(evtParent) {
	delete evtParent;
	}
	throw Exception() << "EvtGen could not decay " << EvtPDL::name(particle->getId())
	<<" with mass "<< particle->mass()
	<<" to decay channel number "<< particle->getChannel() << "\n"
	<< Exception::eventerror;
	}
	assert(decayer !=0 );
	// perform decay
	decayer->makeDecay(particle,recursive);
	}
	// cast the decayer
	EvtDecayAmp damp = dynamic_cast<EvtDecayAmp>(decayer);
	// translate the decay products
	output=decayProducts(particle);
	// set spin information if needed
	tSpinPtr pspin(getSpinInfo(parent));
	// if has spin information translate it
	if(damp) {
	ParticleVector products;
	unsigned int nbeforerad=decayer->getNDaug();
	for(unsigned int ix=0;ix<nbeforerad;++ix) products.push_back(output[ix]);
	constructVertex(parent,products,damp);
	}
	// otherwise
	else {
	pspin->develop();
	RhoDMatrix rhotemp(pspin->iSpin());
	pspin->DMatrix()=rhotemp;
	}
	pspin->decay();
	if(recursive) {
	pspin->developed();
	pspin->DMatrix()=ThePEGSpinDensity(particle->getSpinDensityBackward(),
	ThePEGID(particle->getId()));
	}
	// delete the EvtGen particle
	particle->deleteDaughters();
	delete particle;
	if(evtParent) {
	delete evtParent;
	}
	particle = 0;
	}
	catch ( ... ) {
	// change stream back
	- if(reDirect_) {
	+ if(reDirect_ && ! generator()->useStdOut() ) {
	cout.rdbuf(temp[0]);
	cerr.rdbuf(temp[1]);
	}
	throw;
	}
	// change stream back
	- if(reDirect_) {
	+ if(reDirect_ && ! generator()->useStdOut() ) {
	cout.rdbuf(temp[0]);
	cerr.rdbuf(temp[1]);
	string stemp2=stemp.str();
	if(stemp2.length()>0)
	throw Exception() << "EvtGen report error in EvtGen::decay "
	<< "killing event\n"
	<< "Error was " << stemp2
	<< Exception::eventerror;
	}
	// return the decay products
	return output;
	}

	EvtParticle * EvtGenInterface::EvtGenParticle(const Particle & part) const {
	// convert the momentum
	Lorentz5Momentum inmom(part.momentum());
	EvtVector4R p4(EvtGenMomentum(inmom));
	EvtParticle *evtpart;
	// EvtGen ID and spin type
	EvtId id = EvtGenID(part.id());
	EvtSpinType::spintype thisSpin=EvtPDL::getSpinType(id);
	// boost to the rest frame to get the basis states right
	PPtr decay(const_ptr_cast<PPtr>(new_ptr(part)));
	decay->transform(LorentzRotation(-inmom.boostVector()));
	// get spin information so can transfer to EvtGen
	tcSpinPtr spin(dynamic_ptr_cast<tcSpinPtr>(part.spinInfo()));
	if(spin) spin->decay();
	// don't convert neutrinos, aren't decays so waste of time
	if ( thisSpin == EvtSpinType::NEUTRINO) {
	throw Exception() << "Tried to convert a neutrino to EvtGen in "
	<< "EvtGen::EvtParticle. This should not be needed as"
	<< "EvtGen does not decay neutrinos"
	<< Exception::eventerror;
	}
	// don't convert photons, aren't decays so waste of time
	else if ( thisSpin == EvtSpinType::PHOTON ) {
	throw Exception() << "Tried to convert a photon to EvtGen in "
	<< "EvtGen::EvtParticle. This should not be needed as"
	<< "EvtGen does not decay photons"
	<< Exception::eventerror;
	}
	// scalar particles
	else if ( thisSpin == EvtSpinType::SCALAR ) {
	// create particle
	EvtScalarParticle * myPart = new EvtScalarParticle;
	myPart->init(id, p4);
	// set rho matrix if needed
	tcScalarSpinPtr sp(dynamic_ptr_cast<tcScalarSpinPtr>(spin));
	if(sp) {
	myPart->setSpinDensityForward(EvtGenSpinDensity(sp->rhoMatrix()));
	}
	else {
	EvtSpinDensity rho;
	rho.setDiag(EvtSpinType::getSpinStates(thisSpin));
	myPart->setSpinDensityForward(rho);
	}
	evtpart=myPart;
	}
	else if ( thisSpin == EvtSpinType::DIRAC ) {
	EvtDiracParticle *myPart =new EvtDiracParticle;
	// get the spin info
	tcFermionSpinPtr sp(dynamic_ptr_cast<tcFermionSpinPtr>(spin));
	// has spin info transfer
	if(sp) {
	vector<EvtDiracSpinor> prod,decay;
	// transfer spinors
	for(unsigned int ix=0;ix<2;++ix) {
	prod .push_back(EvtGenSpinor(sp->getProductionBasisState(ix)));
	decay.push_back(EvtGenSpinor(sp->getDecayBasisState (ix)));
	}
	myPart->init(id, p4,prod[0],prod[1],decay[0],decay[1]);
	// set the density matrix
	myPart->setSpinDensityForward(EvtGenSpinDensity(sp->rhoMatrix()));
	}
	// no spin info
	else {
	myPart->init(id, p4);
	EvtSpinDensity rho;
	rho.setDiag(EvtSpinType::getSpinStates(thisSpin));
	myPart->setSpinDensityForward(rho);
	}
	evtpart=myPart;
	}
	// vector particles
	else if ( thisSpin == EvtSpinType::VECTOR ) {
	EvtVectorParticle *myPart=new EvtVectorParticle;
	// get the spin info
	tcVectorSpinPtr sp(dynamic_ptr_cast<tcVectorSpinPtr>(spin));
	// has spin info transfer
	if(sp) {
	// transfer polarization vectors
	vector<EvtVector4C> eps;
	for(unsigned int ix=0;ix<3;++ix) {
	eps.push_back(EvtGenPolarization(sp->getDecayBasisState(ix)));
	}
	myPart->init(id, p4,eps[0],eps[1],eps[2]);
	// set spin density matrix
	myPart->setSpinDensityForward(EvtGenSpinDensity(sp->rhoMatrix()));
	}
	// no spin info
	else {
	myPart->init(id, p4);
	EvtSpinDensity rho;
	rho.setDiag(EvtSpinType::getSpinStates(thisSpin));
	myPart->setSpinDensityForward(rho);
	}
	evtpart=myPart;
	}
	// spin 3/2 particles
	else if ( thisSpin == EvtSpinType::RARITASCHWINGER ) {
	EvtRaritaSchwingerParticle *myPart=new EvtRaritaSchwingerParticle;
	// get the spin info
	tcRSFermionSpinPtr sp(dynamic_ptr_cast<tcRSFermionSpinPtr>(spin));
	// has spin info transfer
	if(sp) {
	// transfer spinors
	vector<EvtRaritaSchwinger> prod,decay;
	for(unsigned int ix=0;ix<4;++ix) {
	prod .push_back(EvtGenRSSpinor(sp->getProductionBasisState(ix)));
	decay.push_back(EvtGenRSSpinor(sp->getDecayBasisState(ix) ));
	}
	myPart->init(id, p4,
	prod[0] ,prod[1] ,prod[2] , prod[3],
	decay[0],decay[1],decay[2],decay[3]);
	// set spin density matrix
	myPart->setSpinDensityForward(EvtGenSpinDensity(sp->rhoMatrix()));
	}
	// no spin info
	else {
	myPart->init(id, p4);
	EvtSpinDensity rho;
	rho.setDiag(EvtSpinType::getSpinStates(thisSpin));
	myPart->setSpinDensityForward(rho);
	}
	evtpart=myPart;
	}
	// spin-2 particles
	else if ( thisSpin == EvtSpinType::TENSOR ) {
	EvtTensorParticle *myPart =new EvtTensorParticle;
	tcTensorSpinPtr sp(dynamic_ptr_cast<tcTensorSpinPtr>(spin));
	// has spin info transfer
	if(sp) {
	// transfer the polarization tensors
	vector<EvtTensor4C> eps;
	for(unsigned int ix=0;ix<5;++ix) {
	eps.push_back(EvtGenTensor(sp->getDecayBasisState(ix)));
	}
	myPart->init(id, p4,eps[0],eps[1],eps[2],eps[3],eps[4]);
	// set spin density matrix
	myPart->setSpinDensityForward(EvtGenSpinDensity(sp->rhoMatrix()));
	}
	// no spin info
	else {
	myPart->init(id, p4);
	EvtSpinDensity rho;
	rho.setDiag(EvtSpinType::getSpinStates(thisSpin));
	myPart->setSpinDensityForward(rho);
	}
	evtpart=myPart;
	}
	// shouldn't be doing this but here for safety
	else if ( thisSpin == EvtSpinType::STRING ) {
	EvtStringParticle *myPart;
	myPart=new EvtStringParticle;
	myPart->init(id, p4);
	EvtSpinDensity rho;
	rho.setDiag(EvtSpinType::getSpinStates(thisSpin));
	myPart->setSpinDensityForward(rho);
	evtpart=myPart;
	}
	// high spin particles
	else if ( thisSpin == EvtSpinType::SPIN3 \|\| thisSpin == EvtSpinType::SPIN5HALF \|\|
	thisSpin == EvtSpinType::SPIN4 \|\| thisSpin == EvtSpinType::SPIN7HALF) {
	EvtHighSpinParticle *myPart;
	myPart=new EvtHighSpinParticle;
	myPart->init(id, p4);
	EvtSpinDensity rho;
	rho.setDiag(EvtSpinType::getSpinStates(thisSpin));
	myPart->setSpinDensityForward(rho);
	evtpart=myPart;
	}
	else {
	throw Exception() << "Can't convert particle to EvtGen "
	<< part << Exception::eventerror;
	}
	// boost particle back and return the particle
	decay->boost(inmom.boostVector());
	return evtpart;
	}

	// convert from ThePEG to EvtGen
	EvtId EvtGenInterface::EvtGenID(int id, bool exception) const {
	EvtId output;
	int absid=abs(id),ispin(absid%10),isgn(id/absid);
	// handle the easy cases
	if(
	//quarks(+excited)
	absid<=8\|\|
	// leptons(+exicted)
	(absid>=11&&absid<=18)\|\|
	// SM(+extensions) gauge bosons and higgs
	(absid>=21&&absid<=25)\|\|(absid>=32&&absid<=37)\|\|
	// 1 1S0, 1 3S1 and 1 3P2 1 3D3 mesons are the same
	(absid>100&&absid<600&&ispin%2==1&&ispin<=9)\|\|
	// 2 1S0, 2 3S1 are the same
	(absid>100100&&absid<100600&&(ispin==1\|\|ispin==3))\|\|
	// 1 3d1 goes to 3s in evtgen
	(absid>30100&&absid<30600&&ispin==3) \|\|
	// 1 3p0 most same
	(absid>10100&&absid<10600&&ispin==1) \|\|
	// mixed kaons and diffractive states
	(absid>=100&&absid<=3000&&ispin==0)) {
	output = EvtPDL::evtIdFromStdHep(id);
	}
	// lowest baryon multiplets and diquarks are almost the same
	else if(absid>1000&&absid<6000&&ispin>=1&&ispin<=4) {
	output = EvtPDL::evtIdFromStdHep(id);
	}
	// 1 1P1 mesons are the same apart from D_s1
	else if(absid>10100&&absid<10600&&(ispin==3)) {
	if(absid==10433) output = EvtPDL::evtIdFromStdHep(isgn*20433);
	else output = EvtPDL::evtIdFromStdHep(id);
	}
	// 1 3P1 mesons are the sameapart from D_s1
	else if(absid>20100&&absid<20600&&(ispin==3)) {
	if(absid==20433) output = EvtPDL::evtIdFromStdHep(isgn*10433);
	else output = EvtPDL::evtIdFromStdHep(id);
	}
	// excited baryons
	else if(((absid>10000&&absid<16000)\|\|(absid>20000&&absid<26000)\|\|
	(absid>30000&&absid<36000))&&(ispin==2\|\|ispin==4)) {
	output=EvtPDL::evtIdFromStdHep(id);
	}
	// special for any bottomonium states
	else if((absid%1000)/10==55) {
	output=EvtPDL::evtIdFromStdHep(id);
	}
	// special meson codes
	else if (absid>9000000) {
	// f_0/a_0(980) goes into f_0/a_0(980)
	if(absid==9000111\|\|absid==9000211\|\|absid==9010221) {
	output=EvtPDL::evtIdFromStdHep(id);
	}
	// sigma
	else if(absid==9000221) {
	output=EvtPDL::evtIdFromStdHep(id);
	}
	// psi(4040)
	else if(absid==9000443) {
	output=EvtPDL::evtIdFromStdHep(50443);
	}
	// f_0(1500)
	else if(absid==9030221) {
	output=EvtPDL::evtIdFromStdHep(9020221);
	}
	}
	// check its O.K.
	if(output.getAlias()==-1&&exception) {
	throw Exception() << "Can't find the EvtGen Id for particle "
	<< getParticleData(id)->PDGName()
	<< " with PDG code = "
	<< id << " in EvtGen::EvtGenID"
	<< Exception::eventerror;
	}
	return output;
	}

	ParticleVector EvtGenInterface::decayProducts(EvtParticle *part, bool boost) const {
	ParticleVector output,temp;
	for(unsigned int ix=0,N=part->getNDaug();ix<N;++ix) {
	// get the daughter particle
	EvtParticle * daug = part->getDaug(ix);
	// may just have been used for mass generation, so check has valid momentum
	if(!daug->hasValidP4()) continue;
	// get the Herwig ParticleData object
	int id=ThePEGID(daug->getId());
	tcPDPtr pd=getParticleData(id);
	// if easy to convert do it
	if(pd) output.push_back(ThePEGParticle(daug,pd));
	// special for EvtGen particles like string we don't need to include
	else if(id==90) {
	// check if needs to be decayed by EvtGen
	if(EvtPDL::getStdHep(daug->getId())==92\|\|daug->isDecayed()) {
	// add the particles
	ParticleVector temp(decayProducts(daug,EvtPDL::getStdHep(daug->getId())!=92));
	if(EvtPDL::getStdHep(daug->getId())==92) {
	Lorentz5Momentum pstring(ThePEGMomentum(daug->getP4(),daug->mass()));
	Boost bv(-pstring.x()/pstring.e(),-pstring.y()/pstring.e(),-pstring.z()/pstring.e());
	for(unsigned int ix=0;ix<temp.size();++ix) {
	temp[ix]->deepBoost(bv);
	}
	}
	for(unsigned int iy=0;iy<temp.size();++iy) output.push_back(temp[iy]);
	}
	else if(!daug->isDecayed()) {
	EvtDecayBase *decayer = EvtDecayTable::getInstance()->getDecayFunc(daug);
	// must be a decayer
	if(!decayer) throw Exception() << "Could find EvtGen decayer in "
	<< "EvtGenInterface::decayProducts()"
	<< Exception::runerror;
	// If there are already daughters, then this step is already done!
	// figure out the masses
	if ( daug->getNDaug() == 0 ) daug->generateMassTree();
	// perform decay
	decayer->makeDecay(daug,false);
	// add the particles
	ParticleVector temp(decayProducts(daug));
	for(unsigned int iy=0;iy<temp.size();++iy) output.push_back(temp[iy]);
	}
	}
	else {
	throw Exception() << "Tried to convert an unknown particle, PDG code = " << id
	<< " from EvtGen in EvtGen::decayproducts which is "
	<< " unknown to Herwig++. Killing event."
	<< Exception::eventerror;
	}
	}
	// boost to lab
	if(output.size()>0) {
	if(boost) {
	Boost bv=ThePEGMomentum(part->getP4(),part->mass()).boostVector();
	for(unsigned int ix=0; ix<output.size(); ++ix) output[ix]->deepBoost(bv);
	}
	}
	return output;
	}

	void EvtGenInterface::ThePEGSpin(PPtr peg,EvtParticle *evt) const {
	// if no corresponding ThePEG spin type return
	if(evt->getSpinType()==EvtSpinType::STRING\|\|
	evt->getSpinType()==EvtSpinType::SPIN3\|\|evt->getSpinType()==EvtSpinType::SPIN4\|\|
	evt->getSpinType()==EvtSpinType::SPIN5HALF\|\|
	evt->getSpinType()==EvtSpinType::SPIN7HALF) return;
	// now deal with the other cases
	// scalars
	unsigned int ix;
	SpinPtr spin;
	if(evt->getSpinType()==EvtSpinType::SCALAR) {
	ScalarSpinPtr sca(new_ptr(ScalarSpinInfo(peg->momentum(),true)));
	spin=sca;
	}
	// massless spin 1/2
	else if(evt->getSpinType()==EvtSpinType::NEUTRINO) {
	FermionSpinPtr fer(new_ptr(FermionSpinInfo(peg->momentum(),true)));
	spin=fer;
	ix=peg->id()<0;
	fer->setBasisState(ix,ThePEGSpinor(evt->spParentNeutrino()));
	ix=peg->id()>0;
	fer->setBasisState(ix,LorentzSpinor<SqrtEnergy>());
	}
	// massive spin-1/2
	else if(evt->getSpinType()==EvtSpinType::DIRAC) {
	FermionSpinPtr fer(new_ptr(FermionSpinInfo(peg->momentum(),true)));
	spin=fer;
	for(ix=0;ix<2;++ix) {
	fer->setBasisState(ix,ThePEGSpinor(evt->spParent(ix)));
	}
	}
	// massless vectors
	else if(evt->getSpinType()==EvtSpinType::PHOTON) {
	VectorSpinPtr vec(new_ptr(VectorSpinInfo(peg->momentum(),true)));
	spin=vec;
	for(ix=0;ix<1;++ix) {
	vec->setBasisState(2*ix,ThePEGPolarization(evt->epsParentPhoton(ix)));}
	vec->setBasisState(1,LorentzVector<Complex>());
	}
	// massive vectors
	else if(evt->getSpinType()==EvtSpinType::VECTOR) {
	VectorSpinPtr vec(new_ptr(VectorSpinInfo(peg->momentum(),true)));
	spin=vec;
	for(ix=0;ix<3;++ix) {
	vec->setBasisState(ix,ThePEGPolarization(evt->epsParent(ix)));
	}
	}
	// massive spin 3/2
	else if(evt->getSpinType()==EvtSpinType::RARITASCHWINGER) {
	RSFermionSpinPtr rs(new_ptr(RSFermionSpinInfo(peg->momentum(),true)));
	spin=rs;
	for(ix=0;ix<4;++ix) {
	rs->setBasisState(ix,ThePEGRSSpinor(evt->spRSParent(ix)));}
	}
	// tensors
	else if(evt->getSpinType()==EvtSpinType::TENSOR) {
	TensorSpinPtr ten(new_ptr(TensorSpinInfo(peg->momentum(),true)));
	spin=ten;
	for(ix=0;ix<5;++ix) {
	ten->setBasisState(ix,ThePEGTensor(evt->epsTensorParent(ix)));
	}
	}
	else {
	throw Exception() << "Unknown EvtGen spin type in EvtGen::ThePEGSpin() "
	<< Exception::runerror;
	}
	// set the spininfo
	peg->spinInfo(spin);
	// final set up if the particle has decayed
	if(evt->getSpinDensityForward().getDim()!=0) {
	spin->rhoMatrix()=ThePEGSpinDensity(evt->getSpinDensityForward(),peg->id());
	}
	if(evt->isDecayed()) {
	spin->decayed(true);
	spin->develop();
	if(evt->getSpinDensityBackward().getDim()!=0) {
	spin->DMatrix()=ThePEGSpinDensity(evt->getSpinDensityBackward(),peg->id());
	}
	}
	}

	// convert from EvtGen to ThePEG
	int EvtGenInterface::ThePEGID(EvtId eid,bool exception) const {
	int output(0);
	int id(EvtPDL::getStdHep(eid)),absid(abs(id)),ispin(absid%10),isgn(id/absid);
	// handle the easy cases
	// special particles like string which we delete and include decay products
	if(absid==92\|\|absid==41\|\|absid==42\|\|absid==43\|\|absid==44\|\|
	absid==30343\|\|absid==30353\|\|
	absid==30363\|\|absid==30373\|\|absid==30383) {
	output=90;
	}
	//quarks(+excited)
	else if(absid<=8\|\|
	// leptons(+excited)
	(absid>=11&&absid<=18)\|\|
	// SM gauge bosons and higgs
	(absid>=21&&absid<=25)\|\|(absid>=32&&absid<=37)\|\|
	// 1 1S0, 1 3S1 and 1 3P2 mesons are the same
	(absid>100&&absid<600&&(ispin%2==1&&ispin<9))\|\|
	// 2 1S0 and 2 3S1 mesons are the same
	(absid>100100&&absid<100600&&(ispin==1\|\|ispin==3))\|\|
	// 1 3p0 are the same
	(absid>10100&&absid<10600&&ispin==1) \|\|
	// 3d1 are the same
	(absid>30100&&absid<30600&&ispin==3) \|\|
	// lowest baryon multiplets and diquarks
	(absid>1000&&absid<6000&&ispin>=1&&ispin<=4)\|\|
	// mixed kaons and diffractive states
	(absid>=100&&absid<=3000&&ispin==0)) {
	output=id;
	}
	// 1 1P1 mesons are the same apart from D_s1
	else if (absid>10100&&absid<10600&&(ispin==3)) {
	if(absid==10433) output=isgn*20433;
	else output=id;
	}
	// 1 3P1 mesons are the same apart from D_s1
	else if(absid>20100&&absid<20600&&(ispin==3)) {
	if(absid==20433) output=isgn*10433;
	else output=id;
	}
	// special for the virtual W (we don't distinguish so return W)
	else if(absid==89) return id/absid*24;
	// excited baryons
	else if(((absid>10000&&absid<16000)\|\|(absid>20000&&absid<26000)\|\|
	(absid>30000&&absid<36000))&&(ispin==2\|\|ispin==4)) {
	output=id;
	}
	// bottomium
	else if((absid%1000)/10==55) {
	output = id;
	}
	// special mesons
	else if(absid>9000000) {
	// f_0/a_0(980) goes into f_0/a_0(980)
	if(absid==9000111\|\|absid==9000211\|\|absid==9010221) {
	output=id;
	}
	// sigma
	else if(absid==9000221) {
	output=id;
	}
	// psi(4040)
	else if(absid==9000443) {
	output=id;
	}
	// f_0(1500)
	else if(absid==9020221) {
	output=9030221;
	}
	}
	// check its O.K.
	if(output==0&&exception) {
	throw Exception() << "Can't find the ThePEG id for particle "
	<< EvtPDL::name(eid)
	<< " with PDG code = "
	<< id << " in EvtGen::ThePEGID"
	<< Exception::eventerror;
	}
	return output;
	}

	RhoDMatrix EvtGenInterface::ThePEGSpinDensity(const EvtSpinDensity & rho, int id) const {
	RhoDMatrix output;
	unsigned int ix,iy;
	// special for neutrinos
	if(abs(id)%2==0&&abs(id)>=12&&abs(id)<=16) {
	output = RhoDMatrix(PDT::Spin1Half, false);
	if(id>0) output(0,0)=ThePEGComplex(rho.get(0,0));
	else output(1,1)=ThePEGComplex(rho.get(0,0));
	}
	// special for photons
	else if(id==ParticleID::gamma) {
	output = RhoDMatrix(PDT::Spin1, false);
	for(ix=0;ix<2;++ix) {
	for(iy=0;iy<2;++iy) output(2ix,2iy)=ThePEGComplex(rho.get(ix,iy));
	}
	}
	// normal behaviour
	else {
	unsigned int ndim(abs(rho.getDim()));
	if(ndim!=0) {
	output = RhoDMatrix(PDT::Spin(ndim));
	for(ix=0;ix<ndim;++ix) {
	for(iy=0;iy<ndim;++iy) output(ix,iy)=ThePEGComplex(rho.get(ix,iy));
	}
	}
	else if(ndim==0) {
	output = RhoDMatrix(PDT::Spin0);
	}
	}
	return output;
	}

	void EvtGenInterface::checkConversion() const {
	// check the translation of particles from ThePEG to EvtGen.
	ParticleMap::const_iterator pit = generator()->particles().begin();
	ParticleMap::const_iterator pend = generator()->particles().end();
	- generator()->log() << "Testing conversion of particles from ThePEG to EvtGen\n";
	+ cout << "Testing conversion of particles from ThePEG to EvtGen\n";
	EvtId etemp;
	for(;pit!=pend;++pit) {
	- generator()->log() << pit->first << " ";
	+ cout << pit->first << " ";
	etemp=EvtGenID(pit->first,false);
	Energy mass = pit->second->mass();
	Energy width = pit->second->width();
	if(etemp.getAlias()>=0) {
	- generator()->log() << pit->second->PDGName() << "\t becomes "
	+ cout << pit->second->PDGName() << "\t becomes "
	<< EvtPDL::name(etemp) << "\t "
	<< EvtPDL::getStdHep(etemp);
	if(ThePEGID(etemp,false)-pit->first!=0) {
	- generator()->log() << " and converting back to ThePEG fails";
	+ cout << " and converting back to ThePEG fails";
	}
	double mass2 = EvtPDL::getMeanMass(etemp);
	double width2 = EvtPDL::getWidth(etemp);
	if(mass!=ZERO) {
	double delta = (mass-mass2*GeV)/mass;
	if(abs(delta)>1e-6)
	- generator()->log() << " Mass Difference " << mass/GeV-mass2;
	+ cout << " Mass Difference " << mass/GeV-mass2;
	}
	if(width>ZERO) {
	double delta = (width-width2*GeV)/width;
	if(abs(delta)>1e-6)
	- generator()->log() << " Width Difference " << width/GeV-width2;
	+ cout << " Width Difference " << width/GeV-width2;
	}
	- generator()->log() << "\n";
	+ cout << "\n";
	}
	else {
	- generator()->log() << pit->second->PDGName()
	+ cout << pit->second->PDGName()
	<< " has no match in EvtGen \n";
	}
	}
	// test the conversion the other way
	- generator()->log() << "Testing conversion of particles from EvtGen to ThePEG\n";
	+ cout << "Testing conversion of particles from EvtGen to ThePEG\n";
	int idtemp;
	tcPDPtr ptemp;
	for(unsigned int ix=0;ix<EvtPDL::entries();++ix) {
	- generator()->log() << EvtPDL::getStdHep(EvtPDL::getEntry(ix)) << " "
	+ cout << EvtPDL::getStdHep(EvtPDL::getEntry(ix)) << " "
	<< EvtPDL::name(EvtPDL::getEntry(ix));
	idtemp=ThePEGID(EvtPDL::getEntry(ix),false);
	ptemp=getParticleData(idtemp);
	if(ptemp) {
	- generator()->log() << " becomes " << ptemp->PDGName() << " "
	+ cout << " becomes " << ptemp->PDGName() << " "
	<< ptemp->id();
	if(EvtGenID(ptemp->id(),false)!=EvtPDL::getEntry(ix)) {
	- generator()->log() << " and converting back to EvtGEN fails ";
	+ cout << " and converting back to EvtGEN fails ";
	}
	- generator()->log() << "\n";
	+ cout << "\n";
	}
	else {
	- generator()->log() << " has no match in ThePEG \n";
	+ cout << " has no match in ThePEG \n";
	}
	}
	}

	void EvtGenInterface::outputEvtGenDecays(long parentid) const {
	// output the decay modes from EvtGen so we can read them in
	EvtId parent(EvtGenID(parentid));
	// get evtids for children
	- generator()->log() << "Outputting decays for "
	+ cout << "Outputting decays for "
	<< getParticleData(parentid)->PDGName() << "\n";
	for(int ix=0;ix<EvtDecayTable::getInstance()->getNMode(parent.getAlias());++ix) {
	EvtDecayBase* decayer=EvtDecayTable::getInstance()->getDecay(parent.getAlias(),ix);
	vector<long> pegid;
	for(int iy=0;iy<decayer->getNDaug();++iy) {
	pegid.push_back(ThePEGID(decayer->getDaug(iy),false));
	}
	for(unsigned int iy=pegid.size();iy<7;++iy) pegid.push_back(0);
	double br=decayer->getBranchingFraction();
	- generator()->log()
	+ cout
	<< "insert into decay_modes (incomingID,BR,decayon,star,outgoingID1,"
	<< "outgoingID2,outgoingID3,outgoingID4,outgoingID5,outgoingID6,outgoingID7,"
	<< "description,decayer) values (" << parentid << "," << br << ",'on','*',";
	for(int iy=0;iy<7;++iy) {
	- generator()->log() << pegid[iy] << ",";
	+ cout << pegid[iy] << ",";
	}
	- generator()->log() << "'Decay of %name% with branching "
	- << "ratio taken from EvtGen.',3);\n";
	+ cout << "'Decay of %name% with branching "
	+ << "ratio taken from EvtGen.',3);\n";
	}
	}

	int EvtGenInterface::EvtGenChannel(const DecayMode & dm) const {
	// get evtid for parent
	EvtId parent(EvtGenID(dm.parent()->id()));
	// get evtids for children
	EvtId daugs[20];
	unsigned int ndaug(0);
	ParticleMSet::const_iterator pit = dm.products().begin();
	ParticleMSet::const_iterator pend = dm.products().end();
	EvtId etemp;
	for( ;pit!=pend&&ndaug<20;++pit) {
	etemp=EvtGenID((**pit).id());
	daugs[ndaug]=etemp;
	ndaug++;
	}
	if(ndaug>=20) throw Exception() << "Too many decay products "
	<< "in EvtGen::EvtGenChannel"
	<< Exception::eventerror;
	// find the channel
	int output=EvtDecayTable::getInstance()->inChannelList(parent,ndaug,daugs);
	if(output<0) {
	string dmode;
	dmode = "decay of " + dm.parent()->PDGName() + " -> ";
	pit = dm.products().begin();
	pend = dm.products().end();
	for( ;pit!=pend&&ndaug<20;++pit) dmode += (**pit).PDGName()+" ";
	throw Exception() << "Can't find EVtGen decay mode in EvtGenChannel for "
	<< dmode << " in EvtGen::EvtGenChannel"
	<< Exception::runerror;
	}
	return output;
	}

	// translate the matrix element to our form
	void EvtGenInterface::constructVertex(const Particle & parent,
	ParticleVector products,
	EvtDecayAmp* damp) const {
	unsigned int ix,iy;
	vector <PDT::Spin> outspin;
	for(ix=0;ix<products.size();++ix) {
	outspin.push_back(products[ix]->dataPtr()->iSpin());
	}
	vector<int> constants(products.size()+2,0);
	unsigned int nstate(1),idtemp;
	for(ix=outspin.size();ix>0;--ix) {
	idtemp=abs(products[ix-1]->id());
	if(idtemp==ParticleID::gamma) {
	nstate*=2;
	constants[ix]=nstate;
	}
	else if(idtemp%2==0&&idtemp>=12&&idtemp<=16) {
	constants[ix]=nstate;
	}
	else {
	nstate*=outspin[ix-1];
	constants[ix]=nstate;
	}
	}
	PDT::Spin inspin(parent.dataPtr()->iSpin());
	nstate*=inspin;
	constants[0]=nstate;
	constants[outspin.size()+1]=1;
	int eind[10]={0,0,0,0,0,0,0,0,0,0};
	unsigned int iloc;
	vector<unsigned int> hind(outspin.size()+1,0);
	DecayMEPtr newME = new_ptr(GeneralDecayMatrixElement(inspin,outspin));
	for(ix=0;ix<nstate;++ix) {
	iloc=0;
	hind[0]=ix/constants[1];
	if(inspin!=PDT::Spin0) {
	eind[iloc]=hind[0];
	++iloc;
	}
	for(iy=0;iy<outspin.size();++iy) {
	if(outspin[iy]==PDT::Spin0) {
	hind[iy+1]=0;
	}
	else if(products[iy]->id()==ParticleID::gamma) {
	hind[iy+1]=2*(ix%constants[iy+1])/constants[iy+2];
	eind[iloc]=hind[iy+1]/2;
	++iloc;
	}
	else if(constants[iy+1]==1) {
	hind[iy+1]=products[iy]->id()<0;
	}
	else {
	hind[iy+1]=(ix%constants[iy+1])/constants[iy+2];
	eind[iloc]=hind[iy+1];++iloc;
	}
	}
	if(iloc==0) eind[0]=0;
	(*newME)(hind)=ThePEGComplex(damp->amplitude().getAmp(eind));
	}
	// create the decay vertex
	VertexPtr vertex(new_ptr(DecayVertex()));
	DVertexPtr Dvertex(dynamic_ptr_cast<DVertexPtr>(vertex));
	// set the incoming particle for the decay vertex
	dynamic_ptr_cast<tcSpinPtr>(parent.spinInfo())->decayVertex(vertex);
	for(ix=0;ix<products.size();++ix) {
	dynamic_ptr_cast<tcSpinPtr>(products[ix]->spinInfo())
	->productionVertex(vertex);
	}
	// set the matrix element
	Dvertex->ME(newME);
	}
	diff --git a/Decay/EvtGen/EvtGenInterface.h b/Decay/EvtGen/EvtGenInterface.h
	--- a/Decay/EvtGen/EvtGenInterface.h
	+++ b/Decay/EvtGen/EvtGenInterface.h
	@@ -1,470 +1,486 @@
	// -- C++ --
	#ifndef Herwig_EvtGenInterface_H
	#define Herwig_EvtGenInterface_H
	//
	// This is the declaration of the EvtGenInterface class.
	//

	#include "EvtGenInterface.fh"
	#include "ThePEG/Interface/Interfaced.h"
	#include "ThePEG/Vectors/Lorentz5Vector.h"
	#include "ThePEG/Helicity/ScalarSpinInfo.h"
	#include "ThePEG/Helicity/FermionSpinInfo.h"
	#include "ThePEG/Helicity/VectorSpinInfo.h"
	#include "ThePEG/Helicity/RSFermionSpinInfo.h"
	#include "ThePEG/Helicity/TensorSpinInfo.h"
	#include "ThePEG/EventRecord/Particle.h"

	#include "EvtGenRandom.h"
	#include "EvtGen/EvtGen.hh"
	#include "EvtGenBase/EvtParticle.hh"
	#include "EvtGenBase/EvtSpinDensity.hh"
	#include "EvtGenBase/EvtVector4R.hh"
	#include "EvtGenBase/EvtVector4C.hh"
	#include "EvtGenBase/EvtTensor4C.hh"
	#include "EvtGenBase/EvtDiracSpinor.hh"
	#include "EvtGenBase/EvtRaritaSchwinger.hh"
	#include "EvtGenBase/EvtDecayAmp.hh"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* The EvtGenInterface class is the main class for the use of the EvtGen decay
	* package with Herwig.
	*
	* @see \ref EvtGenInterfaceInterfaces "The interfaces"
	* defined for EvtGenInterface.
	*/
	class EvtGenInterface: public Interfaced {

	public:

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

	/**
	+ * The copy constructor (explicit as cannot copy streams)
	+ */
	+ EvtGenInterface(const EvtGenInterface &);
	+
	+ /**
	* The destructor.
	*/
	virtual ~EvtGenInterface();
	//@}

	public:

	/**
	* Use EvtGen to perform a decay
	* @param parent The decaying particle
	* @param dm The decaymode
	* @return The decay products
	*/
	ParticleVector decay(const Particle &parent,
	bool recursive, const DecayMode & dm) const;

	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 Functions to convert between EvtGen and Herwig classes */
	//@{
	/**
	* Convert a particle to an EvtGen particle.
	* @param part The particle to be converted.
	*/
	EvtParticle *EvtGenParticle(const Particle & part) const;

	/**
	* Return the decay products of an EvtGen particle in as ThePEG particles
	* @param evtpart The EvtGen particle
	*/
	ParticleVector decayProducts(EvtParticle* evtpart, bool boost=true) const;

	/**
	* Convert a Lorentz5Momentum to a real EvtGen 4-vector
	* @param mom The momentum to be converted
	*/
	EvtVector4R EvtGenMomentum(const Lorentz5Momentum & mom) const {
	return EvtVector4R(mom.t()/GeV,mom.x()/GeV,mom.y()/GeV,mom.z()/GeV);
	}

	/**
	* Convert a PDG code from ThePEG into an EvtGen particle id
	* @param id The PDG code
	* @param exception Whether or not to throw an Exception if fails
	*/
	EvtId EvtGenID(int id,bool exception=true) const;

	/**
	* Convert a LorentzSpinor to an EvtGen one. The spinor is converted to the
	* EvtGen Dirac representation/
	* @param sp The LorentzSpinor
	*/
	EvtDiracSpinor EvtGenSpinor(const LorentzSpinor<SqrtEnergy> & sp) const {
	InvSqrtEnergy norm(sqrt(0.5)/sqrt(GeV));
	EvtDiracSpinor output;
	output.set(EvtGenComplex(-norm*( sp.s1()+sp.s3())),
	EvtGenComplex(-norm*( sp.s2()+sp.s4())),
	EvtGenComplex(-norm*(-sp.s1()+sp.s3())),
	EvtGenComplex(-norm*(-sp.s2()+sp.s4())));
	return output;
	}

	/**
	* Convert a LorentzPolarizationVector to a complex EvtGen 4-vector
	* @param eps The polarization vector to be converted
	*/
	EvtVector4C EvtGenPolarization(const LorentzPolarizationVector & eps) const {
	return EvtVector4C(EvtGenComplex(eps.t()),EvtGenComplex(eps.x()),
	EvtGenComplex(eps.y()),EvtGenComplex(eps.z()));
	}

	/**
	* Convert our Rarita-Schwinger spinor to the EvtGen one
	* @param sp Our RS Spinor
	*/
	EvtRaritaSchwinger EvtGenRSSpinor(const LorentzRSSpinor<SqrtEnergy> & sp) const {
	InvSqrtEnergy norm(sqrt(0.5)/sqrt(GeV));
	complex<double> out[4][4];
	for(unsigned int ix=0;ix<4;++ix) {
	out[ix][0] =-norm*( sp(ix,0)+sp(ix,2));
	out[ix][1] =-norm*( sp(ix,1)+sp(ix,3));
	out[ix][2] =-norm*(-sp(ix,0)+sp(ix,2));
	out[ix][3] =-norm*(-sp(ix,1)+sp(ix,3));
	}
	EvtRaritaSchwinger output;
	unsigned int ix,iy;
	// remember we have vec,spin and evtgen spin,vec
	for(ix=0;ix<4;++ix) {
	for(iy=0;iy<4;++iy) output.set(ix,iy,EvtGenComplex(out[iy][ix]));
	}
	return output;
	}

	/**
	* Convert our tensor to the EvtGen one.
	* @param ten Our tensor
	*/
	EvtTensor4C EvtGenTensor(const LorentzTensor<double> & ten) const {
	EvtTensor4C output;
	unsigned int ix,iy;
	for(ix=0;ix<4;++ix){
	for(iy=0;iy<4;++iy) output.set(ix,iy,EvtGenComplex(ten(ix,iy)));
	}
	return output;
	}

	/**
	* Convert a spin density matrix to an EvtGen spin density matrix.
	* @param rho The spin density matrix to be converted.
	*/
	EvtSpinDensity EvtGenSpinDensity(const RhoDMatrix & rho) const {
	EvtSpinDensity rhoout;
	unsigned int ix,iy,ispin(rho.iSpin());
	rhoout.setDim(ispin);
	for(ix=0;ix<ispin;++ix) {
	for(iy=0;iy<ispin;++iy)
	rhoout.set(ix,iy,EvtGenComplex(rho(ix,iy)));
	}
	return rhoout;
	}

	/**
	* Convert from our complex to the EvtGen one
	*/
	EvtComplex EvtGenComplex(Complex c) const {
	return EvtComplex(c.real(),c.imag());
	}
	//@}

	/**
	* Functions to convert between EvtGen and Herwig classes
	*/
	//@{
	/**
	* Convert a particle from an EvtGen one to ThePEG one.
	* @param part The EvtGen particle.
	* @param pd Pointer to the particle data object of ThePEG for the particle.
	* @param spin Convert the spin information as well
	*/
	PPtr ThePEGParticle(EvtParticle *part, tcPDPtr pd,bool spin=true) const {
	PPtr output(new_ptr(Particle(pd)));
	output->set5Momentum(ThePEGMomentum(part->getP4(),part->mass()));
	if(spin) ThePEGSpin(output,part);
	// make the daughters
	ParticleVector daug(decayProducts(part,spin));
	for(unsigned int ix=0;ix<daug.size();++ix) output->addChild(daug[ix]);
	return output;
	}

	/**
	* Set the SpinInfo for a ThePEG particle using an EvtGen particle
	* @param pegpart ThePEG particle.
	* @param evtpart The EvtGen particle.
	*/
	void ThePEGSpin(PPtr pegpart,EvtParticle *evtpart) const;

	/**
	* Convert an EvtGen EvtId to a PDG code in our conventions
	* @param id The EvtGen ID.
	* @param exception Whether or not to throw an Exception if fails
	*/
	int ThePEGID(EvtId id,bool exception=true) const;

	/**
	* Convert from EvtGen momentum to Lorentz5Momentum
	* @param mom The EvtGen 4-momentum
	* @param mass The mass
	*/
	Lorentz5Momentum ThePEGMomentum(const EvtVector4R & mom,double mass) const {
	return Lorentz5Momentum(mom.get(1)GeV,mom.get(2)GeV,
	mom.get(3)GeV,mom.get(0)GeV,mass*GeV);
	}
	/**
	* Convert from EvtGen complex to ours
	*/
	Complex ThePEGComplex(EvtComplex c) const {
	return Complex(real(c),imag(c));
	}

	/**
	* Convert a spin density to a ThePEG one from an EvtGen one
	* @param rho The spin density matrix to be converted
	* @param id The PDG code of the particle to get special cases right.
	*/
	RhoDMatrix ThePEGSpinDensity(const EvtSpinDensity & rho, int id) const;

	/**
	* Convert an EvtDiracSpinor a LorentzSpinor. This spinor is converted to
	* the default Dirac matrix representation used by ThePEG.
	* @param sp The EvtDiracSpinor
	*/
	LorentzSpinor<SqrtEnergy> ThePEGSpinor(const EvtDiracSpinor & sp) const {
	SqrtEnergy norm(sqrt(0.5)*sqrt(GeV));
	vector<complex<SqrtEnergy> > evtSpin(4);
	for(unsigned int ix=0;ix<4;++ix) evtSpin[ix] = -norm*ThePEGComplex(sp.get_spinor(ix));
	return LorentzSpinor<SqrtEnergy>(evtSpin[0]-evtSpin[2],evtSpin[1]-evtSpin[3],
	evtSpin[0]+evtSpin[2],evtSpin[1]+evtSpin[3]);
	}

	/**
	* Convert an EvtGen complex 4-vector to a LorentzPolarizationVector
	* @param eps The complex 4-vector to be converted.
	*/
	LorentzPolarizationVector ThePEGPolarization(const EvtVector4C & eps) const {
	return LorentzPolarizationVector(conj(ThePEGComplex(eps.get(1))),
	conj(ThePEGComplex(eps.get(2))),
	conj(ThePEGComplex(eps.get(3))),
	conj(ThePEGComplex(eps.get(0))));
	}

	/**
	* Convert an EvtGen Rarita-Schwinger spinor to ours
	* @param sp The EvtGen RS spinor.
	*/
	LorentzRSSpinor<SqrtEnergy> ThePEGRSSpinor(const EvtRaritaSchwinger & sp) const {
	complex<SqrtEnergy> evtSpin[4][4];
	SqrtEnergy norm(sqrt(0.5)*sqrt(GeV));
	// normalisation and swap vec,spin order
	for(unsigned int ix=0;ix<4;++ix) {
	for(unsigned int iy=0;iy<4;++iy) evtSpin[ix][iy]=-norm*ThePEGComplex(sp.get(iy,ix));
	}
	LorentzRSSpinor<SqrtEnergy> output;
	for(unsigned int ix=0;ix<4;++ix) {
	output(ix,0) = evtSpin[ix][0] - evtSpin[ix][2];
	output(ix,1) = evtSpin[ix][1] - evtSpin[ix][3];
	output(ix,2) = evtSpin[ix][0] + evtSpin[ix][2];
	output(ix,3) = evtSpin[ix][1] + evtSpin[ix][3];
	}
	// output.changeRep(Helicity::defaultDRep);
	return output;
	}

	/**
	* Convert an EvtGen tensor to ThePEG
	* @param ten The EvtGen tensor
	*/
	LorentzTensor<double> ThePEGTensor(const EvtTensor4C & ten) const {
	LorentzTensor<double> output;
	unsigned int ix,iy;
	for(ix=0;ix<4;++ix) {
	for(iy=0;iy<4;++iy)output(ix,iy)=conj(ThePEGComplex(ten.get(ix,iy)));
	}
	return output;
	}
	//@}

	/**
	* Check the conversion of particles between Herwig++ and EvtGen
	*/
	void checkConversion() const;

	/**
	* Output the EvtGen decay modes for a given particle
	* @param id The PDG code of the particle to output
	*/
	void outputEvtGenDecays(long id) const;

	/**
	* Find the location in the EvtGen list of decay channels for
	* a given decay mode.
	*/
	int EvtGenChannel(const DecayMode &dm) const;

	/**
	* Check the particle has SpinInfo and if not create it
	* @param part The particle
	*/
	tSpinPtr getSpinInfo(const Particle &part) const {
	// return spin info if exists
	if(part.spinInfo()) {
	return dynamic_ptr_cast<tSpinPtr>(const_ptr_cast<tPPtr>(&part)->spinInfo());
	}
	// otherwise make it
	tPPtr ptemp(const_ptr_cast<tPPtr>(&part));
	PDT::Spin spin(part.dataPtr()->iSpin());
	SpinPtr pspin;
	if(spin==PDT::Spin0) pspin=new_ptr(ScalarSpinInfo(part.momentum(),true));
	else if(spin==PDT::Spin1Half) pspin=new_ptr(FermionSpinInfo(part.momentum(),true));
	else if(spin==PDT::Spin1) pspin=new_ptr(VectorSpinInfo(part.momentum(),true));
	else if(spin==PDT::Spin3Half) pspin=new_ptr(RSFermionSpinInfo(part.momentum(),true));
	else if(spin==PDT::Spin2) pspin=new_ptr(TensorSpinInfo(part.momentum(),true));
	else throw Exception() << "Can't create spinInfo for decaying particle in "
	<< "EvtGen::checkSpinInfo for spin " << spin << "particle "
	<< Exception::eventerror;
	ptemp->spinInfo(pspin);
	return pspin;
	}

	/**
	* Construct the DecayVertex for Herwig using the information from
	* EvtGen
	* @param parent The decaying particle
	* @param products The outgoing particles
	* @param damp Pointer to the EvtGen decayer
	*/
	void constructVertex(const Particle & parent,ParticleVector products,
	EvtDecayAmp* damp) const;

	protected:

	/** @name Clone Methods. */
	//@{
	/**
	* Make a simple clone of this object.
	* @return a pointer to the new object.
	*/
	virtual IBPtr clone() const;

	/** Make a clone of this object, possibly modifying the cloned object
	* to make it sane.
	* @return a pointer to the new object.
	*/
	virtual IBPtr fullclone() const;
	//@}

	protected:

	/** @name Standard Interfaced functions. */
	//@{
	/**
	* Initialize this object. Called in the run phase just before
	* a run begins.
	*/
	virtual void doinitrun();
	+
	+ /**
	+ * Finalize this object. Called in the run phase just after a
	+ * run has ended. Used eg. to write out statistics.
	+ */
	+ virtual void dofinish();
	//@}

	private:

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

	private:

	/**
	* Names of the various EvtGen parameter files
	*/
	//@{
	/**
	* The name of the file containing the decays
	*/
	string decayName_;

	/**
	* The name of the file containing the particle data
	*/
	string pdtName_;

	/**
	* Names of addition user specified decays
	*/
	vector<string> userDecays_;
	//@}

	/**
	* Whether or not to redirect cout and cerr when EvtGen is running
	*/
	bool reDirect_;

	/**
	* Check the conversion of the particles
	*/
	bool checkConv_;

	/**
	* Particles for which to output the EvtGen decays
	*/
	vector<long> convID_;

	/**
	* Location of the PYTHIA8 data directory
	*/
	string p8Data_;

	private:

	/**
	* Pointer to the random number generator for EvtGen
	*/
	EvtRandomEngine * evtrnd_;

	/**
	* Main EvtGen object
	*/
	EvtGen * evtgen_;

	+ /**
	+ * File to output the log info to
	+ */
	+ mutable ofstream logFile_;
	+
	};

	}

	#endif /* Herwig_EvtGenInterface_H */
	diff --git a/Shower/Dipole/Kinematics/FIMassiveKinematics.cc b/Shower/Dipole/Kinematics/FIMassiveKinematics.cc
	--- a/Shower/Dipole/Kinematics/FIMassiveKinematics.cc
	+++ b/Shower/Dipole/Kinematics/FIMassiveKinematics.cc
	@@ -1,301 +1,301 @@
	// -- C++ --
	//
	// FIMassiveKinematics.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2007 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 FIMassiveKinematics class.
	//

	#include "FIMassiveKinematics.h"
	#include "ThePEG/Interface/ClassDocumentation.h"


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

	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "Herwig/Shower/Dipole/Base/DipoleSplittingInfo.h"
	#include "Herwig/Shower/Dipole/Kernels/DipoleSplittingKernel.h"

	using namespace Herwig;

	FIMassiveKinematics::FIMassiveKinematics()
	: DipoleSplittingKinematics() {}

	FIMassiveKinematics::~FIMassiveKinematics() {}

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

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

	pair<double,double> FIMassiveKinematics::kappaSupport(const DipoleSplittingInfo&) const {
	return {0.0,1.0};
	}

	pair<double,double> FIMassiveKinematics::xiSupport(const DipoleSplittingInfo& split) const {
	double c = sqrt(1.-4.*sqr(IRCutoff()/generator()->maximumCMEnergy()));
	if ( split.index().emitterData()->id() == ParticleID::g ) {
	if ( split.emissionData()->id() != ParticleID::g )
	return {0.5(1.-c),0.5(1.+c)};
	double b = log((1.+c)/(1.-c));
	return {-b,b};
	}
	return {-log(0.5(1.+c)),-log(0.5(1.-c))};
	}

	// sbar
	Energy FIMassiveKinematics::dipoleScale(const Lorentz5Momentum& pEmitter,
	const Lorentz5Momentum& pSpectator) const {
	return sqrt(2.(pEmitterpSpectator));
	}

	Energy FIMassiveKinematics::ptMax(Energy dScale,
	double, double specX,
	const DipoleIndex& ind,
	const DipoleSplittingKernel& split) const {
	Energy mi = split.emitter(ind)->mass(), m = split.emission(ind)->mass();
	Energy2 mi2 = sqr(mi), m2 = sqr(m);
	// Energy2 Mi2 = split.emitter(int)->id() + split.emission(int)->id() == 0 ?
	// 0.*GeV2 : mi2;
	Energy2 Mi2 = mi2 == m2 ? 0.*GeV2 : mi2;

	// s^star/x
	Energy2 s = sqr(dScale) * (1.-specX)/specX + Mi2;
	return .5 * sqrt(s) * rootOfKallen( s/s, mi2/s, m2/s );
	}

	// what is this? in FF it is given by y+dScale = sqrt( 2qiq / bar )->max
	Energy FIMassiveKinematics::QMax(Energy dScale,
	double, double specX,
	const DipoleIndex&,
	const DipoleSplittingKernel&) const {
	generator()->log() << "FIMassiveKinematics::QMax called.\n" << flush;
	assert(false && "implementation missing");
	// this is sqrt( 2qi*q ) -> max;
	return dScale * sqrt((1.-specX)/specX);
	}

	Energy FIMassiveKinematics::PtFromQ(Energy scale, const DipoleSplittingInfo& split) const {
	// from Martin's thesis
	double z = split.lastZ();
	Energy mi = split.emitterData()->mass();
	Energy m = split.emissionData()->mass();
	Energy2 pt2 = z(1.-z)sqr(scale) - (1-z)sqr(mi) - zsqr(m);
	assert(pt2 >= ZERO);
	return sqrt(pt2);
	}

	Energy FIMassiveKinematics::QFromPt(Energy pt, const DipoleSplittingInfo& split) const {
	// from Martin's thesis
	double z = split.lastZ();
	Energy mi = split.emitterData()->mass();
	Energy m = split.emissionData()->mass();
	Energy2 Q2 = (sqr(pt) + (1-z)sqr(mi) + zsqr(m))/(z*(1.-z));
	return sqrt(Q2);
	}

	double FIMassiveKinematics::ptToRandom(Energy pt, Energy,
	double,double,
	const DipoleIndex&,
	const DipoleSplittingKernel&) const {
	return log(pt/IRCutoff()) / log(0.5 * generator()->maximumCMEnergy()/IRCutoff());
	}

	bool FIMassiveKinematics::generateSplitting(double kappa, double xi, double rphi,
	DipoleSplittingInfo& info,
	const DipoleSplittingKernel&) {

	if ( info.spectatorX() < xMin() ) {
	jacobian(0.0);
	return false;
	}

	Energy pt = IRCutoff() * pow(0.5 * generator()->maximumCMEnergy()/IRCutoff(),kappa);

	if ( pt > info.hardPt() \|\| pt < IRCutoff() ) {
	jacobian(0.0);
	return false;
	}

	double z;
	double mapZJacobian;

	if ( info.index().emitterData()->id() == ParticleID::g ) {
	if ( info.emissionData()->id() != ParticleID::g ) {
	z = xi;
	mapZJacobian = 1.;
	} else {
	z = exp(xi)/(1.+exp(xi));
	mapZJacobian = z*(1.-z);
	}
	} else {
	z = 1.-exp(-xi);
	mapZJacobian = 1.-z;
	}

	// double s = z*(1.-z);

	// double xs = info.spectatorX();

	// double x = 1. / ( 1. + sqr(pt/info.scale()) / s );

	// double zp = 0.5*(1.+sqrt(1.-sqr(pt/info.hardPt())));
	// double zm = 0.5*(1.-sqrt(1.-sqr(pt/info.hardPt())));

	Energy2 mi2 = sqr(info.emitterData()->mass());
	Energy2 m2 = sqr(info.emissionData()->mass());
	Energy2 Mi2 = info.emitterData()->id()+info.emissionData()->id() == 0 ?
	0.*GeV2 : mi2;

	// s^star/x
	Energy2 s = sqr(info.scale()) * (1.-info.spectatorX())/info.spectatorX() + Mi2;

	double xs = info.spectatorX();
	double x = 1. / ( 1. +
	( sqr(pt) + (1.-z)mi2 + zm2 - z(1.-z)Mi2 ) /
	( z(1.-z)s ) );



	Energy hard=info.hardPt();


	if(openZBoundaries()==1){
	hard=.5 * sqrt(s) * rootOfKallen( s/s, mi2/s, m2/s );
	}
	Energy2 sdip = sqr(info.scale()) + Mi2;
	if(openZBoundaries()==2){
	hard=min(0.5sqrt(s)
	rootOfKallen( s/s, mi2/s, m2/s ) ,
	sqrt(sdip) *
	rootOfKallen( sdip/sdip, mi2/sdip, m2/sdip ));
	}


	double ptRatio = sqrt(1.-sqr(pt/info.hardPt()));




	double zm1 = .5( 1.+(mi2-m2)/s - rootOfKallen(s/s,mi2/s,m2/s) ptRatio);
	double zp1 = .5( 1.+(mi2-m2)/s + rootOfKallen(s/s,mi2/s,m2/s) ptRatio);

	if ( // pt < IRCutoff() \|\|
	// pt > info.hardPt() \|\|
	z > zp1 \|\| z < zm1 \|\|
	x < xs ) {
	jacobian(0.0);
	return false;
	}

	// additional purely kinematic constraints
	double mui2 = x*mi2/sqr(info.scale());
	double mu2 = x*m2/sqr(info.scale());
	double Mui2 = x*Mi2/sqr(info.scale());
	double xp = 1. + Mui2 - sqr(sqrt(mui2)+sqrt(mu2));
	double root = sqr(1.-x+Mui2-mui2-mu2)-4.mui2mu2;
	if( root < 0. && root>-1e-10 )
	root = 0.;
	else if (root <0. ) {
	jacobian(0.0);
	return false;
	}

	root = sqrt(root);
	double zm = .5*( 1.-x+Mui2+mui2-mui2 - root ) / (1.-x+Mui2);
	double zp = .5*( 1.-x+Mui2+mui2-mui2 + root ) / (1.-x+Mui2);
	if (x > xp \|\|
	z > zp \|\| z < zm ) {
	jacobian(0.0);
	return false;
	}

	double phi = 2.Constants::pirphi;

	-// Compute and store the jacobian
	+ // Compute and store the jacobian
	Energy2 pt2 = sqr(pt);
	- double jacPt2 = 1. / ( 1. + sqr(1.-z)mi2/pt2 + zz*m2/pt2 );
	+ double jacPt2 = 1. / ( 1. + (1.-z)mi2/pt2 + zm2/pt2 - z(1.-z)Mi2/pt2 );
	jacobian( jacPt2 * mapZJacobian * 2.log(0.5 generator()->maximumCMEnergy()/IRCutoff()));

	lastPt(pt);
	lastZ(z);
	lastPhi(phi);
	lastSpectatorZ(x);

	if ( theMCCheck )
	theMCCheck->book(1.,info.spectatorX(),info.scale(),info.hardPt(),pt,z,jacobian());

	return true;

	}

	void FIMassiveKinematics::generateKinematics(const Lorentz5Momentum& pEmitter,
	const Lorentz5Momentum& pSpectator,
	const DipoleSplittingInfo& dInfo) {

	Energy pt = dInfo.lastPt();
	double z = dInfo.lastZ();

	Lorentz5Momentum kt =
	getKt (pSpectator, pEmitter, pt, dInfo.lastPhi(),true);

	Energy2 mi2 = sqr(dInfo.emitterData()->mass());
	Energy2 m2 = sqr(dInfo.emissionData()->mass());
	Energy2 Mi2 = dInfo.emitterData()->id() + dInfo.emissionData()->id() == 0 ?
	0.*GeV2 : mi2;

	double xInv = ( 1. +
	(ptpt+(1.-z)mi2+zm2-z(1.-z)*Mi2) /
	(z(1.-z)sqr(dInfo.scale())) );

	Lorentz5Momentum em = z*pEmitter +
	(sqr(pt)+mi2-zzMi2)/(zsqr(dInfo.scale()))pSpectator + kt;
	em.setMass(sqrt(mi2));
	em.rescaleEnergy();

	Lorentz5Momentum emm = (1.-z)*pEmitter +
	(ptpt+m2-sqr(1.-z)Mi2)/((1.-z)sqr(dInfo.scale()))pSpectator - kt;
	emm.setMass(sqrt(m2));
	emm.rescaleEnergy();

	Lorentz5Momentum spe = xInv*pSpectator;
	spe.setMass(ZERO);
	spe.rescaleEnergy();

	emitterMomentum(em);
	emissionMomentum(emm);
	spectatorMomentum(spe);

	}

	// If needed, insert default implementations of function defined
	// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).


	void FIMassiveKinematics::persistentOutput(PersistentOStream & ) const {
	}

	void FIMassiveKinematics::persistentInput(PersistentIStream & , int) {
	}

	ClassDescription<FIMassiveKinematics> FIMassiveKinematics::initFIMassiveKinematics;
	// Definition of the static class description member.

	void FIMassiveKinematics::Init() {

	static ClassDocumentation<FIMassiveKinematics> documentation
	("FIMassiveKinematics implements massless splittings "
	"off a final-initial dipole.");

	}

File Metadata

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

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions