Commit d659378e authored by Pawel Sznajder's avatar Pawel Sznajder
Browse files

add masses, tmin/max evaulation, revert evaluation from gpds+

parent bdcc7fe4
......@@ -79,6 +79,20 @@ const double T2_ELEC_CHARGE = 4. / 9.; ///< Square of elect
//TODO more explicit name
const double POSITRON_CHARGE = 0.30282211985966434; ///< Positron charge \f$e\f$ in natural units (\f$e^2 = 4\pi\alpha\f$).
// Mesons
const double MESON_RHOMINUS_MASS = 775.26e-3; ///< \f$\rho^-\f$ mass in \f$\mathrm{GeV}\f$ Phys. 2020, 083C01 (2020).
const double MESON_RHO0_MASS = 775.26e-3; ///< \f$\rho^0\f$ mass in \f$\mathrm{GeV}\f$ Phys. 2020, 083C01 (2020).
const double MESON_RHOPLUS_MASS = 775.26e-3; ///< \f$\rho^+\f$ mass in \f$\mathrm{GeV}\f$ Phys. 2020, 083C01 (2020).
const double MESON_OMEGA_MASS = 782.65e-3; ///< \f$\omega\f$ mass in \f$\mathrm{GeV}\f$ Phys. 2020, 083C01 (2020).
const double MESON_PHI_MASS = 1019.461e-3; ///< \f$\phi\f$ mass in \f$\mathrm{GeV}\f$ Phys. 2020, 083C01 (2020).
const double MESON_JPSI_MASS = 3096.900e-3; ///< \f$J/\Psi\f$ mass in \f$\mathrm{GeV}\f$ Phys. 2020, 083C01 (2020).
const double MESON_UPSILON_MASS = 9460.30e-3; ///< \f$\Upsilon\f$ mass in \f$\mathrm{GeV}\f$ Phys. 2020, 083C01 (2020).
const double MESON_PIMINUS_MASS = 139.57039e-3; ///< \f$\pi^-\f$ mass in \f$\mathrm{GeV}\f$ Phys. 2020, 083C01 (2020).
const double MESON_PI0_MASS = 134.9768e-3; ///< \f$\pi^0\f$ mass in \f$\mathrm{GeV}\f$ Phys. 2020, 083C01 (2020).
const double MESON_PIPLUS_MASS = 139.57039e-3; ///< \f$\pi^+\f$ mass in \f$\mathrm{GeV}\f$ Phys. 2020, 083C01 (2020).
} /* namespace Constant */
} /* namespace PARTONS */
......
......@@ -76,6 +76,11 @@ public:
*/
List<GPDType> getPossibleGPDTypes() const;
/**
* Get mass associated to set meson type.
*/
double getMass() const;
/**
* Get string representation of type being assigned to a declared object of this class.
* @return String representation of assigned type.
......
......@@ -166,7 +166,6 @@ protected:
double m_tmax; ///< Maximal value of t.
double m_xBmin; ///< Minimal value of xB.
double m_y; ///< Lepton energy fraction.
double m_epsilon; ///< @f$ \epsilon = \frac{2 x_B M}{Q} @f$.
DVMPScalesModule* m_pScaleModule; ///< Pointer to the underlying scale module.
DVMPXiConverterModule* m_pXiConverterModule; ///< Pointer to the underlying xi converter module.
......
#include "../../../include/partons/beans/MesonType.h"
#include <ElementaryUtils/string_utils/Formatter.h>
#include <ElementaryUtils/string_utils/StringUtils.h>
#include <ElementaryUtils/thread/Packet.h>
#include "../../../include/partons/FundamentalPhysicalConstants.h"
namespace PARTONS {
MesonType::MesonType() :
......@@ -45,11 +48,60 @@ List<GPDType> MesonType::getPossibleGPDTypes() const {
list.add(GPDType(GPDType::E));
list.add(GPDType(GPDType::HTrans));
list.add(GPDType(GPDType::ETrans));
} else {
throw ElemUtils::CustomException("MesonType", __func__,
ElemUtils::Formatter() << "Not able to return value for type "
<< toString());
}
return list;
}
double MesonType::getMass() const {
switch (m_type) {
case RHOMINUS:
return Constant::MESON_RHOMINUS_MASS;
break;
case RHO0:
return Constant::MESON_RHO0_MASS;
break;
case RHOPLUS:
return Constant::MESON_RHOPLUS_MASS;
break;
case OMEGA:
return Constant::MESON_OMEGA_MASS;
break;
case PHI:
return Constant::MESON_PHI_MASS;
break;
case JPSI:
return Constant::MESON_JPSI_MASS;
break;
case UPSILON:
return Constant::MESON_UPSILON_MASS;
break;
case PIMINUS:
return Constant::MESON_PIMINUS_MASS;
break;
case PI0:
return Constant::MESON_PI0_MASS;
break;
case PIPLUS:
return Constant::MESON_PIPLUS_MASS;
break;
default:
throw ElemUtils::CustomException("MesonType", __func__,
ElemUtils::Formatter() << "Not able to return value for type "
<< toString());
return 0.;
}
}
std::string MesonType::toString() const {
switch (m_type) {
......
......@@ -449,10 +449,10 @@ double DVMPCFFGK06::getMesonGPDCombination(
return 1. / sqrt(2.)
* (Constant::U_ELEC_CHARGE
* partonDistribution.getQuarkDistribution(
QuarkFlavor::UP).getQuarkDistributionPlus()
QuarkFlavor::UP).getQuarkDistribution()
- Constant::D_ELEC_CHARGE
* partonDistribution.getQuarkDistribution(
QuarkFlavor::DOWN).getQuarkDistributionPlus());
QuarkFlavor::DOWN).getQuarkDistribution());
}
break;
......@@ -462,10 +462,10 @@ double DVMPCFFGK06::getMesonGPDCombination(
return 1. / sqrt(2.)
* (Constant::U_ELEC_CHARGE
* partonDistribution.getQuarkDistribution(
QuarkFlavor::UP).getQuarkDistributionPlus()
QuarkFlavor::UP).getQuarkDistribution()
- Constant::D_ELEC_CHARGE
* partonDistribution.getQuarkDistribution(
QuarkFlavor::DOWN).getQuarkDistributionPlus());
QuarkFlavor::DOWN).getQuarkDistribution());
}
break;
......@@ -611,9 +611,9 @@ double DVMPCFFGK06::convolutionFunction(double *xtaub, size_t dim,
std::complex<double> DVMPCFFGK06::convolution(GPDType::Type gpdType,
size_t twist) const {
double resultRe[2], errorRe[2], resultIm[2], errorIm[2]; //results and errors
double resultRe[3], errorRe[3], resultIm[3], errorIm[3]; //results and errors
double rangeMin[3] = { 0., 0.0, 0.0 }; // lower bounds of the 3D integral: with respect to 1) x, 2) tau, and 3) b
double rangeMin[3] = { -1.0, 0.0, 0.0 }; // lower bounds of the 3D integral: with respect to 1) x, 2) tau, and 3) b
double rangeMax[3] = { 1.0, 1.0, 1.0 / m_cLambdaQCD }; // upper bounds of the 3D integral: with respect to 1) x, 2) tau, and 3) b
gsl_rng_env_setup(); //random generator
......@@ -634,11 +634,22 @@ std::complex<double> DVMPCFFGK06::convolution(GPDType::Type gpdType,
//Re
dvmpCFFGK06IntegrationParameters.m_isReal = true;
for (size_t i = 0; i < 2; i++) {
for (size_t i = 0; i < 3; i++) {
//range
rangeMin[0] = (i == 0) ? (0.) : (m_xi);
rangeMax[0] = (i == 0) ? (m_xi) : (1.);
if (i == 0) {
rangeMin[0] = -1.;
rangeMax[0] = -m_xi;
} else if (i == 1) {
rangeMin[0] = -m_xi;
rangeMax[0] = m_xi;
} else if (i == 2) {
rangeMin[0] = m_xi;
rangeMax[0] = 1.;
} else {
throw ElemUtils::CustomException(getClassName(), __func__,
"Wrong index");
}
//state
gslState = gsl_monte_vegas_alloc(3);
......@@ -656,8 +667,29 @@ std::complex<double> DVMPCFFGK06::convolution(GPDType::Type gpdType,
//integrate
do {
if (resultRe[i] == 0. && errorRe[i] == 0.)
if (errorRe[i] == 0.) {
if (std::isnan(resultRe[i])) {
warn(__func__,
ElemUtils::Formatter() << "Range: " << i
<< " for GPD "
<< GPDType(m_currentGPDComputeType).toString()
<< " gives null error, result is NaN and changed to zero");
resultRe[i] = 0.;
} else {
warn(__func__,
ElemUtils::Formatter() << "Range: " << i
<< " for GPD "
<< GPDType(m_currentGPDComputeType).toString()
<< " gives null error, result is: "
<< resultRe[i]);
}
break;
}
gsl_monte_vegas_integrate(&gslFunction, rangeMin, rangeMax, 3,
m_MCCalls, gslRnd, gslState, &resultRe[i], &errorRe[i]);
......@@ -682,11 +714,22 @@ std::complex<double> DVMPCFFGK06::convolution(GPDType::Type gpdType,
//Im parameters
dvmpCFFGK06IntegrationParameters.m_isReal = false;
for (size_t i = 0; i < 2; i++) {
for (size_t i = 0; i < 3; i++) {
//range
rangeMin[0] = (i == 0) ? (0.) : (m_xi);
rangeMax[0] = (i == 0) ? (m_xi) : (1.);
if (i == 0) {
rangeMin[0] = -1.;
rangeMax[0] = -m_xi;
} else if (i == 1) {
rangeMin[0] = -m_xi;
rangeMax[0] = m_xi;
} else if (i == 2) {
rangeMin[0] = m_xi;
rangeMax[0] = 1.;
} else {
throw ElemUtils::CustomException(getClassName(), __func__,
"Wrong index");
}
//state
gslState = gsl_monte_vegas_alloc(3);
......@@ -704,9 +747,29 @@ std::complex<double> DVMPCFFGK06::convolution(GPDType::Type gpdType,
//Im integrate
do {
if (resultIm[i] == 0. && errorIm[i] == 0.)
break;
if (errorIm[i] == 0.) {
if (std::isnan(resultIm[i])) {
warn(__func__,
ElemUtils::Formatter() << "Range: " << i
<< " for GPD "
<< GPDType(m_currentGPDComputeType).toString()
<< " gives null error, result is NaN and changed to zero");
resultIm[i] = 0.;
} else {
warn(__func__,
ElemUtils::Formatter() << "Range: " << i
<< " for GPD "
<< GPDType(m_currentGPDComputeType).toString()
<< " gives null error, result is: "
<< resultIm[i]);
}
break;
}
gsl_monte_vegas_integrate(&gslFunction, rangeMin, rangeMax, 3,
m_MCCalls, gslRnd, gslState, &resultIm[i], &errorIm[i]);
......@@ -731,8 +794,9 @@ std::complex<double> DVMPCFFGK06::convolution(GPDType::Type gpdType,
gsl_rng_free(gslRnd);
// Return
return (resultRe[0] + resultRe[1])
+ std::complex<double>(0., 1.) * (resultIm[0] + resultIm[1]);
return (resultRe[0] + resultRe[1] + resultRe[2])
+ std::complex<double>(0., 1.)
* (resultIm[0] + resultIm[1] + resultIm[2]);
}
std::complex<double> DVMPCFFGK06::convolutionTwist2(
......@@ -775,7 +839,7 @@ std::complex<double> DVMPCFFGK06::convolutionTwist3B(
gsl_integration_workspace* w; // workspace
double result[2], error[2]; //result
double result[3], error[3]; //result
DVMPCFFGK06IntegrationParameters dvmpCFFGK06IntegrationParameters; // parameters
dvmpCFFGK06IntegrationParameters.m_pDVMPCFFGK06 = this;
......@@ -785,18 +849,39 @@ std::complex<double> DVMPCFFGK06::convolutionTwist3B(
gslFunction.function = &gslWrapper1;
gslFunction.params = &dvmpCFFGK06IntegrationParameters;
for (size_t i = 0; i < 2; i++) {
double rangeMin, rangeMax; // range
for (size_t i = 0; i < 3; i++) {
//range
if (i == 0) {
rangeMin = -1.;
rangeMax = -m_xi;
} else if (i == 1) {
rangeMin = -m_xi;
rangeMax = m_xi;
} else if (i == 2) {
rangeMin = m_xi;
rangeMax = 1.;
} else {
throw ElemUtils::CustomException(getClassName(), __func__,
"Wrong index");
}
w = gsl_integration_workspace_alloc(10000); //workspace
gsl_integration_qag(&gslFunction, (i == 0) ? (0.) : (m_xi),
(i == 0) ? (m_xi) : (1.), 0.0, GSL_INTEG_GAUSS61, 1e-5, 10000,
w, &result[i], &error[i]); //evaluate
gsl_integration_qag(&gslFunction, rangeMin, rangeMax, 0.0, 1e-5, 10000,
GSL_INTEG_GAUSS61, w, &result[i], &error[i]); //evaluate
gsl_integration_workspace_free(w); //free
info(__func__,
ElemUtils::Formatter() << "i: " << i << " CFF "
<< GPDType(m_currentGPDComputeType).toString()
<< " result: " << result[i] << " error: " << error[i]);
}
return std::complex<double>(result[0] + result[1], 0.);
return std::complex<double>(result[0] + result[1] + result[2], 0.);
}
std::complex<double> DVMPCFFGK06::convolutionTwist3C(
......
......@@ -48,6 +48,10 @@ void DVMPProcessGK06::initModule() {
getXiConverterModule()->compute(
DVMPObservableKinematic(m_xB, m_t, m_Q2, m_E, m_phi,
m_mesonType)).getValue();
//reevaluate t_min
m_tmin = -4. * pow(Constant::PROTON_MASS, 2.) * pow(m_xi, 2.)
/ (1. - pow(m_xi, 2.));
}
void DVMPProcessGK06::isModuleWellConfigured() {
......
......@@ -31,8 +31,8 @@ DVMPProcessModule::DVMPProcessModule(const std::string &className) :
0.), m_E(0.), m_phi(0.), m_mesonType(MesonType::UNDEFINED), m_beamHelicity(
0.), m_beamCharge(0.), m_mesonPolarization(
MesonPolarization::UNDEFINED), m_tmin(0.), m_tmax(0.), m_xBmin(
0), m_y(0.), m_epsilon(0.), m_pScaleModule(0), m_pXiConverterModule(
0), m_pConvolCoeffFunctionModule(0) {
0), m_y(0.), m_pScaleModule(0), m_pXiConverterModule(0), m_pConvolCoeffFunctionModule(
0) {
}
DVMPProcessModule::~DVMPProcessModule() {
......@@ -60,9 +60,8 @@ DVMPProcessModule::DVMPProcessModule(const DVMPProcessModule& other) :
other.m_beamCharge), m_targetPolarization(
other.m_targetPolarization), m_mesonPolarization(
other.m_mesonPolarization), m_tmin(other.m_tmin), m_tmax(
other.m_tmax), m_xBmin(other.m_xBmin), m_y(other.m_y), m_epsilon(
other.m_epsilon), m_pScaleModule(0), m_pXiConverterModule(0), m_pConvolCoeffFunctionModule(
0) {
other.m_tmax), m_xBmin(other.m_xBmin), m_y(other.m_y), m_pScaleModule(
0), m_pXiConverterModule(0), m_pConvolCoeffFunctionModule(0) {
m_lastCCFKinematics = other.m_lastCCFKinematics;
m_dvcsConvolCoeffFunctionResult = other.m_dvcsConvolCoeffFunctionResult;
......@@ -403,14 +402,21 @@ void DVMPProcessModule::initModule() {
ProcessModule<DVMPObservableKinematic, DVMPObservableResult>::initModule();
//evaluate internal variables
m_epsilon = 2 * m_xB * Constant::PROTON_MASS / sqrt(m_Q2);
m_y = m_Q2 / (2 * m_xB * Constant::PROTON_MASS * m_E);
double eps2 = m_epsilon * m_epsilon;
double epsroot = sqrt(1 + eps2);
double tfactor = -m_Q2 / (4 * m_xB * (1 - m_xB) + eps2);
m_tmin = tfactor * (2 * (1 - m_xB) * (1 - epsroot) + eps2); //TODO
m_tmax = tfactor * (2 * (1 - m_xB) * (1 + epsroot) + eps2); //TODO
m_xBmin = 2 * m_Q2 * m_E / Constant::PROTON_MASS / (4 * m_E * m_E - m_Q2);
double s = pow(Constant::PROTON_MASS, 2) + 2. * Constant::PROTON_MASS * m_E;
double E1cm = (s - m_Q2 - pow(Constant::PROTON_MASS, 2)) / (2 * sqrt(s));
double p1cm = sqrt(pow(E1cm, 2) + m_Q2);
double E3cm = (s + pow(MesonType(m_mesonType).getMass(), 2)
- pow(Constant::PROTON_MASS, 2)) / (2 * sqrt(s));
double p3cm = sqrt(pow(E3cm, 2) - pow(MesonType(m_mesonType).getMass(), 2));
m_tmin = pow(-m_Q2 - pow(MesonType(m_mesonType).getMass(), 2), 2) / (4 * s)
- pow(p1cm - p3cm, 2);
m_tmax = pow(-m_Q2 - pow(MesonType(m_mesonType).getMass(), 2), 2) / (4 * s)
- pow(p1cm + p3cm, 2);
}
void DVMPProcessModule::isModuleWellConfigured() {
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment