From c2f607ed086226429dc55541393ac541d50055fd Mon Sep 17 00:00:00 2001 From: niccolo Date: Tue, 12 Mar 2024 17:37:20 +0100 Subject: [PATCH] Make integrand of convolutions static functions --- inc/adani/CoefficientFunction.h | 22 +-- inc/adani/Convolutions.h | 41 ++--- inc/adani/ExactCoefficientFunctions.h | 10 +- inc/adani/SplittingFunctions.h | 53 +++++- src/Convolutions.cc | 111 ++++++------- src/ExactCoefficientFunctions.cc | 142 ++++++++-------- src/SplittingFunctions.cc | 223 ++++++++++++++++++-------- 7 files changed, 356 insertions(+), 246 deletions(-) diff --git a/inc/adani/CoefficientFunction.h b/inc/adani/CoefficientFunction.h index 18bdf5f..b2f9ab0 100644 --- a/inc/adani/CoefficientFunction.h +++ b/inc/adani/CoefficientFunction.h @@ -3,10 +3,15 @@ #include "adani/SplittingFunctions.h" -struct result { - double central; - double higher; - double lower; +class Value { + public: + Value(const double& central, const double& higher, const double& lower) ; + + + private: + double central; + double higher; + double lower; }; class CoefficientFunction { @@ -18,6 +23,9 @@ class CoefficientFunction { virtual double fx(const double x, const double m2Q2, const double m2mu2, const int nf) const = 0 ; virtual double MuIndependentTerms(const double x, const double m2Q2, const int nf) const {return fx(x, m2Q2, 1., nf);} ; + virtual double MuDependentTerms(const double x, const double m2Q2, const double m2mu2, const int nf) const { + return fx(x, m2Q2, m2mu2, nf) - fx(x, m2Q2, 1., nf); + } // get methods int GetOrder() const { return order_; } ; @@ -29,17 +37,11 @@ class CoefficientFunction { void SetKind(const char& kind) ; void SetChannel(const char& channel) ; - // virtual void Set_fx() = 0; - private: int order_ ; // order = 1, 2, or 3 char kind_ ; // kind_ = '2' for F2 and 'L' for FL char channel_ ; // channel_ = 'g' for Cg and 'q' for Cq - // TODO: IMPLEMENT POINTER TO THE RIGHT FUNCTION - // protected: - // double (*fx_)(double, double, double, int); - }; #endif diff --git a/inc/adani/Convolutions.h b/inc/adani/Convolutions.h index c596e9c..3d054eb 100644 --- a/inc/adani/Convolutions.h +++ b/inc/adani/Convolutions.h @@ -35,12 +35,12 @@ struct function_params { double x; double m2Q2; int nf; + const AbstractConvolution* conv; }; class AbstractConvolution { public: - // AbstractConvolution(const CoefficientFunction& coefffunc, const SplittingFunction& splitfunc, const double& abserr = 1e-3, const double& relerr = 1e-3, const int& dim = 1000); - AbstractConvolution(CoefficientFunction* coefffunc, SplittingFunction* splitfunc, const double& abserr = 1e-3, const double& relerr = 1e-3, const int& dim = 1000); + AbstractConvolution(CoefficientFunction* coefffunc, AbstractSplittingFunction* splitfunc, const double& abserr = 1e-3, const double& relerr = 1e-3, const int& dim = 1000); ~AbstractConvolution() {} ; double Convolute(double x, double m2Q2, int nf) const ; @@ -54,8 +54,8 @@ class AbstractConvolution { int GetRelerr() const {return relerr_;}; int GetDim() const {return dim_;} ; - // CoefficientFunction* GetCoeffFunc() const {return coefffunc_;}; - // SplittingFunction* GetSplitFunc() const {return splitfunc_;}; + CoefficientFunction* GetCoeffFunc() const {return coefffunc_;}; + AbstractSplittingFunction* GetSplitFunc() const {return splitfunc_;}; // set methods void SetAbserr(const double& abserr) ; @@ -69,38 +69,27 @@ class AbstractConvolution { protected: CoefficientFunction *coefffunc_ ; - SplittingFunction *splitfunc_ ; + AbstractSplittingFunction *splitfunc_ ; }; class Convolution : public AbstractConvolution { public: - Convolution(CoefficientFunction* coefffunc, SplittingFunction* splitfunc, const double& abserr = 1e-3, const double& relerr = 1e-3, const int& dim = 1000) : AbstractConvolution(coefffunc, splitfunc, abserr, relerr, dim) {} ; + Convolution(CoefficientFunction* coefffunc, AbstractSplittingFunction* splitfunc, const double& abserr = 1e-3, const double& relerr = 1e-3, const int& dim = 1000) : AbstractConvolution(coefffunc, splitfunc, abserr, relerr, dim) {} ; ~Convolution() {} ; - // double Convolute(double x, double m2Q2, int nf) const ; - 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; - // friend double regular_integrand(double z, void *p) ; - // friend double singular_integrand(double z, void *p) ; - - private: - double regular_integrand(double z, void *p) const ; - double singular_integrand(double z, void *p) const ; - - // static double static_regular_integrand(double z, void *params) { - // // Call the member function through the params pointer - // return static_cast(params)->regular_integrand(z, nullptr); - // } + static double regular_integrand(double z, void *p) ; + static double singular_integrand(double z, void *p) ; }; class MonteCarloDoubleConvolution : public AbstractConvolution { public: - MonteCarloDoubleConvolution(CoefficientFunction* coefffunc, SplittingFunction* 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& MCcalls = 25000) ; ~MonteCarloDoubleConvolution() ; double RegularPart(double x, double m2Q2, int nf) const override; @@ -118,13 +107,13 @@ class MonteCarloDoubleConvolution : public AbstractConvolution { Convolution* convolution_; - double regular1_integrand(double z[], size_t dim, void *p) const ; - double regular2_integrand(double z[], size_t dim, void *p) const ; - double regular3_integrand(double z, void *p) const ; + 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) ; - double singular1_integrand(double z[], size_t dim, void *p) const ; - double singular2_integrand(double z[], size_t dim, void *p) const ; - double singular3_integrand(double z, void *p) const ; + 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) ; }; //==========================================================================================// diff --git a/inc/adani/ExactCoefficientFunctions.h b/inc/adani/ExactCoefficientFunctions.h index 70cbe46..b738436 100644 --- a/inc/adani/ExactCoefficientFunctions.h +++ b/inc/adani/ExactCoefficientFunctions.h @@ -59,6 +59,8 @@ class ExactCoefficientFunction : public CoefficientFunction { double MuIndependentTerms(const double x, const double m2Q2, const int nf) const override ; double MuDependentTerms(const double x, const double m2Q2, const double m2mu2, const int nf) const ; + void SetFunctions(); + private: int method_flag_ ; @@ -67,7 +69,11 @@ class ExactCoefficientFunction : public CoefficientFunction { int MCcalls_ ; int dim_ ; - std::vector convolutions_; + double (ExactCoefficientFunction::*mu_indep_)(double, double, int) const; + double (ExactCoefficientFunction::*mu_dep_)(double, double, double, int) const; + + std::vector convolutions_lmu1; + std::vector convolutions_lmu2; ExactCoefficientFunction* gluon_lo_ ; ExactCoefficientFunction* gluon_nlo_ ; @@ -136,6 +142,8 @@ class ExactCoefficientFunction : public CoefficientFunction { double C_ps3_MuDep(const double x, const double m2Q2, const double m2mu2, const int nf) const ; double C_g3_MuDep(const double x, const double m2Q2, const double m2mu2, const int nf) const ; + double ZeroFunction(double x, double m2mu2, int nf) const {return 0.;}; + }; #endif diff --git a/inc/adani/SplittingFunctions.h b/inc/adani/SplittingFunctions.h index 0765d8f..8de5e72 100644 --- a/inc/adani/SplittingFunctions.h +++ b/inc/adani/SplittingFunctions.h @@ -47,6 +47,10 @@ class SplittingFunction : public AbstractSplittingFunction{ SplittingFunction operator*(const double& rhs) const; friend SplittingFunction operator*(const double& lhs, const SplittingFunction& rhs); + SplittingFunction operator/(const double& rhs) const; + + void SetFunctions() ; + // get methods double GetOrder() const {return order_ ;} ; char GetEntry1() const {return entry1_;} ; @@ -57,6 +61,11 @@ class SplittingFunction : public AbstractSplittingFunction{ char entry1_; char entry2_; + double (SplittingFunction::*reg_)(double, int) const; + double (SplittingFunction::*sing_)(double, int) const; + double (SplittingFunction::*sing_int_)(double, int) const; + double (SplittingFunction::*loc_)(int) const; + //==========================================================================================// // Splitting functions O(alpha_s) // without color factors @@ -95,25 +104,44 @@ class SplittingFunction : public AbstractSplittingFunction{ double Pgg1sing_integrated(const double x, const int nf) const ; double Pgg1loc(const int nf) const ; + double ZeroFunction_x_nf(double x, int nf) const {return 0.;}; + double ZeroFunction_nf(int nf) const {return 0.;}; + }; -class ConvolutedSplittingFunctions : public SplittingFunction { +class ConvolutedSplittingFunctions : public AbstractSplittingFunction { public: - ConvolutedSplittingFunctions(const int& order, const char& entry1, const char& entry2, const char& entry3); + ConvolutedSplittingFunctions(const int& order, const char& entry1, const char& entry2, const char& entry3, const char& entry4); char GetEntry3() const {return entry3_;} ; double Regular(const double x, const int nf) const ; - double Singular(const double x, const int nf) const {return 0.;} ; - double Local(const int nf) const {return 0.;} ; - double SingularIntegrated(const double x, const int nf) const {return 0.;}; + double Singular(const double x, const int nf) const override {return 0.;} ; + double Local(const int nf) const override {return 0.;} ; + double SingularIntegrated(const double x, const int nf) const override {return 0.;}; ConvolutedSplittingFunctions operator*(const double& rhs) const; friend ConvolutedSplittingFunctions operator*(const double& lhs, const ConvolutedSplittingFunctions& rhs); + ConvolutedSplittingFunctions operator/(const double& rhs) const; + + void SetFunctions() ; + + // get methods + double GetOrder() const {return order_ ;} ; + char GetEntry1() const {return entry1_;} ; + char GetEntry2() const {return entry2_;} ; + char GetEntry3() const {return entry3_;} ; + private: + int order_; + char entry1_; + char entry2_; char entry3_; + char entry4_; + + double (ConvolutedSplittingFunctions::*reg_)(double, int) const; //==========================================================================================// // Convolution between the splitting functions Pgg0/Pqq0 @@ -128,14 +156,25 @@ class ConvolutedSplittingFunctions : public SplittingFunction { //------------------------------------------------------------------------------------------// double Pgq0_x_Pqg0(const double x, const int nf) const ; + + double ZeroFunction_x_nf(double x, int nf) const {return 0.;}; + double ZeroFunction_nf(int nf) const {return 0.;}; }; -class Delta : AbstractSplittingFunction { +class Delta : public AbstractSplittingFunction { public: + Delta() : AbstractSplittingFunction() {}; + ~Delta() {} ; + double Regular(const double x, const int nf) const {return 0.;}; double Singular(const double x, const int nf) const {return 0.;}; - double Local(const int nf) const {return 1.;}; + double Local(const int nf) const {return GetMultFact() * 1.;}; double SingularIntegrated(const double x, const int nf) const {return 0.;}; + + Delta operator*(const double& rhs) const; + friend Delta operator*(const double& lhs, const Delta& rhs); + + Delta operator/(const double& rhs) const; }; #endif diff --git a/src/Convolutions.cc b/src/Convolutions.cc index 8461628..d73429f 100644 --- a/src/Convolutions.cc +++ b/src/Convolutions.cc @@ -15,16 +15,7 @@ using std::endl; // TODO : in all the numerical convolutions the gsl default error is first // switched off and then swithced on again: see if it can be done once for all -// AbstractConvolution::AbstractConvolution(const CoefficientFunction& coefffunc, const SplittingFunction& splitfunc, const double& abserr, const double& relerr, const int& dim) { -// SetAbserr(abserr); -// SetRelerr(relerr); -// SetDim(dim); - -// coefffunc_ = &coefffunc; -// splitfunc_ = &splitfunc; -// } - -AbstractConvolution::AbstractConvolution(CoefficientFunction* coefffunc, SplittingFunction* splitfunc, const double& abserr, const double& relerr, const int& dim) { +AbstractConvolution::AbstractConvolution(CoefficientFunction* coefffunc, AbstractSplittingFunction* splitfunc, const double& abserr, const double& relerr, const int& dim) { abserr_ = abserr; relerr_ = relerr; dim_ = dim; @@ -64,31 +55,7 @@ double AbstractConvolution::Convolute(double x, double m2Q2, int nf) const { return RegularPart(x, m2Q2, nf) + SingularPart(x, m2Q2, nf) + LocalPart(x, m2Q2, nf) ; } - -// double regular_integrand(double z, void *p) { - -// struct function_params *params = (struct function_params *)p; - -// double m2Q2 = (params->m2Q2); -// double x = (params->x); -// int nf = (params->nf); - -// return coefffunc_ -> MuIndependentTerms(z, m2Q2, nf) * splitfunc_ -> Regular(x / z, nf) / z ; -// } - -// double singular_integrand(double z, void *p) { - -// struct function_params *params = (struct function_params *)p; - -// double m2Q2 = (params->m2Q2); -// double x = (params->x); -// int nf = (params->nf); - -// return splitfunc_ -> Regular(z, nf) * (coefffunc_-> MuIndependentTerms(x / z, m2Q2, nf) - coefffunc_ -> MuIndependentTerms(x, m2Q2, nf) ) ; - -// } - -double Convolution::regular_integrand(double z, void *p) const { +double regular_integrand(double z, void *p) { struct function_params *params = (struct function_params *)p; @@ -96,17 +63,17 @@ double Convolution::regular_integrand(double z, void *p) const { double x = (params->x); int nf = (params->nf); - return coefffunc_ -> MuIndependentTerms(z, m2Q2, nf) * splitfunc_ -> Regular(x / z, nf) / z ; + return ((params->conv) -> GetCoeffFunc()) -> MuIndependentTerms(z, m2Q2, nf) * ((params->conv) -> GetSplitFunc()) -> Regular(x / z, nf) / z ; } -double Convolution::singular_integrand(double z, void *p) const { +double singular_integrand(double z, void *p) { struct function_params *params = (struct function_params *)p; double m2Q2 = (params->m2Q2); double x = (params->x); int nf = (params->nf); - return splitfunc_ -> Singular(z, nf) * (coefffunc_-> MuIndependentTerms(x / z, m2Q2, nf) - coefffunc_ -> MuIndependentTerms(x, m2Q2, nf) ) ; + return ((params->conv) -> GetSplitFunc()) -> Singular(z, nf) * (((params->conv) -> GetCoeffFunc())-> MuIndependentTerms(x / z, m2Q2, nf) - ((params->conv) -> GetCoeffFunc()) -> MuIndependentTerms(x, m2Q2, nf) ) ; } @@ -132,7 +99,7 @@ double Convolution::RegularPart(double x, double m2Q2, int nf) const { double regular, error; - struct function_params params = { x, m2Q2, nf }; + struct function_params params = { x, m2Q2, nf, this }; gsl_function F; F.function = ®ular_integrand; @@ -163,7 +130,7 @@ double Convolution::SingularPart(double x, double m2Q2, int nf) const { double singular, error; - struct function_params params = { x, m2Q2, nf }; + struct function_params params = { x, m2Q2, nf, this }; gsl_function F; F.function = &singular_integrand; @@ -195,7 +162,7 @@ double Convolution::LocalPart(double x, double m2Q2, int nf) const { } -MonteCarloDoubleConvolution::MonteCarloDoubleConvolution(CoefficientFunction* coefffunc, SplittingFunction* 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& MCcalls) : AbstractConvolution(coefffunc, splitfunc, abserr, relerr, dim) { SetMCcalls(MCcalls); convolution_ = new Convolution(coefffunc, splitfunc, abserr, relerr, dim); @@ -215,7 +182,7 @@ void MonteCarloDoubleConvolution::SetMCcalls(const int& MCcalls) { MCcalls_ = MCcalls; } -double MonteCarloDoubleConvolution::regular1_integrand(double z[], size_t dim, void *p) const { +double regular1_integrand(double z[], size_t dim, void *p) { struct function_params *params = (struct function_params *)p; @@ -223,34 +190,40 @@ double MonteCarloDoubleConvolution::regular1_integrand(double z[], size_t dim, v double x = (params->x); int nf = (params->nf); + CoefficientFunction* coefffunc = (params->conv) -> GetCoeffFunc(); + AbstractSplittingFunction* splitfunc = (params->conv) -> GetSplitFunc(); + double z1 = z[0], z2 = z[1]; if (z2 > z1) { - return 1. / (z1 * z2) * splitfunc_ -> Regular(x / z1, nf) * splitfunc_ -> Regular(z1 / z2, nf) - * coefffunc_-> MuIndependentTerms(z2, m2Q2, nf); + return 1. / (z1 * z2) * splitfunc -> Regular(x / z1, nf) * splitfunc -> Regular(z1 / z2, nf) + * coefffunc-> MuIndependentTerms(z2, m2Q2, nf); } else { return 0.; } } -double MonteCarloDoubleConvolution::regular2_integrand(double z[], size_t dim, void *p) const { +double regular2_integrand(double z[], size_t dim, void *p) { struct function_params *params = (struct function_params *)p; double m2Q2 = (params->m2Q2); double x = (params->x); int nf = (params->nf); + CoefficientFunction* coefffunc = (params->conv) -> GetCoeffFunc(); + AbstractSplittingFunction* splitfunc = (params->conv) -> GetSplitFunc(); + double z1 = z[0], z2 = z[1]; if (z2 > z1) { - return 1. / (z1 * z2) * splitfunc_ -> Regular(x / z1, nf) * splitfunc_ -> Singular(z1 / z2, nf) - * (coefffunc_-> MuIndependentTerms(z2, m2Q2, nf) - z1 / z2 * coefffunc_-> MuIndependentTerms(z1, m2Q2, nf)); + return 1. / (z1 * z2) * splitfunc -> Regular(x / z1, nf) * splitfunc -> Singular(z1 / z2, nf) + * (coefffunc-> MuIndependentTerms(z2, m2Q2, nf) - z1 / z2 * coefffunc-> MuIndependentTerms(z1, m2Q2, nf)); } else { return 0.; } } -double MonteCarloDoubleConvolution::regular3_integrand(double z, void *p) const { +double regular3_integrand(double z, void *p) { struct function_params *params = (struct function_params *)p; @@ -258,16 +231,19 @@ double MonteCarloDoubleConvolution::regular3_integrand(double z, void *p) const double x = (params->x); int nf = (params->nf); + CoefficientFunction* coefffunc = (params->conv) -> GetCoeffFunc(); + AbstractSplittingFunction* splitfunc = (params->conv) -> GetSplitFunc(); + double x_max = 1. / (1. + 4 * m2Q2); - return -1. / z * splitfunc_ -> Regular(x / z, nf) * coefffunc_-> MuIndependentTerms(z, m2Q2, nf) - * splitfunc_ -> SingularIntegrated(z / x_max, nf); + return -1. / z * splitfunc -> Regular(x / z, nf) * coefffunc-> MuIndependentTerms(z, m2Q2, nf) + * splitfunc -> SingularIntegrated(z / x_max, nf); } double MonteCarloDoubleConvolution::RegularPart(double x, double m2Q2, int nf) const { double x_max = 1. / (1. + 4 * m2Q2); - struct function_params params = { x, m2Q2, nf }; + struct function_params params = { x, m2Q2, nf, this }; double xl[2] = { x, x }; double xu[2] = { x_max, x_max }; @@ -330,7 +306,7 @@ double MonteCarloDoubleConvolution::RegularPart(double x, double m2Q2, int nf) c // monte carlo methods //------------------------------------------------------------------------------------------// -double MonteCarloDoubleConvolution::singular1_integrand(double z[], size_t dim, void *p) const { +double singular1_integrand(double z[], size_t dim, void *p) { struct function_params *params = (struct function_params *)p; @@ -338,19 +314,22 @@ double MonteCarloDoubleConvolution::singular1_integrand(double z[], size_t dim, double x = (params->x); int nf = (params->nf); + CoefficientFunction* coefffunc = (params->conv) -> GetCoeffFunc(); + AbstractSplittingFunction* splitfunc = (params->conv) -> GetSplitFunc(); + double z1 = z[0], z2 = z[1]; double tmp; if (z2 > x / z1) { - tmp = splitfunc_ -> Regular(x / (z1 * z2), nf) / z1; + tmp = splitfunc -> Regular(x / (z1 * z2), nf) / z1; } else { tmp = 0.; } - return splitfunc_ -> Singular(z1, nf) * (tmp - splitfunc_ -> Regular(x / z2, nf)) * coefffunc_ -> MuIndependentTerms(z2, m2Q2, nf) / z2; + return splitfunc -> Singular(z1, nf) * (tmp - splitfunc -> Regular(x / z2, nf)) * coefffunc -> MuIndependentTerms(z2, m2Q2, nf) / z2; } -double MonteCarloDoubleConvolution::singular2_integrand(double z[], size_t dim, void *p) const { +double singular2_integrand(double z[], size_t dim, void *p) { struct function_params *params = (struct function_params *)p; @@ -358,22 +337,25 @@ double MonteCarloDoubleConvolution::singular2_integrand(double z[], size_t dim, double x = (params->x); int nf = (params->nf); + CoefficientFunction* coefffunc = (params->conv) -> GetCoeffFunc(); + AbstractSplittingFunction* splitfunc = (params->conv) -> GetSplitFunc(); + double x_max = 1. / (1. + 4 * m2Q2); double z1 = z[0], z2 = z[1]; double tmp; if (z2 > x / (x_max * z1)) { - tmp = (coefffunc_ -> MuIndependentTerms(x / (z1 * z2), m2Q2, nf) / z2 - coefffunc_ -> MuIndependentTerms(x / z1, m2Q2, nf)) / z1; + tmp = (coefffunc -> MuIndependentTerms(x / (z1 * z2), m2Q2, nf) / z2 - coefffunc -> MuIndependentTerms(x / z1, m2Q2, nf)) / z1; } else { tmp = 0.; } - return splitfunc_ -> Singular(z1, nf) * splitfunc_ -> Singular(z2, nf) - * (tmp - (coefffunc_ -> MuIndependentTerms(x / z2, m2Q2, nf) / z2 - coefffunc_ -> MuIndependentTerms(x, m2Q2, nf))); + return splitfunc -> Singular(z1, nf) * splitfunc -> Singular(z2, nf) + * (tmp - (coefffunc -> MuIndependentTerms(x / z2, m2Q2, nf) / z2 - coefffunc -> MuIndependentTerms(x, m2Q2, nf))); } -double MonteCarloDoubleConvolution::singular3_integrand(double z, void *p) const { +double singular3_integrand(double z, void *p) { struct function_params *params = (struct function_params *)p; @@ -381,19 +363,22 @@ double MonteCarloDoubleConvolution::singular3_integrand(double z, void *p) const double x = (params->x); int nf = (params->nf); + CoefficientFunction* coefffunc = (params->conv) -> GetCoeffFunc(); + AbstractSplittingFunction* splitfunc = (params->conv) -> GetSplitFunc(); + double x_max = 1. / (1. + 4 * m2Q2); return -( - splitfunc_ -> Singular(z, nf) - * (coefffunc_ -> MuIndependentTerms(x / z, m2Q2, nf) * splitfunc_ -> SingularIntegrated(x / (x_max * z), nf) / z - - coefffunc_ -> MuIndependentTerms(x, m2Q2, nf) * splitfunc_ -> SingularIntegrated(x / x_max, nf)) + splitfunc -> Singular(z, nf) + * (coefffunc -> MuIndependentTerms(x / z, m2Q2, nf) * splitfunc -> SingularIntegrated(x / (x_max * z), nf) / z + - coefffunc -> MuIndependentTerms(x, m2Q2, nf) * splitfunc -> SingularIntegrated(x / x_max, nf)) ); } double MonteCarloDoubleConvolution::SingularPart(double x, double m2Q2, int nf) const { double x_max = 1. / (1. + 4 * m2Q2); - struct function_params params = { x, m2Q2, nf }; + struct function_params params = { x, m2Q2, nf, this }; double xl[2] = { x / x_max, x }; double xu[2] = { 1, x_max }; diff --git a/src/ExactCoefficientFunctions.cc b/src/ExactCoefficientFunctions.cc index 30440d8..7a61463 100644 --- a/src/ExactCoefficientFunctions.cc +++ b/src/ExactCoefficientFunctions.cc @@ -30,93 +30,85 @@ ExactCoefficientFunction::ExactCoefficientFunction(const int& order, const char& if (GetOrder() == 2) { if (GetChannel() == 'q') { - convolutions_.push_back( &Convolution(gluon_lo_, &SplittingFunction(0, 'g', 'q'), abserr, relerr, dim) ); + convolutions_lmu1.push_back( &Convolution(gluon_lo_, &SplittingFunction(0, 'g', 'q'), abserr, relerr, dim) ); } else if (GetChannel() == 'g'){ - convolutions_.push_back( new Convolution(gluon_lo_, &SplittingFunction(0, 'g', 'g'), abserr, relerr, dim) ); + convolutions_lmu1.push_back( new Convolution(gluon_lo_, &SplittingFunction(0, 'g', 'g'), abserr, relerr, dim) ); + convolutions_lmu1.push_back( new Convolution(gluon_lo_, &((-1)*Delta())) ); } } else if (GetOrder() == 3) { if (GetChannel() == 'q') { - convolutions_.push_back( new Convolution(gluon_lo_, &SplittingFunction(1, 'g', 'q'), abserr, relerr, dim) ); - convolutions_.push_back( new Convolution(gluon_nlo_, &SplittingFunction(0, 'g', 'q'), abserr, relerr, dim) ); - convolutions_.push_back( new Convolution(quark_nlo_, &SplittingFunction(0, 'q', 'q'), abserr, relerr, dim) ); - convolutions_.push_back( new Convolution(gluon_lo_, &ConvolutedSplittingFunctions(1, 'g', 'q', 'q'), abserr, relerr, dim) ); - convolutions_.push_back( new Convolution(gluon_lo_, &SplittingFunction(0, 'g', 'q'), abserr, relerr, dim) ); + convolutions_lmu1.push_back( new Convolution(gluon_lo_, &SplittingFunction(1, 'g', 'q'), abserr, relerr, dim) ); + convolutions_lmu1.push_back( new Convolution(gluon_nlo_, &SplittingFunction(0, 'g', 'q'), abserr, relerr, dim) ); + convolutions_lmu1.push_back( new Convolution(quark_nlo_, &SplittingFunction(0, 'q', 'q'), abserr, relerr, dim) ); + convolutions_lmu1.push_back( new Convolution(quark_nlo_, &(-2*Delta())) ); + + convolutions_lmu2.push_back( new Convolution(gluon_lo_, &(0.5*ConvolutedSplittingFunctions(1, 'g', 'g', 'g', 'q')), abserr, relerr, dim) ); + convolutions_lmu2.push_back( new Convolution(gluon_lo_, &(0.5*ConvolutedSplittingFunctions(1, 'q', 'q', 'g', 'q')), abserr, relerr, dim) ); + convolutions_lmu2.push_back( new Convolution(gluon_lo_, &(-3./2*SplittingFunction(0, 'g', 'q')), abserr, relerr, dim) ); } else { - convolutions_.push_back( new Convolution(gluon_lo_, &SplittingFunction(1, 'g', 'g'), abserr, relerr, dim) ); - convolutions_.push_back( new Convolution(quark_nlo_, &SplittingFunction(0, 'q', 'g'), abserr, relerr, dim) ); - convolutions_.push_back( new Convolution(gluon_nlo_, &SplittingFunction(0, 'g', 'g'), abserr, relerr, dim) ); - convolutions_.push_back( new MonteCarloDoubleConvolution(gluon_lo_, &SplittingFunction(0, 'g', 'g'), abserr, relerr, dim, MCcalls) ); - convolutions_.push_back( new Convolution(gluon_lo_, &ConvolutedSplittingFunctions(0, 'g', 'q', 'g'), abserr, relerr, dim) ); - convolutions_.push_back( new Convolution(gluon_lo_, &SplittingFunction(0, 'g', 'g'), abserr, relerr, dim) ); + convolutions_lmu1.push_back( new Convolution(gluon_lo_, &SplittingFunction(1, 'g', 'g'), abserr, relerr, dim) ); + convolutions_lmu1.push_back( new Convolution(gluon_lo_, &Delta()) ); + convolutions_lmu1.push_back( new Convolution(quark_nlo_, &SplittingFunction(0, 'q', 'g'), abserr, relerr, dim) ); + convolutions_lmu1.push_back( new Convolution(gluon_nlo_, &SplittingFunction(0, 'g', 'g'), abserr, relerr, dim) ); + convolutions_lmu1.push_back( new Convolution(gluon_nlo_, &Delta()) ); + + convolutions_lmu2.push_back( new MonteCarloDoubleConvolution(gluon_lo_, &SplittingFunction(0, 'g', 'g'), abserr, relerr, dim, MCcalls) ); + convolutions_lmu2.push_back( new Convolution(gluon_lo_, &ConvolutedSplittingFunctions(0, 'g', 'q', 'q', 'g'), abserr, relerr, dim) ); + convolutions_lmu2.push_back( new Convolution(gluon_lo_, &SplittingFunction(0, 'g', 'g'), abserr, relerr, dim) ); + convolutions_lmu1.push_back( new Convolution(gluon_lo_, &Delta()) ); } } + SetFunctions(); + } double ExactCoefficientFunction::fx(const double x, const double m2Q2, const double m2mu2, const int nf) const { + return MuIndependentTerms(x, m2Q2, nf) + MuDependentTerms(x, m2Q2, m2mu2, nf); +} - // O(as) - if (GetOrder() == 1 && GetKind() == '2') return C2_g1(x, m2Q2) ; - if (GetOrder() == 1 && GetKind() == 'L') return CL_g1(x, m2Q2) ; - - // O(as^2) - if (GetOrder() == 2 && GetKind() == '2' && GetChannel() == 'g') return C2_g2(x, m2Q2, m2mu2) ; - if (GetOrder() == 2 && GetKind() == '2' && GetChannel() == 'q') return C2_ps2(x, m2Q2, m2mu2) ; - if (GetOrder() == 2 && GetKind() == 'L' && GetChannel() == 'g') return CL_g2(x, m2Q2, m2mu2) ; - if (GetOrder() == 2 && GetKind() == 'L' && GetChannel() == 'q') return CL_ps2(x, m2Q2, m2mu2) ; +double ExactCoefficientFunction::MuIndependentTerms(const double x, const double m2Q2, const int nf) const { + return (this->*mu_indep_)(x, m2Q2, nf); +} - else { - cout << "Error: something has gone wrong!" << endl; - } +double ExactCoefficientFunction::MuDependentTerms(const double x, const double m2Q2, const double m2mu2, const int nf) const { + return (this->*mu_dep_)(x, m2Q2, m2mu2, nf); } -double ExactCoefficientFunction::MuIndependentTerms(const double x, const double m2Q2, const int nf) const { +void ExactCoefficientFunction::SetFunctions() { if (GetOrder() == 1) { - if (GetKind() == '2') return C2_g1(x, m2Q2) ; - else if (GetKind() == 'L') return CL_g1(x, m2Q2) ; - } else if (GetOrder() == 2) { if (GetKind() == '2') { - if (GetChannel() == 'g') return C2_g20(x, m2Q2) ; - else if (GetChannel() == 'q') return C2_ps20(x, m2Q2) ; - } else if (GetKind() == 'L') { - if (GetChannel() == 'g') return CL_g20(x, m2Q2) ; - else if (GetChannel() == 'q') return CL_ps20(x, m2Q2) ; + mu_indep_ = &ExactCoefficientFunction::C2_g1 ; + } else if (GetKind() =='L') { + mu_indep_ = &ExactCoefficientFunction::CL_g1 ; + } + mu_dep_ = &ExactCoefficientFunction::ZeroFunction ; + } else if (GetOrder() == 2) { + if (GetChannel() == 'q') { + if (GetKind() == '2') { + mu_indep_ = &ExactCoefficientFunction::C2_ps20 ; + } else if (GetKind() == 'L') { + mu_indep_ = &ExactCoefficientFunction::CL_ps20 ; + } + mu_dep_ = &ExactCoefficientFunction::C_ps2_MuDep ; + } else if (GetChannel() == 'g') { + if (GetKind() == '2') { + mu_indep_ = &ExactCoefficientFunction::C2_g20 ; + } else if (GetKind() == 'L') { + mu_indep_ = &ExactCoefficientFunction::CL_g20 ; + } + mu_dep_ = &ExactCoefficientFunction::C_g2_MuDep ; } - } else if (GetOrder() == 3){ - cout << "Error: mu independent terms are not known at O(as^3)!" << endl; - exit(-1); - } -} - -double ExactCoefficientFunction::MuDependentTerms(const double x, const double m2Q2, const double m2mu2, const int nf) const { - if (GetOrder() == 1) return 0.; - else if (GetOrder() == 2) { - if (GetChannel() == 'q') return C2 ; - else return (convolutions_[0] -> Convolute(x, m2Q2, nf) - beta(0, nf) * gluon_lo_ -> fx(x, m2Q2, static_cast(nan("")), nf)) * log(1. / m2mu2) ; } else if (GetOrder() == 3) { if (GetChannel() == 'q') { - + mu_dep_ = &ExactCoefficientFunction::C_ps3_MuDep ; + } else if (GetChannel() == 'g') { + mu_dep_ = &ExactCoefficientFunction::C_g3_MuDep ; } + mu_indep_ = nullptr; } } -// void ExactCoefficientFunction::Set_fx() { - -// // O(as) -// if (GetOrder() == 1 && GetKind() == '2') fx_ = C2_g1 ; -// if (GetOrder() == 1 && GetKind() == 'L') fx_ = CL_g1 ; - -// // O(as^2) -// if (GetOrder() == 2 && GetKind() == '2' && GetChannel() == 'g') fx_ = C2_g2 ; -// if (GetOrder() == 2 && GetKind() == '2' && GetChannel() == 'q') fx_ = C2_q2 ; -// if (GetOrder() == 2 && GetKind() == 'L' && GetChannel() == 'g') fx_ = CL_g2 ; -// if (GetOrder() == 2 && GetKind() == 'L' && GetChannel() == 'q') fx_ = CL_q2 ; - -// else { -// cout << "Error: something has gone wrong!" << endl; -// } -// } - void ExactCoefficientFunction::SetMethodFlag(const int& method_flag) { // check method_flag if (method_flag != 0 && method_flag != 1) { @@ -332,7 +324,7 @@ double ExactCoefficientFunction::C_ps21(const double x, const double m2Q2) const int nf = static_cast(nan("")); - return - convolutions_[0] -> Convolute(x, m2Q2, nf); + return - convolutions_lmu1[0] -> Convolute(x, m2Q2, nf); // The minus sign comes from the fact that in [arXiv:1205.5727] // the expansion is performed in terms of log(m^2/mu^2) // (even if it says the opposite) but we are @@ -396,7 +388,7 @@ double ExactCoefficientFunction::C_g21(const double x, const double m2Q2) const int nf = 1; // Put nf to 1 since the nf contribution cancels for any value of nf - return -(convolutions_[0] -> Convolute(x, m2Q2, nf) - gluon_lo_ -> MuIndependentTerms(x, m2Q2, nf) * beta(0, nf)); + return -(convolutions_lmu1[0] -> Convolute(x, m2Q2, nf) + convolutions_lmu1[1] -> Convolute(x, m2Q2, nf) * beta(0, nf)); } //==========================================================================================// @@ -452,8 +444,8 @@ double ExactCoefficientFunction::C_ps2_MuDep(const double x, const double m2Q2, double ExactCoefficientFunction::C_ps31(const double x, const double m2Q2, const int nf) const { return -( - convolutions_[0] -> Convolute(x, m2Q2, nf) + convolutions_[1] -> Convolute(x, m2Q2, nf) - + convolutions_[2] -> Convolute(x, m2Q2, nf) - 2. * beta(0, nf) * quark_nlo_ -> MuIndependentTerms(x, m2Q2, nf) + convolutions_lmu1[0] -> Convolute(x, m2Q2, nf) + convolutions_lmu1[1] -> Convolute(x, m2Q2, nf) + + convolutions_lmu1[2] -> Convolute(x, m2Q2, nf) + beta(0, nf) * convolutions_lmu1[3] -> Convolute(x, m2Q2, nf) ); } @@ -481,8 +473,8 @@ double ExactCoefficientFunction::C_ps31(const double x, const double m2Q2, const double ExactCoefficientFunction::C_ps32(const double x, const double m2Q2, const int nf) const { - return 0.5 * convolutions_[3] -> Convolute(x, m2Q2, nf) - - 3. / 2 * beta(0, nf) * convolutions_[4] -> Convolute(x, m2Q2, nf); + return convolutions_lmu2[0] -> Convolute(x, m2Q2, nf) + convolutions_lmu2[1] -> Convolute(x, m2Q2, nf) + + beta(0, nf) * convolutions_lmu2[2] -> Convolute(x, m2Q2, nf); } //==========================================================================================// @@ -510,9 +502,9 @@ double ExactCoefficientFunction::C_ps32(const double x, const double m2Q2, const double ExactCoefficientFunction::C_g31(const double x, const double m2Q2, const int nf) const { return -( - convolutions_[0] -> Convolute(x, m2Q2, nf) - beta(1, nf) * gluon_lo_ -> MuIndependentTerms(x, m2Q2, nf) - + convolutions_[1] -> Convolute(x, m2Q2, nf) + convolutions_[2] -> Convolute(x, m2Q2, nf) - - 2. * beta(0, nf) * gluon_nlo_ -> MuIndependentTerms(x, m2Q2, nf) + convolutions_lmu1[0] -> Convolute(x, m2Q2, nf) + beta(1, nf) * convolutions_lmu1[1] -> Convolute(x, m2Q2, nf) + + convolutions_lmu1[2] -> Convolute(x, m2Q2, nf) + convolutions_lmu1[3] -> Convolute(x, m2Q2, nf) + + beta(0, nf) * convolutions_lmu1[4] -> Convolute(x, m2Q2, nf) ); } @@ -555,9 +547,9 @@ double ExactCoefficientFunction::C_g32(const double x, const double m2Q2, const // exit(-1); // } - return 0.5 * convolutions_[3] -> Convolute(x, m2Q2, nf) + 0.5 * convolutions_[4] -> Convolute(x, m2Q2, nf) - - 3. / 2 * beta0 * convolutions_[5] -> Convolute(x, m2Q2, nf) - + beta0 * beta0 * gluon_lo_ -> MuIndependentTerms(x, m2Q2, nf); + return convolutions_lmu2[0] -> Convolute(x, m2Q2, nf) + convolutions_lmu2[1] -> Convolute(x, m2Q2, nf) + + beta0 * convolutions_lmu2[2] -> Convolute(x, m2Q2, nf) + + beta0 * beta0 * convolutions_lmu2[3] -> Convolute(x, m2Q2, nf); } //==========================================================================================// diff --git a/src/SplittingFunctions.cc b/src/SplittingFunctions.cc index 1f2a8aa..68a0e3b 100644 --- a/src/SplittingFunctions.cc +++ b/src/SplittingFunctions.cc @@ -1,24 +1,13 @@ #include "adani/SplittingFunctions.h" #include "adani/Constants.h" #include "adani/SpecialFunctions.h" + #include #include using std::cout; using std::endl; -SplittingFunction SplittingFunction::operator*(const double& rhs) const { - SplittingFunction res(order_, entry1_, entry2_); - res.SetMultFact(GetMultFact() * rhs); - return res; -} - -SplittingFunction operator*(const double& lhs, const SplittingFunction& rhs){ - SplittingFunction res(rhs.order_, res.entry1_, res.entry2_); - res.SetMultFact(rhs.GetMultFact() * lhs); - return res; -} - SplittingFunction::SplittingFunction(const int& order, const char& entry1, const char& entry2) : AbstractSplittingFunction(){ // check order @@ -42,94 +31,200 @@ SplittingFunction::SplittingFunction(const int& order, const char& entry1, const } entry2_ = entry2 ; - // check for not implemented splitting functions - if (order == 1 && entry1_ == 'q') { - cout << "Error: Pq" << entry2_ << " is not implemented at O(as)!" << endl ; - exit(-1); - } - + SetFunctions(); } double SplittingFunction::Regular(const double x, const int nf) const { - if (order_ == 0) { - if (entry1_ == 'g' && entry2_ == 'q') return Pgq0(x); - else if (entry1_ == 'q' && entry2_ == 'g') return Pqg0(x, nf); - else if (entry1_ == 'g' && entry2_ == 'g') return Pgg0reg(x); - else if (entry1_ == 'q' && entry2_ == 'q') return Pqq0reg(x); - } else if (order_ == 1) { - if (entry1_ == 'g' && entry2_ == 'q') return Pgq1(x, nf); - else if (entry1_ == 'g' && entry2_ == 'g') return Pgg1reg(x, nf); - } + return GetMultFact() * (this->*reg_)(x, nf); } double SplittingFunction::Singular(const double x, const int nf) const { - if (order_ == 0) { - if (entry1_ == 'g' && entry2_ == 'q') return 0.; - else if (entry1_ == 'q' && entry2_ == 'g') return 0.; - else if (entry1_ == 'g' && entry2_ == 'g') return Pgg0sing(x); - else if (entry1_ == 'q' && entry2_ == 'q') return Pqq0sing(x); - } else if (order_ == 1) { - if (entry1_ == 'g' && entry2_ == 'q') return 0.; - else if (entry1_ == 'g' && entry2_ == 'g') return Pgg1sing(x, nf); - } + return GetMultFact() * (this->*sing_)(x, nf); +} + +double SplittingFunction::SingularIntegrated(const double x, const int nf) const { + return GetMultFact() * (this->*sing_int_)(x, nf); } double SplittingFunction::Local(const int nf) const { - if (order_ == 0) { - if (entry1_ == 'g' && entry2_ == 'q') return 0.; - else if (entry1_ == 'q' && entry2_ == 'g') return 0.; - else if (entry1_ == 'g' && entry2_ == 'g') return Pgg0loc(nf); - else if (entry1_ == 'q' && entry2_ == 'q') return Pqq0loc(); - } else if (order_ == 1) { - if (entry1_ == 'g' && entry2_ == 'q') return 0.; - else if (entry1_ == 'g' && entry2_ == 'g') return Pgg1loc(nf); - } + return GetMultFact() * (this->*loc_)(nf); } -double SplittingFunction::SingularIntegrated(const double x, const int nf) const { - if (order_ == 0) { - if (entry1_ == 'g' && entry2_ == 'q') return 0.; - else if (entry1_ == 'q' && entry2_ == 'g') return 0.; - else if (entry1_ == 'g' && entry2_ == 'g') return Pgg0sing_integrated(x); - else if (entry1_ == 'q' && entry2_ == 'q') return Pqq0sing_integrated(x); +SplittingFunction SplittingFunction::operator*(const double& rhs) const { + SplittingFunction res(order_, entry1_, entry2_); + res.SetMultFact(GetMultFact() * rhs); + return res; +} + +SplittingFunction operator*(const double& lhs, const SplittingFunction& rhs){ + SplittingFunction res(rhs.order_, res.entry1_, res.entry2_); + res.SetMultFact(rhs.GetMultFact() * lhs); + return res; +} + +SplittingFunction SplittingFunction::operator/(const double& rhs) const { + SplittingFunction res(order_, entry1_, entry2_); + res.SetMultFact(GetMultFact() / rhs); + return res; +} + +void SplittingFunction::SetFunctions() { + + if (order_ == 0) { + if (entry1_ == 'g' && entry2_ == 'q') { + + reg_ = &SplittingFunction::Pgq0; + sing_ = &ZeroFunction_x_nf; + loc_ = &ZeroFunction_nf; + sing_int_ = &ZeroFunction_x_nf; + + } else if (entry1_ == 'q' && entry2_ == 'g') { + + reg_ = &SplittingFunction::Pqg0; + + sing_ = &ZeroFunction_x_nf; + loc_ = &ZeroFunction_nf; + sing_int_ = &ZeroFunction_x_nf; + + } else if (entry1_ == 'g' && entry2_ == 'g') { + + reg_ = &SplittingFunction::Pgg0reg; + sing_ = &SplittingFunction::Pgg0sing; + loc_ = &SplittingFunction::Pgg0loc; + sing_int_ = &SplittingFunction::Pgg0sing_integrated; + + } else if (entry1_ == 'q' && entry2_ == 'q') { + + reg_ = &SplittingFunction::Pqq0reg; + sing_ = &SplittingFunction::Pqq0sing; + loc_ = &SplittingFunction::Pqq0loc; + sing_int_ = &SplittingFunction::Pqq0sing_integrated; + + } } else if (order_ == 1) { - if (entry1_ == 'g' && entry2_ == 'q') return 0.; - else if (entry1_ == 'g' && entry2_ == 'g') return Pgg1sing_integrated(x, nf); + if (entry1_ == 'g' && entry2_ == 'q') { + + reg_ = &SplittingFunction::Pgq1; + sing_ = &ZeroFunction_x_nf; + loc_ = &ZeroFunction_nf; + sing_int_ = &ZeroFunction_x_nf; + + } else if (entry1_ == 'g' && entry2_ == 'g') { + + reg_ = &SplittingFunction::Pgg1reg; + sing_ = &SplittingFunction::Pgg1sing; + loc_ = &SplittingFunction::Pgg1loc; + sing_int_ = &SplittingFunction::Pgg1sing_integrated; + + } + } else { + cout << "Error: P" << entry1_ << entry2_ << order_ << " is not implemented!" << endl ; + exit(-1); } + } -// SplittingFunction SplittingFunction::operator+(const SplittingFunction& splitfunc) const { - -// } +ConvolutedSplittingFunctions::ConvolutedSplittingFunctions(const int& order, const char& entry1, const char& entry2, const char& entry3, const char& entry4) : AbstractSplittingFunction() { + + // check order + if (order != 0 && order !=1) { + cout << "Error: order must be 0 or 1. Got " << order << endl ; + exit(-1) ; + } + order_ = order ; + + // check entry1 + if (entry1 != 'g' && entry2 != 'q') { + cout << "Error: entry1 must be g or q. Got " << entry1 << endl ; + exit(-1) ; + } + entry1_ = entry1 ; + + // check entry2 + if (entry2 != 'g' && entry2 != 'q') { + cout << "Error: entry2 must be g or q. Got " << entry2 << endl ; + exit(-1) ; + } + entry2_ = entry2 ; + + // check entry3 + if (entry3 != 'g' && entry3 != 'q') { + cout << "Error: entry3 must be g or q. Got " << entry3 << endl ; + exit(-1) ; + } + entry3_ = entry3; -ConvolutedSplittingFunctions::ConvolutedSplittingFunctions(const int& order, const char& entry1, const char& entry2, const char& entry3) : SplittingFunction(order, entry1, entry2) { // check entry3 if (entry3 != 'g' && entry3 != 'q') { cout << "Error: entry3 must be g or q. Got " << entry3 << endl ; exit(-1) ; } entry3_ = entry3; + + SetFunctions(); } double ConvolutedSplittingFunctions::Regular(const double x, const int nf) const { - - if (GetEntry1() == 'g' && GetEntry2() == 'q' && GetEntry3() == 'g') return Pgq0_x_Pqg0(x, nf) ; - if (GetEntry1() == 'g' && GetEntry2() == 'q' && GetEntry3() == 'q') return Pgg0_x_Pgq0(x, nf) + Pqq0_x_Pgq0(x) ; - + return GetMultFact() * (this->*reg_)(x, nf); } ConvolutedSplittingFunctions ConvolutedSplittingFunctions::operator*(const double& rhs) const { - ConvolutedSplittingFunctions res(GetOrder(), GetEntry1(), GetEntry2(), entry3_); + ConvolutedSplittingFunctions res(order_, entry1_, entry2_, entry3_, entry4_); res.SetMultFact(GetMultFact() * rhs); return res; } ConvolutedSplittingFunctions operator*(const double& lhs, const ConvolutedSplittingFunctions& rhs){ - ConvolutedSplittingFunctions res(rhs.GetOrder(), res.GetEntry1(), res.GetEntry2(), res.entry3_); + ConvolutedSplittingFunctions res(rhs.order_, res.entry1_, res.entry2_, res.entry3_, res.entry4_); + res.SetMultFact(rhs.GetMultFact() * lhs); + return res; +} + +ConvolutedSplittingFunctions ConvolutedSplittingFunctions::operator/(const double& rhs) const { + ConvolutedSplittingFunctions res(order_, entry1_, entry2_, entry3_, entry4_); + res.SetMultFact(GetMultFact() / rhs); + return res; +} + +void ConvolutedSplittingFunctions::SetFunctions() { + + if (entry1_ == 'g' && entry2_ == 'q' && entry3_ == 'q' && entry4_ == 'g') { + + reg_ = &ConvolutedSplittingFunctions::Pgq0_x_Pqg0; + + } else if (entry1_ == 'g' && entry2_ == 'g' && entry3_ == 'g' && entry4_ == 'q') { + + reg_ = &ConvolutedSplittingFunctions::Pgg0_x_Pgq0; + + } else if (entry1_ == 'q' && entry2_ == 'q' && entry3_ == 'g' && entry4_ == 'q') { + + reg_ = &ConvolutedSplittingFunctions::Pqq0_x_Pgq0; + + } else { + cout << "Error: P" << entry1_ << entry2_ << " x P" << entry3_ << entry4_ << " is not implemented!" << endl ; + exit(-1); + } + +} + +Delta Delta::operator*(const double& rhs) const { + Delta res = Delta(); + res.SetMultFact(GetMultFact() * rhs); + return res; +} + +Delta operator*(const double& lhs, const Delta& rhs){ + Delta res = Delta(); res.SetMultFact(rhs.GetMultFact() * lhs); return res; } +Delta Delta::operator/(const double& rhs) const { + Delta res = Delta(); + res.SetMultFact(GetMultFact() / rhs); + return res; +} + //==========================================================================================// // Gluon-quark splitting functions O(alpha_s) without color factors // @@ -389,4 +484,4 @@ double ConvolutedSplittingFunctions::Pgq0_x_Pqg0(const double x, const int nf) c return 4. * CF * nf * (1. + 4. / 3 / x - x - 4. * x * x / 3 + 2. * (1 + x) * log(x)); -} \ No newline at end of file +}