add new SC module

include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21.h

+/*
+ * GPDSubtractionConstantDLMSTW21.h
+ *
+ *  Created on: Jan 8, 2021
+ *      Author: Pawel Sznajder (NCBJ)
+ */
+
+#ifndef GPD_SUBTRACTION_CONSTANT_DLMSTW21_H
+#define GPD_SUBTRACTION_CONSTANT_DLMSTW21_H
+
+#include <ElementaryUtils/parameters/Parameters.h>
+#include <include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantModule.h>
+#include <include/partons/utils/type/PhysicalType.h>
+#include <stddef.h>
+#include <string>
+#include <vector>
+
+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 */

include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21Evolution.h

+/*
+ * evolution.h
+ *
+ *  Created on: Aug 6, 2019
+ *      Author: partons
+ */
+
+#ifndef INCLUDE_EVOLUTION_H_
+#define INCLUDE_EVOLUTION_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_EVOLUTION_H_ */

include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21Replicas.h

+#include <stddef.h>
+
+namespace PARTONS {
+namespace GPDSubtractionConstantDLMSTW21Replicas {
+
+const size_t c_nReplicas = 1;
+
+const double c_Replicas[c_nReplicas][2] = { { 0., 0. } };
+
+}
+}

src/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21.cpp

+#include <ElementaryUtils/logger/CustomException.h>
+#include <ElementaryUtils/parameters/GenericType.h>
+#include <ElementaryUtils/string_utils/Formatter.h>
+#include <include/partons/BaseObjectRegistry.h>
+#include <include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21.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 <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 = GPDSubtractionConstantDLMSTW21Replicas::c_Replicas[replica][0];
+	m_d1q = GPDSubtractionConstantDLMSTW21Replicas::c_Replicas[replica][1];
+
+	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.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<double> 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<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 */

src/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantDLMSTW21Evolution.cpp

+/*
+ * 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);
+
+	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);