Skip to content
GitLab
Commit 2fd33ad4 authored by Cédric Mezrag's avatar Cédric Mezrag
Browse files

Merge branch 'APFELevolution_check' of... 

Merge branch 'APFELevolution_check' of git@drf-gitlab.cea.fr:partons/core/partons into APFELevolution_check
parents 063d3f08 1e6a0780
Hide whitespace changes
Inline Side-by-side
include/partons/modules/evolution/gpd/GPDEvolutionApfel.h
View file @ 2fd33ad4
......@@ -81,8 +81,10 @@ public:
double getTabLowerBound() const;
double getTabUpperBound() const;
int getTabInterDegree() const;
double getPreviousXi() const;
std::shared_ptr<apfel::Grid> getGrid() const;
std::function<double(double const&)> getAlphas() const;
std::shared_ptr<apfel::TabulateObject<apfel::Set<apfel::Distribution>>> getTabulatedGPDs() const;
protected:
......@@ -110,9 +112,9 @@ private:
int m_tabInterDegree;
double m_xi_prev;
std::unique_ptr<apfel::Grid> m_g;
std::shared_ptr<apfel::Grid> m_g;
std::function<double(double const&)> m_as;
std::unique_ptr<apfel::TabulateObject<apfel::Set<apfel::Distribution>>> m_tabulatedGpds;
std::shared_ptr<apfel::TabulateObject<apfel::Set<apfel::Distribution>>> m_tabulatedGpds;
};
} /* namespace PARTONS */
......
src/partons/modules/evolution/gpd/GPDEvolutionApfel.cpp
View file @ 2fd33ad4
......@@ -19,7 +19,7 @@ const unsigned int GPDEvolutionApfel::classId = BaseObjectRegistry::getInstance(
GPDEvolutionApfel::GPDEvolutionApfel(const std::string &className) :
GPDEvolutionModule(className),
m_subgridNodes({100, 60, 50, 50}),
m_subgridLowerBounds({0.0001, 0.1, 0.6, 0.8}),
m_subgridLowerBounds({0.0000001, 0.1, 0.6, 0.8}),
m_subgridInterDegrees({3, 3, 3, 3}),
m_tabNodes(100),
m_tabLowerBound(1),
......@@ -100,32 +100,35 @@ void GPDEvolutionApfel::prepareSubModules(const std::map<std::string, BaseObject
}
PartonDistribution GPDEvolutionApfel::compute(GPDModule* pGPDModule, const GPDType::Type &type) {
// Set initial scale
m_MuF2_ref = pGPDModule->getMuF2Ref();
// Set current GPD type
setGPDType(type);
// Get thresholds. Set to zero whatever is below one.
std::vector<double> thresholds;
std::vector<ActiveFlavorsThresholds> afts = m_pActiveFlavorsModule->getNfIntervals();
if (afts.size() == 1)
for (int i = 0; i < afts[0].getNf(); i++)
thresholds.push_back(0);
else
for (ActiveFlavorsThresholds aft : afts)
thresholds.push_back(aft.getLowerBound() < 1 ? 0 : sqrt(aft.getLowerBound()));
// Initialize QCD evolution objects
const auto GpdObj = apfel::InitializeGpdObjects(*m_g, thresholds, m_xi);
// Construct the DGLAP evolution operators
const auto EvolvedGpds = apfel::BuildDglap(GpdObj, initialScaleDistributions(pGPDModule), sqrt(m_MuF2_ref), this->getPertOrder() - 1, m_as);
// Set initial scale
m_MuF2_ref = pGPDModule->getMuF2Ref();
// Tabulate evolution Operator
m_tabulatedGpds = std::unique_ptr<apfel::TabulateObject<apfel::Set<apfel::Distribution>>>(new apfel::TabulateObject<apfel::Set<apfel::Distribution>>
{*EvolvedGpds, m_tabNodes, m_tabLowerBound, m_tabUpperBound, m_tabInterDegree});
//if (m_xi != m_xi_prev)
{
// Get thresholds. Set to zero whatever is below one.
std::vector<double> thresholds;
std::vector<ActiveFlavorsThresholds> afts = m_pActiveFlavorsModule->getNfIntervals();
if (afts.size() == 1)
for (int i = 0; i < afts[0].getNf(); i++)
thresholds.push_back(0);
else
for (ActiveFlavorsThresholds aft : afts)
thresholds.push_back(aft.getLowerBound() < 1 ? 0 : sqrt(aft.getLowerBound()));
// Initialize QCD evolution objects
const auto GpdObj = apfel::InitializeGpdObjects(*m_g, thresholds, m_xi);
// Construct the DGLAP evolution operators
const auto EvolvedGpds = apfel::BuildDglap(GpdObj, initialScaleDistributions(pGPDModule), sqrt(m_MuF2_ref), this->getPertOrder() - 1, m_as);
// Tabulate evolution Operator
m_tabulatedGpds = std::shared_ptr<apfel::TabulateObject<apfel::Set<apfel::Distribution>>>(new apfel::TabulateObject<apfel::Set<apfel::Distribution>>
{*EvolvedGpds, m_tabNodes, m_tabLowerBound, m_tabUpperBound, m_tabInterDegree});
m_xi_prev = m_xi;
}
// Get kinematics
const GPDKinematic kin = this->getKinematics();
......@@ -143,8 +146,8 @@ PartonDistribution GPDEvolutionApfel::compute(GPDModule* pGPDModule, const GPDTy
else if (i == 2)
j = 1;
const double q = gpds.at(i).Evaluate(this->m_x);
const double qb = gpds.at(-i).Evaluate(this->m_x);
const double q = gpds.at(i).Evaluate(this->m_x) / this->m_x;
const double qb = gpds.at(-i).Evaluate(this->m_x) / this->m_x;
qd.insert({(QuarkFlavor::Type) j, QuarkDistribution{(QuarkFlavor::Type) j, q, q + qb, q - qb}});
}
......@@ -157,7 +160,6 @@ PartonDistribution GPDEvolutionApfel::compute(GPDModule* pGPDModule, const GPDTy
std::function<std::map<int, double>(double const&, double const&)> GPDEvolutionApfel::initialScaleDistributions(GPDModule* pGPDModule) {
return [=] (double const& x, double const& muF0) -> std::map<int, double> {
// Define kinematics
const double muF02 = pow(muF0, 2);
const GPDKinematic kin{x, m_xi, m_t, muF02, muF02};
......@@ -176,8 +178,8 @@ std::function<std::map<int, double>(double const&, double const&)> GPDEvolutionA
else if ((int) t == 2)
fl = 1;
const double q = pd.getQuarkDistribution(t).getQuarkDistribution();
const double qb = pd.getQuarkDistribution(t).getQuarkDistributionPlus() - q;
const double q = x * pd.getQuarkDistribution(t).getQuarkDistribution();
const double qb = x * pd.getQuarkDistribution(t).getQuarkDistributionPlus() - q;
physd.insert({fl, q});
physd.insert({-fl, qb});
}
......@@ -252,6 +254,18 @@ double GPDEvolutionApfel::getPreviousXi() const {
return m_xi_prev;
}
std::shared_ptr<apfel::Grid> GPDEvolutionApfel::getGrid() const {
return m_g;
}
std::function<double(double const&)> GPDEvolutionApfel::getAlphas() const {
return m_as;
}
std::shared_ptr<apfel::TabulateObject<apfel::Set<apfel::Distribution>>> GPDEvolutionApfel::getTabulatedGPDs() const {
return m_tabulatedGpds;
}
GPDEvolutionApfel::GPDEvolutionApfel(const GPDEvolutionApfel &other) :
GPDEvolutionModule(other) {
m_subgridNodes = other.getSubgridNodes();
......@@ -262,18 +276,22 @@ GPDEvolutionApfel::GPDEvolutionApfel(const GPDEvolutionApfel &other) :
m_tabUpperBound = other.getTabUpperBound();
m_tabInterDegree = other.getTabInterDegree();
m_xi_prev = other.getPreviousXi();
m_g = other.getGrid();
m_as = other.getAlphas();
m_tabulatedGpds = other.getTabulatedGPDs();
}
void GPDEvolutionApfel::initModule() {
GPDEvolutionModule::initModule();
// Silence APFEL
apfel::SetVerbosityLevel(0);
// Setup APFEL++ x-space
std::vector<apfel::SubGrid> vsg;
for (int i = 0; i < (int) m_subgridNodes.size(); i++)
vsg.push_back(apfel::SubGrid{m_subgridNodes[i], m_subgridLowerBounds[i], m_subgridInterDegrees[i]});
m_g = std::unique_ptr<apfel::Grid> (new apfel::Grid{vsg});
m_g = std::shared_ptr<apfel::Grid> (new apfel::Grid{vsg});
// Running coupling
m_as = [=] (double const& mu) -> double{ return getRunningAlphaStrongModule()->compute(mu * mu); };
......
src/partons/modules/evolution/gpd/GPDEvolutionVinnikov.cpp
View file @ 2fd33ad4
......@@ -34,7 +34,7 @@ const double GPDEvolutionVinnikov::C_A = (2.0 * C_F + 1.0 / N_C);
const double GPDEvolutionVinnikov::T_R = (N_F / 2.0);
GPDEvolutionVinnikov::GPDEvolutionVinnikov(const std::string& className) :
GPDEvolutionModule(className), MathIntegratorModule(), m_gridSize(10), m_xiLast(
GPDEvolutionModule(className), MathIntegratorModule(), m_gridSize(200), m_xiLast(
0), m_tLast(0), m_MuR2Last(0), m_nFlavors(0), m_nFlavors_ref(0), m_pGPDModule(
0) {
......
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