D119.id502.diff
No OneTemporary
Actions

Size

62 KB

Subscribers

None

D119.id502.diff
View Options

	diff --git a/EvtGenBase/EvtDalitzReso.hh b/EvtGenBase/EvtDalitzReso.hh
	--- a/EvtGenBase/EvtDalitzReso.hh
	+++ b/EvtGenBase/EvtDalitzReso.hh
	@@ -87,11 +87,11 @@
	WA76
	};

	- EvtDalitzReso() : _typeN( NON_RES ){};
	+ EvtDalitzReso() : m_typeN( NON_RES ){};

	EvtDalitzReso( const EvtDalitzPlot& dp, EvtCyclic3::Pair pairRes,
	NumType typeN, double alpha = 0.0 ) :
	- _dp( dp ), _pairRes( pairRes ), _typeN( typeN ), _alpha( alpha ){};
	+ m_dp( dp ), m_pairRes( pairRes ), m_typeN( typeN ), m_alpha( alpha ){};

	EvtDalitzReso( const EvtDalitzPlot& dp, EvtCyclic3::Pair pairAng,
	EvtCyclic3::Pair pairRes, EvtSpinType::spintype spin,
	@@ -124,106 +124,109 @@

	EvtDalitzReso* clone() const { return new EvtDalitzReso( *this ); }

	- EvtComplex evaluate( const EvtDalitzPoint& p );
	+ EvtComplex evaluate( const EvtDalitzPoint& p ) const;

	- void set_fd( double R ) { _vd.set_f( R ); }
	- void set_fb( double R ) { _vb.set_f( R ); }
	+ void set_fd( double R ) { m_vd.set_f( R ); }
	+ void set_fb( double R ) { m_vb.set_f( R ); }

	void addFlatteParam( const EvtFlatteParam& param )
	{
	- _flatteParams.push_back( param );
	+ m_flatteParams.push_back( param );
	}

	private:
	- EvtComplex psFactor( double& ma, double& mb, double& m );
	- EvtComplex psFactor( double& ma1, double& mb1, double& ma2, double& mb2,
	- double& m );
	- EvtComplex propGauss( const double& m0, const double& s0, const double& m );
	+ EvtComplex psFactor( const double& ma, const double& mb,
	+ const double& m ) const;
	+ EvtComplex psFactor( const double& ma1, const double& mb1, const double& ma2,
	+ const double& mb2, const double& m ) const;
	+ EvtComplex propGauss( const double& m0, const double& s0,
	+ const double& m ) const;
	EvtComplex propBreitWigner( const double& m0, const double& g0,
	- const double& m );
	+ const double& m ) const;
	EvtComplex propBreitWignerRel( const double& m0, const double& g0,
	- const double& m );
	+ const double& m ) const;
	EvtComplex propBreitWignerRel( const double& m0, const EvtComplex& g0,
	- const double& m );
	+ const double& m ) const;
	EvtComplex propBreitWignerRelCoupled( const double& m0, const EvtComplex& g1,
	- const EvtComplex& g2, const double& m );
	+ const EvtComplex& g2,
	+ const double& m ) const;
	EvtComplex propGounarisSakurai( const double& m0, const double& g0,
	const double& k0, const double& m,
	- const double& g, const double& k );
	+ const double& g, const double& k ) const;
	inline double GS_f( const double& m0, const double& g0, const double& k0,
	- const double& m, const double& k );
	- inline double GS_h( const double& m, const double& k );
	- inline double GS_dhods( const double& m0, const double& k0 );
	- inline double GS_d( const double& m0, const double& k0 );
	+ const double& m, const double& k ) const;
	+ inline double GS_h( const double& m, const double& k ) const;
	+ inline double GS_dhods( const double& m0, const double& k0 ) const;
	+ inline double GS_d( const double& m0, const double& k0 ) const;

	EvtComplex numerator( const EvtDalitzPoint& p, const EvtTwoBodyKine& vb,
	- const EvtTwoBodyKine& vd );
	- double angDep( const EvtDalitzPoint& p );
	- EvtComplex mixFactor( EvtComplex prop, EvtComplex prop_mix );
	- EvtComplex Fvector( double s, int index );
	- EvtComplex lass( double s );
	- EvtComplex flatte( const double& m );
	-
	- inline EvtComplex sqrtCplx( double in )
	+ const EvtTwoBodyKine& vd ) const;
	+ double angDep( const EvtDalitzPoint& p ) const;
	+ EvtComplex mixFactor( EvtComplex prop, EvtComplex prop_mix ) const;
	+ EvtComplex Fvector( double s, int index ) const;
	+ EvtComplex lass( double s ) const;
	+ EvtComplex flatte( const double& m ) const;
	+
	+ inline EvtComplex sqrtCplx( double in ) const
	{
	return ( in > 0 ) ? EvtComplex( sqrt( in ), 0 )
	: EvtComplex( 0, sqrt( -in ) );
	}

	// Dalitz plot
	- EvtDalitzPlot _dp;
	+ EvtDalitzPlot m_dp;

	// Pairing indices:
	- EvtCyclic3::Pair _pairAng; // angular
	- EvtCyclic3::Pair _pairRes; // resonance
	+ EvtCyclic3::Pair m_pairAng; // angular
	+ EvtCyclic3::Pair m_pairRes; // resonance

	// Spin
	- EvtSpinType::spintype _spin;
	+ EvtSpinType::spintype m_spin;

	// Numerator type
	- NumType _typeN;
	+ NumType m_typeN;

	// Nominal mass and width
	- double _m0, _g0;
	+ double m_m0, m_g0;

	// Vertices
	- EvtTwoBodyVertex _vb;
	- EvtTwoBodyVertex _vd;
	+ EvtTwoBodyVertex m_vb;
	+ EvtTwoBodyVertex m_vd;

	// Daughter masses
	- double _massFirst, _massSecond;
	+ double m_massFirst, m_massSecond;

	// variables for electromagnetic mass mixing
	- double _m0_mix, _g0_mix, _delta_mix;
	- EvtComplex _amp_mix;
	+ double m_m0_mix, m_g0_mix, m_delta_mix;
	+ EvtComplex m_amp_mix;

	// variables for coupled Breit-Wigner
	- double _g1, _g2;
	- CouplingType _coupling2;
	+ double m_g1, m_g2;
	+ CouplingType m_coupling2;

	// variables for Blatt-Weisskopf form factors
	- double _f_b, _f_d;
	+ double m_f_b, m_f_d;

	// K-matrix
	- int _kmatrix_index;
	- EvtComplex _fr12prod, _fr13prod, _fr14prod, _fr15prod;
	- double _s0prod;
	+ int m_kmatrix_index;
	+ EvtComplex m_fr12prod, m_fr13prod, m_fr14prod, m_fr15prod;
	+ double m_s0prod;

	// LASS
	- double _a;
	- double _r;
	- double _Blass;
	- double _phiB;
	- double _R;
	- double _phiR;
	- double _cutoff;
	- bool _scaleByMOverQ;
	+ double m_a;
	+ double m_r;
	+ double m_Blass;
	+ double m_phiB;
	+ double m_R;
	+ double m_phiR;
	+ double m_cutoff;
	+ bool m_scaleByMOverQ;

	//Nonresonant
	- double _alpha;
	+ double m_alpha;

	// Flatte
	- std::vector<EvtFlatteParam> _flatteParams;
	+ std::vector<EvtFlatteParam> m_flatteParams;
	};

	#endif
	diff --git a/EvtGenModels/EvtD0ToKspipi.hh b/EvtGenModels/EvtD0ToKspipi.hh
	new file mode 100644
	--- /dev/null
	+++ b/EvtGenModels/EvtD0ToKspipi.hh
	@@ -0,0 +1,71 @@
	+#ifndef EVTD0TOKSPIPI_HH
	+#define EVTD0TOKSPIPI_HH
	+
	+#include "EvtGenBase/EvtComplex.hh"
	+#include "EvtGenBase/EvtDalitzPoint.hh"
	+#include "EvtGenBase/EvtDalitzReso.hh"
	+#include "EvtGenBase/EvtDecayAmp.hh"
	+
	+#include <string>
	+#include <utility>
	+#include <vector>
	+
	+class EvtParticle;
	+
	+class EvtD0ToKspipi : public EvtDecayAmp {
	+ public:
	+ std::string getName() override;
	+ EvtDecayBase* clone() override;
	+
	+ void init() override;
	+ void initProbMax() override;
	+ void decay( EvtParticle* parent ) override;
	+
	+ private:
	+ // Calculate the total amplitude given the Dalitz plot point
	+ EvtComplex calcTotAmp( const EvtDalitzPoint& point ) const;
	+
	+ // Set particle IDs and PDG masses
	+ void setPDGValues();
	+
	+ // Setup the Dalitz plot resonances and their amplitude coefficients
	+ void initResonances();
	+
	+ // Daughter IDs (updated according to decay file ordering)
	+ int m_d0 = 0;
	+ int m_d1 = 1;
	+ int m_d2 = 2;
	+
	+ // Resonance lineshape and complex amplitude coefficient pair
	+ typedef std::pair<EvtDalitzReso, EvtComplex> ResAmpPair;
	+
	+ // Vector of (resonance, coeff) pairs
	+ std::vector<ResAmpPair> m_resonances;
	+
	+ // IDs of the relevant particles
	+ EvtId m_BP;
	+ EvtId m_BM;
	+ EvtId m_B0;
	+ EvtId m_B0B;
	+ EvtId m_D0;
	+ EvtId m_D0B;
	+ EvtId m_KM;
	+ EvtId m_KP;
	+ EvtId m_K0;
	+ EvtId m_K0B;
	+ EvtId m_KL;
	+ EvtId m_KS;
	+ EvtId m_PIM;
	+ EvtId m_PIP;
	+
	+ // Flavor of the B parent and D
	+ EvtId m_bFlavor;
	+ EvtId m_dFlavor;
	+
	+ // Masses of the relevant particles
	+ double m_mD0;
	+ double m_mKs;
	+ double m_mPi;
	+ double m_mK;
	+};
	+#endif
	diff --git a/History.md b/History.md
	--- a/History.md
	+++ b/History.md
	@@ -8,6 +8,9 @@

	https://phab.hepforge.org/Dxyz

	+17 Jun 2024 John Back
	+* D119: Add EvtD0ToKspipi DP model, courtesy of Camille Normand (LHCb).
	+
	===
	## R02-02-01

	diff --git a/src/EvtGenBase/EvtDalitzReso.cpp b/src/EvtGenBase/EvtDalitzReso.cpp
	--- a/src/EvtGenBase/EvtDalitzReso.cpp
	+++ b/src/EvtGenBase/EvtDalitzReso.cpp
	@@ -43,47 +43,48 @@
	EvtDalitzReso::EvtDalitzReso( const EvtDalitzPlot& dp, Pair pairAng, Pair pairRes,
	EvtSpinType::spintype spin, double m0, double g0,
	NumType typeN, double f_b, double f_d ) :
	- _dp( dp ),
	- _pairAng( pairAng ),
	- _pairRes( pairRes ),
	- _spin( spin ),
	- _typeN( typeN ),
	- _m0( m0 ),
	- _g0( g0 ),
	- _massFirst( dp.m( first( pairRes ) ) ),
	- _massSecond( dp.m( second( pairRes ) ) ),
	- _m0_mix( -1. ),
	- _g0_mix( 0. ),
	- _delta_mix( 0. ),
	- _amp_mix( 0., 0. ),
	- _g1( -1. ),
	- _g2( -1. ),
	- _coupling2( Undefined ),
	- _f_b( f_b ),
	- _f_d( f_d ),
	- _kmatrix_index( -1 ),
	- _fr12prod( 0., 0. ),
	- _fr13prod( 0., 0. ),
	- _fr14prod( 0., 0. ),
	- _fr15prod( 0., 0. ),
	- _s0prod( 0. ),
	- _a( 0. ),
	- _r( 0. ),
	- _Blass( 0. ),
	- _phiB( 0. ),
	- _R( 0. ),
	- _phiR( 0. ),
	- _cutoff( -1. ),
	- _scaleByMOverQ( false ),
	- _alpha( 0. )
	+ m_dp( dp ),
	+ m_pairAng( pairAng ),
	+ m_pairRes( pairRes ),
	+ m_spin( spin ),
	+ m_typeN( typeN ),
	+ m_m0( m0 ),
	+ m_g0( g0 ),
	+ m_massFirst( dp.m( first( pairRes ) ) ),
	+ m_massSecond( dp.m( second( pairRes ) ) ),
	+ m_m0_mix( -1. ),
	+ m_g0_mix( 0. ),
	+ m_delta_mix( 0. ),
	+ m_amp_mix( 0., 0. ),
	+ m_g1( -1. ),
	+ m_g2( -1. ),
	+ m_coupling2( Undefined ),
	+ m_f_b( f_b ),
	+ m_f_d( f_d ),
	+ m_kmatrix_index( -1 ),
	+ m_fr12prod( 0., 0. ),
	+ m_fr13prod( 0., 0. ),
	+ m_fr14prod( 0., 0. ),
	+ m_fr15prod( 0., 0. ),
	+ m_s0prod( 0. ),
	+ m_a( 0. ),
	+ m_r( 0. ),
	+ m_Blass( 0. ),
	+ m_phiB( 0. ),
	+ m_R( 0. ),
	+ m_phiR( 0. ),
	+ m_cutoff( -1. ),
	+ m_scaleByMOverQ( false ),
	+ m_alpha( 0. )
	{
	- _vb = EvtTwoBodyVertex( _m0, _dp.m( EvtCyclic3::other( _pairRes ) ),
	- _dp.bigM(), _spin );
	- _vd = EvtTwoBodyVertex( _massFirst, _massSecond, _m0, _spin );
	- _vb.set_f( _f_b ); // Default values for Blatt-Weisskopf factors are 0.0 and 1.5.
	- _vd.set_f( _f_d );
	- assert( _typeN != K_MATRIX && _typeN != K_MATRIX_I &&
	- _typeN != K_MATRIX_II ); // single BW cannot be K-matrix
	+ m_vb = EvtTwoBodyVertex( m_m0, m_dp.m( EvtCyclic3::other( m_pairRes ) ),
	+ m_dp.bigM(), m_spin );
	+ m_vd = EvtTwoBodyVertex( m_massFirst, m_massSecond, m_m0, m_spin );
	+ m_vb.set_f(
	+ m_f_b ); // Default values for Blatt-Weisskopf factors are 0.0 and 1.5.
	+ m_vd.set_f( m_f_d );
	+ assert( m_typeN != K_MATRIX && m_typeN != K_MATRIX_I &&
	+ m_typeN != K_MATRIX_II ); // single BW cannot be K-matrix
	}

	// Breit-Wigner with electromagnetic mass mixing
	@@ -91,47 +92,48 @@
	EvtSpinType::spintype spin, double m0, double g0,
	NumType typeN, double m0_mix, double g0_mix,
	double delta_mix, EvtComplex amp_mix ) :
	- _dp( dp ),
	- _pairAng( pairAng ),
	- _pairRes( pairRes ),
	- _spin( spin ),
	- _typeN( typeN ),
	- _m0( m0 ),
	- _g0( g0 ),
	- _massFirst( dp.m( first( pairRes ) ) ),
	- _massSecond( dp.m( second( pairRes ) ) ),
	- _m0_mix( m0_mix ),
	- _g0_mix( g0_mix ),
	- _delta_mix( delta_mix ),
	- _amp_mix( amp_mix ),
	- _g1( -1. ),
	- _g2( -1. ),
	- _coupling2( Undefined ),
	- _f_b( 0.0 ),
	- _f_d( 1.5 ),
	- _kmatrix_index( -1 ),
	- _fr12prod( 0., 0. ),
	- _fr13prod( 0., 0. ),
	- _fr14prod( 0., 0. ),
	- _fr15prod( 0., 0. ),
	- _s0prod( 0. ),
	- _a( 0. ),
	- _r( 0. ),
	- _Blass( 0. ),
	- _phiB( 0. ),
	- _R( 0. ),
	- _phiR( 0. ),
	- _cutoff( -1. ),
	- _scaleByMOverQ( false ),
	- _alpha( 0. )
	+ m_dp( dp ),
	+ m_pairAng( pairAng ),
	+ m_pairRes( pairRes ),
	+ m_spin( spin ),
	+ m_typeN( typeN ),
	+ m_m0( m0 ),
	+ m_g0( g0 ),
	+ m_massFirst( dp.m( first( pairRes ) ) ),
	+ m_massSecond( dp.m( second( pairRes ) ) ),
	+ m_m0_mix( m0_mix ),
	+ m_g0_mix( g0_mix ),
	+ m_delta_mix( delta_mix ),
	+ m_amp_mix( amp_mix ),
	+ m_g1( -1. ),
	+ m_g2( -1. ),
	+ m_coupling2( Undefined ),
	+ m_f_b( 0.0 ),
	+ m_f_d( 1.5 ),
	+ m_kmatrix_index( -1 ),
	+ m_fr12prod( 0., 0. ),
	+ m_fr13prod( 0., 0. ),
	+ m_fr14prod( 0., 0. ),
	+ m_fr15prod( 0., 0. ),
	+ m_s0prod( 0. ),
	+ m_a( 0. ),
	+ m_r( 0. ),
	+ m_Blass( 0. ),
	+ m_phiB( 0. ),
	+ m_R( 0. ),
	+ m_phiR( 0. ),
	+ m_cutoff( -1. ),
	+ m_scaleByMOverQ( false ),
	+ m_alpha( 0. )
	{
	- _vb = EvtTwoBodyVertex( _m0, _dp.m( EvtCyclic3::other( _pairRes ) ),
	- _dp.bigM(), _spin );
	- _vd = EvtTwoBodyVertex( _massFirst, _massSecond, _m0, _spin );
	- _vb.set_f( 0.0 ); // Default values for Blatt-Weisskopf factors.
	- _vd.set_f( 1.5 );
	+ m_vb = EvtTwoBodyVertex( m_m0, m_dp.m( EvtCyclic3::other( m_pairRes ) ),
	+ m_dp.bigM(), m_spin );
	+ m_vd = EvtTwoBodyVertex( m_massFirst, m_massSecond, m_m0, m_spin );
	+ m_vb.set_f( 0.0 ); // Default values for Blatt-Weisskopf factors.
	+ m_vd.set_f( 1.5 );
	// single BW (with electromagnetic mixing) cannot be K-matrix
	- assert( _typeN != K_MATRIX && _typeN != K_MATRIX_I && _typeN != K_MATRIX_II );
	+ assert( m_typeN != K_MATRIX && m_typeN != K_MATRIX_I &&
	+ m_typeN != K_MATRIX_II );
	}

	// coupled Breit-Wigner
	@@ -139,50 +141,50 @@
	Pair pairRes, EvtSpinType::spintype spin,
	double m0, NumType typeN, double g1, double g2,
	CouplingType coupling2 ) :
	- _dp( dp ),
	- _pairAng( pairAng ),
	- _pairRes( pairRes ),
	- _spin( spin ),
	- _typeN( typeN ),
	- _m0( m0 ),
	- _g0( -1. ),
	- _massFirst( dp.m( first( pairRes ) ) ),
	- _massSecond( dp.m( second( pairRes ) ) ),
	- _m0_mix( -1. ),
	- _g0_mix( 0. ),
	- _delta_mix( 0. ),
	- _amp_mix( 0., 0. ),
	- _g1( g1 ),
	- _g2( g2 ),
	- _coupling2( coupling2 ),
	- _f_b( 0.0 ),
	- _f_d( 1.5 ),
	- _kmatrix_index( -1 ),
	- _fr12prod( 0., 0. ),
	- _fr13prod( 0., 0. ),
	- _fr14prod( 0., 0. ),
	- _fr15prod( 0., 0. ),
	- _s0prod( 0. ),
	- _a( 0. ),
	- _r( 0. ),
	- _Blass( 0. ),
	- _phiB( 0. ),
	- _R( 0. ),
	- _phiR( 0. ),
	- _cutoff( -1. ),
	- _scaleByMOverQ( false ),
	- _alpha( 0. )
	+ m_dp( dp ),
	+ m_pairAng( pairAng ),
	+ m_pairRes( pairRes ),
	+ m_spin( spin ),
	+ m_typeN( typeN ),
	+ m_m0( m0 ),
	+ m_g0( -1. ),
	+ m_massFirst( dp.m( first( pairRes ) ) ),
	+ m_massSecond( dp.m( second( pairRes ) ) ),
	+ m_m0_mix( -1. ),
	+ m_g0_mix( 0. ),
	+ m_delta_mix( 0. ),
	+ m_amp_mix( 0., 0. ),
	+ m_g1( g1 ),
	+ m_g2( g2 ),
	+ m_coupling2( coupling2 ),
	+ m_f_b( 0.0 ),
	+ m_f_d( 1.5 ),
	+ m_kmatrix_index( -1 ),
	+ m_fr12prod( 0., 0. ),
	+ m_fr13prod( 0., 0. ),
	+ m_fr14prod( 0., 0. ),
	+ m_fr15prod( 0., 0. ),
	+ m_s0prod( 0. ),
	+ m_a( 0. ),
	+ m_r( 0. ),
	+ m_Blass( 0. ),
	+ m_phiB( 0. ),
	+ m_R( 0. ),
	+ m_phiR( 0. ),
	+ m_cutoff( -1. ),
	+ m_scaleByMOverQ( false ),
	+ m_alpha( 0. )
	{
	- _vb = EvtTwoBodyVertex( _m0, _dp.m( EvtCyclic3::other( _pairRes ) ),
	- _dp.bigM(), _spin );
	- _vd = EvtTwoBodyVertex( _massFirst, _massSecond, _m0, _spin );
	- _vb.set_f( 0.0 ); // Default values for Blatt-Weisskopf factors.
	- _vd.set_f( 1.5 );
	- assert( _coupling2 != Undefined );
	- assert( _typeN != K_MATRIX && _typeN != K_MATRIX_I &&
	- _typeN != K_MATRIX_II ); // coupled BW cannot be K-matrix
	- assert( _typeN != LASS ); // coupled BW cannot be LASS
	- assert( _typeN != NBW ); // for coupled BW, only relativistic BW
	+ m_vb = EvtTwoBodyVertex( m_m0, m_dp.m( EvtCyclic3::other( m_pairRes ) ),
	+ m_dp.bigM(), m_spin );
	+ m_vd = EvtTwoBodyVertex( m_massFirst, m_massSecond, m_m0, m_spin );
	+ m_vb.set_f( 0.0 ); // Default values for Blatt-Weisskopf factors.
	+ m_vd.set_f( 1.5 );
	+ assert( m_coupling2 != Undefined );
	+ assert( m_typeN != K_MATRIX && m_typeN != K_MATRIX_I &&
	+ m_typeN != K_MATRIX_II ); // coupled BW cannot be K-matrix
	+ assert( m_typeN != LASS ); // coupled BW cannot be LASS
	+ assert( m_typeN != NBW ); // for coupled BW, only relativistic BW
	}

	// K-Matrix (A&S)
	@@ -191,52 +193,53 @@
	EvtComplex fr12prod, EvtComplex fr13prod,
	EvtComplex fr14prod, EvtComplex fr15prod,
	double s0prod ) :
	- _dp( dp ),
	- _pairRes( pairRes ),
	- _typeN( typeN ),
	- _m0( 0. ),
	- _g0( 0. ),
	- _massFirst( dp.m( first( pairRes ) ) ),
	- _massSecond( dp.m( second( pairRes ) ) ),
	- _m0_mix( -1. ),
	- _g0_mix( 0. ),
	- _delta_mix( 0. ),
	- _amp_mix( 0., 0. ),
	- _g1( -1. ),
	- _g2( -1. ),
	- _coupling2( Undefined ),
	- _f_b( 0. ),
	- _f_d( 0. ),
	- _kmatrix_index( -1 ),
	- _fr12prod( fr12prod ),
	- _fr13prod( fr13prod ),
	- _fr14prod( fr14prod ),
	- _fr15prod( fr15prod ),
	- _s0prod( s0prod ),
	- _a( 0. ),
	- _r( 0. ),
	- _Blass( 0. ),
	- _phiB( 0. ),
	- _R( 0. ),
	- _phiR( 0. ),
	- _cutoff( -1. ),
	- _scaleByMOverQ( false ),
	- _alpha( 0. )
	+ m_dp( dp ),
	+ m_pairRes( pairRes ),
	+ m_typeN( typeN ),
	+ m_m0( 0. ),
	+ m_g0( 0. ),
	+ m_massFirst( dp.m( first( pairRes ) ) ),
	+ m_massSecond( dp.m( second( pairRes ) ) ),
	+ m_m0_mix( -1. ),
	+ m_g0_mix( 0. ),
	+ m_delta_mix( 0. ),
	+ m_amp_mix( 0., 0. ),
	+ m_g1( -1. ),
	+ m_g2( -1. ),
	+ m_coupling2( Undefined ),
	+ m_f_b( 0. ),
	+ m_f_d( 0. ),
	+ m_kmatrix_index( -1 ),
	+ m_fr12prod( fr12prod ),
	+ m_fr13prod( fr13prod ),
	+ m_fr14prod( fr14prod ),
	+ m_fr15prod( fr15prod ),
	+ m_s0prod( s0prod ),
	+ m_a( 0. ),
	+ m_r( 0. ),
	+ m_Blass( 0. ),
	+ m_phiB( 0. ),
	+ m_R( 0. ),
	+ m_phiR( 0. ),
	+ m_cutoff( -1. ),
	+ m_scaleByMOverQ( false ),
	+ m_alpha( 0. )
	{
	- assert( _typeN == K_MATRIX \|\| _typeN == K_MATRIX_I \|\| _typeN == K_MATRIX_II );
	- _spin = EvtSpinType::SCALAR;
	+ assert( m_typeN == K_MATRIX \|\| m_typeN == K_MATRIX_I \|\|
	+ m_typeN == K_MATRIX_II );
	+ m_spin = EvtSpinType::SCALAR;
	if ( nameIndex == "Pole1" )
	- _kmatrix_index = 1;
	+ m_kmatrix_index = 1;
	else if ( nameIndex == "Pole2" )
	- _kmatrix_index = 2;
	+ m_kmatrix_index = 2;
	else if ( nameIndex == "Pole3" )
	- _kmatrix_index = 3;
	+ m_kmatrix_index = 3;
	else if ( nameIndex == "Pole4" )
	- _kmatrix_index = 4;
	+ m_kmatrix_index = 4;
	else if ( nameIndex == "Pole5" )
	- _kmatrix_index = 5;
	+ m_kmatrix_index = 5;
	else if ( nameIndex == "f11prod" )
	- _kmatrix_index = 6;
	+ m_kmatrix_index = 6;
	else
	assert( 0 );
	}
	@@ -246,197 +249,212 @@
	double g0, double a, double r, double B,
	double phiB, double R, double phiR, double cutoff,
	bool scaleByMOverQ ) :
	- _dp( dp ),
	- _pairRes( pairRes ),
	- _typeN( LASS ),
	- _m0( m0 ),
	- _g0( g0 ),
	- _massFirst( dp.m( first( pairRes ) ) ),
	- _massSecond( dp.m( second( pairRes ) ) ),
	- _m0_mix( -1. ),
	- _g0_mix( 0. ),
	- _delta_mix( 0. ),
	- _amp_mix( 0., 0. ),
	- _g1( -1. ),
	- _g2( -1. ),
	- _coupling2( Undefined ),
	- _f_b( 0.0 ),
	- _f_d( 1.5 ),
	- _kmatrix_index( -1 ),
	- _fr12prod( 0., 0. ),
	- _fr13prod( 0., 0. ),
	- _fr14prod( 0., 0. ),
	- _fr15prod( 0., 0. ),
	- _s0prod( 0. ),
	- _a( a ),
	- _r( r ),
	- _Blass( B ),
	- _phiB( phiB ),
	- _R( R ),
	- _phiR( phiR ),
	- _cutoff( cutoff ),
	- _scaleByMOverQ( scaleByMOverQ ),
	- _alpha( 0. )
	+ m_dp( dp ),
	+ m_pairRes( pairRes ),
	+ m_typeN( LASS ),
	+ m_m0( m0 ),
	+ m_g0( g0 ),
	+ m_massFirst( dp.m( first( pairRes ) ) ),
	+ m_massSecond( dp.m( second( pairRes ) ) ),
	+ m_m0_mix( -1. ),
	+ m_g0_mix( 0. ),
	+ m_delta_mix( 0. ),
	+ m_amp_mix( 0., 0. ),
	+ m_g1( -1. ),
	+ m_g2( -1. ),
	+ m_coupling2( Undefined ),
	+ m_f_b( 0.0 ),
	+ m_f_d( 1.5 ),
	+ m_kmatrix_index( -1 ),
	+ m_fr12prod( 0., 0. ),
	+ m_fr13prod( 0., 0. ),
	+ m_fr14prod( 0., 0. ),
	+ m_fr15prod( 0., 0. ),
	+ m_s0prod( 0. ),
	+ m_a( a ),
	+ m_r( r ),
	+ m_Blass( B ),
	+ m_phiB( phiB ),
	+ m_R( R ),
	+ m_phiR( phiR ),
	+ m_cutoff( cutoff ),
	+ m_scaleByMOverQ( scaleByMOverQ ),
	+ m_alpha( 0. )
	{
	- _spin = EvtSpinType::SCALAR;
	- _vd = EvtTwoBodyVertex( _massFirst, _massSecond, _m0, _spin );
	- _vd.set_f( 1.5 ); // Default values for Blatt-Weisskopf factors.
	+ m_spin = EvtSpinType::SCALAR;
	+ m_vd = EvtTwoBodyVertex( m_massFirst, m_massSecond, m_m0, m_spin );
	+ m_vd.set_f( 1.5 ); // Default values for Blatt-Weisskopf factors.
	}

	//Flatte
	EvtDalitzReso::EvtDalitzReso( const EvtDalitzPlot& dp, EvtCyclic3::Pair pairRes,
	double m0 ) :
	- _dp( dp ),
	- _pairRes( pairRes ),
	- _typeN( FLATTE ),
	- _m0( m0 ),
	- _g0( 0. ),
	- _massFirst( dp.m( first( pairRes ) ) ),
	- _massSecond( dp.m( second( pairRes ) ) ),
	- _m0_mix( -1. ),
	- _g0_mix( 0. ),
	- _delta_mix( 0. ),
	- _amp_mix( 0., 0. ),
	- _g1( -1. ),
	- _g2( -1. ),
	- _coupling2( Undefined ),
	- _f_b( 0. ),
	- _f_d( 0. ),
	- _kmatrix_index( -1 ),
	- _fr12prod( 0., 0. ),
	- _fr13prod( 0., 0. ),
	- _fr14prod( 0., 0. ),
	- _fr15prod( 0., 0. ),
	- _s0prod( 0. ),
	- _a( 0. ),
	- _r( 0. ),
	- _Blass( 0. ),
	- _phiB( 0. ),
	- _R( 0. ),
	- _phiR( 0. ),
	- _cutoff( -1. ),
	- _scaleByMOverQ( false ),
	- _alpha( 0. )
	+ m_dp( dp ),
	+ m_pairRes( pairRes ),
	+ m_typeN( FLATTE ),
	+ m_m0( m0 ),
	+ m_g0( 0. ),
	+ m_massFirst( dp.m( first( pairRes ) ) ),
	+ m_massSecond( dp.m( second( pairRes ) ) ),
	+ m_m0_mix( -1. ),
	+ m_g0_mix( 0. ),
	+ m_delta_mix( 0. ),
	+ m_amp_mix( 0., 0. ),
	+ m_g1( -1. ),
	+ m_g2( -1. ),
	+ m_coupling2( Undefined ),
	+ m_f_b( 0. ),
	+ m_f_d( 0. ),
	+ m_kmatrix_index( -1 ),
	+ m_fr12prod( 0., 0. ),
	+ m_fr13prod( 0., 0. ),
	+ m_fr14prod( 0., 0. ),
	+ m_fr15prod( 0., 0. ),
	+ m_s0prod( 0. ),
	+ m_a( 0. ),
	+ m_r( 0. ),
	+ m_Blass( 0. ),
	+ m_phiB( 0. ),
	+ m_R( 0. ),
	+ m_phiR( 0. ),
	+ m_cutoff( -1. ),
	+ m_scaleByMOverQ( false ),
	+ m_alpha( 0. )
	{
	- _spin = EvtSpinType::SCALAR;
	+ m_spin = EvtSpinType::SCALAR;
	}

	-EvtComplex EvtDalitzReso::evaluate( const EvtDalitzPoint& x )
	+EvtComplex EvtDalitzReso::evaluate( const EvtDalitzPoint& x ) const
	{
	- double m = sqrt( x.q( _pairRes ) );
	+ double m = sqrt( x.q( m_pairRes ) );

	- if ( _typeN == NON_RES )
	+ if ( m_typeN == NON_RES )
	return EvtComplex( 1.0, 0.0 );

	- if ( _typeN == NON_RES_LIN )
	+ if ( m_typeN == NON_RES_LIN )
	return m * m;

	- if ( _typeN == NON_RES_EXP )
	- return exp( -_alpha * m * m );
	+ if ( m_typeN == NON_RES_EXP )
	+ return exp( -m_alpha * m * m );

	// do use always hash table (speed up fitting)
	- if ( _typeN == K_MATRIX \|\| _typeN == K_MATRIX_I \|\| _typeN == K_MATRIX_II )
	- return Fvector( m * m, _kmatrix_index );
	+ if ( m_typeN == K_MATRIX \|\| m_typeN == K_MATRIX_I \|\| m_typeN == K_MATRIX_II )
	+ return Fvector( m * m, m_kmatrix_index );

	- if ( _typeN == LASS )
	+ if ( m_typeN == LASS )
	return lass( m * m );

	- if ( _typeN == FLATTE )
	+ if ( m_typeN == FLATTE )
	return flatte( m );

	EvtComplex amp( 1.0, 0.0 );

	- if ( fabs( _dp.bigM() - x.bigM() ) > 0.000001 ) {
	- _vb = EvtTwoBodyVertex( _m0, _dp.m( EvtCyclic3::other( _pairRes ) ),
	- x.bigM(), _spin );
	- _vb.set_f( _f_b );
	+ if ( fabs( m_dp.bigM() - x.bigM() ) > 0.000001 ) {
	+ EvtGenReport( EVTGEN_WARNING, "EvtGen" )
	+ << "Warning in EvtDalitzReso::evaluate."
	+ << "The mass of the mother has changed from " << m_dp.bigM()
	+ << " to " << x.bigM() << ". " << std::endl;
	}
	- EvtTwoBodyKine vb( m, x.m( EvtCyclic3::other( _pairRes ) ), x.bigM() );
	- EvtTwoBodyKine vd( _massFirst, _massSecond, m );
	+
	+ EvtTwoBodyKine vb( m, x.m( EvtCyclic3::other( m_pairRes ) ), x.bigM() );
	+ EvtTwoBodyKine vd( m_massFirst, m_massSecond, m );

	EvtComplex prop( 0, 0 );
	- if ( _typeN == NBW ) {
	- prop = propBreitWigner( _m0, _g0, m );
	- } else if ( _typeN == GAUSS_CLEO \|\| _typeN == GAUSS_CLEO_ZEMACH ) {
	- prop = propGauss( _m0, _g0, m );
	+ if ( m_typeN == NBW ) {
	+ prop = propBreitWigner( m_m0, m_g0, m );
	+ } else if ( m_typeN == GAUSS_CLEO \|\| m_typeN == GAUSS_CLEO_ZEMACH ) {
	+ prop = propGauss( m_m0, m_g0, m );
	} else {
	- if ( _coupling2 == Undefined ) {
	+ if ( m_coupling2 == Undefined ) {
	// single BW
	- double g = ( _g0 <= 0. \|\| _vd.pD() <= 0. )
	- ? -_g0
	- : _g0 * _vd.widthFactor( vd ); // running width
	- if ( _typeN == GS_CLEO \|\| _typeN == GS_CLEO_ZEMACH ) {
	+ double g = ( m_g0 <= 0. \|\| m_vd.pD() <= 0. )
	+ ? -m_g0
	+ : m_g0 * m_vd.widthFactor( vd ); // running width
	+ if ( m_typeN == GS_CLEO \|\| m_typeN == GS_CLEO_ZEMACH ) {
	// Gounaris-Sakurai (GS)
	- prop = propGounarisSakurai( _m0, fabs( _g0 ), _vd.pD(), m, g,
	+ prop = propGounarisSakurai( m_m0, fabs( m_g0 ), m_vd.pD(), m, g,
	vd.p() );
	} else {
	// standard relativistic BW
	- prop = propBreitWignerRel( _m0, g, m );
	+ prop = propBreitWignerRel( m_m0, g, m );
	}
	} else {
	// coupled width BW
	EvtComplex G1, G2;
	- switch ( _coupling2 ) {
	+ switch ( m_coupling2 ) {
	case PicPic: {
	- G1 = _g1 * _g1 * psFactor( _massFirst, _massSecond, m );
	- static double mPic = EvtPDL::getMass( EvtPDL::getId( "pi+" ) );
	- G2 = _g2 * _g2 * psFactor( mPic, mPic, m );
	+ G1 = m_g1 * m_g1 * psFactor( m_massFirst, m_massSecond, m );
	+ static const double mPic = EvtPDL::getMass(
	+ EvtPDL::getId( "pi+" ) );
	+ G2 = m_g2 * m_g2 * psFactor( mPic, mPic, m );
	break;
	}
	case PizPiz: {
	- G1 = _g1 * _g1 * psFactor( _massFirst, _massSecond, m );
	- static double mPiz = EvtPDL::getMass( EvtPDL::getId( "pi0" ) );
	- G2 = _g2 * _g2 * psFactor( mPiz, mPiz, m );
	+ G1 = m_g1 * m_g1 * psFactor( m_massFirst, m_massSecond, m );
	+ static const double mPiz = EvtPDL::getMass(
	+ EvtPDL::getId( "pi0" ) );
	+ G2 = m_g2 * m_g2 * psFactor( mPiz, mPiz, m );
	break;
	}
	case PiPi: {
	- G1 = _g1 * _g1 * psFactor( _massFirst, _massSecond, m );
	- static double mPic = EvtPDL::getMass( EvtPDL::getId( "pi+" ) );
	- static double mPiz = EvtPDL::getMass( EvtPDL::getId( "pi0" ) );
	- G2 = _g2 * _g2 * psFactor( mPic, mPic, mPiz, mPiz, m );
	+ G1 = m_g1 * m_g1 * psFactor( m_massFirst, m_massSecond, m );
	+ static const double mPic = EvtPDL::getMass(
	+ EvtPDL::getId( "pi+" ) );
	+ static const double mPiz = EvtPDL::getMass(
	+ EvtPDL::getId( "pi0" ) );
	+ G2 = m_g2 * m_g2 * psFactor( mPic, mPic, mPiz, mPiz, m );
	break;
	}
	case KcKc: {
	- G1 = _g1 * _g1 * psFactor( _massFirst, _massSecond, m );
	- static double mKc = EvtPDL::getMass( EvtPDL::getId( "K+" ) );
	- G2 = _g2 * _g2 * psFactor( mKc, mKc, m );
	+ G1 = m_g1 * m_g1 * psFactor( m_massFirst, m_massSecond, m );
	+ static const double mKc = EvtPDL::getMass(
	+ EvtPDL::getId( "K+" ) );
	+ G2 = m_g2 * m_g2 * psFactor( mKc, mKc, m );
	break;
	}
	case KzKz: {
	- G1 = _g1 * _g1 * psFactor( _massFirst, _massSecond, m );
	- static double mKz = EvtPDL::getMass( EvtPDL::getId( "K0" ) );
	- G2 = _g2 * _g2 * psFactor( mKz, mKz, m );
	+ G1 = m_g1 * m_g1 * psFactor( m_massFirst, m_massSecond, m );
	+ static const double mKz = EvtPDL::getMass(
	+ EvtPDL::getId( "K0" ) );
	+ G2 = m_g2 * m_g2 * psFactor( mKz, mKz, m );
	break;
	}
	case KK: {
	- G1 = _g1 * _g1 * psFactor( _massFirst, _massSecond, m );
	- static double mKc = EvtPDL::getMass( EvtPDL::getId( "K+" ) );
	- static double mKz = EvtPDL::getMass( EvtPDL::getId( "K0" ) );
	- G2 = _g2 * _g2 * psFactor( mKc, mKc, mKz, mKz, m );
	+ G1 = m_g1 * m_g1 * psFactor( m_massFirst, m_massSecond, m );
	+ static const double mKc = EvtPDL::getMass(
	+ EvtPDL::getId( "K+" ) );
	+ static const double mKz = EvtPDL::getMass(
	+ EvtPDL::getId( "K0" ) );
	+ G2 = m_g2 * m_g2 * psFactor( mKc, mKc, mKz, mKz, m );
	break;
	}
	case EtaPic: {
	- G1 = _g1 * _g1 * psFactor( _massFirst, _massSecond, m );
	- static double mEta = EvtPDL::getMass( EvtPDL::getId( "eta" ) );
	- static double mPic = EvtPDL::getMass( EvtPDL::getId( "pi+" ) );
	- G2 = _g2 * _g2 * psFactor( mEta, mPic, m );
	+ G1 = m_g1 * m_g1 * psFactor( m_massFirst, m_massSecond, m );
	+ static const double mEta = EvtPDL::getMass(
	+ EvtPDL::getId( "eta" ) );
	+ static const double mPic = EvtPDL::getMass(
	+ EvtPDL::getId( "pi+" ) );
	+ G2 = m_g2 * m_g2 * psFactor( mEta, mPic, m );
	break;
	}
	case EtaPiz: {
	- G1 = _g1 * _g1 * psFactor( _massFirst, _massSecond, m );
	- static double mEta = EvtPDL::getMass( EvtPDL::getId( "eta" ) );
	- static double mPiz = EvtPDL::getMass( EvtPDL::getId( "pi0" ) );
	- G2 = _g2 * _g2 * psFactor( mEta, mPiz, m );
	+ G1 = m_g1 * m_g1 * psFactor( m_massFirst, m_massSecond, m );
	+ static const double mEta = EvtPDL::getMass(
	+ EvtPDL::getId( "eta" ) );
	+ static const double mPiz = EvtPDL::getMass(
	+ EvtPDL::getId( "pi0" ) );
	+ G2 = m_g2 * m_g2 * psFactor( mEta, mPiz, m );
	break;
	}
	case PicPicKK: {
	- static double mPic = EvtPDL::getMass( EvtPDL::getId( "pi+" ) );
	- //G1 = _g1_g1psFactor(mPic,mPic,m);
	- G1 = _g1 * psFactor( mPic, mPic, m );
	- static double mKc = EvtPDL::getMass( EvtPDL::getId( "K+" ) );
	- static double mKz = EvtPDL::getMass( EvtPDL::getId( "K0" ) );
	- //G2 = _g2_g2psFactor(mKc,mKc,mKz,mKz,m);
	- G2 = _g2 * psFactor( mKc, mKc, mKz, mKz, m );
	+ static const double mPic = EvtPDL::getMass(
	+ EvtPDL::getId( "pi+" ) );
	+ G1 = m_g1 * psFactor( mPic, mPic, m );
	+ static const double mKc = EvtPDL::getMass(
	+ EvtPDL::getId( "K+" ) );
	+ static const double mKz = EvtPDL::getMass(
	+ EvtPDL::getId( "K0" ) );
	+ G2 = m_g2 * psFactor( mKc, mKc, mKz, mKz, m );
	break;
	}
	default:
	@@ -447,28 +465,28 @@
	break;
	}
	// calculate standard couple BW propagator
	- if ( _coupling2 != WA76 )
	- prop = _g1 * propBreitWignerRelCoupled( _m0, G1, G2, m );
	+ if ( m_coupling2 != WA76 )
	+ prop = m_g1 * propBreitWignerRelCoupled( m_m0, G1, G2, m );
	}
	}
	amp *= prop;

	// Compute form-factors (Blatt-Weisskopf penetration factor)
	- amp *= _vb.formFactor( vb );
	- amp *= _vd.formFactor( vd );
	+ amp *= m_vb.formFactor( vb );
	+ amp *= m_vd.formFactor( vd );

	// Compute numerator (angular distribution)
	amp *= numerator( x, vb, vd );

	// Compute electromagnetic mass mixing factor
	- if ( _m0_mix > 0. ) {
	+ if ( m_m0_mix > 0. ) {
	EvtComplex prop_mix;
	- if ( _typeN == NBW ) {
	- prop_mix = propBreitWigner( _m0_mix, _g0_mix, m );
	+ if ( m_typeN == NBW ) {
	+ prop_mix = propBreitWigner( m_m0_mix, m_g0_mix, m );
	} else {
	- assert( _g1 < 0. ); // running width only
	- double g_mix = _g0_mix * _vd.widthFactor( vd );
	- prop_mix = propBreitWignerRel( _m0_mix, g_mix, m );
	+ assert( m_g1 < 0. ); // running width only
	+ double g_mix = m_g0_mix * m_vd.widthFactor( vd );
	+ prop_mix = propBreitWignerRel( m_m0_mix, g_mix, m );
	}
	amp *= mixFactor( prop, prop_mix );
	}
	@@ -476,7 +494,8 @@
	return amp;
	}

	-EvtComplex EvtDalitzReso::psFactor( double& ma, double& mb, double& m )
	+EvtComplex EvtDalitzReso::psFactor( const double& ma, const double& mb,
	+ const double& m ) const
	{
	if ( m > ( ma + mb ) ) {
	EvtTwoBodyKine vd( ma, mb, m );
	@@ -490,14 +509,15 @@
	}
	}

	-EvtComplex EvtDalitzReso::psFactor( double& ma1, double& mb1, double& ma2,
	- double& mb2, double& m )
	+EvtComplex EvtDalitzReso::psFactor( const double& ma1, const double& mb1,
	+ const double& ma2, const double& mb2,
	+ const double& m ) const
	{
	return 0.5 * ( psFactor( ma1, mb1, m ) + psFactor( ma2, mb2, m ) );
	}

	EvtComplex EvtDalitzReso::propGauss( const double& m0, const double& s0,
	- const double& m )
	+ const double& m ) const
	{
	// Gaussian
	double gauss = 1. / sqrt( EvtConst::twoPi ) / s0 *
	@@ -506,14 +526,14 @@
	}

	EvtComplex EvtDalitzReso::propBreitWigner( const double& m0, const double& g0,
	- const double& m )
	+ const double& m ) const
	{
	// non-relativistic BW
	return sqrt( g0 / EvtConst::twoPi ) / ( m - m0 - EvtComplex( 0.0, g0 / 2. ) );
	}

	-EvtComplex EvtDalitzReso::propBreitWignerRel( const double& m0,
	- const double& g0, const double& m )
	+EvtComplex EvtDalitzReso::propBreitWignerRel( const double& m0, const double& g0,
	+ const double& m ) const
	{
	// relativistic BW with real width
	return 1. / ( m0 * m0 - m * m - EvtComplex( 0., m0 * g0 ) );
	@@ -521,7 +541,7 @@

	EvtComplex EvtDalitzReso::propBreitWignerRel( const double& m0,
	const EvtComplex& g0,
	- const double& m )
	+ const double& m ) const
	{
	// relativistic BW with complex width
	return 1. / ( m0 * m0 - m * m - EvtComplex( 0., m0 ) * g0 );
	@@ -530,15 +550,16 @@
	EvtComplex EvtDalitzReso::propBreitWignerRelCoupled( const double& m0,
	const EvtComplex& g1,
	const EvtComplex& g2,
	- const double& m )
	+ const double& m ) const
	{
	// relativistic coupled BW
	return 1. / ( m0 * m0 - m * m - ( g1 + g2 ) );
	}

	-EvtComplex EvtDalitzReso::propGounarisSakurai( const double& m0, const double& g0,
	- const double& k0, const double& m,
	- const double& g, const double& k )
	+EvtComplex EvtDalitzReso::propGounarisSakurai( const double& m0,
	+ const double& g0, const double& k0,
	+ const double& m, const double& g,
	+ const double& k ) const
	{
	// Gounaris-Sakurai parameterization of pi+pi- P wave. PRD, Vol61, 112002. PRL, Vol21, 244.
	// Expressions taken from BAD637v4, after fixing the imaginary part of the BW denominator: i M_R Gamma_R(s) --> i sqrt(s) Gamma_R(s)
	@@ -549,7 +570,7 @@

	inline double EvtDalitzReso::GS_f( const double& m0, const double& g0,
	const double& k0, const double& m,
	- const double& k )
	+ const double& k ) const
	{
	// m: sqrt(s)
	// m0: nominal resonance mass
	@@ -560,66 +581,66 @@
	( m0 * m0 - m * m ) * k0 * k0 * GS_dhods( m0, k0 ) );
	}

	-inline double EvtDalitzReso::GS_h( const double& m, const double& k )
	+inline double EvtDalitzReso::GS_h( const double& m, const double& k ) const
	{
	return 2. / EvtConst::pi * k / m *
	- log( ( m + 2. * k ) / ( 2. * _massFirst ) );
	+ log( ( m + 2. * k ) / ( 2. * m_massFirst ) );
	}

	-inline double EvtDalitzReso::GS_dhods( const double& m0, const double& k0 )
	+inline double EvtDalitzReso::GS_dhods( const double& m0, const double& k0 ) const
	{
	return GS_h( m0, k0 ) * ( 0.125 / ( k0 * k0 ) - 0.5 / ( m0 * m0 ) ) +
	0.5 / ( EvtConst::pi * m0 * m0 );
	}

	-inline double EvtDalitzReso::GS_d( const double& m0, const double& k0 )
	+inline double EvtDalitzReso::GS_d( const double& m0, const double& k0 ) const
	{
	- return 3. / EvtConst::pi * _massFirst * _massFirst / ( k0 * k0 ) *
	- log( ( m0 + 2. * k0 ) / ( 2. * _massFirst ) ) +
	+ return 3. / EvtConst::pi * m_massFirst * m_massFirst / ( k0 * k0 ) *
	+ log( ( m0 + 2. * k0 ) / ( 2. * m_massFirst ) ) +
	m0 / ( 2. * EvtConst::pi * k0 ) -
	- _massFirst * _massFirst * m0 / ( EvtConst::pi * k0 * k0 * k0 );
	+ m_massFirst * m_massFirst * m0 / ( EvtConst::pi * k0 * k0 * k0 );
	}

	EvtComplex EvtDalitzReso::numerator( const EvtDalitzPoint& x,
	const EvtTwoBodyKine& vb,
	- const EvtTwoBodyKine& vd )
	+ const EvtTwoBodyKine& vd ) const
	{
	EvtComplex ret( 0., 0. );

	// Non-relativistic Breit-Wigner
	- if ( NBW == _typeN ) {
	+ if ( NBW == m_typeN ) {
	ret = angDep( x );
	}

	// Standard relativistic Zemach propagator
	- else if ( RBW_ZEMACH == _typeN ) {
	- ret = _vd.phaseSpaceFactor( vd, EvtTwoBodyKine::AB ) * angDep( x );
	+ else if ( RBW_ZEMACH == m_typeN ) {
	+ ret = m_vd.phaseSpaceFactor( vd, EvtTwoBodyKine::AB ) * angDep( x );
	}

	// Standard relativistic Zemach propagator
	- else if ( RBW_ZEMACH2 == _typeN ) {
	- ret = _vd.phaseSpaceFactor( vd, EvtTwoBodyKine::AB ) *
	- _vb.phaseSpaceFactor( vb, EvtTwoBodyKine::AB ) * angDep( x );
	- if ( _spin == EvtSpinType::VECTOR ) {
	+ else if ( RBW_ZEMACH2 == m_typeN ) {
	+ ret = m_vd.phaseSpaceFactor( vd, EvtTwoBodyKine::AB ) *
	+ m_vb.phaseSpaceFactor( vb, EvtTwoBodyKine::AB ) * angDep( x );
	+ if ( m_spin == EvtSpinType::VECTOR ) {
	ret *= -4.;
	- } else if ( _spin == EvtSpinType::TENSOR ) {
	+ } else if ( m_spin == EvtSpinType::TENSOR ) {
	ret *= 16. / 3.;
	- } else if ( _spin != EvtSpinType::SCALAR )
	+ } else if ( m_spin != EvtSpinType::SCALAR )
	assert( 0 );
	}

	// Kuehn-Santamaria normalization:
	- else if ( RBW_KUEHN == _typeN ) {
	- ret = _m0 * _m0 * angDep( x );
	+ else if ( RBW_KUEHN == m_typeN ) {
	+ ret = m_m0 * m_m0 * angDep( x );
	}

	// CLEO amplitude
	- else if ( ( RBW_CLEO == _typeN ) \|\| ( GS_CLEO == _typeN ) \|\|
	- ( RBW_CLEO_ZEMACH == _typeN ) \|\| ( GS_CLEO_ZEMACH == _typeN ) \|\|
	- ( GAUSS_CLEO == _typeN ) \|\| ( GAUSS_CLEO_ZEMACH == _typeN ) ) {
	- Index iA = other( _pairAng ); // A = other(BC)
	- Index iB = common( _pairRes, _pairAng ); // B = common(AB,BC)
	- Index iC = other( _pairRes ); // C = other(AB)
	+ else if ( ( RBW_CLEO == m_typeN ) \|\| ( GS_CLEO == m_typeN ) \|\|
	+ ( RBW_CLEO_ZEMACH == m_typeN ) \|\| ( GS_CLEO_ZEMACH == m_typeN ) \|\|
	+ ( GAUSS_CLEO == m_typeN ) \|\| ( GAUSS_CLEO_ZEMACH == m_typeN ) ) {
	+ Index iA = other( m_pairAng ); // A = other(BC)
	+ Index iB = common( m_pairRes, m_pairAng ); // B = common(AB,BC)
	+ Index iC = other( m_pairRes ); // C = other(AB)

	double M = x.bigM();
	double mA = x.m( iA );
	@@ -630,20 +651,20 @@
	double qCA = x.q( combine( iC, iA ) );

	double M2 = M * M;
	- double m02 = ( ( RBW_CLEO_ZEMACH == _typeN ) \|\|
	- ( GS_CLEO_ZEMACH == _typeN ) \|\|
	- ( GAUSS_CLEO_ZEMACH == _typeN ) )
	+ double m02 = ( ( RBW_CLEO_ZEMACH == m_typeN ) \|\|
	+ ( GS_CLEO_ZEMACH == m_typeN ) \|\|
	+ ( GAUSS_CLEO_ZEMACH == m_typeN ) )
	? qAB
	- : _m0 * _m0;
	+ : m_m0 * m_m0;
	double mA2 = mA * mA;
	double mB2 = mB * mB;
	double mC2 = mC * mC;

	- if ( _spin == EvtSpinType::SCALAR )
	+ if ( m_spin == EvtSpinType::SCALAR )
	ret = EvtComplex( 1., 0. );
	- else if ( _spin == EvtSpinType::VECTOR ) {
	+ else if ( m_spin == EvtSpinType::VECTOR ) {
	ret = qCA - qBC + ( M2 - mC2 ) * ( mB2 - mA2 ) / m02;
	- } else if ( _spin == EvtSpinType::TENSOR ) {
	+ } else if ( m_spin == EvtSpinType::TENSOR ) {
	double x1 = qBC - qCA + ( M2 - mC2 ) * ( mA2 - mB2 ) / m02;
	double x2 = M2 - mC2;
	double x3 = qAB - 2 * M2 - 2 * mC2 + x2 * x2 / m02;
	@@ -657,30 +678,30 @@
	return ret;
	}

	-double EvtDalitzReso::angDep( const EvtDalitzPoint& x )
	+double EvtDalitzReso::angDep( const EvtDalitzPoint& x ) const
	{
	// Angular dependece for factorizable amplitudes
	// unphysical cosines indicate we are in big trouble
	double cosTh = x.cosTh(
	- _pairAng, _pairRes ); // angle between common(reso,ang) and other(reso)
	+ m_pairAng, m_pairRes ); // angle between common(reso,ang) and other(reso)
	if ( fabs( cosTh ) > 1. ) {
	EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "cosTh " << cosTh << std::endl;
	assert( 0 );
	}

	// in units of half-spin
	- return EvtdFunction::d( EvtSpinType::getSpin2( _spin ), 2 * 0, 2 * 0,
	+ return EvtdFunction::d( EvtSpinType::getSpin2( m_spin ), 2 * 0, 2 * 0,
	acos( cosTh ) );
	}

	-EvtComplex EvtDalitzReso::mixFactor( EvtComplex prop, EvtComplex prop_mix )
	+EvtComplex EvtDalitzReso::mixFactor( EvtComplex prop, EvtComplex prop_mix ) const
	{
	- double Delta = _delta_mix * ( _m0 + _m0_mix );
	+ double Delta = m_delta_mix * ( m_m0 + m_m0_mix );
	return 1 / ( 1 - Delta * Delta * prop * prop_mix ) *
	- ( 1 + _amp_mix * Delta * prop_mix );
	+ ( 1 + m_amp_mix * Delta * prop_mix );
	}

	-EvtComplex EvtDalitzReso::Fvector( double s, int index )
	+EvtComplex EvtDalitzReso::Fvector( double s, int index ) const
	{
	assert( index >= 1 && index <= 6 );

	@@ -696,11 +717,11 @@
	double g[5][5]; // Coupling constants. The first index is the pole index. The second index is the decay channel
	double ma[5]; // Pole masses. The unit is in GeV

	- int solution = ( _typeN == K_MATRIX )
	+ int solution = ( m_typeN == K_MATRIX )
	? 3
	- : ( ( _typeN == K_MATRIX_I )
	+ : ( ( m_typeN == K_MATRIX_I )
	? 1
	- : ( ( _typeN == K_MATRIX_II ) ? 2 : 0 ) );
	+ : ( ( m_typeN == K_MATRIX_II ) ? 2 : 0 ) );
	if ( solution == 0 ) {
	std::cout << "EvtDalitzReso::Fvector() error. Kmatrix solution incorrectly chosen ! "
	<< std::endl;
	@@ -965,7 +986,7 @@

	//This is not correct!
	//(1-ipK) != (1-iKp)
	- static EvtMatrix<EvtComplex> mat;
	+ static thread_local EvtMatrix<EvtComplex> mat;
	mat.setRange(
	5 ); // Try to do in only the first time. DEFINE ALLOCATION IN CONSTRUCTOR.

	@@ -998,59 +1019,59 @@
	} else {
	//this calculates fprod Factors
	value += U1j[0];
	- value += U1j[1] * _fr12prod;
	- value += U1j[2] * _fr13prod;
	- value += U1j[3] * _fr14prod;
	- value += U1j[4] * _fr15prod;
	+ value += U1j[1] * m_fr12prod;
	+ value += U1j[2] * m_fr13prod;
	+ value += U1j[3] * m_fr14prod;
	+ value += U1j[4] * m_fr15prod;

	- value = ( 1 - _s0prod ) / ( s - _s0prod ) smallTerm;
	+ value = ( 1 - m_s0prod ) / ( s - m_s0prod ) smallTerm;
	}

	return value;
	}

	//replace Breit-Wigner with LASS
	-EvtComplex EvtDalitzReso::lass( double s )
	+EvtComplex EvtDalitzReso::lass( double s ) const
	{
	- EvtTwoBodyKine vd( _massFirst, _massSecond, sqrt( s ) );
	+ EvtTwoBodyKine vd( m_massFirst, m_massSecond, sqrt( s ) );
	double q = vd.p();
	- double GammaM = _g0 * _vd.widthFactor( vd ); // running width;
	+ double GammaM = m_g0 * m_vd.widthFactor( vd ); // running width;

	//calculate the background phase motion
	- double cot_deltaB = 1.0 / ( _a * q ) + 0.5 * _r * q;
	+ double cot_deltaB = 1.0 / ( m_a * q ) + 0.5 * m_r * q;
	double deltaB = atan( 1.0 / cot_deltaB );
	- double totalB = deltaB + _phiB;
	+ double totalB = deltaB + m_phiB;

	//calculate the resonant phase motion
	- double deltaR = atan( ( _m0 * GammaM / ( _m0 * _m0 - s ) ) );
	- double totalR = deltaR + _phiR;
	+ double deltaR = atan( ( m_m0 * GammaM / ( m_m0 * m_m0 - s ) ) );
	+ double totalR = deltaR + m_phiR;

	//sum them up
	EvtComplex bkgB, resT;
	- bkgB = EvtComplex( _Blass * sin( totalB ), 0 ) *
	+ bkgB = EvtComplex( m_Blass * sin( totalB ), 0 ) *
	EvtComplex( cos( totalB ), sin( totalB ) );
	- resT = EvtComplex( _R * sin( deltaR ), 0 ) *
	+ resT = EvtComplex( m_R * sin( deltaR ), 0 ) *
	EvtComplex( cos( totalR ), sin( totalR ) ) *
	EvtComplex( cos( 2 * totalB ), sin( 2 * totalB ) );

	EvtComplex T;
	- if ( _cutoff > 0 && sqrt( s ) > _cutoff )
	+ if ( m_cutoff > 0 && sqrt( s ) > m_cutoff )
	T = resT;
	else
	T = bkgB + resT;

	- if ( _scaleByMOverQ )
	+ if ( m_scaleByMOverQ )
	T *= ( sqrt( s ) / q );

	return T;
	}

	-EvtComplex EvtDalitzReso::flatte( const double& m )
	+EvtComplex EvtDalitzReso::flatte( const double& m ) const
	{
	EvtComplex w;

	- for ( vector<EvtFlatteParam>::const_iterator param = _flatteParams.begin();
	- param != _flatteParams.end(); ++param ) {
	+ for ( vector<EvtFlatteParam>::const_iterator param = m_flatteParams.begin();
	+ param != m_flatteParams.end(); ++param ) {
	double m1 = ( *param ).m1();
	double m2 = ( *param ).m2();
	double g = ( *param ).g();
	@@ -1059,7 +1080,7 @@
	( 1 - ( ( m1 + m2 ) * ( m1 + m2 ) ) / ( m * m ) ) ) );
	}

	- EvtComplex denom = _m0 * _m0 - m * m - EvtComplex( 0, 1 ) * w;
	+ EvtComplex denom = m_m0 * m_m0 - m * m - EvtComplex( 0, 1 ) * w;

	return EvtComplex( 1.0, 0.0 ) / denom;
	}
	diff --git a/src/EvtGenModels/EvtD0ToKspipi.cpp b/src/EvtGenModels/EvtD0ToKspipi.cpp
	new file mode 100644
	--- /dev/null
	+++ b/src/EvtGenModels/EvtD0ToKspipi.cpp
	@@ -0,0 +1,331 @@
	+#include "EvtGenModels/EvtD0ToKspipi.hh"
	+
	+#include "EvtGenBase/EvtDecayTable.hh"
	+#include "EvtGenBase/EvtPDL.hh"
	+#include "EvtGenBase/EvtParticle.hh"
	+#include "EvtGenBase/EvtPatches.hh"
	+#include "EvtGenBase/EvtReport.hh"
	+#include "EvtGenBase/EvtVector4R.hh"
	+
	+#include <iostream>
	+
	+std::string EvtD0ToKspipi::getName()
	+{
	+ return "D0TOKSPIPI";
	+}
	+
	+EvtDecayBase* EvtD0ToKspipi::clone()
	+{
	+ return new EvtD0ToKspipi;
	+}
	+
	+void EvtD0ToKspipi::init()
	+{
	+ // Check that there are 0 arguments
	+ checkNArg( 0 );
	+
	+ // Check parent and daughter types
	+ checkNDaug( 3 );
	+ checkSpinDaughter( 0, EvtSpinType::SCALAR );
	+ checkSpinDaughter( 1, EvtSpinType::SCALAR );
	+ checkSpinDaughter( 2, EvtSpinType::SCALAR );
	+ checkSpinParent( EvtSpinType::SCALAR );
	+
	+ // Set the particle IDs and PDG masses
	+ setPDGValues();
	+
	+ // Set the EvtId of the three D0 daughters
	+ const int nDaug = 3;
	+ std::vector<EvtId> dau;
	+ for ( int index = 0; index < nDaug; index++ ) {
	+ dau.push_back( getDaug( index ) );
	+ }
	+
	+ // Look for K0bar h+ h-. The order will be K[0SL] h+ h-
	+ for ( int index = 0; index < nDaug; index++ ) {
	+ if ( ( dau[index] == m_K0B ) \|\| ( dau[index] == m_KS ) \|\|
	+ ( dau[index] == m_KL ) ) {
	+ m_d0 = index;
	+ } else if ( dau[index] == m_PIP ) {
	+ m_d1 = index;
	+ } else if ( dau[index] == m_PIM ) {
	+ m_d2 = index;
	+ } else {
	+ EvtGenReport( EVTGEN_ERROR, "EvtD0ToKspipi" )
	+ << "Daughter " << index << " has wrong ID" << std::endl;
	+ }
	+ }
	+
	+ // Setup the Dalitz plot resonances and their amplitude coefficients
	+ initResonances();
	+}
	+
	+void EvtD0ToKspipi::initProbMax()
	+{
	+ setProbMax( 6000.0 );
	+}
	+
	+void EvtD0ToKspipi::decay( EvtParticle* p )
	+{
	+ // Phase space generation and 4-momenta
	+ p->initializePhaseSpace( getNDaug(), getDaugs() );
	+ const EvtVector4R p0 = p->getDaug( m_d0 )->getP4();
	+ const EvtVector4R p1 = p->getDaug( m_d1 )->getP4();
	+ const EvtVector4R p2 = p->getDaug( m_d2 )->getP4();
	+
	+ // Squared invariant masses
	+ const double mSq01 = ( p0 + p1 ).mass2();
	+ const double mSq02 = ( p0 + p2 ).mass2();
	+ const double mSq12 = ( p1 + p2 ).mass2();
	+
	+ // For the decay amplitude
	+ EvtComplex amp( 0.0, 0.0 );
	+
	+ // Flavor of the D meson
	+ m_dFlavor = p->getId();
	+
	+ // Direct and conjugated Dalitz points
	+ const EvtDalitzPoint pointD0( m_mKs, m_mPi, m_mPi, mSq02, mSq12, mSq01 );
	+ const EvtDalitzPoint pointD0b( m_mKs, m_mPi, m_mPi, mSq01, mSq12, mSq02 );
	+
	+ // Check if the D is from a B+- -> D0 K+- decay with the appropriate model
	+ EvtParticle* parent = p->getParent();
	+ if ( parent != nullptr &&
	+ EvtDecayTable::getInstance()->getDecayFunc( parent )->getName() ==
	+ "BTODDALITZCPK" ) {
	+ const EvtId parId = parent->getId();
	+ if ( ( parId == m_BP ) \|\| ( parId == m_BM ) \|\| ( parId == m_B0 ) \|\|
	+ ( parId == m_B0B ) ) {
	+ m_bFlavor = parId;
	+
	+ // D0 parent particle is a B meson from the BTODDALITZCPK decay model.
	+ // D0 decay amplitude combines the interference of D0 and D0bar.
	+ // Read the D decay parameters from the B decay model.
	+ // Gamma angle in radians
	+ const double gamma =
	+ EvtDecayTable::getInstance()->getDecayFunc( parent )->getArg( 0 );
	+ // Strong phase in radians
	+ const double delta =
	+ EvtDecayTable::getInstance()->getDecayFunc( parent )->getArg( 1 );
	+ // Ratio between B -> D0 K and B -> D0bar K
	+ const double rB =
	+ EvtDecayTable::getInstance()->getDecayFunc( parent )->getArg( 2 );
	+
	+ // Direct and conjugated amplitudes
	+ const EvtComplex ampD0 = calcTotAmp( pointD0 );
	+ const EvtComplex ampD0b = calcTotAmp( pointD0b );
	+
	+ if ( m_bFlavor == m_BP \|\| m_bFlavor == m_B0 ) {
	+ // B+ or B0
	+ const EvtComplex iPhase( 0.0, delta + gamma );
	+ const EvtComplex expo( exp( iPhase ) );
	+ amp = ampD0b + rB * expo * ampD0;
	+ } else {
	+ // B- or B0bar
	+ const EvtComplex iPhase( 0.0, delta - gamma );
	+ const EvtComplex expo( exp( iPhase ) );
	+ amp = ampD0 + rB * expo * ampD0b;
	+ }
	+ }
	+ } else if ( !parent ) {
	+ // D0 has no parent particle. Use direct or conjugated amplitude
	+ if ( m_dFlavor == m_D0 ) {
	+ amp = calcTotAmp( pointD0 );
	+ } else {
	+ amp = calcTotAmp( pointD0b );
	+ }
	+ }
	+
	+ // Set the decay vertex amplitude
	+ vertex( amp );
	+}
	+
	+EvtComplex EvtD0ToKspipi::calcTotAmp( const EvtDalitzPoint& point ) const
	+{
	+ // Initialise the total amplitude
	+ EvtComplex totAmp( 0.0, 0.0 );
	+ // Check that the Dalitz plot point is OK
	+ if ( point.isValid() == false ) {
	+ return totAmp;
	+ }
	+
	+ // Add the resonance amplitudes by iterating over the (resonance, coeff) pairs.
	+ // This includes the BW and LASS lineshapes, as well as the K-matrix contributions
	+ for ( auto& [res, amp] : m_resonances ) {
	+ // Evaluate the resonance amplitude and multiply by the coeff
	+ totAmp += res.evaluate( point ) * amp;
	+ }
	+ // Return the total amplitude
	+ return totAmp;
	+}
	+
	+void EvtD0ToKspipi::initResonances()
	+{
	+ // Dalitz plot model from combined BaBar and BELLE paper hep-ex/1804.06153
	+
	+ // Define the Dalitz plot axes
	+ const EvtCyclic3::Pair AB = EvtCyclic3::AB;
	+ const EvtCyclic3::Pair AC = EvtCyclic3::AC;
	+ const EvtCyclic3::Pair BC = EvtCyclic3::BC;
	+
	+ // Define the particle spin and lineshape types
	+ const EvtSpinType::spintype vector = EvtSpinType::VECTOR;
	+ const EvtSpinType::spintype tensor = EvtSpinType::TENSOR;
	+ const EvtDalitzReso::NumType RBW = EvtDalitzReso::RBW_CLEO_ZEMACH;
	+ const EvtDalitzReso::NumType KMAT = EvtDalitzReso::K_MATRIX;
	+
	+ // Define the Dalitz plot
	+ const EvtDalitzPlot DP( m_mKs, m_mPi, m_mPi, m_mD0, 0, 0 );
	+
	+ // Clear the internal vector of (resonance, coeff) pairs
	+ m_resonances.clear();
	+
	+ // rho BW
	+ const EvtDalitzReso rhoBW( DP, AB, BC, vector, 0.77155, 0.13469, RBW, 5.0,
	+ 1.5 );
	+ const EvtComplex rhoCoeff( 1.0, 0.0 );
	+ m_resonances.push_back( std::make_pair( rhoBW, rhoCoeff ) );
	+
	+ // Omega BW
	+ const EvtDalitzReso omegaBW( DP, AB, BC, vector, 0.78265, 0.00849, RBW, 5.0,
	+ 1.5 );
	+ const EvtComplex omegaCoeff( -0.019829903319132, 0.033339785741436 );
	+ m_resonances.push_back( std::make_pair( omegaBW, omegaCoeff ) );
	+
	+ // K*(892)- BW
	+ const EvtDalitzReso KstarBW( DP, BC, AB, vector, 0.893709298220334,
	+ 0.047193287094108, RBW, 5.0, 1.5 );
	+ const EvtComplex KstarCoeff( -1.255025021860793, 1.176780750003210 );
	+ m_resonances.push_back( std::make_pair( KstarBW, KstarCoeff ) );
	+
	+ // K*0(1430)- LASS
	+ const double LASS_F = 0.955319683174069;
	+ const double LASS_phi_F = 0.001737032480754;
	+ const double LASS_R = 1.0;
	+ const double LASS_phi_R = -1.914503836666840;
	+ const double LASS_a = 0.112673863011817;
	+ const double LASS_r = -33.799002116066454;
	+ const EvtDalitzReso Kstar0_1430LASS = EvtDalitzReso(
	+ DP, AB, 1.440549945739415, 0.192611512914605, LASS_a, LASS_r, LASS_F,
	+ LASS_phi_F, LASS_R, LASS_phi_R, false, -1.0 );
	+ const EvtComplex Kstar0_1430Coeff( -0.386469884688245, 2.330315087713914 );
	+ m_resonances.push_back( std::make_pair( Kstar0_1430LASS, Kstar0_1430Coeff ) );
	+
	+ // K*2(1430)- BW
	+ const EvtDalitzReso Kstar2_1430BW( DP, BC, AB, tensor, 1.4256, 0.0985, RBW,
	+ 5.0, 1.5 );
	+ const EvtComplex Kstar2_1430Coeff( 0.914470111251261, -0.885129049790117 );
	+ m_resonances.push_back( std::make_pair( Kstar2_1430BW, Kstar2_1430Coeff ) );
	+
	+ // K*(1680)- BW
	+ const EvtDalitzReso Kstar_1680BW( DP, BC, AB, vector, 1.717, 0.322, RBW,
	+ 5.0, 1.5 );
	+ const EvtComplex Kstar_1680Coeff( -1.560837188791231, -2.916210561577914 );
	+ m_resonances.push_back( std::make_pair( Kstar_1680BW, Kstar_1680Coeff ) );
	+
	+ // K*(1410)- BW
	+ const EvtDalitzReso Kstar_1410BW( DP, BC, AB, vector, 1.414, 0.232, RBW,
	+ 5.0, 1.5 );
	+ const EvtComplex Kstar_1410Coeff( -0.046795079734847, 0.283085379985959 );
	+ m_resonances.push_back( std::make_pair( Kstar_1410BW, Kstar_1410Coeff ) );
	+
	+ // K*(892)+ DCS BW
	+ const EvtDalitzReso Kstar_DCSBW( DP, BC, AC, vector, 0.893709298220334,
	+ 0.047193287094108, RBW, 5.0, 1.5 );
	+ const EvtComplex Kstar_DCSCoeff( 0.121693743404499, -0.110206354657867 );
	+ m_resonances.push_back( std::make_pair( Kstar_DCSBW, Kstar_DCSCoeff ) );
	+
	+ // K*0(1430)+ DCS LASS
	+ const EvtDalitzReso Kstar0_1430_DCSLASS = EvtDalitzReso(
	+ DP, AC, 1.440549945739415, 0.192611512914605, LASS_a, LASS_r, LASS_F,
	+ LASS_phi_F, LASS_R, LASS_phi_R, false, -1.0 );
	+ const EvtComplex Kstar0_1430_DCSCoeff( -0.101484805664368, 0.032368302993344 );
	+ m_resonances.push_back(
	+ std::make_pair( Kstar0_1430_DCSLASS, Kstar0_1430_DCSCoeff ) );
	+
	+ // K*2(1430)+ DCS BW
	+ const EvtDalitzReso Kstar2_1430_DCSBW( DP, AB, AC, tensor, 1.4256, 0.0985,
	+ RBW, 5.0, 1.5 );
	+ const EvtComplex Kstar2_1430_DCSCoeff( 0.000699701539252, -0.102571188336701 );
	+ m_resonances.push_back(
	+ std::make_pair( Kstar2_1430_DCSBW, Kstar2_1430_DCSCoeff ) );
	+
	+ // K*(1410)+ DCS BW
	+ const EvtDalitzReso Kstar_1410_DCSBW( DP, BC, AC, vector, 1.414, 0.232, RBW,
	+ 5.0, 1.5 );
	+ const EvtComplex Kstar_1410_DCSCoeff( -0.181330401419455, 0.103990039950039 );
	+ m_resonances.push_back(
	+ std::make_pair( Kstar_1410_DCSBW, Kstar_1410_DCSCoeff ) );
	+
	+ // f2(1270) BW
	+ const EvtDalitzReso f2_1270BW( DP, AB, BC, tensor, 1.2751, 0.1842, RBW, 5.0,
	+ 1.5 );
	+ const EvtComplex f2_1270Coeff( 1.151785277682948, -0.845612891825272 );
	+ m_resonances.push_back( std::make_pair( f2_1270BW, f2_1270Coeff ) );
	+
	+ // rho(1450) BW
	+ const EvtDalitzReso rho_1450BW( DP, AB, BC, vector, 1.465, 0.400, RBW, 5.0,
	+ 1.5 );
	+ const EvtComplex rho_1450Coeff( -0.597963342540235, 2.787903868470057 );
	+ m_resonances.push_back( std::make_pair( rho_1450BW, rho_1450Coeff ) );
	+
	+ // K-matrix pole 1
	+ const double sProd0 = -0.07;
	+ const EvtDalitzReso pole1( DP, BC, "Pole1", KMAT, 0, 0, 0, 0, sProd0 );
	+ const EvtComplex p1Coeff( 3.122415682166643, 7.928823290976309 );
	+ m_resonances.push_back( std::make_pair( pole1, p1Coeff ) );
	+
	+ // K-matrix pole 2
	+ const EvtDalitzReso pole2( DP, BC, "Pole2", KMAT, 0, 0, 0, 0, sProd0 );
	+ const EvtComplex p2Coeff( 11.139907856904129, 4.948420661321371 );
	+ m_resonances.push_back( std::make_pair( pole2, p2Coeff ) );
	+
	+ // K-matrix pole 3
	+ const EvtDalitzReso pole3( DP, BC, "Pole3", KMAT, 0, 0, 0, 0, sProd0 );
	+ const EvtComplex p3Coeff( 29.146102368470210, -0.053588781806890 );
	+ m_resonances.push_back( std::make_pair( pole3, p3Coeff ) );
	+
	+ // K-matrix pole 4
	+ const EvtDalitzReso pole4( DP, BC, "Pole4", KMAT, 0, 0, 0, 0, sProd0 );
	+ const EvtComplex p4Coeff( 6.631556203215280, -8.455370251307063 );
	+ m_resonances.push_back( std::make_pair( pole4, p4Coeff ) );
	+
	+ // K-matrix pole 5 is not included since its amplitude coefficient is zero
	+
	+ // K-matrix slowly varying part
	+ const EvtComplex fProd11( -4.724094278696236, -6.511009103363590 );
	+ const EvtComplex fProd12( -23.289333360304212, -12.215597571354197 );
	+ const EvtComplex fProd13( -1.860311896516422, -32.982507366353126 );
	+ const EvtComplex fProd14( -13.638752211193912, -22.339804683783186 );
	+ const EvtComplex fProd15( 0.0, 0.0 );
	+
	+ const EvtDalitzReso KMSVP( DP, BC, "f11prod", KMAT, fProd12 / fProd11,
	+ fProd13 / fProd11, fProd14 / fProd11,
	+ fProd15 / fProd11, sProd0 );
	+ m_resonances.push_back( std::make_pair( KMSVP, fProd11 ) );
	+}
	+
	+void EvtD0ToKspipi::setPDGValues()
	+{
	+ // Set the EvtIds
	+ m_BP = EvtPDL::getId( "B+" );
	+ m_BM = EvtPDL::getId( "B-" );
	+ m_B0 = EvtPDL::getId( "B0" );
	+ m_B0B = EvtPDL::getId( "anti-B0" );
	+ m_D0 = EvtPDL::getId( "D0" );
	+ m_D0B = EvtPDL::getId( "anti-D0" );
	+ m_KM = EvtPDL::getId( "K-" );
	+ m_KP = EvtPDL::getId( "K+" );
	+ m_K0 = EvtPDL::getId( "K0" );
	+ m_K0B = EvtPDL::getId( "anti-K0" );
	+ m_KL = EvtPDL::getId( "K_L0" );
	+ m_KS = EvtPDL::getId( "K_S0" );
	+ m_PIM = EvtPDL::getId( "pi-" );
	+ m_PIP = EvtPDL::getId( "pi+" );
	+
	+ // Set the particle masses
	+ m_mD0 = EvtPDL::getMeanMass( m_D0 );
	+ m_mKs = EvtPDL::getMeanMass( m_KS );
	+ m_mPi = EvtPDL::getMeanMass( m_PIP );
	+ m_mK = EvtPDL::getMeanMass( m_KP );
	+}
	diff --git a/src/EvtGenModels/EvtModelReg.cpp b/src/EvtGenModels/EvtModelReg.cpp
	--- a/src/EvtGenModels/EvtModelReg.cpp
	+++ b/src/EvtGenModels/EvtModelReg.cpp
	@@ -56,6 +56,7 @@
	#include "EvtGenModels/EvtBtoXsll.hh"
	#include "EvtGenModels/EvtCBTo3piMPP.hh"
	#include "EvtGenModels/EvtCBTo3piP00.hh"
	+#include "EvtGenModels/EvtD0ToKspipi.hh"
	#include "EvtGenModels/EvtD0gammaDalitz.hh"
	#include "EvtGenModels/EvtD0mixDalitz.hh"
	#include "EvtGenModels/EvtDDalitz.hh"
	@@ -286,6 +287,7 @@
	modelist.registerModel( new EvtDMix );
	modelist.registerModel( new EvtD0mixDalitz );
	modelist.registerModel( new EvtD0gammaDalitz );
	+ modelist.registerModel( new EvtD0ToKspipi );

	modelist.registerModel( new EvtbTosllAli );
	modelist.registerModel( new EvtBaryonPCR );

File Metadata

Mime Type: text/plain
Expires: Mon, Nov 18, 7:22 PM (10 h, 30 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3800938
Default Alt Text: D119.id502.diff (62 KB)

D119.id502.diffNo OneTemporaryActions

D119.id502.diffView Options

File Metadata

Event Timeline

D119.id502.diff
No OneTemporary
Actions

D119.id502.diff
View Options