No OneTemporary
Actions

Size

123 KB

Subscribers

None

View Options

	diff --git a/Decay/ResonanceHelpers.h b/Decay/ResonanceHelpers.h
	--- a/Decay/ResonanceHelpers.h
	+++ b/Decay/ResonanceHelpers.h
	@@ -1,238 +1,241 @@
	// -- C++ --
	//
	// ResonanceHelpers.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2018 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_ResonanceHelpers_H
	#define HERWIG_ResonanceHelpers_H

	namespace Resonance {

	/**
	* The velocity squared
	*/
	inline double beta2(const Energy2 & s, const Energy & m1, const Energy & m2) {
	return max(0.,(1.-sqr(m1+m2)/s)*(1.-sqr(m1-m2)/s));
	}

	/**
	* The velocity
	*/
	inline double beta(const Energy2 & s, const Energy & m1, const Energy & m2) {
	return sqrt(beta2(s,m1,m2));
	}

	/**
	* The derivative of the function \f$\hat{H}(s)\f$ for the GS Breit-Wigner evaluated
	* at the resonance mass
	*/
	inline double dHhatds(const Energy & mRes, const Energy & gamma,
	const Energy & m1, const Energy & m2) {
	double v2 = beta2(sqr(mRes),m1,m2);
	double v = sqrt(v2);
	double r = (sqr(m1) + sqr(m2))/sqr(mRes);
	return gamma/Constants::pi/mRes/v2*
	((3.-2.v2- 3.r)log((1.+v)/(1.-v)) + 2.v*(1.- r/(1.-v2)));
	}

	/**
	* The \f$\hat{H}(s)\f$ function for the GS Breit-Wigner
	*/
	inline Energy2 Hhat(const Energy2 & s, const Energy & mRes, const Energy & gamma,
	const Energy & m1, const Energy & m2) {
	double vR = beta(sqr(mRes),m1,m2);
	double v = beta( s ,m1,m2);
	return gamma/mRes/Constants::pispow(v/vR,3)*log((1.+v)/(1.-v));
	}

	/**
	* The \f$H(s)\f$ function for the GS Breit-Wigner (input the derivative term)
	*/
	inline Energy2 H(const Energy2 & s, const Energy & mRes, const Energy & gamma,
	const Energy & m1, const Energy & m2,
	const double & dH, const Energy2 & Hres) {
	if(s!=ZERO)
	return Hhat(s,mRes,gamma,m1,m2) - Hres - (s-sqr(mRes))*dH;
	else
	return -2.sqr(m1+m2)/Constants::pigamma/mRes/pow(beta(sqr(mRes),m1,m2),3) - Hres + sqr(mRes)*dH;
	}

	/**
	* The \f$H(s)\f$ function for the GS Breit-Wigner (with calc of derivative term)
	*/
	inline Energy2 H(const Energy2 & s, const Energy & mRes, const Energy & gamma,
	const Energy & m1, const Energy & m2) {
	double dH = dHhatds(mRes,gamma,m1,m2);
	Energy2 Hres = Hhat(sqr(mRes),mRes,gamma,m1,m2);
	return H(s,mRes,gamma,m1,m2,dH,Hres);
	}

	/**
	* The \f$p\f$-wave runningwidth
	*/
	inline Energy gammaP(const Energy2 & s, const Energy & mRes, const Energy & gamma,
	const Energy & m1, const Energy & m2) {
	double v2 = beta2(s,m1,m2);
	if(v2<=0.) return ZERO;
	double vR2 = beta2(sqr(mRes),m1,m2);
	double rp = sqrt(v2/vR2);
	return sqrt(s)/mRespow(rp,3)gamma;
	}

	/**
	* The \f$p\f$-wave runningwidth
	*/
	inline Energy gammaD(const Energy2 & s, const Energy & mRes, const Energy & gamma,
	const Energy & m1, const Energy & m2) {
	double v2 = beta2(s,m1,m2);
	if(v2<=0.) return ZERO;
	double vR2 = beta2(sqr(mRes),m1,m2);
	double rp = sqrt(v2/vR2);
	return pow(sqrt(s)/mRes,3)pow(rp,5)gamma;
	}

	/**
	* The \f$p\f$-wave runningwidth
	*/
	inline Energy gammaS(const Energy2 & s, const Energy & mRes, const Energy & gamma,
	const Energy & m1, const Energy & m2) {
	double v2 = beta2(s,m1,m2);
	if(v2<=0.) return ZERO;
	double vR2 = beta2(sqr(mRes),m1,m2);
	double rp = sqrt(v2/vR2);
	return mRes/sqrt(s)rpgamma;
	}

	/**
	* The GS form of the Breit-Wigner distribution
	*/
	inline Complex BreitWignerGS(const Energy2 & s, const Energy & mRes, const Energy & gamma,
	const Energy & m1, const Energy & m2,
	const Energy2 & H0, const double &dH, const Energy2 & Hres) {
	Energy2 mR2=sqr(mRes);
	return (mR2+H0)/(mR2-s+H(s,mRes,gamma,m1,m2,dH,Hres)-Complex(0.,1.)sqrt(s)gammaP(s,mRes,gamma,m1,m2));
	}

	/**
	* The GS form of the Breit-Wigner distribution
	*/
	inline Complex BreitWignerGS(const Energy2 & s, const Energy & mRes,
	const Energy & gamma,
	const Energy & m1, const Energy & m2) {
	double dH = dHhatds(mRes,gamma,m1,m2);
	Energy2 Hres = Hhat(sqr(mRes),mRes,gamma,m1,m2);
	Energy2 H0 = H(ZERO,mRes,gamma,m1,m2,dH,Hres);
	return BreitWignerGS(s,mRes,gamma,m1,m2,H0,dH,Hres);
	}

	/**
	* Standard \f$p\f$-wave Breit-Wigner
	*/
	inline Complex BreitWignerPWave(const Energy2 & s, const Energy & mRes, const Energy & gamma,
	const Energy & m1, const Energy & m2) {
	Energy2 mR2=sqr(mRes);
	return mR2/(mR2-s-Complex(0.,1.)sqrt(s)gammaP(s,mRes,gamma,m1,m2));
	}

	/**
	* Standard \f$s\f$-wave Breit-Wigner
	*/
	inline Complex BreitWignerSWave(const Energy2 & s, const Energy & mRes, const Energy & gamma,
	const Energy & m1, const Energy & m2) {
	Energy2 mR2=sqr(mRes);
	return mR2/(mR2-s-Complex(0.,1.)sqrt(s)gammaS(s,mRes,gamma,m1,m2));
	}

	/**
	* Standard \f$p\f$-wave Breit-Wigner
	*/
	inline Complex BreitWignerDWave(const Energy2 & s, const Energy & mRes, const Energy & gamma,
	const Energy & m1, const Energy & m2) {
	Energy2 mR2=sqr(mRes);
	return mR2/(mR2-s-Complex(0.,1.)sqrt(s)gammaD(s,mRes,gamma,m1,m2));
	}

	/**
	* Standard fixed width Breit-Wigner, no width in numerator
	*/
	inline Complex BreitWignerFW(const Energy2 & s, const Energy & mRes, const Energy & gamma) {
	Energy2 mR2=sqr(mRes);
	return mR2/(mR2-s-Complex(0.,1.)mResgamma);
	}

	/**
	* Standard fixed width Breit-Wigner, width in numerator
	*/
	inline Complex BreitWignerFW_GN(const Energy2 & s, const Energy & mRes, const Energy & gamma) {
	Energy2 mR2=sqr(mRes);
	complex<Energy2> fact = mR2 - Complex(0.,1.)mResgamma;
	return fact/(fact-s);
	}

	/**
	* The \f$H\f$ function from 0512180
	*/
	Complex H(const Energy & mass, const Energy & width, const Energy2 & sp, const Energy2 & sm,
	const Energy2 & s0, const Energy & mp, const Energy & m0) {
	return
	Resonance::BreitWignerPWave(sp,mass,width,mp,m0)+
	Resonance::BreitWignerPWave(sm,mass,width,mp,m0)+
	Resonance::BreitWignerPWave(s0,mass,width,mp,mp);
	}

	/**
	* Sum over \f$p\f$-wave resonaces
	*/
	template<typename Value>
	Complex F_rho(const Energy2 & s,
	const vector<Value> weights,
	const vector<Energy> & mass,
	const vector<Energy> & width,
	const Energy & m1, const Energy & m2) {
	Value norm(0.);
	Complex output;
	for(unsigned int ix=0;ix<weights.size();++ix) {
	norm += weights[ix];
	output += weights[ix]*
	BreitWignerPWave(s,mass[ix],width[ix],m1,m2);
	}
	return output/norm;
	}

	double ga1(const Energy2 &s) {
	static const Energy mpi=0.13957*GeV;
	- if(s>0.838968432668*GeV2) {
	+ if(s<9.*sqr(mpi)) {
	+ return 0.;
	+ }
	+ else if(s>0.838968432668*GeV2) {
	double Q2 = s/GeV2;
	return 1.623*Q2+10.38-9.32/Q2+0.65/sqr(Q2);
	}
	else {
	double Q2 = (s-9.*sqr(mpi))/GeV2;
	return 4.1Q2sqr(Q2)(1.-3.3Q2+5.8*sqr(Q2));
	}
	}

	/**
	* GS form of the \f$a_1\f$ Breit-Wigner
	*/
	Complex BreitWignera1(const Energy2 & s, const Energy & mRes,
	const Energy & gamma) {
	Energy2 mR2 = sqr(mRes);
	return mR2/(mR2-s-Complex(0.,1)gammamRes*ga1(s)/ga1(mR2));
	}

	/**
	* Difference between two resonances
	*/
	complex<InvEnergy2> BreitWignerDiff(const Energy2 & s,
	const Energy & mRes1, const Energy & gamma1,
	const Energy & mRes2, const Energy & gamma2,
	const Energy & m1, const Energy & m2) {
	return
	BreitWignerPWave(s,mRes1,gamma1,m1,m2)/sqr(mRes1)-
	BreitWignerPWave(s,mRes2,gamma2,m1,m2)/sqr(mRes2);
	}
	}
	#endif
	diff --git a/Decay/WeakCurrents/OneKaonTwoPionCurrent.h b/Decay/WeakCurrents/OneKaonTwoPionCurrent.h
	--- a/Decay/WeakCurrents/OneKaonTwoPionCurrent.h
	+++ b/Decay/WeakCurrents/OneKaonTwoPionCurrent.h
	@@ -1,401 +1,401 @@
	// -- C++ --
	//
	// OneKaonTwoPionCurrent.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2017 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_OneKaonTwoPionCurrent_H
	#define HERWIG_OneKaonTwoPionCurrent_H
	//
	// This is the declaration of the OneKaonTwoPionCurrent class.
	//

	#include "WeakCurrent.h"
	#include "Herwig/Decay/ResonanceHelpers.h"
	#include <numeric>

	namespace Herwig {

	using namespace ThePEG;

	/**
	* The OneKaonTwoPionCurrent class implements the model of M. Finkemeier
	* and E.~Mirkes, Z. Phys. C 69 (1996) 243 [arXiv:hep-ph/9503474],
	* for the weak current for three mesons where at least one of the mesons is
	* a kaon.
	*
	* \ingroup Decay
	*
	* This is the base class for the three meson decays of the weak current.
	* It is designed so that the currents for the following modes can be implemented
	* in classes inheriting from this
	* - \f$ \pi^- \pi^- \pi^+ \f$, (imode=0)
	* - \f$ \pi^0 \pi^0 \pi^- \f$, (imode=1)
	* - \f$ K^- \pi^- K^+ \f$, (imode=2)
	* - \f$ K^0 \pi^- \bar{K}^0\f$, (imode=3)
	* - \f$ K^- \pi^0 K^0 \f$, (imode=4)
	* - \f$ \pi^0 \pi^0 K^- \f$, (imode=5)
	* - \f$ K^- \pi^- \pi^+ \f$, (imode=6)
	* - \f$ \pi^- \bar{K}^0 \pi^0 \f$, (imode=7)
	* - \f$ \pi^- \pi^0 \eta \f$, (imode=8)
	*
	* obviously there are other modes with three pseudoscalar mesons for the decay
	* of the weak current but this model original came from \f$\tau\f$ decay where
	* these are the only modes. However one case which is important is the inclusion
	* of the mixing in the neutral kaon sector for which we include the additional
	* currents
	* - \f$ K^0_S \pi^- K^0_S\f$, (imode=9)
	* - \f$ K^0_L \pi^- K^0_L\f$, (imode=10)
	* - \f$ K^0_S \pi^- K^0_L\f$, (imode=11)
	*
	* In this case the current is given by
	* \f[ J^\mu = \left(g^{\mu\nu}-\frac{q^\mu q^\nu}{q^2}\right)
	* \left[F_1(p_2-p_3)^\mu +F_2(p_3-p_1)^\mu+F_3(p_1-p_2)^\mu\right]
	* +q^\mu F_4
	* +F_5\epsilon^{\mu\alpha\beta\gamma}p_1^\alpha p_2^\beta p_3^\gamma
	* \f]
	* where
	* - \f$p_{1,2,3}\f$ are the momenta of the mesons in the order given above.
	* - \f$F_1,F_2,F_3,F_4,F_5\f$ are the form factors which must be
	* calculated in the calculateFormFactors member which should be implemented
	* in classes inheriting from this.
	*
	* @see WeakCurrent.
	*
	* \author Peter Richardson
	* @see \ref OneKaonTwoPionCurrentInterfaces "The interfaces"
	* defined for OneKaonTwoPionCurrent.
	*/
	class OneKaonTwoPionCurrent: public WeakCurrent {

	public:

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

	/** @name Methods for the construction of the phase space integrator. */
	//@{
	/**
	* Complete the construction of the decay mode for integration.classes inheriting
	* from this one.
	* This method is purely virtual and must be implemented in the classes inheriting
	* from WeakCurrent.
	* @param icharge The total charge of the outgoing particles in the current.
	* @param resonance If specified only include terms with this particle
	* @param Itotal If specified the total isospin of the current
	* @param I3 If specified the thrid component of isospin
	* @param imode The mode in the current being asked for.
	* @param mode The phase space mode for the integration
	* @param iloc The location of the of the first particle from the current in
	* the list of outgoing particles.
	* @param ires The location of the first intermediate for the current.
	* @param phase The prototype phase space channel for the integration.
	* @param upp The maximum possible mass the particles in the current are
	* allowed to have.
	* @return Whether the current was sucessfully constructed.
	*/
	virtual bool createMode(int icharge, tcPDPtr resonance,
	IsoSpin::IsoSpin Itotal, IsoSpin::I3 i3,
	unsigned int imode,PhaseSpaceModePtr mode,
	unsigned int iloc,int ires,
	PhaseSpaceChannel phase, Energy upp );
	//@}


	/**
	* Hadronic current. This method is purely virtual and must be implemented in
	* all classes inheriting from this one.
	* @param resonance If specified only include terms with this particle
	* @param Itotal If specified the total isospin of the current
	* @param I3 If specified the thrid component of isospin
	* @param imode The mode
	* @param ichan The phase-space channel the current is needed for.
	* @param scale The invariant mass of the particles in the current.
	* @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 current.
	*/
	virtual vector<LorentzPolarizationVectorE>
	current(tcPDPtr resonance,
	IsoSpin::IsoSpin Itotal, IsoSpin::I3 i3,
	const int imode, const int ichan,Energy & scale,
	const tPDVector & outgoing,
	const vector<Lorentz5Momentum> & momenta,
	DecayIntegrator::MEOption meopt) const;

	/**
	* Accept the decay. Checks the mesons against the list.
	* @param id The id's of the particles in the current.
	* @return Can this current have the external particles specified.
	*/
	virtual bool accept(vector<int> id);

	/**
	* Return the decay mode number for a given set of particles in the current.
	* Checks the mesons against the list.
	* @param id The id's of the particles in the current.
	* @return The number of the mode
	*/
	virtual unsigned int decayMode(vector<int> id);

	/**
	* The particles produced by the current. This returns the mesons for the mode.
	* @param icharge The total charge of the particles in the current.
	* @param imode The mode for which the particles are being requested
	* @param iq The PDG code for the quark
	* @param ia The PDG code for the antiquark
	* @return The external particles for the current.
	*/
	virtual tPDVector particles(int icharge, unsigned int imode, int iq, int ia);

	/**
	* 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
	* @param create Whether or not to add a statement creating the object
	*/
	virtual void dataBaseOutput(ofstream & os,bool header,bool create) 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 Clone Methods. */
	//@{
	/**
	* Make a simple clone of this object.
	* @return a pointer to the new object.
	*/
	virtual IBPtr clone() const;

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

	protected:

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

	private:

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

	private:

	/**
	* The \f$\rho\f$ lineshape for the axial-vector terms
	* @param q2 The scale \f$q^2\f$ for the lineshape
	* @param ires Which \f$\rho\f$ multiplet
	*/
	Complex Trho1(Energy2 q2,int ires) const {
	if(ires>=int(_rho1wgts.size())) return 0.;
	double norm = std::accumulate(_rho1wgts.begin(),_rho1wgts.end(),0.);
	unsigned int imin=0,imax=_rho1wgts.size();
	if(ires>0) {
	imin=ires;
	imax=imin+1;
	}
	Complex output(0.);
	for(unsigned int ix=imin;ix<imax;++ix)
	output+=_rho1wgts[ix]*
	Resonance::BreitWignerPWave(q2,_rho1mass[ix],_rho1width[ix],_mpi,_mpi);
	return output/norm;
	}

	/**
	* The \f$K^*\f$ lineshape for the axial-vector terms
	* @param q2 The scale \f$q^2\f$ for the lineshape
	* @param ires Which \f$K^*\f$ multiplet
	*/
	Complex TKstar1(Energy2 q2,int ires) const {
	if(ires>=int(_kstar1wgts.size())) return 0.;
	double norm = std::accumulate(_kstar1wgts.begin(),_kstar1wgts.end(),0.);
	unsigned int imin=0,imax=_kstar1wgts.size();
	if(ires>0) {
	imin=ires;
	imax=imin+1;
	}
	Complex output(0.);
	for(unsigned int ix=imin;ix<imax;++ix)
	output+=_kstar1wgts[ix]*
	Resonance::BreitWignerPWave(q2,_kstar1mass[ix],_kstar1width[ix],_mK,_mpi);
	return output/norm;
	}

	/**
	- * The \f$\rho\f$ lineshape for the vector terms
	+ * The \f$K^*\f$ lineshape for the vector terms
	* @param q2 The scale \f$q^2\f$ for the lineshape
	* @param ires Which \f$K^*\f$ multiplet
	*/
	Complex TKstar2(Energy2 q2,int ires) const {
	if(ires>=int(_kstar2wgts.size())) return 0.;
	double norm = std::accumulate(_kstar2wgts.begin(),_kstar2wgts.end(),0.);
	unsigned int imin=0,imax=_kstar2wgts.size();
	if(ires>0) {
	imin=ires;
	imax=imin+1;
	}
	Complex output(0.);
	for(unsigned int ix=imin;ix<imax;++ix)
	output+=_kstar2wgts[ix]*
	Resonance::BreitWignerPWave(q2,_kstar2mass[ix],_kstar2width[ix],_mK,_mpi);
	return output/norm;
	}

	/**
	* The \f$K_1\f$ line shape
	* @param q2 The scale \f$q^2\f$ for the Breit-Wigner
	* @param iopt Whether this is \f$K^*\pi\f$ or \f$\rho K\f$.
	* @param ires the resonance
	*/
	Complex TK1(Energy2 q2,unsigned int iopt,int ires) const;

	private:

	/**
	* Parameters for the \f$\rho\f$ in the axial-vector terms
	*/
	//@{
	/**
	* Weight for the different resonances
	*/
	vector<double> _rho1wgts;

	/**
	* Masses
	*/
	vector<Energy> _rho1mass;

	/**
	* Widths
	*/
	vector<Energy> _rho1width;
	//@}

	/**
	* Parameters for the \f$K^*\f$ in the axial-vector terms
	*/
	//@{
	/**
	* Weight for the different resonances
	*/
	vector<double> _kstar1wgts;

	/**
	* Masses
	*/
	vector<Energy> _kstar1mass;

	/**
	* Widths
	*/
	vector<Energy> _kstar1width;
	//@}

	/**
	* Parameters for the \f$K^*\f$ in the vector terms
	*/
	//@{
	/**
	* Weight for the different resonances
	*/
	vector<double> _kstar2wgts;

	/**
	* Masses
	*/
	vector<Energy> _kstar2mass;

	/**
	* Widths
	*/
	vector<Energy> _kstar2width;
	//@}

	/**
	* Parameters for the three meson resonances
	*/
	//@{

	/**
	* The masses of the \f$aK1\f$ resonances.
	*/
	vector<Energy> _k1mass;

	/**
	* The widths of the \f$K_1\f$ resonances.
	*/
	vector<Energy> _k1width;

	/**
	* The weights for the different \f$K_1\f$ resonances for \f$K_1\to K^*\pi\f$
	*/
	vector<double> _k1wgta;

	/**
	* The weights for the different \f$K_1\f$ resonaces for \f$K_1\to\rho K\f$.
	*/
	vector<double> _k1wgtb;
	//@}

	/**
	* The pion decay constant, \f$f_\pi\f$.
	*/
	Energy _fpi;

	/**
	* The pion mass
	*/
	Energy _mpi;

	/**
	* The kaon mass
	*/
	Energy _mK;
	//@}
	};

	}

	#endif /* HERWIG_OneKaonTwoPionCurrent_H */
	diff --git a/Decay/WeakCurrents/TwoKaonOnePionCurrent.cc b/Decay/WeakCurrents/TwoKaonOnePionCurrent.cc
	--- a/Decay/WeakCurrents/TwoKaonOnePionCurrent.cc
	+++ b/Decay/WeakCurrents/TwoKaonOnePionCurrent.cc
	@@ -1,1161 +1,988 @@
	// -- C++ --
	//
	// TwoKaonOnePionCurrent.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2017 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 TwoKaonOnePionCurrent class.
	//

	#include "TwoKaonOnePionCurrent.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 "Herwig/PDT/ThreeBodyAllOnCalculator.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "ThePEG/Helicity/epsilon.h"

	using namespace Herwig;

	DescribeClass<TwoKaonOnePionCurrent,WeakCurrent>
	describeHerwigTwoKaonOnePionCurrent("Herwig::TwoKaonOnePionCurrent",
	"HwWeakCurrents.so");
	HERWIG_INTERPOLATOR_CLASSDESC(TwoKaonOnePionCurrent,Energy,Energy2)

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

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

	TwoKaonOnePionCurrent::TwoKaonOnePionCurrent() {
	// the quarks for the different modes
	addDecayMode(2,-1);
	addDecayMode(2,-1);
	addDecayMode(2,-1);
	addDecayMode(2,-1);
	addDecayMode(2,-1);
	addDecayMode(2,-1);
	setInitialModes(7);
	// rho parameters
	// rho parameters for axial-vector pieces
	- _rho1wgts.push_back( 1.0 ); _rho1wgts.push_back(-0.145);
	- _rho1wgts.push_back(0.);
	- _rho1mass.push_back(0.773GeV);_rho1mass.push_back(1.370GeV);
	- _rho1mass.push_back(1.750*GeV);
	- _rho1width.push_back(0.145GeV);_rho1width.push_back(0.510GeV);
	- _rho1width.push_back(0.120*GeV);
	+ _rho1wgts = {1.0,-0.145,0.};
	+ _rho1mass = {0.773GeV,1.370GeV,1.750*GeV};
	+ _rho1width = {0.145GeV,0.510GeV,0.120*GeV};
	// rho parameters for vector pieces
	- _rho2wgts.push_back( 1.0 ); _rho2wgts.push_back(-0.25 );
	- _rho2wgts.push_back(-0.038);
	- _rho2mass.push_back(0.773GeV);_rho2mass.push_back(1.500GeV);
	- _rho2mass.push_back(1.750*GeV);
	- _rho2width.push_back(0.145GeV);_rho2width.push_back(0.220GeV);
	- _rho2width.push_back(0.120*GeV);
	+ _rho2wgts = {1.0,-0.25,-0.038};
	+ _rho2mass = {0.773GeV,1.500GeV,1.750*GeV};
	+ _rho2width = {0.145GeV,0.220GeV,0.120*GeV};
	// K* parameters
	// K* parameters for the axial-vector pieces
	- _kstar1wgts.push_back( 1.0 ); _kstar1wgts.push_back(-0.135);
	- _kstar1wgts.push_back(0.);
	- _kstar1mass.push_back(0.892GeV);_kstar1mass.push_back(1.412GeV);
	- _kstar1mass.push_back(1.714*GeV);
	- _kstar1width.push_back(0.050GeV);_kstar1width.push_back(0.227GeV);
	- _kstar1width.push_back(0.323*GeV);
	- // K* parameters for vector pieces
	- _kstar2wgts.push_back( 1.0 ); _kstar2wgts.push_back(-0.25 );
	- _kstar2wgts.push_back(-0.038);
	- _kstar2mass.push_back(0.892GeV);_kstar2mass.push_back(1.412GeV);
	- _kstar2mass.push_back(1.714*GeV);
	- _kstar2width.push_back(0.050GeV);_kstar2width.push_back(0.227GeV);
	- _kstar2width.push_back(0.323*GeV);
	+ _kstar1wgts = {1.0,-0.135,0.};
	+ _kstar1mass = {0.892GeV,1.412GeV,1.714*GeV};
	+ _kstar1width = {0.050GeV,0.227GeV,0.323*GeV};
	// a_1 parameters
	_initializea1 = false;
	_a1opt = true;
	_a1mass = 1.251*GeV;
	_a1width = 0.475*GeV;
	_a1runwidth = {0GeV, 0GeV, 0GeV, 0GeV, 0*GeV,
	0GeV, 0GeV, 0GeV, 0GeV, 0*GeV,
	0GeV, 0GeV, 1.47729e-06GeV, 1.19209e-05GeV, 3.884e-05*GeV,
	8.83255e-05GeV, 0.00016561GeV, 0.000275439GeV, 0.000422332GeV,
	0.000610773GeV, 0.000845357GeV, 0.00113092GeV, 0.00147264GeV,
	0.00187616GeV, 0.0023477GeV, 0.00289413GeV, 0.00352315GeV,
	0.00424342GeV, 0.0050647GeV, 0.00599808GeV, 0.00705616GeV,
	0.00825335GeV, 0.0096062GeV, 0.0111337GeV, 0.0128579GeV,
	0.0148041GeV, 0.017002GeV, 0.0194858GeV, 0.0222956GeV,
	0.0254781GeV, 0.0290874GeV, 0.0331862GeV, 0.0378467GeV,
	0.0431501GeV, 0.0491862GeV, 0.0560496GeV, 0.0638341GeV,
	0.0726215GeV, 0.0824662GeV, 0.0933765GeV, 0.105297GeV,
	0.118103GeV, 0.131602GeV, 0.145564GeV, 0.159749GeV,
	0.173938GeV, 0.18795GeV, 0.201649GeV, 0.214943GeV,
	0.227773GeV, 0.240109GeV, 0.25194GeV, 0.263268GeV,
	0.274104GeV, 0.284466GeV, 0.294372GeV, 0.303845GeV,
	0.312905GeV, 0.321576GeV, 0.329878GeV, 0.337832GeV,
	0.345456GeV, 0.35277GeV, 0.35979GeV, 0.366532GeV,
	0.373012GeV, 0.379243GeV, 0.38524GeV, 0.391014GeV,
	0.396577GeV, 0.401939GeV, 0.407111GeV, 0.412102GeV,
	0.416923GeV, 0.421577GeV, 0.426078GeV, 0.430427GeV,
	0.434636GeV, 0.43871GeV, 0.442654GeV, 0.446475GeV,
	0.450177GeV, 0.453765GeV, 0.457245GeV, 0.460621GeV,
	0.463899GeV, 0.467077GeV, 0.470164GeV, 0.473162GeV,
	0.476076GeV, 0.478909GeV, 0.481658GeV, 0.484333GeV,
	0.486934GeV, 0.489465GeV, 0.491926GeV, 0.494321GeV,
	0.496651GeV, 0.49892GeV, 0.501128GeV, 0.503277GeV,
	0.505371GeV, 0.507409GeV, 0.509395GeV, 0.511328GeV,
	0.513212GeV, 0.515047GeV, 0.516846GeV, 0.518624GeV,
	0.520285GeV, 0.52194GeV, 0.523553GeV, 0.525124GeV,
	0.526646GeV, 0.52814GeV, 0.529638GeV, 0.531016GeV,
	0.532401GeV, 0.533751GeV, 0.535069GeV, 0.536354GeV,
	0.537608GeV, 0.538831GeV, 0.540039GeV, 0.541194GeV,
	0.542327GeV, 0.543438GeV, 0.544522GeV, 0.545582GeV,
	0.546616GeV, 0.54764GeV, 0.548615GeV, 0.549581GeV,
	0.550525GeV, 0.551449GeV, 0.552351GeV, 0.55324GeV,
	0.554101GeV, 0.554944GeV, 0.555772GeV, 0.556583GeV,
	0.557373GeV, 0.558155GeV, 0.558917GeV, 0.559664GeV,
	0.560396GeV, 0.561114GeV, 0.561849GeV, 0.562508GeV,
	0.563186GeV, 0.563851GeV, 0.564503GeV, 0.565145GeV,
	0.565774GeV, 0.566394GeV, 0.567001GeV, 0.567595GeV,
	0.568182GeV, 0.56876GeV, 0.56933GeV, 0.569886GeV,
	0.570433GeV, 0.570976GeV, 0.571504GeV, 0.572027GeV,
	0.572542GeV, 0.573114GeV, 0.573548GeV, 0.574108GeV,
	0.574524GeV, 0.575002GeV, 0.575473GeV, 0.575937GeV,
	0.576394GeV, 0.576845GeV, 0.57729GeV, 0.57773GeV,
	0.578173GeV, 0.5786GeV, 0.579013GeV, 0.579431GeV,
	0.579834GeV, 0.580246GeV, 0.580649GeV, 0.581045GeV,
	0.581437GeV, 0.581827GeV, 0.582208GeV, 0.582586GeV, 0.582959*GeV};
	_a1runq2 = { 0GeV2 , 0.0158678GeV2, 0.0317356GeV2, 0.0476034GeV2, 0.0634712*GeV2,
	0.079339GeV2, 0.0952068GeV2, 0.111075GeV2, 0.126942GeV2, 0.14281*GeV2,
	0.158678GeV2, 0.174546GeV2, 0.190414GeV2, 0.206281GeV2, 0.222149*GeV2,
	0.238017GeV2, 0.253885GeV2, 0.269753GeV2, 0.285621GeV2, 0.301488*GeV2,
	0.317356GeV2, 0.333224GeV2, 0.349092GeV2, 0.36496GeV2, 0.380827*GeV2,
	0.396695GeV2, 0.412563GeV2, 0.428431GeV2, 0.444299GeV2, 0.460166*GeV2,
	0.476034GeV2, 0.491902GeV2, 0.50777GeV2, 0.523638GeV2, 0.539505*GeV2,
	0.555373GeV2, 0.571241GeV2, 0.587109GeV2, 0.602977GeV2, 0.618844*GeV2,
	0.634712GeV2, 0.65058GeV2, 0.666448GeV2, 0.682316GeV2, 0.698183*GeV2,
	0.714051GeV2, 0.729919GeV2, 0.745787GeV2, 0.761655GeV2, 0.777523*GeV2,
	0.79339GeV2, 0.809258GeV2, 0.825126GeV2, 0.840994GeV2, 0.856862*GeV2,
	0.872729GeV2, 0.888597GeV2, 0.904465GeV2, 0.920333GeV2, 0.936201*GeV2,
	0.952068GeV2, 0.967936GeV2, 0.983804GeV2, 0.999672GeV2, 1.01554*GeV2,
	1.03141GeV2, 1.04728GeV2, 1.06314GeV2, 1.07901GeV2, 1.09488*GeV2,
	1.11075GeV2, 1.12661GeV2, 1.14248GeV2, 1.15835GeV2, 1.17422*GeV2,
	1.19009GeV2, 1.20595GeV2, 1.22182GeV2, 1.23769GeV2, 1.25356*GeV2,
	1.26942GeV2, 1.28529GeV2, 1.30116GeV2, 1.31703GeV2, 1.3329*GeV2,
	1.34876GeV2, 1.36463GeV2, 1.3805GeV2, 1.39637GeV2, 1.41223*GeV2,
	1.4281GeV2, 1.44397GeV2, 1.45984GeV2, 1.47571GeV2, 1.49157*GeV2,
	1.50744GeV2, 1.52331GeV2, 1.53918GeV2, 1.55505GeV2, 1.57091*GeV2,
	1.58678GeV2, 1.60265GeV2, 1.61852GeV2, 1.63438GeV2, 1.65025*GeV2,
	1.66612GeV2, 1.68199GeV2, 1.69786GeV2, 1.71372GeV2, 1.72959*GeV2,
	1.74546GeV2, 1.76133GeV2, 1.77719GeV2, 1.79306GeV2, 1.80893*GeV2,
	1.8248GeV2, 1.84067GeV2, 1.85653GeV2, 1.8724GeV2, 1.88827*GeV2,
	1.90414GeV2, 1.92GeV2, 1.93587GeV2, 1.95174GeV2, 1.96761*GeV2,
	1.98348GeV2, 1.99934GeV2, 2.01521GeV2, 2.03108GeV2, 2.04695*GeV2,
	2.06281GeV2, 2.07868GeV2, 2.09455GeV2, 2.11042GeV2, 2.12629*GeV2,
	2.14215GeV2, 2.15802GeV2, 2.17389GeV2, 2.18976GeV2, 2.20563*GeV2,
	2.22149GeV2, 2.23736GeV2, 2.25323GeV2, 2.2691GeV2, 2.28496*GeV2,
	2.30083GeV2, 2.3167GeV2, 2.33257GeV2, 2.34844GeV2, 2.3643*GeV2,
	2.38017GeV2, 2.39604GeV2, 2.41191GeV2, 2.42777GeV2, 2.44364*GeV2,
	2.45951GeV2, 2.47538GeV2, 2.49125GeV2, 2.50711GeV2, 2.52298*GeV2,
	2.53885GeV2, 2.55472GeV2, 2.57058GeV2, 2.58645GeV2, 2.60232*GeV2,
	2.61819GeV2, 2.63406GeV2, 2.64992GeV2, 2.66579GeV2, 2.68166*GeV2,
	2.69753GeV2, 2.71339GeV2, 2.72926GeV2, 2.74513GeV2, 2.761*GeV2,
	2.77687GeV2, 2.79273GeV2, 2.8086GeV2, 2.82447GeV2, 2.84034*GeV2,
	2.85621GeV2, 2.87207GeV2, 2.88794GeV2, 2.90381GeV2, 2.91968*GeV2,
	2.93554GeV2, 2.95141GeV2, 2.96728GeV2, 2.98315GeV2, 2.99902*GeV2,
	3.01488GeV2, 3.03075GeV2, 3.04662GeV2, 3.06249GeV2, 3.07835*GeV2,
	3.09422GeV2, 3.11009GeV2, 3.12596GeV2, 3.14183GeV2, 3.15769*GeV2};
	// parameters for the T_omega function
	- _epsomega=0.05;
	+ _epsomega = 0.05;
	_omegamass = 0.782*GeV;
	_omegawidth = 0.00843*GeV;
	_phimass = 1.020*GeV;
	_phiwidth = 0.00443*GeV;
	_omegaKstarwgt=1./sqrt(2.);
	// the pion decay constant
	- _fpi=130.7*MeV/sqrt(2.);
	- _mpi=ZERO;
	- _mK=ZERO;
	- _maxmass=ZERO;
	- _maxcalc=ZERO;
	+ _fpi = 130.7*MeV/sqrt(2.);
	+ _mpi = ZERO;
	+ _mK = ZERO;
	+ _maxmass = ZERO;
	+ _maxcalc = ZERO;
	}


	void TwoKaonOnePionCurrent::persistentOutput(PersistentOStream & os) const {
	os << _a1runinter
	<< _rho1wgts << ounit(_rho1mass,GeV) << ounit(_rho1width,GeV)
	<< _rho2wgts << ounit(_rho2mass,GeV) << ounit(_rho2width,GeV)
	- << _kstar1wgts << ounit(_kstar1mass,GeV) << ounit(_kstar1width,GeV)
	- << _kstar2wgts << ounit(_kstar2mass,GeV) << ounit(_kstar2width,GeV)
	+ << _kstar1wgts << ounit(_kstar1mass,GeV) << ounit(_kstar1width,GeV)
	<< ounit(_a1mass,GeV) << ounit(_a1width,GeV)
	<< ounit(_a1runwidth,GeV) << ounit(_a1runq2,GeV2) << _epsomega
	<< ounit(_omegamass,GeV) << ounit(_omegawidth,GeV)
	<< ounit(_phimass,GeV) << ounit(_phiwidth,GeV) << _omegaKstarwgt
	<< ounit(_fpi,GeV) << ounit(_mpi,GeV) << ounit(_mK,GeV)
	<< _initializea1 << _a1opt
	<< ounit(_maxmass,GeV) << ounit(_maxcalc,GeV);
	}

	void TwoKaonOnePionCurrent::persistentInput(PersistentIStream & is, int) {
	is >> _a1runinter
	>> _rho1wgts >> iunit(_rho1mass,GeV) >> iunit(_rho1width,GeV)
	>> _rho2wgts >> iunit(_rho2mass,GeV) >> iunit(_rho2width,GeV)
	- >> _kstar1wgts >> iunit(_kstar1mass,GeV) >> iunit(_kstar1width,GeV)
	- >> _kstar2wgts >> iunit(_kstar2mass,GeV) >> iunit(_kstar2width,GeV)
	+ >> _kstar1wgts >> iunit(_kstar1mass,GeV) >> iunit(_kstar1width,GeV)
	>> iunit(_a1mass,GeV) >> iunit(_a1width,GeV)
	>> iunit(_a1runwidth,GeV) >> iunit(_a1runq2,GeV2) >> _epsomega
	>> iunit(_omegamass,GeV) >> iunit(_omegawidth,GeV)
	>> iunit(_phimass,GeV) >> iunit(_phiwidth,GeV) >> _omegaKstarwgt
	>> iunit(_fpi,GeV) >> iunit(_mpi,GeV) >> iunit(_mK,GeV)
	>> _initializea1 >> _a1opt
	>> iunit(_maxmass,GeV) >> iunit(_maxcalc,GeV);
	}


	void TwoKaonOnePionCurrent::Init() {

	static ClassDocumentation<TwoKaonOnePionCurrent> documentation
	("The TwoKaonOnePionCurrent class implements the model of "
	"Z. Phys. C 69 (1996) 243 [arXiv:hep-ph/9503474]"
	" for the weak current with three "
	"mesons, at least one of which is a kaon",
	"The TwoKaonOnePionCurrent class implements the model of "
	"\\cite{Finkemeier:1995sr} for the weak current with three "
	"mesons, at least one of which is a kaon.",
	"\\bibitem{Finkemeier:1995sr}\n"
	"M.~Finkemeier and E.~Mirkes,\n"
	"Z.\\ Phys.\\ C {\\bf 69} (1996) 243 [arXiv:hep-ph/9503474].\n"
	" %%CITATION = ZEPYA,C69,243;%%\n"

	);

	static Switch<TwoKaonOnePionCurrent,bool> interfaceInitializea1
	("Initializea1",
	"Initialise the calculation of the a_1 running width",
	&TwoKaonOnePionCurrent::_initializea1, false, false, false);
	static SwitchOption interfaceInitializea1Initialization
	(interfaceInitializea1,
	"Yes",
	"Initialize the calculation",
	true);
	static SwitchOption interfaceInitializea1NoInitialization
	(interfaceInitializea1,
	"No",
	"Use the default values",
	false);

	static Parameter<TwoKaonOnePionCurrent,Energy> interfaceA1Width
	("A1Width",
	"The a_1 width if using local values.",
	&TwoKaonOnePionCurrent::_a1width, GeV, 0.599GeV, ZERO, 10.0GeV,
	false, false, false);

	static Parameter<TwoKaonOnePionCurrent,Energy> interfaceA1Mass
	("A1Mass",
	"The a_1 mass if using local values.",
	&TwoKaonOnePionCurrent::_a1mass, GeV, 1.251GeV, ZERO, 10.0GeV,
	false, false, false);

	static Parameter<TwoKaonOnePionCurrent,Energy> interfaceFPi
	("FPi",
	"The pion decay constant",
	&TwoKaonOnePionCurrent::_fpi, MeV, 92.4MeV, ZERO, 200.0MeV,
	false, false, true);

	static ParVector<TwoKaonOnePionCurrent,Energy> interfaceRhoAxialMasses
	("RhoAxialMasses",
	"The masses for the rho resonances if used local values",
	&TwoKaonOnePionCurrent::_rho1mass, GeV, -1, 1.0GeV, ZERO, 10.0GeV,
	false, false, true);

	static ParVector<TwoKaonOnePionCurrent,Energy> interfaceRhoAxialWidths
	("RhoAxialWidths",
	"The widths for the rho resonances if used local values",
	&TwoKaonOnePionCurrent::_rho1width, GeV, -1, 1.0GeV, ZERO, 10.0GeV,
	false, false, true);

	static ParVector<TwoKaonOnePionCurrent,Energy> interfaceRhoVectorMasses
	("RhoVectorMasses",
	"The masses for the rho resonances if used local values",
	&TwoKaonOnePionCurrent::_rho2mass, GeV, -1, 1.0GeV, ZERO, 10.0GeV,
	false, false, true);

	static ParVector<TwoKaonOnePionCurrent,Energy> interfaceRhoVectorWidths
	("RhoVectorWidths",
	"The widths for the rho resonances if used local values",
	&TwoKaonOnePionCurrent::_rho2width, GeV, -1, 1.0GeV, ZERO, 10.0GeV,
	false, false, true);

	static ParVector<TwoKaonOnePionCurrent,Energy> interfaceKstarAxialMasses
	("KstarAxialMasses",
	"The masses for the Kstar resonances if used local values",
	&TwoKaonOnePionCurrent::_kstar1mass, GeV, -1, 1.0GeV, ZERO, 10.0GeV,
	false, false, true);

	static ParVector<TwoKaonOnePionCurrent,Energy> interfaceKstarAxialWidths
	("KstarAxialWidths",
	"The widths for the Kstar resonances if used local values",
	&TwoKaonOnePionCurrent::_kstar1width, GeV, -1, 1.0GeV, ZERO, 10.0GeV,
	false, false, true);

	- static ParVector<TwoKaonOnePionCurrent,Energy> interfaceKstarVectorMasses
	- ("KstarVectorMasses",
	- "The masses for the Kstar resonances if used local values",
	- &TwoKaonOnePionCurrent::_kstar2mass, GeV, -1, 1.0GeV, ZERO, 10.0GeV,
	- false, false, true);
	-
	- static ParVector<TwoKaonOnePionCurrent,Energy> interfaceKstarVectorWidths
	- ("KstarVectorWidths",
	- "The widths for the Kstar resonances if used local values",
	- &TwoKaonOnePionCurrent::_kstar2width, GeV, -1, 1.0GeV, ZERO, 10.0GeV,
	- false, false, true);
	-
	static ParVector<TwoKaonOnePionCurrent,double> interfaceAxialRhoWeight
	("AxialRhoWeight",
	"The weights of the different rho resonances in the F1,2,3 form factor",
	&TwoKaonOnePionCurrent::_rho1wgts,
	0, 0, 0, -1000, 1000, false, false, true);

	static ParVector<TwoKaonOnePionCurrent,double> interfaceAxialKStarWeight
	("AxialKStarWeight",
	"The weights of the different Kstar resonances in the F1,2,3 form factor",
	&TwoKaonOnePionCurrent::_kstar1wgts,
	0, 0, 0, -1000, 1000, false, false, true);

	static ParVector<TwoKaonOnePionCurrent,double> interfaceVectorRhoWeight
	("VectorRhoWeight",
	"The weights of the different rho resonances in the F1,2,3 form factor",
	&TwoKaonOnePionCurrent::_rho2wgts,
	0, 0, 0, -1000, 1000, false, false, true);

	- static ParVector<TwoKaonOnePionCurrent,double> interfaceVectorKStarWeight
	- ("VectorKStarWeight",
	- "The weights of the different Kstar resonances in the F1,2,3 form factor",
	- &TwoKaonOnePionCurrent::_kstar2wgts,
	- 0, 0, 0, -1000, 1000, false, false, true);
	-
	static Switch<TwoKaonOnePionCurrent,bool> interfacea1WidthOption
	("a1WidthOption",
	"Option for the treatment of the a1 width",
	&TwoKaonOnePionCurrent::_a1opt, true, false, false);
	static SwitchOption interfacea1WidthOptionLocal
	(interfacea1WidthOption,
	"Local",
	"Use a calculation of the running width based on the parameters as"
	" interpolation table.",
	true);
	static SwitchOption interfacea1WidthOptionParam
	(interfacea1WidthOption,
	"Kuhn",
	"Use the parameterization of Kuhn and Santamaria for default parameters."
	" This should only be used for testing vs TAUOLA",
	false);

	static ParVector<TwoKaonOnePionCurrent,Energy> interfacea1RunningWidth
	("a1RunningWidth",
	"The values of the a_1 width for interpolation to giving the running width.",
	&TwoKaonOnePionCurrent::_a1runwidth, GeV, -1, 1.0GeV, ZERO, 10.0GeV,
	false, false, true);

	static ParVector<TwoKaonOnePionCurrent,Energy2> interfacea1RunningQ2
	("a1RunningQ2",
	"The values of the q^2 for interpolation to giving the running width.",
	&TwoKaonOnePionCurrent::_a1runq2, GeV2, -1, 1.0GeV2, ZERO, 10.0GeV2,
	false, false, true);

	static Parameter<TwoKaonOnePionCurrent,double> interfaceEpsOmega
	("EpsOmega",
	"The omega-phi mixing ",
	&TwoKaonOnePionCurrent::_epsomega, 0.05, 0.0, 1.0,
	false, false, Interface::limited);

	static Parameter<TwoKaonOnePionCurrent,Energy> interfaceOmegaMass
	("OmegaMass",
	"The mass of the omega meson",
	&TwoKaonOnePionCurrent::_omegamass, GeV, 0.782GeV, ZERO, 10.0GeV,
	false, false, Interface::limited);

	static Parameter<TwoKaonOnePionCurrent,Energy> interfaceOmegaWidth
	("OmegaWidth",
	"The width of the omega meson",
	&TwoKaonOnePionCurrent::_omegawidth, GeV, 0.00843GeV, ZERO, 10.0GeV,
	false, false, Interface::limited);

	static Parameter<TwoKaonOnePionCurrent,Energy> interfacePhiMass
	("PhiMass",
	"The mass of the phi meson",
	&TwoKaonOnePionCurrent::_phimass, GeV, 1.020GeV, ZERO, 10.0GeV,
	false, false, Interface::limited);

	static Parameter<TwoKaonOnePionCurrent,Energy> interfacePhiWidth
	("PhiWidth",
	"The width of the phi meson",
	&TwoKaonOnePionCurrent::_phiwidth, GeV, 0.00443GeV, ZERO, 10.0GeV,
	false, false, Interface::limited);

	static Parameter<TwoKaonOnePionCurrent,double> interfaceOmegaKStarWeight
	("OmegaKStarWeight",
	"The relative weight of the omega-phi and K* terms",
	&TwoKaonOnePionCurrent::_omegaKstarwgt, 1./sqrt(2.), 0.0, 100.0,
	false, false, Interface::limited);

	}

	void TwoKaonOnePionCurrent::inita1Width(int iopt) {
	if(iopt==-1) {
	_maxcalc=_maxmass;
	if(!_initializea1\|\|_maxmass==ZERO) return;
	// parameters for the table of values
	Energy2 step(sqr(_maxmass)/199.);
	// integrator to perform the integral
	vector<double> inweights;inweights.push_back(0.5);inweights.push_back(0.5);
	vector<int> intype;intype.push_back(2);intype.push_back(3);
	Energy mrho(getParticleData(ParticleID::rhoplus)->mass()),
	wrho(getParticleData(ParticleID::rhoplus)->width());
	vector<Energy> inmass(2,mrho),inwidth(2,wrho);
	vector<double> inpow(2,0.0);
	ThreeBodyAllOnCalculator<TwoKaonOnePionCurrent>
	widthgen(inweights,intype,inmass,inwidth,inpow,*this,0,_mpi,_mpi,_mpi);
	// normalisation constant to give physical width if on shell
	double a1const(_a1width/(widthgen.partialWidth(sqr(_a1mass))));
	// loop to give the values
	_a1runq2.clear();_a1runwidth.clear();
	for(Energy2 moff2 = ZERO; moff2<=sqr(_maxmass); moff2+=step) {
	_a1runwidth.push_back(widthgen.partialWidth(moff2)*a1const);
	_a1runq2.push_back(moff2);
	}
	}
	// set up the interpolator
	else if(iopt==0) {
	_a1runinter = make_InterpolatorPtr(_a1runwidth,_a1runq2,3);
	}
	}

	// complete the construction of the decay mode for integration
	bool TwoKaonOnePionCurrent::createMode(int icharge, tcPDPtr resonance,
	IsoSpin::IsoSpin Itotal, IsoSpin::I3 i3,
	unsigned int imode,PhaseSpaceModePtr mode,
	unsigned int iloc,int ires,
	PhaseSpaceChannel phase, Energy upp ) {
	+ // check the charge
	if(abs(icharge)!=3) return false;
	+ // check the total isospin
	+ if(Itotal!=IsoSpin::IUnknown) {
	+ if(Itotal!=IsoSpin::IOne) return false;
	+ }
	+ // check I_3
	+ if(i3!=IsoSpin::I3Unknown) {
	+ switch(i3) {
	+ case IsoSpin::I3Zero:
	+ if(imode<=1) return false;
	+ break;
	+ case IsoSpin::I3One:
	+ if( imode>1 \|\| icharge ==-3) return false;
	+ break;
	+ case IsoSpin::I3MinusOne:
	+ if( imode>1 \|\| icharge == 3) return false;
	+ default:
	+ return false;
	+ }
	+ }
	+ // get the particles and check the mass
	int iq(0),ia(0);
	tPDVector extpart(particles(1,imode,iq,ia));
	Energy min(ZERO);
	for(unsigned int ix=0;ix<extpart.size();++ix) min+=extpart[ix]->massMin();
	if(min>upp) return false;
	// the particles we will use a lot
	tPDPtr a1 = getParticleData(ParticleID::a_1minus);
	- tPDPtr k1[2] = {getParticleData(ParticleID::K_1minus),
	- getParticleData(ParticleID::Kstar_1minus)};
	_maxmass=max(_maxmass,upp);
	// the rho0 resonances
	tPDPtr rho0[3] ={getParticleData( 113),getParticleData( 100113),
	getParticleData( 30113)};
	// the charged rho resonances
	tPDPtr rhoc[3] ={getParticleData(-213),getParticleData(-100213),
	getParticleData(-30213)};
	// the K*0 resonances
	tPDPtr Kstar0[3]={getParticleData( 313),getParticleData( 100313),
	getParticleData( 30313)};
	// the charged K* resonances
	tPDPtr Kstarc[3]={getParticleData(-323),getParticleData(-100323),
	getParticleData(-30323)};
	if(icharge==3) {
	a1 = a1->CC();
	- k1[0] = k1[0]->CC();
	- k1[1] = k1[1]->CC();
	for(unsigned int ix=0;ix<3;++ix) {
	if(rhoc[ix]) rhoc[ix]=rhoc[ix]->CC();
	if(Kstar0[ix]) Kstar0[ix]=Kstar0[ix]->CC();
	if(Kstarc[ix]) Kstarc[ix]=Kstarc[ix]->CC();
	}
	}
	if(imode==0) {
	// channels for K- pi- K+
	for(unsigned int ix=0;ix<3;++ix) {
	- if(Kstar0[ix]) {
	+ if(!resonance \|\| resonance==a1) {
	mode->addChannel((PhaseSpaceChannel(phase),ires,a1,ires+1,Kstar0[ix],ires+1,iloc+1,
	ires+2,iloc+2,ires+2,iloc+3));
	- }
	- if(rho0[ix]) {
	mode->addChannel((PhaseSpaceChannel(phase),ires,a1,ires+1,rho0[ix],ires+1,iloc+2,
	ires+2,iloc+1,ires+2,iloc+3));
	}
	for(unsigned int iy=0;iy<3;++iy) {
	- if(!rhoc[ix]) continue;
	- if(Kstar0[iy]) {
	- mode->addChannel((PhaseSpaceChannel(phase),ires,rhoc[ix],ires+1,Kstar0[iy],ires+1,iloc+1,
	- ires+2,iloc+2,ires+2,iloc+3));
	- }
	+ if(resonance && resonance !=rhoc[ix]) continue;
	+ mode->addChannel((PhaseSpaceChannel(phase),ires,rhoc[ix],ires+1,Kstar0[iy],ires+1,iloc+1,
	+ ires+2,iloc+2,ires+2,iloc+3));
	}
	}
	}
	else if(imode==1) {
	// channels for K0 pi- K0bar
	for(unsigned int ix=0;ix<3;++ix) {
	- if(Kstarc[ix]) {
	+ if(!resonance \|\| resonance==a1) {
	mode->addChannel((PhaseSpaceChannel(phase),ires,a1,ires+1,Kstarc[ix],ires+1,iloc+1,
	ires+2,iloc+2,ires+2,iloc+3));
	- }
	- if(rho0[ix]) {
	mode->addChannel((PhaseSpaceChannel(phase),ires,a1,ires+1,rho0[ix],ires+1,iloc+2,
	ires+2,iloc+1,ires+2,iloc+3));
	}
	for(unsigned int iy=0;iy<3;++iy) {
	- if(!rhoc[ix]) continue;
	- if(Kstarc[iy]) {
	- mode->addChannel((PhaseSpaceChannel(phase),ires,rhoc[ix],ires+1,Kstarc[iy],ires+1,iloc+1,
	- ires+2,iloc+2,ires+2,iloc+3));
	- }
	+ if(resonance && resonance !=rhoc[ix]) continue;
	+ mode->addChannel((PhaseSpaceChannel(phase),ires,rhoc[ix],ires+1,Kstarc[iy],ires+1,iloc+1,
	+ ires+2,iloc+2,ires+2,iloc+3));
	}
	}
	}
	else if(imode==2) {
	// channels for K- pi0 K0
	for(unsigned int ix=0;ix<3;++ix) {
	- if(Kstar0[ix]) {
	+ if(!resonance \|\| resonance==a1) {
	mode->addChannel((PhaseSpaceChannel(phase),ires,a1,ires+1,Kstar0[ix],ires+1,iloc+1,
	ires+2,iloc+2,ires+2,iloc+3));
	- }
	- if(Kstarc[ix]) {
	mode->addChannel((PhaseSpaceChannel(phase),ires,a1,ires+1,Kstarc[ix],ires+1,iloc+3,
	ires+2,iloc+1,ires+2,iloc+2));
	- }
	- if(rhoc[ix]) {
	mode->addChannel((PhaseSpaceChannel(phase),ires,a1,ires+1,rhoc[ix],ires+1,iloc+2,
	ires+2,iloc+1,ires+2,iloc+3));
	}
	for(unsigned int iy=0;iy<3;++iy) {
	- if(!rhoc[ix]) continue;
	- if(Kstar0[iy]) {
	- mode->addChannel((PhaseSpaceChannel(phase),ires,rhoc[ix],ires+1,Kstar0[iy],ires+1,iloc+1,
	- ires+2,iloc+2,ires+2,iloc+3));
	- }
	- if(Kstarc[iy]) {
	- mode->addChannel((PhaseSpaceChannel(phase),ires,rhoc[ix],ires+1,Kstarc[iy],ires+1,iloc+3,
	- ires+2,iloc+1,ires+2,iloc+2));
	- }
	+ if(resonance && resonance !=rhoc[ix]) continue;
	+ mode->addChannel((PhaseSpaceChannel(phase),ires,rhoc[ix],ires+1,Kstar0[iy],ires+1,iloc+1,
	+ ires+2,iloc+2,ires+2,iloc+3));
	+ mode->addChannel((PhaseSpaceChannel(phase),ires,rhoc[ix],ires+1,Kstarc[iy],ires+1,iloc+3,
	+ ires+2,iloc+1,ires+2,iloc+2));
	}
	}
	}
	else if(imode==3\|\|imode==4) {
	// channels for K_S0 pi- K_S0 and K_L0 pi- K_L0
	for(unsigned int ix=0;ix<3;++ix) {
	- if(Kstarc[ix]) {
	+ if(!resonance \|\| resonance==a1) {
	mode->addChannel((PhaseSpaceChannel(phase),ires,a1,ires+1,Kstarc[ix],ires+1,iloc+1,
	ires+2,iloc+2,ires+2,iloc+3));
	mode->addChannel((PhaseSpaceChannel(phase),ires,a1,ires+1,Kstarc[ix],ires+1,iloc+3,
	ires+2,iloc+1,ires+2,iloc+2));
	}
	for(unsigned int iy=0;iy<3;++iy) {
	- if(!rhoc[ix]) continue;
	- if(Kstarc[iy]) {
	- mode->addChannel((PhaseSpaceChannel(phase),ires,rhoc[ix],ires+1,Kstarc[iy],ires+1,iloc+1,
	- ires+2,iloc+2,ires+2,iloc+3));
	- mode->addChannel((PhaseSpaceChannel(phase),ires,rhoc[ix],ires+1,Kstarc[iy],ires+1,iloc+3,
	- ires+2,iloc+1,ires+2,iloc+2));
	- }
	+ if(resonance && resonance !=rhoc[ix]) continue;
	+ mode->addChannel((PhaseSpaceChannel(phase),ires,rhoc[ix],ires+1,Kstarc[iy],ires+1,iloc+1,
	+ ires+2,iloc+2,ires+2,iloc+3));
	+ mode->addChannel((PhaseSpaceChannel(phase),ires,rhoc[ix],ires+1,Kstarc[iy],ires+1,iloc+3,
	+ ires+2,iloc+1,ires+2,iloc+2));
	}
	}
	}
	else if(imode==5) {
	// channels for K_S0 pi- K_L0
	for(unsigned int ix=0;ix<3;++ix) {
	- if(Kstarc[ix]) {
	+ if(!resonance \|\| resonance==a1) {
	mode->addChannel((PhaseSpaceChannel(phase),ires,a1,ires+1,Kstarc[ix],ires+1,iloc+1,
	ires+2,iloc+2,ires+2,iloc+3));
	mode->addChannel((PhaseSpaceChannel(phase),ires,a1,ires+1,Kstarc[ix],ires+1,iloc+3,
	ires+2,iloc+1,ires+2,iloc+2));
	- }
	- if(rho0[ix])
	mode->addChannel((PhaseSpaceChannel(phase),ires,a1,ires+1,rho0[ix],ires+1,iloc+2,
	ires+2,iloc+1,ires+2,iloc+3));
	+ }
	for(unsigned int iy=0;iy<3;++iy) {
	- if(!rhoc[ix]) continue;
	- if(Kstarc[iy]) {
	- mode->addChannel((PhaseSpaceChannel(phase),ires,rhoc[ix],ires+1,Kstarc[ix],ires+1,iloc+1,
	- ires+2,iloc+2,ires+2,iloc+3));
	- mode->addChannel((PhaseSpaceChannel(phase),ires,rhoc[ix],ires+1,Kstarc[ix],ires+1,iloc+3,
	- ires+2,iloc+1,ires+2,iloc+2));
	- }
	+ if(resonance && resonance !=rhoc[ix]) continue;
	+ mode->addChannel((PhaseSpaceChannel(phase),ires,rhoc[ix],ires+1,Kstarc[ix],ires+1,iloc+1,
	+ ires+2,iloc+2,ires+2,iloc+3));
	+ mode->addChannel((PhaseSpaceChannel(phase),ires,rhoc[ix],ires+1,Kstarc[ix],ires+1,iloc+3,
	+ ires+2,iloc+1,ires+2,iloc+2));
	}
	}
	}
	+ // update the integration parameters
	for(unsigned int ix=0;ix<_rho1mass.size();++ix) {
	mode->resetIntermediate(rhoc[ix],_rho1mass[ix],
	_rho1width[ix]);
	mode->resetIntermediate(rho0[ix],_rho1mass[ix],
	_rho1width[ix]);
	}
	- // K star parameters in the base class
	for(unsigned int ix=0;ix<_kstar1mass.size();++ix) {
	mode->resetIntermediate(Kstarc[ix],_kstar1mass[ix],
	_kstar1width[ix]);
	mode->resetIntermediate(Kstar0[ix],_kstar1mass[ix],
	_kstar1width[ix]);
	}
	return true;
	}

	void TwoKaonOnePionCurrent::dataBaseOutput(ofstream & os,
	bool header,bool create) const {
	if(header) os << "update decayers set parameters=\"";
	if(create) os << "create Herwig::TwoKaonOnePionCurrent "
	<< name() << " HwWeakCurrents.so\n";
	for(unsigned int ix=0;ix<_rho1wgts.size();++ix) {
	if(ix<3) {
	os << "newdef " << name() << ":AxialRhoWeight " << ix
	<< " " << _rho1wgts[ix] << "\n";
	}
	else {
	os << "insert " << name() << ":AxialRhoWeight " << ix
	<< " " << _rho1wgts[ix] << "\n";
	}
	}
	for(unsigned int ix=0;ix<_kstar1wgts.size();++ix) {
	if(ix<3) {
	os << "newdef " << name() << ":AxialKStarWeight " << ix
	<< " " << _kstar1wgts[ix] << "\n";}
	else {
	os << "insert " << name() << ":AxialKStarWeight " << ix
	<< " " << _kstar1wgts[ix] << "\n";
	}
	}
	for(unsigned int ix=0;ix<_rho2wgts.size();++ix) {
	if(ix<3) {
	os << "newdef " << name() << ":VectorRhoWeight " << ix
	<< " " << _rho2wgts[ix] << "\n";
	}
	else {
	os << "insert " << name() << ":VectorRhoWeight " << ix
	<< " " << _rho2wgts[ix] << "\n";
	}
	}
	- for(unsigned int ix=0;ix<_kstar2wgts.size();++ix) {
	- if(ix<3) {
	- os << "newdef " << name() << ":VectorKStarWeight " << ix
	- << " " << _kstar2wgts[ix] << "\n";}
	- else {
	- os << "insert " << name() << ":VectorKStarWeight " << ix
	- << " " << _kstar2wgts[ix] << "\n";
	- }
	- }
	os << "newdef " << name() << ":OmegaKStarWeight " << _omegaKstarwgt << "\n";
	os << "newdef " << name() << ":EpsOmega " << _epsomega << "\n";
	os << "newdef " << name() << ":Initializea1 " << _initializea1 << "\n";
	os << "newdef " << name() << ":a1WidthOption " << _a1opt << "\n";
	for(unsigned int ix=0;ix<_a1runwidth.size();++ix) {
	os << "newdef " << name() << ":a1RunningWidth " << ix
	<< " " << _a1runwidth[ix]/GeV << "\n";
	}
	for(unsigned int ix=0;ix<_a1runq2.size();++ix) {
	os << "newdef " << name() << ":a1RunningQ2 " << ix
	<< " " << _a1runq2[ix]/GeV2 << "\n";
	}
	os << "newdef " << name() << ":A1Width " << _a1width/GeV << "\n";
	os << "newdef " << name() << ":A1Mass " << _a1mass/GeV << "\n";
	os << "newdef " << name() << ":OmegaWidth " << _omegawidth/GeV << "\n";
	os << "newdef " << name() << ":OmegaMass " << _omegamass/GeV << "\n";
	os << "newdef " << name() << ":PhiWidth " << _phiwidth/GeV << "\n";
	os << "newdef " << name() << ":PhiMass " << _phimass/GeV << "\n";
	os << "newdef " << name() << ":FPi " << _fpi/MeV << "\n";
	for(unsigned int ix=0;ix<_rho1mass.size();++ix) {
	if(ix<3) os << "newdef " << name() << ":RhoAxialMasses " << ix
	<< " " << _rho1mass[ix]/GeV << "\n";
	else os << "insert " << name() << ": RhoAxialMasses" << ix
	<< " " << _rho1mass[ix]/GeV << "\n";
	}
	for(unsigned int ix=0;ix<_rho1width.size();++ix) {
	if(ix<3) os << "newdef " << name() << ":RhoAxialWidths " << ix
	<< " " << _rho1width[ix]/GeV << "\n";
	else os << "insert " << name() << ":RhoAxialWidths " << ix
	<< " " << _rho1width[ix]/GeV << "\n";
	}
	for(unsigned int ix=0;ix<_rho2mass.size();++ix) {
	if(ix<3) os << "newdef " << name() << ":RhoVectorMasses " << ix
	<< " " << _rho2mass[ix]/GeV << "\n";
	else os << "insert " << name() << ": RhoVectorMasses" << ix
	<< " " << _rho2mass[ix]/GeV << "\n";
	}
	for(unsigned int ix=0;ix<_rho2width.size();++ix) {
	if(ix<3) os << "newdef " << name() << ":RhoVectorWidths " << ix
	<< " " << _rho2width[ix]/GeV << "\n";
	else os << "insert " << name() << ":RhoVectorWidths " << ix
	<< " " << _rho2width[ix]/GeV << "\n";
	}
	for(unsigned int ix=0;ix<_kstar1mass.size();++ix) {
	if(ix<3) os << "newdef " << name() << ":KstarAxialMasses " << ix
	<< " " << _kstar1mass[ix]/GeV << "\n";
	else os << "insert " << name() << ": KstarAxialMasses" << ix
	<< " " << _kstar1mass[ix]/GeV << "\n";
	}
	for(unsigned int ix=0;ix<_kstar1width.size();++ix) {
	if(ix<3) os << "newdef " << name() << ":KstarAxialWidths " << ix
	<< " " << _kstar1width[ix]/GeV << "\n";
	else os << "insert " << name() << ":KstarAxialWidths " << ix
	<< " " << _kstar1width[ix]/GeV << "\n";
	}
	- for(unsigned int ix=0;ix<_kstar2mass.size();++ix) {
	- if(ix<3) os << "newdef " << name() << ":KstarVectorMasses " << ix
	- << " " << _kstar2mass[ix]/GeV << "\n";
	- else os << "insert " << name() << ": KstarVectorMasses" << ix
	- << " " << _kstar2mass[ix]/GeV << "\n";
	- }
	- for(unsigned int ix=0;ix<_kstar2width.size();++ix) {
	- if(ix<3) os << "newdef " << name() << ":KstarVectorWidths " << ix
	- << " " << _kstar2width[ix]/GeV << "\n";
	- else os << "insert " << name() << ":KstarVectorWidths " << ix
	- << " " << _kstar2width[ix]/GeV << "\n";
	- }
	WeakCurrent::dataBaseOutput(os,false,false);
	if(header) os << "\n\" where BINARY ThePEGName=\""
	<< fullName() << "\";" << endl;
	}

	void TwoKaonOnePionCurrent::doinit() {
	WeakCurrent::doinit();
	// masses for the running widths
	_mpi = getParticleData(ParticleID::piplus)->mass();
	_mK = getParticleData(ParticleID::K0) ->mass();
	// initialise the a_1 running width calculation
	inita1Width(-1);
	inita1Width(0);
	}

	-TwoKaonOnePionCurrent::FormFactors
	-TwoKaonOnePionCurrent::calculateFormFactors(const int ichan,const int imode,
	- Energy2 q2,Energy2 s1,
	- Energy2 s2,Energy2 s3) const {
	- useMe();
	- Complex F1, F2, F5;
	- // calculate the K- pi - K+ factor
	- if(imode==0) {
	- Complex a1fact(a1BreitWigner(q2)*sqrt(2.)/3.);
	- if(ichan<0) {
	- F1 = -a1fact*TKstar1(s1,-1);
	- F2 = a1fact*Trho1(s2,-1);
	- F5 = Trho2(q2,-1)TOmegaKStar(s1,s2,-1)sqrt(2.);
	- }
	- else if(ichan%5==0) F1 = -a1fact*TKstar1(s1, ichan/5);
	- else if(ichan%5==1) F2 = a1fact*Trho1( s2,(ichan-1)/5);
	- else if(ichan%5>=2) F5 = Trho2(q2,ichan/5)TOmegaKStar(s1,s2,2((ichan-2)%5))
	- *sqrt(2.);
	- }
	- // calculate the K0 pi- K0bar
	- else if(imode==1) {
	- Complex a1fact(a1BreitWigner(q2)*sqrt(2.)/3.);
	- if(ichan<0) {
	- F1 =-a1fact*TKstar1(s1,-1);
	- F2 = a1fact*Trho1 (s2,-1);
	- F5 =-Trho2(q2,-1)TOmegaKStar(s1,s2,-1)sqrt(2.);
	- }
	- else if(ichan%5==0) F1 = -a1fact*TKstar1(s1, ichan/5);
	- else if(ichan%5==1) F2 = a1fact*Trho1 (s2,(ichan-1)/5);
	- else if(ichan%5>=2) F5 = -Trho2(q2,ichan/5)TOmegaKStar(s1,s2,2((ichan-2)%5))
	- *sqrt(2.);
	- }
	- // calculate the K- pi0 k0
	- else if(imode==2) {
	- Complex a1fact(a1BreitWigner(q2)/3.);
	- if(ichan<0) {
	- F1 = a1fact*( TKstar1(s1,-1)-TKstar1(s3,-1));
	- F2 = -a1fact(2.Trho1(s2,-1)+TKstar1(s3,-1));
	- F5 = Trho2(q2,-1)*(TKstar1(s3,-1)-TKstar1(s1,-1))/(1.+_omegaKstarwgt)/sqrt(2.);
	- }
	- else if(ichan%9==0) F1 = a1fact*TKstar1(s1,ichan/9)/3.;
	- else if(ichan%9==1) {
	- F1 = +a1fact*TKstar1(s3,(ichan-1)/9)/3.;
	- F2 = -a1fact*TKstar1(s3,(ichan-1)/9)/3.;
	- }
	- else if(ichan%9==2) F2 = -a1fact2.Trho1(s2,(ichan-2)/9)/3.;
	- else if(ichan%9<6) F5 =-Trho2(q2,ichan/9)*TKstar1(s1,(ichan-3)%9)
	- /(1.+_omegaKstarwgt)/sqrt(2.);
	- else F5 = Trho2(q2,ichan/9)*TKstar1(s3,(ichan-6)%9)
	- /(1.+_omegaKstarwgt)/sqrt(2.);
	- }
	- // calculate the K_S0 pi- K_S0 or K_L0 pi- K_L0
	- else if(imode==3\|\|imode==4) {
	- Complex a1fact(a1BreitWigner(q2)/6.);
	- if(ichan<0) {
	- F1 = a1fact*(TKstar1(s1,-1)+TKstar1(s3,-1));
	- F2 = a1fact*TKstar1(s3,-1);
	- F5 = 0.5Trho2(q2,-1)(TOmegaKStar(s1,s2,-1)-TOmegaKStar(s3,s2,-1));
	- }
	- else if(ichan%8==0) F1=a1fact*TKstar1(s1,ichan/8);
	- else if(ichan%8==1) {
	- F1 = a1fact*TKstar1(s3,ichan/8);
	- F2 = a1fact*TKstar1(s3,ichan/8);
	- }
	- else if(ichan%8<5 ) F5 = -Trho2(q2,ichan/8)*TKstar1(s1,(ichan-2)%8)
	- /(1.+_omegaKstarwgt)/2.;
	- else F5 = Trho2(q2,ichan/8)*TKstar1(s3,(ichan-5)%8)
	- /(1.+_omegaKstarwgt)/2.;
	- }
	- else if(imode==5) {
	- Complex a1fact(a1BreitWigner(q2)/3./sqrt(2.));
	- if(ichan<0) {
	- F1 = -a1fact*(TKstar1(s1,-1)-TKstar1(s3,-1));
	- F2 = a1fact(2.Trho1(s2,-1)+TKstar1(s3,-1));
	- F5 = -Trho2(q2,-1)*(TOmegaKStar(s1,s2,-1)+TOmegaKStar(s3,s2,-1))/sqrt(2.);
	- }
	- else if(ichan%9==0) F1 =- a1fact*TKstar1(s1,ichan/9);
	- else if(ichan%9==1) {
	- F1 = a1fact*TKstar1(s3,ichan/9);
	- F2 = a1fact*TKstar1(s3,ichan/9);
	- }
	- else if(ichan%9==2) F2 = 2.a1factTrho1( s2,ichan/9);
	- else if(ichan%9<6 ) F5 = -sqrt(0.5)TKstar2(q2,ichan/9)
	- TOmegaKStar(s1,s2,2*((ichan-3)%9))/sqrt(2.);
	- else F5 = -sqrt(0.5)TKstar2(q2,ichan/9)
	- TOmegaKStar(s3,s2,2*((ichan-6)%9))/sqrt(2.);
	- }
	- return FormFactors(F1 / _fpi,
	- F2 / _fpi,
	- InvEnergy(),
	- InvEnergy(),
	- -F5 / sqr(Constants::twopi) / pow<3,1>(_fpi));
	-}
	-
	void TwoKaonOnePionCurrent::doinitrun() {
	// set up the running a_1 width
	inita1Width(0);
	WeakCurrent::doinitrun();
	}

	void TwoKaonOnePionCurrent::doupdate() {
	WeakCurrent::doupdate();
	// update running width if needed
	if ( !touched() ) return;
	if(_maxmass!=_maxcalc) inita1Width(-1);
	}

	double TwoKaonOnePionCurrent::
	threeBodyMatrixElement(const int , const Energy2 q2,
	const Energy2 s3, const Energy2 s2,
	const Energy2 s1, const Energy ,
	const Energy , const Energy ) const {
	Energy2 mpi2(sqr(_mpi));
	Complex propb(Trho1(s1,-1)),propa(Trho1(s2,-1));
	// the matrix element
	Energy2 output(ZERO);
	// first resonance
	output+= ((s1-4.mpi2)+0.25(s3-s2)(s3-s2)/q2)real(propb*conj(propb));
	// second resonance
	output+= ((s2-4.mpi2)+0.25(s3-s1)(s3-s1)/q2)real(propa*conj(propa));
	// the interference term
	output+= (0.5q2-s3-0.5mpi2+0.25(s3-s2)(s3-s1)/q2)real(propaconj(propb)+
	propb*conj(propa));
	return output / sqr(_rho1mass[0]);
	}
	-
	-Complex TwoKaonOnePionCurrent::Trho1(Energy2 q2,int ires) const {
	- Complex output(0.);
	- double norm(0.);
	- for(unsigned int ix=0,N=_rho1wgts.size();ix<N;++ix) norm+=_rho1wgts[ix];
	- if(ires<0) {
	- for(unsigned int ix=0,N=_rho1wgts.size();ix<N;++ix) {
	- output+=_rho1wgts[ix]*BWrho1(q2,ix);
	- }
	- }
	- else {
	- unsigned int temp(ires);
	- if(temp<_rho1wgts.size()) output=_rho1wgts[temp]*BWrho1(q2,temp);
	- }
	- return output/norm;
	-}
	-
	-Complex TwoKaonOnePionCurrent::Trho2(Energy2 q2,int ires) const {
	- Complex output(0.);
	- double norm(0.);
	- for(unsigned int ix=0,N=_rho2wgts.size();ix<N;++ix) norm+=_rho2wgts[ix];
	- if(ires<0) {
	- for(unsigned int ix=0,N=_rho2wgts.size();ix<N;++ix) {
	- output+=_rho2wgts[ix]*BWrho2(q2,ix);
	- }
	- }
	- else {
	- unsigned int temp(ires);
	- if(temp<_rho2wgts.size()) output=_rho2wgts[temp]*BWrho2(q2,temp);
	- }
	- return output/norm;
	-}
	-
	-Complex TwoKaonOnePionCurrent::TKstar1(Energy2 q2,int ires) const {
	- Complex output(0.);
	- double norm(0.);
	- for(unsigned int ix=0,N=_kstar1wgts.size();ix<N;++ix) norm+=_kstar1wgts[ix];
	- if(ires<0) {
	- for(unsigned int ix=0,N=_kstar1wgts.size();ix<N;++ix) {
	- output+=_kstar1wgts[ix]*BWKstar1(q2,ix);
	- }
	- }
	- else {
	- unsigned int temp(ires);
	- if(temp<_kstar1wgts.size()) output=_kstar1wgts[temp]*BWKstar1(q2,temp);
	- }
	- return output/norm;
	-}
	-
	-Complex TwoKaonOnePionCurrent::TKstar2(Energy2 q2,int ires) const {
	- Complex output(0.);
	- double norm(0.);
	- for(unsigned int ix=0,N=_kstar2wgts.size();ix<N;++ix) norm+=_kstar2wgts[ix];
	- if(ires<0) {
	- for(unsigned int ix=0,N=_kstar2wgts.size();ix<N;++ix) {
	- output+=_kstar2wgts[ix]*BWKstar2(q2,ix);
	- }
	- }
	- else {
	- unsigned int temp(ires);
	- if(temp<_kstar2wgts.size()) output=_kstar2wgts[temp]*BWKstar2(q2,temp);
	- }
	- return output/norm;
	-}
	-
	-Complex TwoKaonOnePionCurrent::BWrho1(Energy2 q2, unsigned int ires) const {
	- if(ires>=_rho1mass.size()) return 0.;
	- Energy mass = _rho1mass [ires];
	- Energy width = _rho1width[ires];
	- Energy q=sqrt(q2);
	- Energy pcm0 = Kinematics::pstarTwoBodyDecay(mass,_mpi,_mpi);
	- Energy pcm = q<=2.*_mpi ? ZERO : Kinematics::pstarTwoBodyDecay(q,_mpi,_mpi);
	- double ratio = Math::powi(pcm/pcm0, 3);
	- Energy gam(widthmassratio/q);
	- return sqr(mass)/(sqr(mass)-q2-Complex(0.,1.)massgam);
	-}
	-
	-Complex TwoKaonOnePionCurrent::BWrho2(Energy2 q2, unsigned int ires) const {
	- if(ires>=_rho2mass.size()) return 0.;
	- Energy mass = _rho2mass [ires];
	- Energy width = _rho2width[ires];
	- Energy q=sqrt(q2);
	- Energy pcm0 = Kinematics::pstarTwoBodyDecay(mass,_mpi,_mpi);
	- Energy pcm = q<=2.*_mpi ? ZERO : Kinematics::pstarTwoBodyDecay(q,_mpi,_mpi);
	- double ratio(pcm/pcm0);ratio=ratioratio;
	- Energy gam(widthmassratio/q);
	- return sqr(mass)/(sqr(mass)-q2-Complex(0.,1.)massgam);
	-}
	-
	-Complex TwoKaonOnePionCurrent::BWKstar1(Energy2 q2, unsigned int ires) const {
	- if(ires>=_kstar1mass.size()) return 0.;
	- Energy mass = _kstar1mass [ires];
	- Energy width = _kstar1width[ires];
	- Energy q=sqrt(q2);
	- Energy pcm0 = Kinematics::pstarTwoBodyDecay(mass,_mK,_mpi);
	- Energy pcm = q<=_mpi+_mK ? ZERO : Kinematics::pstarTwoBodyDecay(q,_mK,_mpi);
	- double ratio(pcm/pcm0);ratio=ratioratio;
	- Energy gam(widthmassratio/q);
	- return sqr(mass)/(sqr(mass)-q2-Complex(0.,1.)massgam);
	-}
	-
	-Complex TwoKaonOnePionCurrent::BWKstar2(Energy2 q2, unsigned int ires) const {
	- if(ires>=_kstar2mass.size()) return 0.;
	- Energy mass = _kstar2mass [ires];
	- Energy width = _kstar2width[ires];
	- Energy q=sqrt(q2);
	- Energy pcm0 = Kinematics::pstarTwoBodyDecay(mass,_mK,_mpi);
	- Energy pcm = q<=_mpi+_mK ? ZERO : Kinematics::pstarTwoBodyDecay(q,_mK,_mpi);
	- double ratio(pcm/pcm0);ratio=ratioratio;
	- Energy gam(widthmassratio/q);
	- return sqr(mass)/(sqr(mass)-q2-Complex(0.,1.)massgam);
	-}
	-
	-Complex TwoKaonOnePionCurrent::a1BreitWigner(Energy2 q2) const {
	- Complex ii(0.,1.);
	- Energy2 m2(_a1mass*_a1mass);
	- Energy q(sqrt(q2));
	- return m2/(m2-q2-iiqa1Width(q2));
	-}
	-
	-Energy TwoKaonOnePionCurrent::a1Width(Energy2 q2) const {
	- if(!_a1opt) return _a1mass_a1widthg(q2)/g(_a1mass*_a1mass)/sqrt(q2);
	- else return (*_a1runinter)(q2);
	-}
	-
	-double TwoKaonOnePionCurrent::g(Energy2 q2) const {
	- double output;
	- if(q2<9._mpi_mpi) {
	- output=0.;
	- }
	- else if(q2<sqr(_rho1mass[0]+_mpi)) {
	- double diff=(q2-9._mpi_mpi)/GeV2;
	-
	- output=4.1sqr(diff)diff(1.-3.3diff+5.8*sqr(diff));
	- }
	- else {
	- double ratio = q2/GeV2;
	- output = ratio(1.623+10.38/ratio-9.32/sqr(ratio)+0.65/(ratiosqr(ratio)));
	- }
	- return output;
	-}

	Complex TwoKaonOnePionCurrent::Tomega(Energy2 q2, int ires) const {
	double denom=(1.+_epsomega);
	Complex num(0.);
	if(ires<0) num=OmegaPhiBreitWigner(q2,0)+_epsomega*OmegaPhiBreitWigner(q2,1);
	else if(ires==0) num=OmegaPhiBreitWigner(q2,0);
	else num=OmegaPhiBreitWigner(q2,1);
	return num/denom;
	}

	-Complex TwoKaonOnePionCurrent::OmegaPhiBreitWigner(Energy2 q2, unsigned int ires) const {
	- Energy2 m2,mg;
	- if(ires==0) {
	- m2=sqr(_omegamass);
	- mg=_omegamass*_omegawidth;
	- }
	- else {
	- m2=sqr(_phimass);
	- mg=_phimass*_phiwidth;
	- }
	- return (-m2+Complex(0.,1.)mg)/(q2-m2+Complex(0.,1.)mg);
	-}
	-
	Complex TwoKaonOnePionCurrent::TOmegaKStar(Energy2 s1,Energy2 s2,int ires) const {
	Complex output;
	if(ires<0) output = _omegaKstarwgt*TKstar1(s1,-1)+Tomega(s2,-1);
	else if(ires%2==0) output = _omegaKstarwgt*TKstar1(s1,ires/2);
	else if(ires%2==1) output = Tomega(s2,ires/2);
	return output/(1.+_omegaKstarwgt);
	}


	// the hadronic currents
	vector<LorentzPolarizationVectorE>
	TwoKaonOnePionCurrent::current(tcPDPtr resonance,
	IsoSpin::IsoSpin Itotal, IsoSpin::I3 i3,
	const int imode, const int ichan, Energy & scale,
	const tPDVector & ,
	const vector<Lorentz5Momentum> & momenta,
	DecayIntegrator::MEOption) const {
	+ // check the isospin
	+ if(Itotal!=IsoSpin::IUnknown && Itotal!=IsoSpin::IOne)
	+ return vector<LorentzPolarizationVectorE>();
	+ // check I_3
	+ if(i3!=IsoSpin::I3Unknown) {
	+ switch(i3) {
	+ case IsoSpin::I3Zero:
	+ return vector<LorentzPolarizationVectorE>();
	+ break;
	+ case IsoSpin::I3One: case IsoSpin::I3MinusOne:
	+ break;
	+ default:
	+ return vector<LorentzPolarizationVectorE>();
	+ }
	+ }
	+ // check the resonance
	+ int ires1=-1;
	+ if(resonance) {
	+ switch(abs(resonance->id())/1000) {
	+ case 0:
	+ ires1=0; break;
	+ case 100:
	+ ires1=1; break;
	+ case 30:
	+ ires1=2; break;
	+ case 10:
	+ ires1=3; break;
	+ default:
	+ assert(false);
	+ }
	+ }
	+ useMe();
	// calculate q2,s1,s2,s3
	Lorentz5Momentum q;
	for(unsigned int ix=0;ix<momenta.size();++ix)
	q+=momenta[ix];
	q.rescaleMass();
	scale=q.mass();
	Energy2 q2=q.mass2();
	Energy2 s1 = (momenta[1]+momenta[2]).m2();
	Energy2 s2 = (momenta[0]+momenta[2]).m2();
	Energy2 s3 = (momenta[0]+momenta[1]).m2();
	- FormFactors F = calculateFormFactors(ichan,imode,q2,s1,s2,s3);
	- //if(inpart.id()==ParticleID::tauplus){F.F5=conj(F.F5);}
	+ // calculate the form factors
	+ useMe();
	+ Complex F1(0.), F2(0.), F5(0.);
	+ Complex a1fact = ires1<0 \|\| ires1==3 ? a1BreitWigner(q2) : 0.;
	+ // calculate the K- pi - K+ factor
	+ if(imode==0) {
	+ a1fact *= sqrt(2.)/3.;
	+ if(ichan<0) {
	+ F1 = -a1fact*TKstar1(s1,-1);
	+ F2 = a1fact*Trho1(s2,-1);
	+ if(ires1<0)
	+ F5 = Trho2(q2, -1)TOmegaKStar(s1,s2,-1)sqrt(2.);
	+ else if(ires1<3)
	+ F5 = Trho2(q2,ires1)TOmegaKStar(s1,s2,-1)sqrt(2.);
	+ else
	+ F5 = 0.;
	+ }
	+ else if(ichan%5==0) F1 = -a1fact*TKstar1(s1, ichan/5);
	+ else if(ichan%5==1) F2 = a1fact*Trho1( s2,(ichan-1)/5);
	+ else if(ichan%5>=2) F5 = Trho2(q2,ichan/5)TOmegaKStar(s1,s2,2((ichan-2)%5))
	+ *sqrt(2.);
	+ }
	+ // calculate the K0 pi- K0bar
	+ else if(imode==1) {
	+ a1fact *= sqrt(2.)/3.;
	+ if(ichan<0) {
	+ F1 =-a1fact*TKstar1(s1,-1);
	+ F2 = a1fact*Trho1 (s2,-1);
	+ if(ires1<0)
	+ F5 =-Trho2(q2, -1)TOmegaKStar(s1,s2,-1)sqrt(2.);
	+ else if(ires1<3)
	+ F5 =-Trho2(q2,ires1)TOmegaKStar(s1,s2,-1)sqrt(2.);
	+ else
	+ F5 = 0.;
	+ }
	+ else if(ichan%5==0) F1 = -a1fact*TKstar1(s1, ichan/5);
	+ else if(ichan%5==1) F2 = a1fact*Trho1 (s2,(ichan-1)/5);
	+ else if(ichan%5>=2) F5 = -Trho2(q2,ichan/5)TOmegaKStar(s1,s2,2((ichan-2)%5))
	+ *sqrt(2.);
	+ }
	+ // calculate the K- pi0 k0
	+ else if(imode==2) {
	+ a1fact /= 3.;
	+ if(ichan<0) {
	+ F1 = a1fact*( TKstar1(s1,-1)-TKstar1(s3,-1));
	+ F2 = -a1fact(2.Trho1(s2,-1)+TKstar1(s3,-1));
	+ if(ires1<0)
	+ F5 = Trho2(q2, -1)*(TKstar1(s3,-1)-TKstar1(s1,-1))/(1.+_omegaKstarwgt)/sqrt(2.);
	+ else if(ires1<3)
	+ F5 = Trho2(q2,ires1)*(TKstar1(s3,-1)-TKstar1(s1,-1))/(1.+_omegaKstarwgt)/sqrt(2.);
	+ else
	+ F5 = 0.;
	+ }
	+ else if(ichan%9==0) F1 = a1fact*TKstar1(s1,ichan/9)/3.;
	+ else if(ichan%9==1) {
	+ F1 = +a1fact*TKstar1(s3,(ichan-1)/9)/3.;
	+ F2 = -a1fact*TKstar1(s3,(ichan-1)/9)/3.;
	+ }
	+ else if(ichan%9==2) F2 = -a1fact2.Trho1(s2,(ichan-2)/9)/3.;
	+ else if(ichan%9<6) F5 =-Trho2(q2,ichan/9)*TKstar1(s1,(ichan-3)%9)
	+ /(1.+_omegaKstarwgt)/sqrt(2.);
	+ else F5 = Trho2(q2,ichan/9)*TKstar1(s3,(ichan-6)%9)
	+ /(1.+_omegaKstarwgt)/sqrt(2.);
	+ }
	+ // calculate the K_S0 pi- K_S0 or K_L0 pi- K_L0
	+ else if(imode==3\|\|imode==4) {
	+ a1fact /=6;
	+ if(ichan<0) {
	+ F1 = a1fact*(TKstar1(s1,-1)+TKstar1(s3,-1));
	+ F2 = a1fact*TKstar1(s3,-1);
	+ if(ires1<0)
	+ F5 = 0.5Trho2(q2, -1)(TOmegaKStar(s1,s2,-1)-TOmegaKStar(s3,s2,-1));
	+ else if(ires1<3)
	+ F5 = 0.5Trho2(q2,ires1)(TOmegaKStar(s1,s2,-1)-TOmegaKStar(s3,s2,-1));
	+ else
	+ F5 = 0.;
	+ }
	+ else if(ichan%8==0) F1=a1fact*TKstar1(s1,ichan/8);
	+ else if(ichan%8==1) {
	+ F1 = a1fact*TKstar1(s3,ichan/8);
	+ F2 = a1fact*TKstar1(s3,ichan/8);
	+ }
	+ else if(ichan%8<5 ) F5 = -Trho2(q2,ichan/8)*TKstar1(s1,(ichan-2)%8)
	+ /(1.+_omegaKstarwgt)/2.;
	+ else F5 = Trho2(q2,ichan/8)*TKstar1(s3,(ichan-5)%8)
	+ /(1.+_omegaKstarwgt)/2.;
	+ }
	+ else if(imode==5) {
	+ a1fact *= 1./3./sqrt(2.);
	+ if(ichan<0) {
	+ F1 = -a1fact*(TKstar1(s1,-1)-TKstar1(s3,-1));
	+ F2 = a1fact(2.Trho1(s2,-1)+TKstar1(s3,-1));
	+ if(ires1<0)
	+ F5 = -Trho2(q2, -1)*(TOmegaKStar(s1,s2,-1)+TOmegaKStar(s3,s2,-1))/sqrt(2.);
	+ else if(ires1<3)
	+ F5 = -Trho2(q2,ires1)*(TOmegaKStar(s1,s2,-1)+TOmegaKStar(s3,s2,-1))/sqrt(2.);
	+ else
	+ F5 = 0.;
	+ }
	+ else if(ichan%9==0) F1 =- a1fact*TKstar1(s1,ichan/9);
	+ else if(ichan%9==1) {
	+ F1 = a1fact*TKstar1(s3,ichan/9);
	+ F2 = a1fact*TKstar1(s3,ichan/9);
	+ }
	+ else if(ichan%9==2) F2 = 2.a1factTrho1( s2,ichan/9);
	+ else if(ichan%9<6 ) F5 = -sqrt(0.5)Trho2(q2,ichan/9)
	+ TOmegaKStar(s1,s2,2*((ichan-3)%9))/sqrt(2.);
	+ else F5 = -sqrt(0.5)Trho2(q2,ichan/9)
	+ TOmegaKStar(s3,s2,2*((ichan-6)%9))/sqrt(2.);
	+ }
	// the first three form-factors
	- LorentzPolarizationVector vect;
	- vect = (F.F2-F.F1)*momenta[2]
	- +(F.F1-F.F3)*momenta[1]
	- +(F.F3-F.F2)*momenta[0];
	+ LorentzPolarizationVectorE vect = (F2-F1)momenta[2] + F1momenta[1] - F2*momenta[0];
	// multiply by the transverse projection operator
	- complex<InvEnergy> dot=(vect*q)/q2;
	+ Complex dot=(vect*q)/q2;
	// scalar and parity violating terms
	- vect += (F.F4-dot)*q;
	- if(F.F5!=complex<InvEnergy3>())
	- vect += Complex(0.,1.)F.F5Helicity::epsilon(momenta[0],
	- momenta[1],
	- momenta[2]);
	+ vect -= dot*q;
	+ if(F5!=0.)
	+ vect -= Complex(0.,1.)F5/sqr(Constants::twopi)/sqr(_fpi)
	+ Helicity::epsilon(momenta[0],momenta[1],momenta[2]);
	// factor to get dimensions correct
	- return vector<LorentzPolarizationVectorE>(1,q.mass()*vect);
	+ return vector<LorentzPolarizationVectorE>(1,q.mass()/_fpi*vect);
	}

	bool TwoKaonOnePionCurrent::accept(vector<int> id) {
	if(id.size()!=3) return false;
	int npip(0),npim(0),nkp(0),nkm(0);
	int npi0(0),nk0(0),nk0bar(0),nks(0),nkl(0);
	for(unsigned int ix=0;ix<id.size();++ix) {
	if(id[ix]==ParticleID::piplus) ++npip;
	else if(id[ix]==ParticleID::piminus) ++npim;
	else if(id[ix]==ParticleID::Kplus) ++nkp;
	else if(id[ix]==ParticleID::Kminus) ++nkm;
	else if(id[ix]==ParticleID::pi0) ++npi0;
	else if(id[ix]==ParticleID::K0) ++nk0;
	else if(id[ix]==ParticleID::Kbar0) ++nk0bar;
	else if(id[ix]==ParticleID::K_S0) ++nks;
	else if(id[ix]==ParticleID::K_L0) ++nkl;
	}
	if ( (nkp==1&&nkm==1&&npip==1) \|\|
	(nkp==1&&nkm==1&&npim==1)) return true;
	else if( (nk0==1&&nk0bar==1&&npip==1) \|\|
	(nk0==1&&nk0bar==1&&npim==1)) return true;
	else if( (nkp==1&&nk0bar==1&&npi0==1) \|\|
	(nkm==1&&npi0==1&&nk0==1)) return true;
	else if( nks==2 && (npip==1\|\|npim==1) ) return true;
	else if( nkl==2 && (npip==1\|\|npim==1) ) return true;
	else if( nks==1&&nkl==1 && (npip==1\|\|npim==1) ) return true;
	return false;
	}

	unsigned int TwoKaonOnePionCurrent::decayMode(vector<int> id) {
	assert(id.size()==3);
	int npip(0),npim(0),nkp(0),nkm(0),
	npi0(0),nk0(0),nk0bar(0),neta(0),nks(0),nkl(0);
	for(unsigned int ix=0;ix<id.size();++ix) {
	if(id[ix]==ParticleID::piplus) ++npip;
	else if(id[ix]==ParticleID::piminus) ++npim;
	else if(id[ix]==ParticleID::Kplus) ++nkp;
	else if(id[ix]==ParticleID::Kminus) ++nkm;
	else if(id[ix]==ParticleID::pi0) ++npi0;
	else if(id[ix]==ParticleID::K0) ++nk0;
	else if(id[ix]==ParticleID::Kbar0) ++nk0bar;
	else if(id[ix]==ParticleID::eta) ++neta;
	else if(id[ix]==ParticleID::K_S0) ++nks;
	else if(id[ix]==ParticleID::K_L0) ++nkl;
	}
	if ( (nkp==1&&nkm==1&&npip==1) \|\|
	(nkp==1&&nkm==1&&npim==1)) return 0;
	else if( (nk0==1&&nk0bar==1&&npip==1) \|\|
	(nk0==1&&nk0bar==1&&npim==1)) return 1;
	else if( (nkp==1&&nk0bar==1&&npi0==1) \|\|
	(nkm==1&&npi0==1&&nk0==1)) return 2;
	else if( nks==2 && (npip==1\|\|npim==1) ) return 3;
	else if( nkl==2 && (npip==1\|\|npim==1) ) return 4;
	else if( nks==1&&nkl==1 && (npip==1\|\|npim==1) ) return 5;
	assert(false);
	}


	tPDVector TwoKaonOnePionCurrent::particles(int icharge, unsigned int imode,int,int) {
	tPDVector extpart(3);
	if(imode==0) {
	extpart[0]=getParticleData(ParticleID::Kminus);
	extpart[1]=getParticleData(ParticleID::piminus);
	extpart[2]=getParticleData(ParticleID::Kplus);
	}
	else if(imode==1) {
	extpart[0]=getParticleData(ParticleID::K0);
	extpart[1]=getParticleData(ParticleID::piminus);
	extpart[2]=getParticleData(ParticleID::Kbar0);
	}
	else if(imode==2) {
	extpart[0]=getParticleData(ParticleID::Kminus);
	extpart[1]=getParticleData(ParticleID::pi0);
	extpart[2]=getParticleData(ParticleID::K0);
	}
	else if(imode==3) {
	extpart[0]=getParticleData(ParticleID::K_S0);
	extpart[1]=getParticleData(ParticleID::piminus);
	extpart[2]=getParticleData(ParticleID::K_S0);
	}
	else if(imode==4) {
	extpart[0]=getParticleData(ParticleID::K_L0);
	extpart[1]=getParticleData(ParticleID::piminus);
	extpart[2]=getParticleData(ParticleID::K_L0);
	}
	else if(imode==5) {
	extpart[0]=getParticleData(ParticleID::K_S0);
	extpart[1]=getParticleData(ParticleID::piminus);
	extpart[2]=getParticleData(ParticleID::K_L0);
	}
	// conjugate the particles if needed
	if(icharge==3) {
	for(unsigned int ix=0;ix<3;++ix) {
	if(extpart[ix]->CC()) extpart[ix]=extpart[ix]->CC();
	}
	}
	// return the answer
	return extpart;
	}
	diff --git a/Decay/WeakCurrents/TwoKaonOnePionCurrent.h b/Decay/WeakCurrents/TwoKaonOnePionCurrent.h
	--- a/Decay/WeakCurrents/TwoKaonOnePionCurrent.h
	+++ b/Decay/WeakCurrents/TwoKaonOnePionCurrent.h
	@@ -1,607 +1,536 @@
	// -- C++ --
	//
	// TwoKaonOnePionCurrent.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2017 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_TwoKaonOnePionCurrent_H
	#define HERWIG_TwoKaonOnePionCurrent_H
	//
	// This is the declaration of the TwoKaonOnePionCurrent class.
	//

	#include "WeakCurrent.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* The TwoKaonOnePionCurrent class implements the model of M. Finkemeier
	* and E.~Mirkes, Z. Phys. C 69 (1996) 243 [arXiv:hep-ph/9503474],
	* for the weak current for three mesons where at least one of the mesons is
	* a kaon.
	*
	* \ingroup Decay
	*
	* This is the base class for the three meson decays of the weak current.
	* It is designed so that the currents for the following modes can be implemented
	* in classes inheriting from this
	* - \f$ K^- \pi^- K^+ \f$, (imode=0)
	* - \f$ K^0 \pi^- \bar{K}^0\f$, (imode=1)
	* - \f$ K^- \pi^0 K^0 \f$, (imode=2)
	* - \f$ \pi^0 \pi^0 K^- \f$, (imode=3)
	* - \f$ K^- \pi^- \pi^+ \f$, (imode=4)
	* - \f$ \pi^- \bar{K}^0 \pi^0 \f$, (imode=5)
	* - \f$ \pi^- \pi^0 \eta \f$, (imode=6)
	*
	* obviously there are other modes with three pseudoscalar mesons for the decay
	* of the weak current but this model original came from \f$\tau\f$ decay where
	* these are the only modes. However one case which is important is the inclusion
	* of the mixing in the neutral kaon sector for which we include the additional
	* currents
	* - \f$ K^0_S \pi^- K^0_S\f$, (imode=9)
	* - \f$ K^0_L \pi^- K^0_L\f$, (imode=10)
	* - \f$ K^0_S \pi^- K^0_L\f$, (imode=11)
	*
	* In this case the current is given by
	* \f[ J^\mu = \left(g^{\mu\nu}-\frac{q^\mu q^\nu}{q^2}\right)
	* \left[F_1(p_2-p_3)^\mu +F_2(p_3-p_1)^\mu+F_3(p_1-p_2)^\mu\right]
	* +q^\mu F_4
	* +F_5\epsilon^{\mu\alpha\beta\gamma}p_1^\alpha p_2^\beta p_3^\gamma
	* \f]
	* where
	* - \f$p_{1,2,3}\f$ are the momenta of the mesons in the order given above.
	* - \f$F_1,F_2,F_3,F_4,F_5\f$ are the form factors which must be
	* calculated in the calculateFormFactors member which should be implemented
	* in classes inheriting from this.
	*
	* @see WeakCurrent.
	*
	* \author Peter Richardson
	* @see \ref TwoKaonOnePionCurrentInterfaces "The interfaces"
	* defined for TwoKaonOnePionCurrent.
	*/
	class TwoKaonOnePionCurrent: public WeakCurrent {

	public:

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

	/** @name Methods for the construction of the phase space integrator. */
	//@{
	/**
	* Complete the construction of the decay mode for integration.classes inheriting
	* from this one.
	* This method is purely virtual and must be implemented in the classes inheriting
	* from WeakCurrent.
	* @param icharge The total charge of the outgoing particles in the current.
	* @param resonance If specified only include terms with this particle
	* @param Itotal If specified the total isospin of the current
	* @param I3 If specified the thrid component of isospin
	* @param imode The mode in the current being asked for.
	* @param mode The phase space mode for the integration
	* @param iloc The location of the of the first particle from the current in
	* the list of outgoing particles.
	* @param ires The location of the first intermediate for the current.
	* @param phase The prototype phase space channel for the integration.
	* @param upp The maximum possible mass the particles in the current are
	* allowed to have.
	* @return Whether the current was sucessfully constructed.
	*/
	virtual bool createMode(int icharge, tcPDPtr resonance,
	IsoSpin::IsoSpin Itotal, IsoSpin::I3 i3,
	unsigned int imode,PhaseSpaceModePtr mode,
	unsigned int iloc,int ires,
	PhaseSpaceChannel phase, Energy upp );
	//@}


	/**
	* Hadronic current. This method is purely virtual and must be implemented in
	* all classes inheriting from this one.
	* @param resonance If specified only include terms with this particle
	* @param Itotal If specified the total isospin of the current
	* @param I3 If specified the thrid component of isospin
	* @param imode The mode
	* @param ichan The phase-space channel the current is needed for.
	* @param scale The invariant mass of the particles in the current.
	* @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 current.
	*/
	virtual vector<LorentzPolarizationVectorE>
	current(tcPDPtr resonance,
	IsoSpin::IsoSpin Itotal, IsoSpin::I3 i3,
	const int imode, const int ichan,Energy & scale,
	const tPDVector & outgoing,
	const vector<Lorentz5Momentum> & momenta,
	DecayIntegrator::MEOption meopt) const;

	/**
	* Accept the decay. Checks the mesons against the list.
	* @param id The id's of the particles in the current.
	* @return Can this current have the external particles specified.
	*/
	virtual bool accept(vector<int> id);

	/**
	* Return the decay mode number for a given set of particles in the current.
	* Checks the mesons against the list.
	* @param id The id's of the particles in the current.
	* @return The number of the mode
	*/
	virtual unsigned int decayMode(vector<int> id);

	/**
	* The particles produced by the current. This returns the mesons for the mode.
	* @param icharge The total charge of the particles in the current.
	* @param imode The mode for which the particles are being requested
	* @param iq The PDG code for the quark
	* @param ia The PDG code for the antiquark
	* @return The external particles for the current.
	*/
	virtual tPDVector particles(int icharge, unsigned int imode, int iq, int ia);

	/**
	* 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
	* @param create Whether or not to add a statement creating the object
	*/
	virtual void dataBaseOutput(ofstream & os,bool header,bool create) const;

	/**
	* the matrix element for the \f$a_1\f$ decay to calculate the running width
	* @param imode The mode for which the matrix element is needed.
	* @param q2 The mass of the decaying off-shell \f$a_1\f$, \f$q^2\f$.
	* @param s3 The invariant mass squared of particles 1 and 2, \f$s_3=m^2_{12}\f$.
	* @param s2 The invariant mass squared of particles 1 and 3, \f$s_2=m^2_{13}\f$.
	* @param s1 The invariant mass squared of particles 2 and 3, \f$s_1=m^2_{23}\f$.
	* @param m1 The mass of the first outgoing particle.
	* @param m2 The mass of the second outgoing particle.
	* @param m3 The mass of the third outgoing particle.
	* @return The matrix element squared summed over spins.
	*/
	double threeBodyMatrixElement(const int imode, const Energy2 q2,
	const Energy2 s3, const Energy2 s2,
	const Energy2 s1, const Energy m1,
	const Energy m2, const Energy m3) const;
	-
	-protected:
	-
	- /**
	- * Helper class for form factors
	- */
	- struct FormFactors {
	-
	- /**
	- * @param F1 The \f$F_1\f$ form factor
	- */
	- complex<InvEnergy> F1;
	-
	- /**
	- * @param F2 The \f$F_2\f$ form factor
	- */
	- complex<InvEnergy> F2;
	-
	- /**
	- * @param F3 The \f$F_3\f$ form factor
	- */
	- complex<InvEnergy> F3;
	-
	- /**
	- * @param F4 The \f$F_4\f$ form factor
	- */
	- complex<InvEnergy> F4;
	-
	- /**
	- * @param F5 The \f$F_5\f$ form factor
	- */
	- complex<InvEnergy3> F5;
	-
	- /**
	- * Constructor
	- * @param f1 The \f$F_1\f$ form factor
	- * @param f2 The \f$F_2\f$ form factor
	- * @param f3 The \f$F_3\f$ form factor
	- * @param f4 The \f$F_4\f$ form factor
	- * @param f5 The \f$F_5\f$ form factor
	- */
	- FormFactors(complex<InvEnergy> f1 = InvEnergy(),
	- complex<InvEnergy> f2 = InvEnergy(),
	- complex<InvEnergy> f3 = InvEnergy(),
	- complex<InvEnergy> f4 = InvEnergy(),
	- complex<InvEnergy3> f5 = InvEnergy3())
	- : F1(f1), F2(f2), F3(f3), F4(f4), F5(f5) {}
	- };
	-
	- /**
	- * Calculate the form factor for the current. Implements the form factors
	- * described above.
	- * @param ichan The phase space channel
	- * @param imode The mode
	- * @param q2 The scale \f$q^2\f$ for the current.
	- * @param s1 The invariant mass squared of particles 2 and 3, \f$s_1=m^2_{23}\f$.
	- * @param s2 The invariant mass squared of particles 1 and 3, \f$s_2=m^2_{13}\f$.
	- * @param s3 The invariant mass squared of particles 1 and 2, \f$s_3=m^2_{12}\f$.
	- */
	- virtual FormFactors calculateFormFactors(const int ichan,const int imode,
	- Energy2 q2,Energy2 s1,Energy2 s2,Energy2 s3) 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 Clone Methods. */
	//@{
	/**
	* Make a simple clone of this object.
	* @return a pointer to the new object.
	*/
	virtual IBPtr clone() const;

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

	protected:

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

	/**
	* Initialize this object to the begining of the run phase.
	*/
	virtual void doinitrun();

	/**
	* Check sanity of the object during the setup phase.
	*/
	virtual void doupdate();
	//@}

	private:

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

	private:

	/**
	* The \f$\rho\f$ lineshape for the axial-vector terms
	* @param q2 The scale \f$q^2\f$ for the lineshape
	* @param ires Which \f$\rho\f$ multiplet
	*/
	- Complex Trho1(Energy2 q2,int ires) const;
	+ Complex Trho1(Energy2 q2,int ires) const {
	+ if(ires>=int(_rho1wgts.size())) return 0.;
	+ double norm = std::accumulate(_rho1wgts.begin(),_rho1wgts.end(),0.);
	+ unsigned int imin=0,imax=_rho1wgts.size();
	+ if(ires>0) {
	+ imin=ires;
	+ imax=imin+1;
	+ }
	+ Complex output(0.);
	+ for(unsigned int ix=imin;ix<imax;++ix)
	+ output+=_rho1wgts[ix]*
	+ Resonance::BreitWignerPWave(q2,_rho1mass[ix],_rho1width[ix],_mpi,_mpi);
	+ return output/norm;
	+ }

	/**
	* The \f$\rho\f$ lineshape for the vector terms
	* @param q2 The scale \f$q^2\f$ for the lineshape
	* @param ires Which \f$\rho\f$ multiplet
	*/
	- Complex Trho2(Energy2 q2,int ires) const;
	+ Complex Trho2(Energy2 q2,int ires) const {
	+ if(ires>=int(_rho2wgts.size())) return 0.;
	+ double norm = std::accumulate(_rho2wgts.begin(),_rho2wgts.end(),0.);
	+ unsigned int imin=0,imax=_rho2wgts.size();
	+ if(ires>0) {
	+ imin=ires;
	+ imax=imin+1;
	+ }
	+ Complex output(0.);
	+ for(unsigned int ix=imin;ix<imax;++ix)
	+ output+=_rho2wgts[ix]*
	+ Resonance::BreitWignerPWave(q2,_rho2mass[ix],_rho2width[ix],_mpi,_mpi);
	+ return output/norm;
	+ }

	/**
	* The \f$K^*\f$ lineshape for the axial-vector terms
	* @param q2 The scale \f$q^2\f$ for the lineshape
	* @param ires Which \f$K^*\f$ multiplet
	*/
	- Complex TKstar1(Energy2 q2,int ires) const;
	-
	- /**
	- * The \f$\rho\f$ lineshape for the vector terms
	- * @param q2 The scale \f$q^2\f$ for the lineshape
	- * @param ires Which \f$K^*\f$ multiplet
	- */
	- Complex TKstar2(Energy2 q2,int ires) const;
	-
	- /**
	- * The \f$\rho\f$ Breit-Wigner for the axial-vector terms
	- * @param q2 The scale \f$q^2\f$ for the Breit-Wigner
	- * @param ires Which \f$\rho\f$ multiplet
	- */
	- Complex BWrho1(Energy2 q2, unsigned int ires) const;
	-
	- /**
	- * The \f$\rho\f$ Breit-Wigner for the vector terms
	- * @param q2 The scale \f$q^2\f$ for the Breit-Wigner
	- * @param ires Which \f$\rho\f$ multiplet
	- */
	- Complex BWrho2(Energy2 q2, unsigned int ires) const;
	-
	- /**
	- * The \f$K^*\f$ Breit-Wigner for the axial-vector terms
	- * @param q2 The scale \f$q^2\f$ for the Breit-Wigner
	- * @param ires Which \f$K^*\f$ multiplet
	- */
	- Complex BWKstar1(Energy2 q2, unsigned int ires) const;
	-
	- /**
	- * The \f$K^*\f$ Breit-Wigner for the vector terms
	- * @param q2 The scale \f$q^2\f$ for the Breit-Wigner
	- * @param ires Which \f$K^*\f$ multiplet
	- */
	- Complex BWKstar2(Energy2 q2, unsigned int ires) const;
	+ Complex TKstar1(Energy2 q2,int ires) const {
	+ if(ires>=int(_kstar1wgts.size())) return 0.;
	+ double norm = std::accumulate(_kstar1wgts.begin(),_kstar1wgts.end(),0.);
	+ unsigned int imin=0,imax=_kstar1wgts.size();
	+ if(ires>0) {
	+ imin=ires;
	+ imax=imin+1;
	+ }
	+ Complex output(0.);
	+ for(unsigned int ix=imin;ix<imax;++ix)
	+ output+=_kstar1wgts[ix]*
	+ Resonance::BreitWignerPWave(q2,_kstar1mass[ix],_kstar1width[ix],_mK,_mpi);
	+ return output/norm;
	+ }

	/**
	* \f$a_1\f$ Breit-Wigner
	* @param q2 The scale \f$q^2\f$ for the Breit-Wigner
	* @return The Breit-Wigner
	*/
	- Complex a1BreitWigner(Energy2 q2) const;
	-
	- /**
	- * The \f$a_1\f$ running width
	- * @param q2 The scale \f$q^2\f$ for the Breit-Wigner
	- * @return The \f$a_1\f$ running width.
	- */
	- Energy a1Width(Energy2 q2) const;
	-
	- /**
	- * The \f$g(Q^2)\f$ function of Kuhn and Santamaria
	- */
	- double g(Energy2 q2) const;
	+ Complex a1BreitWigner(Energy2 q2) const {
	+ Complex ii(0.,1.);
	+ Energy2 m2(_a1mass*_a1mass);
	+ Energy q(sqrt(q2));
	+ Energy gamma = !_a1opt ?
	+ _a1mass_a1widthResonance::ga1(q2)/Resonance::ga1(_a1mass_a1mass)/sqrt(q2) : (_a1runinter)(q2);
	+ return m2/(m2-q2-iiqgamma);
	+ }

	/**
	* Initialize the \f$a_1\f$ running width
	* @param iopt Initialization option (-1 full calculation, 0 set up the interpolation)
	*/
	void inita1Width(int iopt);

	/**
	* The \f$T_\omega\f$ function
	* @param q2 The scale
	* @param ires the resonance
	*/
	Complex Tomega(Energy2 q2, int ires) const;

	/**
	* The \f$\omega\f$ and \f$\phi\f$ Breit-Wigner
	* @param q2 The scale
	* @param ires the resonance
	*/
	- Complex OmegaPhiBreitWigner(Energy2 q2, unsigned int ires) const;
	+ Complex OmegaPhiBreitWigner(Energy2 q2, unsigned int ires) const {
	+ Energy2 m2,mg;
	+ if(ires==0) {
	+ m2=sqr(_omegamass);
	+ mg=_omegamass*_omegawidth;
	+ }
	+ else {
	+ m2=sqr(_phimass);
	+ mg=_phimass*_phiwidth;
	+ }
	+ return (-m2+Complex(0.,1.)mg)/(q2-m2+Complex(0.,1.)mg);
	+ }

	/**
	* The \f$\omega-\phi\f$ \f$K^*\f$ form-factor for the \f$F_5\f$ form-factor
	* @param s1 The scale \f$s_1\f$.
	* @param s2 The scale \f$s_2\f$.
	* @param ires Which resonances to use
	* @return The mixed Breit-Wigner
	*/
	Complex TOmegaKStar(Energy2 s1,Energy2 s2,int ires) const;

	private:

	/**
	* Parameters for the \f$\rho\f$ in the axial-vector terms
	*/
	//@{
	/**
	* Weight for the different resonances
	*/
	vector<double> _rho1wgts;

	/**
	* Masses
	*/
	vector<Energy> _rho1mass;

	/**
	* Widths
	*/
	vector<Energy> _rho1width;
	//@}

	/**
	* Parameters for the \f$\rho\f$ in the vector terms
	*/
	//@{
	/**
	* Weight for the different resonances
	*/
	vector<double> _rho2wgts;

	/**
	* Masses
	*/
	vector<Energy> _rho2mass;

	/**
	* Widths
	*/
	vector<Energy> _rho2width;
	//@}

	/**
	* Parameters for the \f$K^*\f$ in the axial-vector terms
	*/
	//@{
	/**
	* Weight for the different resonances
	*/
	vector<double> _kstar1wgts;

	/**
	* Masses
	*/
	vector<Energy> _kstar1mass;

	/**
	* Widths
	*/
	vector<Energy> _kstar1width;
	//@}

	/**
	- * Parameters for the \f$K^*\f$ in the vector terms
	- */
	- //@{
	- /**
	- * Weight for the different resonances
	- */
	- vector<double> _kstar2wgts;
	-
	- /**
	- * Masses
	- */
	- vector<Energy> _kstar2mass;
	-
	- /**
	- * Widths
	- */
	- vector<Energy> _kstar2width;
	- //@}
	-
	- /**
	* Parameters for the three meson resonances
	*/
	//@{
	/**
	* The mass of the \f$a_1\f$ resonances.
	*/
	Energy _a1mass;

	/**
	* The width of the \f$a_1\f$ resonances.
	*/
	Energy _a1width;

	/**
	* The \f$a_1\f$ width for the running \f$a_1\f$ width calculation.
	*/
	vector<Energy> _a1runwidth;

	/**
	* The \f$q^2\f$ for the running \f$a_1\f$ width calculation.
	*/
	vector<Energy2> _a1runq2;

	/**
	* The interpolator for the running \f$a_1\f$ width calculation.
	*/
	Interpolator<Energy,Energy2>::Ptr _a1runinter;
	//@}

	/**
	* Parameters for the \f$T_\omega\f$ function
	*/
	//@{
	/**
	* Mixing parameter
	*/
	double _epsomega;

	/**
	* Mass of the \f$\omega\f$
	*/
	Energy _omegamass;

	/**
	* Width of the \f$\omega\f$
	*/
	Energy _omegawidth;

	/**
	* Mass of the \f$\phi\f$
	*/
	Energy _phimass;

	/**
	* Width of the \f$\phi\f$
	*/
	Energy _phiwidth;
	//@}

	/**
	* The relative weight of the \f$\omega-\phi\f$ and \f$K^*\f$ where needed.
	*/
	double _omegaKstarwgt;

	/**
	* The pion decay constant, \f$f_\pi\f$.
	*/
	Energy _fpi;

	/**
	* The pion mass
	*/
	Energy _mpi;

	/**
	* The kaon mass
	*/
	Energy _mK;

	/**
	* Initialization switches
	*/
	//@{
	/**
	* Initialize the running \f$a_1\f$ width.
	*/
	bool _initializea1;

	/**
	* Option for the \f$a_1\f$ width
	*/
	bool _a1opt;
	//@}

	/**
	* The maximum mass of the hadronic system
	*/
	Energy _maxmass;

	/**
	* The maximum mass when the running width was calculated
	*/
	Energy _maxcalc;
	};

	}

	#endif /* HERWIG_TwoKaonOnePionCurrent_H */
	diff --git a/Decay/WeakCurrents/TwoKaonOnePionDefaultCurrent.cc b/Decay/WeakCurrents/TwoKaonOnePionDefaultCurrent.cc
	--- a/Decay/WeakCurrents/TwoKaonOnePionDefaultCurrent.cc
	+++ b/Decay/WeakCurrents/TwoKaonOnePionDefaultCurrent.cc
	@@ -1,754 +1,754 @@
	// -- C++ --
	//
	// TwoKaonOnePionDefaultCurrent.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2017 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 TwoKaonOnePionDefaultCurrent class.
	//

	#include "TwoKaonOnePionDefaultCurrent.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 "Herwig/PDT/ThreeBodyAllOnCalculator.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "ThePEG/Helicity/epsilon.h"

	using namespace Herwig;
	using namespace ThePEG;

	DescribeClass<TwoKaonOnePionDefaultCurrent,WeakCurrent>
	describeHerwigTwoKaonOnePionDefaultCurrent("Herwig::TwoKaonOnePionDefaultCurrent",
	"HwWeakCurrents.so");
	HERWIG_INTERPOLATOR_CLASSDESC(TwoKaonOnePionDefaultCurrent,Energy,Energy2)


	TwoKaonOnePionDefaultCurrent::TwoKaonOnePionDefaultCurrent() {
	// the quarks for the different modes
	addDecayMode(2,-1);
	addDecayMode(2,-1);
	addDecayMode(2,-1);
	setInitialModes(3);
	// the pion decay constant
	_fpi=130.7*MeV/sqrt(2.);
	_mpi = ZERO;
	_mK = ZERO;
	// set the initial weights for the resonances
	// the rho weights
	_rhoF123wgts = { 1.0,-0.145,0.};
	_rhoF5wgts = {-26., 6.5,1.};
	// the Kstar weights
	_kstarF123wgts = {1.};
	// relative rho/Kstar weights
	_rhoKstarwgt=-0.2;
	// local values of the a_1 parameters
	_a1mass = 1.251*GeV;
	_a1width = 0.599*GeV;
	_a1opt=true;
	// local values of the rho parameters
	_rhoF123masses = {0.773GeV,1.370GeV,1.750*GeV};
	_rhoF123widths = {0.145GeV,0.510GeV,0.120*GeV};
	_rhoF5masses = {0.773GeV,1.500GeV,1.750*GeV};
	_rhoF5widths = {0.145GeV,0.220GeV,0.120*GeV};
	// local values for the Kstar parameters
	_kstarF123masses = {0.8921*GeV};
	_kstarF123widths = {0.0513*GeV};
	// initialization of the a_1 running width
	_initializea1=false;
	double a1q2in[200]={0,15788.6,31577.3,47365.9,63154.6,78943.2,94731.9,110521,
	126309,142098,157886,173675,189464,205252,221041,236830,
	252618,268407,284196,299984,315773,331562,347350,363139,
	378927,394716,410505,426293,442082,457871,473659,489448,
	505237,521025,536814,552603,568391,584180,599969,615757,
	631546,647334,663123,678912,694700,710489,726278,742066,
	757855,773644,789432,805221,821010,836798,852587,868375,
	884164,899953,915741,931530,947319,963107,978896,994685,
	1.01047e+06,1.02626e+06,1.04205e+06,1.05784e+06,1.07363e+06,
	1.08942e+06,1.10521e+06,1.12099e+06,1.13678e+06,1.15257e+06,
	1.16836e+06,1.18415e+06,1.19994e+06,1.21573e+06,1.23151e+06,
	1.2473e+06,1.26309e+06,1.27888e+06,1.29467e+06,1.31046e+06,
	1.32625e+06,1.34203e+06,1.35782e+06,1.37361e+06,1.3894e+06,
	1.40519e+06,1.42098e+06,1.43677e+06,1.45256e+06,1.46834e+06
	,1.48413e+06,1.49992e+06,1.51571e+06,1.5315e+06,1.54729e+06,
	1.56308e+06,1.57886e+06,1.59465e+06,1.61044e+06,1.62623e+06,
	1.64202e+06,1.65781e+06,1.6736e+06,1.68939e+06,1.70517e+06,
	1.72096e+06,1.73675e+06,1.75254e+06,1.76833e+06,1.78412e+06,
	1.79991e+06,1.81569e+06,1.83148e+06,1.84727e+06,1.86306e+06,
	1.87885e+06,1.89464e+06,1.91043e+06,1.92621e+06,1.942e+06,
	1.95779e+06,1.97358e+06,1.98937e+06,2.00516e+06,2.02095e+06,
	2.03674e+06,2.05252e+06,2.06831e+06,2.0841e+06,2.09989e+06,
	2.11568e+06,2.13147e+06,2.14726e+06,2.16304e+06,2.17883e+06,
	2.19462e+06,2.21041e+06,2.2262e+06,2.24199e+06,2.25778e+06,
	2.27356e+06,2.28935e+06,2.30514e+06,2.32093e+06,2.33672e+06,
	2.35251e+06,2.3683e+06,2.38409e+06,2.39987e+06,2.41566e+06,
	2.43145e+06,2.44724e+06,2.46303e+06,2.47882e+06,2.49461e+06,
	2.51039e+06,2.52618e+06,2.54197e+06,2.55776e+06,2.57355e+06,
	2.58934e+06,2.60513e+06,2.62092e+06,2.6367e+06,2.65249e+06,
	2.66828e+06,2.68407e+06,2.69986e+06,2.71565e+06,2.73144e+06,
	2.74722e+06,2.76301e+06,2.7788e+06,2.79459e+06,2.81038e+06,
	2.82617e+06,2.84196e+06,2.85774e+06,2.87353e+06,2.88932e+06,
	2.90511e+06,2.9209e+06,2.93669e+06,2.95248e+06,2.96827e+06,
	2.98405e+06,2.99984e+06,3.01563e+06,3.03142e+06,3.04721e+06,
	3.063e+06,3.07879e+06,3.09457e+06,3.11036e+06,3.12615e+06,
	3.14194e+06};
	double a1widthin[200]={0,0,0,0,0,0,0,0,0,0,0,0,0.00153933,0.0136382,0.0457614,
	0.105567,0.199612,0.333825,0.513831,0.745192,1.0336,1.38501,
	1.80581,2.30295,2.88403,3.5575,4.33278,5.22045,6.23243,
	7.38223,8.68521,10.1589,11.8234,13.7018,15.8206,18.2107,
	20.9078,23.9533,27.3954,31.2905,35.7038,40.7106,46.3984,
	52.8654,60.2207,68.581,78.0637,88.7754,100.794,114.145,
	128.783,144.574,161.299,178.683,196.426,214.248,231.908,
	249.221,266.059,282.336,298.006,313.048,327.46,341.254,
	354.448,367.066,379.133,390.677,401.726,412.304,422.439,
	432.155,441.474,450.419,459.01,467.267,475.207,482.847,
	490.203,497.29,504.121,510.71,517.068,523.207,529.138,
	534.869,540.411,545.776,550.961,556.663,560.851,565.566,
	570.137,574.569,578.869,583.041,587.091,591.023,594.843,
	598.553,602.16,605.664,609.072,612.396,615.626,618.754,
	621.796,624.766,627.656,630.47,633.21,635.878,638.5,
	641.006,643.471,645.873,648.213,650.493,652.715,654.88,
	656.99,659.047,661.052,663.007,664.963,666.771,668.6,
	670.351,672.075,673.828,675.397,676.996,678.567,680.083,
	681.589,683.023,684.457,685.825,687.18,688.499,689.789,
	691.058,692.284,693.501,694.667,695.82,696.947,698.05,
	699.129,700.186,701.221,702.234,703.226,704.198,705.158,
	706.085,707.001,707.899,708.78,709.644,710.474,711.334,
	712.145,712.943,713.727,714.505,715.266,716.015,716.751,
	717.474,718.183,718.88,719.645,720.243,720.91,721.565,
	722.211,722.851,723.473,724.094,724.697,725.296,725.886,
	726.468,727.041,727.608,728.166,728.718,729.262,729.808,
	730.337,730.856,731.374,731.883,732.386,732.884,733.373,
	733.859,734.339,734.813};

	vector<double> tmp1(a1widthin,a1widthin+200);
	_a1runwidth.clear();
	std::transform(tmp1.begin(), tmp1.end(),
	back_inserter(_a1runwidth),
	[](double x){return x*MeV;});

	vector<double> tmp2(a1q2in,a1q2in+200);
	_a1runq2.clear();
	std::transform(tmp2.begin(), tmp2.end(),
	back_inserter(_a1runq2),
	[](double x){return x*MeV2;});
	_maxmass=ZERO;
	_maxcalc=ZERO;
	}

	void TwoKaonOnePionDefaultCurrent::doinit() {
	WeakCurrent::doinit();
	// masses for the running widths
	_mpi = getParticleData(ParticleID::piplus)->mass();
	_mK = getParticleData(ParticleID::Kminus)->mass();
	// initialise the a_1 running width calculation
	inita1Width(-1);
	}

	void TwoKaonOnePionDefaultCurrent::persistentOutput(PersistentOStream & os) const {
	os << _rhoF123wgts << _kstarF123wgts << _rhoF5wgts
	<< _rhoKstarwgt << ounit(_a1runwidth,GeV)<< ounit(_a1runq2,GeV2)
	<< _initializea1 << ounit(_a1mass,GeV)<< ounit(_a1width,GeV)
	<< ounit(_fpi,GeV) << ounit(_mpi,GeV)<< ounit(_mK,GeV)
	<< ounit(_rhoF123masses,GeV) << ounit(_rhoF5masses,GeV)
	<< ounit(_rhoF123widths,GeV)
	<< ounit(_rhoF5widths,GeV) << ounit(_kstarF123masses,GeV)
	<< ounit(_kstarF123widths,GeV)
	<< _a1opt << ounit(_maxmass,GeV) << ounit(_maxcalc,GeV) << _a1runinter;
	}

	void TwoKaonOnePionDefaultCurrent::persistentInput(PersistentIStream & is, int) {
	is >> _rhoF123wgts >> _kstarF123wgts >> _rhoF5wgts
	>> _rhoKstarwgt >> iunit(_a1runwidth,GeV) >> iunit(_a1runq2,GeV2)
	>> _initializea1 >> iunit(_a1mass,GeV) >> iunit(_a1width,GeV)
	>> iunit(_fpi,GeV) >> iunit(_mpi,GeV) >> iunit(_mK,GeV)
	>> iunit(_rhoF123masses,GeV) >> iunit(_rhoF5masses,GeV)
	>> iunit(_rhoF123widths,GeV)
	>> iunit(_rhoF5widths,GeV) >> iunit(_kstarF123masses,GeV)
	>> iunit(_kstarF123widths,GeV)
	>> _a1opt >> iunit(_maxmass,GeV) >> iunit(_maxcalc,GeV) >> _a1runinter;
	}

	void TwoKaonOnePionDefaultCurrent::Init() {

	static ClassDocumentation<TwoKaonOnePionDefaultCurrent> documentation
	("The TwoKaonOnePionDefaultCurrent class is designed to implement "
	"the three meson decays of the tau, ie pi- pi- pi+, pi0 pi0 pi-, "
	"K- pi- K+, K0 pi- Kbar0, K- pi0 K0,pi0 pi0 K-, K- pi- pi+, "
	"pi- Kbar0 pi0, pi- pi0 eta. It uses the same currents as those in TAUOLA.",
	"The three meson decays of the tau, ie pi- pi- pi+, pi0 pi0 pi-, "
	"K- pi- K+, K0 pi- Kbar0, K- pi0 K0,pi0 pi0 K-, K- pi- pi+, "
	"and pi- Kbar0 pi0, pi- pi0 eta "
	"use the same currents as \\cite{Jadach:1993hs,Kuhn:1990ad,Decker:1992kj}.",
	"%\\cite{Jadach:1993hs}\n"
	"\\bibitem{Jadach:1993hs}\n"
	" S.~Jadach, Z.~Was, R.~Decker and J.~H.~Kuhn,\n"
	" %``The Tau Decay Library Tauola: Version 2.4,''\n"
	" Comput.\\ Phys.\\ Commun.\\ {\\bf 76}, 361 (1993).\n"
	" %%CITATION = CPHCB,76,361;%%\n"
	"%\\cite{Kuhn:1990ad}\n"
	"\\bibitem{Kuhn:1990ad}\n"
	" J.~H.~Kuhn and A.~Santamaria,\n"
	" %``Tau decays to pions,''\n"
	" Z.\\ Phys.\\ C {\\bf 48}, 445 (1990).\n"
	" %%CITATION = ZEPYA,C48,445;%%\n"
	"%\\cite{Decker:1992kj}\n"
	"\\bibitem{Decker:1992kj}\n"
	" R.~Decker, E.~Mirkes, R.~Sauer and Z.~Was,\n"
	" %``Tau decays into three pseudoscalar mesons,''\n"
	" Z.\\ Phys.\\ C {\\bf 58}, 445 (1993).\n"
	" %%CITATION = ZEPYA,C58,445;%%\n"
	);

	static ParVector<TwoKaonOnePionDefaultCurrent,double> interfaceF123RhoWgt
	("F123RhoWeight",
	"The weights of the different rho resonances in the F1,2,3 form factor",
	&TwoKaonOnePionDefaultCurrent::_rhoF123wgts,
	0, 0, 0, -1000, 1000, false, false, true);

	static ParVector<TwoKaonOnePionDefaultCurrent,double> interfaceF123KstarWgt
	("F123KstarWeight",
	"The weights of the different Kstar resonances in the F1,2,3 form factor",
	&TwoKaonOnePionDefaultCurrent::_kstarF123wgts,
	0, 0, 0, -1000, 1000, false, false, true);

	static ParVector<TwoKaonOnePionDefaultCurrent,double> interfaceF5RhoWgt
	("F5RhoWeight",
	"The weights of the different rho resonances in the F1,2,3 form factor",
	&TwoKaonOnePionDefaultCurrent::_rhoF5wgts,
	0, 0, 0, -1000, 1000, false, false, true);

	static Parameter<TwoKaonOnePionDefaultCurrent,double> interfaceRhoKstarWgt
	("RhoKstarWgt",
	"The relative weights of the rho and K* in the F5 form factor",
	&TwoKaonOnePionDefaultCurrent::_rhoKstarwgt, -0.2, -10., 10.,
	false, false, false);

	static Switch<TwoKaonOnePionDefaultCurrent,bool> interfaceInitializea1
	("Initializea1",
	"Initialise the calculation of the a_1 running width",
	&TwoKaonOnePionDefaultCurrent::_initializea1, false, false, false);
	static SwitchOption interfaceInitializea1Initialization
	(interfaceInitializea1,
	"Yes",
	"Initialize the calculation",
	true);
	static SwitchOption interfaceInitializea1NoInitialization
	(interfaceInitializea1,
	"No",
	"Use the default values",
	false);

	static Switch<TwoKaonOnePionDefaultCurrent,bool> interfacea1WidthOption
	("a1WidthOption",
	"Option for the treatment of the a1 width",
	&TwoKaonOnePionDefaultCurrent::_a1opt, true, false, false);
	static SwitchOption interfacea1WidthOptionLocal
	(interfacea1WidthOption,
	"Local",
	"Use a calculation of the running width based on the parameters as"
	" interpolation table.",
	true);
	static SwitchOption interfacea1WidthOptionParam
	(interfacea1WidthOption,
	"Kuhn",
	"Use the parameterization of Kuhn and Santamaria for default parameters."
	" This should only be used for testing vs TAUOLA",
	false);

	static ParVector<TwoKaonOnePionDefaultCurrent,Energy> interfacea1RunningWidth
	("a1RunningWidth",
	"The values of the a_1 width for interpolation to giving the running width.",
	&TwoKaonOnePionDefaultCurrent::_a1runwidth, GeV, -1, 1.0GeV, ZERO, 10.0GeV,
	false, false, true);

	static ParVector<TwoKaonOnePionDefaultCurrent,Energy2> interfacea1RunningQ2
	("a1RunningQ2",
	"The values of the q^2 for interpolation to giving the running width.",
	&TwoKaonOnePionDefaultCurrent::_a1runq2, GeV2, -1, 1.0GeV2, ZERO, 10.0GeV2,
	false, false, true);

	static Parameter<TwoKaonOnePionDefaultCurrent,Energy> interfaceA1Width
	("A1Width",
	"The a_1 width if using local values.",
	&TwoKaonOnePionDefaultCurrent::_a1width, GeV, 0.599GeV, ZERO, 10.0GeV,
	false, false, false);

	static Parameter<TwoKaonOnePionDefaultCurrent,Energy> interfaceA1Mass
	("A1Mass",
	"The a_1 mass if using local values.",
	&TwoKaonOnePionDefaultCurrent::_a1mass, GeV, 1.251GeV, ZERO, 10.0GeV,
	false, false, false);

	static ParVector<TwoKaonOnePionDefaultCurrent,Energy> interfacerhoF123masses
	("rhoF123masses",
	"The masses for the rho resonances if used local values",
	&TwoKaonOnePionDefaultCurrent::_rhoF123masses, GeV, -1, 1.0GeV, ZERO, 10.0GeV,
	false, false, true);

	static ParVector<TwoKaonOnePionDefaultCurrent,Energy> interfacerhoF123widths
	("rhoF123widths",
	"The widths for the rho resonances if used local values",
	&TwoKaonOnePionDefaultCurrent::_rhoF123widths, GeV, -1, 1.0GeV, ZERO, 10.0GeV,
	false, false, true);

	static ParVector<TwoKaonOnePionDefaultCurrent,Energy> interfacerhoF5masses
	("rhoF5masses",
	"The masses for the rho resonances if used local values",
	&TwoKaonOnePionDefaultCurrent::_rhoF5masses, GeV, -1, 1.0GeV, ZERO, 10.0GeV,
	false, false, true);

	static ParVector<TwoKaonOnePionDefaultCurrent,Energy> interfacerhoF5widths
	("rhoF5widths",
	"The widths for the rho resonances if used local values",
	&TwoKaonOnePionDefaultCurrent::_rhoF5widths, GeV, -1, 1.0GeV, ZERO, 10.0GeV,
	false, false, true);

	static ParVector<TwoKaonOnePionDefaultCurrent,Energy> interfaceKstarF123masses
	("KstarF123masses",
	"The masses for the Kstar resonances if used local values",
	&TwoKaonOnePionDefaultCurrent::_kstarF123masses, GeV, -1, 1.0GeV, ZERO, 10.0GeV,
	false, false, true);

	static ParVector<TwoKaonOnePionDefaultCurrent,Energy> interfaceKstarF123widths
	("KstarF123widths",
	"The widths for the Kstar resonances if used local values",
	&TwoKaonOnePionDefaultCurrent::_kstarF123widths, GeV, -1, 1.0GeV, ZERO, 10.0GeV,
	false, false, true);

	static Parameter<TwoKaonOnePionDefaultCurrent,Energy> interfaceFPi
	("FPi",
	"The pion decay constant",
	&TwoKaonOnePionDefaultCurrent::_fpi, MeV, 92.4MeV, ZERO, 200.0MeV,
	false, false, true);
	}

	// complete the construction of the decay mode for integration
	bool TwoKaonOnePionDefaultCurrent::createMode(int icharge, tcPDPtr resonance,
	IsoSpin::IsoSpin Itotal, IsoSpin::I3 i3,
	unsigned int imode,PhaseSpaceModePtr mode,
	unsigned int iloc,int ires,
	PhaseSpaceChannel phase, Energy upp ) {
	// check the charge
	if(abs(icharge)!=3) return false;
	// check the total isospin
	if(Itotal!=IsoSpin::IUnknown) {
	if(Itotal!=IsoSpin::IOne) return false;
	}
	// check I_3
	if(i3!=IsoSpin::I3Unknown) {
	switch(i3) {
	case IsoSpin::I3Zero:
	if(imode<=1) return false;
	break;
	case IsoSpin::I3One:
	if( imode>1 \|\| icharge ==-3) return false;
	break;
	case IsoSpin::I3MinusOne:
	if( imode>1 \|\| icharge == 3) return false;
	default:
	return false;
	}
	}
	// get the particles and check the mass
	int iq(0),ia(0);
	tPDVector extpart(particles(1,imode,iq,ia));
	Energy min(ZERO);
	for(unsigned int ix=0;ix<extpart.size();++ix) min+=extpart[ix]->massMin();
	if(min>upp) return false;
	// the particles we will use a lot
	tPDPtr a1=getParticleData(ParticleID::a_1minus);
	if(icharge==3) a1=a1->CC();
	_maxmass=max(_maxmass,upp);
	// the rho0 resonances
	tPDPtr rho0[3] = { getParticleData(113), getParticleData(100113), getParticleData(30113)};
	tPDPtr rhoc[3] = {getParticleData(-213),getParticleData(-100213),getParticleData(-30213)};
	tPDPtr Kstar0[3] = { getParticleData(313), getParticleData(100313), getParticleData(30313)};
	tPDPtr Kstarc[3] = {getParticleData(-323),getParticleData(-100323),getParticleData(-30323)};
	if(icharge==3) {
	for(unsigned int ix=0;ix<3;++ix) {
	rhoc [ix] = rhoc[ix]->CC();
	Kstar0[ix] = Kstar0[ix]->CC();
	Kstarc[ix] = Kstarc[ix]->CC();
	}
	}
	if(imode==0) {
	// channels for K- pi- K+
	for(unsigned int ix=0;ix<3;++ix) {
	if(!resonance \|\| resonance==a1) {
	mode->addChannel((PhaseSpaceChannel(phase),ires,a1,ires+1,iloc+1,ires+1,Kstar0[ix],
	ires+2,iloc+2,ires+2,iloc+3));
	mode->addChannel((PhaseSpaceChannel(phase),ires,a1,ires+1,iloc+2,ires+1,rho0[ix],
	ires+2,iloc+1,ires+2,iloc+3));
	}
	for(unsigned int iy=0;iy<3;++iy) {
	if(resonance && resonance !=rhoc[ix]) continue;
	mode->addChannel((PhaseSpaceChannel(phase),ires,rhoc[ix],ires+1,iloc+1,ires+1,Kstar0[iy],
	ires+2,iloc+2,ires+2,iloc+3));
	mode->addChannel((PhaseSpaceChannel(phase),ires,rhoc[ix],ires+1,iloc+2,ires+1,rho0[iy],
	ires+2,iloc+1,ires+2,iloc+3));
	}
	}
	}
	else if(imode==1) {
	// channels for K0 pi- K0bar
	for(unsigned int ix=0;ix<3;++ix) {
	if(!resonance \|\| resonance==a1) {
	mode->addChannel((PhaseSpaceChannel(phase),ires,a1,ires+1,iloc+1,ires+1,Kstarc[ix],
	ires+2,iloc+2,ires+2,iloc+3));
	mode->addChannel((PhaseSpaceChannel(phase),ires,a1,ires+1,iloc+2,ires+1,rho0[ix],
	ires+2,iloc+1,ires+2,iloc+3));
	}
	for(unsigned int iy=0;iy<3;++iy) {
	if(resonance && resonance !=rhoc[ix]) continue;
	mode->addChannel((PhaseSpaceChannel(phase),ires,rhoc[ix],ires+1,iloc+1,ires+1,Kstarc[iy],
	ires+2,iloc+2,ires+2,iloc+3));
	mode->addChannel((PhaseSpaceChannel(phase),ires,rhoc[ix],ires+1,iloc+2,ires+1,rho0[iy],
	ires+2,iloc+1,ires+2,iloc+3));
	}
	}
	}
	else if(imode==2) {
	// channels for K- pi0 K0
	for(unsigned int ix=0;ix<3;++ix) {
	if(!resonance \|\| resonance==a1) {
	mode->addChannel((PhaseSpaceChannel(phase),ires,a1,ires+1,iloc+2,ires+1,rhoc[ix],
	ires+2,iloc+1,ires+2,iloc+3));
	}
	}
	}
	for(unsigned int ix=0;ix<_rhoF123masses.size();++ix) {
	mode->resetIntermediate(rhoc[ix],_rhoF123masses[ix],
	_rhoF123widths[ix]);
	mode->resetIntermediate(rho0[ix],_rhoF123masses[ix],
	_rhoF123widths[ix]);
	}
	// K star parameters in the base class
	for(unsigned int ix=0;ix<_kstarF123masses.size();++ix) {
	mode->resetIntermediate(Kstarc[ix],_kstarF123masses[ix],
	_kstarF123widths[ix]);
	mode->resetIntermediate(Kstar0[ix],_kstarF123masses[ix],
	_kstarF123widths[ix]);
	}
	return true;
	}

	// initialisation of the a_1 width
	// (iopt=-1 initialises, iopt=0 starts the interpolation)
	void TwoKaonOnePionDefaultCurrent::inita1Width(int iopt) {
	if(iopt==-1) {
	_maxcalc=_maxmass;
	if(!_initializea1\|\|_maxmass==ZERO) return;
	// parameters for the table of values
	Energy2 step(sqr(_maxcalc)/199.);
	// integrator to perform the integral
	vector<double> inweights;inweights.push_back(0.5);inweights.push_back(0.5);
	vector<int> intype;intype.push_back(2);intype.push_back(3);
	Energy mrho(getParticleData(ParticleID::rhoplus)->mass()),
	wrho(getParticleData(ParticleID::rhoplus)->width());
	vector<Energy> inmass(2,mrho),inwidth(2,wrho);
	vector<double> inpow(2,0.0);
	ThreeBodyAllOnCalculator<TwoKaonOnePionDefaultCurrent>
	widthgen(inweights,intype,inmass,inwidth,inpow,*this,0,_mpi,_mpi,_mpi);
	// normalisation constant to give physical width if on shell
	double a1const(_a1width/(widthgen.partialWidth(sqr(_a1mass))));
	// loop to give the values
	_a1runq2.clear(); _a1runwidth.clear();
	for(Energy2 moff2(ZERO); moff2<=sqr(_maxcalc); moff2+=step) {
	_a1runwidth.push_back(widthgen.partialWidth(moff2)*a1const);
	_a1runq2.push_back(moff2);
	}
	}
	// set up the interpolator
	else if(iopt==0) {
	_a1runinter = make_InterpolatorPtr(_a1runwidth,_a1runq2,3);
	}
	}

	void TwoKaonOnePionDefaultCurrent::dataBaseOutput(ofstream & output,bool header,
	bool create) const {
	if(header) output << "update decayers set parameters=\"";
	if(create) output << "create Herwig::TwoKaonOnePionDefaultCurrent "
	<< name() << " HwWeakCurrents.so\n";
	for(unsigned int ix=0;ix<_rhoF123wgts.size();++ix) {
	if(ix<3) output << "newdef ";
	else output << "insert ";
	output << name() << ":F123RhoWeight " << ix << " " << _rhoF123wgts[ix] << "\n";
	}
	for(unsigned int ix=0;ix<_kstarF123wgts.size();++ix) {
	if(ix<1) output << "newdef ";
	else output << "insert ";
	output << name() << ":F123KstarWeight " << ix << " "
	<< _kstarF123wgts[ix] << "\n";
	}
	for(unsigned int ix=0;ix<_rhoF5wgts.size();++ix) {
	if(ix<3) output << "newdef ";
	else output << "insert ";
	output << name() << ":F5RhoWeight " << ix << " " << _rhoF5wgts[ix] << "\n";
	}
	output << "newdef " << name() << ":RhoKstarWgt " << _rhoKstarwgt << "\n";
	output << "newdef " << name() << ":Initializea1 " << _initializea1 << "\n";
	output << "newdef " << name() << ":a1WidthOption " << _a1opt << "\n";
	for(unsigned int ix=0;ix<_a1runwidth.size();++ix) {
	output << "newdef " << name() << ":a1RunningWidth " << ix
	<< " " << _a1runwidth[ix]/GeV << "\n";
	}
	for(unsigned int ix=0;ix<_a1runq2.size();++ix) {
	output << "newdef " << name() << ":a1RunningQ2 " << ix
	<< " " << _a1runq2[ix]/GeV2 << "\n";
	}
	output << "newdef " << name() << ":A1Width " << _a1width/GeV << "\n";
	output << "newdef " << name() << ":A1Mass " << _a1mass/GeV << "\n";
	output << "newdef " << name() << ":FPi " << _fpi/MeV << "\n";
	for(unsigned int ix=0;ix<_rhoF123masses.size();++ix) {
	if(ix<3) output << "newdef ";
	else output << "insert ";
	output << name() << ":rhoF123masses " << ix
	<< " " << _rhoF123masses[ix]/GeV << "\n";
	}
	for(unsigned int ix=0;ix<_rhoF123widths.size();++ix) {
	if(ix<3) output << "newdef ";
	else output << "insert ";
	output << name() << ":rhoF123widths " << ix << " "
	<< _rhoF123widths[ix]/GeV << "\n";
	}
	for(unsigned int ix=0;ix<_rhoF5masses.size();++ix) {
	if(ix<3) output << "newdef ";
	else output << "insert ";
	output << name() << ":rhoF5masses " << ix << " "
	<< _rhoF5masses[ix]/GeV << "\n";
	}
	for(unsigned int ix=0;ix<_rhoF5widths.size();++ix) {
	if(ix<3) output << "newdef ";
	else output << "insert ";
	output << name() << ":rhoF5widths " << ix << " "
	<< _rhoF5widths[ix]/GeV << "\n";
	}
	for(unsigned int ix=0;ix<_kstarF123masses.size();++ix) {
	if(ix<1) output << "newdef ";
	else output << "insert ";
	output << name() << ":KstarF123masses " << ix << " "
	<< _kstarF123masses[ix]/GeV << "\n";
	}
	for(unsigned int ix=0;ix<_kstarF123widths.size();++ix) {
	if(ix<1) output << "newdef ";
	else output << "insert ";
	output << name() << ":KstarF123widths " << ix << " "
	<< _kstarF123widths[ix]/GeV << "\n";
	}
	WeakCurrent::dataBaseOutput(output,false,false);
	if(header) output << "\n\" where BINARY ThePEGName=\""
	<< fullName() << "\";" << endl;
	}

	void TwoKaonOnePionDefaultCurrent::doinitrun() {
	// set up the running a_1 width
	inita1Width(0);
	WeakCurrent::doinitrun();
	}

	void TwoKaonOnePionDefaultCurrent::doupdate() {
	WeakCurrent::doupdate();
	// update running width if needed
	if ( !touched() ) return;
	if(_maxmass!=_maxcalc) inita1Width(-1);
	}

	double TwoKaonOnePionDefaultCurrent::
	threeBodyMatrixElement(const int , const Energy2 q2,
	const Energy2 s3, const Energy2 s2,
	const Energy2 s1, const Energy ,
	const Energy , const Energy ) const {
	Energy2 mpi2(sqr(_mpi));
	Complex propb(BrhoF123(s1,-1)),propa(BrhoF123(s2,-1));
	// the matrix element
	Energy2 output(ZERO);
	// first resonance
	output += ((s1-4.mpi2) + 0.25(s3-s2)(s3-s2)/q2) real(propb*conj(propb));
	// second resonance
	output += ((s2-4.mpi2) + 0.25(s3-s1)(s3-s1)/q2) real(propa*conj(propa));
	// the interference term
	output += (0.5q2-s3-0.5mpi2+0.25(s3-s2)(s3-s1)/q2)real(propaconj(propb)+
	propb*conj(propa));
	return output/sqr(_rhoF123masses[0]);
	}

	// the hadronic currents
	vector<LorentzPolarizationVectorE>
	TwoKaonOnePionDefaultCurrent::current(tcPDPtr resonance,
	- IsoSpin::IsoSpin Itotal, IsoSpin::I3 i3,
	- const int imode, const int ichan, Energy & scale,
	- const tPDVector & ,
	- const vector<Lorentz5Momentum> & momenta,
	- DecayIntegrator::MEOption) const {
	+ IsoSpin::IsoSpin Itotal, IsoSpin::I3 i3,
	+ const int imode, const int ichan, Energy & scale,
	+ const tPDVector & ,
	+ const vector<Lorentz5Momentum> & momenta,
	+ DecayIntegrator::MEOption) const {
	// check the isospin
	if(Itotal!=IsoSpin::IUnknown && Itotal!=IsoSpin::IOne)
	return vector<LorentzPolarizationVectorE>();
	// check I_3
	if(i3!=IsoSpin::I3Unknown) {
	switch(i3) {
	case IsoSpin::I3Zero:
	return vector<LorentzPolarizationVectorE>();
	break;
	case IsoSpin::I3One: case IsoSpin::I3MinusOne:
	break;
	default:
	return vector<LorentzPolarizationVectorE>();
	}
	}
	// check the resonance
	int ires1=-1;
	if(resonance) {
	switch(abs(resonance->id())/1000) {
	case 0:
	ires1=0; break;
	case 100:
	ires1=1; break;
	case 30:
	ires1=2; break;
	case 10:
	ires1=3; break;
	default:
	assert(false);
	}
	}
	useMe();
	// calculate q2,s1,s2,s3
	Lorentz5Momentum q;
	for(unsigned int ix=0;ix<momenta.size();++ix)
	q+=momenta[ix];
	q.rescaleMass();
	scale=q.mass();
	Energy2 q2=q.mass2();
	Energy2 s1 = (momenta[1]+momenta[2]).m2();
	Energy2 s2 = (momenta[0]+momenta[2]).m2();
	Complex F1(0.), F2(0.), F5(0.);
	Complex a1fact = ires1<0 \|\| ires1==3 ? a1BreitWigner(q2) : 0.;
	// calculate the K- pi - K+ factor
	if(imode==0) {
	a1fact *= sqrt(2.)/3.;
	if(ichan<0) {
	F1 =-a1fact*BKstarF123(s1,-1);
	F2 = a1fact*BrhoF123(s2,-1);
	if(ires1<0)
	F5 = BrhoF5(q2,-1)FKrho(s1,s2,-1)sqrt(2.);
	else if(ires1<3)
	F5 = BrhoF5(q2,ires1)FKrho(s1,s2,-1)sqrt(2.);
	else
	F5 = 0.;
	}
	else if(ichan%8==0) F1 =-a1fact*BKstarF123(s1,ichan/8);
	else if(ichan%8==1) F2 = a1fact*BrhoF123(s2,(ichan-1)/8);
	else if(ichan%8>=2) F5 = BrhoF5(q2,ichan/8)FKrho(s1,s2,(ichan-2)%8)sqrt(2.);
	}
	// calculate the K0 pi- K0bar
	else if(imode==1) {
	a1fact *= sqrt(2.)/3.;
	if(ichan<0) {
	F1 =-a1fact*BKstarF123(s1,-1);
	F2 = a1fact*BrhoF123(s2,-1);
	if(ires1<0)
	F5 =-BrhoF5(q2,-1)FKrho(s1,s2,-1)sqrt(2.);
	else if(ires1<3)
	F5 =-BrhoF5(q2,ires1)FKrho(s1,s2,-1)sqrt(2.);
	else
	F5 = 0.;
	}
	else if(ichan%8==0) F1 = -a1fact*BKstarF123(s1,ichan/8);
	else if(ichan%8==1) F2 = a1fact*BrhoF123(s2,(ichan-1)/8);
	else if(ichan%8>=2) F5 = -BrhoF5(q2,ichan/8)FKrho(s1,s2,(ichan-2)%8)sqrt(2.);
	}
	// calculate the K- pi0 k0
	else if(imode==2) {
	if(ichan<0) F2 =-a1fact*BrhoF123(s2,-1);
	else F2 =-a1fact*BrhoF123(s2,ichan);
	}
	// the first three form-factors
	LorentzPolarizationVectorE vect =
	(F2-F1)momenta[2] + F1momenta[1] - F2*momenta[0];
	// multiply by the transverse projection operator
	Complex dot=(vect*q)/q2;
	// scalar and parity violating terms
	vect -= dot*q;
	if(F5!=0.)
	vect -= Complex(0.,1.)F5/sqr(Constants::twopi)/sqr(_fpi)
	Helicity::epsilon(momenta[0],momenta[1],momenta[2]);
	// factor to get dimensions correct
	return vector<LorentzPolarizationVectorE>(1,q.mass()/_fpi*vect);
	}

	bool TwoKaonOnePionDefaultCurrent::accept(vector<int> id) {
	int npip(0),npim(0),nkp(0),nkm(0);
	int npi0(0),nk0(0),nk0bar(0);
	for(unsigned int ix=0;ix<id.size();++ix) {
	if(id[ix]==ParticleID::piplus) ++npip;
	else if(id[ix]==ParticleID::piminus) ++npim;
	else if(id[ix]==ParticleID::Kplus) ++nkp;
	else if(id[ix]==ParticleID::Kminus) ++nkm;
	else if(id[ix]==ParticleID::pi0) ++npi0;
	else if(id[ix]==ParticleID::K0) ++nk0;
	else if(id[ix]==ParticleID::Kbar0) ++nk0bar;
	}
	if ( (nkp==1&&nkm==1&&npip==1) \|\|
	(nkp==1&&nkm==1&&npim==1)) return true;
	else if( (nk0==1&&nk0bar==1&&npip==1) \|\|
	(nk0==1&&nk0bar==1&&npim==1)) return true;
	else if( (nkp==1&&nk0bar==1&&npi0==1) \|\|
	(nkm==1&&npi0==1&&nk0==1)) return true;
	return false;
	}

	unsigned int TwoKaonOnePionDefaultCurrent::decayMode(vector<int> id) {
	int npip(0),npim(0),nkp(0),nkm(0);
	int npi0(0),nk0(0),nk0bar(0);
	for(unsigned int ix=0;ix<id.size();++ix) {
	if(id[ix]==ParticleID::piplus) ++npip;
	else if(id[ix]==ParticleID::piminus) ++npim;
	else if(id[ix]==ParticleID::Kplus) ++nkp;
	else if(id[ix]==ParticleID::Kminus) ++nkm;
	else if(id[ix]==ParticleID::pi0) ++npi0;
	else if(id[ix]==ParticleID::K0) ++nk0;
	else if(id[ix]==ParticleID::Kbar0) ++nk0bar;
	}
	if ( (nkp==1&&nkm==1&&npip==1) \|\|
	(nkp==1&&nkm==1&&npim==1)) return 0;
	else if( (nk0==1&&nk0bar==1&&npip==1) \|\|
	(nk0==1&&nk0bar==1&&npim==1)) return 1;
	else if( (nkp==1&&nk0bar==1&&npi0==1) \|\|
	(nkm==1&&npi0==1&&nk0==1)) return 2;
	assert(false);
	}

	tPDVector TwoKaonOnePionDefaultCurrent::particles(int icharge, unsigned int imode,int,int) {
	tPDVector extpart(3);
	if(imode==0) {
	extpart[0]=getParticleData(ParticleID::Kminus);
	extpart[1]=getParticleData(ParticleID::piminus);
	extpart[2]=getParticleData(ParticleID::Kplus);
	}
	else if(imode==1) {
	extpart[0]=getParticleData(ParticleID::K0);
	extpart[1]=getParticleData(ParticleID::piminus);
	extpart[2]=getParticleData(ParticleID::Kbar0);
	}
	else if(imode==2) {
	extpart[0]=getParticleData(ParticleID::Kminus);
	extpart[1]=getParticleData(ParticleID::pi0);
	extpart[2]=getParticleData(ParticleID::K0);
	}
	// conjugate the particles if needed
	if(icharge==3) {
	for(unsigned int ix=0;ix<3;++ix) {
	if(extpart[ix]->CC()) extpart[ix]=extpart[ix]->CC();
	}
	}
	// return the answer
	return extpart;
	}

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 5:20 PM (1 d, 14 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3805369
Default Alt Text: (123 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions