Skip to content
Snippets Groups Projects
Commit be7a070a authored by Gloria Corti's avatar Gloria Corti
Browse files

Merge branch 'LHCBGAUSS-1584.EvtSLDiBaryonAmpFix' into 'Sim09' 

Use correct B reference frame in EvtSLDiBaryonAmp: LHCBGAUSS-1584

See merge request lhcb/Gauss!415
parents 3f754c14 1efbfc29
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Gen/EvtGen/EvtGenModels/EvtBToDiBaryonlnupQCD.hh
+ 2
− 2
View file @ be7a070a
......@@ -48,8 +48,8 @@ public:
private:
EvtBToDiBaryonlnupQCDFF* ffModel;
EvtSLDiBaryonAmp* calcAmp;
EvtBToDiBaryonlnupQCDFF* ffModel_;
EvtSLDiBaryonAmp* calcAmp_;
};
......
Gen/EvtGen/EvtGenModels/EvtBToDiBaryonlnupQCDFF.hh
+ 4
− 4
View file @ be7a070a
......@@ -17,7 +17,7 @@
//
// Mark Smith July 17, 2017 Module created
// John B Oct 2018 Added FormFactors class
//
//
//------------------------------------------------------------------------
#ifndef EVTBTODIBARYONLNUPQCDFF_HH
......@@ -46,13 +46,13 @@ public:
EvtBToDiBaryonlnupQCDFF(std::vector<double>& DParameters);
void getDiracFF(EvtParticle* parent, double dibaryonMass,
EvtBToDiBaryonlnupQCDFF::FormFactors& FF);
EvtBToDiBaryonlnupQCDFF::FormFactors& FF) const;
void getRaritaFF(EvtParticle* parent, double dibaryonMass,
EvtBToDiBaryonlnupQCDFF::FormFactors& FF);
EvtBToDiBaryonlnupQCDFF::FormFactors& FF) const;
void getFF(EvtParticle* parent, double dibaryonMass,
EvtBToDiBaryonlnupQCDFF::FormFactors& FF);
EvtBToDiBaryonlnupQCDFF::FormFactors& FF) const;
private:
......
Gen/EvtGen/EvtGenModels/EvtSLDiBaryonAmp.hh
+ 20
− 6
View file @ be7a070a
......@@ -23,25 +23,39 @@
#define EVTSLDIBARYONAMP_HH
#include "EvtGenBase/EvtAmp.hh"
#include "EvtGenBase/EvtDiracSpinor.hh"
#include "EvtGenBase/EvtId.hh"
#include "EvtGenModels/EvtBToDiBaryonlnupQCDFF.hh"
#include "EvtGenBase/EvtVector4R.hh"
#include "EvtGenBase/EvtVector4C.hh"
#include <vector>
class EvtParticle;
class EvtBToDiBaryonlnupQCDFF;
class EvtSLDiBaryonAmp {
public:
EvtSLDiBaryonAmp(EvtBToDiBaryonlnupQCDFF*);
EvtSLDiBaryonAmp(const EvtBToDiBaryonlnupQCDFF&);
virtual ~EvtSLDiBaryonAmp() {;}
void CalcAmp(EvtParticle *parent, EvtAmp& amp) const;
private:
protected:
int checkDibaryonParity(const EvtId& id1, const EvtId& id2,
const int J1, const int J2) const;
EvtBToDiBaryonlnupQCDFF* ffModel;
int getBaryonParity(const EvtId& id) const;
std::vector<EvtVector4C> getHadronicCurrents(const EvtDiracSpinor& u, const EvtDiracSpinor& v,
const EvtVector4R& p, const EvtVector4R& gMtmTerms,
const EvtVector4R& fMtmTerms) const;
private:
EvtBToDiBaryonlnupQCDFF ffModel_;
};
#endif
......
Gen/EvtGen/src/EvtBToDiBaryonlnupQCD.cpp
+ 19
− 22
View file @ be7a070a
//
// Environment:
// This software is part of the EvtGen package developed jointly
......@@ -42,16 +40,16 @@
using std::endl;
EvtBToDiBaryonlnupQCD::EvtBToDiBaryonlnupQCD() :
ffModel(0),
calcAmp(0)
ffModel_(0),
calcAmp_(0)
{
}
EvtBToDiBaryonlnupQCD::~EvtBToDiBaryonlnupQCD() {
delete ffModel;
ffModel = 0;
delete calcAmp;
calcAmp = 0;
delete ffModel_;
ffModel_ = 0;
delete calcAmp_;
calcAmp_ = 0;
}
std::string EvtBToDiBaryonlnupQCD::getName() {
......@@ -70,7 +68,7 @@ void EvtBToDiBaryonlnupQCD::decay(EvtParticle *p) {
p->initializePhaseSpace(getNDaug(), getDaugs(), true);
calcAmp->CalcAmp(p, _amp2);
calcAmp_->CalcAmp(p, _amp2);
}
......@@ -101,7 +99,7 @@ void EvtBToDiBaryonlnupQCD::init() {
if (parentType != EvtSpinType::SCALAR) {
report(ERROR,"EvtGen") << "EvtBToDiBaryonlnupQCD model expected "
report(ERROR,"EvtGen") << "EvtBToDiBaryonlnupQCD model expected "
<< " a SCALAR parent, found:"
<< EvtPDL::name(getParentId()) << endl;
report(ERROR,"EvtGen") << "Will terminate execution!" << endl;
......@@ -133,7 +131,7 @@ void EvtBToDiBaryonlnupQCD::init() {
}
// Form factor model
ffModel = new EvtBToDiBaryonlnupQCDFF(DPars);
ffModel_ = new EvtBToDiBaryonlnupQCDFF(DPars);
// Set amplitude calculation pointer.
// Accomodate for spin 1/2 (DIRAC) or 3/2 (RARITASCHWINGER) baryons
......@@ -144,15 +142,14 @@ void EvtBToDiBaryonlnupQCD::init() {
(baryon1Type == EvtSpinType::RARITASCHWINGER && baryon2Type == EvtSpinType::DIRAC) ||
(baryon1Type == EvtSpinType::DIRAC && baryon2Type == EvtSpinType::DIRAC) ) {
calcAmp = new EvtSLDiBaryonAmp(ffModel);
calcAmp_ = new EvtSLDiBaryonAmp(*ffModel_);
} else {
report(ERROR,"EvtGen") << "Wrong baryon spin type in EvtBToDiBaryonlnupQCD model. "
<< "Expected spin type " << EvtSpinType::DIRAC
<< " or " << EvtSpinType::RARITASCHWINGER
<< ", found spin types " << baryon1Type
<< " and " << baryon2Type << endl;
<< ", found spin types " << baryon1Type << " and " << baryon2Type << endl;
report(ERROR,"EvtGen") << "Will terminate execution!" << endl;
::abort();
}
......@@ -201,7 +198,7 @@ void EvtBToDiBaryonlnupQCD::initProbMax() {
if (Delta.contains(bar1Id) || Delta.contains(bar2Id)) {
// Delta
setProbMax(2.2e6);
setProbMax(1e7);
} else if (LambdaC.contains(bar1Id) || LambdaC.contains(bar2Id)) {
// Lambda_c+
......@@ -217,31 +214,31 @@ void EvtBToDiBaryonlnupQCD::initProbMax() {
} else if (N1440.contains(bar1Id) || N1440.contains(bar2Id)) {
// N(1440)
setProbMax(1.25e6);
setProbMax(8e5);
} else if (N1520.contains(bar1Id) || N1520.contains(bar2Id)) {
// N(1520)
setProbMax(1.25e6);
setProbMax(8e6);
} else if (N1535.contains(bar1Id) || N1535.contains(bar2Id)) {
// N(1535)
setProbMax(1.25e6);
setProbMax(8e5);
} else if (N1650.contains(bar1Id) || N1650.contains(bar2Id)) {
// N(1650)
setProbMax(1.25e6);
setProbMax(8e5);
} else if (N1700.contains(bar1Id) || N1700.contains(bar2Id)) {
// N(1700)
setProbMax(1.25e6);
setProbMax(4e6);
} else if (N1710.contains(bar1Id) || N1710.contains(bar2Id)) {
// N(1710)
setProbMax(1.25e6);
setProbMax(5e5);
} else if (N1720.contains(bar1Id) || N1720.contains(bar2Id)) {
// N(1720)
setProbMax(1.25e6);
setProbMax(4e6);
} // Baryon combinations
......
Gen/EvtGen/src/EvtBToDiBaryonlnupQCDFF.cpp
+ 3
− 3
View file @ be7a070a
......@@ -41,7 +41,7 @@ EvtBToDiBaryonlnupQCDFF::EvtBToDiBaryonlnupQCDFF(std::vector<double>& DParameter
}
void EvtBToDiBaryonlnupQCDFF::getFF(EvtParticle*, double dibaryonMass,
EvtBToDiBaryonlnupQCDFF::FormFactors& FF)
EvtBToDiBaryonlnupQCDFF::FormFactors& FF) const
{
if (nDPars == 6 && dibaryonMass > 0.0) {
......@@ -73,7 +73,7 @@ void EvtBToDiBaryonlnupQCDFF::getFF(EvtParticle*, double dibaryonMass,
}
void EvtBToDiBaryonlnupQCDFF::getDiracFF(EvtParticle* parent, double dibaryonMass,
EvtBToDiBaryonlnupQCDFF::FormFactors &FF)
EvtBToDiBaryonlnupQCDFF::FormFactors &FF) const
{
this->getFF(parent, dibaryonMass, FF);
......@@ -81,7 +81,7 @@ void EvtBToDiBaryonlnupQCDFF::getDiracFF(EvtParticle* parent, double dibaryonMas
}
void EvtBToDiBaryonlnupQCDFF::getRaritaFF(EvtParticle* parent, double dibaryonMass,
EvtBToDiBaryonlnupQCDFF::FormFactors &FF)
EvtBToDiBaryonlnupQCDFF::FormFactors &FF) const
{
this->getFF(parent, dibaryonMass, FF);
......
Gen/EvtGen/src/EvtSLDiBaryonAmp.cpp
+ 221
− 123
View file @ be7a070a
......@@ -25,23 +25,20 @@
#include "EvtGenBase/EvtPatches.hh"
#include "EvtGenBase/EvtDiracSpinor.hh"
#include "EvtGenBase/EvtGammaMatrix.hh"
#include "EvtGenBase/EvtId.hh"
#include "EvtGenBase/EvtIdSet.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtRaritaSchwinger.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtVector4C.hh"
#include "EvtGenBase/EvtVector4R.hh"
#include "EvtGenBase/EvtTensor4C.hh"
#include "EvtGenModels/EvtBToDiBaryonlnupQCDFF.hh"
#include "EvtGenModels/EvtSLDiBaryonAmp.hh"
using std::endl;
EvtSLDiBaryonAmp::EvtSLDiBaryonAmp(EvtBToDiBaryonlnupQCDFF* formFactors) :
ffModel(formFactors)
EvtSLDiBaryonAmp::EvtSLDiBaryonAmp(const EvtBToDiBaryonlnupQCDFF& formFactors) :
ffModel_(formFactors)
{
}
......@@ -54,27 +51,56 @@ void EvtSLDiBaryonAmp::CalcAmp(EvtParticle *parent, EvtAmp& amp) const {
static EvtId MUP = EvtPDL::getId("mu+");
static EvtId TAUP = EvtPDL::getId("tau+");
// The amplitude assumes B- -> p+ p- l- nubar ordering
// i.e. the B- decay is the "particle" mode
// B charge (x3) to check for antiparticle mode and baryon daughter ordering
EvtId BId = parent->getId();
int qB3 = EvtPDL::chg3(BId);
bool particleMode(true);
// Check if we have B+ instead (antiparticle mode)
if (qB3 > 0) {particleMode = false;}
// The baryon, charged lepton and neutrino daughters
// Make sure the first baryon has a charge opposite to the B, since the
// amplitude expressions assume this order
EvtParticle* baryon1 = parent->getDaug(0);
EvtParticle* baryon2 = parent->getDaug(1);
// Check if we need to reverse the baryon ordering
if (EvtPDL::chg3(baryon1->getId()) == qB3) {
baryon1 = parent->getDaug(1);
baryon2 = parent->getDaug(0);
}
EvtParticle* lepton = parent->getDaug(2);
EvtParticle* neutrino = parent->getDaug(3);
// 4-momenta
// 4-momenta in B rest frame
EvtVector4R p0(parent->mass(), 0.0, 0.0, 0.0);
EvtVector4R p1 = baryon1->getP4();
EvtVector4R p2 = baryon2->getP4();
EvtVector4R pSum = p1 + p2;
EvtVector4R p = parent->getP4() - pSum;
EvtVector4R p = p0 - pSum;
EvtVector4R pDiff = p2 - p1;
// Particle id's
// Particle id's: retrieve 1st baryon again in case order has changed
EvtId Id1 = baryon1->getId();
EvtId Id2 = baryon2->getId();
EvtId l_num = lepton->getId();
EvtSpinType::spintype type1 = EvtPDL::getSpinType(Id1);
EvtSpinType::spintype type2 = EvtPDL::getSpinType(Id2);
// Parity of B+- = -1. Check if the parity of the dibaryon state is the same.
// If so, set the sameParity integer to 1. Otherwise set it to -1,
// i.e. the dibaryon system has opposite parity to the B meson
int J1 = EvtSpinType::getSpin2(type1);
int J2 = EvtSpinType::getSpin2(type2);
int sameParity = this->checkDibaryonParity(Id1, Id2, J1, J2);
// Number of chiral components of the baryon spinors
int N1 = EvtSpinType::getSpinStates(type1);
......@@ -91,45 +117,54 @@ void EvtSLDiBaryonAmp::CalcAmp(EvtParticle *parent, EvtAmp& amp) const {
if (l_num == EM || l_num == MUM || l_num == TAUM) {
// B-
l1 = EvtLeptonVACurrent(lepton->spParent(0), neutrino->spParentNeutrino());
l2 = EvtLeptonVACurrent(lepton->spParent(1), neutrino->spParentNeutrino());
} else if (l_num == EP || l_num == MUP || l_num == TAUP) {
// B+
l1 = EvtLeptonVACurrent(neutrino->spParentNeutrino(), lepton->spParent(0));
l2 = EvtLeptonVACurrent(neutrino->spParentNeutrino(), lepton->spParent(1));
l2 = EvtLeptonVACurrent(neutrino->spParentNeutrino(), lepton->spParent(1));
} else {
report(ERROR,"EvtSLDiBaryonAmp") << "Wrong lepton number"<<endl;
}
// Sigma summation terms. Store in a vector to avoid recalculating this for each i,j loop
std::vector<EvtGammaMatrix> sigmaVect;
for (int k = 0; k < 4; k++) {
EvtGammaMatrix sigmaSum;
for (int index = 0; index < 4; index++) {
sigmaSum += EvtGammaMatrix::sigmaLower(k, index)*p.get(index);
}
sigmaVect.push_back(sigmaSum);
}
// Vector and axial-vector terms; these get reset within the loop over k
EvtVector4C term1;
EvtVector4C term2;
// Parity multiplication factors for the antiparticle mode hadronic currents
double sign1 = (particleMode == true) ? 1.0: 1.0*sameParity;
double sign2 = (particleMode == true) ? 1.0: 1.0*sameParity;
double sign3 = (particleMode == true) ? 1.0: -1.0*sameParity;
double sign4 = (particleMode == true) ? 1.0: -1.0*sameParity;
double sign5 = (particleMode == true) ? 1.0: -1.0*sameParity;
double sign6 = (particleMode == true) ? 1.0: 1.0*sameParity;
// Handle case of two Dirac-type daughters, e.g. ppbar, pN(1440)
// Define form factor coeff variables
double f1(0.0), f2(0.0), f3(0.0), f4(0.0), f5(0.0);
double g1(0.0), g2(0.0), g3(0.0), g4(0.0), g5(0.0);
// Handle case of two Dirac-type daughters, e.g. p pbar, p N(1440)
if (type1 == EvtSpinType::DIRAC && type2 == EvtSpinType::DIRAC) {
// Form factor parameters
EvtBToDiBaryonlnupQCDFF::FormFactors FF;
ffModel->getDiracFF(parent, m_dibaryon, FF);
ffModel_.getDiracFF(parent, m_dibaryon, FF);
if (sameParity == 1) {
f1 = FF.F1; f2 = FF.F2; f3 = FF.F3; f4 = FF.F4; f5 = FF.F5;
g1 = FF.G1; g2 = FF.G2; g3 = FF.G3; g4 = FF.G4; g5 = FF.G5;
} else {
// Swap coeffs: f_i <--> g_i
f1 = FF.G1; f2 = FF.G2; f3 = FF.G3; f4 = FF.G4; f5 = FF.G5;
g1 = FF.F1; g2 = FF.F2; g3 = FF.F3; g4 = FF.F4; g5 = FF.F5;
}
EvtVector4R gMtmTerms = g3*p + g4*pSum + g5*pDiff;
EvtVector4R fMtmTerms = f3*p + f4*pSum + f5*pDiff;
// First baryon
for (int i = 0; i < N1; i++) {
......@@ -141,33 +176,27 @@ void EvtSLDiBaryonAmp::CalcAmp(EvtParticle *parent, EvtAmp& amp) const {
// Second baryon
for(int j = 0; j < N2; j++) {
EvtDiracSpinor v = baryon2->spParent(j);
EvtDiracSpinor g5v = EvtGammaMatrix::g5()*v;
// Calculate and set the 4 complex amplitude components
for (int k = 0; k < 4; k++) {
EvtDiracSpinor v = baryon2->spParent(j);
EvtGammaMatrix sigmaSum = sigmaVect[k];
// Need two terms owing to the requirements of using the product operators
EvtDiracSpinor vgTermA = (FF.G1*EvtGammaMatrix::g(k) + I*FF.G2*sigmaSum)*g5v;
EvtDiracSpinor vgTermB = (FF.G3*p.get(k) + FF.G4*pSum.get(k) + FF.G5*pDiff.get(k))*g5v;
// S current = ubar adjoint*vTerms
term1.set(k, EvtLeptonSCurrent(u, vgTermA+vgTermB));
// Hadronic currents
std::vector<EvtVector4C> hadCurrents = this->getHadronicCurrents(u, v, p, gMtmTerms,
fMtmTerms);
EvtDiracSpinor vfTermA = (FF.F1*EvtGammaMatrix::g(k) + I*FF.F2*sigmaSum)*v;
EvtDiracSpinor vfTermB = (FF.F3*p.get(k) + FF.F4*pSum.get(k) + FF.F5*pDiff.get(k))*v;
// First amplitude terms: 3rd current already has the form factor coeffs applied (gMtmTerms)
EvtVector4C amp1 = g1*sign1*hadCurrents[0] + g2*sign2*hadCurrents[1] + sign3*hadCurrents[2];
// S current = ubar adjoint*vTerms
term2.set(k, EvtLeptonSCurrent(u, vfTermA+vfTermB));
// Second amplitude terms: 6th current already has the form factor coeffs applied (fMtmTerms)
EvtVector4C amp2 = f1*sign4*hadCurrents[3] + f2*sign5*hadCurrents[4] + sign6*hadCurrents[5];
EvtVector4C hadAmp;
if (sameParity == 1) {
hadAmp = amp1 - amp2;
} else {
hadAmp = amp2 - amp1;
}
// Set the decay amplitude element
EvtVector4C term = term1 - term2;
amp.vertex(i,j,0,l1*term);
amp.vertex(i,j,1,l2*term);
amp.vertex(i, j, 0, l1*hadAmp);
amp.vertex(i, j, 1, l2*hadAmp);
} // j
......@@ -181,7 +210,19 @@ void EvtSLDiBaryonAmp::CalcAmp(EvtParticle *parent, EvtAmp& amp) const {
// Form factor parameters
EvtBToDiBaryonlnupQCDFF::FormFactors FF;
ffModel->getRaritaFF(parent, m_dibaryon, FF);
ffModel_.getRaritaFF(parent, m_dibaryon, FF);
if (sameParity == 1) {
f1 = FF.F1; f2 = FF.F2; f3 = FF.F3; f4 = FF.F4; f5 = FF.F5;
g1 = FF.G1; g2 = FF.G2; g3 = FF.G3; g4 = FF.G4; g5 = FF.G5;
} else {
// Swap coeffs: f_i <--> g_i
f1 = FF.G1; f2 = FF.G2; f3 = FF.G3; f4 = FF.G4; f5 = FF.G5;
g1 = FF.F1; g2 = FF.F2; g3 = FF.F3; g4 = FF.F4; g5 = FF.F5;
}
EvtVector4R gMtmTerms = g3*p + g4*pSum + g5*pDiff;
EvtVector4R fMtmTerms = f3*p + f4*pSum + f5*pDiff;
if (type1 == EvtSpinType::DIRAC) {
......@@ -195,48 +236,33 @@ void EvtSLDiBaryonAmp::CalcAmp(EvtParticle *parent, EvtAmp& amp) const {
// Second baryon is RS-type
for (int j = 0; j < N2; j++) {
EvtRaritaSchwinger vRS = baryon2->spRSParent(j);
// Store products of g5 with the spinor components as well as
// EvtDiracSpinors to limit constructor calls inside k loop
std::vector<EvtDiracSpinor> g5vVect, vVect;
for (int index = 0; index < 4; index++) {
EvtDiracSpinor v = vRS.getSpinor(index);
EvtDiracSpinor g5v = EvtGammaMatrix::g5()*v;
g5vVect.push_back(g5v);
vVect.push_back(v);
EvtRaritaSchwinger vRS = baryon2->spRSParent(j);
EvtDiracSpinor v;
for (int k = 0; k < 4; k++) {
v.set_spinor(k, vRS.getVector(k)*p0);
}
// Calculate and set the 4 complex amplitude components
for (int k = 0; k < 4; k++) {
EvtGammaMatrix sigmaSum = sigmaVect[k];
EvtDiracSpinor g5v = g5vVect[k];
EvtDiracSpinor v = vVect[k];
// Hadronic currents
std::vector<EvtVector4C> hadCurrents = this->getHadronicCurrents(u, v, p, gMtmTerms,
fMtmTerms);
// Need two terms owing to the requirements of using the product operators
EvtDiracSpinor vgTermA = (FF.G1*EvtGammaMatrix::g(k) + I*FF.G2*sigmaSum)*g5v;
EvtDiracSpinor vgTermB = (FF.G3*p.get(k) + FF.G4*pSum.get(k) + FF.G5*pDiff.get(k))*g5v;
// First amplitude terms: 3rd current already has the form factor coeffs applied (gMtmTerms)
EvtVector4C amp1 = g1*sign1*hadCurrents[0] + g2*sign2*hadCurrents[1] + sign3*hadCurrents[2];
// S current = ubar adjoint*vTerms
term1.set(k, EvtLeptonSCurrent(u, vgTermA+vgTermB));
EvtDiracSpinor vfTermA = (FF.F1*EvtGammaMatrix::g(k) + I*FF.F2*sigmaSum)*v;
EvtDiracSpinor vfTermB = (FF.F3*p.get(k) + FF.F4*pSum.get(k) + FF.F5*pDiff.get(k))*v;
// Second amplitude terms: 6th current already has the form factor coeffs applied (fMtmTerms)
EvtVector4C amp2 = f1*sign4*hadCurrents[3] + f2*sign5*hadCurrents[4] + sign6*hadCurrents[5];
// S current = ubar adjoint*vTerms
term2.set(k, EvtLeptonSCurrent(u, vfTermA+vfTermB));
EvtVector4C hadAmp;
if (sameParity == 1) {
hadAmp = amp1 - amp2;
} else {
hadAmp = amp2 - amp1;
}
// Set the decay amplitude element
EvtVector4C term = term1 - term2;
amp.vertex(i,j,0,l1*term);
amp.vertex(i,j,1,l2*term);
amp.vertex(i, j, 0, l1*hadAmp);
amp.vertex(i, j, 1, l2*hadAmp);
} // j
} // i
......@@ -250,50 +276,37 @@ void EvtSLDiBaryonAmp::CalcAmp(EvtParticle *parent, EvtAmp& amp) const {
// Get the baryon spinor in the B rest frame
EvtRaritaSchwinger uRS = baryon1->spRSParent(i);
// Store EvtDiracSpinors to reduce constructor calls inside k loop
std::vector<EvtDiracSpinor> uVect;
for (int index = 0; index < 4; index++) {
// Just use u and not i*u, since the imaginary constant factor is not needed
// for the probability
EvtDiracSpinor u = uRS.getSpinor(index);
uVect.push_back(u);
EvtDiracSpinor u;
for (int k = 0; k < 4; k++) {
u.set_spinor(k, uRS.getVector(k)*p0);
}
// Second baryon is Dirac
for (int j = 0; j < N2; j++) {
EvtDiracSpinor v = baryon2->spParent(j);
EvtDiracSpinor g5v = EvtGammaMatrix::g5()*v;
// Calculate and set the 4 complex amplitude components
for (int k = 0; k < 4; k++) {
EvtGammaMatrix sigmaSum = sigmaVect[k];
EvtDiracSpinor u = uVect[k];
// Need two terms owing to the requirements of using the product operators
EvtDiracSpinor vgTermA = (FF.G1*EvtGammaMatrix::g(k) + I*FF.G2*sigmaSum)*g5v;
EvtDiracSpinor vgTermB = (FF.G3*p.get(k) + FF.G4*pSum.get(k) + FF.G5*pDiff.get(k))*g5v;
// S current = ubar adjoint*vTerms
term1.set(k, EvtLeptonSCurrent(u, vgTermA+vgTermB));
EvtDiracSpinor vfTermA = (FF.F1*EvtGammaMatrix::g(k) + I*FF.F2*sigmaSum)*v;
EvtDiracSpinor vfTermB = (FF.F3*p.get(k) + FF.F4*pSum.get(k) + FF.F5*pDiff.get(k))*v;
// S current = ubar adjoint*vTerms
term2.set(k, EvtLeptonSCurrent(u, vfTermA+vfTermB));
// Hadronic currents
std::vector<EvtVector4C> hadCurrents = this->getHadronicCurrents(u, v, p, gMtmTerms,
fMtmTerms);
// First amplitude terms: 3rd current already has the form factor coeffs applied (gMtmTerms)
EvtVector4C amp1 = g1*sign1*hadCurrents[0] + g2*sign2*hadCurrents[1] + sign3*hadCurrents[2];
// Second amplitude terms: 6th current already has the form factor coeffs applied (fMtmTerms)
EvtVector4C amp2 = f1*sign4*hadCurrents[3] + f2*sign5*hadCurrents[4] + sign6*hadCurrents[5];
EvtVector4C hadAmp;
if (sameParity == 1) {
hadAmp = amp1 - amp2;
} else {
hadAmp = amp2 - amp1;
}
// Set the decay amplitude element
EvtVector4C term = term1 - term2;
amp.vertex(i,j,0,l1*term);
amp.vertex(i,j,1,l2*term);
amp.vertex(i, j, 0, l1*hadAmp);
amp.vertex(i, j, 1, l2*hadAmp);
} // j
} // i
......@@ -303,3 +316,88 @@ void EvtSLDiBaryonAmp::CalcAmp(EvtParticle *parent, EvtAmp& amp) const {
} // Have Dirac and RS baryons
}
std::vector<EvtVector4C> EvtSLDiBaryonAmp::getHadronicCurrents(const EvtDiracSpinor& u, const EvtDiracSpinor& v,
const EvtVector4R& p, const EvtVector4R& gMtmTerms,
const EvtVector4R& fMtmTerms) const {
// Store the currents used in Eq 6 (in order of appearance)
std::vector<EvtVector4C> currents;
currents.reserve(6);
EvtDiracSpinor g5v = EvtGammaMatrix::g5()*v;
// ubar*gamma*gamma5*v
EvtVector4C current1 = EvtLeptonACurrent(u, v);
currents.push_back(current1);
// ubar*sigma*p*gamma5*v -> [ubar*sigma*(gamma5*v)]*p
EvtTensor4C TC1 = EvtLeptonTCurrent(u, g5v);
// Contract tensor with 4-momentum
EvtVector4C current2 = TC1.cont2(p);
currents.push_back(current2);
// ubar*p*gamma5*v; "p" = p, pSum and pDiff
EvtComplex PC1 = EvtLeptonPCurrent(u, v);
EvtVector4C current3 = PC1*gMtmTerms;
currents.push_back(current3);
// ubar*gamma*v
EvtVector4C current4 = EvtLeptonVCurrent(u, v);
currents.push_back(current4);
// ubar*sigma*p*v -> [ubar*sigma*v]*p
EvtTensor4C TC2 = EvtLeptonTCurrent(u, v);
// Contract tensor with 4-momentum
EvtVector4C current5 = TC2.cont2(p);
currents.push_back(current5);
// ubar*p*v; "p" = p, pSum and pDiff
EvtComplex S1 = EvtLeptonSCurrent(u, v);
EvtVector4C current6 = S1*fMtmTerms;
currents.push_back(current6);
return currents;
}
int EvtSLDiBaryonAmp::checkDibaryonParity(const EvtId& id1, const EvtId& id2,
const int J1, const int J2) const {
// Get intrisic parities of the two baryons, then multiply by (-1)^|J1 - J2|.
// Note here that the J1 and J2 function arguments = 2*spin
int par1 = this->getBaryonParity(id1);
int par2 = this->getBaryonParity(id2);
// mult should be either 0 or 1 for allowed Dirac/RS baryon pairs
int mult = static_cast<int>(pow(-1.0, 0.5*fabs(J1 - J2)));
int dbParity = par1*par2*mult;
// Initialise result to 1, i.e. dibaryon parity = B parity = -1
int result(1);
// Dibaryon parity is opposite to the negative B parity
if (dbParity > 0) {result = -1;}
return result;
}
int EvtSLDiBaryonAmp::getBaryonParity(const EvtId& id) const {
// Initialise parity to +1
int parity(1);
// List of baryons with parity = +1
static EvtIdSet posParity("p+", "Delta+", "Lambda_c+", "anti-Lambda_c(2593)-",
"anti-Lambda_c(2625)-", "N(1440)+", "anti-N(1520)-",
"anti-N(1535)-", "anti-N(1650)-", "anti-N(1700)-",
"N(1710)+", "N(1720)+");
// If the baryon id is not in the list, set the parity to -1
if (!posParity.contains(id)) {parity = -1;}
return parity;
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment