diff --git a/doc/bib/citations.bib b/doc/bib/citations.bib index 381486613271a849e16684e1372fa595fd034d58..51dc52da6d8eb5b877f86d4f05b04e6708cb7553 100644 --- a/doc/bib/citations.bib +++ b/doc/bib/citations.bib @@ -264,3 +264,45 @@ primaryClass = "hep-ph", SLACcitation = "%%CITATION = ARXIV:1512.06174;%%" } + +@article{Jakob:1997wg, + author = "Jakob, R. and Mulders, P. J. and Rodrigues, J.", + title = "{Modeling quark distribution and fragmentation + functions}", + journal = "Nucl. Phys.", + volume = "A626", + year = "1997", + pages = "937-965", + doi = "10.1016/S0375-9474(97)00588-5", + eprint = "hep-ph/9704335", + archivePrefix = "arXiv", + primaryClass = "hep-ph", + reportNumber = "NIKHEF-97-018, VUTH-97-7", + SLACcitation = "%%CITATION = HEP-PH/9704335;%%" +} + +@article{Hwang:2007tb, + author = "Hwang, D. S. and Mueller, Dieter", + title = "{Implication of the overlap representation for modelling + generalized parton distributions}", + journal = "Phys. Lett.", + volume = "B660", + year = "2008", + pages = "350-359", + doi = "10.1016/j.physletb.2008.01.014", + eprint = "0710.1567", + archivePrefix = "arXiv", + primaryClass = "hep-ph", + SLACcitation = "%%CITATION = ARXIV:0710.1567;%%" +} + +@article{Muller:2014tqa, + author = "Müller, Dieter and Hwang, Dae Sung", + title = "{The concept of phenomenological light-front wave + functions -- Regge improved diquark model predictions}", + year = "2014", + eprint = "1407.1655", + archivePrefix = "arXiv", + primaryClass = "hep-ph", + SLACcitation = "%%CITATION = ARXIV:1407.1655;%%" +} \ No newline at end of file diff --git a/include/partons/modules/gpd/GPDHM18.h b/include/partons/modules/gpd/GPDHM18.h new file mode 100644 index 0000000000000000000000000000000000000000..cc02bcfe782ff7cc1900bfdf80b692d1289cd1ab --- /dev/null +++ b/include/partons/modules/gpd/GPDHM18.h @@ -0,0 +1,180 @@ +#ifndef GPD_HM18_H +#define GPD_HM18_H + +/** + * @file GPDHM18.h + * @author Arkadiusz P. Trawiński (CEA) + * @date July 2, 2018 + * @version 1.0 + */ + +#include +#include +#include + +#include "../MathIntegratorModule.h" +#include "GPDModule.h" + +namespace PARTONS { + +/** + * @class GPDHM18 + * + * The proposition of the scalar di-quark model (SDQM) in the context + * of GPDs back to the Jakob-Mulders-Rodrigues paper hep-ph/9704335 + * @cite Jakob:1997wg. + * + * Later, the model has been studied by Dae Sung Hwang and Dieter Mueller. + * The short description of the GPDs derived from this model is given in + * arXiv:0710.1567 @cite Hwang:2007tb. In the Ref. arXiv:1407.16655 + * @cite Muller:2014tqa full description is given. + * + * Available GPD types: H, Ht, E, Et. + */ + +class GPDHM18: public GPDModule, public MathIntegratorModule { +public: + static const std::string PARAMETER_NAME_HM18MODEL_M; ///< Name of parameter to set #m_M via configure() + static const std::string PARAMETER_NAME_HM18MODEL_m; ///< Name of parameter to set #m_m via configure() + static const std::string PARAMETER_NAME_HM18MODEL_lambda; ///< Name of parameter to set #m_lambda via configure() + static const std::string PARAMETER_NAME_HM18MODEL_p; ///< Name of parameter to set #m_p via configure() + /** + * ID assigned by BaseObjectRegistry. + */ + static const unsigned int classId; + /** + * Default constructor. + * @param className Name of class. + */ + GPDHM18(const std::string& className); + /** + * Destructor. + */ + virtual ~GPDHM18(); + virtual GPDHM18* clone() const; + virtual void resolveObjectDependencies(); + virtual void configure(const ElemUtils::Parameters ¶meters); +protected: + /** + * Copy constructor. + * @param other Object to be copied. + */ + GPDHM18(const GPDHM18& other); + virtual void isModuleWellConfigured(); + virtual void initModule(); + virtual PartonDistribution computeH(); + virtual PartonDistribution computeHt(); + virtual PartonDistribution computeE(); + virtual PartonDistribution computeEt(); + +private: + double m_M; ///< Mass of the proton; + double m_m; ///< Mass of the active quark. + double m_lambda; ///< Mass of the spectator + double m_p; ///< Parameter controlling the power. + double m_N; ///< Normalization of the GPDs. + + /** + * Function setting the normalization of the GPDs @param m_N + * according to the Eq. (22) Ref. @cite Hwang:2007tb. * + */ + void Normalize(); + + /** + * Integral given in Eq. (22) Ref. @cite Hwang:2007tb. + * + * @param y + * @param par optional parameters of integration, not used. + * @return 1/m_N + */ + double IntNorm(double y, std::vector par); + NumA::FunctionType1D* m_pint_IntNorm; ///< Functor related to IntNorm. + + /** %Double distribution function for GPD E. Definition can be found + * in Eqs. (18) and (19) in Ref. @cite Hwang:2007tb. + * + * @param y %Double distribution parameter y. + * @param z %Double distribution parameter z. + * @param t %Double distribution parameter t. + * @return GPD E + */ + double DD_E(double y, double z, double t); + double IntE(double y, std::vector par); ///< Integrand of GPD E for xi <> 0. + double IntE0(double z, std::vector par); ///< Integrand of GPD E for xi == 0. + NumA::FunctionType1D* m_pint_IntE; ///< Functor related to IntE. + NumA::FunctionType1D* m_pint_IntE0; ///< Functor related to IntE0. + + /** + * %Double distribution function for GPD H. Definition can be found + * in Eqs. (21) and (19) in Ref. @cite Hwang:2007tb. + * + * @param y %Double distribution parameter y. + * @param z %Double distribution parameter z. + * @param t %Double distribution parameter t. + * @return GPD H + */ + double DD_H(double y, double z, double t); + double IntH(double y, std::vector par); ///< Integrand of GPD H for xi <> 0. + double IntH0(double z, std::vector par); ///< Integrand of GPD H for xi == 0. + NumA::FunctionType1D* m_pint_IntH; ///< Functor related to IntH. + NumA::FunctionType1D* m_pint_IntH0; ///< Functor related to IntH0. + + /** + * %Double distribution function for GPD Ht. Definition can be found + * in Eqs. (4.23) and (4.26) in Ref. @cite Muller:2014tqa. + * + * @param y %Double distribution parameter y. + * @param z %Double distribution parameter z. + * @param t %Double distribution parameter t. + * @return GPD Ht + */ + double DD_Ht(double y, double z, double t); + double IntHt(double y, std::vector par); ///< Integrand of GPD Ht for xi <> 0. + double IntHt0(double z, std::vector par); ///< Integrad of GPD Ht for xi == 0. + NumA::FunctionType1D* m_pint_IntHt; ///< Functor related to IntHt. + NumA::FunctionType1D* m_pint_IntHt0; ///< Functor related to IntHt0. + + /** %Double distribution function for GPD Et. Definition can be found + * in Eqs. (4.23) and (4.27) in Ref. @cite Muller:2014tqa. + * + * @param y %Double distribution parameter y. + * @param z %Double distribution parameter z. + * @param t %Double distribution parameter t. + * @return GPD Et + */ + + double DD_Et(double y, double z, double t); + double IntEt(double y, std::vector par); ///< Integrand of GPD Et for xi <> 0. + NumA::FunctionType1D* m_pint_IntEt; ///< Functor related to IntE for xi <> 0. + + /** + * Evaluates GPDs using %double distribution function + * by integrating p_fun0 or p_fun function depending on kinematics. + * + * @param x + * @param p_fun0 integrate this function if xi == 0 + * @param p_fun integrate this function if xi <> 0 + * @return computed GPD + */ + double evaluate(double x, NumA::FunctionType1D* p_fun0, + NumA::FunctionType1D* p_fun); + + /** + * Computes GPDs using GPDHM18::evaluate function + * and creates appropriate parton distributions. + * + * @param p_fun0 integrate this function if xi == 0 + * @param p_fun integrate this function if xi <> 0 + * @return parton distributions + */ + PartonDistribution compute(NumA::FunctionType1D* p_fun0, + NumA::FunctionType1D* p_fun); + + void initializeFunctorsForIntegrations(); ///< Initialize functors. + void deleteFunctorsForIntegrations(); ///< Delete functors. + +}; + +} /* namespace PARTONS */ + +#endif /* GPD_HM18_H */ diff --git a/src/partons/modules/gpd/GPDHM18.cpp b/src/partons/modules/gpd/GPDHM18.cpp new file mode 100644 index 0000000000000000000000000000000000000000..81f034753ad95c5dc21fb3ad26ce29f5d8c076b0 --- /dev/null +++ b/src/partons/modules/gpd/GPDHM18.cpp @@ -0,0 +1,376 @@ +#include "../../../../include/partons/modules/gpd/GPDHM18.h" + +#include +#include +#include +#include +#include +#include +#include + +#include "../../../../include/partons/beans/gpd/GPDType.h" +#include "../../../../include/partons/beans/parton_distribution/GluonDistribution.h" +#include "../../../../include/partons/beans/parton_distribution/QuarkDistribution.h" +#include "../../../../include/partons/beans/QuarkFlavor.h" +#include "../../../../include/partons/BaseObjectRegistry.h" +#include "../../../../include/partons/FundamentalPhysicalConstants.h" + +namespace PARTONS { + +const unsigned int GPDHM18::classId = + BaseObjectRegistry::getInstance()->registerBaseObject( + new GPDHM18("GPDHM18")); + +const std::string GPDHM18::PARAMETER_NAME_HM18MODEL_M = "HM18MODEL_M"; +const std::string GPDHM18::PARAMETER_NAME_HM18MODEL_m = "HM18MODEL_m"; +const std::string GPDHM18::PARAMETER_NAME_HM18MODEL_lambda = "HM18MODEL_lambda"; +const std::string GPDHM18::PARAMETER_NAME_HM18MODEL_p = "HM18MODEL_p"; + +GPDHM18::GPDHM18(const std::string &className) : + GPDModule(className), MathIntegratorModule() { + + initializeFunctorsForIntegrations(); + + m_MuF2_ref = 4.; + + m_M = Constant::PROTON_MASS; + m_m = 0.45; + m_lambda = 0.75; + m_p = 1.; + + /** + * Value of m_N is calculated in GPDHM18::configure using Normalize() function. + * In this place Normalization() function cann't be run, since it is using + * integration mechanism, which will be initiated not sooner that in GPDHM18::configure. + * In order to not to leave the value of m_N not initialize the value returned + * by Normalization() with initial values of model parameters has been used. + */ + m_N = 0.048900116463331; + + //relate a specific GPD type with the appropriate function + m_listGPDComputeTypeAvailable.insert( + std::make_pair(GPDType::H, &GPDModule::computeH)); + m_listGPDComputeTypeAvailable.insert( + std::make_pair(GPDType::Ht, &GPDModule::computeHt)); + m_listGPDComputeTypeAvailable.insert( + std::make_pair(GPDType::E, &GPDModule::computeE)); + m_listGPDComputeTypeAvailable.insert( + std::make_pair(GPDType::Et, &GPDModule::computeEt)); +} + +GPDHM18::GPDHM18(const GPDHM18& other) : + GPDModule(other), MathIntegratorModule(other) { + + initializeFunctorsForIntegrations(); + + m_M = other.m_M; + m_m = other.m_m; + m_lambda = other.m_lambda; + m_p = other.m_p; + m_N = other.m_N; + +} + +GPDHM18::~GPDHM18() { + deleteFunctorsForIntegrations(); +} + +void GPDHM18::initializeFunctorsForIntegrations() { + + m_pint_IntNorm = NumA::Integrator1D::newIntegrationFunctor(this, + &GPDHM18::IntNorm); + + m_pint_IntE = NumA::Integrator1D::newIntegrationFunctor(this, + &GPDHM18::IntE); + + m_pint_IntE0 = NumA::Integrator1D::newIntegrationFunctor(this, + &GPDHM18::IntE0); + + m_pint_IntH = NumA::Integrator1D::newIntegrationFunctor(this, + &GPDHM18::IntH); + + m_pint_IntH0 = NumA::Integrator1D::newIntegrationFunctor(this, + &GPDHM18::IntH0); + + m_pint_IntHt = NumA::Integrator1D::newIntegrationFunctor(this, + &GPDHM18::IntHt); + + m_pint_IntHt0 = NumA::Integrator1D::newIntegrationFunctor(this, + &GPDHM18::IntHt0); + + m_pint_IntEt = NumA::Integrator1D::newIntegrationFunctor(this, + &GPDHM18::IntEt); +} + +void GPDHM18::deleteFunctorsForIntegrations() { + + delete m_pint_IntNorm; + delete m_pint_IntE; + delete m_pint_IntE0; + delete m_pint_IntH; + delete m_pint_IntH0; + delete m_pint_IntHt; + delete m_pint_IntHt0; + delete m_pint_IntEt; +} + +GPDHM18* GPDHM18::clone() const { + return new GPDHM18(*this); +} +void GPDHM18::resolveObjectDependencies() { + setIntegrator(NumA::IntegratorType1D::DEXP); +} +void GPDHM18::configure(const ElemUtils::Parameters ¶meters) { + MathIntegratorModule::configureIntegrator(parameters); + GPDModule::configure(parameters); + + if (parameters.isAvailable(GPDHM18::PARAMETER_NAME_HM18MODEL_M)) { + + m_M = parameters.getLastAvailable().toDouble(); + info(__func__, + ElemUtils::Formatter() << "Parameter " + << GPDHM18::PARAMETER_NAME_HM18MODEL_M << " changed to " + << m_M); + } + + if (parameters.isAvailable(GPDHM18::PARAMETER_NAME_HM18MODEL_m)) { + + m_m = parameters.getLastAvailable().toDouble(); + info(__func__, + ElemUtils::Formatter() << "Parameter " + << GPDHM18::PARAMETER_NAME_HM18MODEL_m << " changed to " + << m_m); + } + + if (parameters.isAvailable(GPDHM18::PARAMETER_NAME_HM18MODEL_lambda)) { + m_lambda = parameters.getLastAvailable().toDouble(); + info(__func__, + ElemUtils::Formatter() << "Parameter " + << GPDHM18::PARAMETER_NAME_HM18MODEL_lambda + << " changed to " << m_lambda); + } + + if (parameters.isAvailable(GPDHM18::PARAMETER_NAME_HM18MODEL_p)) { + m_p = parameters.getLastAvailable().toDouble(); + info(__func__, + ElemUtils::Formatter() << "Parameter " + << GPDHM18::PARAMETER_NAME_HM18MODEL_p << " changed to " + << m_p); + } + + Normalize(); +} +void GPDHM18::isModuleWellConfigured() { + GPDModule::isModuleWellConfigured(); +} +void GPDHM18::initModule() { + GPDModule::initModule(); +} + +double GPDHM18::IntNorm(double y, std::vector par) { + + double m2 = pow(m_m, 2); + double M2 = pow(m_M, 2); + double lambda2 = pow(m_lambda, 2); + double coeff = sqrt(Constant::PI) / (4 * M2 * m_p) * tgamma(1 + m_p) + / tgamma(1.5 + m_p); + + double Num = m2 * (1 + 2 * m_p - y) + 4 * m_m * m_M * m_p * y + + y * (M2 * (y - 1 + 2 * m_p * y) + lambda2); + double Den = m2 / M2 - y + lambda2 / M2 * y / (1 - y); + + return coeff * Num / pow(Den, 2 * m_p + 1); +} + +void GPDHM18::Normalize() { + //set variables for integrations + std::vector parameters; + + double result = integrate(m_pint_IntNorm, 0, 1, parameters); + m_N = 1 / result; + + info(__func__, ElemUtils::Formatter() << "Normalization set to " << m_N); +} + +double GPDHM18::DD_E(double y, double z, double t) { + double M2 = pow(m_M, 2); + double m2 = pow(m_m, 2); + double Num = (m_m / m_M + y) * pow(pow(1 - y, 2) - pow(z, 2), m_p); + double Den = (1 - y) * m2 / M2 + y * pow(m_lambda, 2) / M2 - y * (1 - y) + - (pow(1 - y, 2) - pow(z, 2)) * t / 4 / M2; + return m_N * Num / pow(Den, 2 * m_p + 1); + +} + +double GPDHM18::IntE0(double z, std::vector par) { + double x = par[0]; + return (1 - x) * DD_E(x, z, m_t); + +} + +double GPDHM18::IntE(double y, std::vector par) { + double x = par[0]; + return (1 - x) / m_xi * DD_E(y, (x - y) / m_xi, m_t); + +} + +double GPDHM18::DD_H(double y, double z, double t) { + double M2 = pow(m_M, 2); + double m2 = pow(m_m, 2); + double Num1 = pow(pow(1 - y, 2) - pow(z, 2), m_p); + double Num2 = (-y + y * pow(m_lambda, 2) / M2 + (2 - y) * m2 / M2) * Num1; + double Den = (1 - y) * m2 / M2 + y * pow(m_lambda, 2) / M2 - y * (1 - y) + - (pow(1 - y, 2) - pow(z, 2)) * t / 4 / M2; + return m_N * (1 - 2 * m_p) / (4 * m_p) * Num1 / pow(Den, 2 * m_p) + + m_N * Num2 / 2 / pow(Den, 2 * m_p + 1); +} + +double GPDHM18::IntH0(double z, std::vector par) { + double x = par[0]; + return DD_H(x, z, m_t) + x * DD_E(x, z, m_t); +} + +double GPDHM18::IntH(double y, std::vector par) { + double x = par[0]; + return 1 / m_xi * DD_H(y, (x - y) / m_xi, m_t) + + x / m_xi * DD_E(y, (x - y) / m_xi, m_t); +} + +double GPDHM18::DD_Ht(double y, double z, double t) { + double M2 = pow(m_M, 2); + double m2 = pow(m_m, 2); + double lambda2 = pow(m_lambda, 2); + double Num = pow(pow(1 - y, 2) - pow(z, 2), m_p); + double Den = (1 - y) * m2 / M2 + y * lambda2 / M2 - y * (1 - y) + - (pow(1 - y, 2) - pow(z, 2)) * t / 4 / M2; + return m_N / 2 * (m2 / M2 * y - y * lambda2 / M2 + y+2 * y * m_m / m_M ) + * Num / pow(Den, 2 * m_p + 1) + - m_N * (1 - 2 * m_p) / (4 * m_p) * Num / pow(Den, 2 * m_p); +} + +double GPDHM18::IntHt0(double z, std::vector par) { + double x = par[0]; + return DD_Ht(x, z, m_t); +} + +double GPDHM18::IntHt(double y, std::vector par) { + double x = par[0]; + return 1 / m_xi * DD_Ht(y, (x - y) / m_xi, m_t); +} + +double GPDHM18::DD_Et(double y, double z, double t) { + double M2 = pow(m_M, 2); + double m2 = pow(m_m, 2); + double Bracket = pow(1 - y, 2) - pow(z, 2) + + (1 - y - z / m_xi) * (m_m / m_M + y); + double Num = Bracket * pow(pow(1 - y, 2) - pow(z, 2), m_p); + double Den = (1 - y) * m2 / M2 + y * pow(m_lambda, 2) / M2 - y * (1 - y) + - (pow(1 - y, 2) - pow(z, 2)) * t / 4 / M2; + return m_N * Num / pow(Den, 2 * m_p + 1); +} + +double GPDHM18::IntEt(double y, std::vector par) { + double x = par[0]; + return 1 / m_xi * DD_Et(y, (x - y) / m_xi, m_t); +} + +double GPDHM18::evaluate(double x, NumA::FunctionType1D* p_fun0, + NumA::FunctionType1D* p_fun) { + //set variables for integrations + std::vector parameters; + parameters.push_back(x); + + //checking kinetic region + if (fabs(x) > 1 || fabs(m_xi) > 1) + throw ElemUtils::CustomException(getClassName(), __func__, + ElemUtils::Formatter() + << "Calculation not in the kinetic domain."); + + //calculate GPD for x < - xi + if (x < -fabs(m_xi)) + return 0; + + //calculate GPD for xi == 0 + if (m_xi == 0) + if (p_fun0 == 0) + throw ElemUtils::CustomException(getClassName(), __func__, + ElemUtils::Formatter() << "Integrand not set."); + else + return integrate(p_fun0, -1 + x, 1 - x, parameters); + + //calculate GPD for xi <> 0 and x > xi + if (x > fabs(m_xi)) + if (p_fun == 0) + throw ElemUtils::CustomException(getClassName(), __func__, + ElemUtils::Formatter() << "Integrand not set."); + else + return integrate(p_fun, (x - m_xi) / (1 - m_xi), + (x + m_xi) / (1 + m_xi), parameters); + + //calculate GPD for xi <> 0 and - xi < x < xi + if (x <= fabs(m_xi) && x >= -fabs(m_xi)) + if (p_fun == 0) + throw ElemUtils::CustomException(getClassName(), __func__, + ElemUtils::Formatter() << "Integrand not set."); + else + return integrate(p_fun, 0., (x + m_xi) / (1 + m_xi), parameters); + + return 0; +} + +PartonDistribution GPDHM18::compute(NumA::FunctionType1D* p_fun0, + NumA::FunctionType1D* p_fun) { + //variables + double aValPx = GPDHM18::evaluate(m_x, p_fun0, p_fun); + double aValMx = GPDHM18::evaluate(-m_x, p_fun0, p_fun); + double Sea = 0; + double g = 0; + + //store + PartonDistribution partonDistribution; + QuarkDistribution quarkDistribution(QuarkFlavor::UNDEFINED); + GluonDistribution gluonDistribution(g); + + quarkDistribution.setQuarkDistribution(aValPx + Sea); + quarkDistribution.setQuarkDistributionPlus(aValPx - aValMx + 2 * Sea); + quarkDistribution.setQuarkDistributionMinus(aValPx + aValMx); + + partonDistribution.setGluonDistribution(gluonDistribution); + partonDistribution.addQuarkDistribution(quarkDistribution); + + //return + return partonDistribution; +} + +PartonDistribution GPDHM18::computeE() { + //compute and return parton distribution for GPH E + return GPDHM18::compute(m_pint_IntE0, m_pint_IntE); +} + +PartonDistribution GPDHM18::computeH() { + //compute and return parton distribution for GPH H + return GPDHM18::compute(m_pint_IntH0, m_pint_IntH); +} + +PartonDistribution GPDHM18::computeHt() { + //compute and return parton distribution for GPH Ht + return GPDHM18::compute(m_pint_IntHt0, m_pint_IntHt); +} + +PartonDistribution GPDHM18::computeEt() { + //No result of GPH Et for xi == 0. + PartonDistribution no_result; + if (m_xi == 0) { + warn(__func__, + ElemUtils::Formatter() + << "GPD Et is divergent for xi==0 in the SDQM by Hwang-Mueller."); + warn(__func__, + ElemUtils::Formatter() + << "No quark GPD Et result will be returned. "); + return no_result; + } else + //compute and return parton distribution for GPH Et for xi <> 0. + return GPDHM18::compute(0, m_pint_IntEt); +} + +} /* namespace PARTONS */