D124.diff
No OneTemporary
Actions

Size

38 KB

Subscribers

None

D124.diff
View Options

	diff --git a/EvtGenBase/EvtDalitzReso.hh b/EvtGenBase/EvtDalitzReso.hh
	--- a/EvtGenBase/EvtDalitzReso.hh
	+++ b/EvtGenBase/EvtDalitzReso.hh
	@@ -124,7 +124,7 @@

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

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

	void set_fd( double R ) { m_vd.set_f( R ); }
	void set_fb( double R ) { m_vb.set_f( R ); }
	@@ -135,36 +135,39 @@
	}

	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 ) );
	diff --git a/EvtGenModels/EvtD0ToKspipi.hh b/EvtGenModels/EvtD0ToKspipi.hh
	new file mode 100644
	--- /dev/null
	+++ b/EvtGenModels/EvtD0ToKspipi.hh
	@@ -0,0 +1,67 @@
	+#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;
	+
	+ // 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
	@@ -11,6 +11,9 @@
	===
	## R02-0X-00

	+24 June 2024 John Back
	+* D124: Add EvtD0ToKspipi DP model, courtesy of Camille Normand (LHCb).
	+
	21 June 2024 Thomas Latham
	* D123: Removed `EvtPatches.hh`. Removed unused `clock_t` typedef and redundant `UNUSED` macro.

	diff --git a/src/EvtGenBase/EvtDalitzReso.cpp b/src/EvtGenBase/EvtDalitzReso.cpp
	--- a/src/EvtGenBase/EvtDalitzReso.cpp
	+++ b/src/EvtGenBase/EvtDalitzReso.cpp
	@@ -323,7 +323,7 @@
	m_spin = EvtSpinType::SCALAR;
	}

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

	@@ -349,10 +349,12 @@
	EvtComplex amp( 1.0, 0.0 );

	if ( fabs( m_dp.bigM() - x.bigM() ) > 0.000001 ) {
	- m_vb = EvtTwoBodyVertex( m_m0, m_dp.m( EvtCyclic3::other( m_pairRes ) ),
	- x.bigM(), m_spin );
	- m_vb.set_f( m_f_b );
	+ 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( m_pairRes ) ), x.bigM() );
	EvtTwoBodyKine vd( m_massFirst, m_massSecond, m );

	@@ -381,61 +383,76 @@
	switch ( m_coupling2 ) {
	case PicPic: {
	G1 = m_g1 * m_g1 * psFactor( m_massFirst, m_massSecond, m );
	- static double mPic = EvtPDL::getMass( EvtPDL::getId( "pi+" ) );
	+ static const double mPic = EvtPDL::getMass(
	+ EvtPDL::getId( "pi+" ) );
	G2 = m_g2 * m_g2 * psFactor( mPic, mPic, m );
	break;
	}
	case PizPiz: {
	G1 = m_g1 * m_g1 * psFactor( m_massFirst, m_massSecond, m );
	- static double mPiz = EvtPDL::getMass( EvtPDL::getId( "pi0" ) );
	+ static const double mPiz = EvtPDL::getMass(
	+ EvtPDL::getId( "pi0" ) );
	G2 = m_g2 * m_g2 * psFactor( mPiz, mPiz, m );
	break;
	}
	case PiPi: {
	G1 = m_g1 * m_g1 * psFactor( m_massFirst, m_massSecond, m );
	- static double mPic = EvtPDL::getMass( EvtPDL::getId( "pi+" ) );
	- static double mPiz = EvtPDL::getMass( EvtPDL::getId( "pi0" ) );
	+ 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 = m_g1 * m_g1 * psFactor( m_massFirst, m_massSecond, m );
	- static double mKc = EvtPDL::getMass( EvtPDL::getId( "K+" ) );
	+ static const double mKc = EvtPDL::getMass(
	+ EvtPDL::getId( "K+" ) );
	G2 = m_g2 * m_g2 * psFactor( mKc, mKc, m );
	break;
	}
	case KzKz: {
	G1 = m_g1 * m_g1 * psFactor( m_massFirst, m_massSecond, m );
	- static double mKz = EvtPDL::getMass( EvtPDL::getId( "K0" ) );
	+ static const double mKz = EvtPDL::getMass(
	+ EvtPDL::getId( "K0" ) );
	G2 = m_g2 * m_g2 * psFactor( mKz, mKz, m );
	break;
	}
	case KK: {
	G1 = m_g1 * m_g1 * psFactor( m_massFirst, m_massSecond, m );
	- static double mKc = EvtPDL::getMass( EvtPDL::getId( "K+" ) );
	- static double mKz = EvtPDL::getMass( EvtPDL::getId( "K0" ) );
	+ 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 = m_g1 * m_g1 * psFactor( m_massFirst, m_massSecond, m );
	- static double mEta = EvtPDL::getMass( EvtPDL::getId( "eta" ) );
	- static double mPic = EvtPDL::getMass( EvtPDL::getId( "pi+" ) );
	+ 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 = m_g1 * m_g1 * psFactor( m_massFirst, m_massSecond, m );
	- static double mEta = EvtPDL::getMass( EvtPDL::getId( "eta" ) );
	- static double mPiz = EvtPDL::getMass( EvtPDL::getId( "pi0" ) );
	+ 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+" ) );
	+ static const double mPic = EvtPDL::getMass(
	+ EvtPDL::getId( "pi+" ) );
	G1 = m_g1 * psFactor( mPic, mPic, m );
	- static double mKc = EvtPDL::getMass( EvtPDL::getId( "K+" ) );
	- static double mKz = EvtPDL::getMass( EvtPDL::getId( "K0" ) );
	+ 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;
	}
	@@ -476,7 +493,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 +508,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 +525,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 +540,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 +549,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 +569,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,19 +580,19 @@
	( 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. * 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 * m_massFirst * m_massFirst / ( k0 * k0 ) *
	log( ( m0 + 2. * k0 ) / ( 2. * m_massFirst ) ) +
	@@ -582,7 +602,7 @@

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

	@@ -657,7 +677,7 @@
	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
	@@ -673,14 +693,14 @@
	acos( cosTh ) );
	}

	-EvtComplex EvtDalitzReso::mixFactor( EvtComplex prop, EvtComplex prop_mix )
	+EvtComplex EvtDalitzReso::mixFactor( EvtComplex prop, EvtComplex prop_mix ) const
	{
	double Delta = m_delta_mix * ( m_m0 + m_m0_mix );
	return 1 / ( 1 - Delta * Delta * prop * 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 );

	@@ -965,7 +985,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.

	@@ -1010,7 +1030,7 @@
	}

	//replace Breit-Wigner with LASS
	-EvtComplex EvtDalitzReso::lass( double s )
	+EvtComplex EvtDalitzReso::lass( double s ) const
	{
	EvtTwoBodyKine vd( m_massFirst, m_massSecond, sqrt( s ) );
	double q = vd.p();
	@@ -1045,7 +1065,7 @@
	return T;
	}

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

	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,325 @@
	+#include "EvtGenModels/EvtD0ToKspipi.hh"
	+
	+#include "EvtGenBase/EvtDecayTable.hh"
	+#include "EvtGenBase/EvtPDL.hh"
	+#include "EvtGenBase/EvtParticle.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 );
	+
	+ // 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();
	+ EvtDecayBase* decayFun = ( parent != nullptr )
	+ ? EvtDecayTable::getInstance()->getDecayFunc(
	+ parent )
	+ : nullptr;
	+ if ( parent != nullptr && decayFun != nullptr &&
	+ decayFun->getName() == "BTODDALITZCPK" ) {
	+ const EvtId parId = parent->getId();
	+ if ( ( parId == m_BP ) \|\| ( parId == m_BM ) \|\| ( parId == m_B0 ) \|\|
	+ ( parId == m_B0B ) ) {
	+ // 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 = decayFun->getArg( 0 );
	+ // Strong phase in radians
	+ const double delta = decayFun->getArg( 1 );
	+ // Ratio between B -> D0 K and B -> D0bar K
	+ const double rB = decayFun->getArg( 2 );
	+
	+ // Direct and conjugated amplitudes
	+ const EvtComplex ampD0 = calcTotAmp( pointD0 );
	+ const EvtComplex ampD0b = calcTotAmp( pointD0b );
	+
	+ if ( parId == m_BP \|\| parId == 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 ( p->getId() == 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 ( const 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, -1.0, false );
	+ 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, -1.0, false );
	+ 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
	@@ -55,6 +55,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"
	@@ -273,6 +274,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 );
	diff --git a/test/jsonFiles/BTODDALITZCPK=D0TOKSPIPI__Bu_D0K_KS0pipi.json b/test/jsonFiles/BTODDALITZCPK=D0TOKSPIPI__Bu_D0K_KS0pipi.json
	new file mode 100644
	--- /dev/null
	+++ b/test/jsonFiles/BTODDALITZCPK=D0TOKSPIPI__Bu_D0K_KS0pipi.json
	@@ -0,0 +1,29 @@
	+{
	+ "parent" : "B-",
	+ "daughters" : ["D0", "K-"],
	+ "grand_daughters" : [["K_S0", "pi+", "pi-"], []],
	+ "models" : ["BTODDALITZCPK", "D0TOKSPIPI", ""],
	+ "parameters" : [["1.22", "2.27", "0.10"], [], []],
	+ "extras" : ["noFSR"],
	+ "events" : 5000,
	+ "histograms" : [
	+ {"variable" : "prob_daug1", "title" : "Prob(D0TOKSPIPI)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 1.0},
	+ {"variable" : "parMass", "title" : "m(B^{-})", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 5.25, "xmax" : 5.3},
	+ {"variable" : "cosTheta", "title" : "cosTheta(D^{0})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosTheta", "title" : "cosTheta(K^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "mass_daug1", "title" : "m(K^{0}_{S} #pi^{+})", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 0.5, "xmax" : 1.8},
	+ {"variable" : "mass_daug1", "title" : "m(K^{0}_{S} #pi^{-})", "d1" : 1, "d2" : 3, "nbins" : 100, "xmin" : 0.5, "xmax" : 1.8},
	+ {"variable" : "mass_daug1", "title" : "m(#pi^{+} #pi^{-})", "d1" : 2, "d2" : 3, "nbins" : 100, "xmin" : 0.2, "xmax" : 1.5},
	+ {"variable" : "massSq_daug1", "title" : "m^{2}(K^{0}_{S} #pi^{+})", "d1" : 1, "d2" : 2, "nbins" : 150, "xmin" : 0.0, "xmax" : 3.0},
	+ {"variable" : "massSq_daug1", "title" : "m^{2}(K^{0}_{S} #pi^{-})", "d1" : 1, "d2" : 3, "nbins" : 150, "xmin" : 0.0, "xmax" : 3.0},
	+ {"variable" : "massSq_daug1", "title" : "m^{2}(#pi^{+} #pi^{-})", "d1" : 2, "d2" : 3, "nbins" : 100, "xmin" : 0.0, "xmax" : 2.0},
	+ {"variable" : "massSq_daug1", "title" : "Dalitz m^{2}(K^{0}_{S} #pi^{+}) vs m^{2}(K^{0}_{S} #pi^{-})", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 3.1, "variableY" : "massSq_daug1", "d1Y" : 1, "d2Y" : 3, "nbinsY" : 50, "ymin" : 0.0, "ymax" : 3.1},
	+ {"variable" : "massSq_daug1", "title" : "Dalitz m^{2}(#pi^{+} #pi^{-}) vs m^{2}(K^{0}_{S} #pi^{+})", "d1" : 2, "d2" : 3, "nbins" : 50, "xmin" : 0.0, "xmax" : 2.0, "variableY" : "massSq_daug1", "d1Y" : 1, "d2Y" : 2, "nbinsY" : 50, "ymin" : 0.0, "ymax" : 3.1},
	+ {"variable" : "mPrime_daug1", "title" : "sqDalitz", "d1" : 1, "d2" : 2, "nbins" : 25, "xmin" : 0.0, "xmax" : 1.00, "variableY" : "thetaPrime_daug1", "d1Y" : 1, "d2Y" : 2, "nbinsY" : 25, "ymin" : 0.0, "ymax" : 1.0},
	+ {"variable" : "cosHel", "title" : "cosHel(D^{0},K^{0}_{S})", "d1" : 1, "d2" : 0, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosHel_daug1", "title" : "cosHel12", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosHel_daug1", "title" : "cosHel23", "d1" : 2, "d2" : 3, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosHel_daug1", "title" : "cosHel13", "d1" : 1, "d2" : 3, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0}
	+ ],
	+ "outfile" : "BTODDALITZCPK=D0TOKSPIPI__Bu_D0K_KS0pipi.root"
	+}
	diff --git a/test/jsonFiles/D0TOKSPIPI__D0_KS0pipi.json b/test/jsonFiles/D0TOKSPIPI__D0_KS0pipi.json
	new file mode 100644
	--- /dev/null
	+++ b/test/jsonFiles/D0TOKSPIPI__D0_KS0pipi.json
	@@ -0,0 +1,29 @@
	+{
	+ "parent" : "D0",
	+ "daughters" : ["K_S0", "pi+", "pi-"],
	+ "models" : ["D0TOKSPIPI"],
	+ "parameters" : [[]],
	+ "extras" : ["noFSR"],
	+ "events" : 5000,
	+ "histograms" : [
	+ {"variable" : "prob", "title" : "Prob", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 1.0},
	+ {"variable" : "parMass", "title" : "m(D^{0})", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 1.8, "xmax" : 2.0},
	+ {"variable" : "cosTheta", "title" : "cosTheta(K^{0}_{S})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosTheta", "title" : "cosTheta(#pi^{+})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosTheta", "title" : "cosTheta(#pi^{-})", "d1" : 3, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "mass", "title" : "m(K^{0}_{S} #pi^{+})", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 0.5, "xmax" : 1.8},
	+ {"variable" : "mass", "title" : "m(K^{0}_{S} #pi^{-})", "d1" : 1, "d2" : 3, "nbins" : 100, "xmin" : 0.5, "xmax" : 1.8},
	+ {"variable" : "mass", "title" : "m(#pi^{+} #pi^{-})", "d1" : 2, "d2" : 3, "nbins" : 100, "xmin" : 0.2, "xmax" : 1.5},
	+ {"variable" : "massSq", "title" : "m^{2}(K^{0}_{S} #pi^{+})", "d1" : 1, "d2" : 2, "nbins" : 150, "xmin" : 0.0, "xmax" : 3.0},
	+ {"variable" : "massSq", "title" : "m^{2}(K^{0}_{S} #pi^{-})", "d1" : 1, "d2" : 3, "nbins" : 150, "xmin" : 0.0, "xmax" : 3.0},
	+ {"variable" : "massSq", "title" : "m^{2}(#pi^{+} #pi^{-})", "d1" : 2, "d2" : 3, "nbins" : 100, "xmin" : 0.0, "xmax" : 2.0},
	+ {"variable" : "massSq", "title" : "Dalitz m^{2}(K^{0}_{S} #pi^{+}) vs m^{2}(K^{0}_{S} #pi^{-})", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 3.1, "variableY" : "massSq", "d1Y" : 1, "d2Y" : 3, "nbinsY" : 50, "ymin" : 0.0, "ymax" : 3.1},
	+ {"variable" : "massSq", "title" : "Dalitz m^{2}(#pi^{+} #pi^{-}) vs m^{2}(K^{0}_{S} #pi^{+})", "d1" : 2, "d2" : 3, "nbins" : 50, "xmin" : 0.0, "xmax" : 2.0, "variableY" : "massSq", "d1Y" : 1, "d2Y" : 2, "nbinsY" : 50, "ymin" : 0.0, "ymax" : 3.1},
	+ {"variable" : "mPrime", "title" : "sqDalitz", "d1" : 1, "d2" : 2, "nbins" : 25, "xmin" : 0.0, "xmax" : 1.00, "variableY" : "thetaPrime", "d1Y" : 1, "d2Y" : 2, "nbinsY" : 25, "ymin" : 0.0, "ymax" : 1.0},
	+ {"variable" : "cosHel", "title" : "cosHel(D^{0},K^{0}_{S})", "d1" : 1, "d2" : 0, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosHel", "title" : "cosHel12", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosHel", "title" : "cosHel23", "d1" : 2, "d2" : 3, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosHel", "title" : "cosHel13", "d1" : 1, "d2" : 3, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0}
	+ ],
	+ "outfile" : "D0TOKSPIPI__D0_KS0pipi.root"
	+}

File Metadata

Mime Type: text/plain
Expires: Mon, Nov 18, 6:47 PM (16 h, 27 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3805669
Default Alt Text: D124.diff (38 KB)

D124.diffNo OneTemporaryActions

D124.diffView Options

File Metadata

Event Timeline

D124.diff
No OneTemporary
Actions

D124.diff
View Options