Commit 94896f30 authored by Pawel Sznajder's avatar Pawel Sznajder
Browse files

Merge branch 'new_channels_this_is_it_mesons' of...

Merge branch 'new_channels_this_is_it_mesons' of git@drf-gitlab.cea.fr:partons/core/partons.git into new_channels_this_is_it_mesons
parents dc6c1565 19899274
......@@ -33,7 +33,7 @@ namespace PARTONS {
* 4) \bar{E}_T: Table 1 in arxiv:1901.11380
* with PDFs (forward-limits) parametrizations ABM: arXiv:1202.2281 and DSSV: arXiv:0904.3821
*
* Parameterization of GPDs H and E as in GPDGK19
* Parameterization of GPDs H and E as in GPDGK16
*/
class GPDGK19: public GPDGK16 {
......
......@@ -68,6 +68,15 @@ private:
double CrossSectionLT(); ///< Partial cross-section LT
double CrossSectionTT(); ///< Partial cross-section TT
double poleResidue(); ///< The residue of the pole in Pi+ production
double poleAmplitude0p0p(); ///< Computation of the pion pole amplitude M_{0+0+}.
double poleAmplitude0m0p(); ///< Computation of the pion pole amplitude M_{0-0+}.
double poleAmplitude0mpp(); ///< Computation of the pion pole amplitude M_{0-++}.
double poleAmplitude0ppp(); ///< Computation of the pion pole amplitude M_{0+++}.
double poleAmplitude0pmp(); ///< Computation of the pion pole amplitude M_{0+-+}.
double poleAmplitude0mmp(); ///< Computation of the pion pole amplitude M_{0--+}.
std::complex<double> Amplitude0p0p(); ///< Computation of the amplitude M_{0+0+}.
std::complex<double> Amplitude0m0p(); ///< Computation of the amplitude M_{0-0+}.
std::complex<double> Amplitude0mpp(); ///< Computation of the amplitude M_{0-++}.
......
......@@ -211,10 +211,11 @@ std::complex<double> DVMPCFFGK06::computeCFF() {
//switch over mesons
switch (m_currentGPDComputeType) {
//Ht
//Ht
case GPDType::Ht: {
return convolutionTwist2(m_currentGPDComputeType);
}
break;
......@@ -223,6 +224,7 @@ std::complex<double> DVMPCFFGK06::computeCFF() {
case GPDType::Et: {
return convolutionTwist2(m_currentGPDComputeType);
}
break;
......@@ -233,6 +235,7 @@ std::complex<double> DVMPCFFGK06::computeCFF() {
return convolutionTwist3A(m_currentGPDComputeType)
+ convolutionTwist3B(m_currentGPDComputeType)
+ convolutionTwist3C(m_currentGPDComputeType);
}
break;
......@@ -243,6 +246,7 @@ std::complex<double> DVMPCFFGK06::computeCFF() {
return convolutionTwist3A(m_currentGPDComputeType)
+ convolutionTwist3B(m_currentGPDComputeType)
+ convolutionTwist3C(m_currentGPDComputeType);
}
break;
......@@ -359,7 +363,26 @@ std::pair<double, double> DVMPCFFGK06::mesonWFParameters(size_t twist) const {
switch (m_mesonType) {
//pi0, pi+
case MesonType::PI0:
case MesonType::PI0: {
//twist 2
if (twist == 2) {
return std::pair<double, double>(0.132,
1. / (8. * pow(M_PI, 2.0) * pow(0.132, 2.)));
}
//twist 3
else if (twist == 3) {
return std::pair<double, double>(0.132, 1.8);
}
//???
else {
throw ElemUtils::CustomException(getClassName(), __func__,
ElemUtils::Formatter() << "Twist " << twist
<< " not implemented here for meson "
<< MesonType(m_mesonType).toString());
}
}
case MesonType::PIPLUS: {
//twist 2
......@@ -412,10 +435,10 @@ double DVMPCFFGK06::mesonWF(double tau, double b, size_t twist) const {
return 4. * M_PI * par.first / sqrt(2. * m_Nc) * m_muPi
* pow(par.second, 2.)
* exp(
-1.0 * pow(b, 2.) / (8. * pow(par.second, 2.0))
-1.0 * pow(b, 2.) / (8. * pow(par.second, 2.0)))
* gsl_sf_bessel_In(0,
pow(b, 2.)
/ (8. * pow(par.second, 2.0))));
/ (8. * pow(par.second, 2.0)));
}
//???
else {
......@@ -446,20 +469,16 @@ double DVMPCFFGK06::getMesonGPDCombination(
break;
//pi+
//pi+
case MesonType::PIPLUS: {
return 1. / sqrt(2.)
* (Constant::U_ELEC_CHARGE
* partonDistribution.getQuarkDistribution(
QuarkFlavor::UP).getQuarkDistribution()
- Constant::D_ELEC_CHARGE
* partonDistribution.getQuarkDistribution(
QuarkFlavor::DOWN).getQuarkDistribution());
return partonDistribution.getQuarkDistribution(QuarkFlavor::UP).getQuarkDistribution()
- partonDistribution.getQuarkDistribution(QuarkFlavor::DOWN).getQuarkDistribution();
}
break;
//???
//???
default: {
throw ElemUtils::CustomException(getClassName(), __func__,
ElemUtils::Formatter() << "No implementation for meson: "
......@@ -475,41 +494,40 @@ double DVMPCFFGK06::getMesonGPDCombination(
std::complex<double> DVMPCFFGK06::subProcess(double x, double tau, double b,
size_t twist) const {
std::complex<double> Ts, Tu;
if (x >= m_xi) {
Ts = -1. * std::complex<double>(0., 1.) / 4.
* (gsl_sf_bessel_J0(
sqrt((1. - tau) * (x - m_xi) / (2. * m_xi)) * b
* sqrt(m_Q2))
+ std::complex<double>(0., 1.)
* gsl_sf_bessel_Y0(
sqrt(
(1. - tau) * (x - m_xi)
/ (2. * m_xi))
* b * sqrt(m_Q2)));
} else {
Ts = -1. / (2. * M_PI)
* gsl_sf_bessel_K0(
sqrt((1. - tau) * (m_xi - x) / (2. * m_xi)) * b
* sqrt(m_Q2));
}
Tu = -1. / (2. * M_PI)
* gsl_sf_bessel_K0(
sqrt(tau * (x + m_xi) / (2. * m_xi)) * b
* sqrt(m_Q2));
//switch over mesons
switch (m_mesonType) {
//pi0 pi+
case MesonType::PI0:
case MesonType::PIPLUS: {
std::complex<double> Ts, Tu;
case MesonType::PI0: {
//twist 2
if (twist == 2) {
if (x >= m_xi) {
Ts = -1. * std::complex<double>(0., 1.) / 4.
* (gsl_sf_bessel_J0(
sqrt((1. - tau) * (x - m_xi) / (2. * m_xi)) * b
* sqrt(m_Q2))
+ std::complex<double>(0., 1.)
* gsl_sf_bessel_Y0(
sqrt(
(1. - tau) * (x - m_xi)
/ (2. * m_xi))
* b * sqrt(m_Q2)));
} else {
Ts = -1. / (2. * M_PI)
* gsl_sf_bessel_K0(
sqrt((1. - tau) * (m_xi - x) / (2. * m_xi)) * b
* sqrt(m_Q2));
}
Tu = -1. / (2. * M_PI)
* gsl_sf_bessel_K0(
sqrt(tau * (x + m_xi) / (2. * m_xi)) * b
* sqrt(m_Q2));
return m_Cf * sqrt(2. / m_Nc) * m_Q2 / m_xi * 2. * M_PI * b
* mesonWF(tau, b, twist)
* m_pRunningAlphaStrongModule->compute(
......@@ -519,34 +537,42 @@ std::complex<double> DVMPCFFGK06::subProcess(double x, double tau, double b,
//twist 3
else if (twist == 3) {
if (x >= m_xi) {
Ts = -1. * std::complex<double>(0., 1.) / 4.
* (gsl_sf_bessel_J0(
sqrt((1. - tau) * (x - m_xi) / (2. * m_xi)) * b
* sqrt(m_Q2))
+ std::complex<double>(0., 1.)
* gsl_sf_bessel_Y0(
sqrt(
(1. - tau) * (x - m_xi)
/ (2. * m_xi))
* b * sqrt(m_Q2)));
} else {
Ts = -1. / (2. * M_PI)
* gsl_sf_bessel_K0(
sqrt((1. - tau) * (m_xi - x) / (2. * m_xi)) * b
* sqrt(m_Q2));
}
return 4.0 * m_Cf / sqrt(2. * m_Nc) * m_Q2 / m_xi * 2. * M_PI * b
* mesonWF(tau, b, twist)
* m_pRunningAlphaStrongModule->compute(
pow(computeMuR(tau, b), 2))
* expSudakovFactor(tau, b) * ((1. - tau) * Ts + tau * Tu);
}
//???
else {
throw ElemUtils::CustomException(getClassName(), __func__,
ElemUtils::Formatter() << "Twist " << twist
<< " not implemented here for meson "
<< MesonType(m_mesonType).toString());
}
}
case MesonType::PIPLUS: {
//twist 2
if (twist == 2) {
Tu = -1. / (2. * M_PI)
* gsl_sf_bessel_K0(
sqrt(tau * (x + m_xi) / (2. * m_xi)) * b
* sqrt(m_Q2));
return m_Cf * sqrt(2. / m_Nc) * m_Q2 / m_xi * 2. * M_PI * b
* mesonWF(tau, b, twist)
* m_pRunningAlphaStrongModule->compute(
pow(computeMuR(tau, b), 2))
* expSudakovFactor(tau, b) * (Constant::U_ELEC_CHARGE * Ts - Constant::D_ELEC_CHARGE * Tu);
}
//twist 3
else if (twist == 3) {
return 4.0 * m_Cf / sqrt(2. * m_Nc) * m_Q2 / m_xi * 2. * M_PI * b
* mesonWF(tau, b, twist)
* m_pRunningAlphaStrongModule->compute(
pow(computeMuR(tau, b), 2))
* expSudakovFactor(tau, b) * ((1. - tau) * Ts + tau * Tu);
* expSudakovFactor(tau, b) * ((1. - tau) * Constant::U_ELEC_CHARGE * Ts + tau * Constant::D_ELEC_CHARGE * Tu);
}
//???
......@@ -818,17 +844,58 @@ double DVMPCFFGK06::convolutionTwist3BFunction(double x, void * params) const {
// - getMesonGPDCombination(m_gpdResultXiXi))
// / (x - m_xi);
double convolution = 1. / (x + m_xi)
* (getMesonGPDCombination(gpdResult)
- getMesonGPDCombination(m_gpdResultmXiXi))
- 1.
//switch over mesons
switch (m_mesonType) {
//pi0
case MesonType::PI0: {
double convolution = 1. / (x + m_xi)
* (getMesonGPDCombination(gpdResult)
- getMesonGPDCombination(m_gpdResultXiXi))
/ (x - m_xi);
- getMesonGPDCombination(m_gpdResultmXiXi))
- 1.
* (getMesonGPDCombination(gpdResult)
- getMesonGPDCombination(m_gpdResultXiXi))
/ (x - m_xi);
return 16. * M_PI * m_Cf / m_Nc
* m_pRunningAlphaStrongModule->compute(m_Q2 / 2.) * par.first
* m_muPi * pow(par.second, 2.) * convolution;
}
break;
//pi+
case MesonType::PIPLUS: {
double convolution = Constant::D_ELEC_CHARGE / (x + m_xi)
* (getMesonGPDCombination(gpdResult)
- getMesonGPDCombination(m_gpdResultmXiXi))
- Constant::U_ELEC_CHARGE
* (getMesonGPDCombination(gpdResult)
- getMesonGPDCombination(m_gpdResultXiXi))
/ (x - m_xi);
return 16. * M_PI * m_Cf / m_Nc
* m_pRunningAlphaStrongModule->compute(m_Q2 / 2.) * par.first
* m_muPi * pow(par.second, 2.) * convolution;
}
break;
//???
default: {
throw ElemUtils::CustomException(getClassName(), __func__,
ElemUtils::Formatter() << "No implementation for meson: "
<< MesonType(m_mesonType).toString());
}
break;
}
return 16. * M_PI * m_Cf / m_Nc
* m_pRunningAlphaStrongModule->compute(m_Q2 / 2.) * par.first
* m_muPi * pow(par.second, 2.) * convolution;
}
std::complex<double> DVMPCFFGK06::convolutionTwist3B(
......@@ -894,15 +961,50 @@ std::complex<double> DVMPCFFGK06::convolutionTwist3C(
// * (std::complex<double>(0., 1.) * M_PI
// - log((1. - m_xi) / (2. * m_xi))));
return 16. * M_PI * m_Cf / m_Nc
* m_pRunningAlphaStrongModule->compute(m_Q2 / 2.) * par.first
* m_muPi * pow(par.second, 2.)
* (getMesonGPDCombination(m_gpdResultXiXi)
* (std::complex<double>(0., 1.) * M_PI
+ log((1. + m_xi) / (1. - m_xi)))
+ getMesonGPDCombination(m_gpdResultmXiXi)
//switch over mesons
switch (m_mesonType) {
//pi0
case MesonType::PI0: {
return 16. * M_PI * m_Cf / m_Nc
* m_pRunningAlphaStrongModule->compute(m_Q2 / 2.) * par.first
* m_muPi * pow(par.second, 2.)
* (getMesonGPDCombination(m_gpdResultXiXi)
* (std::complex<double>(0., 1.) * M_PI
+ log((1. + m_xi) / (1. - m_xi)))
+ getMesonGPDCombination(m_gpdResultmXiXi)
* (std::complex<double>(0., 1.) * M_PI
+ log((1. + m_xi) / (1. - m_xi))));
}
break;
//pi+
case MesonType::PIPLUS: {
return 16. * M_PI * m_Cf / m_Nc
* m_pRunningAlphaStrongModule->compute(m_Q2 / 2.) * par.first
* m_muPi * pow(par.second, 2.)
* (Constant::U_ELEC_CHARGE * getMesonGPDCombination(m_gpdResultXiXi)
* (std::complex<double>(0., 1.) * M_PI
+ log((1. + m_xi) / (1. - m_xi))));
+ log((1. + m_xi) / (1. - m_xi)))
+ Constant::D_ELEC_CHARGE * getMesonGPDCombination(m_gpdResultmXiXi)
* (std::complex<double>(0., 1.) * M_PI
+ log((1. + m_xi) / (1. - m_xi))));
}
break;
//???
default: {
throw ElemUtils::CustomException(getClassName(), __func__,
ElemUtils::Formatter() << "No implementation for meson: "
<< MesonType(m_mesonType).toString());
}
break;
}
}
......
......@@ -189,58 +189,330 @@ double DVMPProcessGK06::lambdaFunction(double a, double b, double c) const {
- 2.0 * (a * b + a * c + b * c);
}
double DVMPProcessGK06::poleResidue() {
//the residue of the pole. See Eq. (8) in https://arxiv.org/pdf/0906.0460.pdf
//The fit parameter to be used in pion form factor
double FitParameter = 0.5;
//A parameter to be used in pion-nucleon vertex
double LambdaN = 0.44;
//Pion-Nucleon coupling constant
double PionNucleonCoupling = 13.1;
//Pion form factor
double FF = 1.0 / (1.0 + m_Q2 / FitParameter);
//Pion-Nucleon vertex
double PionNucleonVertex = (pow(LambdaN, 2.0) - pow(Constant::MESON_PIPLUS_MASS, 2.0)) / (pow(LambdaN, 2.0) - m_t);
//The residue of the pole
double residue = sqrt(2.0) * PionNucleonCoupling * FF * PionNucleonVertex;
return residue;
}
double DVMPProcessGK06::poleAmplitude0p0p() {
double poleAmplitude = - Constant::POSITRON_CHARGE * 2.0 * Constant::PROTON_MASS * m_xi * sqrt(m_Q2) / sqrt(1.0 - pow(m_xi, 2.0)) *
poleResidue() / (m_t - pow(Constant::MESON_PIPLUS_MASS, 2.0));
return poleAmplitude;
}
double DVMPProcessGK06::poleAmplitude0m0p() {
double poleAmplitude = Constant::POSITRON_CHARGE * sqrt(m_Q2) * sqrt(-(m_t - m_tmin)) * poleResidue() / (m_t - pow(Constant::MESON_PIPLUS_MASS, 2.0));
return poleAmplitude;
}
double DVMPProcessGK06::poleAmplitude0ppp() {
double poleAmplitude = Constant::POSITRON_CHARGE * 2.0 * sqrt(2.0) * m_xi * Constant::PROTON_MASS * sqrt(-(m_t - m_tmin)) * poleResidue() /
(m_t - pow(Constant::MESON_PIPLUS_MASS, 2.0));
return poleAmplitude;
}
double DVMPProcessGK06::poleAmplitude0pmp() {
double poleAmplitude = - Constant::POSITRON_CHARGE * 2.0 * sqrt(2.0) * m_xi * Constant::PROTON_MASS * sqrt(-(m_t - m_tmin)) * poleResidue() /
(m_t - pow(Constant::MESON_PIPLUS_MASS, 2.0));
return poleAmplitude;
}
double DVMPProcessGK06::poleAmplitude0mpp() {
double poleAmplitude = Constant::POSITRON_CHARGE * sqrt(2.0) * (m_t - m_tmin) * sqrt(1.0 - pow(m_xi, 2.0)) * poleResidue() /
(m_t - pow(Constant::MESON_PIPLUS_MASS, 2.0));
return poleAmplitude;
}
double DVMPProcessGK06::poleAmplitude0mmp() {
double poleAmplitude = - Constant::POSITRON_CHARGE * sqrt(2.0) * (m_t - m_tmin) * sqrt(1.0 - pow(m_xi, 2.0)) * poleResidue() /
(m_t - pow(Constant::MESON_PIPLUS_MASS, 2.0));
return poleAmplitude;
}
std::complex<double> DVMPProcessGK06::Amplitude0p0p() {
//the handbag amplitude
std::complex<double> amplitude0p0p = sqrt(1. - pow(m_xi, 2.))
* Constant::POSITRON_CHARGE / sqrt(m_Q2)
* (getConvolCoeffFunctionValue(GPDType::Ht)
- pow(m_xi, 2.) / (1. - pow(m_xi, 2.))
* getConvolCoeffFunctionValue(GPDType::Et));
return amplitude0p0p;
//switch over mesons
switch (m_mesonType) {
//pi0
case MesonType::PI0: {
return amplitude0p0p;
}
break;
//pi+
case MesonType::PIPLUS: {
return amplitude0p0p + poleAmplitude0p0p();
}
break;
//???
default: {
throw ElemUtils::CustomException(getClassName(), __func__,
ElemUtils::Formatter() << "No implementation for meson: "
<< MesonType(m_mesonType).toString());
}
break;
}
return std::complex<double>(0., 0.);
}
std::complex<double> DVMPProcessGK06::Amplitude0m0p() {
//the handbag amplitude
std::complex<double> amplitude0m0p = Constant::POSITRON_CHARGE / sqrt(m_Q2)
* sqrt(-(m_t - m_tmin)) * m_xi / (2. * Constant::PROTON_MASS)
* getConvolCoeffFunctionValue(GPDType::Et);
return amplitude0m0p;
}
//switch over mesons
switch (m_mesonType) {
std::complex<double> DVMPProcessGK06::Amplitude0mpp() {
//pi0
case MesonType::PI0: {
std::complex<double> amplitude0mpp = Constant::POSITRON_CHARGE
* sqrt(1. - pow(m_xi, 2.))
* getConvolCoeffFunctionValue(GPDType::HTrans);
return amplitude0m0p;
}
break;
//pi+
case MesonType::PIPLUS: {
return amplitude0m0p + poleAmplitude0m0p();
}
break;
//???
default: {
throw ElemUtils::CustomException(getClassName(), __func__,
ElemUtils::Formatter() << "No implementation for meson: "
<< MesonType(m_mesonType).toString());
}
break;
}
return amplitude0mpp;
return std::complex<double>(0., 0.);
}
std::complex<double> DVMPProcessGK06::Amplitude0ppp() {
//the handbag amplitude
std::complex<double> amplitude0ppp = -1.0 * Constant::POSITRON_CHARGE
* sqrt(-(m_t - m_tmin)) / (4. * Constant::PROTON_MASS)
* getConvolCoeffFunctionValue(GPDType::EbarTrans);
return amplitude0ppp;
//switch over mesons
switch (m_mesonType) {
//pi0
case MesonType::PI0: {
return amplitude0ppp;
}
break;
//pi+
case MesonType::PIPLUS: {
return amplitude0ppp + poleAmplitude0ppp();