Merge pull request #49 from niclaurenti/PggxPgg

Implement analytical convolution Pgg x Pgg
niclaurenti · Apr 25, 2024 · 2362d9c · 2362d9c
2 parents f4f05d7 + b14bc92
commit 2362d9c
Show file tree

Hide file tree

Showing 13 changed files with 269 additions and 44 deletions.
diff --git a/examples/test.cpp b/examples/test.cpp
@@ -19,16 +19,16 @@ int main() {
     int nf = 4;
 
     ApproximateCoefficientFunction F2g(
-        3, '2', 'g', true, "abmp", 1e-3, 1e-3, 1000, false, 25000
+        3, '2', 'g', true, "abmp", 1e-3, 1e-3, 1000, "analytical", 25000
     );
     ApproximateCoefficientFunction FLg(
-        3, 'L', 'g', true, "abmp", 1e-3, 1e-3, 1000, false, 25000
+        3, 'L', 'g', true, "abmp", 1e-3, 1e-3, 1000, "analytical", 25000
     );
     ApproximateCoefficientFunction F2q(
-        3, '2', 'q', true, "exact", 1e-3, 1e-3, 1000, false, 25000
+        3, '2', 'q', true, "exact", 1e-3, 1e-3, 1000
     );
     ApproximateCoefficientFunction FLq(
-        3, 'L', 'q', true, "exact", 1e-3, 1e-3, 1000, false, 25000
+        3, 'L', 'q', true, "exact", 1e-3, 1e-3, 1000
     );
 
     for (int i = 0; i < N; i++) {

diff --git a/examples/test.py b/examples/test.py
@@ -4,10 +4,10 @@
 nf = 4
 m = 4.92
 
-F2g = ad.ApproximateCoefficientFunction(3, '2', 'g', True, "abmp", 1e-3, 1e-3, 1000, False, 25000)
-FLg = ad.ApproximateCoefficientFunction(3, 'L', 'g', True, "abmp", 1e-3, 1e-3, 1000, False, 25000)
-F2q = ad.ApproximateCoefficientFunction(3, '2', 'q', True, "exact", 1e-3, 1e-3, 1000, False, 25000)
-FLq = ad.ApproximateCoefficientFunction(3, 'L', 'q', True, "exact", 1e-3, 1e-3, 1000, False, 25000)
+F2g = ad.ApproximateCoefficientFunction(3, '2', 'g', True, "abmp", 1e-3, 1e-3, 1000, "analytical", 25000)
+FLg = ad.ApproximateCoefficientFunction(3, 'L', 'g', True, "abmp", 1e-3, 1e-3, 1000, "analytical", 25000)
+F2q = ad.ApproximateCoefficientFunction(3, '2', 'q', True, "exact", 1e-3, 1e-3, 1000)
+FLq = ad.ApproximateCoefficientFunction(3, 'L', 'q', True, "exact", 1e-3, 1e-3, 1000)
 
 for x in np.geomspace(1e-4, 1, 5):
     for Q in np.geomspace(10, 100, 5):

diff --git a/inc/adani/ApproximateCoefficientFunction.h b/inc/adani/ApproximateCoefficientFunction.h
@@ -64,7 +64,8 @@ class AbstractApproximate : public CoefficientFunction {
         AbstractApproximate(
             const int &order, const char &kind, const char &channel,
             const double &abserr = 1e-3, const double &relerr = 1e-3,
-            const int &dim = 1000, const bool &MCintegral = false,
+            const int &dim = 1000,
+            const string &double_int_method = "analytical",
             const int &MCcalls = 25000
         );
         ~AbstractApproximate();
@@ -91,7 +92,8 @@ class ApproximateCoefficientFunction : public AbstractApproximate {
             const int &order, const char &kind, const char &channel,
             const bool &NLL = true, const string &highscale_version = "klmv",
             const double &abserr = 1e-3, const double &relerr = 1e-3,
-            const int &dim = 1000, const bool &MCintegral = false,
+            const int &dim = 1000,
+            const string &double_int_method = "analytical",
             const int &MCcalls = 25000
         );
         ~ApproximateCoefficientFunction() override;
@@ -124,7 +126,8 @@ class ApproximateCoefficientFunctionKLMV : public AbstractApproximate {
             const int &order, const char &kind, const char &channel,
             const string &highscale_version = "klmv", const bool &lowxi = false,
             const double &abserr = 1e-3, const double &relerr = 1e-3,
-            const int &dim = 1000, const bool &MCintegral = false,
+            const int &dim = 1000,
+            const string &double_int_method = "analytical",
             const int &MCcalls = 25000
         );
         ~ApproximateCoefficientFunctionKLMV() override;

diff --git a/inc/adani/ExactCoefficientFunction.h b/inc/adani/ExactCoefficientFunction.h
@@ -22,8 +22,11 @@
 #include "adani/Convolution.h"
 #include "adani/SplittingFunction.h"
 
+#include <string>
 #include <vector>
 
+using std::string;
+
 //==========================================================================================//
 //  forward declaration of class ThresholdCoefficientFunction to avoid circular
 //  dependencies
@@ -41,7 +44,8 @@ class ExactCoefficientFunction : public CoefficientFunction {
         ExactCoefficientFunction(
             const int &order, const char &kind, const char &channel,
             const double &abserr = 1e-3, const double &relerr = 1e-3,
-            const int &dim = 1000, const bool &MCintegral = false,
+            const int &dim = 1000,
+            const string &double_int_method = "analytical",
             const int &MCcalls = 25000
         );
         ~ExactCoefficientFunction() override;
@@ -84,6 +88,7 @@ class ExactCoefficientFunction : public CoefficientFunction {
         SplittingFunction *Pgg1_;
         SplittingFunction *Pqg0_;
         ConvolutedSplittingFunctions *Pgq0Pqg0_;
+        ConvolutedSplittingFunctions *Pgg0Pgg0_;
 
         Delta *delta_;
 

diff --git a/inc/adani/SplittingFunction.h b/inc/adani/SplittingFunction.h
@@ -140,7 +140,7 @@ class ConvolutedSplittingFunctions : public AbstractSplittingFunction {
             const int &order1, const char &entry1, const char &entry2,
             const int &order2, const char &entry3, const char &entry4
         );
-        ~ConvolutedSplittingFunctions() override{};
+        ~ConvolutedSplittingFunctions() override;
 
         // overloading operators
         ConvolutedSplittingFunctions operator*(const double &rhs) const;
@@ -151,12 +151,10 @@ class ConvolutedSplittingFunctions : public AbstractSplittingFunction {
 
         // Components of the Convoluted Splitting Function
         double Regular(double x, int nf) const override;
-        double Singular(double /*x*/, int /*nf*/) const override { return 0.; };
-        double Local(int /*nf*/) const override { return 0.; };
+        double Singular(double x, int nf) const override;
+        double Local(int nf) const override;
         // Integral from 0 to x of the Singular part
-        double SingularIntegrated(double /*x*/, int /*nf*/) const override {
-            return 0.;
-        };
+        double SingularIntegrated(double x, int nf) const override;
 
         void SetFunctions();
 
@@ -177,6 +175,11 @@ class ConvolutedSplittingFunctions : public AbstractSplittingFunction {
         const char entry4_;
 
         double (ConvolutedSplittingFunctions::*reg_)(double, int) const;
+        double (ConvolutedSplittingFunctions::*sing_)(double, int) const;
+        double (ConvolutedSplittingFunctions::*sing_int_)(double, int) const;
+        double (ConvolutedSplittingFunctions::*loc_)(int) const;
+
+        SplittingFunction *Pgg0_;
 
         //==========================================================================================//
         //  Convolution between the splitting functions Pgg0/Pqq0
@@ -192,6 +195,22 @@ class ConvolutedSplittingFunctions : public AbstractSplittingFunction {
 
         double Pgq0_x_Pqg0(double x, int nf) const;
 
+        //==========================================================================================//
+        //  Convolution between the splitting functions Pgg0 and Pgg0
+        //------------------------------------------------------------------------------------------//
+
+        double Pgg0reg_x_Pgg0reg(double x) const;
+        double Pgg0reg_x_Pgg0sing(double x) const;
+        double Pgg0sing_x_Pgg0sing_reg(double x) const;
+        double Pgg0sing_x_Pgg0sing_sing(double x) const;
+        double Pgg0sing_x_Pgg0sing_sing_integrated(double x) const;
+        double Pgg0sing_x_Pgg0sing_loc() const;
+
+        double Pgg0_x_Pgg0_reg(double x, int nf) const;
+        double Pgg0_x_Pgg0_sing(double x, int nf) const;
+        double Pgg0_x_Pgg0_sing_integrated(double x, int nf) const;
+        double Pgg0_x_Pgg0_loc(int nf) const;
+
         //==========================================================================================//
         //  Function needed to return a zero function
         //------------------------------------------------------------------------------------------//

diff --git a/output/output_grid.cpp b/output/output_grid.cpp
@@ -53,6 +53,7 @@ int main(int argc, char **argv) {
         inputQ >> Qtmp;
     }
 
+    // clang-format off
     cout << "                                                                                      " << endl;
     cout << "                    888       888          888                                        " << endl;
     cout << "                    888   o   888          888                                        " << endl;
@@ -116,6 +117,7 @@ int main(int argc, char **argv) {
     cout << "                                  @@@@@@@@@@@@@@                   .+-  =@@@@@@@                     " << endl;
     cout << endl;
     cout << endl;
+    // clang-format on
 
     cout << "Computation of the grid for the coefficient function C" << channel
          << " for m = " << m << " GeV, nf = " << nf << " and µ/Q = " << mufrac

diff --git a/output/output_grid.py b/output/output_grid.py
@@ -31,7 +31,7 @@
     1e-3,
     1e-3,
     1000,
-    False,
+    "analytical",
     25000,
 )
 

diff --git a/pywrap/pywrap.cc b/pywrap/pywrap.cc
@@ -56,11 +56,11 @@ PYBIND11_MODULE(_core, m) {
             py::init<
                 const int &, const char &, const char &, const bool &,
                 const string &, const double &, const double &, const int &,
-                const bool &, const int &>(),
+                const string &, const int &>(),
             py::arg("order"), py::arg("kind"), py::arg("channel"),
             py::arg("NLL") = true, py::arg("highscale_version") = "klmv",
             py::arg("abserr") = 1e-3, py::arg("relerr") = 1e-3,
-            py::arg("dim") = 1000, py::arg("MCintegral") = false,
+            py::arg("dim") = 1000, py::arg("double_int_method") = "analytical",
             py::arg("MCcalls") = 25000
         )
         .def(
@@ -100,11 +100,11 @@ PYBIND11_MODULE(_core, m) {
             py::init<
                 const int &, const char &, const char &, const string &,
                 const bool &, const double &, const double &, const int &,
-                const bool &, const int &>(),
+                const string &, const int &>(),
             py::arg("order"), py::arg("kind"), py::arg("channel"),
             py::arg("highscale_version") = "klmv", py::arg("lowxi") = false,
             py::arg("abserr") = 1e-3, py::arg("relerr") = 1e-3,
-            py::arg("dim") = 1000, py::arg("MCintegral") = false,
+            py::arg("dim") = 1000, py::arg("double_int_method") = "analytical",
             py::arg("MCcalls") = 25000
         )
         .def(
@@ -183,10 +183,10 @@ PYBIND11_MODULE(_core, m) {
         .def(
             py::init<
                 const int &, const char &, const char &, const double &,
-                const double &, const int &, const bool &, const int &>(),
+                const double &, const int &, const string &, const int &>(),
             py::arg("order"), py::arg("kind"), py::arg("channel"),
             py::arg("abserr") = 1e-3, py::arg("relerr") = 1e-3,
-            py::arg("dim") = 1000, py::arg("MCintegral") = false,
+            py::arg("dim") = 1000, py::arg("double_int_method") = "analytical",
             py::arg("MCcalls") = 25000
         )
         .def(

diff --git a/src/ApproximateCoefficientFunction.cc b/src/ApproximateCoefficientFunction.cc
@@ -20,12 +20,12 @@ using std::vector;
 AbstractApproximate::AbstractApproximate(
     const int &order, const char &kind, const char &channel,
     const double &abserr, const double &relerr, const int &dim,
-    const bool &MCintegral, const int &MCcalls
+    const string &double_int_method, const int &MCcalls
 )
     : CoefficientFunction(order, kind, channel) {
 
     muterms_ = new ExactCoefficientFunction(
-        order, kind, channel, abserr, relerr, dim, MCintegral, MCcalls
+        order, kind, channel, abserr, relerr, dim, double_int_method, MCcalls
     );
 }
 
@@ -35,6 +35,10 @@ AbstractApproximate::AbstractApproximate(
 
 AbstractApproximate::~AbstractApproximate() { delete muterms_; }
 
+//==========================================================================================//
+//  AbstractApproximate: central value of the mu-independent terms
+//------------------------------------------------------------------------------------------//
+
 double AbstractApproximate::MuIndependentTerms(
     double x, double m2Q2, int nf
 ) const {
@@ -96,10 +100,10 @@ struct variation_parameters CL_var = { 0.2, 2. };
 ApproximateCoefficientFunction::ApproximateCoefficientFunction(
     const int &order, const char &kind, const char &channel, const bool &NLL,
     const string &highscale_version, const double &abserr, const double &relerr,
-    const int &dim, const bool &MCintegral, const int &MCcalls
+    const int &dim, const string &double_int_method, const int &MCcalls
 )
     : AbstractApproximate(
-          order, kind, channel, abserr, relerr, dim, MCintegral, MCcalls
+          order, kind, channel, abserr, relerr, dim, double_int_method, MCcalls
       ) {
 
     threshold_ = new ThresholdCoefficientFunction(order, kind, channel);
@@ -262,11 +266,11 @@ struct klmv_params klmv_C2g3B_lowxi = { 0.8, 10.7, 0.055125, 2, 0.3825 };
 ApproximateCoefficientFunctionKLMV::ApproximateCoefficientFunctionKLMV(
     const int &order, const char &kind, const char &channel,
     const string &highscale_version, const bool &lowxi, const double &abserr,
-    const double &relerr, const int &dim, const bool &MCintegral,
+    const double &relerr, const int &dim, const string &double_int_method,
     const int &MCcalls
 )
     : AbstractApproximate(
-          order, kind, channel, abserr, relerr, dim, MCintegral, MCcalls
+          order, kind, channel, abserr, relerr, dim, double_int_method, MCcalls
       ) {
     if (GetOrder() == 1) {
         cout << "Error: KLMV approximation is not implemented at O(as)!"

diff --git a/src/ExactCoefficientFunction.cc b/src/ExactCoefficientFunction.cc
@@ -16,7 +16,7 @@ using std::endl;
 ExactCoefficientFunction::ExactCoefficientFunction(
     const int &order, const char &kind, const char &channel,
     const double &abserr, const double &relerr, const int &dim,
-    const bool &MCintegral, const int &MCcalls
+    const string &double_int_method, const int &MCcalls
 )
     : CoefficientFunction(order, kind, channel) {
 
@@ -33,9 +33,20 @@ ExactCoefficientFunction::ExactCoefficientFunction(
     Pgg1_ = nullptr;
     Pqg0_ = nullptr;
     Pgq0Pqg0_ = nullptr;
+    Pgg0Pgg0_ = nullptr;
 
     delta_ = nullptr;
 
+    // check double_int_method
+    if (double_int_method != "analytical"
+        && double_int_method != "double_numerical"
+        && double_int_method != "monte_carlo") {
+        cout << "Error in ExactCoefficientFunction: double_int_method must be "
+                "'analytical', 'double_numerical' or 'monte_carlo'! Got "
+             << double_int_method << endl;
+        exit(-1);
+    }
+
     if (GetOrder() == 2) {
         asy_ = new AsymptoticCoefficientFunction(order, kind, channel);
         thr_ = new ThresholdCoefficientFunction(order, kind, channel);
@@ -62,6 +73,9 @@ ExactCoefficientFunction::ExactCoefficientFunction(
         Pgg1_ = new SplittingFunction(1, 'g', 'g');
         Pqg0_ = new SplittingFunction(0, 'q', 'g');
         Pgq0Pqg0_ = new ConvolutedSplittingFunctions(0, 'g', 'q', 0, 'q', 'g');
+        if (double_int_method == "analytical")
+            Pgg0Pgg0_ =
+                new ConvolutedSplittingFunctions(0, 'g', 'g', 0, 'g', 'g');
     }
 
     if (GetOrder() == 2) {
@@ -110,9 +124,19 @@ ExactCoefficientFunction::ExactCoefficientFunction(
             );
             convolutions_lmu1_.push_back(new Convolution(gluon_as2_, delta_));
 
-            convolutions_lmu2_.push_back(new DoubleConvolution(
-                gluon_as1_, Pgg0_, abserr, relerr, dim, MCintegral, MCcalls
-            ));
+            if (double_int_method == "monte_carlo") {
+                convolutions_lmu2_.push_back(new DoubleConvolution(
+                    gluon_as1_, Pgg0_, abserr, relerr, dim, true, MCcalls
+                ));
+            } else if (double_int_method == "double_numerical") {
+                convolutions_lmu2_.push_back(new DoubleConvolution(
+                    gluon_as1_, Pgg0_, abserr, relerr, dim, false, MCcalls
+                ));
+            } else {
+                convolutions_lmu2_.push_back(
+                    new Convolution(gluon_as1_, Pgg0Pgg0_, abserr, relerr, dim)
+                );
+            }
             convolutions_lmu2_.push_back(
                 new Convolution(gluon_as1_, Pgq0Pqg0_, abserr, relerr, dim)
             );
@@ -153,6 +177,7 @@ ExactCoefficientFunction::~ExactCoefficientFunction() {
     delete Pgg1_;
     delete Pqg0_;
     delete Pgq0Pqg0_;
+    delete Pgg0Pgg0_;
 
     delete delta_;