diff --git a/EvtGenModels/EvtEtaLLPiPi.hh b/EvtGenModels/EvtEtaLLPiPi.hh
new file mode 100644
--- /dev/null
+++ b/EvtGenModels/EvtEtaLLPiPi.hh
@@ -0,0 +1,81 @@
+
+/***********************************************************************
+* Copyright 1998-2020 CERN for the benefit of the EvtGen authors *
+* *
+* This file is part of EvtGen. *
+* *
+* EvtGen is free software: you can redistribute it and/or modify *
+* it under the terms of the GNU General Public License as published by *
+* the Free Software Foundation, either version 3 of the License, or *
+* (at your option) any later version. *
+* *
+* EvtGen is distributed in the hope that it will be useful, *
+* but WITHOUT ANY WARRANTY; without even the implied warranty of *
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+* GNU General Public License for more details. *
+* *
+* You should have received a copy of the GNU General Public License *
+* along with EvtGen. If not, see . *
+***********************************************************************/
+
+#ifndef EVT_ETALLPIPI_HH
+#define EVT_ETALLPIPI_HH
+
+#include "EvtGenBase/EvtConst.hh"
+#include "EvtGenBase/EvtDecayProb.hh"
+
+#include
+
+class EvtParticle;
+
+// eta' -> mu+ mu- pi+ pi- or e+ e- pi+ pi-
+// From Zhang Zhen-Yu et al, Chinese Phys. C 36, p926, 2012
+
+class EvtEtaLLPiPi : public EvtDecayProb {
+ public:
+ EvtEtaLLPiPi() = default;
+
+ void init() override;
+ void initProbMax() override;
+
+ std::string getName() override;
+ EvtDecayBase* clone() override;
+
+ void decay( EvtParticle* p ) override;
+
+ private:
+ void updateMassPars( double mLep, double mPi );
+
+ double rhoWidth( double s, double m ) const;
+
+ double F0( double sLL, double sPiPi ) const;
+
+ double lambda( double a, double b, double c ) const;
+
+ double ampSquared( EvtParticle* p ) const;
+
+ double m_alpha{ 1.0 / 137.0 };
+ double m_eSq{ 4.0 * EvtConst::pi * m_alpha };
+ double m_fPi{ 0.0924 };
+ double m_f8{ 1.3 * m_fPi };
+ double m_f0{ 1.04 * m_fPi };
+ double m_thetaMix{ 20.0 * EvtConst::pi / 180.0 };
+ double m_mixSq{ 0.0 };
+ double m_c1{ 1.0 };
+ double m_c2{ 0.0 };
+ double m_c3{ m_c1 - m_c2 }; // Eq 9
+ double m_par1{ 1.0 - ( 3.0 * ( m_c1 - m_c2 + m_c3 ) / 4.0 ) };
+ double m_parLL{ 3.0 * ( m_c1 - m_c2 - m_c3 ) / 4.0 };
+ double m_parPiPi{ 3.0 * m_c3 / 2.0 };
+ double m_rhoMass{ 0.775 }; // updated in init()
+ double m_rhoMassSq{ m_rhoMass * m_rhoMass };
+ double m_rhoGamma{ 0.149 }; // updated in init()
+ double m_lepMass{ 0.106 }; // modified in updateMassPars()
+ double m_lepMassSq{ m_lepMass * m_lepMass };
+ double m_piMass{ 0.140 }; // modified in updateMassPars()
+ double m_piMassSq{ m_piMass * m_piMass };
+ double m_4LepMassSq{ 4.0 * m_lepMassSq };
+ double m_4PiMassSq{ 4.0 * m_piMassSq };
+};
+
+#endif
diff --git a/History.txt b/History.txt
--- a/History.txt
+++ b/History.txt
@@ -1,9 +1,25 @@
//==========================================================================
//
-// History file for EvtGen
+// History file for EvtGen.
+// From version 2.0.0, Tabc labels refer to Maniphest tasks while
+// Dxy labels refer to Differential code reviews on HepForge:
+// https://phab.hepforge.org/Tabc
+// https://phab.hepforge.org/Dxy
//
//===========================================================================
+21st Aug 2020 John Back
+ T109: Add EvtEtaLLPiPi model for eta' -> l l pi pi decays (l = e or mu),
+ courtesy of Aleksei Luchinsky (LHCb).
+
+29th Jun 2020 Michal Kreps
+ T103: Add missing include in EvtGenBase/EvtMatrix.hh.
+
+15th May 2020 Michal Kreps
+ D28: Add EvtThreeBodyPhsp (rectangular Dalitz plot selection) and
+ EvtPhspDecaytimeCut (minimum decay time requirement) models.
+ D27: Fix initialisation of constants in EvtBTo3hCP model.
+
R02-00-00-------------------------------------------------------------------
24th Apr 2020 Michal Kreps
Update particle properties file evt.pdl to 2019 version of RPP by PDG.
diff --git a/src/EvtGenModels/EvtEtaLLPiPi.cpp b/src/EvtGenModels/EvtEtaLLPiPi.cpp
new file mode 100644
--- /dev/null
+++ b/src/EvtGenModels/EvtEtaLLPiPi.cpp
@@ -0,0 +1,219 @@
+
+/***********************************************************************
+* Copyright 1998-2020 CERN for the benefit of the EvtGen authors *
+* *
+* This file is part of EvtGen. *
+* *
+* EvtGen is free software: you can redistribute it and/or modify *
+* it under the terms of the GNU General Public License as published by *
+* the Free Software Foundation, either version 3 of the License, or *
+* (at your option) any later version. *
+* *
+* EvtGen is distributed in the hope that it will be useful, *
+* but WITHOUT ANY WARRANTY; without even the implied warranty of *
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+* GNU General Public License for more details. *
+* *
+* You should have received a copy of the GNU General Public License *
+* along with EvtGen. If not, see . *
+***********************************************************************/
+
+#include "EvtGenModels/EvtEtaLLPiPi.hh"
+
+#include "EvtGenBase/EvtId.hh"
+#include "EvtGenBase/EvtKine.hh"
+#include "EvtGenBase/EvtPDL.hh"
+#include "EvtGenBase/EvtParticle.hh"
+#include "EvtGenBase/EvtSpinType.hh"
+#include "EvtGenBase/EvtVector4R.hh"
+
+#include
+
+// eta' -> mu+ mu- pi+ pi- or e+ e- pi+ pi-
+// From Zhang Zhen-Yu et al, Chinese Phys. C 36, p926, 2012
+
+void EvtEtaLLPiPi::init()
+{
+ // Check for 0 or 1 (maxProb) arguments
+ checkNArg( 0, 1 );
+
+ // Check particle types
+ checkSpinParent( EvtSpinType::SCALAR );
+ checkSpinDaughter( 0, EvtSpinType::DIRAC );
+ checkSpinDaughter( 1, EvtSpinType::DIRAC );
+ checkSpinDaughter( 2, EvtSpinType::SCALAR );
+ checkSpinDaughter( 3, EvtSpinType::SCALAR );
+
+ // Mass and width of rho0 from particle properties file
+ m_rhoMass = EvtPDL::getMeanMass( EvtPDL::getId( "rho0" ) );
+ m_rhoMassSq = m_rhoMass * m_rhoMass;
+ m_rhoGamma = EvtPDL::getWidth( EvtPDL::getId( "rho0" ) );
+
+ // Mixing parameter squared, using Eq 6
+ const double denom = 8.0 * pow( EvtConst::pi * m_fPi, 2 );
+ const double factor = m_eSq / ( denom * denom * 3.0 );
+ const double fTerm8 = sin( m_thetaMix ) / m_f8;
+ const double fTerm0 = 2.0 * sqrt( 2.0 ) * cos( m_thetaMix ) / m_f0;
+ m_mixSq = factor * pow( fTerm8 + fTerm0, 2 );
+}
+
+void EvtEtaLLPiPi::initProbMax()
+{
+ if ( getNArg() == 1 ) {
+ setProbMax( getArg( 0 ) );
+
+ } else {
+ int lepId = getDaug( 0 ).getId();
+ if ( lepId == EvtPDL::getId( "e-" ).getId() ||
+ lepId == EvtPDL::getId( "e+" ).getId() ) {
+ setProbMax( 27500.0 );
+
+ } else if ( lepId == EvtPDL::getId( "mu-" ).getId() ||
+ lepId == EvtPDL::getId( "mu+" ).getId() ) {
+ setProbMax( 20.0 );
+ }
+ }
+}
+
+std::string EvtEtaLLPiPi::getName()
+{
+ return "ETA_LLPIPI";
+}
+
+EvtDecayBase* EvtEtaLLPiPi::clone()
+{
+ return new EvtEtaLLPiPi();
+}
+
+void EvtEtaLLPiPi::decay( EvtParticle* p )
+{
+ p->initializePhaseSpace( getNDaug(), getDaugs() );
+
+ const double mLep = p->getDaug( 0 )->mass();
+ const double mPi = p->getDaug( 2 )->mass();
+
+ updateMassPars( mLep, mPi );
+
+ const double prob = ampSquared( p );
+ setProb( prob );
+}
+
+void EvtEtaLLPiPi::updateMassPars( double mLep, double mPi )
+{
+ // Update mass parameters used in various functions
+ m_lepMass = mLep;
+ m_lepMassSq = mLep * mLep;
+ m_4LepMassSq = 4.0 * m_lepMassSq;
+
+ m_piMass = mPi;
+ m_piMassSq = mPi * mPi;
+ m_4PiMassSq = 4.0 * m_piMassSq;
+}
+
+double EvtEtaLLPiPi::rhoWidth( double s, double m ) const
+{
+ // Define width of rho using modified vector meson dynamics, Eq 8
+ double gamma( 0.0 );
+
+ if ( s >= m_4PiMassSq ) {
+ const double mSq = m * m;
+ const double num = 1.0 - ( 4.0 * mSq / s );
+ const double denom = 1.0 - ( 4.0 * mSq / m_rhoMassSq );
+ const double ratio = denom > 0.0 ? num / denom : 0.0;
+ gamma = m_rhoGamma * ( s / m_rhoMassSq ) * pow( ratio, 1.5 );
+ }
+
+ return gamma;
+}
+
+double EvtEtaLLPiPi::F0( double sLL, double sPiPi ) const
+{
+ // Equation 7
+ double ampSq( 0.0 );
+
+ const double rhoWidthL = rhoWidth( sLL, m_lepMass );
+ const double rhoWidthPi = rhoWidth( sPiPi, m_piMass );
+
+ const double mSqDiffL = m_rhoMassSq - sLL;
+ const double mRhoWidthL = m_rhoMass * rhoWidthL;
+
+ const double mSqDiffPi = m_rhoMassSq - sPiPi;
+ const double mRhoWidthPi = m_rhoMass * rhoWidthPi;
+
+ const double denomLL = mSqDiffL * mSqDiffL + mRhoWidthL * mRhoWidthL;
+ const double denomPiPi = mSqDiffPi * mSqDiffPi + mRhoWidthPi * mRhoWidthPi;
+
+ if ( denomLL > 0.0 && denomPiPi > 0.0 ) {
+ const double mRho4 = m_rhoMassSq * m_rhoMassSq;
+ const double denomProd = denomLL * denomPiPi;
+
+ double realAmp = m_par1 + m_parLL * ( m_rhoMassSq * mSqDiffL / denomLL );
+ realAmp += m_parPiPi * mRho4 *
+ ( ( mSqDiffPi * mSqDiffL ) - mRhoWidthL * mRhoWidthPi ) /
+ denomProd;
+
+ double imagAmp = m_parLL * ( m_rhoMassSq * mRhoWidthL / denomLL );
+ imagAmp += m_parPiPi * mRho4 *
+ ( mRhoWidthPi * mSqDiffL + mRhoWidthL * mSqDiffPi ) /
+ denomProd;
+
+ ampSq = realAmp * realAmp + imagAmp * imagAmp;
+ }
+
+ return ampSq;
+}
+
+double EvtEtaLLPiPi::lambda( double a, double b, double c ) const
+{
+ const double sumSq = a * a + b * b + c * c;
+ const double prod = a * b + b * c + c * a;
+ const double L = sumSq - 2.0 * prod;
+ return L;
+}
+
+double EvtEtaLLPiPi::ampSquared( EvtParticle* p ) const
+{
+ // Equation 3
+ const double zeroProb( 0.0 );
+
+ // Mass of eta' meson
+ const double mEta = p->mass();
+
+ const EvtVector4R pL1 = p->getDaug( 0 )->getP4();
+ const EvtVector4R pL2 = p->getDaug( 1 )->getP4();
+ const EvtVector4R pPi1 = p->getDaug( 2 )->getP4();
+ const EvtVector4R pPi2 = p->getDaug( 3 )->getP4();
+
+ const EvtVector4R pLL = pL1 + pL2;
+ const double sLL = pLL.mass2();
+ const EvtVector4R pPiPi = pPi1 + pPi2;
+ const double sPiPi = pPiPi.mass2();
+
+ if ( sLL < 1e-4 || sPiPi < m_4PiMassSq || sLL < m_4LepMassSq ) {
+ // To avoid negative square roots etc
+ return zeroProb;
+ }
+
+ // Angles theta_p, theta_k and phi defined in Fig 1
+ const EvtVector4R pSum = pLL + pPiPi;
+ // Helicity angle of first lepton
+ const double cosThp = -EvtDecayAngle( pSum, pLL, pL1 );
+ const double sinThp = sqrt( 1.0 - cosThp * cosThp );
+ // Helicity angle of first pion
+ const double cosThk = -EvtDecayAngle( pSum, pPiPi, pPi2 );
+ const double sinThk = sqrt( 1.0 - cosThk * cosThk );
+ // Angle between the dilepton and dipion decay planes
+ const double phi = EvtDecayAngleChi( pSum, pL1, pL2, pPi1, pPi2 );
+ const double sinPhi = sin( phi );
+
+ const double betaLL = sqrt( 1.0 - ( m_4LepMassSq / sLL ) );
+ const double betaPiPi = sqrt( 1.0 - ( m_4PiMassSq / sPiPi ) );
+
+ const double betaProd = ( 1.0 - pow( betaLL * sinThp * sinPhi, 2 ) ) *
+ sPiPi * pow( betaPiPi * sinThk, 2 );
+ const double L = lambda( mEta * mEta, sLL, sPiPi );
+ const double ampSq = m_eSq * F0( sLL, sPiPi ) * m_mixSq * L * betaProd /
+ ( 8.0 * sLL );
+
+ return ampSq;
+}
diff --git a/src/EvtGenModels/EvtModelReg.cpp b/src/EvtGenModels/EvtModelReg.cpp
--- a/src/EvtGenModels/EvtModelReg.cpp
+++ b/src/EvtGenModels/EvtModelReg.cpp
@@ -160,6 +160,7 @@
#include "EvtGenModels/EvtbsToLLLL.hh"
#include "EvtGenModels/EvtbsToLLLLHyperCP.hh"
#include "EvtGenModels/EvtThreeBodyPhsp.hh"
+#include "EvtGenModels/EvtEtaLLPiPi.hh"
#include
#include
@@ -341,4 +342,5 @@
modelist.registerModel( new EvtPsi2JpsiPiPi );
modelist.registerModel( new EvtThreeBodyPhsp );
+ modelist.registerModel( new EvtEtaLLPiPi );
}