From eb1b4e97aab1cba4ee326ed6d10ab5ec1dcd1d9c Mon Sep 17 00:00:00 2001 From: Valerio Bertone Date: Tue, 3 Dec 2019 11:51:56 +0100 Subject: [PATCH] Starting reorganisation --- inc/NangaParbat/tmdgrid.h | 39 ++++ inc/NangaParbat/twobodyphasespace.h | 78 -------- src/CMakeLists.txt | 7 +- src/fastinterface/CMakeLists.txt | 2 - src/fastinterface/fastinterface.cc | 4 +- src/fastinterface/twobodyphasespace.cc | 166 ------------------ src/tmdgrid/CMakeLists.txt | 6 + .../createtmdgrid.cc | 4 +- src/tmdgrid/tmdgrid.cc | 62 +++++++ tests/IntegrationTest.cc | 3 +- tests/qTintegration.cc | 3 +- 11 files changed, 118 insertions(+), 256 deletions(-) create mode 100644 inc/NangaParbat/tmdgrid.h delete mode 100644 inc/NangaParbat/twobodyphasespace.h delete mode 100644 src/fastinterface/twobodyphasespace.cc create mode 100644 src/tmdgrid/CMakeLists.txt rename src/{fastinterface => tmdgrid}/createtmdgrid.cc (99%) create mode 100644 src/tmdgrid/tmdgrid.cc diff --git a/inc/NangaParbat/tmdgrid.h b/inc/NangaParbat/tmdgrid.h new file mode 100644 index 00000000..a37bed30 --- /dev/null +++ b/inc/NangaParbat/tmdgrid.h @@ -0,0 +1,39 @@ +// +// Author: Valerio Bertone: valerio.bertone@cern.ch +// + +#pragma once + +#include +#include + +namespace NangaParbat +{ + /** + * @brief Class for the interpolation of a single TMD grid + */ + class TMDGrid + { + public: + /** + * @brief The "TMDGrid" constructor + * @param grid: the YAML node containing the grid + */ + TMDGrid(YAML::Node const& grid); + + /** + * @brief Function that returns the value of one of the functions. + * @param x: momentum fraction + * @param qT: transverse momentum + * @param Q: renormalisation scale (assumed to be equal to the square root of zeta) + * @return It returns the value of all the flavour TMD distributions + */ + std::map Evaluate(double const& x, double const& qT, double const& Q) const; + + private: + std::unique_ptr> _xg; + std::unique_ptr> _qToQg; + std::unique_ptr> _Qg; + std::map>>> const _tmds; + }; +} diff --git a/inc/NangaParbat/twobodyphasespace.h b/inc/NangaParbat/twobodyphasespace.h deleted file mode 100644 index b631316f..00000000 --- a/inc/NangaParbat/twobodyphasespace.h +++ /dev/null @@ -1,78 +0,0 @@ -// -// Authors: Valerio Bertone: valerio.bertone@cern.ch -// Fulvio Piacenza: fulvio.piacenza01@universitadipavia.it -// - -#pragma once - -namespace NangaParbat -{ - /** - * @brief Class for the calculation of the phase-space reduction - * factor due to cuts on the single outgoing lepton in Drell-Yan - * production. The relevant process is:
- * γ(q) → l+(k1) + l-(k2)
- * with:
- * kT,1(2) > pT,min
- * |η1(2)| < ηmax
- */ - class TwoBodyPhaseSpace - { - public: - /** - * @brief The "TwoBodyPhaseSpace" constructor. - * @param pTmin: the minimum cut in the pT of the single lepton - * @param etamin: the minimum cut in the η of the single lepton - * @param etamax: the maximum cut in the η of the single lepton - * @param eps: the integration accuracy (default = 10-9) - */ - TwoBodyPhaseSpace(double const& pTmin, double const& etamin, double const& etamax, double const& eps = 1e-9); - - /** - * @brief Function that returns the phase-space reduction factor. - * @param M: invariant mass of the leptonic pair - * @param y: rapidity of the leptonic pair - * @param qT: transverse momentum of the leptonic pair - * @return the phase-space reduction factor as a function of the - * invariant mass, transverse momentum and rapidity of the lepton - * pair. - */ - double PhaseSpaceReduction(double const& Q, double const& y, double const& qT); - - /** - * @brief Function that returns the derivative w.r.t. qT of the - * phase-space reduction factor. - * @param M: invariant mass of the leptonic pair - * @param y: rapidity of the leptonic pair - * @param qT: transverse momentum of the leptonic pair - * @return the derivative of the phase-space reduction factor as a - * function of the invariant mass, transverse momentum and rapidity - * of the lepton pair. - */ - double DerivePhaseSpaceReduction(double const& Q, double const& y, double const& qT); - - /** - * @brief Function that returns the phase-space reduction factor - * associated to the parity violating contribution. - * @param M: invariant mass of the leptonic pair - * @param y: rapidity of the leptonic pair - * @param qT: transverse momentum of the leptonic pair - * @return the phase-space reduction factor as a function of the - * invariant mass, transverse momentum and rapidity of the lepton - * pair. - */ - double ParityViolatingPhaseSpaceReduction(double const& Q, double const& y, double const& qT); - - private: - /** - * @name Cut variables - * Cut Variables that define the cuts on the single leptons - */ - ///@{ - double _pTmin; - double _etamin; - double _etamax; - double _eps; - ///@} - }; -} diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 4bf62e4f..8c379967 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -4,16 +4,19 @@ add_subdirectory(chi2) add_subdirectory(minimisation) add_subdirectory(parameterisation) add_subdirectory(preprocessing) +add_subdirectory(tmdgrid) option(SHARED "Build shared-libray instead of static-libray" ON) if(SHARED) add_library(NangaParbat SHARED $ $ $ - $ $ $) + $ $ + $ $) else(SHARED) add_library(NangaParbat STATIC $ $ $ - $ $ $) + $ $ + $ $) endif(SHARED) target_link_libraries(NangaParbat ${YAML_LDFLAGS} ${APFELXX_LIBRARIES} ${ROOT_LIBRARIES} diff --git a/src/fastinterface/CMakeLists.txt b/src/fastinterface/CMakeLists.txt index 1e446e9f..37eadbd6 100644 --- a/src/fastinterface/CMakeLists.txt +++ b/src/fastinterface/CMakeLists.txt @@ -3,8 +3,6 @@ set(fastinterface_source convolutiontable.cc utilities.cc bstar.cc - twobodyphasespace.cc - createtmdgrid.cc ) add_library(fastinterface OBJECT ${fastinterface_source}) diff --git a/src/fastinterface/fastinterface.cc b/src/fastinterface/fastinterface.cc index 779eeb62..4746ba02 100644 --- a/src/fastinterface/fastinterface.cc +++ b/src/fastinterface/fastinterface.cc @@ -268,7 +268,7 @@ namespace NangaParbat const std::pair etaRange = kin.etaRange; // Allowed range in eta of the final-state leptons // Initialise two-particle phase-space object - TwoBodyPhaseSpace ps{pTMin, etaRange.first, etaRange.second}; + apfel::TwoBodyPhaseSpace ps{pTMin, etaRange.first, etaRange.second}; // Ogata-quadrature object of degree one or zero according to // weather the cross sections have to be integrated over the @@ -553,7 +553,7 @@ namespace NangaParbat const std::pair etaRange = kin.etaRange; // Allowed range in eta of the final-state leptons // Initialise two-particle phase-space object - TwoBodyPhaseSpace ps{pTMin, etaRange.first, etaRange.second}; + apfel::TwoBodyPhaseSpace ps{pTMin, etaRange.first, etaRange.second}; // Q integration bounds const double Qmin = kin.var1b.first; diff --git a/src/fastinterface/twobodyphasespace.cc b/src/fastinterface/twobodyphasespace.cc deleted file mode 100644 index 1f1eef91..00000000 --- a/src/fastinterface/twobodyphasespace.cc +++ /dev/null @@ -1,166 +0,0 @@ -// -// Authors: Valerio Bertone: valerio.bertone@cern.ch -// Fulvio Piacenza: fulvio.piacenza01@universitadipavia.it -// - -#include "NangaParbat/twobodyphasespace.h" - -#include -#include -#include - -namespace NangaParbat -{ - //_________________________________________________________________________ - TwoBodyPhaseSpace::TwoBodyPhaseSpace(double const& pTmin, double const& etamin, double const& etamax, double const& eps): - _pTmin(pTmin), - _etamin(etamin), - _etamax(etamax), - _eps(eps) - { - } - - //_________________________________________________________________________ - double TwoBodyPhaseSpace::PhaseSpaceReduction(double const& Q, double const& y, double const& qT) - { - // Return automatically zero if "y" is larger that "_etamax" - if (y <= _etamin || y >= _etamax) - return 0; - - // Useful definitions - const double Q2 = Q * Q; - const double qT2 = qT * qT; - const double qT4 = qT2 * qT2; - const double pTmin2 = _pTmin * _pTmin; - const double M2 = Q2 + qT2; - const double M = sqrt(M2); - const double ctghmax = 1 / tanh(y - _etamax); - const double ctghmin = 1 / tanh(y - _etamin); - - // Integrand function - const auto IntegrandP = [&] (double const& eta) -> double - { - // More useful definitions - const double ch = cosh(eta - y); - const double sh2 = ch * ch - 1; - const double Eq = M * ch; - const double Eq2 = Eq * Eq; - const double Eq4 = Eq2 * Eq2; - const double EmqT2 = Eq2 - qT2; - const double EmqT4 = EmqT2 * EmqT2; - const double EmqT = sqrt(EmqT2); - - // Auxiliary functions - const double f2 = ( 2 * _pTmin * Eq - Q2 ) / 2 / _pTmin / qT; - const double f3max = Eq / qT - Q2 * ( sinh(eta - y) * ctghmax + ch ) / 2 / qT / M; - const double f3min = Eq / qT - Q2 * ( sinh(eta - y) * ctghmin + ch ) / 2 / qT / M; - const double f4 = ( Eq * ( Q2 - 2 * pTmin2 + 2 * qT2 ) - Q2 * sqrt( Eq2 - M2 + pTmin2 ) ) / 2 / qT / ( M2 - pTmin2 ); - - // Integration limits in cos(phi) - const double x1 = std::max(f2, -1.); - const double x2 = std::min(std::min(f4, std::min(f3max, f3min)), 1.); - - // Return zero if x2 < x1 - if (x2 <= x1) - return 0; - - // Primitive function of the intregration in cos(phi) (up to a - // factor that can be compute externally) when contractin the - // leptonic tensor with g^{\mu\nu}_\perp. - const auto Fbar = [&] (double const& x) -> double - { - const double x2 = x * x; - const double xp = sqrt( 1 - x2 ); - const double atanfact = Eq * ( atan( ( qT - x * Eq ) / EmqT / xp ) - atan( ( qT + x * Eq ) / EmqT / xp ) ) / EmqT; - const double Fi = qT2 * x * xp / ( x2 * qT2 - Eq2 ) - atanfact; - const double Gi = Q2 * sh2 - * ( qT2 * xp * ( ( ( 11 * Eq2 * qT2 + 4 * qT4 ) * x2 + 3 * Eq * qT * ( 9 * Eq2 + qT2 ) * x + 18 * Eq4 - 5 * Eq2 * qT2 + 2 * qT4 ) / pow(x * qT + Eq, 3) + - ( ( 11 * Eq2 * qT2 + 4 * qT4 ) * x2 - 3 * Eq * qT * ( 9 * Eq2 + qT2 ) * x + 18 * Eq4 - 5 * Eq2 * qT2 + 2 * qT4 ) / pow(x * qT - Eq, 3) ) - - 6 * ( 2 * Eq2 + 3 * qT2 ) * atanfact ) / EmqT4 / 4; - - return 3 * Fi + Gi; - }; - return ( Fbar(x2) - Fbar(x1) ) / EmqT2; - }; - - // Return integral (symmetrise with respect to the center of the - // rapidity range to make sure that the result is symmetric). This - // "misbehaviour" is due to the fact that the integrand function - // is pieceswise and thus the numerical integrator struggles. - const apfel::Integrator Ieta{IntegrandP}; - return Q2 * ( y > ( _etamin + _etamax ) / 2 ? Ieta.integrate(_etamin, _etamax, _eps) : - Ieta.integrate(_etamax, _etamin, _eps) ) / 16 / M_PI; - } - - //_________________________________________________________________________ - double TwoBodyPhaseSpace::DerivePhaseSpaceReduction(double const& Q, double const& y, double const& qT) - { - const double eps = 1e-3; - return ( PhaseSpaceReduction(Q, y, qT * ( 1 + eps )) - PhaseSpaceReduction(Q, y, qT * ( 1 - eps )) ) / 2 / eps / qT; - } - - //_________________________________________________________________________ - double TwoBodyPhaseSpace::ParityViolatingPhaseSpaceReduction(double const& Q, double const& y, double const& qT) - { - // Return automatically zero if "y" is larger that "_etamax" - if (y <= _etamin || y >= _etamax) - return 0; - - // Useful definitions - const double Q2 = Q * Q; - const double qT2 = qT * qT; - const double pTmin2 = _pTmin * _pTmin; - const double M2 = Q2 + qT2; - const double M = sqrt(M2); - const double ctghmax = 1 / tanh(y - _etamax); - const double ctghmin = 1 / tanh(y - _etamin); - - // Integrand function - const auto IntegrandP = [&] (double const& eta) -> double - { - // More useful definitions - const double ch = cosh(eta - y); - const double sh = sinh(y - eta); - const double Eq = M * ch; - const double Eq2 = Eq * Eq; - const double EmqT2 = Eq2 - qT2; - const double EmqT = sqrt(EmqT2); - - // Auxiliary functions - const double f2 = ( 2 * _pTmin * Eq - Q2 ) / 2 / _pTmin / qT; - const double f3max = Eq / qT - Q2 * ( sinh(eta - y) * ctghmax + ch ) / 2 / qT / M; - const double f3min = Eq / qT - Q2 * ( sinh(eta - y) * ctghmin + ch ) / 2 / qT / M; - const double f4 = ( Eq * ( Q2 - 2 * pTmin2 + 2 * qT2 ) - Q2 * sqrt( Eq2 - M2 + pTmin2 ) ) / 2 / qT / ( M2 - pTmin2 ); - - // Integration limits in cos(phi) - const double x1 = std::max(f2, -1.); - const double x2 = std::min(std::min(f4, std::min(f3max, f3min)), 1.); - - // Return zero if x2 < x1 - if (x2 <= x1) - return 0; - - // Primitive function of the intregration in cos(phi) (up to a - // factor that can be compute externally) when contractin the - // leptonic tensor with g^{\mu\nu}_\perp. - const auto Hbar = [&] (double const& x) -> double - { - const double x2 = x * x; - const double xp = sqrt( 1 - x2 ); - const double Hi = sh - * ( qT2 * xp * ( ( 3 * Eq * qT * x - ( 4 * Eq2 - qT2 ) ) / pow(x * qT - Eq, 2) + - ( 3 * Eq * qT * x + ( 4 * Eq2 - qT2 ) ) / pow(x * qT + Eq, 2) ) - - ( 2 * Eq2 + qT2 ) - * ( atan( ( qT - x * Eq ) / EmqT / xp ) - atan( ( qT + x * Eq ) / EmqT / xp ) ) / EmqT ) / EmqT2 / EmqT2; - return Hi; - }; - return ( Hbar(x2) - Hbar(x1) ) / EmqT2; - }; - - // Return integral (symmetrise with respect to the center of the - // rapidity range to make sure that the result is symmetric). This - // "misbehaviour" is due to the fact that the integrand function - // is pieceswise and thus the numerical integrator struggles. - const apfel::Integrator Ieta{IntegrandP}; - return 3 * Q2 * Q * ( y > ( _etamin + _etamax ) / 2 ? Ieta.integrate(_etamin, _etamax, _eps) : - Ieta.integrate(_etamax, _etamin, _eps) ) / 128 / M_PI; - } -} diff --git a/src/tmdgrid/CMakeLists.txt b/src/tmdgrid/CMakeLists.txt new file mode 100644 index 00000000..49e85999 --- /dev/null +++ b/src/tmdgrid/CMakeLists.txt @@ -0,0 +1,6 @@ +set(tmdgrid_source + createtmdgrid.cc + tmdgrid.cc + ) + +add_library(tmdgrid OBJECT ${tmdgrid_source}) diff --git a/src/fastinterface/createtmdgrid.cc b/src/tmdgrid/createtmdgrid.cc similarity index 99% rename from src/fastinterface/createtmdgrid.cc rename to src/tmdgrid/createtmdgrid.cc index d780e945..6e62c7b6 100644 --- a/src/fastinterface/createtmdgrid.cc +++ b/src/tmdgrid/createtmdgrid.cc @@ -124,8 +124,8 @@ namespace NangaParbat std::unique_ptr out = std::unique_ptr(new YAML::Emitter); // Dump grids to emitter - out->SetFloatPrecision(8); - out->SetDoublePrecision(8); + out->SetFloatPrecision(6); + out->SetDoublePrecision(6); (*out) << YAML::BeginMap; (*out) << YAML::Key << "Qg" << YAML::Value << YAML::Flow << Qg; (*out) << YAML::Key << "xg" << YAML::Value << YAML::Flow << xg; diff --git a/src/tmdgrid/tmdgrid.cc b/src/tmdgrid/tmdgrid.cc new file mode 100644 index 00000000..235db561 --- /dev/null +++ b/src/tmdgrid/tmdgrid.cc @@ -0,0 +1,62 @@ +// +// Authors: Valerio Bertone: valerio.bertone@cern.ch +// + +#include "NangaParbat/tmdgrid.h" + +namespace NangaParbat +{ + //_________________________________________________________________________________ + TMDGrid::TMDGrid(YAML::Node const& grid): + _xg(std::unique_ptr>(new apfel::QGrid {grid["xg"].as>(), 3})), + _qToQg(std::unique_ptr>(new apfel::QGrid {grid["qToQg"].as>(), 3})), + _Qg(std::unique_ptr>(new apfel::QGrid {grid["Qg"].as>(), 3})), + _tmds(grid["TMDs"].as>>>>()) + { + } + + //_________________________________________________________________________________ + std::map TMDGrid::Evaluate(double const& x, double const& qT, double const& Q) const + { + // Get summation bounds + const std::tuple xbounds = _xg->SumBounds(x); + const std::tuple qToQbounds = _qToQg->SumBounds(qT/Q); + const std::tuple Qbounds = _Qg->SumBounds(Q); + + // Compute interpolating functions singularly to minimise the + // number of calls due to the the nesting of the interpolation. + // Q grid + const int inQ = std::get<1>(Qbounds); + const int nQ = std::get<2>(Qbounds) - inQ; + std::vector IQ(nQ); + for (int iQ = 0; iQ < nQ; iQ++) + IQ[iQ] = _Qg->Interpolant(std::get<0>(Qbounds), inQ + iQ, Q); + + // x grid + const int inx = std::get<1>(xbounds); + const int nx = std::get<2>(xbounds) - inx; + std::vector Ix(nx); + for (int ix = 0; ix < nx; ix++) + Ix[ix] = _xg->Interpolant(std::get<0>(xbounds), inx + ix, x); + + const int inqT = std::get<1>(qToQbounds); + const int nqT = std::get<2>(qToQbounds) - inqT; + std::vector IqT(nqT); + for (int iqT = 0; iqT < nqT; iqT++) + IqT[iqT] = _qToQg->Interpolant(std::get<0>(qToQbounds), inqT + iqT, qT/Q); + + // Initialise output map to zero + std::map result; + for (auto const& tmd : _tmds) + result.insert({tmd.first, 0}); + + // Do the interpolation + for (auto const& tmd : _tmds) + for (int iQ = 0; iQ < nQ; iQ++) + for (int ix = 0; ix < nx; ix++) + for (int iqT = 0; iqT < nqT; iqT++) + result[tmd.first] += IQ[iQ] * Ix[ix] * IqT[iqT] * tmd.second[inQ + iQ][inx + ix][inqT + iqT]; + + return result; + } +} diff --git a/tests/IntegrationTest.cc b/tests/IntegrationTest.cc index 408c5860..19c40a9c 100644 --- a/tests/IntegrationTest.cc +++ b/tests/IntegrationTest.cc @@ -6,7 +6,6 @@ #include "NangaParbat/fastinterface.h" #include "NangaParbat/convolutiontable.h" -#include "NangaParbat/twobodyphasespace.h" #include "NangaParbat/bstar.h" #include @@ -111,7 +110,7 @@ int main() apfel::OgataQuadrature OgataObj{}; // Phase-space reduction factor - NangaParbat::TwoBodyPhaseSpace ps{pTMin, etaRange.first, etaRange.second}; + apfel::TwoBodyPhaseSpace ps{pTMin, etaRange.first, etaRange.second}; // Loop over the qT-bin bounds const auto qTintegrand = [&] (double const& qT) -> double diff --git a/tests/qTintegration.cc b/tests/qTintegration.cc index a0df1eea..b4f1b2bb 100644 --- a/tests/qTintegration.cc +++ b/tests/qTintegration.cc @@ -7,7 +7,6 @@ #include "NangaParbat/fastinterface.h" #include "NangaParbat/convolutiontable.h" #include "NangaParbat/bstar.h" -#include "NangaParbat/twobodyphasespace.h" #include #include @@ -135,7 +134,7 @@ int main() apfel::OgataQuadrature bintegrand{}; // Phase-space reduction factor - NangaParbat::TwoBodyPhaseSpace ps{pTMin, etaRange.first, etaRange.second}; + apfel::TwoBodyPhaseSpace ps{pTMin, etaRange.first, etaRange.second}; // Define qT-distribution function using a DoubleObject const auto Lumi = [=] (double const& bs) -> apfel::DoubleObject