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

Merge branch 'subtraction_constant_0121' into 'master'

Subtraction constant 0121

See merge request !36
parents 865213f4 794cd55b
/*
* GPDSubtractionConstantDLMSTW21.h
*
* Created on: Jan 8, 2021
* Author: Pawel Sznajder (NCBJ)
*/
#ifndef GPD_SUBTRACTION_CONSTANT_DLMSTW21_H
#define GPD_SUBTRACTION_CONSTANT_DLMSTW21_H
#include <stddef.h>
#include <string>
#include <vector>
#include <ElementaryUtils/parameters/Parameters.h>
#include "../../modules/gpd_subtraction_constant/GPDSubtractionConstantModule.h"
#include "../../utils/type/PhysicalType.h"
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 &parameters);
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<double>& v, double& mean,
double& unc) const;
protected:
/** Copy constructor.
* @param other Object to be copied.
*/
GPDSubtractionConstantDLMSTW21(const GPDSubtractionConstantDLMSTW21& other);
virtual PhysicalType<double> computeSubtractionConstant();
private:
/**
* Evaluate mean from a given vector.
*/
double getMean(const std::vector<double>& v) const;
/**
* Evaluate sigma from a given vector.
*/
double getSigma(const std::vector<double>& v) const;
/**
* Remove outliers from a given vector using 3sigma rule.
*/
size_t removeOutliers(std::vector<double>& 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 */
/*
* evolution.h
*
* Created on: Aug 6, 2019
* Author: partons
*/
#ifndef INCLUDE_GPDSUBTRACTIONCONSTANTDLMSTW21EVOLUTION_H_
#define INCLUDE_GPDSUBTRACTIONCONSTANTDLMSTW21EVOLUTION_H_
#include <stddef.h>
#include <cmath>
#include <vector>
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<double> evolveQuarkGluon(double mu2, double mu20,
std::vector<double> 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<double> evolveQuarkGluon(double mu2, double mu20,
std::vector<double> dPartons, size_t n);
}
}
#endif /* INCLUDE_GPDSUBTRACTIONCONSTANTDLMSTW21EVOLUTION_H_ */
#include <stddef.h>
namespace PARTONS {
namespace GPDSubtractionConstantDLMSTW21Replicas {
const size_t c_nReplicas = 101;
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
};
}
}
#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 <ElementaryUtils/logger/CustomException.h>
#include <ElementaryUtils/parameters/GenericType.h>
#include <ElementaryUtils/string_utils/Formatter.h>
#include <cmath>
#include <iterator>
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 &parameters) {
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 = 0.;
m_d1q = GPDSubtractionConstantDLMSTW21Replicas::c_Replicas[replica];
if (printInfo) {
info(__func__,
ElemUtils::Formatter() << "Parameters set for replica "
<< replica);
}
}
PhysicalType<double> GPDSubtractionConstantDLMSTW21::computeSubtractionConstant() {
//charges u,d,s,c,b,t
std::vector<double> charges(6);
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<double> d1(4);
d1[0] = m_d1g; //g
d1[1] = m_d1q; //u
d1[2] = m_d1q; //d
d1[3] = m_d1q; //s
//evolve
std::vector<double> 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<double>(sCdTermFactor * dTermFormFactor * tDep,
PhysicalUnit::NONE);
}
double GPDSubtractionConstantDLMSTW21::getMean(
const std::vector<double>& 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<double>& 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<double>& 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<double> result;
std::vector<double>::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<double>& v, double& mean, double& unc) const {
std::vector<double> vOutlierFree = v;
removeOutliers(vOutlierFree);
mean = getMean(vOutlierFree);
unc = getSigma(vOutlierFree);
}
} /* namespace PARTONS */
/*
* evolution.cpp
*
* Created on: Aug 6, 2019
* Author: Paweł Sznajder and Arkadiusz P. Trawiński
*/
#include "../../../../include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21Evolution.h"
#include <cstdlib>
#include <iostream>
#include <iterator>
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);