From c74f8696e87c1195468ce2a77b1bd1defc4004f9 Mon Sep 17 00:00:00 2001 From: Valerio Bertone Date: Tue, 11 Aug 2020 14:15:58 +0200 Subject: [PATCH 1/5] Update CMakeLists.txt --- CMakeLists.txt | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index 72a1cca8..8d09be10 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -56,6 +56,9 @@ find_package(ElementaryUtils REQUIRED) # find libraries: NumA++ find_package(NumA++ REQUIRED) +# find libraries: GSL +find_package(GSL REQUIRED) + # directories containing headers include_directories(${QT_INCLUDE_DIRS} ${SFML_INCLUDE_DIR} ${CLN_INCLUDE_DIR} ${ELEMENTARY_UTILS_INCLUDE_DIR} ${NUMA_INCLUDE_DIR}) include_directories(include) @@ -99,6 +102,8 @@ target_link_libraries( ${ELEMENTARY_UTILS_LIBRARIES} ${NUMA_LIBRARIES} + + ${GSL_LIBRARIES} ) # install -- GitLab From 87056a4db54e6162059777c64e145bcf2f08e35a Mon Sep 17 00:00:00 2001 From: Pawel Sznajder Date: Fri, 8 Jan 2021 16:29:19 +0100 Subject: [PATCH 2/5] add new SC module --- .../GPDSubtractionConstantDLMSTW21.h | 102 ++++++ .../GPDSubtractionConstantDLMSTW21Evolution.h | 111 +++++++ .../GPDSubtractionConstantDLMSTW21Replicas.h | 11 + .../GPDSubtractionConstantDLMSTW21.cpp | 199 +++++++++++ ...PDSubtractionConstantDLMSTW21Evolution.cpp | 311 ++++++++++++++++++ 5 files changed, 734 insertions(+) create mode 100644 include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21.h create mode 100644 include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21Evolution.h create mode 100644 include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21Replicas.h create mode 100644 src/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21.cpp create mode 100644 src/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21Evolution.cpp diff --git a/include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21.h b/include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21.h new file mode 100644 index 00000000..b7b0ad97 --- /dev/null +++ b/include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21.h @@ -0,0 +1,102 @@ +/* + * GPDSubtractionConstantDLMSTW21.h + * + * Created on: Jan 8, 2021 + * Author: Pawel Sznajder (NCBJ) + */ + +#ifndef GPD_SUBTRACTION_CONSTANT_DLMSTW21_H +#define GPD_SUBTRACTION_CONSTANT_DLMSTW21_H + +#include +#include +#include +#include +#include +#include + +namespace PARTONS { + +/** + * @class GPDSubtractionConstantDLMSTW21 + * + * @brief Implementation of DLMSTW subtraction constant. + * + * This module implements the subtraction constant based on tripole Ansatz with parameters + * extracted as described in 2021 paper by Dutrieux, Lorce, Moutarde, Sznajder, Trawinski and Wagner. + */ +class GPDSubtractionConstantDLMSTW21: public GPDSubtractionConstantModule { + +public: + + static const std::string PARAMETER_NAME_REPLICA; ///< Name of parameter to set replica index via xml file. + + /** + * Unique ID to automatically register the class in the registry. + */ + static const unsigned int classId; + + /** + * Constructor. + * @param className Name of this class. + */ + GPDSubtractionConstantDLMSTW21(const std::string &className); + + /** + * Destructor. + */ + virtual ~GPDSubtractionConstantDLMSTW21(); + + virtual void configure(const ElemUtils::Parameters ¶meters); + virtual GPDSubtractionConstantDLMSTW21* clone() const; + + /** + * Load parameters for a given replica index. + */ + void loadParameters(size_t replica, bool printInfo = true); + + /** + * Evaluate mean and uncertainty for a given vector of numbers. The procedure include removing of outliers. + */ + void getMeanAndUncertainty(const std::vector& v, double& mean, + double& unc) const; + +protected: + + /** Copy constructor. + * @param other Object to be copied. + */ + GPDSubtractionConstantDLMSTW21(const GPDSubtractionConstantDLMSTW21& other); + + virtual PhysicalType computeSubtractionConstant(); + +private: + + /** + * Evaluate mean from a given vector. + */ + double getMean(const std::vector& v) const; + + /** + * Evaluate sigma from a given vector. + */ + double getSigma(const std::vector& v) const; + + /** + * Remove outliers from a given vector using 3sigma rule. + */ + size_t removeOutliers(std::vector& v) const; + + size_t m_replica; ///< Replica index. + + double m_d1g; ///< d1 for gluons at initial scale. + double m_d1q; ///< d1 for light quarks at initial scale. + + double m_muF20; ///< Initial factorization scale. + double m_M; ///< Parameter of multipole Ansatz. + double m_alpha; ///< Parameter of multipole Ansatz. +}; + +} /* namespace PARTONS */ + +#endif /* GPD_SUBTRACTION_CONSTANT_DLMSTW21_H */ diff --git a/include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21Evolution.h b/include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21Evolution.h new file mode 100644 index 00000000..b43f6214 --- /dev/null +++ b/include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21Evolution.h @@ -0,0 +1,111 @@ +/* + * evolution.h + * + * Created on: Aug 6, 2019 + * Author: partons + */ + +#ifndef INCLUDE_EVOLUTION_H_ +#define INCLUDE_EVOLUTION_H_ + +#include +#include +#include + +namespace PARTONS { +namespace GPDSubtractionConstantDLMSTW21Evolution { + +/** + * Color factors. + */ +const double c_cf = 4. / 3.; +const double c_tf = 0.5; +const double c_ca = 3.; + +/** + * Quark masses (squared). + */ +const double c_m_u2 = pow(2.2 / 1.E3, 2); +const double c_m_d2 = pow(4.7 / 1.E3, 2); +const double c_m_s2 = pow(96. / 1.E3, 2); +const double c_m_c2 = pow(1.28, 2); +const double c_m_b2 = pow(4.18, 2); +const double c_m_t2 = pow(173.1, 2); + +/** + * Number of active flavors. + */ +size_t getNActiveFlavors(double muF2); + +/** + * Threshold for a given number of active flavors. + */ +double getThreshold2(size_t nf); + +/** + * Lambda_QCD. + */ +double lambdaQCD(size_t nf); + +/** + * Lambda_QCD. + */ +double lambdaQCD(double muF2); + +/** + * Beta coefficient. + */ +double beta(size_t i, size_t nf); + +/** + * Running coupling. + */ +double alphaS(double muF2); + +/** + * Anomalous dimension. + */ +double evolutionGamma(double n, int sign, size_t nf); + +/** + * Gluon evolution coefficients denoted as a_n^\pm. + */ +double evolutionA(size_t n, int sign, size_t nf); + +/** + * Evolve quark non-singlet ignoring quark's thresholds. + * dQuarkDiff is a diffrence of twho d^q. + */ +double evolveQuarkDiff(double muF2, double muF20, double dQuarkDiff, size_t n, + size_t nf); +/** + * Evolve quark singlet ignoring quark's thresholds. + * dQuarkSum is a sum of d^q. + */ +double evolveQuarkSum(double muF2, double muF20, double dQuarkSum, + double dGluon, size_t n, size_t nf); + +/** + * Evolve gluons ignoring quark's thresholds. + * dQuarkSum is a sum of d^q. + */ +double evolveGluon(double muF2, double muF20, double dQuarkSum, double dGluon, + size_t n, size_t nf); + +/** + * Evolve quarks and gluons ignoring quark's thresholds. + * dPartons is a vector of dGluon, dQuark1, dQuark2, ... + */ +std::vector evolveQuarkGluon(double mu2, double mu20, + std::vector dPartons, size_t n, size_t nf); + +/** + * Evolve quarks and gluons including quark's thresholds. + * dPartons is a vector of dGluon, dQuark1, dQuark2, ... + */ +std::vector evolveQuarkGluon(double mu2, double mu20, + std::vector dPartons, size_t n); +} +} + +#endif /* INCLUDE_EVOLUTION_H_ */ diff --git a/include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21Replicas.h b/include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21Replicas.h new file mode 100644 index 00000000..f8533d2e --- /dev/null +++ b/include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21Replicas.h @@ -0,0 +1,11 @@ +#include + +namespace PARTONS { +namespace GPDSubtractionConstantDLMSTW21Replicas { + +const size_t c_nReplicas = 1; + +const double c_Replicas[c_nReplicas][2] = { { 0., 0. } }; + +} +} diff --git a/src/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21.cpp b/src/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21.cpp new file mode 100644 index 00000000..88e00910 --- /dev/null +++ b/src/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21.cpp @@ -0,0 +1,199 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace PARTONS { + +const std::string GPDSubtractionConstantDLMSTW21::PARAMETER_NAME_REPLICA = + "replica"; + +const unsigned int GPDSubtractionConstantDLMSTW21::classId = + BaseObjectRegistry::getInstance()->registerBaseObject( + new GPDSubtractionConstantDLMSTW21( + "GPDSubtractionConstantDLMSTW21")); + +GPDSubtractionConstantDLMSTW21::GPDSubtractionConstantDLMSTW21( + const std::string& className) : + GPDSubtractionConstantModule(className), m_muF20(0.1), m_M(0.8), m_alpha( + 3.), m_replica(0), m_d1g(0.), m_d1q(0.) { + + //by default load for 0-th replica + loadParameters(0, false); +} + +GPDSubtractionConstantDLMSTW21::~GPDSubtractionConstantDLMSTW21() { +} + +GPDSubtractionConstantDLMSTW21* GPDSubtractionConstantDLMSTW21::clone() const { + return new GPDSubtractionConstantDLMSTW21(*this); +} + +GPDSubtractionConstantDLMSTW21::GPDSubtractionConstantDLMSTW21( + const GPDSubtractionConstantDLMSTW21& other) : + GPDSubtractionConstantModule(other), m_muF20(other.m_muF20), m_M( + other.m_M), m_alpha(other.m_alpha), m_replica(other.m_replica), m_d1g( + other.m_d1g), m_d1q(other.m_d1q) { +} + +void GPDSubtractionConstantDLMSTW21::configure( + const ElemUtils::Parameters ¶meters) { + + GPDSubtractionConstantModule::configure(parameters); + + if (parameters.isAvailable(PARAMETER_NAME_REPLICA)) { + loadParameters(parameters.getLastAvailable().toUInt()); + } +} + +void GPDSubtractionConstantDLMSTW21::loadParameters(size_t replica, + bool printInfo) { + + if (replica >= GPDSubtractionConstantDLMSTW21Replicas::c_nReplicas) { + throw ElemUtils::CustomException(getClassName(), __func__, + ElemUtils::Formatter() << "Illegal replica " << replica); + } + + m_d1g = GPDSubtractionConstantDLMSTW21Replicas::c_Replicas[replica][0]; + m_d1q = GPDSubtractionConstantDLMSTW21Replicas::c_Replicas[replica][1]; + + if (printInfo) { + info(__func__, + ElemUtils::Formatter() << "Parameters set for replica " + << replica); + } +} + +PhysicalType GPDSubtractionConstantDLMSTW21::computeSubtractionConstant() { + + //charges u,d,s,c,b,t + std::vector charges(6); + + charges.push_back(2 / 3.); //u + charges.push_back(-1 / 3.); //d + charges.push_back(-1 / 3.); //s + charges.push_back(2 / 3.); //c + charges.push_back(-1 / .3); //b + charges.push_back(2 / 3.); //t + + //d1 at initial scale + std::vector d1(4); + + d1.push_back(m_d1g); //g + d1.push_back(m_d1q); //u + d1.push_back(m_d1q); //d + d1.push_back(m_d1q); //s + + //evolve + std::vector d1Evolved = + GPDSubtractionConstantDLMSTW21Evolution::evolveQuarkGluon(m_MuF2, + m_muF20, d1, 1); + + //evaluate (gluons are skipped) + double dTermFormFactor = 0.; + + for (size_t i = 1; i < d1Evolved.size(); i++) + dTermFormFactor += pow(charges.at(i - 1), 2) * d1Evolved.at(i); + + //t-dependence + double tDep = pow(1. - m_t / pow(m_M, 2), -1 * m_alpha); + + //this takes into account that SC = 4 * dTermFormFactor + double sCdTermFactor = 4.; + + //return + return PhysicalType(sCdTermFactor * dTermFormFactor * tDep, + PhysicalUnit::NONE); +} + +double GPDSubtractionConstantDLMSTW21::getMean( + const std::vector& v) const { + + if (v.size() == 0) { + throw ElemUtils::CustomException(getClassName(), __func__, + "vector size is 0"); + } + + double mean = 0.; + + for (int i = 0; i < v.size(); i++) { + mean += v.at(i); + } + + return mean / v.size(); +} + +double GPDSubtractionConstantDLMSTW21::getSigma( + const std::vector& v) const { + + if (v.size() == 0) { + throw ElemUtils::CustomException(getClassName(), __func__, + "vector size is 0"); + } + + double mean = getMean(v); + + double sigma = 0.; + + for (int i = 0; i < v.size(); i++) { + sigma += pow(mean - v.at(i), 2); + } + + return sqrt(sigma / double(v.size())); +} + +size_t GPDSubtractionConstantDLMSTW21::removeOutliers( + std::vector& v) const { + + if (v.size() == 0) { + throw ElemUtils::CustomException(getClassName(), __func__, + "vector size is 0"); + } + + double meanData = getMean(v); + double sigmaData = getSigma(v); + + if (sigmaData == 0.) { + + warn(__func__, "sigma size is 0"); + return 0; + } + + std::vector result; + std::vector::iterator it; + size_t nRemoved = 0; + + for (it = v.begin(); it != v.end(); it++) { + + if (fabs((*it) - meanData) / sigmaData > 3.) { + nRemoved++; + } else { + result.push_back(*it); + } + } + + v = result; + + if (nRemoved != 0) + nRemoved += removeOutliers(v); + + return nRemoved; +} + +void GPDSubtractionConstantDLMSTW21::getMeanAndUncertainty( + const std::vector& v, double& mean, double& unc) const { + + std::vector vOutlierFree = v; + removeOutliers(vOutlierFree); + + mean = getMean(vOutlierFree); + unc = getSigma(vOutlierFree); +} + +} /* namespace PARTONS */ diff --git a/src/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21Evolution.cpp b/src/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21Evolution.cpp new file mode 100644 index 00000000..d613f74b --- /dev/null +++ b/src/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21Evolution.cpp @@ -0,0 +1,311 @@ +/* + * evolution.cpp + * + * Created on: Aug 6, 2019 + * Author: Paweł Sznajder and Arkadiusz P. Trawiński + */ + +#include +#include +#include +#include + +namespace PARTONS { +namespace GPDSubtractionConstantDLMSTW21Evolution { + +size_t getNActiveFlavors(double mu2) { + + if (mu2 >= c_m_u2 && mu2 < c_m_d2) { + return 1; + } else if (mu2 >= c_m_d2 && mu2 < c_m_s2) { + return 2; + } else if (mu2 >= c_m_s2 && mu2 < c_m_c2) { + return 3; + } else if (mu2 >= c_m_c2 && mu2 < c_m_b2) { + return 4; + } else if (mu2 >= c_m_b2 && mu2 < c_m_t2) { + return 5; + } else if (mu2 >= c_m_t2) { + return 6; + } + + return 0; +} + +double getThreshold2(size_t nf) { + + if (nf == 1) { + return c_m_u2; + } else if (nf == 2) { + return c_m_d2; + } else if (nf == 3) { + return c_m_s2; + } else if (nf == 4) { + return c_m_c2; + } else if (nf == 5) { + return c_m_b2; + } else if (nf == 6) { + return c_m_t2; + } + + return 0.; +} + +double lambdaQCD(size_t nf) { + + if (nf < 2) { + std::cout << "lambdaQCD(): nf cann't be smaller than 2" << std::endl; + exit(0); + } else if (nf > 6) { + std::cout << "lambdaQCD(): nf cann't be larger than 6" << std::endl; + exit(0); + } else if (nf == 2) { + return 0.14442; + } else { + double m = sqrt(getThreshold2(nf)); + double lamdaQCDsmall = lambdaQCD(nf - 1); + return lamdaQCDsmall * pow(lamdaQCDsmall / m, 2. / (33. - 2. * nf)); + } + + return 0; +} + +double lambdaQCD(double mu2) { + + size_t nf = getNActiveFlavors(mu2); + + if (nf < 2) { + std::cout << "lambdaQCD(): mu2 smaller than u quark mass" << std::endl; + exit(0); + } else + return lambdaQCD(nf); + + return 0.; +} + +double beta(size_t i, size_t nf) { + + switch (i) { + case 0: + return 11. - 2. * nf / 3.; + break; + + default: + std::cout << "beta(): beta_" << i << " undefined" << std::endl; + exit(0); + break; + } + + return 0.; +} + +double alphaS(double mu2) { + + if (mu2 <= 0.) { + std::cout << "alphaS(): illegal mu2 = " << mu2 << std::endl; + exit(0); + } + + size_t nf = getNActiveFlavors(mu2); + + return 4 * M_PI / (beta(0, nf) * log(mu2 / (pow(lambdaQCD(mu2), 2)))); +} + +double evolutionGamma(size_t n, int sign, size_t nf) { + + if (abs(sign) > 1) { + std::cout << "evolutionGamma(): illegal sign = " << sign << std::endl; + exit(0); + } + + double sum = 0.; + + for (size_t k = 2; k <= (n + 1); k++) { + sum += 1. / k; + } + + double gammaQQ = c_cf * (0.5 - (1. / ((n + 1.) * (n + 2.))) + 2 * sum); + + if (sign == 0) + return gammaQQ; + + double gammaQG = -1. * nf * c_tf * (n * n + 3 * n + 4.) + / (n * (n + 1.) * (n + 2.)); + double gammaGQ = -2. * c_cf * (n * n + 3 * n + 4.) + / ((n + 1.) * (n + 2.) * (n + 3.)); + double gammaGG = + c_ca + * (1. / 6. - 2. / (n * (n + 1)) - 2. / ((n + 2) * (n + 3)) + + 2 * sum) + (2. / 3.) * nf * c_tf; + return 0.5 + * (gammaQQ + gammaGG + + sign + * sqrt( + pow(gammaQQ - gammaGG, 2) + + 4 * gammaQG * gammaGQ)); +} + +double evolutionA(size_t n, int sign, size_t nf) { + + if (abs(sign) != 1) { + std::cout << "evolutionA(): illegal sign = " << sign << std::endl; + exit(0); + } + + double sum = 0.; + + for (size_t k = 2; k <= (n + 1); k++) { + sum += 1. / k; + } + + double gammaQQ = c_cf * (0.5 - (1. / ((n + 1.) * (n + 2.))) + 2 * sum); + double gammaQG = -1. * nf * c_tf * (n * n + 3 * n + 4.) + / (n * (n + 1.) * (n + 2.)); + + return 2. * nf / n * (evolutionGamma(n, sign, nf) - gammaQQ) / gammaQG; +} + +double evolveQuarkDiff(double mu2, double mu20, double dQuarkDiff, size_t n, + size_t nf) { + + return dQuarkDiff + * pow(alphaS(mu2) / alphaS(mu20), + 2 * evolutionGamma(n, 0, nf) / beta(0, nf)); +} + +double evolveQuarkSum(double mu2, double mu20, double dQuarkSum, double dGluon, + size_t n, size_t nf) { + + double aPlus = evolutionA(n, +1, nf); + double aMinus = evolutionA(n, -1, nf); + + double dPlus = (dGluon - aMinus * dQuarkSum / nf) / (aPlus - aMinus); + double dMinus = (-dGluon + aPlus * dQuarkSum / nf) / (aPlus - aMinus); + + return dPlus * nf + * pow(alphaS(mu2) / alphaS(mu20), + 2 * evolutionGamma(n, +1, nf) / beta(0, nf)) + + dMinus * nf + * pow(alphaS(mu2) / alphaS(mu20), + 2 * evolutionGamma(n, -1, nf) / beta(0, nf)); +} + +double evolveGluon(double mu2, double mu20, double dQuarkSum, double dGluon, + size_t n, size_t nf) { + + double aPlus = evolutionA(n, +1, nf); + double aMinus = evolutionA(n, -1, nf); + + double dPlus = (dGluon - aMinus * dQuarkSum / nf) / (aPlus - aMinus); + double dMinus = (-dGluon + aPlus * dQuarkSum / nf) / (aPlus - aMinus); + + return aPlus * dPlus + * pow(alphaS(mu2) / alphaS(mu20), + 2 * evolutionGamma(n, +1, nf) / beta(0, nf)) + + aMinus * dMinus + * pow(alphaS(mu2) / alphaS(mu20), + 2 * evolutionGamma(n, -1, nf) / beta(0, nf)); +} + +std::vector evolveQuarkGluon(double mu2, double mu20, + std::vector dPartons, size_t n, size_t nf) { + + if (nf < 2) { + std::cout << "evolveQuarkGluon(): nf cann't be smaller than 2" + << std::endl; + exit(0); + } else if (nf > 6) { + std::cout << "evolveQuarkGluon(): nf cann't be larger than 6" + << std::endl; + exit(0); + } else if (dPartons.size() != nf + 1) { + std::cout << "evolveQuarkGluon(): nf does't match size of dPartons " + << std::endl; + exit(0); + } + + /* Gluons are the first parton */ + double dGluon = dPartons.front(); + + /* Computing sum over quarks */ + double dQuarkSum = 0; + for (std::vector::iterator it = dPartons.begin() + 1; + it != dPartons.end(); it++) + dQuarkSum += *it; + + /* Vector to return */ + std::vector dResult; + + /* Computing evolved gluons and returing it */ + double dGluonEvolved = evolveGluon(mu2, mu20, dQuarkSum, dGluon, n, nf); + dResult.push_back(dGluonEvolved); + + /* Computing evolved sum over all quarks */ + double dQuarkSumEvolved = evolveQuarkSum(mu2, mu20, dQuarkSum, dGluon, n, + nf); + + /* Computing evolved difference of quarks in respect to the first one */ + std::vector dQuarkDiffEvolved; + double dQuark1 = *(dPartons.begin() + 1); + + for (std::vector::iterator it = dPartons.begin() + 2; + it != dPartons.end(); it++) + dQuarkDiffEvolved.push_back( + evolveQuarkDiff(mu2, mu20, *it - dQuark1, n, nf)); + + /* Computing evolved first quark and returing it*/ + double dQuark1Evolved = dQuarkSumEvolved; + for (std::vector::iterator it = dQuarkDiffEvolved.begin(); + it != dQuarkDiffEvolved.end(); it++) + dQuark1Evolved -= *it; + dQuark1Evolved /= nf; + + dResult.push_back(dQuark1Evolved); + + /* Returning of other quarks */ + for (std::vector::iterator it = dQuarkDiffEvolved.begin(); + it != dQuarkDiffEvolved.end(); it++) + dResult.push_back(dQuark1Evolved + *it); + + return dResult; +} + +std::vector evolveQuarkGluon(double mu2, double mu20, + std::vector dPartons, size_t n) { + size_t nf_mu2 = getNActiveFlavors(mu2); + size_t nf_mu20 = getNActiveFlavors(mu20); + + if (dPartons.size() != nf_mu20 + 1) { + std::cout << "evolveQuarkGluon(): At scale mu20 = " << mu20 + << " it is required of " << nf_mu20 + << " active quarks, but provided " << dPartons.size() - 1 + << " in dPartons." << std::endl; + exit(0); + } + + if (n % 2 == 0) { + std::cout << "evolveQuarkGluon(): n has to be odd" << std::endl; + exit(0); + } + + if (nf_mu2 > nf_mu20) { + double m2 = getThreshold2(nf_mu20 + 1); + std::vector dPartonsEvolved = evolveQuarkGluon(m2, mu20, + dPartons, n, nf_mu20); + dPartonsEvolved.push_back(0); /* contribution form heavier quark is generated radiatively */ + return evolveQuarkGluon(mu2, m2, dPartonsEvolved, n); + } + + if (nf_mu2 < nf_mu20) { + double m2 = getThreshold2(nf_mu2 + 1); + std::vector dPartonsEvolved = evolveQuarkGluon(m2, mu20, + dPartons, n); + dPartonsEvolved.pop_back(); /* removing contribution from the heaviest quark */ + return evolveQuarkGluon(mu2, m2, dPartonsEvolved, n, nf_mu2); + } + /* nf_mu2 == nf_mu20 */ + return evolveQuarkGluon(mu2, mu20, dPartons, n, nf_mu20); +} + +} +} -- GitLab From 348a6b236da0f4fd4a4677384c0c3bcba086179c Mon Sep 17 00:00:00 2001 From: Pawel Sznajder Date: Fri, 8 Jan 2021 22:48:13 +0100 Subject: [PATCH 3/5] add correct parameters for new SC module --- .../GPDSubtractionConstantDLMSTW21Replicas.h | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21Replicas.h b/include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21Replicas.h index f8533d2e..cdc54b8f 100644 --- a/include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21Replicas.h +++ b/include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21Replicas.h @@ -3,9 +3,23 @@ namespace PARTONS { namespace GPDSubtractionConstantDLMSTW21Replicas { -const size_t c_nReplicas = 1; +const size_t c_nReplicas = 101; -const double c_Replicas[c_nReplicas][2] = { { 0., 0. } }; +const double c_Replicas[c_nReplicas] = { -1.32998, -1.36677, 1.57145, -0.848952, + 2.51711, -3.17167, -2.65668, -0.193518, -3.60872, -1.2062, -1.6528, + -1.20646, -1.41473, -0.35459, -1.46504, -0.85388, 0.787904, -2.84513, + 0.99233, 9.5808, 0.988474, 1.55902, 0.426887, 1.30152, -1.81678, + 0.304013, -0.116567, 2.55685, -2.44012, -1.44651, -0.894461, 3.32256, + 1.2449, -3.42956, 0.989982, -0.00999596, 9.29666, -1.35486, 6.44742, + 1.63883, 0.828096, -3.28229, -1.20501, -0.848672, -2.61778, -3.8984, + -1.83257, -1.4171, 0.907111, 2.81445, 0.725484, 0.0538566, 2.78283, + 0.693457, -0.933821, -0.383909, -0.642304, -0.294465, 3.14944, -3.03039, + -1.12463, -1.96843, 1.66537, 6.48931, 5.55679, 8.89635, 1.32581, + -1.37508, -4.15628, -1.37918, -0.647732, -2.84596, -2.1652, -0.349579, + -2.25791, 0.936877, 2.19065, -1.9533, 0.380512, -1.73229, -0.2101, + -4.66385, -0.387866, 1.53858, -1.9789, -1.93822, -2.18833, -0.912989, + -0.213293, -1.49808, -1.31773, -1.69803, -1.02708, 6.28663, -0.734018, + -2.68933, 0.52901, -1.57246, -1.00882, 0.0262994, -0.945482 }; } } -- GitLab From 286a0df9c708c6a6dd5feb668ca969af5d9ab5f8 Mon Sep 17 00:00:00 2001 From: Pawel Sznajder Date: Fri, 8 Jan 2021 23:55:29 +0100 Subject: [PATCH 4/5] update --- .../GPDSubtractionConstantDLMSTW21.h | 8 ++-- .../GPDSubtractionConstantDLMSTW21Evolution.h | 6 +-- .../GPDSubtractionConstantDLMSTW21.cpp | 37 ++++++++++--------- ...PDSubtractionConstantDLMSTW21Evolution.cpp | 3 +- 4 files changed, 30 insertions(+), 24 deletions(-) diff --git a/include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21.h b/include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21.h index b7b0ad97..9058818e 100644 --- a/include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21.h +++ b/include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21.h @@ -8,13 +8,15 @@ #ifndef GPD_SUBTRACTION_CONSTANT_DLMSTW21_H #define GPD_SUBTRACTION_CONSTANT_DLMSTW21_H -#include -#include -#include #include #include #include +#include + +#include "../../modules/gpd_subtraction_constant/GPDSubtractionConstantModule.h" +#include "../../utils/type/PhysicalType.h" + namespace PARTONS { /** diff --git a/include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21Evolution.h b/include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21Evolution.h index b43f6214..7de22fa9 100644 --- a/include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21Evolution.h +++ b/include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21Evolution.h @@ -5,8 +5,8 @@ * Author: partons */ -#ifndef INCLUDE_EVOLUTION_H_ -#define INCLUDE_EVOLUTION_H_ +#ifndef INCLUDE_GPDSUBTRACTIONCONSTANTDLMSTW21EVOLUTION_H_ +#define INCLUDE_GPDSUBTRACTIONCONSTANTDLMSTW21EVOLUTION_H_ #include #include @@ -108,4 +108,4 @@ std::vector evolveQuarkGluon(double mu2, double mu20, } } -#endif /* INCLUDE_EVOLUTION_H_ */ +#endif /* INCLUDE_GPDSUBTRACTIONCONSTANTDLMSTW21EVOLUTION_H_ */ diff --git a/src/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21.cpp b/src/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21.cpp index 88e00910..17461273 100644 --- a/src/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21.cpp +++ b/src/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21.cpp @@ -1,11 +1,14 @@ +#include "../../../../include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21.h" + +#include "../../../../include/partons/BaseObjectRegistry.h" +#include "../../../../include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21Evolution.h" +#include "../../../../include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21Replicas.h" +#include "../../../../include/partons/utils/type/PhysicalUnit.h" + #include #include #include -#include -#include -#include -#include -#include + #include #include @@ -60,8 +63,8 @@ void GPDSubtractionConstantDLMSTW21::loadParameters(size_t replica, ElemUtils::Formatter() << "Illegal replica " << replica); } - m_d1g = GPDSubtractionConstantDLMSTW21Replicas::c_Replicas[replica][0]; - m_d1q = GPDSubtractionConstantDLMSTW21Replicas::c_Replicas[replica][1]; + m_d1g = 0.; + m_d1q = GPDSubtractionConstantDLMSTW21Replicas::c_Replicas[replica]; if (printInfo) { info(__func__, @@ -75,20 +78,20 @@ PhysicalType GPDSubtractionConstantDLMSTW21::computeSubtractionConstant( //charges u,d,s,c,b,t std::vector charges(6); - charges.push_back(2 / 3.); //u - charges.push_back(-1 / 3.); //d - charges.push_back(-1 / 3.); //s - charges.push_back(2 / 3.); //c - charges.push_back(-1 / .3); //b - charges.push_back(2 / 3.); //t + charges[0] = 2 / 3.; //u + charges[1] = -1 / 3.; //d + charges[2] = -1 / 3.; //s + charges[3] = 2 / 3.; //c + charges[4] = -1 / 3.; //b + charges[5] = 2 / 3.; //t //d1 at initial scale std::vector d1(4); - d1.push_back(m_d1g); //g - d1.push_back(m_d1q); //u - d1.push_back(m_d1q); //d - d1.push_back(m_d1q); //s + d1[0] = m_d1g; //g + d1[1] = m_d1q; //u + d1[2] = m_d1q; //d + d1[3] = m_d1q; //s //evolve std::vector d1Evolved = diff --git a/src/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21Evolution.cpp b/src/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21Evolution.cpp index d613f74b..6d755f5f 100644 --- a/src/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21Evolution.cpp +++ b/src/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21Evolution.cpp @@ -5,7 +5,8 @@ * Author: Paweł Sznajder and Arkadiusz P. Trawiński */ -#include +#include "../../../../include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21Evolution.h" + #include #include #include -- GitLab From 3d2a85481a69b6be1fc05462da37fb9061730d78 Mon Sep 17 00:00:00 2001 From: Pawel Sznajder Date: Sat, 9 Jan 2021 15:50:43 +0100 Subject: [PATCH 5/5] final parameters --- .../GPDSubtractionConstantDLMSTW21Replicas.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21Replicas.h b/include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21Replicas.h index cdc54b8f..64730961 100644 --- a/include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21Replicas.h +++ b/include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21Replicas.h @@ -19,7 +19,8 @@ const double c_Replicas[c_nReplicas] = { -1.32998, -1.36677, 1.57145, -0.848952, -2.25791, 0.936877, 2.19065, -1.9533, 0.380512, -1.73229, -0.2101, -4.66385, -0.387866, 1.53858, -1.9789, -1.93822, -2.18833, -0.912989, -0.213293, -1.49808, -1.31773, -1.69803, -1.02708, 6.28663, -0.734018, - -2.68933, 0.52901, -1.57246, -1.00882, 0.0262994, -0.945482 }; + -2.68933, 0.52901, -1.57246, -1.00882, 0.0262994, -0.945482 +}; } } -- GitLab