Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F19245247
D100.1759117826.diff
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Award Token
Flag For Later
Size
71 KB
Referenced Files
None
Subscribers
None
D100.1759117826.diff
View Options
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,6 +11,12 @@
===
## R02-0X-00
+10 Oct 2023 John Back
+* T230: Add EvtBcVPPHad model for Bc to Jpsi p pbar pi decays and generalised the 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).
+
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.
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/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
@@ -14,7 +14,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"
File Metadata
Details
Attached
Mime Type
text/plain
Expires
Mon, Sep 29, 4:50 AM (9 h, 54 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
6564837
Default Alt Text
D100.1759117826.diff (71 KB)
Attached To
Mode
D100: Add EvtBcVPPHad model and updated EvtPhiDalitz to use helicity amplitudes.
Attached
Detach File
Event Timeline
Log In to Comment