D100.1759161376.diff
No OneTemporary
Actions

Size

84 KB

Referenced Files

None

Subscribers

None

D100.1759161376.diff
View Options

	diff --git a/EvtGenBase/EvtVector4R.hh b/EvtGenBase/EvtVector4R.hh
	--- a/EvtGenBase/EvtVector4R.hh
	+++ b/EvtGenBase/EvtVector4R.hh
	@@ -53,7 +53,7 @@
	void applyRotateEuler( double alpha, double beta, double gamma );
	void applyBoostTo( const EvtVector4R& p4, bool inverse = false );
	void applyBoostTo( const EvtVector3R& boost, bool inverse = false );
	- EvtVector4R cross( const EvtVector4R& v2 );
	+ EvtVector4R cross( const EvtVector4R& v2 ) const;
	double dot( const EvtVector4R& v2 ) const;
	double d3mag() const;

	diff --git a/EvtGenModels/EvtBcVHad.hh b/EvtGenModels/EvtBcVHad.hh
	--- a/EvtGenModels/EvtBcVHad.hh
	+++ b/EvtGenModels/EvtBcVHad.hh
	@@ -27,6 +27,7 @@
	#include "EvtGenModels/EvtBCVFF2.hh"
	#include "EvtGenModels/EvtWHad.hh"

	+#include <array>
	#include <memory>
	#include <string>

	@@ -40,37 +41,44 @@
	EvtDecayBase* clone() override;
	void initProbMax() override;
	void init() override;
	- void decay( EvtParticle* p ) override;
	+ void decay( EvtParticle* parent ) override;

	protected:
	// Hadronic current function
	- EvtVector4C hardCurr( EvtParticle* root_particle ) const;
	+ EvtVector4C hardCurr( EvtParticle* parent ) const;
	+ void parseDecay();

	private:
	- // whichfit --- code of the Bc -> VW formfactor set:
	- // 1 - SR
	- // 2 - PM
	- int whichfit;
	+ // Code of the Bc -> VW formfactor set:
	+ // 1 - SR
	+ // 2 - PM
	+ int m_whichFit;

	- // final vector particle code
	- int idVector;
	+ // Final vector particle code
	+ int m_idVector;

	- // out_code: code of the hadronic final state
	- // 1 - pi+
	- // 2 - pi+ pi0
	- // 3 - pi+ pi+ pi-
	- // 4 - 4pi
	- // 5 - pi+ pi+ pi- pi- pi+
	- // 6 - K+ K- pi+
	- // 7 - K+ pi+ pi-
	- // 8 - K_S0 K+
	- // 9 - K+ K- pi+ pi+ pi-
	- // 10 - pi+ pi+ pi+ pi+ pi- pi- pi- (7 charged pi)
	- // 11 - K+ pi+ pi+ pi- pi-
	- int out_code;
	+ // Code of the hadronic final state
	+ // 1: B_c+ -> V pi+
	+ // 2: B_c+ -> V pi+ pi0
	+ // 3: B_c+ -> V 2pi+ pi-
	+ // 4: B_c+ -> V 2pi+ pi- pi0 (not implemented)
	+ // 5: B_c+ -> V 3pi+ 2pi-
	+ // 6: B_c+ -> V K+ K- pi+
	+ // 7: B_c+ -> V K+ pi+ pi-
	+ // 8: B_c+ -> V K_S0 K+
	+ // 9: B_c+ -> V K+ K- 2pi+ pi-
	+ // 10: B_c+ -> V 4pi+ 3pi-
	+ // 11: B_c+ -> V K+ 2pi+ 2pi-
	+ int m_outCode;

	- std::unique_ptr<EvtBCVFF2> ffmodel;
	- std::unique_ptr<EvtWHad> wcurr;
	+ std::unique_ptr<EvtBCVFF2> m_FFModel;
	+ std::unique_ptr<EvtWHad> m_WCurr;
	+
	+ std::array<int, 4> m_iPiPlus = { { -1, -1, -1, -1 } };
	+ std::array<int, 4> m_iPiMinus = { { -1, -1, -1, -1 } };
	+ std::array<int, 4> m_iPiZero = { { -1, -1, -1, -1 } };
	+ std::array<int, 4> m_iKPlus = { { -1, -1, -1, -1 } };
	+ std::array<int, 4> m_iKMinus = { { -1, -1, -1, -1 } };
	};

	#endif
	diff --git a/EvtGenModels/EvtPhiDalitz.hh b/EvtGenModels/EvtBcVPPHad.hh
	copy from EvtGenModels/EvtPhiDalitz.hh
	copy to EvtGenModels/EvtBcVPPHad.hh
	--- a/EvtGenModels/EvtPhiDalitz.hh
	+++ b/EvtGenModels/EvtBcVPPHad.hh
	@@ -1,6 +1,5 @@
	-
	/***********************************************************************
	-* Copyright 1998-2020 CERN for the benefit of the EvtGen authors *
	+* Copyright 1998-2023 CERN for the benefit of the EvtGen authors *
	* *
	* This file is part of EvtGen. *
	* *
	@@ -18,32 +17,48 @@
	* along with EvtGen. If not, see <https://www.gnu.org/licenses/>. *
	***********************************************************************/

	-#ifndef EVTPHIDALITZ_HH
	-#define EVTPHIDALITZ_HH
	+#ifndef EvtBcVPPHad_HH
	+#define EvtBcVPPHad_HH

	#include "EvtGenBase/EvtDecayAmp.hh"
	+#include "EvtGenBase/EvtVector4C.hh"
	+
	+#include "EvtGenModels/EvtBCVFF2.hh"
	+#include "EvtGenModels/EvtWHad.hh"
	+
	+#include <string>

	class EvtParticle;

	-class EvtPhiDalitz : public EvtDecayAmp {
	+// Description: Module to implement Bc -> psi + p + pbar + pi decays
	+
	+class EvtBcVPPHad : public EvtDecayAmp {
	public:
	std::string getName() override;
	EvtDecayBase* clone() override;
	-
	- void init() override;
	void initProbMax() override;
	- void decay( EvtParticle* p ) override;
	+ void init() override;
	+ void decay( EvtParticle* parent ) override;
	+
	+ protected:
	+ // Hadronic current function
	+ EvtVector4C hardCurrPP( EvtParticle* parent, int i1, int i2 ) const;

	private:
	- double _mRho;
	- double _gRho;
	- double _aD;
	- double _phiD;
	- double _aOmega;
	- double _phiOmega;
	- int _locPip;
	- int _locPim;
	- int _locPi0;
	+ // Code of the Bc -> VW formfactor set:
	+ // 1 - SR
	+ // 2 - PM
	+ int m_whichFit;
	+
	+ // Final vector particle code
	+ int m_idVector;
	+
	+ // Code of the hadronic final state
	+ // 1: p+ p- pi+
	+ int m_outCode;
	+
	+ std::unique_ptr<EvtBCVFF2> m_FFModel;
	+ std::unique_ptr<EvtWHad> m_WCurr;
	};

	#endif
	diff --git a/EvtGenModels/EvtPhiDalitz.hh b/EvtGenModels/EvtPhiDalitz.hh
	--- a/EvtGenModels/EvtPhiDalitz.hh
	+++ b/EvtGenModels/EvtPhiDalitz.hh
	@@ -35,15 +35,17 @@
	void decay( EvtParticle* p ) override;

	private:
	- double _mRho;
	- double _gRho;
	- double _aD;
	- double _phiD;
	- double _aOmega;
	- double _phiOmega;
	- int _locPip;
	- int _locPim;
	- int _locPi0;
	+ double calc_q( double, double, double ) const;
	+
	+ double m_mRho;
	+ double m_gRho;
	+ double m_aD;
	+ double m_phiD;
	+ double m_aOmega;
	+ double m_phiOmega;
	+ int m_locPip;
	+ int m_locPim;
	+ int m_locPi0;
	};

	#endif
	diff --git a/EvtGenModels/EvtWHad.hh b/EvtGenModels/EvtWHad.hh
	--- a/EvtGenModels/EvtWHad.hh
	+++ b/EvtGenModels/EvtWHad.hh
	@@ -22,6 +22,7 @@
	#define EvtWHad_HH

	#include "EvtGenBase/EvtComplex.hh"
	+#include "EvtGenBase/EvtDiracSpinor.hh"
	#include "EvtGenBase/EvtVector4C.hh"
	#include "EvtGenBase/EvtVector4R.hh"

	@@ -72,11 +73,16 @@
	const EvtVector4R& p5, const EvtVector4R& p6,
	const EvtVector4R& p7 ) const;

	- // 1=K+ 2 = pi+ 3 = pi+ 4 = pi- 5 = pi- with symmetrizatiom
	+ // 1=K+ 2 = pi+ 3 = pi+ 4 = pi- 5 = pi- with symmetrization
	EvtVector4C WCurrent_K4pi( const EvtVector4R& p1, const EvtVector4R& p2,
	const EvtVector4R& p3, const EvtVector4R& p4,
	const EvtVector4R& p5 ) const;

	+ // a1 -> p+ p- pi+
	+ EvtVector4C WCurrent_ppPi( const EvtVector4R& p1, const EvtDiracSpinor& sp1,
	+ const EvtVector4R& p2, const EvtDiracSpinor& sp2,
	+ const EvtVector4R& k ) const;
	+
	protected:
	EvtVector4C JB( const EvtVector4R& q1, const EvtVector4R& q2,
	const EvtVector4R& q3, const EvtVector4R& q4,
	diff --git a/History.md b/History.md
	--- a/History.md
	+++ b/History.md
	@@ -11,11 +11,17 @@
	===
	## R02-0X-00

	+3 Feb 2024 John Back
	+* T230: Add EvtBcVPPHad model for Bc to Jpsi p pbar pi decays and generalised particle ordering
	+ for EvtBcVHad decay files, courtesy of Aleksei Luchinsky, Dmitrii Pereima & Vanya Belyaev (LHCb).
	+ Updated EvtPhiDalitz model to use helicity amplitudes and fixed the indices used in
	+ EvtVector3R::dot(), courtesy of Arnau Brossa Gonzalo & Antonio Romero Vidal (LHCb).
	+
	31 Oct 2023 Fernando Abudinen
	* D102: Bugfix probmax issue for X38722-+\_PSI\_GAMMA model
	- Calculation of amplitude moved to dedicated calcAmp function and fixed initialisation.
	- Introduced weighting to compensate for phase-space changes due to resonance width.
	- Fixed bugs in rho couplings and loops over photon and vector states.
	+ - Calculation of amplitude moved to dedicated calcAmp function and fixed initialisation.
	+ - Introduced weighting to compensate for phase-space changes due to resonance width.
	+ - Fixed bugs in rho couplings and loops over photon and vector states.

	16 Oct 2023 Fernando Abudinen
	* D99: Add tests for PHOTOS
	@@ -34,7 +40,7 @@

	22 Aug 2023 Fernando Abudinen
	* D97: Bugfix probmax issue and introduced pole compensation for VTOSLL model
	- Calculation of amplitude moved to dedicated calcAmp function and fixed initialisation.
	+ - Calculation of amplitude moved to dedicated calcAmp function and fixed initialisation.

	22 Aug 2023 Tom Latham
	* CMake updates
	diff --git a/src/EvtGenBase/EvtVector3R.cpp b/src/EvtGenBase/EvtVector3R.cpp
	--- a/src/EvtGenBase/EvtVector3R.cpp
	+++ b/src/EvtGenBase/EvtVector3R.cpp
	@@ -87,24 +87,23 @@
	}

	double EvtVector3R::d3mag() const
	-
	-// returns the 3 momentum mag.
	{
	+ // Returns the 3 momentum mag
	double temp;

	temp = v[0] * v[0] + v[1] * v[1] + v[2] * v[2];
	temp = sqrt( temp );

	return temp;
	-} // r3mag
	+}

	double EvtVector3R::dot( const EvtVector3R& p2 )
	{
	double temp;

	temp = v[0] * p2.v[0];
	- temp += v[0] * p2.v[0];
	- temp += v[0] * p2.v[0];
	+ temp += v[1] * p2.v[1];
	+ temp += v[2] * p2.v[2];

	return temp;
	-} //dot
	+}
	diff --git a/src/EvtGenBase/EvtVector4R.cpp b/src/EvtGenBase/EvtVector4R.cpp
	--- a/src/EvtGenBase/EvtVector4R.cpp
	+++ b/src/EvtGenBase/EvtVector4R.cpp
	@@ -170,7 +170,7 @@
	}
	}

	-EvtVector4R EvtVector4R::cross( const EvtVector4R& p2 )
	+EvtVector4R EvtVector4R::cross( const EvtVector4R& p2 ) const
	{
	//Calcs the cross product. Added by djl on July 27, 1995.
	//Modified for real vectros by ryd Aug 28-96
	diff --git a/src/EvtGenModels/EvtBcVHad.cpp b/src/EvtGenModels/EvtBcVHad.cpp
	--- a/src/EvtGenModels/EvtBcVHad.cpp
	+++ b/src/EvtGenModels/EvtBcVHad.cpp
	@@ -47,6 +47,19 @@

	void EvtBcVHad::init()
	{
	+ // The following decay m_outCodes are supported:
	+ // 1: B_c+ -> V pi+
	+ // 2: B_c+ -> V pi+ pi0
	+ // 3: B_c+ -> V 2pi+ pi-
	+ // 4: B_c+ -> V 2pi+ pi- pi0 (not implemented)
	+ // 5: B_c+ -> V 3pi+ 2pi-
	+ // 6: B_c+ -> V K+ K- pi+
	+ // 7: B_c+ -> V K+ pi+ pi-
	+ // 8: B_c+ -> V K_S0 K+
	+ // 9: B_c+ -> V K+ K- 2pi+ pi-
	+ // 10: B_c+ -> V 4pi+ 3pi-
	+ // 11: B_c+ -> V K+ 2pi+ 2pi-
	+
	checkNArg( 1 );
	checkSpinParent( EvtSpinType::SCALAR );
	checkSpinDaughter( 0, EvtSpinType::VECTOR );
	@@ -54,277 +67,323 @@
	checkSpinDaughter( i, EvtSpinType::SCALAR );
	}

	- idVector = getDaug( 0 ).getId();
	- whichfit = int( getArg( 0 ) + 0.1 );
	- ffmodel = std::make_unique<EvtBCVFF2>( idVector, whichfit );
	-
	- wcurr = std::make_unique<EvtWHad>();
	-
	- // determine the code of final hadronic state
	- EvtIdSet thePis{ "pi+", "pi-", "pi0" };
	- EvtIdSet theK{ "K+", "K-", "K_S0" };
	- if ( getNDaug() == 2 && thePis.contains( getDaug( 1 ) ) ) {
	- out_code = 1; // pi+
	- } else if ( getNDaug() == 3 && thePis.contains( getDaug( 1 ) ) &&
	- thePis.contains( getDaug( 2 ) ) ) {
	- out_code = 2; // pi+ pi0
	- } else if ( getNDaug() == 4 && thePis.contains( getDaug( 1 ) ) &&
	- thePis.contains( getDaug( 2 ) ) &&
	- thePis.contains( getDaug( 3 ) ) ) {
	- out_code = 3; // pi+ pi+ pi-
	- } else if ( getNDaug() == 5 && thePis.contains( getDaug( 1 ) ) &&
	- thePis.contains( getDaug( 2 ) ) &&
	- thePis.contains( getDaug( 3 ) ) &&
	- thePis.contains( getDaug( 4 ) ) ) {
	- out_code = 4; //4pi
	- } else if ( getNDaug() == 6 && thePis.contains( getDaug( 1 ) ) &&
	- thePis.contains( getDaug( 2 ) ) &&
	- thePis.contains( getDaug( 3 ) ) &&
	- thePis.contains( getDaug( 4 ) ) &&
	- thePis.contains( getDaug( 5 ) ) ) {
	- out_code = 5; //5pi
	- } else if ( getNDaug() == 4 && theK.contains( getDaug( 1 ) ) &&
	- theK.contains( getDaug( 2 ) ) &&
	- thePis.contains( getDaug( 3 ) ) ) {
	- out_code = 6; //KKpi
	- } else if ( getNDaug() == 4 && theK.contains( getDaug( 1 ) ) &&
	- thePis.contains( getDaug( 2 ) ) &&
	- thePis.contains( getDaug( 3 ) ) ) {
	- out_code = 7; // Kpi pi
	+ m_idVector = getDaug( 0 ).getId();
	+ m_whichFit = int( getArg( 0 ) + 0.1 );
	+ m_FFModel = std::make_unique<EvtBCVFF2>( m_idVector, m_whichFit );
	+
	+ m_WCurr = std::make_unique<EvtWHad>();
	+
	+ // Determine the code of the final hadronic state
	+ parseDecay();
	+}
	+
	+//======================================================
	+
	+void EvtBcVHad::parseDecay()
	+{
	+ const EvtIdSet BcPlusID{ "B_c+" }, BcMinusID{ "B_c-" };
	+ const EvtIdSet theK{ "K+", "K-", "K_S0" };
	+ const int cMode = BcMinusID.contains( getParentId() );
	+
	+ const EvtIdSet PiPlusID{ cMode == 0 ? "pi+" : "pi-" };
	+ const EvtIdSet PiMinusID{ cMode == 0 ? "pi-" : "pi+" };
	+ const EvtIdSet PiZeroID{ "pi0" };
	+ const EvtIdSet KPlusID{ cMode == 0 ? "K+" : "K-" };
	+ const EvtIdSet KMinusID{ cMode == 0 ? "K-" : "K+" };
	+ EvtGenReport( EVTGEN_INFO, "EvtBcVHad" )
	+ << "parentId = " << getParentId() << std::endl;
	+
	+ int PiPlusFound = 0, PiMinusFound = 0, PiZeroFound = 0, KPlusFound = 0,
	+ KMinusFound = 0;
	+ for ( int iDaughter = 0; iDaughter < getNDaug(); ++iDaughter ) {
	+ const EvtId daugId = getDaug( iDaughter );
	+ EvtGenReport( EVTGEN_INFO, "EvtBcVHad" )
	+ << "iDaughter = " << iDaughter
	+ << " id = " << getDaug( iDaughter ).getName() << std::endl;
	+ if ( PiPlusID.contains( daugId ) && PiPlusFound < 4 ) {
	+ m_iPiPlus[PiPlusFound] = iDaughter;
	+ PiPlusFound++;
	+ } else if ( PiMinusID.contains( daugId ) && PiMinusFound < 4 ) {
	+ m_iPiMinus[PiMinusFound] = iDaughter;
	+ PiMinusFound++;
	+ } else if ( PiZeroID.contains( daugId ) && PiZeroFound < 4 ) {
	+ m_iPiZero[PiZeroFound] = iDaughter;
	+ PiZeroFound++;
	+ } else if ( KPlusID.contains( daugId ) && KPlusFound < 4 ) {
	+ m_iKPlus[KPlusFound] = iDaughter;
	+ KPlusFound++;
	+ } else if ( KMinusID.contains( daugId ) && KMinusFound < 4 ) {
	+ m_iKMinus[KMinusFound] = iDaughter;
	+ KMinusFound++;
	+ }
	+ }
	+
	+ if ( getNDaug() == 2 && PiPlusFound == 1 ) {
	+ m_outCode = 1; // pi+
	+ } else if ( getNDaug() == 3 && PiPlusFound == 1 && PiZeroFound == 1 ) {
	+ m_outCode = 2; // pi+ pi0
	+ } else if ( getNDaug() == 4 && PiPlusFound == 2 && PiMinusFound == 1 ) {
	+ m_outCode = 3; // pi+ pi+ pi-
	+ } else if ( getNDaug() == 5 && PiPlusFound == 2 && PiMinusFound == 1 &&
	+ PiZeroFound == 1 ) {
	+ m_outCode = 4; // pi+ pi+ pi- pi0
	+ } else if ( getNDaug() == 6 && PiPlusFound == 3 && PiMinusFound == 2 ) {
	+ m_outCode = 5; // 5pi
	+ } else if ( getNDaug() == 4 && KPlusFound == 1 && KMinusFound == 1 &&
	+ PiPlusFound == 1 ) {
	+ m_outCode = 6; // KKpi
	+ } else if ( getNDaug() == 4 && KPlusFound == 1 && PiPlusFound == 1 &&
	+ PiMinusFound == 1 ) {
	+ m_outCode = 7; // K+ pi+ pi-
	} else if ( getNDaug() == 3 && theK.contains( getDaug( 1 ) ) &&
	theK.contains( getDaug( 2 ) ) ) {
	- out_code = 8; // KK
	- } else if ( getNDaug() == 6 && theK.contains( getDaug( 1 ) ) &&
	- theK.contains( getDaug( 2 ) ) && thePis.contains( getDaug( 3 ) ) &&
	- thePis.contains( getDaug( 4 ) ) &&
	- thePis.contains( getDaug( 5 ) ) ) {
	- out_code = 9; // KK+3pi
	- } else if ( getNDaug() == 8 && thePis.contains( getDaug( 1 ) ) &&
	- thePis.contains( getDaug( 2 ) ) &&
	- thePis.contains( getDaug( 3 ) ) &&
	- thePis.contains( getDaug( 4 ) ) &&
	- thePis.contains( getDaug( 5 ) ) &&
	- thePis.contains( getDaug( 6 ) ) &&
	- thePis.contains( getDaug( 7 ) ) ) {
	- out_code = 10; //7pi
	- } else if ( getNDaug() == 6 && theK.contains( getDaug( 1 ) ) &&
	- thePis.contains( getDaug( 2 ) ) &&
	- thePis.contains( getDaug( 3 ) ) &&
	- thePis.contains( getDaug( 4 ) ) &&
	- thePis.contains( getDaug( 5 ) ) ) {
	- out_code = 11; // K+ pi+ pi+ pi- pi-
	+ m_outCode = 8; // KK
	+ } else if ( getNDaug() == 6 && KPlusFound == 1 && KMinusFound == 1 &&
	+ PiPlusFound == 2 && PiMinusFound == 1 ) {
	+ m_outCode = 9; // K+ K- pi+ pi+ pi-
	+ } else if ( getNDaug() == 8 && PiPlusFound == 4 && PiMinusFound == 3 ) {
	+ m_outCode = 10; // 7pi
	+ } else if ( getNDaug() == 6 && KPlusFound == 1 && PiPlusFound == 2 &&
	+ PiMinusFound == 2 ) {
	+ m_outCode = 11; // K+ pi+ pi+ pi- pi-
	} else {
	- EvtGenReport( EVTGEN_ERROR, "EvtBcHad" )
	+ EvtGenReport( EVTGEN_ERROR, "EvtBcVHad" )
	<< "Init: unknown decay" << std::endl;
	- };
	+ }

	- EvtGenReport( EVTGEN_INFO, "EvtBcHad" )
	- << "out_code = " << out_code << ", whichfit = " << whichfit << std::endl;
	+ EvtGenReport( EVTGEN_INFO, "EvtBcVHad" )
	+ << "m_outCode = " << m_outCode << ", m_whichFit = " << m_whichFit
	+ << std::endl;
	+ for ( int i = 0; i < 4; ++i ) {
	+ EvtGenReport( EVTGEN_INFO, "EvtBcVHad" )
	+ << " i = " << i << ", m_iPiPlus = " << m_iPiPlus[i]
	+ << ", m_iPiMinus = " << m_iPiMinus[i]
	+ << ", m_iPiZero = " << m_iPiZero[i] << ", m_iKPlus = " << m_iKPlus[i]
	+ << ", m_iKMinus = " << m_iKMinus[i] << std::endl;
	+ }
	+ EvtGenReport( EVTGEN_INFO, "EvtBcVHad" )
	+ << "PiPlusFound = " << PiPlusFound << ", PiMinusFound = " << PiMinusFound
	+ << ", PiZeroFound = " << PiZeroFound << ", KPlusFound = " << KPlusFound
	+ << ", KMinusFound = " << KMinusFound << std::endl;
	}

	//======================================================

	void EvtBcVHad::initProbMax()
	{
	- if ( out_code == 1 ) {
	- if ( idVector == EvtPDL::getId( "J/psi" ).getId() && whichfit == 1 &&
	- getNDaug() == 2 )
	+ if ( m_outCode == 1 ) {
	+ if ( m_idVector == EvtPDL::getId( "J/psi" ).getId() &&
	+ m_whichFit == 1 && getNDaug() == 2 )
	setProbMax( 500. );
	- else if ( idVector == EvtPDL::getId( "J/psi" ).getId() &&
	- whichfit == 2 && getNDaug() == 2 )
	+ else if ( m_idVector == EvtPDL::getId( "J/psi" ).getId() &&
	+ m_whichFit == 2 && getNDaug() == 2 )
	setProbMax( 300. );
	- else if ( idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
	- whichfit == 1 && getNDaug() == 2 )
	+ else if ( m_idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
	+ m_whichFit == 1 && getNDaug() == 2 )
	setProbMax( 17. );
	- else if ( idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
	- whichfit == 2 && getNDaug() == 2 )
	+ else if ( m_idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
	+ m_whichFit == 2 && getNDaug() == 2 )
	setProbMax( 40. );
	- } else if ( out_code == 2 ) {
	- if ( idVector == EvtPDL::getId( "J/psi" ).getId() && whichfit == 1 &&
	- getNDaug() == 3 )
	+ } else if ( m_outCode == 2 ) {
	+ if ( m_idVector == EvtPDL::getId( "J/psi" ).getId() &&
	+ m_whichFit == 1 && getNDaug() == 3 )
	setProbMax( 10950. );
	- else if ( idVector == EvtPDL::getId( "J/psi" ).getId() &&
	- whichfit == 2 && getNDaug() == 3 )
	+ else if ( m_idVector == EvtPDL::getId( "J/psi" ).getId() &&
	+ m_whichFit == 2 && getNDaug() == 3 )
	setProbMax( 4200. );
	- else if ( idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
	- whichfit == 1 && getNDaug() == 3 )
	+ else if ( m_idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
	+ m_whichFit == 1 && getNDaug() == 3 )
	setProbMax( 500. );
	- else if ( idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
	- whichfit == 2 && getNDaug() == 3 )
	+ else if ( m_idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
	+ m_whichFit == 2 && getNDaug() == 3 )
	setProbMax( 700. );
	- } else if ( out_code == 3 ) {
	- if ( idVector == EvtPDL::getId( "J/psi" ).getId() && whichfit == 1 &&
	- getNDaug() == 4 )
	+ } else if ( m_outCode == 3 ) {
	+ if ( m_idVector == EvtPDL::getId( "J/psi" ).getId() &&
	+ m_whichFit == 1 && getNDaug() == 4 )
	setProbMax( 42000. );
	- else if ( idVector == EvtPDL::getId( "J/psi" ).getId() &&
	- whichfit == 2 && getNDaug() == 4 )
	+ else if ( m_idVector == EvtPDL::getId( "J/psi" ).getId() &&
	+ m_whichFit == 2 && getNDaug() == 4 )
	setProbMax( 90000. );
	- else if ( idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
	- whichfit == 1 && getNDaug() == 4 )
	+ else if ( m_idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
	+ m_whichFit == 1 && getNDaug() == 4 )
	setProbMax( 1660. );
	- else if ( idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
	- whichfit == 2 && getNDaug() == 4 )
	+ else if ( m_idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
	+ m_whichFit == 2 && getNDaug() == 4 )
	setProbMax( 2600. );
	- } else if ( out_code == 5 ) {
	- if ( idVector == EvtPDL::getId( "J/psi" ).getId() && whichfit == 1 &&
	- getNDaug() == 6 )
	+ } else if ( m_outCode == 5 ) {
	+ if ( m_idVector == EvtPDL::getId( "J/psi" ).getId() &&
	+ m_whichFit == 1 && getNDaug() == 6 )
	setProbMax( 720000. );
	- else if ( idVector == EvtPDL::getId( "J/psi" ).getId() &&
	- whichfit == 2 && getNDaug() == 6 )
	+ else if ( m_idVector == EvtPDL::getId( "J/psi" ).getId() &&
	+ m_whichFit == 2 && getNDaug() == 6 )
	setProbMax( 519753. );
	- else if ( idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
	- whichfit == 1 && getNDaug() == 6 )
	+ else if ( m_idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
	+ m_whichFit == 1 && getNDaug() == 6 )
	setProbMax( 40000. );
	- else if ( idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
	- whichfit == 2 && getNDaug() == 6 )
	+ else if ( m_idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
	+ m_whichFit == 2 && getNDaug() == 6 )
	setProbMax( 30000. );
	- } else if ( out_code == 6 ) {
	- if ( idVector == EvtPDL::getId( "J/psi" ).getId() && whichfit == 1 )
	+ } else if ( m_outCode == 6 ) {
	+ if ( m_idVector == EvtPDL::getId( "J/psi" ).getId() && m_whichFit == 1 )
	setProbMax( 50000. );
	- else if ( idVector == EvtPDL::getId( "J/psi" ).getId() && whichfit == 2 )
	+ else if ( m_idVector == EvtPDL::getId( "J/psi" ).getId() &&
	+ m_whichFit == 2 )
	setProbMax( 22000.0 );
	- else if ( idVector == EvtPDL::getId( "psi(2S)" ).getId() && whichfit == 1 )
	+ else if ( m_idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
	+ m_whichFit == 1 )
	setProbMax( 2300.0 );
	- else if ( idVector == EvtPDL::getId( "psi(2S)" ).getId() && whichfit == 2 )
	+ else if ( m_idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
	+ m_whichFit == 2 )
	setProbMax( 1700.00 );
	- } else if ( out_code == 7 ) {
	- if ( idVector == EvtPDL::getId( "J/psi" ).getId() && whichfit == 1 )
	+ } else if ( m_outCode == 7 ) {
	+ if ( m_idVector == EvtPDL::getId( "J/psi" ).getId() && m_whichFit == 1 )
	setProbMax( 2.2e+06 );
	- else if ( idVector == EvtPDL::getId( "J/psi" ).getId() && whichfit == 2 )
	+ else if ( m_idVector == EvtPDL::getId( "J/psi" ).getId() &&
	+ m_whichFit == 2 )
	setProbMax( 930000 );
	- else if ( idVector == EvtPDL::getId( "psi(2S)" ).getId() && whichfit == 1 )
	+ else if ( m_idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
	+ m_whichFit == 1 )
	setProbMax( 92000.0 );
	- else if ( idVector == EvtPDL::getId( "psi(2S)" ).getId() && whichfit == 2 )
	+ else if ( m_idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
	+ m_whichFit == 2 )
	setProbMax( 93000.0 );
	- } else if ( out_code == 8 ) {
	- if ( idVector == EvtPDL::getId( "J/psi" ).getId() && whichfit == 1 )
	+ } else if ( m_outCode == 8 ) {
	+ if ( m_idVector == EvtPDL::getId( "J/psi" ).getId() && m_whichFit == 1 )
	setProbMax( 2e2 );
	- else if ( idVector == EvtPDL::getId( "J/psi" ).getId() && whichfit == 2 )
	+ else if ( m_idVector == EvtPDL::getId( "J/psi" ).getId() &&
	+ m_whichFit == 2 )
	setProbMax( 80 );
	- else if ( idVector == EvtPDL::getId( "psi(2S)" ).getId() && whichfit == 1 )
	+ else if ( m_idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
	+ m_whichFit == 1 )
	setProbMax( 10 );
	- else if ( idVector == EvtPDL::getId( "psi(2S)" ).getId() && whichfit == 2 )
	+ else if ( m_idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
	+ m_whichFit == 2 )
	setProbMax( 10 );
	- } else if ( out_code == 9 ) {
	- if ( idVector == EvtPDL::getId( "J/psi" ).getId() && whichfit == 1 )
	+ } else if ( m_outCode == 9 ) {
	+ if ( m_idVector == EvtPDL::getId( "J/psi" ).getId() && m_whichFit == 1 )
	setProbMax( 3e4 );
	- else if ( idVector == EvtPDL::getId( "J/psi" ).getId() && whichfit == 2 )
	+ else if ( m_idVector == EvtPDL::getId( "J/psi" ).getId() &&
	+ m_whichFit == 2 )
	setProbMax( 18540 );
	- else if ( idVector == EvtPDL::getId( "psi(2S)" ).getId() && whichfit == 1 )
	+ else if ( m_idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
	+ m_whichFit == 1 )
	setProbMax( 0.15 * 1e4 );
	- else if ( idVector == EvtPDL::getId( "psi(2S)" ).getId() && whichfit == 2 )
	+ else if ( m_idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
	+ m_whichFit == 2 )
	setProbMax( 2 * 500 );
	- } else if ( out_code == 10 ) {
	- if ( idVector == EvtPDL::getId( "J/psi" ).getId() && whichfit == 1 )
	+ } else if ( m_outCode == 10 ) {
	+ if ( m_idVector == EvtPDL::getId( "J/psi" ).getId() && m_whichFit == 1 )
	setProbMax( 2e6 );
	- else if ( idVector == EvtPDL::getId( "J/psi" ).getId() && whichfit == 2 )
	+ else if ( m_idVector == EvtPDL::getId( "J/psi" ).getId() &&
	+ m_whichFit == 2 )
	setProbMax( 5e6 );
	- else if ( idVector == EvtPDL::getId( "psi(2S)" ).getId() && whichfit == 1 )
	+ else if ( m_idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
	+ m_whichFit == 1 )
	setProbMax( 1.5e5 );
	- else if ( idVector == EvtPDL::getId( "psi(2S)" ).getId() && whichfit == 2 )
	+ else if ( m_idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
	+ m_whichFit == 2 )
	setProbMax( 1e5 );
	- } else if ( out_code == 11 ) {
	- if ( idVector == EvtPDL::getId( "J/psi" ).getId() && whichfit == 1 )
	+ } else if ( m_outCode == 11 ) {
	+ if ( m_idVector == EvtPDL::getId( "J/psi" ).getId() && m_whichFit == 1 )
	setProbMax( 2.5e8 );
	- else if ( idVector == EvtPDL::getId( "J/psi" ).getId() && whichfit == 2 )
	+ else if ( m_idVector == EvtPDL::getId( "J/psi" ).getId() &&
	+ m_whichFit == 2 )
	setProbMax( 1.4e7 );
	- else if ( idVector == EvtPDL::getId( "psi(2S)" ).getId() && whichfit == 1 )
	+ else if ( m_idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
	+ m_whichFit == 1 )
	setProbMax( 2e6 );
	- else if ( idVector == EvtPDL::getId( "psi(2S)" ).getId() && whichfit == 2 )
	+ else if ( m_idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
	+ m_whichFit == 2 )
	setProbMax( 8e4 );
	} else {
	EvtGenReport( EVTGEN_ERROR, "EvtBcHad" )
	- << "probmax: Have not yet implemented this final state in BC_VHAD model, out_code = "
	- << out_code << std::endl;
	- ::abort();
	+ << "probmax: Have not yet implemented this final state in BC_VHAD model, m_outCode = "
	+ << m_outCode << std::endl;
	+ for ( int id = 0; id < ( getNDaug() - 1 ); id++ ) {
	+ EvtGenReport( EVTGEN_ERROR, "EvtBcVHad" )
	+ << "Daug " << id << " " << EvtPDL::name( getDaug( id ) ).c_str()
	+ << std::endl;
	+ }
	}
	}

	//======================================================

	-EvtVector4C EvtBcVHad::hardCurr( EvtParticle* root_particle ) const
	+EvtVector4C EvtBcVHad::hardCurr( EvtParticle* parent ) const
	{
	EvtVector4C hardCur;

	- if ( out_code == 1 ) {
	+ if ( m_outCode == 1 ) {
	// pi+
	- hardCur = wcurr->WCurrent( root_particle->getDaug( 1 )->getP4() );
	+ hardCur = m_WCurr->WCurrent( parent->getDaug( m_iPiPlus[0] )->getP4() );

	- } else if ( out_code == 2 ) {
	+ } else if ( m_outCode == 2 ) {
	// pi+ pi0
	- hardCur = wcurr->WCurrent( root_particle->getDaug( 1 )->getP4(),
	- root_particle->getDaug( 2 )->getP4() );
	+ hardCur = m_WCurr->WCurrent( parent->getDaug( m_iPiPlus[0] )->getP4(),
	+ parent->getDaug( m_iPiZero[0] )->getP4() );

	- } else if ( out_code == 3 ) {
	+ } else if ( m_outCode == 3 ) {
	// pi+ pi+ pi-
	- hardCur = wcurr->WCurrent( root_particle->getDaug( 1 )->getP4(),
	- root_particle->getDaug( 2 )->getP4(),
	- root_particle->getDaug( 3 )->getP4() );
	-
	- } else if ( out_code == 5 ) {
	- // Bc -> psi pi+ pi+ pi- pi- pi+ from Kuhn, Was, hep-ph/0602162
	- hardCur = wcurr->WCurrent( root_particle->getDaug( 1 )->getP4(),
	- root_particle->getDaug( 2 )->getP4(),
	- root_particle->getDaug( 3 )->getP4(),
	- root_particle->getDaug( 4 )->getP4(),
	- root_particle->getDaug( 5 )->getP4() );
	-
	- } else if ( out_code == 6 ) {
	+ hardCur = m_WCurr->WCurrent( parent->getDaug( m_iPiPlus[0] )->getP4(),
	+ parent->getDaug( m_iPiPlus[1] )->getP4(),
	+ parent->getDaug( m_iPiMinus[0] )->getP4() );
	+ } else if ( m_outCode == 5 ) {
	+ // Bc -> psi pi+ pi+ pi- pi- pi+ from Kuhn & Was, hep-ph/0602162
	+ hardCur =
	+ m_WCurr->WCurrent_5pi( parent->getDaug( m_iPiPlus[0] )->getP4(),
	+ parent->getDaug( m_iPiPlus[1] )->getP4(),
	+ parent->getDaug( m_iPiPlus[2] )->getP4(),
	+ parent->getDaug( m_iPiMinus[0] )->getP4(),
	+ parent->getDaug( m_iPiMinus[1] )->getP4() );
	+
	+ } else if ( m_outCode == 6 ) {
	// K+ K- pi+
	- hardCur = wcurr->WCurrent_KKP( root_particle->getDaug( 1 )->getP4(),
	- root_particle->getDaug( 2 )->getP4(),
	- root_particle->getDaug( 3 )->getP4() );
	+ hardCur =
	+ m_WCurr->WCurrent_KKP( parent->getDaug( m_iKPlus[0] )->getP4(),
	+ parent->getDaug( m_iKMinus[0] )->getP4(),
	+ parent->getDaug( m_iPiPlus[0] )->getP4() );

	- } else if ( out_code == 7 ) {
	+ } else if ( m_outCode == 7 ) {
	// K+ pi+ pi-
	- hardCur = wcurr->WCurrent_KPP( root_particle->getDaug( 1 )->getP4(),
	- root_particle->getDaug( 2 )->getP4(),
	- root_particle->getDaug( 3 )->getP4() );
	+ hardCur =
	+ m_WCurr->WCurrent_KPP( parent->getDaug( m_iKPlus[0] )->getP4(),
	+ parent->getDaug( m_iPiPlus[0] )->getP4(),
	+ parent->getDaug( m_iPiMinus[0] )->getP4() );

	- } else if ( out_code == 8 ) {
	+ } else if ( m_outCode == 8 ) {
	// K_S0 K+
	- hardCur = wcurr->WCurrent_KSK( root_particle->getDaug( 1 )->getP4(),
	- root_particle->getDaug( 2 )->getP4() );
	+ hardCur = m_WCurr->WCurrent_KSK( parent->getDaug( 1 )->getP4(),
	+ parent->getDaug( 2 )->getP4() );

	- } else if ( out_code == 9 ) {
	+ } else if ( m_outCode == 9 ) {
	// K+ K- pi+ pi+ pi-
	- hardCur = wcurr->WCurrent_KKPPP(
	- root_particle->getDaug( 1 )->getP4(), // K+
	- root_particle->getDaug( 2 )->getP4(), // K-
	- root_particle->getDaug( 3 )->getP4(), // pi+
	- root_particle->getDaug( 4 )->getP4(), // pi+
	- root_particle->getDaug( 5 )->getP4() // pi-
	+ hardCur = m_WCurr->WCurrent_KKPPP(
	+ parent->getDaug( m_iKPlus[0] )->getP4(), // K+
	+ parent->getDaug( m_iKMinus[0] )->getP4(), // K-
	+ parent->getDaug( m_iPiPlus[0] )->getP4(), // pi+
	+ parent->getDaug( m_iPiPlus[1] )->getP4(), // pi+
	+ parent->getDaug( m_iPiMinus[0] )->getP4() // pi-
	);
	- } else if ( out_code == 10 ) {
	- // 1=pi+ 2=pi+ 3=pi+ 4=pi+ 5=pi- 6=pi- 7=pi- with symmetrization of the identical particles
	- hardCur =
	- wcurr->WCurrent_7pi( root_particle->getDaug( 1 )->getP4(), // pi+
	- root_particle->getDaug( 2 )->getP4(), // pi+
	- root_particle->getDaug( 3 )->getP4(), // pi+
	- root_particle->getDaug( 4 )->getP4(), // pi+
	- root_particle->getDaug( 5 )->getP4(), // pi-
	- root_particle->getDaug( 6 )->getP4(), // pi-
	- root_particle->getDaug( 7 )->getP4() // pi-
	- );
	- } else if ( out_code == 11 ) {
	- // 1=K+ 2 = pi+ 3 = pi+ 4 = pi- 5 = pi- with symmetrizatiom
	- hardCur = wcurr->WCurrent_K4pi(
	- root_particle->getDaug( 1 )->getP4(), // K+
	- root_particle->getDaug( 2 )->getP4(), // pi+
	- root_particle->getDaug( 3 )->getP4(), // pi+
	- root_particle->getDaug( 4 )->getP4(), // pi-
	- root_particle->getDaug( 5 )->getP4() // pi-
	+ } else if ( m_outCode == 10 ) {
	+ // 1=pi+ 2=pi+ 3=pi+ 4=pi+ 5=pi- 6=pi- 7=pi- with symmetrization of identical particles
	+ hardCur = m_WCurr->WCurrent_7pi(
	+ parent->getDaug( m_iPiPlus[0] )->getP4(), // pi+
	+ parent->getDaug( m_iPiPlus[1] )->getP4(), // pi+
	+ parent->getDaug( m_iPiPlus[2] )->getP4(), // pi+
	+ parent->getDaug( m_iPiPlus[3] )->getP4(), // pi+
	+ parent->getDaug( m_iPiMinus[0] )->getP4(), // pi-
	+ parent->getDaug( m_iPiMinus[1] )->getP4(), // pi-
	+ parent->getDaug( m_iPiMinus[2] )->getP4() // pi-
	+ );
	+ } else if ( m_outCode == 11 ) {
	+ // 1=K+ 2 = pi+ 3 = pi+ 4 = pi- 5 = pi- with symmetrization
	+ hardCur = m_WCurr->WCurrent_K4pi(
	+ parent->getDaug( m_iKPlus[0] )->getP4(), // K+
	+ parent->getDaug( m_iPiPlus[0] )->getP4(), // pi+
	+ parent->getDaug( m_iPiPlus[1] )->getP4(), // pi+
	+ parent->getDaug( m_iPiMinus[0] )->getP4(), // pi-
	+ parent->getDaug( m_iPiMinus[1] )->getP4() // pi-
	);
	- } else {
	- EvtGenReport( EVTGEN_ERROR, "EvtBcHad" )
	- << "hardCurr: Have not yet implemented this final state in BC_VHAD model"
	- << std::endl;
	- // EvtGenReport(EVTGEN_ERROR, "EvtGen") << "Have not yet implemented this final state in BC_VHAD model" << std::endl;
	- // EvtGenReport(EVTGEN_ERROR, "EvtGen") << "Ndaug=" << getNDaug() << std::endl;
	- for ( int id = 0; id < ( getNDaug() - 1 ); id++ ) {
	- // EvtGenReport(EVTGEN_ERROR, "EvtGen") << "Daug " << id << " " << EvtPDL::name(getDaug(id)).c_str() << std::endl;
	- }
	- ::abort();
	}

	return hardCur;
	@@ -332,42 +391,43 @@

	//======================================================

	-void EvtBcVHad::decay( EvtParticle* root_particle )
	+void EvtBcVHad::decay( EvtParticle* parent )
	{
	- root_particle->initializePhaseSpace( getNDaug(), getDaugs() );
	+ parent->initializePhaseSpace( getNDaug(), getDaugs() );

	// Calculate hadronic current
	- const EvtVector4C hardCur = hardCurr( root_particle );
	+ const EvtVector4C hardCur = hardCurr( parent );

	- EvtParticle* Jpsi = root_particle->getDaug( 0 );
	+ EvtParticle* Jpsi = parent->getDaug( 0 );
	+ //std::cout<<"4th comp: "<<Jpsi->eps(3)<<std::endl;

	- const EvtVector4R p4b( root_particle->mass(), 0., 0., 0. ), // Bc momentum
	- p4meson = Jpsi->getP4(), // J/psi momenta
	+ const EvtVector4R p4b( parent->mass(), 0., 0., 0. ), // Bc momentum
	+ p4meson = Jpsi->getP4(), // J/psi momenta
	Q = p4b - p4meson, p4Sum = p4meson + p4b;
	const double Q2 = Q.mass2();

	// Calculate Bc -> V W form-factors
	double a1f( 0.0 ), a2f( 0.0 ), vf( 0.0 ), a0f( 0.0 );

	- const double m_meson = Jpsi->mass();
	- const double m_b = root_particle->mass();
	- const double mVar = m_b + m_meson;
	+ const double mMeson = Jpsi->mass();
	+ const double mB = parent->mass();
	+ const double mVar = mB + mMeson;

	- ffmodel->getvectorff( root_particle->getId(), Jpsi->getId(), Q2, m_meson,
	- &a1f, &a2f, &vf, &a0f );
	+ m_FFModel->getvectorff( parent->getId(), Jpsi->getId(), Q2, mMeson, &a1f,
	+ &a2f, &vf, &a0f );

	- const double a3f = ( mVar / ( 2.0 * m_meson ) ) * a1f -
	- ( ( m_b - m_meson ) / ( 2.0 * m_meson ) ) * a2f;
	+ const double a3f = ( mVar / ( 2.0 * mMeson ) ) * a1f -
	+ ( ( mB - mMeson ) / ( 2.0 * mMeson ) ) * a2f;

	// Calculate Bc -> V W current
	EvtTensor4C H = a1f * mVar * EvtTensor4C::g();
	H.addDirProd( ( -a2f / mVar ) * p4b, p4Sum );
	H += EvtComplex( 0.0, vf / mVar ) *
	dual( EvtGenFunctions::directProd( p4Sum, Q ) );
	- H.addDirProd( ( a0f - a3f ) * 2.0 * ( m_meson / Q2 ) * p4b, Q );
	+ H.addDirProd( ( a0f - a3f ) * 2.0 * ( mMeson / Q2 ) * p4b, Q );
	const EvtVector4C Heps = H.cont2( hardCur );

	- for ( int i = 0; i < 4; i++ ) {
	+ for ( int i = 0; i < 3; i++ ) {
	const EvtVector4C eps =
	Jpsi->epsParent( i ).conj(); // psi-meson polarization vector
	const EvtComplex amp = eps * Heps;
	diff --git a/src/EvtGenModels/EvtBcVPPHad.cpp b/src/EvtGenModels/EvtBcVPPHad.cpp
	new file mode 100644
	--- /dev/null
	+++ b/src/EvtGenModels/EvtBcVPPHad.cpp
	@@ -0,0 +1,166 @@
	+/***********************************************************************
	+* Copyright 1998-2023 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/EvtBcVPPHad.hh"
	+
	+#include "EvtGenBase/EvtDiracSpinor.hh"
	+#include "EvtGenBase/EvtIdSet.hh"
	+#include "EvtGenBase/EvtPDL.hh"
	+#include "EvtGenBase/EvtParticle.hh"
	+#include "EvtGenBase/EvtReport.hh"
	+#include "EvtGenBase/EvtSpinType.hh"
	+#include "EvtGenBase/EvtTensor4C.hh"
	+#include "EvtGenBase/EvtVector4R.hh"
	+
	+#include <iostream>
	+
	+std::string EvtBcVPPHad::getName()
	+{
	+ return "BC_VPPHAD";
	+}
	+
	+EvtDecayBase* EvtBcVPPHad::clone()
	+{
	+ return new EvtBcVPPHad;
	+}
	+
	+//======================================================
	+
	+void EvtBcVPPHad::init()
	+{
	+ checkNArg( 1 );
	+ checkSpinParent( EvtSpinType::SCALAR );
	+ checkSpinDaughter( 0, EvtSpinType::VECTOR );
	+ checkSpinDaughter( 1, EvtSpinType::DIRAC );
	+ checkSpinDaughter( 2, EvtSpinType::DIRAC );
	+ for ( int i = 3; i <= ( getNDaug() - 1 ); i++ ) {
	+ checkSpinDaughter( i, EvtSpinType::SCALAR );
	+ }
	+
	+ m_idVector = getDaug( 0 ).getId();
	+ m_whichFit = int( getArg( 0 ) + 0.1 );
	+ m_FFModel = std::make_unique<EvtBCVFF2>( m_idVector, m_whichFit );
	+
	+ m_WCurr = std::make_unique<EvtWHad>();
	+
	+ // determine the code of final hadronic state
	+ const EvtIdSet thePis{ "pi+", "pi-", "pi0" };
	+ const EvtIdSet theK{ "K+", "K-", "K_S0" };
	+ const EvtIdSet thePs{ "p+", "anti-p-" };
	+
	+ if ( getNDaug() == 4 && thePs.contains( getDaug( 1 ) ) &&
	+ thePs.contains( getDaug( 2 ) ) && thePis.contains( getDaug( 3 ) ) ) {
	+ m_outCode = 1; // p+ p- pi+
	+ } else {
	+ EvtGenReport( EVTGEN_ERROR, "EvtBcPPHad::init has unknown decay mode" );
	+ }
	+ EvtGenReport( EVTGEN_INFO, "EvtBcVPPHad" )
	+ << ": m_outCode = " << m_outCode << ", m_whichFit = " << m_whichFit
	+ << std::endl;
	+}
	+
	+//======================================================
	+
	+void EvtBcVPPHad::initProbMax()
	+{
	+ if ( m_outCode == 1 ) {
	+ // p+ p- pi+
	+ if ( m_idVector == EvtPDL::getId( "J/psi" ).getId() ) {
	+ if ( m_whichFit == 2 ) {
	+ setProbMax( 550.0 );
	+ } else {
	+ setProbMax( 1100.0 );
	+ }
	+ } else if ( m_idVector == EvtPDL::getId( "psi(2S)" ).getId() ) {
	+ if ( m_whichFit == 2 ) {
	+ setProbMax( 7.0 );
	+ } else {
	+ setProbMax( 50.0 );
	+ }
	+ }
	+ } else {
	+ EvtGenReport(
	+ EVTGEN_ERROR,
	+ "probmax: Have not yet implemented this final state in BC_VPPHAD model, m_outCode = " )
	+ << m_outCode << std::endl;
	+ }
	+}
	+
	+//======================================================
	+EvtVector4C EvtBcVPPHad::hardCurrPP( EvtParticle* parent, int i1, int i2 ) const
	+{
	+ EvtVector4C hardCur;
	+ const EvtVector4R p1 = parent->getDaug( 1 )->getP4();
	+ const EvtDiracSpinor sp1 = parent->getDaug( 1 )->spParent( i1 );
	+
	+ const EvtVector4R p2 = parent->getDaug( 2 )->getP4();
	+ const EvtDiracSpinor sp2 = parent->getDaug( 2 )->spParent( i2 );
	+
	+ if ( m_outCode == 1 ) {
	+ const EvtVector4R k = parent->getDaug( 3 )->getP4();
	+ hardCur = m_WCurr->WCurrent_ppPi( p1, sp1, p2, sp2, k );
	+ }
	+ return hardCur;
	+}
	+
	+//======================================================
	+void EvtBcVPPHad::decay( EvtParticle* parent )
	+{
	+ parent->initializePhaseSpace( getNDaug(), getDaugs() );
	+
	+ EvtParticle* Jpsi = parent->getDaug( 0 );
	+
	+ const EvtVector4R p4b( parent->mass(), 0., 0., 0. ), // Bc momentum
	+ p4meson = Jpsi->getP4(), // J/psi momenta
	+ Q = p4b - p4meson, p4Sum = p4meson + p4b;
	+ const double Q2 = Q.mass2();
	+
	+ // Calculate Bc -> V W form-factors
	+ double a1f( 0.0 ), a2f( 0.0 ), vf( 0.0 ), a0f( 0.0 );
	+
	+ const double mMeson = Jpsi->mass();
	+ const double mB = parent->mass();
	+ const double mVar = mB + mMeson;
	+
	+ m_FFModel->getvectorff( parent->getId(), Jpsi->getId(), Q2, mMeson, &a1f,
	+ &a2f, &vf, &a0f );
	+
	+ const double a3f = ( mVar / ( 2.0 * mMeson ) ) * a1f -
	+ ( ( mB - mMeson ) / ( 2.0 * mMeson ) ) * a2f;
	+
	+ // Calculate Bc -> V W current
	+ EvtTensor4C H = a1f * mVar * EvtTensor4C::g();
	+ H.addDirProd( ( -a2f / mVar ) * p4b, p4Sum );
	+ H += EvtComplex( 0.0, vf / mVar ) *
	+ dual( EvtGenFunctions::directProd( p4Sum, Q ) );
	+ H.addDirProd( ( a0f - a3f ) * 2.0 * ( mMeson / Q2 ) * p4b, Q );
	+
	+ for ( int i1 = 0; i1 < 2; ++i1 ) {
	+ for ( int i2 = 0; i2 < 2; ++i2 ) {
	+ const EvtVector4C hardCur = hardCurrPP( parent, i1, i2 );
	+ const EvtVector4C Heps = H.cont2( hardCur );
	+ for ( int i = 0; i < 3; i++ ) {
	+ // psi-meson polarization vector
	+ const EvtVector4C eps = Jpsi->epsParent( i ).conj();
	+ const EvtComplex amp = eps * Heps;
	+ vertex( i, i1, i2, amp );
	+ }
	+ }
	+ }
	+}
	diff --git a/src/EvtGenModels/EvtModelReg.cpp b/src/EvtGenModels/EvtModelReg.cpp
	--- a/src/EvtGenModels/EvtModelReg.cpp
	+++ b/src/EvtGenModels/EvtModelReg.cpp
	@@ -46,6 +46,7 @@
	#include "EvtGenModels/EvtBcVHad.hh"
	#include "EvtGenModels/EvtBcVMuNu.hh"
	#include "EvtGenModels/EvtBcVNpi.hh"
	+#include "EvtGenModels/EvtBcVPPHad.hh"
	#include "EvtGenModels/EvtBsMuMuKK.hh"
	#include "EvtGenModels/EvtBsquark.hh"
	#include "EvtGenModels/EvtBto2piCPiso.hh"
	@@ -297,7 +298,7 @@
	modelist.registerModel( new EvtPVVCPLH );
	modelist.registerModel( new EvtSSD_DirectCP );

	- modelist.registerModel( new EvtBcToNPi( true ) ); // true = print author info
	+ modelist.registerModel( new EvtBcToNPi );
	modelist.registerModel( new EvtBcPsiNPi );
	modelist.registerModel( new EvtBcBsNPi );
	modelist.registerModel( new EvtBcBsStarNPi );
	@@ -318,9 +319,10 @@
	modelist.registerModel( new EvtVtoSll );

	modelist.registerModel( new EvtBsMuMuKK );
	- modelist.registerModel( new EvtGenericDalitz() );
	+ modelist.registerModel( new EvtGenericDalitz );

	modelist.registerModel( new EvtBcVHad );
	+ modelist.registerModel( new EvtBcVPPHad );

	modelist.registerModel( new Evtbs2llGammaMNT );
	modelist.registerModel( new Evtbs2llGammaISRFSR );
	diff --git a/src/EvtGenModels/EvtPhiDalitz.cpp b/src/EvtGenModels/EvtPhiDalitz.cpp
	--- a/src/EvtGenModels/EvtPhiDalitz.cpp
	+++ b/src/EvtGenModels/EvtPhiDalitz.cpp
	@@ -25,6 +25,9 @@
	#include "EvtGenBase/EvtParticle.hh"
	#include "EvtGenBase/EvtPatches.hh"
	#include "EvtGenBase/EvtReport.hh"
	+#include "EvtGenBase/EvtVector3C.hh"
	+#include "EvtGenBase/EvtVector3R.hh"
	+#include "EvtGenBase/EvtVector4C.hh"
	#include "EvtGenBase/EvtVector4R.hh"

	#include <math.h>
	@@ -57,26 +60,27 @@
	checkSpinDaughter( 1, EvtSpinType::SCALAR );
	checkSpinDaughter( 2, EvtSpinType::SCALAR );

	- _mRho = 0.7758;
	- _gRho = 0.1439;
	- _aD = 0.78;
	- _phiD = -2.47;
	- _aOmega = 0.0071;
	- _phiOmega = -0.22;
	+ // results taken from KLOE results: arxiv.org/abs/hep-ex/0303016
	+ m_mRho = 0.7758;
	+ m_gRho = 0.1439;
	+ m_aD = 0.78;
	+ m_phiD = -2.47;
	+ m_aOmega = 0.0071;
	+ m_phiOmega = -0.22;

	- _locPip = -1;
	- _locPim = -1;
	- _locPi0 = -1;
	+ m_locPip = -1;
	+ m_locPim = -1;
	+ m_locPi0 = -1;

	for ( int i = 0; i < 3; i++ ) {
	if ( getDaug( i ) == EvtPDL::getId( "pi+" ) )
	- _locPip = i;
	+ m_locPip = i;
	if ( getDaug( i ) == EvtPDL::getId( "pi-" ) )
	- _locPim = i;
	+ m_locPim = i;
	if ( getDaug( i ) == EvtPDL::getId( "pi0" ) )
	- _locPi0 = i;
	+ m_locPi0 = i;
	}
	- if ( _locPip == -1 \|\| _locPim == -1 \|\| _locPi0 == -1 ) {
	+ if ( m_locPip == -1 \|\| m_locPim == -1 \|\| m_locPi0 == -1 ) {
	EvtGenReport( EVTGEN_ERROR, "EvtGen" )
	<< getModelName()
	<< "generator expects daughters to be pi+ pi- pi0\n";
	@@ -89,80 +93,118 @@

	void EvtPhiDalitz::initProbMax()
	{
	- setProbMax( 150.0 );
	+ setProbMax( 300.0 );
	}

	void EvtPhiDalitz::decay( EvtParticle* p )
	{
	- EvtId PIP = EvtPDL::getId( "pi+" );
	- EvtId PIM = EvtPDL::getId( "pi-" );
	- EvtId PIZ = EvtPDL::getId( "pi0" );
	- EvtId OMEGA = EvtPDL::getId( "omega" );
	+ const EvtId PIP = EvtPDL::getId( "pi+" );
	+ const EvtId PIM = EvtPDL::getId( "pi-" );
	+ const EvtId PIZ = EvtPDL::getId( "pi0" );
	+ const EvtId OMEGA = EvtPDL::getId( "omega" );

	p->initializePhaseSpace( getNDaug(), getDaugs() );

	- EvtVector4R Ppip = p->getDaug( _locPip )->getP4();
	- EvtVector4R Ppim = p->getDaug( _locPim )->getP4();
	- EvtVector4R Ppi0 = p->getDaug( _locPi0 )->getP4();
	- EvtVector4R Qp = ( Ppim + Ppi0 );
	- EvtVector4R Qm = ( Ppip + Ppi0 );
	- EvtVector4R Q0 = ( Ppip + Ppim );
	- double m2_pip = pow( EvtPDL::getMeanMass( PIP ), 2 );
	- double m2_pim = pow( EvtPDL::getMeanMass( PIM ), 2 );
	- double m2_pi0 = pow( EvtPDL::getMeanMass( PIZ ), 2 );
	- double M2rhop = pow( _mRho, 2 );
	- double M2rhom = pow( _mRho, 2 );
	- double M2rho0 = pow( _mRho, 2 );
	- double M2omega = pow( EvtPDL::getMeanMass( OMEGA ), 2 );
	-
	- double Wrhop = _gRho;
	- double Wrhom = _gRho;
	- double Wrho0 = _gRho;
	- double Womega = EvtPDL::getWidth( OMEGA );
	-
	- EvtComplex Atot( 0, 0 );
	-
	- //Rho+ Risonance Amplitude
	- double Gp = Wrhop *
	- pow( ( ( Qp.mass2() - m2_pim - m2_pi0 ) / 2 - M2rhop / 4 ) /
	- ( M2rhop / 4 - ( m2_pim + m2_pi0 ) / 2 ),
	- 3 / 2 ) *
	- ( M2rhop / Qp.mass2() );
	- EvtComplex Drhop( ( Qp.mass2() - M2rhop ), Qp.mass() * Gp );
	- EvtComplex A1( M2rhop / Drhop );
	-
	- //Rho- Risonance Amplitude
	- double Gm = Wrhom *
	- pow( ( ( Qm.mass2() - m2_pip - m2_pi0 ) / 2 - M2rhom / 4 ) /
	- ( M2rhom / 4 - ( m2_pip + m2_pi0 ) / 2 ),
	- 3 / 2 ) *
	- ( M2rhom / Qm.mass2() );
	- EvtComplex Drhom( ( Qm.mass2() - M2rhom ), Qm.mass() * Gm );
	- EvtComplex A2( M2rhom / Drhom );
	-
	- //Rho0 Risonance Amplitude
	- double G0 = Wrho0 *
	- pow( ( ( Q0.mass2() - m2_pip - m2_pim ) / 2 - M2rho0 / 4 ) /
	- ( M2rho0 / 4 - ( m2_pip + m2_pim ) / 2 ),
	- 3 / 2 ) *
	- ( M2rho0 / Q0.mass2() );
	- EvtComplex Drho0( ( Q0.mass2() - M2rho0 ), Q0.mass() * G0 );
	- EvtComplex A3( M2rho0 / Drho0 );
	-
	- //Omega Risonance Amplitude
	- EvtComplex OmegaPhase( 0, _phiOmega );
	- EvtComplex DOmega( ( Q0.mass2() - M2omega ), Q0.mass() * Womega );
	- EvtComplex A4( _aOmega * M2omega * exp( OmegaPhase ) / DOmega );
	+ const EvtVector4R Ppip = p->getDaug( m_locPip )->getP4();
	+ const EvtVector4R Ppim = p->getDaug( m_locPim )->getP4();
	+ const EvtVector4R Ppi0 = p->getDaug( m_locPi0 )->getP4();
	+ const EvtVector4R Qm = ( Ppim + Ppi0 );
	+ const EvtVector4R Qp = ( Ppip + Ppi0 );
	+ const EvtVector4R Q0 = ( Ppip + Ppim );
	+ const double m_pip = EvtPDL::getMeanMass( PIP );
	+ const double m_pim = EvtPDL::getMeanMass( PIM );
	+ const double m_pi0 = EvtPDL::getMeanMass( PIZ );
	+ const double Mrhop = m_mRho;
	+ const double Mrhom = m_mRho;
	+ const double Mrho0 = m_mRho;
	+ const double M2rhop = pow( Mrhop, 2 );
	+ const double M2rhom = pow( Mrhom, 2 );
	+ const double M2rho0 = pow( Mrho0, 2 );
	+ const double M2omega = pow( EvtPDL::getMeanMass( OMEGA ), 2 );
	+
	+ const double Wrhop = m_gRho;
	+ const double Wrhom = m_gRho;
	+ const double Wrho0 = m_gRho;
	+ const double Womega = EvtPDL::getWidth( OMEGA );
	+
	+ const double QmM = Qm.mass();
	+ const double QmM2 = Qm.mass2();
	+ const double QpM = Qp.mass();
	+ const double QpM2 = Qp.mass2();
	+ const double Q0M = Q0.mass();
	+ const double Q0M2 = Q0.mass2();
	+
	+ //Rho- Resonance Amplitude
	+ const double qm = calc_q( QmM, m_pim, m_pi0 );
	+ const double qm_0 = calc_q( Mrhom, m_pim, m_pi0 );
	+ const double Gm = Wrhom * pow( qm / qm_0, 3 ) * ( M2rhom / QmM2 );
	+ const EvtComplex Drhom( ( QmM2 - M2rhom ), QmM * Gm );
	+ const EvtComplex A1( M2rhom / Drhom );
	+
	+ //Rho+ Resonance Amplitude
	+ const double qp = calc_q( QpM, m_pip, m_pi0 );
	+ const double qp_0 = calc_q( Mrhop, m_pip, m_pi0 );
	+ const double Gp = Wrhop * pow( qp / qp_0, 3 ) * ( M2rhop / QpM2 );
	+ const EvtComplex Drhop( ( QpM2 - M2rhop ), QpM * Gp );
	+ const EvtComplex A2( M2rhop / Drhop );
	+
	+ //Rho0 Resonance Amplitude
	+ const double q0 = calc_q( Q0M, m_pip, m_pim );
	+ const double q0_0 = calc_q( Mrho0, m_pip, m_pim );
	+ const double G0 = Wrho0 * pow( q0 / q0_0, 3 ) * ( M2rho0 / Q0M2 );
	+ const EvtComplex Drho0( ( Q0M2 - M2rho0 ), Q0M * G0 );
	+ const EvtComplex A3( M2rho0 / Drho0 );
	+
	+ //Omega Resonance Amplitude
	+ const EvtComplex OmegaA( m_aOmega * cos( m_phiOmega ),
	+ m_aOmega * sin( m_phiOmega ) );
	+ const EvtComplex DOmega( ( Q0M2 - M2omega ), Q0M * Womega );
	+ const EvtComplex A4( OmegaA * M2omega / DOmega );

	//Direct Decay Amplitude
	- EvtComplex DirPhase( 0, _phiD );
	- EvtComplex A5( _aD * exp( DirPhase ) );
	+ const EvtComplex A5( m_aD * cos( m_phiD ), m_aD * sin( m_phiD ) );

	- Atot = A1 + A2 + A3 + A4 + A5;
	+ const EvtComplex Atot = A1 + A2 + A3 + A4 + A5;

	- vertex( 0, Atot );
	- vertex( 1, Atot );
	- vertex( 2, Atot );
	+ // Polarization
	+
	+ const EvtVector4C ep0 = p->eps( 0 );
	+ const EvtVector4C ep1 = p->eps( 1 );
	+ const EvtVector4C ep2 = p->eps( 2 );
	+
	+ const EvtVector3R p1( Ppip.get( 1 ), Ppip.get( 2 ), Ppip.get( 3 ) );
	+ const EvtVector3R p2( Ppim.get( 1 ), Ppim.get( 2 ), Ppim.get( 3 ) );
	+ const EvtVector3R q = cross( p1, p2 );
	+
	+ const EvtVector3C e1( ep0.get( 1 ), ep0.get( 2 ), ep0.get( 3 ) );
	+ const EvtVector3C e2( ep1.get( 1 ), ep1.get( 2 ), ep1.get( 3 ) );
	+ const EvtVector3C e3( ep2.get( 1 ), ep2.get( 2 ), ep2.get( 3 ) );
	+
	+ //This is an approximate formula of the maximum value that
	+ //\|q\| can have.
	+ const double pM = p->mass();
	+ const double mSum = Ppip.mass() + Ppim.mass() + Ppi0.mass();
	+ const double norm = 10.26 / ( pM * pM - mSum * mSum );
	+
	+ vertex( 0, norm * e1 * q * Atot );
	+ vertex( 1, norm * e2 * q * Atot );
	+ vertex( 2, norm * e3 * q * Atot );

	return;
	}
	+
	+double EvtPhiDalitz::calc_q( double M, double m1, double m2 ) const
	+{
	+ const double m12Sum = m1 + m2;
	+
	+ if ( M > m12Sum ) {
	+ const double MSq = M * M;
	+ const double m12Diff = m1 - m2;
	+
	+ return sqrt( ( MSq - ( m12Sum ) * ( m12Sum ) ) *
	+ ( MSq - ( m12Diff ) * ( m12Diff ) ) ) /
	+ ( 2.0 * M );
	+ } else {
	+ return 0.0;
	+ }
	+}
	diff --git a/src/EvtGenModels/EvtWHad.cpp b/src/EvtGenModels/EvtWHad.cpp
	--- a/src/EvtGenModels/EvtWHad.cpp
	+++ b/src/EvtGenModels/EvtWHad.cpp
	@@ -276,7 +276,7 @@
	p6 ); // pi+ pi+ pi+ pi- pi-
	const EvtVector4R pf0 = p4 + p7;
	return eps1 * BWa( qTot ) * BWf( pf0 );
	-};
	+}

	// hadronic current W+ -> K+ pi+ pi-

	@@ -457,3 +457,59 @@
	dual( EvtGenFunctions::directProd( epsKstar, epsA1 ) ).cont2( pKstar - pa1 );
	return eps;
	}
	+
	+EvtVector4C EvtWHad::WCurrent_ppPi( const EvtVector4R& p1,
	+ const EvtDiracSpinor& sp1,
	+ const EvtVector4R& p2,
	+ const EvtDiracSpinor& sp2,
	+ const EvtVector4R& k ) const
	+{
	+ const EvtVector4R q = p1 + p2 + k;
	+ const double q2 = q.mass2();
	+ const double mp = p1.mass(), mp2 = mp * mp, mpi = k.mass(), mpi2 = mpi * mpi;
	+ const double mn = EvtPDL::getMeanMass( EvtPDL::getId( "n0" ) );
	+ const double mn2 = mn * mn;
	+ const EvtComplex II( 0, 1 );
	+
	+ const double kp2 = k * p2;
	+ const double p1p2 = p1 * p2;
	+
	+ const double f1 = 1.0;
	+ const double f2 = 3.7 / ( 2.0 + mp );
	+ const double g1 = 1.25;
	+ const double g3 = 2.0 * mp * g1 / ( p1p2 + mp2 );
	+
	+ const EvtComplex curS = EvtLeptonSCurrent( sp1, sp2 );
	+ const EvtComplex curP = EvtLeptonPCurrent( sp1, sp2 );
	+ const EvtVector4C curV = EvtLeptonVCurrent( sp1, sp2 );
	+ const EvtVector4C curA = EvtLeptonACurrent( sp1, sp2 );
	+ const EvtTensor4C curT = EvtLeptonTCurrent( sp1, sp2 );
	+
	+ const double D1 = 1 / ( 2 * kp2 + mp2 - mn2 + mpi2 ); //(k+p2)^2-mn^2
	+
	+ // Amplitude: ~U(p1).GA5.Vpp(alpha).prop(-p2-k).V(p2) + permutations
	+ // U() and V() are proton and antiproton spinors, prop() is the proton's propagator,
	+ // and Vpp(alpha) is the W->pp vertex.
	+ // Expand terms to use basic spinor currents (Scalar, Vector, Axial, etc)
	+
	+ EvtVector4C current;
	+ current += curA * ( -( f2 * II ) );
	+ current += ( D1 * g1 * II ) * curT.cont2( k );
	+ current += -( -0.5 * ( D1 * f1 ) ) * dual( curT ).cont2( k );
	+ current += -( D1 * f2 * II * mp ) * dual( curT ).cont2( k );
	+ current += k * ( -( D1 * g1 ) ) * curS;
	+ current += k * ( -( D1 * f1 ) ) * curP;
	+ current += k * ( 2 * D1 * f2 * II * mp ) * curP;
	+ current += k * ( curV * k ) * ( D1 * g3 );
	+ current += k * ( curA * k ) * ( D1 * f2 * II );
	+ current += p1 * ( curV * k ) * ( D1 * g3 );
	+ current += p1 * ( curA * k ) * ( -( D1 * f2 * II ) );
	+ current += p2 * ( curV * k ) * ( D1 * g3 );
	+ current += p2 * ( curA * k ) * ( D1 * f2 * II );
	+
	+ const EvtTensor4C JT = ( 1.0 / q2 ) * EvtGenFunctions::directProd( q, q ) -
	+ EvtTensor4C::g();
	+ current = JT.cont2( current );
	+
	+ return BWa( q ) * current;
	+}
	diff --git a/test/jsonFiles/BC_VPPHAD1=VLL__Bc_Jpsippbarpi_mumu.json b/test/jsonFiles/BC_VPPHAD1=VLL__Bc_Jpsippbarpi_mumu.json
	new file mode 100644
	--- /dev/null
	+++ b/test/jsonFiles/BC_VPPHAD1=VLL__Bc_Jpsippbarpi_mumu.json
	@@ -0,0 +1,38 @@
	+{
	+ "parent" : "B_c+",
	+ "daughters" : ["J/psi", "p+", "anti-p-", "pi+"],
	+ "grand_daughters" : [["mu+", "mu-"], [], [], []],
	+ "models" : ["BC_VPPHAD", "VLL", "", "", ""],
	+ "parameters" : [["1"], [], [], [], []],
	+ "extras" : ["noPhotos"],
	+ "events" : 10000,
	+ "histograms" : [
	+ {"variable" : "prob", "title" : "Prob(BC_VPPHAD)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 1.0},
	+ {"variable" : "prob_daug1", "title" : "Prob(VLL)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 1.0},
	+ {"variable" : "mass", "title" : "m(Jpsi)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 3.0, "xmax" : 3.2},
	+ {"variable" : "mass", "title" : "m(Jpsi p)", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 3.8, "xmax" : 5.4},
	+ {"variable" : "mass", "title" : "m(Jpsi pbar)", "d1" : 1, "d2" : 3, "nbins" : 100, "xmin" : 3.8, "xmax" : 5.4},
	+ {"variable" : "mass", "title" : "m(Jpsi pi)", "d1" : 1, "d2" : 4, "nbins" : 100, "xmin" : 3.0, "xmax" : 4.5},
	+ {"variable" : "mass", "title" : "m(p pbar)", "d1" : 2, "d2" : 3, "nbins" : 100, "xmin" : 1.8, "xmax" : 3.2},
	+ {"variable" : "mass", "title" : "m(p pi)", "d1" : 2, "d2" : 4, "nbins" : 100, "xmin" : 1.0, "xmax" : 2.4},
	+ {"variable" : "mass", "title" : "m(pbar pi)", "d1" : 3, "d2" : 4, "nbins" : 100, "xmin" : 1.0, "xmax" : 2.4},
	+ {"variable" : "cosHel", "title" : "cosHel(Jpsi, p)", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosHel", "title" : "cosHel(Jpsi, pbar)", "d1" : 1, "d2" : 3, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosHel", "title" : "cosHel(Jpsi, pi)", "d1" : 1, "d2" : 4, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosHel", "title" : "cosHel(p, pbar)", "d1" : 2, "d2" : 3, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosHel", "title" : "cosHel(p, pi)", "d1" : 2, "d2" : 4, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosHel", "title" : "cosHel(pbar, pi)", "d1" : 3, "d2" : 4, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosHel_daug1", "title" : "cosHel(mu+,Jpsi)", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosHel_daug1", "title" : "cosHel(mu-,Jpsi)", "d1" : 2, "d2" : 1, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosTheta", "title" : "cosTheta(Jpsi)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "phi", "title" : "phi(Jpsi)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -180.0, "xmax" : 180.0},
	+ {"variable" : "pLab_daug1", "title" : "pLab(mu+)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 0.5, "xmax" : 3.5},
	+ {"variable" : "pLab_daug1", "title" : "pLab(mu-)", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 0.5, "xmax" : 3.5},
	+ {"variable" : "cosTheta_daug1", "title" : "cosTheta(mu+)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "phi_daug1", "title" : "phi(mu+)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -180.0, "xmax" : 180.0},
	+ {"variable" : "cosTheta_daug1", "title" : "cosTheta(mu-)", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "phi_daug1", "title" : "phi(mu-)", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -180.0, "xmax" : 180.0}
	+ ],
	+ "outfile" : "BC_VPPHAD1=VLL__Bc_Jpsippbarpi_mumu.root",
	+ "reference" : "RefBC_VPPHAD1=VLL__Bc_Jpsippbarpi_mumu.root"
	+}
	diff --git a/test/jsonFiles/BC_VPPHAD1=VLL__Bc_psi2Sppbarpi_mumu.json b/test/jsonFiles/BC_VPPHAD1=VLL__Bc_psi2Sppbarpi_mumu.json
	new file mode 100644
	--- /dev/null
	+++ b/test/jsonFiles/BC_VPPHAD1=VLL__Bc_psi2Sppbarpi_mumu.json
	@@ -0,0 +1,38 @@
	+{
	+ "parent" : "B_c+",
	+ "daughters" : ["psi(2S)", "p+", "anti-p-", "pi+"],
	+ "grand_daughters" : [["mu+", "mu-"], [], [], []],
	+ "models" : ["BC_VPPHAD", "VLL", "", "", ""],
	+ "parameters" : [["1"], [], [], [], []],
	+ "extras" : ["noPhotos"],
	+ "events" : 10000,
	+ "histograms" : [
	+ {"variable" : "prob", "title" : "Prob(BC_VPPHAD)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 1.0},
	+ {"variable" : "prob_daug1", "title" : "Prob(VLL)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 1.0},
	+ {"variable" : "mass", "title" : "m(psi2S)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 3.684, "xmax" : 3.69},
	+ {"variable" : "mass", "title" : "m(psi2S p)", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 4.5, "xmax" : 5.3},
	+ {"variable" : "mass", "title" : "m(psi2S pbar)", "d1" : 1, "d2" : 3, "nbins" : 100, "xmin" : 4.5, "xmax" : 5.3},
	+ {"variable" : "mass", "title" : "m(psi2S pi)", "d1" : 1, "d2" : 4, "nbins" : 100, "xmin" : 3.7, "xmax" : 4.5},
	+ {"variable" : "mass", "title" : "m(p pbar)", "d1" : 2, "d2" : 3, "nbins" : 100, "xmin" : 1.8, "xmax" : 2.5},
	+ {"variable" : "mass", "title" : "m(p pi)", "d1" : 2, "d2" : 4, "nbins" : 100, "xmin" : 1.0, "xmax" : 1.7},
	+ {"variable" : "mass", "title" : "m(pbar pi)", "d1" : 3, "d2" : 4, "nbins" : 100, "xmin" : 1.0, "xmax" : 1.7},
	+ {"variable" : "cosHel", "title" : "cosHel(psi2S, p)", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosHel", "title" : "cosHel(psi2S, pbar)", "d1" : 1, "d2" : 3, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosHel", "title" : "cosHel(psi2S, pi)", "d1" : 1, "d2" : 4, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosHel", "title" : "cosHel(p, pbar)", "d1" : 2, "d2" : 3, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosHel", "title" : "cosHel(p, pi)", "d1" : 2, "d2" : 4, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosHel", "title" : "cosHel(pbar, pi)", "d1" : 3, "d2" : 4, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosHel_daug1", "title" : "cosHel(mu+,psi2S)", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosHel_daug1", "title" : "cosHel(mu-,psi2S)", "d1" : 2, "d2" : 1, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosTheta", "title" : "cosTheta(psi2S)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "phi", "title" : "phi(psi2S)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -180.0, "xmax" : 180.0},
	+ {"variable" : "pLab_daug1", "title" : "pLab(mu+)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 0.5, "xmax" : 3.5},
	+ {"variable" : "pLab_daug1", "title" : "pLab(mu-)", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 0.5, "xmax" : 3.5},
	+ {"variable" : "cosTheta_daug1", "title" : "cosTheta(mu+)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "phi_daug1", "title" : "phi(mu+)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -180.0, "xmax" : 180.0},
	+ {"variable" : "cosTheta_daug1", "title" : "cosTheta(mu-)", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "phi_daug1", "title" : "phi(mu-)", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -180.0, "xmax" : 180.0}
	+ ],
	+ "outfile" : "BC_VPPHAD1=VLL__Bc_psi2Sppbarpi_mumu.root",
	+ "reference" : "RefBC_VPPHAD1=VLL__Bc_psi2Sppbarpi_mumu.root"
	+}
	diff --git a/test/jsonFiles/BC_VPPHAD2=VLL__Bc_Jpsippbarpi_mumu.json b/test/jsonFiles/BC_VPPHAD2=VLL__Bc_Jpsippbarpi_mumu.json
	new file mode 100644
	--- /dev/null
	+++ b/test/jsonFiles/BC_VPPHAD2=VLL__Bc_Jpsippbarpi_mumu.json
	@@ -0,0 +1,38 @@
	+{
	+ "parent" : "B_c+",
	+ "daughters" : ["J/psi", "p+", "anti-p-", "pi+"],
	+ "grand_daughters" : [["mu+", "mu-"], [], [], []],
	+ "models" : ["BC_VPPHAD", "VLL", "", "", ""],
	+ "parameters" : [["2"], [], [], [], []],
	+ "extras" : ["noPhotos"],
	+ "events" : 10000,
	+ "histograms" : [
	+ {"variable" : "prob", "title" : "Prob(BC_VPPHAD)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 1.0},
	+ {"variable" : "prob_daug1", "title" : "Prob(VLL)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 1.0},
	+ {"variable" : "mass", "title" : "m(Jpsi)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 3.0, "xmax" : 3.2},
	+ {"variable" : "mass", "title" : "m(Jpsi p)", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 3.8, "xmax" : 5.4},
	+ {"variable" : "mass", "title" : "m(Jpsi pbar)", "d1" : 1, "d2" : 3, "nbins" : 100, "xmin" : 3.8, "xmax" : 5.4},
	+ {"variable" : "mass", "title" : "m(Jpsi pi)", "d1" : 1, "d2" : 4, "nbins" : 100, "xmin" : 3.0, "xmax" : 4.5},
	+ {"variable" : "mass", "title" : "m(p pbar)", "d1" : 2, "d2" : 3, "nbins" : 100, "xmin" : 1.8, "xmax" : 3.2},
	+ {"variable" : "mass", "title" : "m(p pi)", "d1" : 2, "d2" : 4, "nbins" : 100, "xmin" : 1.0, "xmax" : 2.4},
	+ {"variable" : "mass", "title" : "m(pbar pi)", "d1" : 3, "d2" : 4, "nbins" : 100, "xmin" : 1.0, "xmax" : 2.4},
	+ {"variable" : "cosHel", "title" : "cosHel(Jpsi, p)", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosHel", "title" : "cosHel(Jpsi, pbar)", "d1" : 1, "d2" : 3, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosHel", "title" : "cosHel(Jpsi, pi)", "d1" : 1, "d2" : 4, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosHel", "title" : "cosHel(p, pbar)", "d1" : 2, "d2" : 3, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosHel", "title" : "cosHel(p, pi)", "d1" : 2, "d2" : 4, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosHel", "title" : "cosHel(pbar, pi)", "d1" : 3, "d2" : 4, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosHel_daug1", "title" : "cosHel(mu+,Jpsi)", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosHel_daug1", "title" : "cosHel(mu-,Jpsi)", "d1" : 2, "d2" : 1, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosTheta", "title" : "cosTheta(Jpsi)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "phi", "title" : "phi(Jpsi)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -180.0, "xmax" : 180.0},
	+ {"variable" : "pLab_daug1", "title" : "pLab(mu+)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 0.5, "xmax" : 3.5},
	+ {"variable" : "pLab_daug1", "title" : "pLab(mu-)", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 0.5, "xmax" : 3.5},
	+ {"variable" : "cosTheta_daug1", "title" : "cosTheta(mu+)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "phi_daug1", "title" : "phi(mu+)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -180.0, "xmax" : 180.0},
	+ {"variable" : "cosTheta_daug1", "title" : "cosTheta(mu-)", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "phi_daug1", "title" : "phi(mu-)", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -180.0, "xmax" : 180.0}
	+ ],
	+ "outfile" : "BC_VPPHAD2=VLL__Bc_Jpsippbarpi_mumu.root",
	+ "reference" : "RefBC_VPPHAD2=VLL__Bc_Jpsippbarpi_mumu.root"
	+}
	diff --git a/test/jsonFiles/BC_VPPHAD2=VLL__Bc_psi2Sppbarpi_mumu.json b/test/jsonFiles/BC_VPPHAD2=VLL__Bc_psi2Sppbarpi_mumu.json
	new file mode 100644
	--- /dev/null
	+++ b/test/jsonFiles/BC_VPPHAD2=VLL__Bc_psi2Sppbarpi_mumu.json
	@@ -0,0 +1,38 @@
	+{
	+ "parent" : "B_c+",
	+ "daughters" : ["psi(2S)", "p+", "anti-p-", "pi+"],
	+ "grand_daughters" : [["mu+", "mu-"], [], [], []],
	+ "models" : ["BC_VPPHAD", "VLL", "", "", ""],
	+ "parameters" : [["2"], [], [], [], []],
	+ "extras" : ["noPhotos"],
	+ "events" : 10000,
	+ "histograms" : [
	+ {"variable" : "prob", "title" : "Prob(BC_VPPHAD)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 1.0},
	+ {"variable" : "prob_daug1", "title" : "Prob(VLL)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 1.0},
	+ {"variable" : "mass", "title" : "m(psi2S)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 3.684, "xmax" : 3.69},
	+ {"variable" : "mass", "title" : "m(psi2S p)", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 4.5, "xmax" : 5.3},
	+ {"variable" : "mass", "title" : "m(psi2S pbar)", "d1" : 1, "d2" : 3, "nbins" : 100, "xmin" : 4.5, "xmax" : 5.3},
	+ {"variable" : "mass", "title" : "m(psi2S pi)", "d1" : 1, "d2" : 4, "nbins" : 100, "xmin" : 3.7, "xmax" : 4.5},
	+ {"variable" : "mass", "title" : "m(p pbar)", "d1" : 2, "d2" : 3, "nbins" : 100, "xmin" : 1.8, "xmax" : 2.5},
	+ {"variable" : "mass", "title" : "m(p pi)", "d1" : 2, "d2" : 4, "nbins" : 100, "xmin" : 1.0, "xmax" : 1.7},
	+ {"variable" : "mass", "title" : "m(pbar pi)", "d1" : 3, "d2" : 4, "nbins" : 100, "xmin" : 1.0, "xmax" : 1.7},
	+ {"variable" : "cosHel", "title" : "cosHel(psi2S, p)", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosHel", "title" : "cosHel(psi2S, pbar)", "d1" : 1, "d2" : 3, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosHel", "title" : "cosHel(psi2S, pi)", "d1" : 1, "d2" : 4, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosHel", "title" : "cosHel(p, pbar)", "d1" : 2, "d2" : 3, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosHel", "title" : "cosHel(p, pi)", "d1" : 2, "d2" : 4, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosHel", "title" : "cosHel(pbar, pi)", "d1" : 3, "d2" : 4, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosHel_daug1", "title" : "cosHel(mu+,psi2S)", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosHel_daug1", "title" : "cosHel(mu-,psi2S)", "d1" : 2, "d2" : 1, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosTheta", "title" : "cosTheta(psi2S)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "phi", "title" : "phi(psi2S)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -180.0, "xmax" : 180.0},
	+ {"variable" : "pLab_daug1", "title" : "pLab(mu+)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 0.5, "xmax" : 3.5},
	+ {"variable" : "pLab_daug1", "title" : "pLab(mu-)", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 0.5, "xmax" : 3.5},
	+ {"variable" : "cosTheta_daug1", "title" : "cosTheta(mu+)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "phi_daug1", "title" : "phi(mu+)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -180.0, "xmax" : 180.0},
	+ {"variable" : "cosTheta_daug1", "title" : "cosTheta(mu-)", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "phi_daug1", "title" : "phi(mu-)", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -180.0, "xmax" : 180.0}
	+ ],
	+ "outfile" : "BC_VPPHAD2=VLL__Bc_psi2Sppbarpi_mumu.root",
	+ "reference" : "RefBC_VPPHAD2=VLL__Bc_psi2Sppbarpi_mumu.root"
	+}
	diff --git a/test/jsonFiles/PHI_DALITZ__phi_pipipi0.json b/test/jsonFiles/PHI_DALITZ__phi_pipipi0.json
	--- a/test/jsonFiles/PHI_DALITZ__phi_pipipi0.json
	+++ b/test/jsonFiles/PHI_DALITZ__phi_pipipi0.json
	@@ -15,7 +15,13 @@
	{"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" : "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}
	+ {"variable" : "cosHel", "title" : "cosHel13", "d1" : 1, "d2" : 3, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosTheta", "title" : "cosTheta1", "d1" : 1, "d2" : 0, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosTheta", "title" : "cosTheta2", "d1" : 2, "d2" : 0, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosTheta", "title" : "cosTheta3", "d1" : 3, "d2" : 0, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "phi", "title" : "phi1", "d1" : 1, "d2" : 0, "nbins" : 50, "xmin" : -180.0, "xmax" : 180.0},
	+ {"variable" : "phi", "title" : "phi2", "d1" : 2, "d2" : 0, "nbins" : 50, "xmin" : -180.0, "xmax" : 180.0},
	+ {"variable" : "phi", "title" : "phi3", "d1" : 3, "d2" : 0, "nbins" : 50, "xmin" : -180.0, "xmax" : 180.0}
	],
	"outfile" : "PHI_DALITZ__phi_pipipi0.root",
	"reference" : "RefPHI_DALITZ__phi_pipipi0.root"
	diff --git a/test/jsonFiles/SVS_PHI_DALITZ__Ds_phipi_pipipi0.json b/test/jsonFiles/SVS_PHI_DALITZ__Ds_phipi_pipipi0.json
	new file mode 100644
	--- /dev/null
	+++ b/test/jsonFiles/SVS_PHI_DALITZ__Ds_phipi_pipipi0.json
	@@ -0,0 +1,38 @@
	+{
	+ "parent" : "D_s+",
	+ "daughters" : ["phi", "pi+"],
	+ "grand_daughters" : [["pi+", "pi-", "pi0"], []],
	+ "models" : ["SVS", "PHI_DALITZ", ""],
	+ "parameters" : [[], [], []],
	+ "events" : 10000,
	+ "histograms" : [
	+ {"variable" : "prob", "title" : "Prob(SVS)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 1.0},
	+ {"variable" : "prob_daug1", "title" : "Prob(PHI_DALITZ)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 1.0},
	+ {"variable" : "parMass", "title" : "m(D_s^{+})", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 1.95, "xmax" : 2.0},
	+ {"variable" : "mass", "title" : "m(#phi)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 0.2, "xmax" : 1.0},
	+ {"variable" : "mass_daug1", "title" : "m(#pi^{+} #pi^{-}) from #phi", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 0.2, "xmax" : 1.0},
	+ {"variable" : "mass_daug1", "title" : "m(#pi^{+} #pi^{0}) from #phi", "d1" : 1, "d2" : 3, "nbins" : 100, "xmin" : 0.2, "xmax" : 1.0},
	+ {"variable" : "mass_daug1", "title" : "m(#pi^{-} #pi^{0}) from #phi", "d1" : 2, "d2" : 3, "nbins" : 100, "xmin" : 0.2, "xmax" : 1.0},
	+ {"variable" : "massSq_daug1", "title" : "Dalitz m^{2}(#pi^{+} #pi^{-}) vs m^{2}(#pi^{+} #pi^{0})", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 0.8, "variableY" : "massSq_daug1", "d1Y" : 1, "d2Y" : 3, "nbinsY" : 50, "ymin" : 0.0, "ymax" : 0.8},
	+ {"variable" : "massSq_daug1", "title" : "Dalitz m^{2}(#pi^{-} #pi^{0}) vs m^{2}(#pi^{+} #pi^{-})", "d1" : 2, "d2" : 3, "nbins" : 50, "xmin" : 0.0, "xmax" : 0.8, "variableY" : "massSq_daug1", "d1Y" : 1, "d2Y" : 2, "nbinsY" : 50, "ymin" : 0.0, "ymax" : 0.8},
	+ {"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" : "cosThetaResNorm", "title" : "cosThetaResNorm", "d1" : 0, "d2" : 0, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosBetaRes_daug1", "title" : "cosBetaRes_daug1", "d1" : 0, "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_daug1", "title" : "cosHel1_12", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosHel_daug1", "title" : "cosHel1_23", "d1" : 2, "d2" : 3, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosHel_daug1", "title" : "cosHel1_13", "d1" : 1, "d2" : 3, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosTheta", "title" : "cosTheta1", "d1" : 1, "d2" : 0, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosTheta", "title" : "cosTheta1", "d1" : 2, "d2" : 0, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosTheta_daug1", "title" : "cosTheta1_1", "d1" : 1, "d2" : 0, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosTheta_daug1", "title" : "cosTheta1_2", "d1" : 2, "d2" : 0, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "cosTheta_daug1", "title" : "cosTheta1_3", "d1" : 3, "d2" : 0, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0},
	+ {"variable" : "phi", "title" : "phi1", "d1" : 1, "d2" : 0, "nbins" : 50, "xmin" : -180.0, "xmax" : 180.0},
	+ {"variable" : "phi", "title" : "phi2", "d1" : 2, "d2" : 0, "nbins" : 50, "xmin" : -180.0, "xmax" : 180.0},
	+ {"variable" : "phi_daug1", "title" : "phi1_1", "d1" : 1, "d2" : 0, "nbins" : 50, "xmin" : -180.0, "xmax" : 180.0},
	+ {"variable" : "phi_daug1", "title" : "phi1_2", "d1" : 2, "d2" : 0, "nbins" : 50, "xmin" : -180.0, "xmax" : 180.0},
	+ {"variable" : "phi_daug1", "title" : "phi1_3", "d1" : 3, "d2" : 0, "nbins" : 50, "xmin" : -180.0, "xmax" : 180.0}
	+ ],
	+ "outfile" : "SVS_PHI_DALITZ__Ds_phipi_pipipi0.root",
	+ "reference" : "RefSVS_PHI_DALITZ__Ds_phipi_pipipi0.root"
	+}
	diff --git a/test/testDecayModel.cc b/test/testDecayModel.cc
	--- a/test/testDecayModel.cc
	+++ b/test/testDecayModel.cc
	@@ -976,6 +976,88 @@
	EvtDecayAngleChi( p_parent, p_daug1, p_daug2, p_daug3, p_daug4 ) };
	value = chi * 180.0 / EvtConst::pi;

	+ } else if ( !selectedVarName.compare( "cosThetaResNorm" ) ) {
	+ // P -> R p4 -> (p1 p2 p3) p4, where the resonance R decays to p1 p2 p3.
	+ // Theta is the angle between the normal of the plane containing p1, p2 & p3
	+ // and the bachelor particle p4 in the rest frame of resonance R.
	+ // The normal vector is given by the cross product p3 x p1
	+
	+ if ( sel_NDaugMax > 1 ) {
	+ const EvtParticle* res{ selectedParent->getDaug( 0 ) };
	+ const EvtParticle* bac{ selectedParent->getDaug( 1 ) };
	+
	+ // Check resonance has 3 daughters
	+ if ( res != nullptr && res->getNDaug() == 3 ) {
	+ const EvtParticle* daug1 = res->getDaug( 0 );
	+ const EvtParticle* daug3 = res->getDaug( 2 );
	+
	+ // 4-momenta in base parent P lab frame
	+ const EvtVector4R pRes{ res->getP4Lab() };
	+ const EvtVector4R p4{ bac != nullptr ? bac->getP4Lab()
	+ : EvtVector4R() };
	+ const EvtVector4R p1{ daug1 != nullptr ? daug1->getP4Lab()
	+ : EvtVector4R() };
	+ const EvtVector4R p3{ daug3 != nullptr ? daug3->getP4Lab()
	+ : EvtVector4R() };
	+
	+ // Boost 4-vector for resonance frame
	+ const EvtVector4R boost{ pRes.get( 0 ), -pRes.get( 1 ),
	+ -pRes.get( 2 ), -pRes.get( 3 ) };
	+
	+ // Momentum of p1 and p3 in resonance frame
	+ const EvtVector4R p1Res{ boostTo( p1, boost ) };
	+ const EvtVector4R p3Res{ boostTo( p3, boost ) };
	+
	+ // Plane normal vector (just uses 3-momentum components)
	+ const EvtVector4R norm{ p3Res.cross( p1Res ) };
	+
	+ // Momentum of p4 in resonance frame
	+ const EvtVector4R p4Res{ boostTo( p4, boost ) };
	+
	+ // Cosine of the angle between the normal and p4 in the resonance frame
	+ const double normMag{ norm.d3mag() };
	+ const double p4ResMag{ p4Res.d3mag() };
	+ if ( normMag > 0.0 && p4ResMag > 0.0 ) {
	+ value = norm.dot( p4Res ) / ( normMag * p4ResMag );
	+ }
	+ }
	+ }
	+
	+ } else if ( !selectedVarName.compare( "cosBetaRes" ) ) {
	+ // For resonance P (parent) -> p1 p2 p3, beta is the
	+ // angle between p1 & p3 in the (p1 + p2) rest frame
	+ if ( sel_NDaugMax > 2 ) {
	+ const EvtParticle* daug1 = selectedParent->getDaug( 0 );
	+ const EvtParticle* daug2 = selectedParent->getDaug( 1 );
	+ const EvtParticle* daug3 = selectedParent->getDaug( 2 );
	+
	+ // 4-momenta in base parent frame
	+ const EvtVector4R p1{ daug1 != nullptr ? daug1->getP4Lab()
	+ : EvtVector4R() };
	+ const EvtVector4R p2{ daug2 != nullptr ? daug2->getP4Lab()
	+ : EvtVector4R() };
	+ const EvtVector4R p3{ daug3 != nullptr ? daug3->getP4Lab()
	+ : EvtVector4R() };
	+
	+ // p1 + p2
	+ const EvtVector4R p12{ p1 + p2 };
	+
	+ // Boost 4-vector for p12 frame
	+ const EvtVector4R boost{ p12.get( 0 ), -p12.get( 1 ), -p12.get( 2 ),
	+ -p12.get( 3 ) };
	+
	+ // Momentum of p1 & p3 in p12 frame
	+ const EvtVector4R p1_12{ boostTo( p1, boost ) };
	+ const EvtVector4R p3_12{ boostTo( p3, boost ) };
	+
	+ // Cosine of angle between p1 & p3 in p12 frame
	+ const double p1_12Mag{ p1_12.d3mag() };
	+ const double p3_12Mag{ p3_12.d3mag() };
	+ if ( p1_12Mag > 0.0 && p3_12Mag > 0.0 ) {
	+ value = p1_12.dot( p3_12 ) / ( p1_12Mag * p3_12Mag );
	+ }
	+ }
	+
	} else if ( !selectedVarName.compare( "E" ) ) {
	// Energy of first daughter in lab frame
	value = p1_lab.get( 0 );
	@@ -1165,17 +1247,18 @@
	: daughter1->getP4Lab() + daughter2->getP4Lab() };

	const EvtVector4R daughter4Vector2{
	- sel_NDaugMax == 2 ? daughter2->getP4Lab()
	- : sel_NDaugMax == 3 ? daughter2->getP4Lab() + daughter3->getP4Lab()
	- : daughter3->getP4Lab() + daughter4->getP4Lab() };
	+ sel_NDaugMax == 2
	+ ? daughter2->getP4Lab()
	+ : sel_NDaugMax == 3 ? daughter2->getP4Lab() + daughter3->getP4Lab()
	+ : daughter3->getP4Lab() + daughter4->getP4Lab() };

	const EvtVector4R grandDaughter4Vector1{
	sel_NDaugMax <= 3 ? grandDaughter1->getP4Lab() : daughter1->getP4Lab() };

	const EvtVector4R grandDaughter4Vector2{
	- sel_NDaugMax == 2 ? grandDaughter2->getP4Lab()
	- : sel_NDaugMax == 3 ? daughter2->getP4Lab()
	- : daughter3->getP4Lab() };
	+ sel_NDaugMax == 2 ? grandDaughter2->getP4Lab()
	+ : sel_NDaugMax == 3 ? daughter2->getP4Lab()
	+ : daughter3->getP4Lab() };

	const EvtVector4R parentBoost{ parent4Vector.get( 0 ),
	-parent4Vector.get( 1 ),
	@@ -1193,15 +1276,15 @@
	-daughter4Vector2.get( 3 ) };

	// Boosting daughters to reference frame of the mother
	- EvtVector4R daughter4Vector1_boosted{
	+ const EvtVector4R daughter4Vector1_boosted{
	boostTo( daughter4Vector1, parentBoost ) };
	- EvtVector4R daughter4Vector2_boosted{
	+ const EvtVector4R daughter4Vector2_boosted{
	boostTo( daughter4Vector2, parentBoost ) };

	// Boosting each granddaughter to reference frame of its mother
	- EvtVector4R grandDaughter4Vector1_boosted{
	+ const EvtVector4R grandDaughter4Vector1_boosted{
	boostTo( grandDaughter4Vector1, daughter1Boost ) };
	- EvtVector4R grandDaughter4Vector2_boosted{
	+ const EvtVector4R grandDaughter4Vector2_boosted{
	boostTo( grandDaughter4Vector2, daughter2Boost ) };

	// We calculate the normal vectors of the decay two planes

File Metadata

Mime Type: text/plain
Expires: Mon, Sep 29, 4:56 PM (18 h, 33 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 6545124
Default Alt Text: D100.1759161376.diff (84 KB)

D100.1759161376.diffNo OneTemporaryActions

D100.1759161376.diffView Options

File Metadata

Event Timeline

D100.1759161376.diff
No OneTemporary
Actions

D100.1759161376.diff
View Options