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 <https://www.gnu.org/licenses/>.     *
+***********************************************************************/
+
+#ifndef EVT_ETALLPIPI_HH
+#define EVT_ETALLPIPI_HH
+
+#include "EvtGenBase/EvtDecayProb.hh"
+#include <string>
+
+class EvtParticle;
+
+class EvtEtaLLPiPi : public EvtDecayProb {
+public:
+
+    EvtEtaLLPiPi();
+
+    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;
+    double m_eSq;
+    double m_fPi;
+    double m_f8;
+    double m_f0;
+    double m_thetaMix;
+    double m_mixSq;
+    double m_c1;
+    double m_c2;
+    double m_c3;
+    double m_par1;
+    double m_parLL;
+    double m_parPiPi;
+
+    double m_rhoMass;
+    double m_rhoMassSq;
+    double m_rhoGamma;
+
+    double m_lepMass;
+    double m_lepMassSq;
+    double m_piMass;
+    double m_piMassSq;
+    double m_4LepMassSq;
+    double m_4PiMassSq;
+    
+};
+
+#endif
diff --git a/History.txt b/History.txt
--- a/History.txt
+++ b/History.txt
@@ -4,6 +4,15 @@
 //
 //===========================================================================
 
+21st Aug 2020 John Back
+    Add EvtEtaLLPiPi model for eta' -> l l pi pi decays (l = e or mu),
+    courtesy of Aleksei Luchinsky (LHCb).
+
+15th May 2020 Michal Kreps
+    Add EvtThreeBodyPhsp (rectangular Dalitz plot selection) and
+    EvtPhspDecaytimeCut (minimum decay time requirement) models.
+    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,250 @@
+
+/***********************************************************************
+* 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 <https://www.gnu.org/licenses/>.     *
+***********************************************************************/
+
+#include "EvtGenModels/EvtEtaLLPiPi.hh"
+
+#include "EvtGenBase/EvtConst.hh"
+#include "EvtGenBase/EvtId.hh"
+#include "EvtGenBase/EvtKine.hh"
+#include "EvtGenBase/EvtParticle.hh"
+#include "EvtGenBase/EvtPDL.hh"
+#include "EvtGenBase/EvtSpinType.hh"
+#include "EvtGenBase/EvtVector4R.hh"
+
+#include <cmath>
+
+EvtEtaLLPiPi::EvtEtaLLPiPi() :
+    m_alpha( 1.0 / 137.0 ),
+    m_eSq( 4.0 * EvtConst::pi * m_alpha ),
+    m_fPi( 0.0924 ),
+    m_f8( 1.3 * m_fPi ),
+    m_f0( 1.04 * m_fPi ),
+    m_thetaMix( 20.0 * EvtConst::pi / 180.0 ),
+    m_mixSq( 0.0 ),
+    m_c1( 1.0 ),
+    m_c2( 0.0 ),
+    m_c3( m_c1 - m_c2 ), // Eq 9
+    m_par1( 1.0 - ( 3.0 * ( m_c1 - m_c2 + m_c3 ) / 4.0 ) ),
+    m_parLL( 3.0 * ( m_c1 - m_c2 - m_c3 ) / 4.0 ),
+    m_parPiPi( 3.0 * m_c3 / 2.0 ),
+    m_rhoMass( 0.775 ),  // updated in init()
+    m_rhoMassSq( m_rhoMass * m_rhoMass ),
+    m_rhoGamma( 0.149 ), // updated in init()
+    m_lepMass( 0.106 ),  // modified in updateMassPars()
+    m_lepMassSq ( m_lepMass * m_lepMass ),
+    m_piMass( 0.140 ),   // modified in updateMassPars()
+    m_piMassSq( m_piMass * m_piMass ),
+    m_4LepMassSq( 4.0 * m_lepMassSq ),
+    m_4PiMassSq( 4.0 * m_piMassSq )
+{
+    // 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 <assert.h>
 #include <ctype.h>
@@ -341,4 +342,5 @@
     modelist.registerModel( new EvtPsi2JpsiPiPi );
 
     modelist.registerModel( new EvtThreeBodyPhsp );
+    modelist.registerModel( new EvtEtaLLPiPi );
 }