Implement mathod_flag=0

inc/adani/AsymptoticCoefficientFunctions.h

@@ -24,7 +24,6 @@
 class AsymptoticCoefficientFunction : public AbstractHighEnergyCoefficientFunction {
     public:
         AsymptoticCoefficientFunction(const int& order, const char& kind, const char& channel, const bool& NLL = true);
-        // AsymptoticCoefficientFunction() : AsymptoticCoefficientFunction(1, '2', 'g', true) {} ;
         ~AsymptoticCoefficientFunction() ;
 
         double fx(const double x, const double m2Q2, const double m2mu2, const int nf) const ;

inc/adani/CoefficientFunction.h

@@ -18,7 +18,8 @@ class CoefficientFunction {
 
     public:
         CoefficientFunction(const int& order, const char& kind, const char& channel) ;
-        // CoefficientFunction() : CoefficientFunction(1, '2', 'g') {} ;
+        CoefficientFunction(CoefficientFunction* coeff) : CoefficientFunction(coeff -> GetOrder(), coeff -> GetKind(), coeff -> GetChannel()) {};
+
         virtual ~CoefficientFunction() = 0 ;
 
         virtual double fx(const double x, const double m2Q2, const double m2mu2, const int nf) const = 0 ;

inc/adani/Convolutions.h

@@ -31,12 +31,12 @@
 #include <gsl/gsl_monte.h>
 #include <gsl/gsl_monte_vegas.h>
 
-struct function_params {
-    double x;
-    double m2Q2;
-    int nf;
-    const AbstractConvolution* conv;
-};
+// struct function_params {
+//     double x;
+//     double m2Q2;
+//     int nf;
+//     const AbstractConvolution* conv;
+// };
 
 class AbstractConvolution {
     public:
@@ -89,31 +89,49 @@ class Convolution : public AbstractConvolution {
 
 class MonteCarloDoubleConvolution : public AbstractConvolution {
     public:
-        MonteCarloDoubleConvolution(CoefficientFunction* coefffunc, AbstractSplittingFunction* splitfunc, const double& abserr = 1e-3, const double& relerr = 1e-3, const int& dim = 1000, const int& MCcalls = 25000) ;
+        MonteCarloDoubleConvolution(CoefficientFunction* coefffunc, AbstractSplittingFunction* splitfunc, const double& abserr = 1e-3, const double& relerr = 1e-3, const int& dim = 1000, const int& method_flag = 1, const int& MCcalls = 25000) ;
         ~MonteCarloDoubleConvolution() ;
 
         double RegularPart(double x, double m2Q2, int nf) const override;
         double SingularPart(double x, double m2Q2, int nf) const override;
         double LocalPart(double x, double m2Q2, int nf) const override;
 
         // get methods
+        int GetMethodFlag() const {return method_flag_;};
         int GetMCcalls() const {return MCcalls_;};
 
         // set methods
+        void SetMethodFlag(const int& method_flag) ;
         void SetMCcalls(const int& MCcalls) ;
 
-    private:
-        int MCcalls_;
-
-        Convolution* convolution_;
-
         static double regular1_integrand(double z[], size_t dim, void *p) ;
         static double regular2_integrand(double z[], size_t dim, void *p) ;
         static double regular3_integrand(double z, void *p) ;
 
         static double singular1_integrand(double z[], size_t dim, void *p) ;
         static double singular2_integrand(double z[], size_t dim, void *p) ;
         static double singular3_integrand(double z, void *p) ;
+
+    private:
+        int method_flag_;
+        int MCcalls_;
+
+        Convolution* convolution_;
+        ConvolutedCoefficientFunction* conv_coeff_;
+
+};
+
+class ConvolutedCoefficientFunction : public CoefficientFunction {
+    public:
+        ConvolutedCoefficientFunction(Convolution* conv) : CoefficientFunction(conv -> GetCoeffFunc()) {conv_ = conv;};
+        ~ConvolutedCoefficientFunction() {} ;
+
+        double MuIndependentTerms(const double x, const double m2Q2, const int nf) const override {return conv_ -> Convolute(x, m2Q2, nf);};
+        // double MuDependentTerms(const double x, const double m2Q2, const double m2mu2, const int nf) const {return 0.;};
+        double fx(const double x, const double m2Q2, const double m2mu2, const int nf) const {return MuIndependentTerms(x, m2Q2, nf);};
+    private:
+        Convolution* conv_;
+
 };
 
 //==========================================================================================//

inc/adani/HighEnergyCoefficientFunctions.h

@@ -32,7 +32,6 @@
 class AbstractHighEnergyCoefficientFunction : public CoefficientFunction {
     public:
         AbstractHighEnergyCoefficientFunction(const int& order, const char& kind, const char& channel, const bool& NLL = true) ;
-        // AbstractHighEnergyCoefficientFunction() : AbstractHighEnergyCoefficientFunction(1, '2', 'g', true) {} ;
         ~AbstractHighEnergyCoefficientFunction() {};
 
         // get methods

inc/adani/HighScaleCoefficientFunctions.h

@@ -32,7 +32,6 @@
 class HighScaleCoefficientFunction : public CoefficientFunction {
     public:
         HighScaleCoefficientFunction(const int& order, const char& kind, const char& channel) ;
-        // HighScaleCoefficientFunction() : CoefficientFunction() {} ;
         ~HighScaleCoefficientFunction();
 
         double fx(const double x, const double m2Q2, const double m2mu2, const int nf) const ;

inc/adani/MasslessCoefficientFunctions.h

@@ -25,7 +25,6 @@ class MasslessCoefficientFunction : public CoefficientFunction {
 
     public:
         MasslessCoefficientFunction(const int& order, const char& kind, const char& channel) : CoefficientFunction(order, kind, channel) {} ;
-        // MasslessCoefficientFunction() : CoefficientFunction() {} ;
         ~MasslessCoefficientFunction() {} ;
 
         double fx(const double x, const double m2Q2, const double m2mu2, const int nf) const ;

inc/adani/SplittingFunctions.h

@@ -112,6 +112,7 @@ class SplittingFunction : public AbstractSplittingFunction{
 class ConvolutedSplittingFunctions : public AbstractSplittingFunction {
     public:
         ConvolutedSplittingFunctions(const int& order, const char& entry1, const char& entry2, const char& entry3, const char& entry4);
+        ~ConvolutedSplittingFunctions() {};
 
         char GetEntry3() const {return entry3_;} ;
 

src/Convolutions.cc

@@ -162,15 +162,33 @@ double Convolution::LocalPart(double x, double m2Q2, int nf) const {
 
 }
 
-MonteCarloDoubleConvolution::MonteCarloDoubleConvolution(CoefficientFunction* coefffunc, AbstractSplittingFunction* splitfunc, const double& abserr, const double& relerr, const int& dim, const int& MCcalls) : AbstractConvolution(coefffunc, splitfunc, abserr, relerr, dim) {
+MonteCarloDoubleConvolution::MonteCarloDoubleConvolution(CoefficientFunction* coefffunc, AbstractSplittingFunction* splitfunc, const double& abserr, const double& relerr, const int& dim, const int& method_flag, const int& MCcalls) : AbstractConvolution(coefffunc, splitfunc, abserr, relerr, dim) {
 
+    SetMethodFlag(method_flag);
     SetMCcalls(MCcalls);
-    convolution_ = new Convolution(coefffunc, splitfunc, abserr, relerr, dim);
+
+    if (method_flag == 1) {
+        convolution_ = new Convolution(coefffunc, splitfunc, abserr, relerr, dim);
+        conv_coeff_ = nullptr;
+    } else {
+        conv_coeff_ = new ConvolutedCoefficientFunction( new Convolution(coefffunc, splitfunc, abserr, relerr, dim) );
+        convolution_ = new Convolution(conv_coeff_, splitfunc, abserr, relerr, dim);
+    }
 }
 
 MonteCarloDoubleConvolution::~MonteCarloDoubleConvolution() {
 
     delete convolution_;
+    delete conv_coeff_;
+}
+
+void MonteCarloDoubleConvolution::SetMethodFlag(const int& method_flag) {
+    // check dim
+    if (method_flag != 0 && method_flag != 1) {
+        cout << "Error: method_flag must be 0 or 1. Got " << method_flag << endl;
+        exit(-1);
+    }
+    method_flag_ = method_flag;
 }
 
 void MonteCarloDoubleConvolution::SetMCcalls(const int& MCcalls) {
@@ -242,6 +260,8 @@ double regular3_integrand(double z, void *p) {
 
 double MonteCarloDoubleConvolution::RegularPart(double x, double m2Q2, int nf) const {
 
+    if (method_flag_ == 0) return convolution_ -> RegularPart(x, m2Q2, nf);
+
     double x_max = 1. / (1. + 4 * m2Q2);
     struct function_params params = { x, m2Q2, nf, this };
 
@@ -377,6 +397,9 @@ double singular3_integrand(double z, void *p) {
 
 double MonteCarloDoubleConvolution::SingularPart(double x, double m2Q2, int nf) const {
 
+    if (method_flag_ == 0) return convolution_ -> SingularPart(x, m2Q2, nf);
+
+
     double x_max = 1. / (1. + 4 * m2Q2);
     struct function_params params = { x, m2Q2, nf, this };
 
@@ -440,6 +463,8 @@ double MonteCarloDoubleConvolution::SingularPart(double x, double m2Q2, int nf)
 
 double MonteCarloDoubleConvolution::LocalPart(double x, double m2Q2, int nf) const {
 
+    if (method_flag_ == 0) return convolution_ -> LocalPart(x, m2Q2, nf);
+
     double x_max = 1. / (1. + 4 * m2Q2);
 
     return convolution_ -> Convolute(x, m2Q2, nf) * (splitfunc_ -> Local(nf) - splitfunc_ -> SingularIntegrated(x / x_max, nf));

src/ExactCoefficientFunctions.cc

@@ -385,10 +385,10 @@ double ExactCoefficientFunction::C2_g20(const double x, const double m2Q2) const
 
 double ExactCoefficientFunction::C_g21(const double x, const double m2Q2) const {
 
-    int nf = 1;
+    int nf_one = 1;
     // Put nf to 1 since the nf contribution cancels for any value of nf
 
-    return -(convolutions_lmu1[0] -> Convolute(x, m2Q2, nf) + convolutions_lmu1[1] -> Convolute(x, m2Q2, nf) * beta(0, nf));
+    return -(convolutions_lmu1[0] -> Convolute(x, m2Q2, nf_one) + convolutions_lmu1[1] -> Convolute(x, m2Q2, nf_one) * beta(0, nf_one));
 }
 
 //==========================================================================================//
@@ -579,3 +579,18 @@ double ExactCoefficientFunction::C_g32(const double x, const double m2Q2, const
 //            - 3. / 2 * beta0 * CL_g1_x_Pgg0(x, m2Q2, nf)
 //            + beta0 * beta0 * CL_g1(x, m2Q2);
 // }
+
+
+double ExactCoefficientFunction::C_g3_MuDep(const double x, const double m2Q2, const double m2mu2, int nf) const {
+
+    double lmu = log(1. / m2mu2) ;
+
+    return C_g31(x, m2Q2, nf) * lmu + C_g32(x, m2Q2, nf) * lmu * lmu ;
+}
+
+double ExactCoefficientFunction::C_ps3_MuDep(const double x, const double m2Q2, const double m2mu2, int nf) const {
+
+    double lmu = log(1. / m2mu2) ;
+
+    return C_ps31(x, m2Q2, nf) * lmu + C_ps32(x, m2Q2, nf) * lmu * lmu ;
+}