No OneTemporary
Actions

Size

105 KB

Subscribers

None

View Options

	diff --git a/Decay/FormFactors/KMatrix.cc b/Decay/FormFactors/KMatrix.cc
	--- a/Decay/FormFactors/KMatrix.cc
	+++ b/Decay/FormFactors/KMatrix.cc
	@@ -1,127 +1,133 @@
	// -- C++ --
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the KMatrix class.
	//

	#include "KMatrix.h"
	#include "ThePEG/Interface/ClassDocumentation.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 <boost/numeric/ublas/banded.hpp>
	#include <boost/numeric/ublas/lu.hpp>

	using namespace Herwig;

	KMatrix::KMatrix(FlavourInfo flavour, vector<Channels> channels,
	vector<Energy2> poles) : flavour_(flavour), channels_(channels), poles_(poles)
	{}

	void KMatrix::persistentOutput(PersistentOStream & os) const {
	os << ounit(poles_,GeV2) << ounit(mPiPlus_,GeV) << ounit(mPi0_,GeV)
	<< ounit(mKPlus_,GeV) << ounit(mK0_,GeV) << ounit(mEta_,GeV)
	<< ounit(mEtaPrime_,GeV);
	}

	void KMatrix::persistentInput(PersistentIStream & is, int) {
	is >> iunit(poles_,GeV2) >> iunit(mPiPlus_,GeV) >> iunit(mPi0_,GeV)
	>> iunit(mKPlus_,GeV) >> iunit(mK0_,GeV) >> iunit(mEta_,GeV)
	>> iunit(mEtaPrime_,GeV);
	}


	// The following static variable is needed for the type
	// description system in ThePEG.
	DescribeAbstractClass<KMatrix,Interfaced>
	describeHerwigKMatrix("Herwig::KMatrix", "Herwig.so");

	void KMatrix::Init() {

	static ClassDocumentation<KMatrix> documentation
	("The KMatrix class provides a base class for the implementation of "
	"K-matrix parameterizations in Herwig");

	}

	void KMatrix::doinit() {
	Interfaced::doinit();
	// The charged pion mass
	mPiPlus_=getParticleData(ParticleID::piplus)->mass();
	// The neutral pion mass
	mPi0_=getParticleData(ParticleID::pi0)->mass();
	// The charged kaon mass
	mKPlus_=getParticleData(ParticleID::Kplus)->mass();
	// The neutral kaon mass
	mK0_=getParticleData(ParticleID::K0)->mass();
	// The eta mass
	mEta_=getParticleData(ParticleID::eta)->mass();
	// The eta' mass
	mEtaPrime_=getParticleData(ParticleID::etaprime)->mass();
	}

	namespace {

	double kallen(const Energy2 &s, const Energy &m1, const Energy & m2) {
	return (1.-sqr(m1+m2)/s)*(1.-sqr(m1-m2)/s);
	}
	}

	ublas::matrix<Complex> KMatrix::rho(Energy2 s) const {
	size_t msize = channels_.size();
	ublas::diagonal_matrix<Complex> rho(msize,msize);
	for(unsigned int iChan=0;iChan<msize;++iChan) {
	double val(0);
	switch (channels_[iChan]) {
	case PiPi:
	val=kallen(s,mPiPlus_,mPiPlus_);
	break;
	case KPi:
	val=kallen(s,mKPlus_,mPiPlus_);
	break;
	case KEta:
	val=kallen(s,mKPlus_,mEta_);
	break;
	case KEtaPrime:
	val=kallen(s,mKPlus_,mEtaPrime_);
	break;
	default:
	assert(false);
	}
	if(val>=0)
	rho(iChan,iChan) = sqrt(val);
	else
	rho(iChan,iChan) = Complex(0.,1.)*sqrt(-val);
	}
	return rho;
	}

	ublas::vector<Complex> KMatrix::
	amplitudes(Energy2 s, ublas::vector<Complex> pVector, bool multiplyByPoles) const {
	static const Complex ii(0.,1.);
	const ublas::identity_matrix<Complex> I(channels_.size());
	double fact(1.);
	if(multiplyByPoles) {
	for (Energy2 pole : poles_ ) fact *= 1.-s/pole;
	}
	// matrix for which we need the inverse
	ublas::matrix<Complex> m = factI-iiprod(K(s,multiplyByPoles),rho(s));
	- // create a permutation matrix for the LU-factorization
	- ublas::permutation_matrix<std::size_t> pm(m.size1());
	- // perform LU-factorization
	- int res = ublas::lu_factorize(m,pm);
	- if( res != 0 ) {
	- cerr << "problem with factorization\n";
	- exit(1);
	+ // inverse matrix
	+ ublas::matrix<Complex> inverse = ublas::identity_matrix<Complex>(m.size1());
	+ // just a number
	+ if(m.size1()==1) {
	+ inverse(0,0) = 1./m(0,0);
	}
	- // matrix for the output
	- ublas::matrix<Complex> inverse = ublas::identity_matrix<Complex>(m.size1());
	- // backsubstitute to get the inverse
	- ublas::lu_substitute(m, pm, inverse);
	+ else {
	+ // create a permutation matrix for the LU-factorization
	+ ublas::permutation_matrix<std::size_t> pm(m.size1());
	+ // perform LU-factorization
	+ int res = ublas::lu_factorize(m,pm);
	+ if( res != 0 ) {
	+ cerr << "problem with factorization\n";
	+ exit(1);
	+ }
	+ // backsubstitute to get the inverse
	+ ublas::lu_substitute(m, pm, inverse);
	+ }
	// compute the amplitudes
	return prod(inverse,pVector);
	}
	diff --git a/Decay/ScalarMeson/DtoKPiPiE691.cc b/Decay/ScalarMeson/DtoKPiPiE691.cc
	--- a/Decay/ScalarMeson/DtoKPiPiE691.cc
	+++ b/Decay/ScalarMeson/DtoKPiPiE691.cc
	@@ -1,705 +1,705 @@
	// -- C++ --
	//
	// DtoKPiPiE691.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2019 The Herwig Collaboration
	//
	// Herwig is licenced under version 3 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 DtoKPiPiE691 class.
	//

	#include "DtoKPiPiE691.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/ParVector.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h"
	#include "Herwig/Decay/GeneralDecayMatrixElement.h"

	using namespace Herwig;

	using ThePEG::Helicity::ScalarWaveFunction;
	using ThePEG::Helicity::incoming;
	using ThePEG::Helicity::outgoing;

	DtoKPiPiE691::DtoKPiPiE691() {
	// amplitudes and phases for D+ -> K-pi+pi+
	_a1NR = 1. ; _phi1NR = 0.;
	_a1K892 = 0.78; _phi1K892 = - 60.;
	_a1K1430 = 0.53; _phi1K1430 = 132.;
	_a1K1680 = 0.47; _phi1K1680 = - 51.;
	// amplitudes and phases for D0 -> K-pi+pi0
	_a2NR = 1.00; _phi2NR = 0.;
	_a2K8920 = 3.19; _phi2K8920 = 167.;
	_a2K892m = 2.96; _phi2K892m = -112.;
	_a2rho = 8.56; _phi2rho = 40.;
	// amplitudes and phases for D0 -> Kbar0 pi+pi-
	_a3NR = 1.00; _phi3NR = 0.;
	_a3K892 = 2.31; _phi3K892 = 109.;
	_a3rho = 1.59; _phi3rho = -123.;
	// masses and widths
	_localparameters=true;
	_mK8920 = 0.8961 GeV; _wK8920 = 0.0505GeV;
	_mK892m = 0.89159GeV; _wK892m = 0.0498GeV;
	_mK1680 = 1.714 GeV; _wK1680 = 0.323 GeV;
	_mK1430 = 1.429 GeV; _wK1430 = 0.287 GeV;
	_mrho0 = 0.7681 GeV; _wrho0 = 0.1515GeV;
	_mrhop = 0.7681 GeV; _wrhop = 0.1515GeV;
	// generate the intermediate particles
	generateIntermediates(true);
	}

	void DtoKPiPiE691::doinit() {
	DecayIntegrator::doinit();
	// complex amplitudes calculated from magnitudes and phases
	double fact = Constants::pi/180.;
	// D+ -> K-pi+pi+
	_c1NR = _a1NR Complex(cos(_phi1NR fact),sin(_phi1NR *fact));
	_c1K892 = _a1K892 Complex(cos(_phi1K892 fact),sin(_phi1K892 *fact));
	_c1K1430= _a1K1430Complex(cos(_phi1K1430fact),sin(_phi1K1430*fact));
	_c1K1680= _a1K1680Complex(cos(_phi1K1680fact),sin(_phi1K1680*fact));
	// D0 -> K-pi+pi0
	_c2NR = _a2NR Complex(cos(_phi2NR fact),sin(_phi2NR *fact));
	_c2K8920= _a2K8920Complex(cos(_phi2K8920fact),sin(_phi2K8920*fact));
	_c2K892m= _a2K892mComplex(cos(_phi2K892mfact),sin(_phi2K892m*fact));
	_c2rho = _a2rho Complex(cos(_phi2rho fact),sin(_phi2rho *fact));
	// D0 -> Kbar0pi+pi-
	_c3NR = _a3NR Complex(cos(_phi3NR fact),sin(_phi3NR *fact));
	_c3K892 = _a3K892 Complex(cos(_phi3K892 fact),sin(_phi3K892 *fact));
	_c3rho = _a3rho Complex(cos(_phi3rho fact),sin(_phi3rho *fact));
	// intermediate resonances
	tPDPtr k8920 = getParticleData(ParticleID::Kstarbar0 );
	tPDPtr k892m = getParticleData(ParticleID::Kstarminus);
	tPDPtr k1430 = getParticleData(ParticleID::Kstar_0bar0);
	tPDPtr k1680 = getParticleData(-30313);
	tPDPtr rho0 = getParticleData(ParticleID::rho0);
	tPDPtr rhop = getParticleData(ParticleID::rhoplus);
	// D+ -> K-pi+pi+
	tPDPtr in = getParticleData(ParticleID::Dplus);
	tPDVector out= {getParticleData(ParticleID::Kminus),
	getParticleData(ParticleID::piplus),
	getParticleData(ParticleID::piplus)};
	if(_maxwgt.empty()) _maxwgt.push_back(1.);
	PhaseSpaceModePtr mode = new_ptr(PhaseSpaceMode(in,out,_maxwgt[0]));
	vector<PhaseSpaceChannel> channels;
	if(k8920) {
	channels.push_back((PhaseSpaceChannel(mode),0,k8920,0,3,1,1,1,2));
	channels.push_back((PhaseSpaceChannel(mode),0,k8920,0,2,1,1,1,3));
	}
	if(k1430) {
	channels.push_back((PhaseSpaceChannel(mode),0,k1430,0,3,1,1,1,2));
	channels.push_back((PhaseSpaceChannel(mode),0,k1430,0,2,1,1,1,3));
	}
	if(k1680) {
	channels.push_back((PhaseSpaceChannel(mode),0,k1680,0,3,1,1,1,2));
	channels.push_back((PhaseSpaceChannel(mode),0,k1680,0,2,1,1,1,3));
	}
	// add the mode
	vector<double> wtemp;
	if(channels.size()<=_weights.size()) {
	vector<double>::const_iterator wit=_weights.begin();
	wtemp=vector<double>(wit,wit+channels.size());
	}
	else {
	wtemp=vector<double>(channels.size(),1./double(channels.size()));
	}
	for(unsigned int ix=0;ix<channels.size();++ix) {
	channels[ix].weight(wtemp[ix]);
	mode->addChannel(channels[ix]);
	}
	addMode(mode);
	// D0 -> K-pi+pi0
	if(_maxwgt.size()<2) _maxwgt.push_back(1.);
	in = getParticleData(ParticleID::D0);
	out = {getParticleData(ParticleID::Kminus),
	getParticleData(ParticleID::piplus),
	getParticleData(ParticleID::pi0)};
	mode = new_ptr(PhaseSpaceMode(in,out,_maxwgt[1]));
	int iy = channels.size();
	channels.clear();
	if(k8920) {
	channels.push_back((PhaseSpaceChannel(mode),0,k8920,0,3,1,1,1,2));
	}
	if(k892m) {
	channels.push_back((PhaseSpaceChannel(mode),0,k892m,0,2,1,1,1,3));
	}
	if(rhop) {
	channels.push_back((PhaseSpaceChannel(mode),0,rhop,0,1,1,2,1,3));
	}
	// add the mode
	if(iy+channels.size()<=_weights.size()) {
	vector<double>::const_iterator wit=_weights.begin();
	wtemp=vector<double>(wit+iy,wit+iy+channels.size());
	}
	else {
	wtemp=vector<double>(channels.size(),1./double(channels.size()));
	}
	iy+=channels.size();
	for(unsigned int ix=0;ix<channels.size();++ix) {
	channels[ix].weight(wtemp[ix]);
	mode->addChannel(channels[ix]);
	}
	addMode(mode);
	// D0 -> Kbar0 pi+ pi-
	in = getParticleData(ParticleID::D0);
	out = {getParticleData(ParticleID::Kbar0),
	getParticleData(ParticleID::piplus),
	getParticleData(ParticleID::piminus)};
	if(_maxwgt.size()<3) _maxwgt.push_back(1.);
	mode = new_ptr(PhaseSpaceMode(in,out,_maxwgt[2]));
	channels.clear();
	if(rho0) {
	channels.push_back((PhaseSpaceChannel(mode),0,rho0,0,1,1,2,1,3));
	}
	if(k892m) {
	channels.push_back((PhaseSpaceChannel(mode),0,k892m,0,2,1,1,1,3));
	}
	// add the mode
	if(iy+channels.size()<=_weights.size()) {
	vector<double>::const_iterator wit=_weights.begin();
	wtemp=vector<double>(wit+iy,wit+iy+channels.size());
	}
	else {
	wtemp=vector<double>(channels.size(),1./double(channels.size()));
	}
	for(unsigned int ix=0;ix<channels.size();++ix) {
	channels[ix].weight(wtemp[ix]);
	mode->addChannel(channels[ix]);
	}
	addMode(mode);
	// reset the resonance parameters in the integration if needed
	if(_localparameters) {
	resetIntermediate(k8920,_mK8920,_wK8920);
	resetIntermediate(k892m,_mK892m,_wK892m);
	resetIntermediate(k1430,_mK1680,_wK1680);
	resetIntermediate(k1680,_mK1430,_wK1430);
	resetIntermediate(rho0 ,_mrho0 ,_wrho0 );
	resetIntermediate(rhop ,_mrhop ,_wrhop );
	}
	// get values from the ParticleData objects if needed
	else {
	_mK8920 = k8920->mass();
	_mK892m = k892m->mass();
	_mK1680 = k1430->mass();
	_mK1430 = k1680->mass();
	_mrho0 = rho0 ->mass();
	_mrhop = rhop ->mass();
	_wK8920 = k8920->width();
	_wK892m = k892m->width();
	_wK1680 = k1430->width();
	_wK1430 = k1680->width();
	_wrho0 = rho0 ->width();
	_wrhop = rhop ->width();
	}
	}

	void DtoKPiPiE691::persistentOutput(PersistentOStream & os) const {
	os << _a1NR << _phi1NR << _a1K892 << _phi1K892 << _a1K1430 << _phi1K1430
	<< _a1K1680 << _phi1K1680 << _a2NR << _phi2NR << _a2K8920 << _phi2K8920
	<< _a2K892m << _phi2K892m << _a2rho << _phi2rho << _a3NR << _phi3NR
	<< _a3K892 << _phi3K892 << _a3rho << _phi3rho << _c1NR << _c1K892 << _c1K1430
	<< _c1K1680 << _c2NR << _c2K8920 << _c2K892m << _c2rho << _c3NR << _c3K892
	<< _c3rho << _localparameters << ounit(_mK8920,GeV) << ounit(_wK8920,GeV)
	<< ounit(_mK892m,GeV) << ounit(_wK892m,GeV) << ounit(_mK1680,GeV)
	<< ounit(_wK1680,GeV) << ounit(_mK1430,GeV) << ounit(_wK1430,GeV)
	<< ounit(_mrho0 ,GeV) << ounit(_wrho0 ,GeV) << ounit(_mrhop ,GeV)
	<< ounit(_wrhop ,GeV) << _maxwgt << _weights;
	}

	void DtoKPiPiE691::persistentInput(PersistentIStream & is, int) {
	is >> _a1NR >> _phi1NR >> _a1K892 >> _phi1K892 >> _a1K1430 >> _phi1K1430
	>> _a1K1680 >> _phi1K1680 >> _a2NR >> _phi2NR >> _a2K8920 >> _phi2K8920
	>> _a2K892m >> _phi2K892m >> _a2rho >> _phi2rho >> _a3NR >> _phi3NR
	>> _a3K892 >> _phi3K892 >> _a3rho >> _phi3rho >> _c1NR >> _c1K892 >> _c1K1430
	>> _c1K1680 >> _c2NR >> _c2K8920 >> _c2K892m >> _c2rho >> _c3NR >> _c3K892
	>> _c3rho >> _localparameters >> iunit(_mK8920,GeV) >> iunit(_wK8920,GeV)
	>> iunit(_mK892m,GeV) >> iunit(_wK892m,GeV) >> iunit(_mK1680,GeV)
	>> iunit(_wK1680,GeV) >> iunit(_mK1430,GeV) >> iunit(_wK1430,GeV)
	>> iunit(_mrho0 ,GeV) >> iunit(_wrho0 ,GeV) >> iunit(_mrhop ,GeV)
	>> iunit(_wrhop ,GeV) >> _maxwgt >> _weights;
	}

	// The following static variable is needed for the type
	// description system in ThePEG.
	DescribeClass<DtoKPiPiE691,DecayIntegrator>
	describeHerwigDtoKPiPiE691("Herwig::DtoKPiPiE691", "HwSMDecay.so");

	void DtoKPiPiE691::Init() {

	static ClassDocumentation<DtoKPiPiE691> documentation
	("The DtoKPiPiE691 class implements the model of E691 for D -> K pi pi"
	"Dalitz decays",
	"The model of \\cite{Anjos:1992kb} for $D\\to K\\pi\\pi$ decays was used",
	"\\bibitem{Anjos:1992kb} J.~C.~Anjos {\\it et al.} [E691 Collaboration],"
	"Phys.\\ Rev.\\ D {\\bf 48} (1993) 56.");

	static Parameter<DtoKPiPiE691,double> interfaceKmPipPipNonResonantMagnitude
	("KmPipPipNonResonantMagnitude",
	"The magnitude of the non-resonant component for D+ -> K- pi+ pi+",
	&DtoKPiPiE691::_a1NR, 1.0, 0.0, 10.0,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiE691,double> interfaceKmPipPipNonResonantPhase
	("KmPipPipNonResonantPhase",
	"The phase of the non-resonant component for D+ -> K- pi+ pi+",
	&DtoKPiPiE691::_phi1NR, 0., -180., 180.,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiE691,double> interfaceKmPipPipK892Magnitude
	("KmPipPipK892Magnitude",
	"The magnitude of the K*(892) component for D+ -> K- pi+ pi+",
	&DtoKPiPiE691::_a1K892, 0.78, 0.0, 10.0,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiE691,double> interfaceKmPipPipK892Phase
	("KmPipPipK892Phase",
	"The phase of the K*(892) component for D+ -> K- pi+ pi+",
	&DtoKPiPiE691::_phi1K892, -60., -180., 180.,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiE691,double> interfaceKmPipPipK1430Magnitude
	("KmPipPipK1430Magnitude",
	"The magnitude of the K*0(1430) component for D+ -> K- pi+ pi+",
	&DtoKPiPiE691::_a1K1430, 0.53, 0.0, 10.0,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiE691,double> interfaceKmPipPipK1430Phase
	("KmPipPipK1430Phase",
	"The phase of the K*0(1430) component for D+ -> K- pi+ pi+",
	&DtoKPiPiE691::_phi1K1430, 132., -180., 180.,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiE691,double> interfaceKmPipPipK1680Magnitude
	("KmPipPipK1680Magnitude",
	"The magnitude of the K*(1680) component for D+ -> K- pi+ pi+",
	&DtoKPiPiE691::_a1K1680, 0.47, 0.0, 10.0,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiE691,double> interfaceKmPipPipK1680Phase
	("KmPipPipK1680Phase",
	"The phase of the K*(1680) component for D+ -> K- pi+ pi+",
	&DtoKPiPiE691::_phi1K1680, - 51., -180., 180.,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiE691,double> interfaceKmPipPi0NonResonantMagnitude
	("KmPipPi0NonResonantMagnitude",
	"The magnitude of the non-resonant component for D0 -> K- pi+ pi0",
	&DtoKPiPiE691::_a2NR, 1.0, 0.0, 10.0,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiE691,double> interfaceKmPipPi0NonResonantPhase
	("KmPipPi0NonResonantPhase",
	"The phase of the non-resonant component for D0 -> K- pi+ pi0",
	&DtoKPiPiE691::_phi2NR, 0., -180., 180.,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiE691,double> interfaceKmPipPi0K8920Magnitude
	("KmPipPi0K8920Magnitude",
	"The magnitude of the K*(892)0 component for D0 -> K- pi+ pi0",
	&DtoKPiPiE691::_a2K8920, 3.19, 0.0, 10.0,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiE691,double> interfaceKmPipPi0K8920Phase
	("KmPipPi0K8920Phase",
	"The phase of the K*(892)0 component for D0 -> K- pi+ pi0",
	&DtoKPiPiE691::_phi2K8920, 167., -180., 180.,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiE691,double> interfaceKmPipPi0K892mMagnitude
	("KmPipPi0K892mMagnitude",
	"The magnitude of the K*(892)- component for D0 -> K- pi+ pi0",
	&DtoKPiPiE691::_a2K892m, 2.96, 0.0, 10.0,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiE691,double> interfaceKmPipPi0K892mPhase
	("KmPipPi0K892mPhase",
	"The phase of the K*(892)- component for D0 -> K- pi+ pi0",
	&DtoKPiPiE691::_phi2K892m, -112., -180., 180.,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiE691,double> interfaceKmPipPi0RhoMagnitude
	("KmPipPi0RhoMagnitude",
	"The magnitude of the rho component for D0 -> K- pi+ pi0",
	&DtoKPiPiE691::_a2rho, 8.56, 0.0, 10.0,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiE691,double> interfaceKmPipPi0RhoPhase
	("KmPipPi0RhoPhase",
	"The phase of the rho component for D0 -> K- pi+ pi0",
	&DtoKPiPiE691::_phi2rho, 40., -180., 180.,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiE691,double> interfaceK0PipPimNonResonantMagnitude
	("K0PipPimNonResonantMagnitude",
	"The magnitude of the non-resonant component for D0 -> Kbar0 pi+pi-",
	&DtoKPiPiE691::_a3NR, 1.0, 0.0, 10.0,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiE691,double> interfaceK0PipPimNonResonantPhase
	("K0PipPimNonResonantPhase",
	"The phase of the non-resonant component for D0 -> Kbar0 pi+pi-",
	&DtoKPiPiE691::_phi3NR, 0., -180., 180.,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiE691,double> interfaceK0PipPimK892Magnitude
	("K0PipPimK892Magnitude",
	"The magnitude of the K*(892) component for D0 -> Kbar0 pi+pi-",
	&DtoKPiPiE691::_a3K892, 2.31, 0.0, 10.0,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiE691,double> interfaceK0PipPimK892Phase
	("K0PipPimK892Phase",
	"The phase of the K*(892)0 component for D0 -> Kbar0 pi+pi-",
	&DtoKPiPiE691::_phi3K892, 109., -180., 180.,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiE691,double> interfaceK0PipPimRhoMagnitude
	("K0PipPimRhoMagnitude",
	"The magnitude of the rho component for D0 -> Kbar0 pi+pi-",
	&DtoKPiPiE691::_a3rho, 1.59, 0.0, 10.0,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiE691,double> interfaceK0PipPimRhoPhase
	("K0PipPimRhoPhase",
	"The phase of the rho component for D0 -> Kbar0 pi+pi-",
	&DtoKPiPiE691::_phi3rho, -123., -180., 180.,
	false, false, Interface::limited);

	static Switch<DtoKPiPiE691,bool> interfaceLocalParameters
	("LocalParameters",
	"Whether to use local values for the masses and widths or"
	" those from the ParticleData objects",
	&DtoKPiPiE691::_localparameters, true, false, false);
	static SwitchOption interfaceLocalParametersLocal
	(interfaceLocalParameters,
	"Local",
	"Use local values",
	true);
	static SwitchOption interfaceLocalParametersParticleData
	(interfaceLocalParameters,
	"ParticleData",
	"Use the values from the ParticleData objects",
	false);

	static Parameter<DtoKPiPiE691,Energy> interfaceK8920Mass
	("K8920Mass",
	"The mass of the K*(892)0",
	&DtoKPiPiE691::_mK8920, GeV, 0.8961 GeV, ZERO, 10.0GeV,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiE691,Energy> interfaceK8920Width
	("K8920Width",
	"The width of the K*(892)0",
	&DtoKPiPiE691::_wK8920, GeV, 0.0505GeV, ZERO, 10.0GeV,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiE691,Energy> interfaceK892MinusMass
	("K892MinusMass",
	"The mass of the K*(892)-",
	&DtoKPiPiE691::_mK892m, GeV, 0.89159GeV, ZERO, 10.0GeV,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiE691,Energy> interfaceK892MinusWidth
	("K892MinusWidth",
	"The width of the K*(892)-",
	&DtoKPiPiE691::_wK892m, GeV, 0.0498GeV, ZERO, 10.0GeV,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiE691,Energy> interfaceK1680Mass
	("K1680Mass",
	"The mass of the K*(1680)",
	&DtoKPiPiE691::_mK1680, GeV, 1.714 GeV, ZERO, 10.0GeV,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiE691,Energy> interfaceK1680Width
	("K1680Width",
	"The width of the K*(1680)",
	&DtoKPiPiE691::_wK1680, GeV, 0.323 GeV, ZERO, 10.0GeV,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiE691,Energy> interfaceK1430Mass
	("K1430Mass",
	"The mass of the K*(1430)",
	&DtoKPiPiE691::_mK1430, GeV, 1.429 GeV, ZERO, 10.0GeV,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiE691,Energy> interfaceK1430Width
	("K1430Width",
	"The width of the K*(1430)",
	&DtoKPiPiE691::_wK1430, GeV, 0.287 GeV, ZERO, 10.0GeV,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiE691,Energy> interfaceRho0Mass
	("Rho0Mass",
	"The mass of the rho0",
	&DtoKPiPiE691::_mrho0, GeV, 0.7681 GeV, ZERO, 10.0GeV,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiE691,Energy> interfaceRho0Width
	("Rho0Width",
	"The width of the rho0",
	&DtoKPiPiE691::_wrho0, GeV, 0.1515GeV, ZERO, 10.0GeV,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiE691,Energy> interfaceRhoPlusMass
	("RhoPlusMass",
	"The mass of the rho+",
	&DtoKPiPiE691::_mrhop, GeV, 0.7681 GeV, ZERO, 10.0GeV,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiE691,Energy> interfaceRhoPlusWidth
	("RhoPlusWidth",
	"The width of the rho+",
	&DtoKPiPiE691::_wrhop, GeV, 0.1515GeV, ZERO, 10.0GeV,
	false, false, Interface::limited);

	static ParVector<DtoKPiPiE691,double> interfaceMaximumWeights
	("MaximumWeights",
	"The maximum weights for the unweighting of the decays",
	&DtoKPiPiE691::_maxwgt, -1, 1.0, 0.0, 1.0e11,
	false, false, Interface::limited);

	static ParVector<DtoKPiPiE691,double> interfaceWeights
	("Weights",
	"The weights for the different channels for the phase-space integration",
	&DtoKPiPiE691::_weights, -1, 1.0, 0.0, 1.0,
	false, false, Interface::limited);
	}

	int DtoKPiPiE691::modeNumber(bool & cc,tcPDPtr parent,
	const tPDVector & children) const {
	int id0(parent->id());
	// incoming particle must be D0 or D+
	if(abs(id0)!=ParticleID::D0&&abs(id0)!=ParticleID::Dplus) return -1;
	cc = id0<0;
	// must be three decay products
	if(children.size()!=3) return -1;
	tPDVector::const_iterator pit = children.begin();
	unsigned int npip(0),npim(0),nkm(0),nk0(0),npi0(0);
	int id;
	for( ;pit!=children.end();++pit) {
	id=(**pit).id();
	if(id ==ParticleID::piplus) ++npip;
	else if(id ==ParticleID::pi0) ++npi0;
	else if(id ==ParticleID::piminus) ++npim;
	else if(abs(id)==ParticleID::K0) ++nk0;
	else if(id ==ParticleID::K_L0) ++nk0;
	else if(id ==ParticleID::K_S0) ++nk0;
	else if(abs(id)==ParticleID::Kplus) ++nkm;
	}
	if(abs(id0)==ParticleID::Dplus) {
	if((id0==ParticleID::Dplus &&(nkm==1&&npip==2))\|\|
	(id0==ParticleID::Dminus&&(nkm==1&&npim==2))) return 0;
	else return -1;
	}
	else {
	if(npim==1&&npip==1&&nk0==1) return 2;
	else if(nkm==1&&(npip+npim)==1&&npi0==1) return 1;
	else return -1;
	}
	}

	void DtoKPiPiE691::
	constructSpinInfo(const Particle & part, ParticleVector decay) const {
	// set up the spin information for the decay products
	ScalarWaveFunction::constructSpinInfo(const_ptr_cast<tPPtr>(&part),
	incoming,true);
	for(unsigned int ix=0;ix<3;++ix)
	ScalarWaveFunction::constructSpinInfo(decay[ix],outgoing,true);
	}

	double DtoKPiPiE691::me2(const int ichan, const Particle & part,
	- const tPDVector & ,
	- const vector<Lorentz5Momentum> & momenta,
	- MEOption meopt) const {
	+ const tPDVector & ,
	+ const vector<Lorentz5Momentum> & momenta,
	+ MEOption meopt) const {
	if(!ME())
	ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin0,PDT::Spin0,PDT::Spin0,PDT::Spin0)));
	useMe();
	if(meopt==Initialize) {
	ScalarWaveFunction::
	calculateWaveFunctions(_rho,const_ptr_cast<tPPtr>(&part),incoming);
	}
	Complex amp;
	// D+ -> K-pi+pi+
	if(imode()==0) {
	Lorentz5Momentum pres1=momenta[0]+momenta[1];
	pres1.rescaleMass();
	double ct1 =-decayAngle(part.momentum(),pres1,momenta[0]);
	Lorentz5Momentum pres2=momenta[0]+momenta[2];
	pres2.rescaleMass();
	double ct2 =-decayAngle(part.momentum(),pres2,momenta[0]);
	if(ichan<0) {
	amp = _c1NR*sqrt(2.)
	+_c1K892 *amplitude(1,ct1,pres1.mass(),_wK8920,_mK8920)
	+_c1K892 *amplitude(1,ct2,pres2.mass(),_wK8920,_mK8920)
	+_c1K1430*amplitude(0,ct1,pres1.mass(),_wK1430,_mK1430)
	+_c1K1430*amplitude(0,ct2,pres2.mass(),_wK1430,_mK1430)
	+_c1K1680*amplitude(1,ct1,pres1.mass(),_wK1680,_mK1680)
	+_c1K1680*amplitude(1,ct2,pres2.mass(),_wK1680,_mK1680);
	}
	else if(ichan==0) {
	amp=_c1K892 *amplitude(1,ct1,pres1.mass(),_wK8920,_mK8920);
	}
	else if(ichan==1) {
	amp=_c1K892 *amplitude(1,ct2,pres2.mass(),_wK8920,_mK8920);
	}
	else if(ichan==2) {
	amp=_c1K1430*amplitude(1,ct1,pres1.mass(),_wK1430,_mK1430);
	}
	else if(ichan==3) {
	amp=_c1K1430*amplitude(1,ct2,pres2.mass(),_wK1430,_mK1430);
	}
	else if(ichan==4) {
	amp=_c1K1680*amplitude(1,ct1,pres1.mass(),_wK1680,_mK1680);
	}
	else if(ichan==5) {
	amp=_c1K1680*amplitude(1,ct2,pres2.mass(),_wK1680,_mK1680);
	}
	}
	// D0 -> K-pi+pi0
	else if(imode()==1) {
	Lorentz5Momentum pres1=momenta[0]+momenta[1];
	pres1.rescaleMass();
	double ct1 = decayAngle(part.momentum(),pres1,momenta[0]);
	Lorentz5Momentum pres2=momenta[0]+momenta[2];
	pres2.rescaleMass();
	double ct2 = decayAngle(part.momentum(),pres2,momenta[2]);
	Lorentz5Momentum pres3=momenta[1]+momenta[2];
	pres3.rescaleMass();
	double ct3 = decayAngle(part.momentum(),pres3,momenta[2]);
	if(ichan<0) {
	amp = _c2NR
	+_c2K8920*amplitude(1,ct1,pres1.mass(),_wK8920,_mK8920)
	+_c2K892m*amplitude(1,ct2,pres2.mass(),_wK892m,_mK892m)
	+_c2rho *amplitude(1,ct3,pres3.mass(),_wrhop ,_mrhop );
	}
	else if(ichan==0) {
	amp = _c2K8920*amplitude(1,ct1,pres1.mass(),_wK8920,_mK8920);
	}
	else if(ichan==1) {
	amp = _c2K892m*amplitude(1,ct2,pres2.mass(),_wK892m,_mK892m);
	}
	else if(ichan==2) {
	amp = _c2rho *amplitude(1,ct3,pres3.mass(),_wrhop ,_mrhop);
	}
	}
	// D0 -> Kbar0pi+pi-
	else if(imode()==2) {
	Lorentz5Momentum pres1=momenta[0]+momenta[2];
	pres1.rescaleMass();
	double ct1 = decayAngle(part.momentum(),pres1,momenta[0]);
	Lorentz5Momentum pres2=momenta[1]+momenta[2];
	pres2.rescaleMass();
	double ct2 = decayAngle(part.momentum(),pres2,momenta[1]);
	if(ichan<0) {
	amp = _c3NR
	+_c3K892*amplitude(1,ct1,pres1.mass(),_wK892m,_mK892m)
	+_c3rho *amplitude(1,ct2,pres2.mass(),_wrho0 ,_mrho0 );
	}
	else if(ichan==0) {
	amp = _c3rho *amplitude(1,ct2,pres2.mass(),_wrho0 ,_mrho0 );
	}
	else if(ichan==1) {
	amp = _c3K892*amplitude(1,ct1,pres1.mass(),_wK892m,_mK892m);
	}
	}
	// now compute the matrix element
	(*ME())(0,0,0,0)=amp;
	return norm(amp);
	}

	void DtoKPiPiE691::dataBaseOutput(ofstream & output, bool header) const {
	if(header) output << "update decayers set parameters=\"";
	// parameters for the DecayIntegrator base class
	DecayIntegrator::dataBaseOutput(output,false);
	// parameters
	output << "newdef " << name() << ":KmPipPipNonResonantMagnitude "
	<< _a1NR << "\n";
	output << "newdef " << name() << ":KmPipPipNonResonantPhase "
	<< _phi1NR << "\n";
	output << "newdef " << name() << ":KmPipPipK892Magnitude "
	<< _a1K892 << "\n";
	output << "newdef " << name() << ":KmPipPipK892Phase "
	<< _phi1K892 << "\n";
	output << "newdef " << name() << ":KmPipPipK1430Magnitude "
	<< _a1K1430 << "\n";
	output << "newdef " << name() << ":KmPipPipK1430Phase "
	<< _phi1K1430 << "\n";
	output << "newdef " << name() << ":KmPipPipK1680Magnitude "
	<< _a1K1680 << "\n";
	output << "newdef " << name() << ":KmPipPipK1680Phase "
	<< _phi1K1680 << "\n";
	output << "newdef " << name() << ":KmPipPi0NonResonantMagnitude "
	<< _a2NR << "\n";
	output << "newdef " << name() << ":KmPipPi0NonResonantPhase "
	<< _phi2NR << "\n";
	output << "newdef " << name() << ":KmPipPi0K8920Magnitude "
	<< _a2K8920 << "\n";
	output << "newdef " << name() << ":KmPipPi0K8920Phase "
	<< _phi2K8920 << "\n";
	output << "newdef " << name() << ":KmPipPi0K892mMagnitude "
	<< _a2K892m << "\n";
	output << "newdef " << name() << ":KmPipPi0K892mPhase "
	<< _phi2K892m << "\n";
	output << "newdef " << name() << ":KmPipPi0RhoMagnitude "
	<< _a2rho << "\n";
	output << "newdef " << name() << ":KmPipPi0RhoPhase "
	<< _phi2rho << "\n";
	output << "newdef " << name() << ":K0PipPimNonResonantMagnitude "
	<< _a3NR << "\n";
	output << "newdef " << name() << ":K0PipPimNonResonantPhase "
	<< _phi3NR << "\n";
	output << "newdef " << name() << ":K0PipPimK892Magnitude "
	<< _a3K892 << "\n";
	output << "newdef " << name() << ":K0PipPimK892Phase "
	<< _phi3K892 << "\n";
	output << "newdef " << name() << ":K0PipPimRhoMagnitude "
	<< _a3rho << "\n";
	output << "newdef " << name() << ":K0PipPimRhoPhase "
	<< _phi3rho << "\n";
	output << "newdef " << name() << ":LocalParameters " << _localparameters << "\n";
	output << "newdef " << name() << ":K8920Mass " << _mK8920/GeV << "\n";
	output << "newdef " << name() << ":K8920Width " << _wK8920/GeV << "\n";
	output << "newdef " << name() << ":K892MinusMass " << _mK892m/GeV << "\n";
	output << "newdef " << name() << ":K892MinusWidth " << _wK892m/GeV << "\n";
	output << "newdef " << name() << ":K1680Mass " << _mK1680/GeV << "\n";
	output << "newdef " << name() << ":K1680Width " << _wK1680/GeV << "\n";
	output << "newdef " << name() << ":K1430Mass " << _mK1430/GeV << "\n";
	output << "newdef " << name() << ":K1430Width " << _wK1430/GeV << "\n";
	output << "newdef " << name() << ":Rho0Mass " << _mrho0 /GeV << "\n";
	output << "newdef " << name() << ":Rho0Width " << _wrho0 /GeV << "\n";
	output << "newdef " << name() << ":RhoPlusMass " << _mrhop /GeV << "\n";
	output << "newdef " << name() << ":RhoPlusWidth " << _wrhop /GeV << "\n";
	for(unsigned int ix=0;ix<_maxwgt.size();++ix) {
	output << "insert " << name() << ":MaximumWeights "
	<< ix << " " << _maxwgt[ix] << "\n";
	}
	for(unsigned int ix=0;ix<_weights.size();++ix) {
	output << "insert " << name() << ":Weights "
	<< ix << " " << _weights[ix] << "\n";
	}
	if(header) {
	output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl;
	}
	}

	void DtoKPiPiE691::doinitrun() {
	DecayIntegrator::doinitrun();
	_weights.resize(mode(0)->channels().size()+mode(1)->channels().size()+
	mode(2)->channels().size());
	_maxwgt.resize(3);
	unsigned int iy=0;
	for(unsigned int ix=0;ix<3;++ix) {
	_maxwgt[ix]=mode(ix)->maxWeight();
	for(unsigned int iz=0;iz<mode(ix)->channels().size();++iz) {
	_weights[iy]=mode(ix)->channels()[iz].weight();
	++iy;
	}
	}
	}
	diff --git a/Decay/ScalarMeson/DtoKPiPiE691.h b/Decay/ScalarMeson/DtoKPiPiE691.h
	--- a/Decay/ScalarMeson/DtoKPiPiE691.h
	+++ b/Decay/ScalarMeson/DtoKPiPiE691.h
	@@ -1,458 +1,459 @@
	// -- C++ --
	//
	// DtoKPiPiE691.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2019 The Herwig Collaboration
	//
	// Herwig is licenced under version 3 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef HERWIG_DtoKPiPiE691_H
	#define HERWIG_DtoKPiPiE691_H
	//
	// This is the declaration of the DtoKPiPiE691 class.
	//

	#include "Herwig/Decay/DecayIntegrator.h"
	#include "Herwig/Decay/PhaseSpaceMode.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	- * Here is the documentation of the DtoKPiPiE691 class.
	+ * The DtoKPiPiE691 class implements the model of E691 (Phys. Rev. D48 (1993) 56)
	+ * for D -> K pi pi Dalitz decays.
	*
	* @see \ref DtoKPiPiE691Interfaces "The interfaces"
	* defined for DtoKPiPiE691.
	*/
	class DtoKPiPiE691: public DecayIntegrator {

	public:
	/**
	* The default constructor.
	*/
	DtoKPiPiE691();

	/**
	* Which of the possible decays is required
	* @param cc Is this mode the charge conjugate
	* @param parent The decaying particle
	* @param children The decay products
	*/
	virtual int modeNumber(bool & cc, tcPDPtr parent,
	const tPDVector & children) const;

	/**
	* Return the matrix element squared for a given mode and phase-space channel.
	* @param ichan The channel we are calculating the matrix element for.
	* @param part The decaying Particle.
	* @param outgoing The particles produced in the decay
	* @param momenta The momenta of the particles produced in the decay
	* @param meopt Option for the calculation of the matrix element
	* @return The matrix element squared for the phase-space configuration.
	*/
	double me2(const int ichan,const Particle & part,
	const tPDVector & outgoing,
	const vector<Lorentz5Momentum> & momenta,
	MEOption meopt) const;

	/**
	* Construct the SpinInfos for the particles produced in the decay
	*/
	virtual void constructSpinInfo(const Particle & part,
	ParticleVector outgoing) const;

	/**
	* Output the setup information for the particle database
	* @param os The stream to output the information to
	* @param header Whether or not to output the information for MySQL
	*/
	virtual void dataBaseOutput(ofstream & os,bool header) 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:

	/**
	* Methods to calculate the amplitudes for a given channel
	*/
	//@{
	/**
	* Calculate the decay angle for the amplitude, the angle is the
	* angle between the 2 and 3 for the decay \f$D\to(12)3\f$ in the rest frame of
	* the resonance which decays to 1 and 2.
	* @param pparent The momentum of the parent
	* @param pres The momentum of the resonance
	* @param p1 The momentum of the first decay product of the resonance
	*/
	double decayAngle(const Lorentz5Momentum & pparent,
	const Lorentz5Momentum & pres,
	const Lorentz5Momentum & p1) const {
	Energy2 dot = pparentp1, mREp = prespparent;
	Energy2 mRE1 = pres*p1, mp2 = pparent.mass2();
	Energy2 mres2 = pres.mass2(), m12 = p1.mass2();
	return (dotmres2-mREpmRE1)/
	sqrt((mREpmREp-mres2mp2)(mRE1mRE1-mres2*m12));
	}

	/**
	* Calculate the amplitude
	* @param ispin The spin of the resonance
	* @param costheta The decay angle
	* @param mAB The off-shell mass of the resonance
	* @param wres The width of the resonance
	* @param mres The on-shell mass of the resonance
	*/
	Complex amplitude(int ispin, double costheta,Energy mAB,
	Energy wres, Energy mres) const {
	double s = 0.;
	switch(ispin) {
	case 0: s = 1.; break;
	case 1: s = costheta; break;
	case 2: s = 1.5*sqr(costheta)-0.5; break;
	default: assert(false);
	}
	Complex bw = sqrt(0.5wres/GeV/Constants::pi)GeV/
	(mAB-mres-complex<Energy>(ZERO,0.5*wres));
	return s*bw;
	}
	//@}

	protected:

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

	/** 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 {return new_ptr(*this);}
	//@}

	protected:

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

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

	private:

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

	private:

	/**
	* Amplitudes and phases for the different components
	*/
	//@{
	/**
	* Amplitude of the non-resonant component for \f$D^+\to K^-\pi^+\pi^+\f$
	*/
	double _a1NR;

	/**
	* Phase of the non-resonant component for \f$D^+\to K^-\pi^+\pi^+\f$
	*/
	double _phi1NR;

	/**
	* Amplitude of the \f$\bar{K}^*(892)^0\f$ component for \f$D^+\to K^-\pi^+\pi^+\f$
	*/
	double _a1K892;

	/**
	* Phase of the \f$\bar{K}^*(892)^0\f$ component for \f$D^+\to K^-\pi^+\pi^+\f$
	*/
	double _phi1K892;

	/**
	* Amplitude of the \f$\bar{K}^*_0(1430)^0\f$ component for \f$D^+\to K^-\pi^+\pi^+\f$
	*/
	double _a1K1430;

	/**
	* Phase of the \f$\bar{K}^*_0(1430)^0\f$ component for \f$D^+\to K^-\pi^+\pi^+\f$
	*/
	double _phi1K1430;

	/**
	* Amplitude of the \f$\bar{K}^*_0(1680)^0\f$ component for \f$D^+\to K^-\pi^+\pi^+\f$
	*/
	double _a1K1680;

	/**
	* Phase of the \f$\bar{K}^*_0(1680)^0\f$ component for \f$D^+\to K^-\pi^+\pi^+\f$
	*/
	double _phi1K1680;

	/**
	* Amplitude of the non-resonant component for \f$D^0\to K^-\pi^+\pi^0\f$
	*/
	double _a2NR;

	/**
	* Phase of the non-resonant component for \f$D^0\to K^-\pi^+\pi^0\f$
	*/
	double _phi2NR;

	/**
	* Amplitude of the \f$\bar{K}^*(892)^0\f$ component for \f$D^0\to K^-\pi^+\pi^0\f$
	*/
	double _a2K8920;

	/**
	* Phase of the \f$\bar{K}^*(892)^0\f$ component for \f$D^0\to K^-\pi^+\pi^0\f$
	*/
	double _phi2K8920;

	/**
	* Amplitude of the \f$K^*(892)^-\f$ component for \f$D^0\to K^-\pi^+\pi^0\f$
	*/
	double _a2K892m;

	/**
	* Phase of the \f$K^*(892)^-\f$ component for \f$D^0\to K^-\pi^+\pi^0\f$
	*/
	double _phi2K892m;

	/**
	* Amplitude of the \f$\rho^+\f$ component for \f$D^0\to K^-\pi^+\pi^0\f$
	*/
	double _a2rho;

	/**
	* Phase of the \f$\rho^+\f$ component for \f$D^0\to K^-\pi^+\pi^0\f$
	*/
	double _phi2rho;

	/**
	* Amplitude of the non-resonant component for \f$D^0\to \bar{K}^0\pi^+\pi^-\f$
	*/
	double _a3NR;

	/**
	* Phase of the non-resonant component for \f$D^0\to \bar{K}^0\pi^+\pi^-\f$
	*/
	double _phi3NR;

	/**
	* Amplitude of the \f$K^*(892)^-\f$ component for \f$D^0\to \bar{K}^0\pi^+\pi^-\f$
	*/
	double _a3K892;

	/**
	* Phase of the \f$K^*(892)^-\f$ component for \f$D^0\to \bar{K}^0\pi^+\pi^-\f$
	*/
	double _phi3K892;

	/**
	* Amplitude of the \f$\rho^0\f$ component for \f$D^0\to \bar{K}^0\pi^+\pi^-\f$
	*/
	double _a3rho;

	/**
	* Phase of the \f$\rho^0\f$ component for \f$D^0\to \bar{K}^0\pi^+\pi^-\f$
	*/
	double _phi3rho;
	//@}

	/**
	* Complex amplitudes for use in the matrix element
	*/
	//@{
	/**
	* Amplitude of the non-resonant component for \f$D^+\to K^-\pi^+\pi^+\f$
	*/
	Complex _c1NR;

	/**
	* Amplitude of the \f$\bar{K}^*(892)^0\f$ component for \f$D^+\to K^-\pi^+\pi^+\f$
	*/
	Complex _c1K892;

	/**
	* Amplitude of the \f$\bar{K}^*_0(1430)^0\f$ component for \f$D^+\to K^-\pi^+\pi^+\f$
	*/
	Complex _c1K1430;

	/**
	* Amplitude of the \f$\bar{K}^*_0(1680)^0\f$ component for \f$D^+\to K^-\pi^+\pi^+\f$
	*/
	Complex _c1K1680;

	/**
	* Amplitude of the non-resonant component for \f$D^0\to K^-\pi^+\pi^0\f$
	*/
	Complex _c2NR;

	/**
	* Amplitude of the \f$\bar{K}^*(892)^0\f$ for \f$D^0\to K^-\pi^+\pi^0\f$
	*/
	Complex _c2K8920;

	/**
	* Amplitude of the \f$K^*(892)^-\f$ for \f$D^0\to K^-\pi^+\pi^0\f$
	*/
	Complex _c2K892m;

	/**
	* Amplitude of the \f$\rho^+\f$ for \f$D^0\to K^-\pi^+\pi^0\f$
	*/
	Complex _c2rho;

	/**
	* Amplitude of the non-resonant component for \f$D^0\to \bar{K}^0\pi^+\pi^-\f$
	*/
	Complex _c3NR;

	/**
	* Amplitude of the \f$K^*(892)^-\f$ component for \f$D^0\to \bar{K}^0\pi^+\pi^-\f$
	*/
	Complex _c3K892;

	/**
	* Amplitude of the \f$\rho^0\f$ component for \f$D^0\to \bar{K}^0\pi^+\pi^-\f$
	*/
	Complex _c3rho;
	//@}

	/**
	* Masses and widths of the various resonances
	*/
	//@{
	/**
	* Use local values for the masses and widths
	*/
	bool _localparameters;

	/**
	* Mass of the \f$K^*(892)^0\f$
	*/
	Energy _mK8920;

	/**
	* Width of the \f$K^*(892)^0\f$
	*/
	Energy _wK8920;

	/**
	* Mass of the \f$K^*(892)^-\f$
	*/
	Energy _mK892m;

	/**
	* Width of the \f$K^*(892)^-\f$
	*/
	Energy _wK892m;

	/**
	* Mass of the \f$K^*(1680)^0\f$
	*/
	Energy _mK1680;

	/**
	* Width of the \f$K^*(1680)^0\f$
	*/
	Energy _wK1680;

	/**
	* Mass of the \f$K^*_0(1430)^0\f$
	*/
	Energy _mK1430;

	/**
	* Width of the \f$K^*_0(1430)^0\f$
	*/
	Energy _wK1430;

	/**
	* Mass of the \f$\rho^0\f$
	*/
	Energy _mrho0;

	/**
	* Width of the \f$\rho^0\f$
	*/
	Energy _wrho0;

	/**
	* Mass of the \f$\rho^+\f$
	*/
	Energy _mrhop;

	/**
	* Width of the \f$\rho^+\f$
	*/
	Energy _wrhop;
	//@}

	/**
	* Parameters for the phase-space integration
	*/
	//@{
	/**
	* Maximum weights for the various modes
	*/
	vector<double> _maxwgt;

	/**
	* Weights for the different integration channels
	*/
	vector<double> _weights;
	//@}

	/**
	* Spin density matrix
	*/
	mutable RhoDMatrix _rho;
	};

	}

	#endif /* HERWIG_DtoKPiPiE691_H */
	diff --git a/Decay/ScalarMeson/DtoKPiPiFOCUS.cc b/Decay/ScalarMeson/DtoKPiPiFOCUS.cc
	new file mode 100644
	--- /dev/null
	+++ b/Decay/ScalarMeson/DtoKPiPiFOCUS.cc
	@@ -0,0 +1,391 @@
	+// -- C++ --
	+//
	+// This is the implementation of the non-inlined, non-templated member
	+// functions of the DtoKPiPiFOCUS class.
	+//
	+
	+#include "DtoKPiPiFOCUS.h"
	+#include "ThePEG/Interface/ClassDocumentation.h"
	+#include "ThePEG/Interface/Reference.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"
	+
	+using namespace Herwig;
	+
	+DtoKPiPiFOCUS::DtoKPiPiFOCUS() : rD0_(5./GeV), rRes_(1.5/GeV),
	+ g1_(0.31072GeV), g2_(-0.02323GeV),
	+ beta_(3.389*GeV), theta_(286), gamma_({304,126,211}),
	+ c1_({1.655,0.780,-0.954}), c2_(17.182), c3_(0.734),
	+ maxWgt_(1.),weights_({0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1}),
	+ mRes_({0.896 GeV, 1.414GeV, 1.717GeV, 1.4324GeV}),
	+ wRes_({0.0503GeV, 0.232GeV, 0.322GeV, 0.109 GeV})
	+{}
	+
	+IBPtr DtoKPiPiFOCUS::clone() const {
	+ return new_ptr(*this);
	+}
	+
	+IBPtr DtoKPiPiFOCUS::fullclone() const {
	+ return new_ptr(*this);
	+}
	+
	+void DtoKPiPiFOCUS::persistentOutput(PersistentOStream & os) const {
	+ os << ounit(beta_,GeV) << theta_ << gamma_
	+ << ounit(g1_,GeV) << ounit(g2_,GeV) << c1_ << c2_ << c3_
	+ << KHalf_ << KThreeHalf_ << maxWgt_ << weights_
	+ << ounit(rD0_,1./GeV) << ounit(rRes_,1./GeV)
	+ << ounit(mRes_,GeV) << ounit(wRes_,GeV);
	+}
	+
	+void DtoKPiPiFOCUS::persistentInput(PersistentIStream & is, int) {
	+ is >> iunit(beta_,GeV) >> theta_ >> gamma_
	+ >> iunit(g1_,GeV) >> iunit(g2_,GeV) >> c1_ >> c2_ >> c3_
	+ >> KHalf_ >> KThreeHalf_ >> maxWgt_ >> weights_
	+ >> iunit(rD0_,1./GeV) >> iunit(rRes_,1./GeV)
	+ >> iunit(mRes_,GeV) >> iunit(wRes_,GeV);
	+}
	+
	+
	+// The following static variable is needed for the type
	+// description system in ThePEG.
	+DescribeClass<DtoKPiPiFOCUS,DecayIntegrator>
	+describeHerwigDtoKPiPiFOCUS("Herwig::DtoKPiPiFOCUS", "HwSMDecay.so");
	+
	+void DtoKPiPiFOCUS::Init() {
	+
	+ static ClassDocumentation<DtoKPiPiFOCUS> documentation
	+ ("There is no documentation for the DtoKPiPiFOCUS class");
	+
	+ static Reference<DtoKPiPiFOCUS,KMatrix> interfaceKHalf
	+ ("KHalf",
	+ "The K-matrix for the I=1/2 s-wave component",
	+ &DtoKPiPiFOCUS::KHalf_, false, false, true, false, false);
	+
	+ static Reference<DtoKPiPiFOCUS,KMatrix> interfaceKThreeHalf
	+ ("KThreeHalf",
	+ "The K-matrix for the I=1/2 s-wave component",
	+ &DtoKPiPiFOCUS::KThreeHalf_, false, false, true, false, false);
	+
	+ static Parameter<DtoKPiPiFOCUS,Energy> interfaceg1
	+ ("g1",
	+ "The coupling of the first channel to K*0",
	+ &DtoKPiPiFOCUS::g1_, 0.31072GeV, -10.GeV, 10.*GeV,
	+ false, false, Interface::limited);
	+
	+ static Parameter<DtoKPiPiFOCUS,Energy> interfaceg2
	+ ("g2",
	+ "The coupling of the first channel to K*0",
	+ &DtoKPiPiFOCUS::g2_, -0.02323GeV, -10.GeV, 10.*GeV,
	+ false, false, Interface::limited);
	+
	+ static Parameter<DtoKPiPiFOCUS,Energy> interfacebeta
	+ ("beta",
	+ "The beta coefficient for the P-vector",
	+ &DtoKPiPiFOCUS::beta_, 3.389GeV, -10.GeV, 10.*GeV,
	+ false, false, Interface::limited);
	+
	+ static Parameter<DtoKPiPiFOCUS,double> interfacetheta
	+ ("theta",
	+ "The theta phase (in degrees) for the P-vector",
	+ &DtoKPiPiFOCUS::theta_, 286, 0.0, 360.,
	+ false, false, Interface::limited);
	+
	+ static ParVector<DtoKPiPiFOCUS,double> interfacegamma
	+ ("gamma",
	+ "The gamma phases in degress",
	+ &DtoKPiPiFOCUS::gamma_, 3, 0., 0.0, 360.0,
	+ false, false, Interface::limited);
	+
	+ static ParVector<DtoKPiPiFOCUS,double> interfacec1
	+ ("c1",
	+ "The c1 coefficients for the P-vector",
	+ &DtoKPiPiFOCUS::c1_, -1, 0., -100., 100.,
	+ false, false, Interface::limited);
	+
	+ static ParVector<DtoKPiPiFOCUS,double> interfacec2
	+ ("c2",
	+ "The c2 coefficients for the P-vector",
	+ &DtoKPiPiFOCUS::c2_, -1, 0., -100., 100.,
	+ false, false, Interface::limited);
	+
	+ static ParVector<DtoKPiPiFOCUS,double> interfacec3
	+ ("c3",
	+ "The c3 coefficients for the P-vector",
	+ &DtoKPiPiFOCUS::c3_, -1, 0., -100., 100.,
	+ false, false, Interface::limited);
	+
	+ static Parameter<DtoKPiPiFOCUS,InvEnergy> interfaceDRadius
	+ ("DRadius",
	+ "The radius parameter for the Blatt-Weisskopf form-factor for the D",
	+ &DtoKPiPiFOCUS::rD0_, 1./GeV, 5./GeV, ZERO, 10./GeV,
	+ false, false, Interface::limited);
	+
	+ static Parameter<DtoKPiPiFOCUS,InvEnergy> interfaceResonanceRadius
	+ ("ResonanceRadius",
	+ "The radius parameter for the Blatt-Weisskopf form-factor for the"
	+ "intermediate resonances",
	+ &DtoKPiPiFOCUS::rRes_, 1./GeV, 1.5/GeV, ZERO, 10./GeV,
	+ false, false, Interface::limited);
	+
	+ static ParVector<DtoKPiPiFOCUS,Energy> interfaceMasses
	+ ("Masses",
	+ "The masses of the resonances",
	+ &DtoKPiPiFOCUS::mRes_, GeV, -1, 1.0GeV, 0.0GeV, 10.0*GeV,
	+ false, false, Interface::limited);
	+
	+ static ParVector<DtoKPiPiFOCUS,Energy> interfaceWidths
	+ ("Widths",
	+ "The widths of the resonances",
	+ &DtoKPiPiFOCUS::wRes_, GeV, -1, 1.0GeV, 0.0GeV, 10.0*GeV,
	+ false, false, Interface::limited);
	+
	+}
	+
	+int DtoKPiPiFOCUS::modeNumber(bool & cc,tcPDPtr parent,
	+ const tPDVector & children) const {
	+ int id0(parent->id());
	+ // incoming particle must be D+
	+ if(abs(id0)!=ParticleID::Dplus) return -1;
	+ cc = id0<0;
	+ // must be three decay products
	+ if(children.size()!=3) return -1;
	+ unsigned int npip(0),npim(0),nkm(0),nk0(0),npi0(0);
	+ int id;
	+ for(tPDPtr child : children) {
	+ id=child->id();
	+ if(id ==ParticleID::piplus) ++npip;
	+ else if(id ==ParticleID::pi0) ++npi0;
	+ else if(id ==ParticleID::piminus) ++npim;
	+ else if(abs(id)==ParticleID::K0) ++nk0;
	+ else if(id ==ParticleID::K_L0) ++nk0;
	+ else if(id ==ParticleID::K_S0) ++nk0;
	+ else if(abs(id)==ParticleID::Kplus) ++nkm;
	+ }
	+ if((id0==ParticleID::Dplus &&(nkm==1&&npip==2))\|\|
	+ (id0==ParticleID::Dminus&&(nkm==1&&npim==2))) return 0;
	+ else return -1;
	+}
	+
	+void DtoKPiPiFOCUS::
	+constructSpinInfo(const Particle & part, ParticleVector decay) const {
	+ // set up the spin information for the decay products
	+ ScalarWaveFunction::constructSpinInfo(const_ptr_cast<tPPtr>(&part),
	+ incoming,true);
	+ for(unsigned int ix=0;ix<3;++ix)
	+ ScalarWaveFunction::constructSpinInfo(decay[ix],outgoing,true);
	+}
	+
	+double DtoKPiPiFOCUS::me2(const int ichan, const Particle & part,
	+ const tPDVector & ,
	+ const vector<Lorentz5Momentum> & momenta,
	+ MEOption meopt) const {
	+ if(!ME())
	+ ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin0,PDT::Spin0,PDT::Spin0,PDT::Spin0)));
	+ useMe();
	+ if(meopt==Initialize) {
	+ ScalarWaveFunction::
	+ calculateWaveFunctions(rho_,const_ptr_cast<tPPtr>(&part),incoming);
	+ }
	+ Complex amp(0.);
	+ for(unsigned int ipi=1;ipi<3;++ipi) {
	+ unsigned int iother = ipi==1 ? 2 : 1;
	+ Lorentz5Momentum pres=momenta[0]+momenta[ipi];
	+ pres.rescaleMass();
	+ // momentum for the D Blatt-Weisskopf factor
	+ Energy pD = Kinematics::pstarTwoBodyDecay(part.mass(),pres.mass(),momenta[iother].mass());
	+ // momentum for the resonance Blatt-Weisskopf factor
	+ Energy pAB = Kinematics::pstarTwoBodyDecay(pres.mass(),momenta[0].mass(),momenta[ipi].mass());
	+ // loop over the resonances
	+ double rD = rD0_pD, rR = rRes_pAB;
	+ // most resonaces vectors so calculate that now
	+ double Fd = 1./sqrt(1.+sqr(rD));
	+ double Fr = 1./sqrt(1.+sqr(rR));
	+ unsigned int ipow=3;
	+ for(unsigned int ires=0;ires<4;++ires) {
	+ if(ichan>=0 && ichan+1!=int(ires*2+ipi)) continue;
	+ Energy pR = Kinematics::pstarTwoBodyDecay(mRes_[ires],momenta[0].mass(),momenta[ipi].mass());
	+ // last resonance is the K_2
	+ if(ires==3) {
	+ Fd = 1./sqrt(9.+3.*sqr(rD)+Math::powi(rD,4));
	+ Fr = 1./sqrt(9.+3.*sqr(rR)+Math::powi(rR,4));
	+ ipow = 5;
	+ }
	+ double rB = rRes_*pR;
	+ double Fb = ires<3 ? 1./sqrt(1.+sqr(rB)) : 1./sqrt(9.+3.*sqr(rB)+Math::powi(rB,4));
	+ Energy gamma = wRes_[ires]mRes_[ires]/pres.mass()sqr(Fr/Fb)*Math::powi(pAB/pR,ipow);
	+ complex<InvEnergy2> BW = FdFr/(sqr(mRes_[ires])-pres.mass2()-Complex(0.,1.)gamma*mRes_[ires]);
	+ // // and the decay momenta
	+ //
	+ // Energy pR = Kinematics::pstarTwoBodyDecay(mres,mA,mB);
	+ // double Fd(1.),Fr(1.),s(1.);
	+ // switch(ispin) {
	+ // case 0:
	+ // // default values of parameters are correct
	+ // break;
	+ // case 1:
	+ // Fr = sqrt((1.+sqr(_rrespR ))/(1.+sqr(_rrespAB )));
	+ // Fd =
	+ // s = ((mAC-mBC)(mAC+mBC)+(mD-mC)(mD+mC)(mB-mA)(mB+mA)/sqr(mres))/GeV2;
	+ // break;
	+ // case 2:
	+ // Fr = sqrt((9.+3.sqr(_rrespR )+Math::powi(_rres*pR ,4))/
	+ // (9.+3.sqr(_rrespAB )+Math::powi(_rres*pAB ,4)));
	+ // Fd =
	+ // s = sqr(((mBC-mAC)(mBC+mAC)+(mD-mC)(mD+mC)(mA-mB)(mA+mB)/sqr(mres))/GeV2)
	+ // -(mABmAB-2.mDmD-2.mCmC+sqr((mD-mC)(mD+mC))/sqr(mres))*
	+ // (mABmAB-2.mAmA-2.mBmB+sqr((mA-mB)(mA+mB))/sqr(mres))/3./GeV2/GeV2;
	+ // break;
	+ // default:
	+ // throw Exception() << "D0toKPiPiCLEO::amplitude spin is too high ispin = "
	+ // << ispin << Exception::runerror;
	+ // }
	+ // }
	+ }
	+ if(ichan>=0 && ichan!=int(7+ipi)) continue;
	+ // finally the scalar piece
	+ static const Energy2 sc=2.*GeV2;
	+ ublas::vector<Complex> pVector(2);
	+ double fact = 1.-pres.mass2()/KHalf_->poles()[0];
	+ Complex f1=exp(Complex(0.,1.)gamma_[0]/180.Constants::pi)*fact;
	+ double shat = (pres.mass2()-sc)/GeV2;
	+ pVector[0]=0.;
	+ for(unsigned int ix=0;ix<c1_.size();++ix) {
	+ pVector[0] += f1*c1_[ix];
	+ f1 *= shat;
	+ }
	+ pVector[1]=0.;
	+ Complex f2=exp(Complex(0.,1.)gamma_[1]/180.Constants::pi)*fact;
	+ for(unsigned int ix=0;ix<c2_.size();++ix) {
	+ pVector[1] += f2*c2_[ix];
	+ f2 *= shat;
	+ }
	+ ublas::vector<Complex> p3Vector(1);
	+ p3Vector[0]=0.;
	+ Complex f3=exp(Complex(0.,1.)gamma_[2]/180.Constants::pi);
	+ for(unsigned int ix=0;ix<c3_.size();++ix) {
	+ p3Vector[0] += f2*c3_[ix];
	+ f3 *= shat;
	+ }
	+ Complex phase = exp(Complex(0.,1.)theta_/180.Constants::pi);
	+ pVector[0] += Complex(beta_phaseg1_/KHalf_->poles()[0]);
	+ pVector[1] += beta_phaseg2_/KHalf_->poles()[0];
	+ amp += KHalf_->amplitudes(pres.mass2(), pVector,true)[0];
	+ amp += KThreeHalf_->amplitudes(pres.mass2(),p3Vector,true)[0];
	+ }
	+ // now compute the matrix element
	+ (*ME())(0,0,0,0)=amp;
	+ return norm(amp);
	+}
	+
	+void DtoKPiPiFOCUS::dataBaseOutput(ofstream & output, bool header) const {
	+ if(header) output << "update decayers set parameters=\"";
	+ // parameters for the DecayIntegrator base class
	+ DecayIntegrator::dataBaseOutput(output,false);
	+ // // parameters
	+ // output << "newdef " << name() << ":KmPipPipNonResonantMagnitude "
	+ // << _a1NR << "\n";
	+ // output << "newdef " << name() << ":KmPipPipNonResonantPhase "
	+ // << _phi1NR << "\n";
	+ // output << "newdef " << name() << ":KmPipPipK892Magnitude "
	+ // << _a1K892 << "\n";
	+ // output << "newdef " << name() << ":KmPipPipK892Phase "
	+ // << _phi1K892 << "\n";
	+ // output << "newdef " << name() << ":KmPipPipK1430Magnitude "
	+ // << _a1K1430 << "\n";
	+ // output << "newdef " << name() << ":KmPipPipK1430Phase "
	+ // << _phi1K1430 << "\n";
	+ // output << "newdef " << name() << ":KmPipPipK1680Magnitude "
	+ // << _a1K1680 << "\n";
	+ // output << "newdef " << name() << ":KmPipPipK1680Phase "
	+ // << _phi1K1680 << "\n";
	+ // output << "newdef " << name() << ":KmPipPi0NonResonantMagnitude "
	+ // << _a2NR << "\n";
	+ // output << "newdef " << name() << ":KmPipPi0NonResonantPhase "
	+ // << _phi2NR << "\n";
	+ // output << "newdef " << name() << ":KmPipPi0K8920Magnitude "
	+ // << _a2K8920 << "\n";
	+ // output << "newdef " << name() << ":KmPipPi0K8920Phase "
	+ // << _phi2K8920 << "\n";
	+ // output << "newdef " << name() << ":KmPipPi0K892mMagnitude "
	+ // << _a2K892m << "\n";
	+ // output << "newdef " << name() << ":KmPipPi0K892mPhase "
	+ // << _phi2K892m << "\n";
	+ // output << "newdef " << name() << ":KmPipPi0RhoMagnitude "
	+ // << _a2rho << "\n";
	+ // output << "newdef " << name() << ":KmPipPi0RhoPhase "
	+ // << _phi2rho << "\n";
	+ // output << "newdef " << name() << ":K0PipPimNonResonantMagnitude "
	+ // << _a3NR << "\n";
	+ // output << "newdef " << name() << ":K0PipPimNonResonantPhase "
	+ // << _phi3NR << "\n";
	+ // output << "newdef " << name() << ":K0PipPimK892Magnitude "
	+ // << _a3K892 << "\n";
	+ // output << "newdef " << name() << ":K0PipPimK892Phase "
	+ // << _phi3K892 << "\n";
	+ // output << "newdef " << name() << ":K0PipPimRhoMagnitude "
	+ // << _a3rho << "\n";
	+ // output << "newdef " << name() << ":K0PipPimRhoPhase "
	+ // << _phi3rho << "\n";
	+ // output << "newdef " << name() << ":LocalParameters " << _localparameters << "\n";
	+ // output << "newdef " << name() << ":K8920Mass " << _mK8920/GeV << "\n";
	+ // output << "newdef " << name() << ":K8920Width " << _wK8920/GeV << "\n";
	+ // output << "newdef " << name() << ":K892MinusMass " << _mK892m/GeV << "\n";
	+ // output << "newdef " << name() << ":K892MinusWidth " << _wK892m/GeV << "\n";
	+ // output << "newdef " << name() << ":K1680Mass " << _mK1680/GeV << "\n";
	+ // output << "newdef " << name() << ":K1680Width " << _wK1680/GeV << "\n";
	+ // output << "newdef " << name() << ":K1430Mass " << _mK1430/GeV << "\n";
	+ // output << "newdef " << name() << ":K1430Width " << _wK1430/GeV << "\n";
	+ // output << "newdef " << name() << ":Rho0Mass " << _mrho0 /GeV << "\n";
	+ // output << "newdef " << name() << ":Rho0Width " << _wrho0 /GeV << "\n";
	+ // output << "newdef " << name() << ":RhoPlusMass " << _mrhop /GeV << "\n";
	+ // output << "newdef " << name() << ":RhoPlusWidth " << _wrhop /GeV << "\n";
	+ // for(unsigned int ix=0;ix<_maxwgt.size();++ix) {
	+ // output << "insert " << name() << ":MaximumWeights "
	+ // << ix << " " << _maxwgt[ix] << "\n";
	+ // }
	+ // for(unsigned int ix=0;ix<_weights.size();++ix) {
	+ // output << "insert " << name() << ":Weights "
	+ // << ix << " " << _weights[ix] << "\n";
	+ // }
	+ if(header) {
	+ output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl;
	+ }
	+}
	+
	+void DtoKPiPiFOCUS::doinit() {
	+ DecayIntegrator::doinit();
	+ // set up the phase space
	+ // intermediate resonances
	+ tPDPtr resonances[5]={getParticleData(-313 ),
	+ getParticleData(-100313),
	+ getParticleData(-30313 ),
	+ getParticleData(-315 ),
	+ getParticleData(-10311 )};
	+ // D+ -> K-pi+pi+
	+ tPDPtr in = getParticleData(ParticleID::Dplus);
	+ tPDVector out= {getParticleData(ParticleID::Kminus),
	+ getParticleData(ParticleID::piplus),
	+ getParticleData(ParticleID::piplus)};
	+ PhaseSpaceModePtr mode = new_ptr(PhaseSpaceMode(in,out,maxWgt_));
	+ vector<PhaseSpaceChannel> channels;
	+ for(tPDPtr res : resonances) {
	+ channels.push_back((PhaseSpaceChannel(mode),0,res,0,3,1,1,1,2));
	+ channels.push_back((PhaseSpaceChannel(mode),0,res,0,2,1,1,1,3));
	+ }
	+ // add the mode
	+ for(unsigned int ix=0;ix<channels.size();++ix) {
	+ channels[ix].weight(weights_[ix]);
	+ mode->addChannel(channels[ix]);
	+ }
	+ addMode(mode);
	+}
	+
	+void DtoKPiPiFOCUS::doinitrun() {
	+ DecayIntegrator::doinitrun();
	+}
	diff --git a/Decay/ScalarMeson/DtoKPiPiFOCUS.h b/Decay/ScalarMeson/DtoKPiPiFOCUS.h
	new file mode 100644
	--- /dev/null
	+++ b/Decay/ScalarMeson/DtoKPiPiFOCUS.h
	@@ -0,0 +1,246 @@
	+// -- C++ --
	+#ifndef Herwig_DtoKPiPiFOCUS_H
	+#define Herwig_DtoKPiPiFOCUS_H
	+//
	+// This is the declaration of the DtoKPiPiFOCUS class.
	+//
	+
	+#include "Herwig/Decay/DecayIntegrator.h"
	+#include "Herwig/Decay/FormFactors/KMatrix.h"
	+
	+namespace Herwig {
	+
	+using namespace ThePEG;
	+
	+/**
	+ * Here is the documentation of the DtoKPiPiFOCUS class.
	+ *
	+ * @see \ref DtoKPiPiFOCUSInterfaces "The interfaces"
	+ * defined for DtoKPiPiFOCUS.
	+ */
	+class DtoKPiPiFOCUS: public DecayIntegrator {
	+
	+public:
	+
	+ /**
	+ * The default constructor.
	+ */
	+ DtoKPiPiFOCUS();
	+
	+ /**
	+ * Which of the possible decays is required
	+ * @param cc Is this mode the charge conjugate
	+ * @param parent The decaying particle
	+ * @param children The decay products
	+ */
	+ virtual int modeNumber(bool & cc, tcPDPtr parent,
	+ const tPDVector & children) const;
	+
	+ /**
	+ * Return the matrix element squared for a given mode and phase-space channel.
	+ * @param ichan The channel we are calculating the matrix element for.
	+ * @param part The decaying Particle.
	+ * @param outgoing The particles produced in the decay
	+ * @param momenta The momenta of the particles produced in the decay
	+ * @param meopt Option for the calculation of the matrix element
	+ * @return The matrix element squared for the phase-space configuration.
	+ */
	+ double me2(const int ichan,const Particle & part,
	+ const tPDVector & outgoing,
	+ const vector<Lorentz5Momentum> & momenta,
	+ MEOption meopt) const;
	+
	+ /**
	+ * Construct the SpinInfos for the particles produced in the decay
	+ */
	+ virtual void constructSpinInfo(const Particle & part,
	+ ParticleVector outgoing) const;
	+
	+ /**
	+ * Output the setup information for the particle database
	+ * @param os The stream to output the information to
	+ * @param header Whether or not to output the information for MySQL
	+ */
	+ virtual void dataBaseOutput(ofstream & os,bool header) 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 Standard Interfaced functions. */
	+ //@{
	+ /**
	+ * Initialize this object after the setup phase before saving an
	+ * EventGenerator to disk.
	+ * @throws InitException if object could not be initialized properly.
	+ */
	+ virtual void doinit();
	+
	+ /**
	+ * Initialize this object. Called in the run phase just before
	+ * a run begins.
	+ */
	+ virtual void doinitrun();
	+ //@}
	+
	+protected:
	+
	+ /** @name Clone Methods. */
	+ //@{
	+ /**
	+ * Make a simple clone of this object.
	+ * @return a pointer to the new object.
	+ */
	+ virtual IBPtr clone() const;
	+
	+ /** Make a clone of this object, possibly modifying the cloned object
	+ * to make it sane.
	+ * @return a pointer to the new object.
	+ */
	+ virtual IBPtr fullclone() const;
	+ //@}
	+
	+private:
	+
	+ /**
	+ * The assignment operator is private and must never be called.
	+ * In fact, it should not even be implemented.
	+ */
	+ DtoKPiPiFOCUS & operator=(const DtoKPiPiFOCUS &);
	+
	+private:
	+
	+ /**
	+ * Parameters for the Blatt-Weisskopf form-factors
	+ */
	+ //@{
	+ /**
	+ * Radial size for the \f$D^0\f$
	+ */
	+ InvEnergy rD0_;
	+
	+ /**
	+ * Radial size for the light resonances
	+ */
	+ InvEnergy rRes_;
	+ //@}
	+
	+ /**
	+ * The K-matrices for the s-wave component
	+ */
	+ //@{
	+ /**
	+ * The \f$I=\frac12\f$ component
	+ */
	+ KMatrixPtr KHalf_;
	+
	+ /**
	+ * The \f$I=\frac32\f$ component
	+ */
	+ KMatrixPtr KThreeHalf_;
	+ //@}
	+
	+ /**
	+ * Parameters for the f$pf$-vector
	+ */
	+ //@{
	+ /**
	+ * Pole couplings
	+ */
	+ Energy g1_,g2_;
	+
	+ /**
	+ * \f$\beta\f$
	+ */
	+ Energy beta_;
	+
	+ /**
	+ * \f$\theta\f$
	+ */
	+ double theta_;
	+
	+ /**
	+ * \f$\gamma\f$
	+ */
	+ vector<double> gamma_;
	+
	+ /**
	+ * \f$c_{1i}\f$
	+ */
	+ vector<double> c1_;
	+
	+ /**
	+ * \f$c_{2i}\f$
	+ */
	+ vector<double> c2_;
	+
	+ /**
	+ * \f$c_{3i}\f$
	+ */
	+ vector<double> c3_;
	+ //@}
	+
	+ /**
	+ * Parameters for the phase-space integration
	+ */
	+ //@{
	+ /**
	+ * Maximum weights for the various modes
	+ */
	+ double maxWgt_;
	+
	+ /**
	+ * Weights for the different integration channels
	+ */
	+ vector<double> weights_;
	+ //@}
	+
	+ /**
	+ * Masses and widths of the resonances
	+ */
	+ //@{
	+ /**
	+ * Masses
	+ */
	+ vector<Energy> mRes_;
	+
	+ /**
	+ * Widths
	+ */
	+ vector<Energy> wRes_;
	+ //@}
	+
	+ /**
	+ * Spin density matrix
	+ */
	+ mutable RhoDMatrix rho_;
	+
	+};
	+
	+}
	+
	+#endif /* Herwig_DtoKPiPiFOCUS_H */
	diff --git a/Decay/ScalarMeson/DtoKPiPiMarkIII.cc b/Decay/ScalarMeson/DtoKPiPiMarkIII.cc
	--- a/Decay/ScalarMeson/DtoKPiPiMarkIII.cc
	+++ b/Decay/ScalarMeson/DtoKPiPiMarkIII.cc
	@@ -1,680 +1,680 @@
	// -- C++ --
	//
	// DtoKPiPiMarkIII.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2019 The Herwig Collaboration
	//
	// Herwig is licenced under version 3 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 DtoKPiPiMarkIII class.
	//

	#include "DtoKPiPiMarkIII.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/ParVector.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h"
	#include "Herwig/Decay/GeneralDecayMatrixElement.h"

	using namespace Herwig;

	using ThePEG::Helicity::ScalarWaveFunction;
	using ThePEG::Helicity::incoming;
	using ThePEG::Helicity::outgoing;

	Complex DtoKPiPiMarkIII::amplitude(bool rho, Energy mD,
	- Energy mA , Energy mB , Energy mC ,
	- Energy mAB, Energy mAC, Energy mBC,
	- Energy mres, Energy wres) const {
	+ Energy mA , Energy mB , Energy mC ,
	+ Energy mAB, Energy mAC, Energy mBC,
	+ Energy mres, Energy wres) const {
	InvEnergy radius = rho ? _rrho : _rKstar;
	Energy pAB = Kinematics::pstarTwoBodyDecay(mAB ,mA,mB);
	Energy pR = Kinematics::pstarTwoBodyDecay(mres,mA,mB);
	Energy2 mgam = wressqr(mres)/mABMath::powi(pAB/pR,3)*
	(1.+sqr(radiuspR))/(1.+sqr(radiuspAB));
	Energy2 s = (sqr(mAC)-sqr(mBC)-(sqr(mD)-sqr(mC))(sqr(mA)-sqr(mB))/sqr(mres))
	sqrt((1.+sqr(radiuspR))/(1.+sqr(radiuspAB)));
	Complex output=s/((sqr(mres)- sqr(mAB)-complex<Energy2>(ZERO,mgam)));
	return output;
	}

	DtoKPiPiMarkIII::DtoKPiPiMarkIII() {
	// Amplitudes and phases for D0 -> K- pi+ pi0
	_a1rho = 1.0000; _phi1rho = 0.;
	_a1Kstarm = 0.4018; _phi1Kstarm = 154.;
	_a1Kstar0 = 0.4244; _phi1Kstar0 = 7.;
	_a1NR = 2.0693; _phi1NR = 52.;
	// Amplitudes and phases for D0 -> Kbar0 pi+ pi-
	_a2rho = 0.0975; _phi2rho = 93.;
	_a2Kstar = 0.2225; _phi2Kstar = 0.;
	_a2NR = 1.0000; _phi2NR = 0.;
	// Amplitudes and phases for D+ -> Kbar0 pi+ pi0
	_a3rho = 1.0000; _phi3rho = 0.;
	_a3Kstar = 0.5617; _phi3Kstar = 43.;
	_a3NR = 2.7250; _phi3NR = 250.;
	// Amplitudes and phases for D+ -> K- pi+ pi+
	_a4Kstar = 0.04749; _phi4Kstar = 105.;
	_a4NR = 1.00000; _phi4NR = 0.;
	// masses and widths of the resonances
	_localparameters = true;
	_mrhop = 0.770GeV; _wrhop = 0.1533GeV;
	_mrho0 = 0.770GeV; _wrho0 = 0.1533GeV;
	_mKstarm = 0.8921GeV; _wKstarm = 0.0511GeV;
	_mKstar0 = 0.8695GeV; _wKstar0 = 0.0502GeV;
	// radii of the mesons
	_rrho = 5.mm1e-12/hbarc;
	_rKstar = 2.mm1e-12/hbarc;
	// intermediates
	generateIntermediates(true);
	}


	void DtoKPiPiMarkIII::persistentOutput(PersistentOStream & os) const {
	os << _a1rho << _phi1rho << _a1Kstarm << _phi1Kstarm << _a1Kstar0 << _phi1Kstar0
	<< _a1NR << _phi1NR << _c1rho << _c1Kstarm << _c1Kstar0 << _c1NR << _a2rho
	<< _phi2rho << _a2Kstar << _phi2Kstar << _a2NR << _phi2NR << _c2rho << _c2Kstar
	<< _c2NR << _a3rho << _phi3rho << _a3Kstar << _phi3Kstar << _a3NR << _phi3NR
	<< _c3rho << _c3Kstar << _c3NR << _a4Kstar << _phi4Kstar << _a4NR << _phi4NR
	<< _c4Kstar << _c4NR << _localparameters << ounit(_mrhop,GeV) << ounit(_wrhop,GeV)
	<< ounit(_mrho0,GeV) << ounit(_wrho0,GeV) << ounit(_mKstarm,GeV)
	<< ounit(_wKstarm,GeV) << ounit(_mKstar0,GeV) << ounit(_wKstar0,GeV)
	<< ounit(_rrho,1./GeV) << ounit(_rKstar,1./GeV) << _maxwgt << _weights;
	}

	void DtoKPiPiMarkIII::persistentInput(PersistentIStream & is, int) {
	is >> _a1rho >> _phi1rho >> _a1Kstarm >> _phi1Kstarm >> _a1Kstar0 >> _phi1Kstar0
	>> _a1NR >> _phi1NR >> _c1rho >> _c1Kstarm >> _c1Kstar0 >> _c1NR >> _a2rho
	>> _phi2rho >> _a2Kstar >> _phi2Kstar >> _a2NR >> _phi2NR >> _c2rho >> _c2Kstar
	>> _c2NR >> _a3rho >> _phi3rho >> _a3Kstar >> _phi3Kstar >> _a3NR >> _phi3NR
	>> _c3rho >> _c3Kstar >> _c3NR >> _a4Kstar >> _phi4Kstar >> _a4NR >> _phi4NR
	>> _c4Kstar >> _c4NR >> _localparameters >> iunit(_mrhop,GeV) >> iunit(_wrhop,GeV)
	>> iunit(_mrho0,GeV) >> iunit(_wrho0,GeV) >> iunit(_mKstarm,GeV)
	>> iunit(_wKstarm,GeV) >> iunit(_mKstar0,GeV) >> iunit(_wKstar0,GeV)
	>> iunit(_rrho,1./GeV) >> iunit(_rKstar,1./GeV) >> _maxwgt >> _weights;
	}

	// The following static variable is needed for the type
	// description system in ThePEG.
	DescribeClass<DtoKPiPiMarkIII,DecayIntegrator>
	describeHerwigDtoKPiPiMarkIII("Herwig::DtoKPiPiMarkIII", "HwSMDecay.so");

	void DtoKPiPiMarkIII::Init() {

	static ClassDocumentation<DtoKPiPiMarkIII> documentation
	("The DtoKPiPiMarkIII class performs the D -> K pi pi decays using the fit"
	"of the MarkIII collaboration",
	"The fit of the Mark III collaboration \\cite{Adler:1987sd} was used for "
	"$D\\to K\\pi\\pi$ decays",
	"\\bibitem{Adler:1987sd} J.~Adler {\\it et al.} [MARK-III Collaboration], "
	"Phys.\\ Lett.\\ B {\\bf 196} (1987) 107.");

	static Parameter<DtoKPiPiMarkIII,double> interfaceKmPipPi0RhoMagnitude
	("KmPipPi0RhoMagnitude",
	"The magnitude of the rho component for D0 -> K- pi+ pi0",
	&DtoKPiPiMarkIII::_a1rho, 1.00, 0.0, 10.0,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiMarkIII,double> interfaceKmPipPi0RhoPhase
	("KmPipPi0RhoPhase",
	"The phase of the rho component for D0 -> K- pi+ pi0",
	&DtoKPiPiMarkIII::_phi1rho, 0.,0., 360.,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiMarkIII,double> interfaceKmPipPi0KstarmMagnitude
	("KmPipPi0KstarmMagnitude",
	"The magnitude of the K*(892)- component for D0 -> K- pi+ pi0",
	&DtoKPiPiMarkIII::_a1Kstarm, 0.371, 0.0, 10.0,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiMarkIII,double> interfaceKmPipPi0KstarmPhase
	("KmPipPi0KstarmPhase",
	"The phase of the K*(892)- component for D0 -> K- pi+ pi0",
	&DtoKPiPiMarkIII::_phi1Kstarm, 154.,0., 360.,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiMarkIII,double> interfaceKmPipPi0Kstar0Magnitude
	("KmPipPi0Kstar0Magnitude",
	"The magnitude of the K*(892)0 component for D0 -> K- pi+ pi0",
	&DtoKPiPiMarkIII::_a1Kstar0, 0.391, 0.0, 10.0,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiMarkIII,double> interfaceKmPipPi0Kstar0Phase
	("KmPipPi0Kstar0Phase",
	"The phase of the K*(892)0 component for D0 -> K- pi+ pi0",
	&DtoKPiPiMarkIII::_phi1Kstar0, 7.,0., 360.,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiMarkIII,double> interfaceKmPipPi0NonResonantMagnitude
	("KmPipPi0NonResonantMagnitude",
	"The magnitude of the non-resonant component for D0 -> K- pi+ pi0",
	&DtoKPiPiMarkIII::_a1NR, 1.889, 0.0, 10.0,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiMarkIII,double> interfaceKmPipPi0NonResonantPhase
	("KmPipPi0NonResonantPhase",
	"The phase of the non-resonant component for D0 -> K- pi+ pi0",
	&DtoKPiPiMarkIII::_phi1NR, 52.,0., 360.,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiMarkIII,double> interfaceK0PipPimRhoMagnitude
	("K0PipPimRhoMagnitude",
	"The magnitude of the rho component for D0 -> Kbar0 pi+pi-",
	&DtoKPiPiMarkIII::_a2rho, 1.000, 0.0, 10.0,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiMarkIII,double> interfaceK0PipPimRhoPhase
	("K0PipPimRhoPhase",
	"The phase of the rho component for D0 -> Kbar0 pi+pi-",
	&DtoKPiPiMarkIII::_phi2rho, 93.,0., 360.,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiMarkIII,double> interfaceK0PipPimKstarMagnitude
	("K0PipPimKstarMagnitude",
	"The magnitude of the K*(892) component for D0 -> Kbar0 pi+pi-",
	&DtoKPiPiMarkIII::_a2Kstar, 2.106, 0.0, 10.0,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiMarkIII,double> interfaceK0PipPimKstarPhase
	("K0PipPimKstarPhase",
	"The phase of the K*(892)0 component for D0 -> Kbar0 pi+pi-",
	&DtoKPiPiMarkIII::_phi2Kstar, 0.,0., 360.,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiMarkIII,double> interfaceK0PipPimNonResonantMagnitude
	("K0PipPimNonResonantMagnitude",
	"The magnitude of the non-resonant component for D0 -> Kbar0 pi+pi-",
	&DtoKPiPiMarkIII::_a2NR, 9.379, 0.0, 10.0,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiMarkIII,double> interfaceK0PipPimNonResonantPhase
	("K0PipPimNonResonantPhase",
	"The phase of the non-resonant component for D0 -> Kbar0 pi+pi-",
	&DtoKPiPiMarkIII::_phi2NR, 0.,0., 360.,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiMarkIII,double> interfaceK0PipPi0RhoMagnitude
	("K0PipPi0RhoMagnitude",
	"The magnitude of the rho component for D+ -> Kbar0 pi+ pi0",
	&DtoKPiPiMarkIII::_a3rho, 1.0, 0.0, 10.0,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiMarkIII,double> interfaceK0PipPi0RhoPhase
	("K0PipPi0RhoPhase",
	"The phase of the rho component for D+ -> Kbar0 pi+ pi0",
	&DtoKPiPiMarkIII::_phi3rho, 0.0, 0.0, 360.0,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiMarkIII,double> interfaceK0PipPi0KstarMagnitude
	("K0PipPi0KstarMagnitude",
	"The magnitude of the K* component for D+ -> Kbar0 pi+ pi0",
	&DtoKPiPiMarkIII::_a3Kstar, 0.517, 0.0, 10.0,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiMarkIII,double> interfaceK0PipPi0KstarPhase
	("K0PipPi0KstarPhase",
	"The phase of the K* component for D+ -> Kbar0 pi+ pi0",
	&DtoKPiPiMarkIII::_phi3Kstar, 43.0, 0.0, 360.0,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiMarkIII,double> interfaceK0PipPi0NonResonantMagnitude
	("K0PipPi0NonResonantMagnitude",
	"The magnitude of the non-resonant component for D0 -> Kbar0 pi+pi-",
	&DtoKPiPiMarkIII::_a3NR, 2.490, 0.0, 10.0,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiMarkIII,double> interfaceK0PipPi0NonResonantPhase
	("K0PipPi0NonResonantPhase",
	"The phase of the non-resonant component for D0 -> Kbar0 pi+pi-",
	&DtoKPiPiMarkIII::_phi3NR, 250.,0., 360.,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiMarkIII,double> interfaceKmPipPipNonResonantMagnitude
	("KmPipPipNonResonantMagnitude",
	"The magnitude of the non-resonant component for D+ -> K- pi+ pi+",
	&DtoKPiPiMarkIII::_a4NR, 1.00, 0.0, 10.0,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiMarkIII,double> interfaceKmPipPipNonResonantPhase
	("KmPipPipNonResonantPhase",
	"The phase of the non-resonant component for D+ -> K- pi+ pi+",
	&DtoKPiPiMarkIII::_phi4NR, 0., 0., 360.,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiMarkIII,double> interfaceKmPipPipKstarMagnitude
	("KmPipPipKstarMagnitude",
	"The magnitude of the K*(892) component for D+ -> K- pi+ pi+",
	&DtoKPiPiMarkIII::_a4Kstar, 0.049, 0.0, 10.0,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiMarkIII,double> interfaceKmPipPipKstarPhase
	("KmPipPipKstarPhase",
	"The phase of the K*(892) component for D+ -> K- pi+ pi+",
	&DtoKPiPiMarkIII::_phi4Kstar, 105.,0., 360.,
	false, false, Interface::limited);

	static Switch<DtoKPiPiMarkIII,bool> interfaceLocalParameters
	("LocalParameters",
	"Whether to use local values for the masses and widths or"
	" those from the ParticleData objects",
	&DtoKPiPiMarkIII::_localparameters, true, false, false);
	static SwitchOption interfaceLocalParametersLocal
	(interfaceLocalParameters,
	"Local",
	"Use local values",
	true);
	static SwitchOption interfaceLocalParametersParticleData
	(interfaceLocalParameters,
	"ParticleData",
	"Use the values from the ParticleData objects",
	false);

	static Parameter<DtoKPiPiMarkIII,Energy> interfaceKstar0Mass
	("Kstar0Mass",
	"The mass of the K*(892)0",
	&DtoKPiPiMarkIII::_mKstar0, GeV, 0.8965 GeV, ZERO, 10.0GeV,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiMarkIII,Energy> interfaceKstar0Width
	("Kstar0Width",
	"The width of the K*(892)0",
	&DtoKPiPiMarkIII::_wKstar0, GeV, 0.0502GeV, ZERO, 10.0GeV,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiMarkIII,Energy> interfaceKstarMinusMass
	("KstarMinusMass",
	"The mass of the K*(892)-",
	&DtoKPiPiMarkIII::_mKstarm, GeV, 0.8921GeV, ZERO, 10.0GeV,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiMarkIII,Energy> interfaceKstarMinusWidth
	("KstarMinusWidth",
	"The width of the K*(892)-",
	&DtoKPiPiMarkIII::_wKstarm, GeV, 0.0511GeV, ZERO, 10.0GeV,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiMarkIII,Energy> interfaceRho0Mass
	("Rho0Mass",
	"The mass of the rho0",
	&DtoKPiPiMarkIII::_mrho0, GeV, 0.770 GeV, ZERO, 10.0GeV,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiMarkIII,Energy> interfaceRho0Width
	("Rho0Width",
	"The width of the rho0",
	&DtoKPiPiMarkIII::_wrho0, GeV, 0.1533GeV, ZERO, 10.0GeV,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiMarkIII,Energy> interfaceRhoPlusMass
	("RhoPlusMass",
	"The mass of the rho+",
	&DtoKPiPiMarkIII::_mrhop, GeV, 0.770 GeV, ZERO, 10.0GeV,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiMarkIII,Energy> interfaceRhoPlusWidth
	("RhoPlusWidth",
	"The width of the rho+",
	&DtoKPiPiMarkIII::_wrhop, GeV, 0.1533GeV, ZERO, 10.0GeV,
	false, false, Interface::limited);

	static Parameter<DtoKPiPiMarkIII,InvEnergy> interfaceRhoRadius
	("RhoRadius",
	"The radius of the rho for the Blatt-Weisskopf factor",
	&DtoKPiPiMarkIII::_rrho, mm1e-12/hbarc, 5.0mm1e-12/hbarc, 0.0mm1e-12/hbarc, 10.0mm*1e-12/hbarc,
	true, false, Interface::limited);

	static Parameter<DtoKPiPiMarkIII,InvEnergy> interfaceKstarRadius
	("KstarRadius",
	"The radius of the K* for the Blatt-Weisskopf factor",
	&DtoKPiPiMarkIII::_rKstar, mm1e-12/hbarc, 2.0mm1e-12/hbarc, 0.0mm1e-12/hbarc, 10.0mm*1e-12/hbarc,
	true, false, Interface::limited);

	static ParVector<DtoKPiPiMarkIII,double> interfaceMaximumWeights
	("MaximumWeights",
	"The maximum weights for the unweighting of the decays",
	&DtoKPiPiMarkIII::_maxwgt, -1, 1.0, 0.0, 1.0e11,
	false, false, Interface::limited);

	static ParVector<DtoKPiPiMarkIII,double> interfaceWeights
	("Weights",
	"The weights for the different channels for the phase-space integration",
	&DtoKPiPiMarkIII::_weights, -1, 1.0, 0.0, 1.0,
	false, false, Interface::limited);
	}

	void DtoKPiPiMarkIII::doinit() {
	DecayIntegrator::doinit();
	double fact = Constants::pi/180.;
	// amplitudes for D0 -> K- pi+ pi0
	_c1rho = _a1rho Complex(cos(_phi1rho fact),sin(_phi1rho *fact));
	_c1Kstarm = _a1KstarmComplex(cos(_phi1Kstarmfact),sin(_phi1Kstarm*fact));
	_c1Kstar0 = _a1Kstar0Complex(cos(_phi1Kstar0fact),sin(_phi1Kstar0*fact));
	_c1NR = _a1NR Complex(cos(_phi1NR fact),sin(_phi1NR *fact));
	// amplitudes for D0 -> Kbar0 pi+ pi-
	_c2rho = _a2rho Complex(cos(_phi2rho fact),sin(_phi2rho *fact));
	_c2Kstar = _a2Kstar Complex(cos(_phi2Kstar fact),sin(_phi2Kstar *fact));
	_c2NR = _a2NR Complex(cos(_phi2NR fact),sin(_phi2NR *fact));
	// amplitudes for D+ -> Kbar0 pi+ pi0
	_c3rho = _a3rho Complex(cos(_phi3rho fact),sin(_phi3rho *fact));
	_c3Kstar = _a3Kstar Complex(cos(_phi3Kstar fact),sin(_phi3Kstar *fact));
	_c3NR = _a3NR Complex(cos(_phi3NR fact),sin(_phi3NR *fact));
	// amplitudes for D+ -> K- pi+ pi+
	_c4Kstar = _a4Kstar Complex(cos(_phi4Kstar fact),sin(_phi4Kstar *fact));
	_c4NR = _a4NR Complex(cos(_phi4NR fact),sin(_phi4NR *fact));
	// intermediate resonances
	tPDPtr k8920 = getParticleData(ParticleID::Kstarbar0 );
	tPDPtr k892m = getParticleData(ParticleID::Kstarminus);
	tPDPtr rho0 = getParticleData(ParticleID::rho0 );
	tPDPtr rhop = getParticleData(ParticleID::rhoplus );
	// D0 -> K- pi+ pi0
	tPDPtr in =getParticleData(ParticleID::D0);
	tPDVector out={getParticleData(ParticleID::Kminus),
	getParticleData(ParticleID::piplus),
	getParticleData(ParticleID::pi0)};
	if(_maxwgt.empty()) _maxwgt.push_back(1.);
	PhaseSpaceModePtr mode1 = new_ptr(PhaseSpaceMode(in,out,_maxwgt[0]));
	if(rhop)
	mode1->addChannel((PhaseSpaceChannel(mode1),0,rhop,0,1,1,2,1,3));
	if(k892m)
	mode1->addChannel((PhaseSpaceChannel(mode1),0,k892m,0,2,1,1,1,3));
	if(k8920)
	mode1->addChannel((PhaseSpaceChannel(mode1),0,k8920,0,3,1,1,1,2));
	// add the mode
	vector<double> wtemp;
	if(mode1->channels().size()<=_weights.size()) {
	vector<double>::const_iterator wit=_weights.begin();
	wtemp=vector<double>(wit,wit+mode1->channels().size());
	}
	else {
	wtemp=vector<double>(mode1->channels().size(),1./double(mode1->channels().size()));
	}
	mode1->setWeights(wtemp);
	addMode(mode1);
	unsigned int icount = mode1->channels().size();
	// D0 -> Kbar0 pi+ pi-
	in = getParticleData(ParticleID::D0);
	out = {getParticleData(ParticleID::Kbar0),
	getParticleData(ParticleID::piplus),
	getParticleData(ParticleID::piminus)};
	if(_maxwgt.size()<2) _maxwgt.push_back(1.);
	PhaseSpaceModePtr mode2 = new_ptr(PhaseSpaceMode(in,out,_maxwgt[1]));
	if(rho0)
	mode2->addChannel((PhaseSpaceChannel(mode2),0,rho0,0,1,1,2,1,3));
	if(k892m)
	mode2->addChannel((PhaseSpaceChannel(mode2),0,k892m,0,2,1,1,1,3));
	// add the mode
	if(icount+mode2->channels().size()<=_weights.size()) {
	vector<double>::const_iterator wit=_weights.begin()+icount;
	wtemp=vector<double>(wit,wit+mode2->channels().size());
	}
	else {
	wtemp=vector<double>(mode2->channels().size(),1./double(mode2->channels().size()));
	}
	icount += mode2->channels().size();
	mode2->setWeights(wtemp);
	addMode(mode2);
	// D+ -> Kbar0pi+pi0
	in = getParticleData(ParticleID::Dplus);
	out = {getParticleData(ParticleID::Kbar0),
	getParticleData(ParticleID::piplus),
	getParticleData(ParticleID::pi0)};
	if(_maxwgt.size()<3) _maxwgt.push_back(1.);
	PhaseSpaceModePtr mode3 = new_ptr(PhaseSpaceMode(in,out,_maxwgt[2]));
	if(rhop)
	mode3->addChannel((PhaseSpaceChannel(mode3),0,rhop,0,1,1,2,1,3));
	if(k8920)
	mode3->addChannel((PhaseSpaceChannel(mode3),0,k8920,0,2,1,1,1,3));
	// add the mode
	if(icount+mode3->channels().size()<=_weights.size()) {
	vector<double>::const_iterator wit=_weights.begin()+icount;
	wtemp=vector<double>(wit,wit+mode3->channels().size());
	}
	else {
	wtemp=vector<double>(mode3->channels().size(),1./double(mode3->channels().size()));
	}
	icount += mode3->channels().size();
	addMode(mode3);
	// D+ -> K-pi+pi+
	in = getParticleData(ParticleID::Dplus);
	out = {getParticleData(ParticleID::Kminus),
	getParticleData(ParticleID::piplus),
	getParticleData(ParticleID::piplus)};
	if(_maxwgt.size()<4) _maxwgt.push_back(1.);
	PhaseSpaceModePtr mode4 = new_ptr(PhaseSpaceMode(in,out,_maxwgt[3]));
	if(k8920) {
	mode4->addChannel((PhaseSpaceChannel(mode4),0,k8920,0,3,1,1,1,2));
	mode4->addChannel((PhaseSpaceChannel(mode4),0,k8920,0,2,1,1,1,3));
	}
	// add the mode
	if(icount+mode4->channels().size()<=_weights.size()) {
	vector<double>::const_iterator wit=_weights.begin()+icount;
	wtemp=vector<double>(wit,wit+mode4->channels().size());
	}
	else {
	wtemp=vector<double>(mode4->channels().size(),1./double(mode4->channels().size()));
	}
	addMode(mode4);
	// reset the resonance parameters in the integration if needed
	if(_localparameters) {
	resetIntermediate(k8920,_mKstar0,_wKstar0);
	resetIntermediate(k892m,_mKstarm,_wKstarm);
	resetIntermediate(rho0 ,_mrho0 ,_wrho0 );
	resetIntermediate(rhop ,_mrhop ,_wrhop );
	}
	else {
	_mKstar0 = k8920->mass();
	_mKstarm = k892m->mass();
	_mrho0 = rho0 ->mass();
	_mrhop = rhop ->mass();
	_wKstar0 = k8920->width();
	_wKstarm = k892m->width();
	_wrho0 = rho0 ->width();
	_wrhop = rhop ->width();
	}
	}

	int DtoKPiPiMarkIII::modeNumber(bool & cc, tcPDPtr parent,
	const tPDVector & children) const {
	int id0(parent->id());
	// incoming particle must be D0 or D+
	if(abs(id0)!=ParticleID::D0&&abs(id0)!=ParticleID::Dplus) return -1;
	cc = id0<0;
	// must be three decay products
	if(children.size()!=3) return -1;
	tPDVector::const_iterator pit = children.begin();
	unsigned int npip(0),npim(0),nkm(0),nk0(0),npi0(0);
	int id;
	for( ;pit!=children.end();++pit) {
	id=(**pit).id();
	if(id ==ParticleID::piplus) ++npip;
	else if(id ==ParticleID::pi0) ++npi0;
	else if(id ==ParticleID::piminus) ++npim;
	else if(abs(id)==ParticleID::K0) ++nk0;
	else if(id ==ParticleID::K_L0) ++nk0;
	else if(id ==ParticleID::K_S0) ++nk0;
	else if(abs(id)==ParticleID::Kplus) ++nkm;
	}
	if(abs(id0)==ParticleID::Dplus) {
	if((id0==ParticleID::Dplus &&(nkm==1&&npip==2))\|\|
	(id0==ParticleID::Dminus&&(nkm==1&&npim==2)))
	return 3;
	else if((id0==ParticleID::Dplus &&nk0==1&&npip==1&&npi0==1)\|\|
	(id0==ParticleID::Dminus&&nk0==1&&npim==1&&npi0==1))
	return 2;
	else
	return -1;
	}
	else {
	if(npim==1&&npip==1&&nk0==1) return 1;
	else if(nkm==1&&(npip+npim)==1&&npi0==1) return 0;
	else return -1;
	}
	}

	void DtoKPiPiMarkIII::
	constructSpinInfo(const Particle & part, ParticleVector decay) const {
	// set up the spin information for the decay products
	ScalarWaveFunction::constructSpinInfo(const_ptr_cast<tPPtr>(&part),
	incoming,true);
	for(unsigned int ix=0;ix<3;++ix)
	ScalarWaveFunction::constructSpinInfo(decay[ix],outgoing,true);
	}

	double DtoKPiPiMarkIII::me2(const int ichan, const Particle & part,
	const tPDVector & ,
	const vector<Lorentz5Momentum> & momenta,
	MEOption meopt) const {
	if(!ME())
	ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin0,PDT::Spin0,PDT::Spin0,PDT::Spin0)));
	useMe();
	if(meopt==Initialize) {
	ScalarWaveFunction::
	calculateWaveFunctions(_rho,const_ptr_cast<tPPtr>(&part),incoming);
	}
	// compute the invariant masses needed to calulate the amplitudes
	Energy mA = momenta[0].mass();
	Energy mB = momenta[1].mass();
	Energy mC = momenta[2].mass();
	Energy mD = part.mass();
	Energy mAB = (momenta[0]+momenta[1]).m();
	Energy mAC = (momenta[0]+momenta[2]).m();
	Energy mBC = (momenta[1]+momenta[2]).m();
	// compute the amplitudes for the resonaces present in both models
	Complex amp;
	// D0 -> K- pi+ pi0
	if(imode()==0) {
	if(ichan<0) {
	amp = _c1NR
	-_c1rho *amplitude(true ,mD,mB,mC,mA,mBC,mAB,mAC,_mrhop ,_wrhop )
	-_c1Kstarm*amplitude(false,mD,mA,mC,mB,mAC,mAB,mBC,_mKstarm,_wKstarm)
	-_c1Kstar0*amplitude(false,mD,mA,mB,mC,mAB,mAC,mBC,_mKstar0,_wKstar0);
	}
	else if(ichan==0) {
	amp = _c1rho *amplitude(true ,mD,mB,mC,mA,mBC,mAB,mAC,_mrhop ,_wrhop );
	}
	else if(ichan==1) {
	amp = _c1Kstarm*amplitude(false,mD,mA,mC,mB,mAC,mAB,mBC,_mKstarm,_wKstarm);
	}
	else {
	amp = _c1Kstar0*amplitude(false,mD,mA,mB,mC,mAB,mAC,mBC,_mKstar0,_wKstar0);
	}
	}
	// D0 -> Kbar0 pi+pi-
	else if(imode()==1) {
	if(ichan<0) {
	amp = _c2NR(-1.+2.UseRandom::rndbool())
	-_c2rho *amplitude(true ,mD,mB,mC,mA,mBC,mAB,mAC,_mrho0 ,_wrho0 )
	+_c2Kstar*amplitude(false,mD,mA,mC,mB,mAC,mAB,mBC,_mKstarm,_wKstarm);
	}
	else if(ichan==0) {
	amp =_c2rho *amplitude(true ,mD,mB,mC,mA,mBC,mAB,mAC,_mrho0 ,_wrho0 );
	}
	else {
	amp =_c2Kstar*amplitude(false,mD,mA,mC,mB,mAC,mAB,mBC,_mKstarm,_wKstarm);
	}
	}
	// D+ -> Kbar0 pi+ pi0
	else if(imode()==2) {
	if(ichan<0) {
	amp = _c3NR
	-_c3rho *amplitude(true ,mD,mB,mC,mA,mBC,mAB,mAC,_mrhop ,_wrhop )
	-_c3Kstar*amplitude(false,mD,mA,mC,mB,mAC,mAB,mBC,_mKstar0,_wKstar0);
	}
	else if(ichan==0) {
	amp = _c3rho *amplitude(true ,mD,mB,mC,mA,mBC,mAB,mAC,_mrhop,_wrhop);
	}
	else {
	amp = _c3Kstar*amplitude(false,mD,mA,mC,mB,mAC,mAB,mBC,_mKstar0,_wKstar0);
	}
	}
	// D+ -> K- pi+ pi+
	else {
	if(ichan<0) {
	amp = _c4NR
	+_c4Kstar*amplitude(false,mD,mA,mB,mC,mAB,mAC,mBC,_mKstar0,_wKstar0)
	+_c4Kstar*amplitude(false,mD,mA,mC,mB,mAC,mAB,mBC,_mKstar0,_wKstar0);
	}
	else if(ichan==0) {
	amp =_c4Kstar*amplitude(false,mD,mA,mB,mC,mAB,mAC,mBC,_mKstar0,_wKstar0);
	}
	else {
	amp =_c4Kstar*amplitude(false,mD,mA,mC,mB,mAC,mAB,mBC,_mKstar0,_wKstar0);
	}
	}
	// now compute the matrix element
	(*ME())(0,0,0,0)=amp;
	return norm(amp);
	}

	void DtoKPiPiMarkIII::dataBaseOutput(ofstream & output, bool header) const {
	if(header) output << "update decayers set parameters=\"";
	// parameters for the DecayIntegrator base class
	DecayIntegrator::dataBaseOutput(output,false);
	// parameters
	output << "newdef " << name() << ":KmPipPi0RhoMagnitude " << _a1rho << "\n";
	output << "newdef " << name() << ":KmPipPi0RhoPhase " << _phi1rho << "\n";
	output << "newdef " << name() << ":KmPipPi0KstarmMagnitude " << _a1Kstarm << "\n";
	output << "newdef " << name() << ":KmPipPi0KstarmPhase " << _phi1Kstarm << "\n";
	output << "newdef " << name() << ":KmPipPi0Kstar0Magnitude " << _a1Kstar0 << "\n";
	output << "newdef " << name() << ":KmPipPi0Kstar0Phase " << _phi1Kstar0 << "\n";
	output << "newdef " << name() << ":KmPipPi0NonResonantMagnitude " << _a1NR << "\n";
	output << "newdef " << name() << ":KmPipPi0NonResonantPhase " << _phi1NR << "\n";
	output << "newdef " << name() << ":K0PipPimRhoMagnitude " << _a2rho << "\n";
	output << "newdef " << name() << ":K0PipPimRhoPhase " << _phi2rho << "\n";
	output << "newdef " << name() << ":K0PipPimKstarMagnitude " << _a2Kstar << "\n";
	output << "newdef " << name() << ":K0PipPimKstarPhase " << _phi2Kstar << "\n";
	output << "newdef " << name() << ":K0PipPimNonResonantMagnitude " << _a2NR << "\n";
	output << "newdef " << name() << ":K0PipPimNonResonantPhase " << _phi2NR << "\n";
	output << "newdef " << name() << ":K0PipPi0RhoMagnitude " << _a3rho << "\n";
	output << "newdef " << name() << ":K0PipPi0RhoPhase " << _phi3rho << "\n";
	output << "newdef " << name() << ":K0PipPi0KstarMagnitude " << _a3Kstar << "\n";
	output << "newdef " << name() << ":K0PipPi0KstarPhase " << _phi3Kstar << "\n";
	output << "newdef " << name() << ":K0PipPi0NonResonantMagnitude " << _a3NR << "\n";
	output << "newdef " << name() << ":K0PipPi0NonResonantPhase " << _phi3NR << "\n";
	output << "newdef " << name() << ":KmPipPipNonResonantMagnitude " << _a4NR << "\n";
	output << "newdef " << name() << ":KmPipPipNonResonantPhase " << _phi4NR << "\n";
	output << "newdef " << name() << ":KmPipPipKstarMagnitude " << _a4Kstar << "\n";
	output << "newdef " << name() << ":KmPipPipKstarPhase " << _phi4Kstar << "\n";
	output << "newdef " << name() << ":LocalParameters " << _localparameters << "\n";
	output << "newdef " << name() << ":Kstar0Mass " << _mKstar0/GeV << "\n";
	output << "newdef " << name() << ":Kstar0Width " << _wKstar0/GeV << "\n";
	output << "newdef " << name() << ":KstarMinusMass " << _mKstarm/GeV << "\n";
	output << "newdef " << name() << ":KstarMinusWidth " << _wKstarm/GeV << "\n";
	output << "newdef " << name() << ":Rho0Mass " << _mrho0/GeV << "\n";
	output << "newdef " << name() << ":Rho0Width " << _wrho0/GeV << "\n";
	output << "newdef " << name() << ":RhoPlusMass " << _mrhop/GeV << "\n";
	output << "newdef " << name() << ":RhoPlusWidth " << _wrhop/GeV << "\n";
	output << "newdef " << name() << ":RhoRadius " << _rrho/(mm*1e-12/hbarc) << "\n";
	output << "newdef " << name() << ":KstarRadius " << _rKstar/(mm*1e-12/hbarc) << "\n";
	for(unsigned int ix=0;ix<_maxwgt.size();++ix) {
	output << "insert " << name() << ":MaximumWeights "
	<< ix << " " << _maxwgt[ix] << "\n";
	}
	for(unsigned int ix=0;ix<_weights.size();++ix) {
	output << "insert " << name() << ":Weights "
	<< ix << " " << _weights[ix] << "\n";
	}
	if(header) {
	output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl;
	}
	}

	void DtoKPiPiMarkIII::doinitrun() {
	DecayIntegrator::doinitrun();
	_weights.resize(mode(0)->channels().size()+mode(1)->channels().size()+
	mode(2)->channels().size()+mode(3)->channels().size());
	_maxwgt.resize(4);
	unsigned int iy=0;
	for(unsigned int ix=0;ix<4;++ix) {
	_maxwgt[ix]=mode(ix)->maxWeight();
	for(unsigned int iz=0;iz<mode(ix)->channels().size();++iz) {
	_weights[iy]=mode(ix)->channels()[ix].weight();
	++iy;
	}
	}
	}

	diff --git a/Decay/ScalarMeson/DtoKPiPiMarkIII.h b/Decay/ScalarMeson/DtoKPiPiMarkIII.h
	--- a/Decay/ScalarMeson/DtoKPiPiMarkIII.h
	+++ b/Decay/ScalarMeson/DtoKPiPiMarkIII.h
	@@ -1,452 +1,452 @@
	// -- C++ --
	//
	// DtoKPiPiMarkIII.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2019 The Herwig Collaboration
	//
	// Herwig is licenced under version 3 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef HERWIG_DtoKPiPiMarkIII_H
	#define HERWIG_DtoKPiPiMarkIII_H
	//
	// This is the declaration of the DtoKPiPiMarkIII class.
	//

	#include "Herwig/Decay/DecayIntegrator.h"
	#include "Herwig/Decay/PhaseSpaceMode.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* Here is the documentation of the DtoKPiPiMarkIII class.
	*
	* @see \ref DtoKPiPiMarkIIIInterfaces "The interfaces"
	* defined for DtoKPiPiMarkIII.
	*/
	class DtoKPiPiMarkIII: public DecayIntegrator {

	public:

	/**
	* The default constructor.
	*/
	DtoKPiPiMarkIII();

	/**
	* Which of the possible decays is required
	* @param cc Is this mode the charge conjugate
	* @param parent The decaying particle
	* @param children The decay products
	*/
	virtual int modeNumber(bool & cc, tcPDPtr parent,
	const tPDVector & children) const;

	/**
	* Return the matrix element squared for a given mode and phase-space channel.
	* @param ichan The channel we are calculating the matrix element for.
	* @param part The decaying Particle.
	* @param outgoing The particles produced in the decay
	* @param momenta The momenta of the particles produced in the decay
	* @param meopt Option for the calculation of the matrix element
	* @return The matrix element squared for the phase-space configuration.
	*/
	double me2(const int ichan,const Particle & part,
	const tPDVector & outgoing,
	const vector<Lorentz5Momentum> & momenta,
	MEOption meopt) const;

	/**
	* Construct the SpinInfos for the particles produced in the decay
	*/
	virtual void constructSpinInfo(const Particle & part,
	ParticleVector outgoing) const;

	/**
	* Output the setup information for the particle database
	* @param os The stream to output the information to
	* @param header Whether or not to output the information for MySQL
	*/
	virtual void dataBaseOutput(ofstream & os,bool header) 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:

	/**
	* Calculate the amplitude for a resonance
	* @param rho True for rho resonances and false for \f$K^*\f$
	* @param mD The mass of the decaying particle
	* @param mA The mass of the first decay product
	* @param mB The mass of the second decay product
	* @param mC The mass of the third decay product
	* @param mAB The mass of the pair AB
	* @param mAC The mass of the pair AC
	* @param mBC The mass of the pair BC
	* @param mres The on-shell mass of the intermediate resonance
	* @param wres The width of the intermediate resonance
	*/
	Complex amplitude(bool rho, Energy mD,
	- Energy mA , Energy mB , Energy mC ,
	- Energy mAB, Energy mAC, Energy mBC,
	- Energy mres, Energy wres) const;
	+ Energy mA , Energy mB , Energy mC ,
	+ Energy mAB, Energy mAC, Energy mBC,
	+ Energy mres, Energy wres) const;

	protected:

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

	/** 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 {return new_ptr(*this);}
	//@}

	protected:

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

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

	private:

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

	private:

	/**
	* Amplitudes and phases for \f$D^0\to K^-\pi^+\pi^0\f$
	*/
	//@{
	/**
	* Magnitude of the \f$\rho\f$ component
	*/
	double _a1rho;

	/**
	* Phase of the \f$\rho\f$ component
	*/
	double _phi1rho;

	/**
	* Magnitude of the \f$K^{*-}\f$ component
	*/
	double _a1Kstarm;

	/**
	* Phase of the \f$K^{*-}\f$ component
	*/
	double _phi1Kstarm;

	/**
	* Magnitude of the \f$\bar{K}^{*0}\f$ component
	*/
	double _a1Kstar0;

	/**
	* Phase of the \f$\bar{K}^{*0}\f$ component
	*/
	double _phi1Kstar0;

	/**
	* Magnitude of the non-resonant component
	*/
	double _a1NR;

	/**
	* Phase of the non-resonant component
	*/
	double _phi1NR;

	/**
	* Magnitude of the \f$\rho\f$ component
	*/
	Complex _c1rho;

	/**
	* Magnitude of the \f$K^{*-}\f$ component
	*/
	Complex _c1Kstarm;

	/**
	* Magnitude of the \f$\bar{K}^{*0}\f$ component
	*/
	Complex _c1Kstar0;

	/**
	* Magnitude of the non-resonant component
	*/
	Complex _c1NR;
	//@}

	/**
	* Amplitudes and phases for \f$D^0\to \bar{K}^0\pi^+\pi^-\f$
	*/
	//@{
	/**
	* Magnitude of the \f$\rho^0\f$ component
	*/
	double _a2rho;

	/**
	* Phase of the \f$\rho^0\f$ component
	*/
	double _phi2rho;

	/**
	* Magnitude of the \f$K^{*-}\f$ component
	*/
	double _a2Kstar;

	/**
	* Phase of the \f$K^{*-}\f$ component
	*/
	double _phi2Kstar;

	/**
	* Magnitude of the non-resonant component
	*/
	double _a2NR;

	/**
	* Phase of the non-resonant component
	*/
	double _phi2NR;

	/**
	* Amplitude of the \f$\rho^0\f$ component
	*/
	Complex _c2rho;

	/**
	* Amplitude of the \f$K^{*-}\f$ component
	*/
	Complex _c2Kstar;

	/**
	* Amplitude of the non-resonant component
	*/
	Complex _c2NR;
	//@}

	/**
	* Amplitudes and phases for \f$D^+\to \bar{K}^0\pi^+\pi^0\f$
	*/
	//@{
	/**
	* Magnitude of the \f$\rho\f$ component
	*/
	double _a3rho;

	/**
	* Phase of the \f$\rho\f$ component
	*/
	double _phi3rho;

	/**
	* Magnitude of the \f$\bar{K}^{*0}\f$ component
	*/
	double _a3Kstar;

	/**
	* Phase of the \f$\bar{K}^{*0}\f$ component
	*/
	double _phi3Kstar;

	/**
	* Magnitude of the non-resonant component
	*/
	double _a3NR;

	/**
	* Phase of the non-resonant component
	*/
	double _phi3NR;

	/**
	* Amplitude of the \f$\rho\f$ component
	*/
	Complex _c3rho;

	/**
	* Amplitude of the \f$\bar{K}^{*0}\f$ component
	*/
	Complex _c3Kstar;

	/**
	* Amplitude of the non-resonant component
	*/
	Complex _c3NR;
	//@}

	/**
	* Amplitudes and phases for \f$D^+\to K^-\pi^+\pi^+\f$
	*/
	//@{
	/**
	* Magnitude of the \f$\bar{K}^{*0}\f$ component
	*/
	double _a4Kstar;

	/**
	* Phase of the \f$\bar{K}^{*0}\f$ component
	*/
	double _phi4Kstar;

	/**
	* Magnitude of the non-resonant component
	*/
	double _a4NR;

	/**
	* Phase of the non-resonant component
	*/
	double _phi4NR;

	/**
	* Amplitude of the \f$\bar{K}^{*0}\f$ component
	*/
	Complex _c4Kstar;

	/**
	* Amplitude of the non-resonant component
	*/
	Complex _c4NR;
	//@}

	/**
	* Masses and Widths of the resonances
	*/
	//@{
	/**
	* Use local values of the masses and widths
	*/
	bool _localparameters;

	/**
	* Mass of the \f$\rho^+\f$
	*/
	Energy _mrhop;

	/**
	* Width of the \f$\rho^+\f$
	*/
	Energy _wrhop;

	/**
	* Mass of the \f$\rho^0\f$
	*/
	Energy _mrho0;

	/**
	* Width of the \f$\rho^0\f$
	*/
	Energy _wrho0;

	/**
	* Mass of the \f$K^{*-}\f$
	*/
	Energy _mKstarm;

	/**
	* Width of the \f$K^{*-}\f$
	*/
	Energy _wKstarm;

	/**
	* Mass of the \f$\bar{K}^{*0}\f$
	*/
	Energy _mKstar0;

	/**
	* Width of the \f$\bar{K}^{*0}\f$
	*/
	Energy _wKstar0;
	//@}

	/**
	* The radii of the mesons for the form-factors
	*/
	//@{
	/**
	* \f$\rho\f$ radius
	*/
	InvEnergy _rrho;

	/**
	* \f$K^*\f$ radius
	*/
	InvEnergy _rKstar;
	//@}

	/**
	* Parameters for the phase-space integration
	*/
	//@{
	/**
	* Maximum weights for the different modes
	*/
	vector<double> _maxwgt;

	/**
	* Weights for the different phase-space channels
	*/
	vector<double> _weights;
	//@}

	/**
	* Spin density matrix
	*/
	mutable RhoDMatrix _rho;
	};

	}

	#endif /* HERWIG_DtoKPiPiMarkIII_H */
	diff --git a/Decay/ScalarMeson/Makefile.am b/Decay/ScalarMeson/Makefile.am
	--- a/Decay/ScalarMeson/Makefile.am
	+++ b/Decay/ScalarMeson/Makefile.am
	@@ -1,46 +1,48 @@
	BUILT_SOURCES = SMDecayer__all.cc
	CLEANFILES = SMDecayer__all.cc

	SMDecayer__all.cc : $(DIR_H_FILES) $(DIR_CC_FILES) Makefile
	@echo "Concatenating .cc files into $@"
	@$(top_srcdir)/cat_with_cpplines $(DIR_CC_FILES) > $@

	EXTRA_DIST = $(ALL_H_FILES) $(ALL_CC_FILES)

	DIR_H_FILES = $(addprefix $(srcdir)/,$(ALL_H_FILES))
	ALL_H_FILES = \
	EtaPiGammaGammaDecayer.h\
	EtaPiPiGammaDecayer.h \
	EtaPiPiPiDecayer.h \
	PScalar4FermionsDecayer.h\
	PScalarLeptonNeutrinoDecayer.h\
	PScalarPScalarVectorDecayer.h \
	PScalarVectorFermionsDecayer.h\
	PScalarVectorVectorDecayer.h\
	ScalarMesonTensorScalarDecayer.h\
	ScalarScalarScalarDecayer.h \
	SemiLeptonicScalarDecayer.h \
	ScalarMesonFactorizedDecayer.h \
	ScalarVectorVectorDecayer.h \
	DtoKPiPiCLEO.h \
	DtoKPiPiE691.h \
	+DtoKPiPiFOCUS.h \
	DtoKPiPiMarkIII.h

	DIR_CC_FILES = $(addprefix $(srcdir)/,$(ALL_CC_FILES))
	ALL_CC_FILES = \
	EtaPiGammaGammaDecayer.cc \
	EtaPiPiGammaDecayer.cc \
	EtaPiPiPiDecayer.cc \
	PScalar4FermionsDecayer.cc \
	PScalarLeptonNeutrinoDecayer.cc \
	PScalarPScalarVectorDecayer.cc \
	PScalarVectorFermionsDecayer.cc \
	PScalarVectorVectorDecayer.cc \
	ScalarMesonTensorScalarDecayer.cc \
	ScalarScalarScalarDecayer.cc \
	SemiLeptonicScalarDecayer.cc \
	ScalarMesonFactorizedDecayer.cc \
	ScalarVectorVectorDecayer.cc \
	DtoKPiPiCLEO.cc \
	DtoKPiPiE691.cc\
	+DtoKPiPiFOCUS.cc\
	DtoKPiPiMarkIII.cc

File Metadata

Mime Type: text/x-diff
Expires: Sat, May 3, 6:39 AM (1 d, 21 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4980222
Default Alt Text: (105 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions