From 5d5f0b842bd13bf117ae6598c5747c8a290edb0e Mon Sep 17 00:00:00 2001 From: "Tanjona R. Rabemananjara" Date: Sat, 5 Jun 2021 21:49:15 +0200 Subject: [PATCH] fix some inconsistencies in formatting --- .clang-format | 79 ++-- .github/workflows/clang-format.yml | 7 +- dSigmadpt.cpp | 19 +- dSigmadptN.cpp | 19 +- include/AnomalousDim.h | 12 +- include/CombinedRes.h | 10 +- include/ComplexDefs.h | 9 +- include/HarmonicSum.h | 22 +- include/MellinFunc.h | 226 ++++++------ include/SmallptExp.h | 10 +- include/ThresExp.h | 10 +- meson.build | 4 +- runcards/Mellin-HpT-as-N.yaml | 4 +- src/AnomalousDim.cpp | 81 +++-- src/CombinedRes.cpp | 16 +- src/ComplexDefs.cpp | 121 +++---- src/HarmonicSum.cpp | 562 ++++++++++++++--------------- src/MellinFunc.cpp | 104 +++--- src/SmallptExp.cpp | 170 ++++----- src/ThresExp.cpp | 281 ++++++++------- 20 files changed, 912 insertions(+), 854 deletions(-) diff --git a/.clang-format b/.clang-format index 03ce0b4..70da6ad 100644 --- a/.clang-format +++ b/.clang-format @@ -1,13 +1,13 @@ --- Language: Cpp -# BasedOnStyle: LLVM -AccessModifierOffset: -2 +# BasedOnStyle: Google +AccessModifierOffset: -1 AlignAfterOpenBracket: Align AlignConsecutiveMacros: false AlignConsecutiveAssignments: false AlignConsecutiveBitFields: false AlignConsecutiveDeclarations: false -AlignEscapedNewlines: Right +AlignEscapedNewlines: Left AlignOperands: Align AlignTrailingComments: true AllowAllArgumentsOnNextLine: true @@ -18,12 +18,12 @@ AllowShortBlocksOnASingleLine: Never AllowShortCaseLabelsOnASingleLine: false AllowShortFunctionsOnASingleLine: All AllowShortLambdasOnASingleLine: All -AllowShortIfStatementsOnASingleLine: Never -AllowShortLoopsOnASingleLine: false +AllowShortIfStatementsOnASingleLine: WithoutElse +AllowShortLoopsOnASingleLine: true AlwaysBreakAfterDefinitionReturnType: None AlwaysBreakAfterReturnType: None -AlwaysBreakBeforeMultilineStrings: false -AlwaysBreakTemplateDeclarations: MultiLine +AlwaysBreakBeforeMultilineStrings: true +AlwaysBreakTemplateDeclarations: Yes BinPackArguments: true BinPackParameters: true BraceWrapping: @@ -57,12 +57,12 @@ BreakStringLiterals: true ColumnLimit: 80 CommentPragmas: '^ IWYU pragma:' CompactNamespaces: false -ConstructorInitializerAllOnOneLineOrOnePerLine: false +ConstructorInitializerAllOnOneLineOrOnePerLine: true ConstructorInitializerIndentWidth: 4 ContinuationIndentWidth: 4 Cpp11BracedListStyle: true DeriveLineEnding: true -DerivePointerAlignment: false +DerivePointerAlignment: true DisableFormat: false ExperimentalAutoDetectBinPacking: false FixNamespaceComments: true @@ -70,20 +70,23 @@ ForEachMacros: - foreach - Q_FOREACH - BOOST_FOREACH -IncludeBlocks: Preserve +IncludeBlocks: Regroup IncludeCategories: - - Regex: '^"(llvm|llvm-c|clang|clang-c)/' + - Regex: '^' Priority: 2 SortPriority: 0 - - Regex: '^(<|"(gtest|gmock|isl|json)/)' - Priority: 3 + - Regex: '^<.*\.h>' + Priority: 1 + SortPriority: 0 + - Regex: '^<.*' + Priority: 2 SortPriority: 0 - Regex: '.*' - Priority: 1 + Priority: 3 SortPriority: 0 -IncludeIsMainRegex: '(Test)?$' +IncludeIsMainRegex: '([-_](test|unittest))?$' IncludeIsMainSourceRegex: '' -IndentCaseLabels: false +IndentCaseLabels: true IndentCaseBlocks: false IndentGotoLabels: true IndentPPDirectives: None @@ -93,25 +96,55 @@ IndentWrappedFunctionNames: false InsertTrailingCommas: None JavaScriptQuotes: Leave JavaScriptWrapImports: true -KeepEmptyLinesAtTheStartOfBlocks: true +KeepEmptyLinesAtTheStartOfBlocks: false MacroBlockBegin: '' MacroBlockEnd: '' MaxEmptyLinesToKeep: 1 NamespaceIndentation: None -ObjCBinPackProtocolList: Auto +ObjCBinPackProtocolList: Never ObjCBlockIndentWidth: 2 ObjCBreakBeforeNestedBlockParam: true ObjCSpaceAfterProperty: false ObjCSpaceBeforeProtocolList: true PenaltyBreakAssignment: 2 -PenaltyBreakBeforeFirstCallParameter: 19 +PenaltyBreakBeforeFirstCallParameter: 1 PenaltyBreakComment: 300 PenaltyBreakFirstLessLess: 120 PenaltyBreakString: 1000 PenaltyBreakTemplateDeclaration: 10 PenaltyExcessCharacter: 1000000 -PenaltyReturnTypeOnItsOwnLine: 60 -PointerAlignment: Right +PenaltyReturnTypeOnItsOwnLine: 200 +PointerAlignment: Left +RawStringFormats: + - Language: Cpp + Delimiters: + - cc + - CC + - cpp + - Cpp + - CPP + - 'c++' + - 'C++' + CanonicalDelimiter: '' + BasedOnStyle: google + - Language: TextProto + Delimiters: + - pb + - PB + - proto + - PROTO + EnclosingFunctions: + - EqualsProto + - EquivToProto + - PARSE_PARTIAL_TEXT_PROTO + - PARSE_TEST_PROTO + - PARSE_TEXT_PROTO + - ParseTextOrDie + - ParseTextProtoOrDie + - ParseTestProto + - ParsePartialTestProto + CanonicalDelimiter: '' + BasedOnStyle: google ReflowComments: true SortIncludes: true SortUsingDeclarations: true @@ -126,7 +159,7 @@ SpaceBeforeParens: ControlStatements SpaceBeforeRangeBasedForLoopColon: true SpaceInEmptyBlock: false SpaceInEmptyParentheses: false -SpacesBeforeTrailingComments: 1 +SpacesBeforeTrailingComments: 2 SpacesInAngles: false SpacesInConditionalStatement: false SpacesInContainerLiterals: true @@ -134,7 +167,7 @@ SpacesInCStyleCastParentheses: false SpacesInParentheses: false SpacesInSquareBrackets: false SpaceBeforeSquareBrackets: false -Standard: Latest +Standard: Auto StatementMacros: - Q_UNUSED - QT_REQUIRE_VERSION diff --git a/.github/workflows/clang-format.yml b/.github/workflows/clang-format.yml index 398a786..a776e5a 100644 --- a/.github/workflows/clang-format.yml +++ b/.github/workflows/clang-format.yml @@ -1,7 +1,10 @@ name: Run clang-format Linter on: - pull_request: + push: + branches: + - main + - develop workflow_dispatch: jobs: @@ -16,7 +19,7 @@ jobs: extensions: 'h,cpp' clangFormatVersion: 12 inplace: True - style: llvm + style: google - uses: EndBug/add-and-commit@v4 with: default_author: github_actions diff --git a/dSigmadpt.cpp b/dSigmadpt.cpp index 18d3c72..27c5906 100644 --- a/dSigmadpt.cpp +++ b/dSigmadpt.cpp @@ -17,13 +17,14 @@ * ===================================================================================== */ +#include +#include + #include #include #include #include #include -#include -#include #include #include @@ -76,10 +77,8 @@ int main(int argc, char *argv[]) { "_combined_aspt.dat"}; try { - if (order < 0 || order > 1) - throw err_message(); - if (channel < 0 || channel > 5) - throw err_message(); + if (order < 0 || order > 1) throw err_message(); + if (channel < 0 || channel > 5) throw err_message(); filename += ord_fixod[order]; filename += par_chanl[channel]; filename += matsch[scheme]; @@ -90,14 +89,14 @@ int main(int argc, char *argv[]) { // Factors for Born cross-section double factor; - double gf = 1.16637e-5; // Fermi Constant - double gevpb = 3.8937966e8; // GeV to pb + double gf = 1.16637e-5; // Fermi Constant + double gevpb = 3.8937966e8; // GeV to pb if (inorm == 1) { std::cout << "ERROR in inorm!" << std::endl; - exit(EXIT_FAILURE); // TODO: complete implementation! + exit(EXIT_FAILURE); // TODO: complete implementation! } else if (inorm == 0) { - factor = gf / 288. / M_PI / sqrt(2.); // large-top mass limit + factor = gf / 288. / M_PI / sqrt(2.); // large-top mass limit } else { std::cout << "ERROR in inorm!" << std::endl; exit(EXIT_FAILURE); diff --git a/dSigmadptN.cpp b/dSigmadptN.cpp index 3351f2f..d6f6781 100644 --- a/dSigmadptN.cpp +++ b/dSigmadptN.cpp @@ -17,13 +17,14 @@ * ===================================================================================== */ +#include +#include + #include #include #include #include #include -#include -#include #include #include @@ -73,10 +74,8 @@ int main(int argc, char *argv[]) { "_combined_asN.dat"}; try { - if (order < 0 || order > 1) - throw err_message(); - if (channel < 0 || channel > 5) - throw err_message(); + if (order < 0 || order > 1) throw err_message(); + if (channel < 0 || channel > 5) throw err_message(); filename += ord_fixod[order]; filename += par_chanl[channel]; filename += matsch[scheme]; @@ -87,14 +86,14 @@ int main(int argc, char *argv[]) { // Factors for Born cross-section double factor; - double gf = 1.16637e-5; // Fermi Constant - double gevpb = 3.8937966e8; // GeV to pb + double gf = 1.16637e-5; // Fermi Constant + double gevpb = 3.8937966e8; // GeV to pb if (inorm == 1) { std::cout << "ERROR in inorm!" << std::endl; - exit(EXIT_FAILURE); // TODO: complete implementation! + exit(EXIT_FAILURE); // TODO: complete implementation! } else if (inorm == 0) { - factor = gf / 288. / M_PI / sqrt(2.); // large-top mass limit + factor = gf / 288. / M_PI / sqrt(2.); // large-top mass limit } else { std::cout << "ERROR in inorm!" << std::endl; exit(EXIT_FAILURE); diff --git a/include/AnomalousDim.h b/include/AnomalousDim.h index ea35a9e..7e41c70 100644 --- a/include/AnomalousDim.h +++ b/include/AnomalousDim.h @@ -1,15 +1,15 @@ #pragma once -#include -#include #include #include -#include #include +#include +#include +#include + #include "./ComplexDefs.h" #include "./HarmonicSum.h" - #include "higgs-fo/params.h" using namespace std; @@ -26,7 +26,7 @@ struct gamma1sums { }; class AnomDimensions { -public: + public: AnomDimensions(void *params); virtual ~AnomDimensions(); @@ -43,7 +43,7 @@ class AnomDimensions { long double _NN; }; -private: + private: HSum HAP; gamma1sums g1s; diff --git a/include/CombinedRes.h b/include/CombinedRes.h index 88205fb..4b0d838 100644 --- a/include/CombinedRes.h +++ b/include/CombinedRes.h @@ -18,22 +18,22 @@ #pragma once -#include -#include #include #include #include + +#include +#include #include #include #include "./ComplexDefs.h" #include "./SmallptExp.h" #include "./ThresExp.h" - #include "higgs-fo/partonic.h" class CombinedRes { -public: + public: CombinedRes(int order, int channel, std::string pdfname, void *params); virtual ~CombinedRes(); @@ -41,7 +41,7 @@ class CombinedRes { std::complex CombinedResExpr(std::complex N, long double pt, int scheme); -private: + private: int ORD, EXACT_ORD; /* long double LF, LR, LQ; */ long double QS2, MH2, MUF2, MUR2; diff --git a/include/ComplexDefs.h b/include/ComplexDefs.h index 7eab075..3ead81b 100644 --- a/include/ComplexDefs.h +++ b/include/ComplexDefs.h @@ -1,6 +1,5 @@ #pragma once -#include #include #include #include @@ -8,6 +7,8 @@ #include #include #include + +#include #include #include #include @@ -68,9 +69,9 @@ std::complex operator*(const unsigned int n, const std::complex &z); std::complex operator/(const unsigned int n, const std::complex &z); -std::complex -operator*(const std::vector> &c1, - const std::vector> &c2); +std::complex operator*( + const std::vector> &c1, + const std::vector> &c2); bool operator==(const std::complex &z, const int n); bool operator!=(const std::complex &z, const int n); diff --git a/include/HarmonicSum.h b/include/HarmonicSum.h index b82fd55..0362008 100644 --- a/include/HarmonicSum.h +++ b/include/HarmonicSum.h @@ -1,12 +1,13 @@ #pragma once -#include -#include -#include #include #include #include #include + +#include +#include +#include #include #include #include @@ -17,7 +18,7 @@ using namespace std; class HSum { -public: + public: HSum(bool verbose = false, bool testinterfun = false, bool testharmsums = false); virtual ~HSum(); @@ -29,7 +30,7 @@ class HSum { std::complex HS(int i, int j, int k, int m, std::complex N); -private: + private: bool _verbose; bool _testinterpolatedfunction; bool _testharmonicsums; @@ -232,8 +233,7 @@ class HSum { std::complex sum, ris; z -= 1; sum = p[0]; - for (int i = 1; i < (g + 2); i++) - sum += p[i] / (z + i); + for (int i = 1; i < (g + 2); i++) sum += p[i] / (z + i); ris = 0.5 * gsl_sf_log(2 * M_PI) + (z + 0.5) * log(z + g + 0.5) - (z + g + 0.5) + log(sum); return ris; @@ -248,10 +248,10 @@ class HSum { if (i == 0) { std::complex SUB = 0.; std::complex ZZ = z; - if (std::abs(std::imag(ZZ)) < 10.) { // if too close to the real axis... + if (std::abs(std::imag(ZZ)) < 10.) { // if too close to the real axis... label1: if (std::real(ZZ) < - 10.) { // ...use recurrence relation to push real(z) large enough + 10.) { // ...use recurrence relation to push real(z) large enough SUB = SUB - 1. / ZZ; ZZ = ZZ + 1.; goto label1; @@ -266,14 +266,14 @@ class HSum { int K1, K2; std::complex SUB = 0., SUBM; std::complex ZZ = z; - if (std::abs(std::imag(ZZ)) < 10.) { // if too close to the real axis... + if (std::abs(std::imag(ZZ)) < 10.) { // if too close to the real axis... label2: SUBM = -1. / ZZ; for (K1 = 1; K1 <= i; K1++) { SUBM = -SUBM * K1 / ZZ; } if (std::real(ZZ) < - 10.) { // ...use recurrence relation to push real(z) large enough + 10.) { // ...use recurrence relation to push real(z) large enough SUB = SUB + SUBM; ZZ = ZZ + 1.; goto label2; diff --git a/include/MellinFunc.h b/include/MellinFunc.h index 1dc46b3..266c67f 100644 --- a/include/MellinFunc.h +++ b/include/MellinFunc.h @@ -1,9 +1,10 @@ #pragma once -#include #include #include #include + +#include #include #include "ComplexDefs.h" @@ -12,124 +13,125 @@ using namespace std; class MellinFunc { -public: + public: MellinFunc(); virtual ~MellinFunc(); - std::complex D0(std::complex N); // [1/(1-z)]_+ - std::complex - D1(std::complex N); // [Log(1-z)/(1-z)]_+ - std::complex - D2(std::complex N); // [Log(1-z)^2/(1-z)]_+ - std::complex - D3(std::complex N); // [Log(1-z)^3/(1-z)]_+ + std::complex D0(std::complex N); // [1/(1-z)]_+ + std::complex D1( + std::complex N); // [Log(1-z)/(1-z)]_+ + std::complex D2( + std::complex N); // [Log(1-z)^2/(1-z)]_+ + std::complex D3( + std::complex N); // [Log(1-z)^3/(1-z)]_+ - std::complex Li2z(std::complex N); // Li2(z) - std::complex Li2mz(std::complex N); // Li2(-z) - std::complex S12z(std::complex N); // S12(z) - std::complex S12mz(std::complex N); // S12(mz) - std::complex S12z2(std::complex N); // S12(z^2) - std::complex Li3z(std::complex N); // Li3(z) - std::complex Li3mz(std::complex N); // Li3(mz) - std::complex - Li2zLogz(std::complex N); // Li2(z)Log(z) - std::complex Logminus(std::complex N); // Log(1-z) - std::complex Logplus(std::complex N); // Log(1+z) - std::complex Logz(std::complex N); // Log(z) - std::complex Logz2(std::complex N); // Log(z)^2 - std::complex Logz3(std::complex N); // Log(z)^3 - std::complex - LogzLogminus(std::complex N); // Log(z)Log(1-z) - std::complex - LogzLogminus2(std::complex N); // Log(z)Log(1-z)^2 - std::complex - Logz2Logminus(std::complex N); // Log(z)^2Log(1-z) - std::complex - Logz2minus(std::complex N); // log(z)^2/(1-z) - std::complex - Logminus2(std::complex N); // Log(1-z)^2 - std::complex - Logminus3(std::complex N); // Log(1-z)^3 - std::complex - LogzLogminusminus(std::complex N); // Log(z)Log(1-z)/(1-z) - std::complex - Logzminus(std::complex N); // log(z)/(1-z) - std::complex - Logzminusplus(std::complex N); // log(z)/((1-z)(1+z)) - std::complex plus(std::complex N); // 1/(1+z) - std::complex - Li2minusminus(std::complex N); // Li2(1-z)/(1-z) - std::complex - S12zregminus(std::complex N); // (S12(z)-Zeta(3))/(1-z) - std::complex - Li2zregminus(std::complex N); // (Li2(z)-Zeta(2))/(1-z) - std::complex - Li3zregminus(std::complex N); // (Li3(z)-Zeta(3))/(1-z) + std::complex Li2z(std::complex N); // Li2(z) + std::complex Li2mz(std::complex N); // Li2(-z) + std::complex S12z(std::complex N); // S12(z) + std::complex S12mz(std::complex N); // S12(mz) + std::complex S12z2(std::complex N); // S12(z^2) + std::complex Li3z(std::complex N); // Li3(z) + std::complex Li3mz(std::complex N); // Li3(mz) + std::complex Li2zLogz( + std::complex N); // Li2(z)Log(z) + std::complex Logminus(std::complex N); // Log(1-z) + std::complex Logplus(std::complex N); // Log(1+z) + std::complex Logz(std::complex N); // Log(z) + std::complex Logz2(std::complex N); // Log(z)^2 + std::complex Logz3(std::complex N); // Log(z)^3 + std::complex LogzLogminus( + std::complex N); // Log(z)Log(1-z) + std::complex LogzLogminus2( + std::complex N); // Log(z)Log(1-z)^2 + std::complex Logz2Logminus( + std::complex N); // Log(z)^2Log(1-z) + std::complex Logz2minus( + std::complex N); // log(z)^2/(1-z) + std::complex Logminus2( + std::complex N); // Log(1-z)^2 + std::complex Logminus3( + std::complex N); // Log(1-z)^3 + std::complex LogzLogminusminus( + std::complex N); // Log(z)Log(1-z)/(1-z) + std::complex Logzminus( + std::complex N); // log(z)/(1-z) + std::complex Logzminusplus( + std::complex N); // log(z)/((1-z)(1+z)) + std::complex plus(std::complex N); // 1/(1+z) + std::complex Li2minusminus( + std::complex N); // Li2(1-z)/(1-z) + std::complex S12zregminus( + std::complex N); // (S12(z)-Zeta(3))/(1-z) + std::complex Li2zregminus( + std::complex N); // (Li2(z)-Zeta(2))/(1-z) + std::complex Li3zregminus( + std::complex N); // (Li3(z)-Zeta(3))/(1-z) std::complex Li2zregLogminusminus( - std::complex N); //(Li2(z)-Zeta(2))Log(1-z)/(1-z) - std::complex - Li3zplus(std::complex N); // Li3(z)/(1+z) - std::complex - S12z2plus(std::complex N); // S12(z^2)/(1+z) - std::complex - S12mzplus(std::complex N); // S12(-z)/(1+z) - std::complex - S12zplus(std::complex N); // S12(z)/(1+z) - std::complex - Li3mzplus(std::complex N); // Li3(-z)/(1+z) - std::complex - Li2zLogminus(std::complex N); // Li2(z)Log(1-z) - std::complex - Li2zLogplus(std::complex N); // Li2(z)Log(1+z) - std::complex - Li2mzLogplus(std::complex N); // Li2(-z)Log(1+z) - std::complex - Li2mzLogz(std::complex N); // Li2(-z)Log(z) - std::complex - Li2zLogzplusminus(std::complex N); // Li2(z)Log(z)/((1+z)(1-z)) - std::complex - Li2mzLogzplus(std::complex N); // Li2(-z)Log(z)/(1+z) - std::complex - Li2mzLogminus2plus(std::complex N); // Li2(-z)Log[1-z]^2/(1+z) - std::complex - Li2mzLogplusplus(std::complex N); // Li2(-z)Log[1+z]/(1+z) - std::complex - Logz2Logminusminus(std::complex N); // Log(z)^2Log(1-z)/(1-z) - std::complex - LogzLogminus2minus(std::complex N); // Log(z)Log(1-z)^2/(1-z) - std::complex - Logz3minusplus(std::complex N); // Log(z)^3 /((1-z)(1+z)) - std::complex - Logplusplus(std::complex N); // Log(1+z)/(1+z) - std::complex - Li2zLogplusplus(std::complex N); // Li2(z)Log(1+z)/(1+z) - std::complex - Logz2Logplusplus(std::complex N); // Log(1+z)Log(z)^2/(1+z) - std::complex - LogzLogplus(std::complex N); // Log(z)Log(1+z) - std::complex - Logz2Logplus(std::complex N); // Log(z)^2Log(1+z) - std::complex - LogzLogplus2(std::complex N); // Log(z)Log(1+z)^2 - std::complex - LogzLogplus2plus(std::complex N); // Log(z)Log(1+z)^2/(1+z) + std::complex N); //(Li2(z)-Zeta(2))Log(1-z)/(1-z) + std::complex Li3zplus( + std::complex N); // Li3(z)/(1+z) + std::complex S12z2plus( + std::complex N); // S12(z^2)/(1+z) + std::complex S12mzplus( + std::complex N); // S12(-z)/(1+z) + std::complex S12zplus( + std::complex N); // S12(z)/(1+z) + std::complex Li3mzplus( + std::complex N); // Li3(-z)/(1+z) + std::complex Li2zLogminus( + std::complex N); // Li2(z)Log(1-z) + std::complex Li2zLogplus( + std::complex N); // Li2(z)Log(1+z) + std::complex Li2mzLogplus( + std::complex N); // Li2(-z)Log(1+z) + std::complex Li2mzLogz( + std::complex N); // Li2(-z)Log(z) + std::complex Li2zLogzplusminus( + std::complex N); // Li2(z)Log(z)/((1+z)(1-z)) + std::complex Li2mzLogzplus( + std::complex N); // Li2(-z)Log(z)/(1+z) + std::complex Li2mzLogminus2plus( + std::complex N); // Li2(-z)Log[1-z]^2/(1+z) + std::complex Li2mzLogplusplus( + std::complex N); // Li2(-z)Log[1+z]/(1+z) + std::complex Logz2Logminusminus( + std::complex N); // Log(z)^2Log(1-z)/(1-z) + std::complex LogzLogminus2minus( + std::complex N); // Log(z)Log(1-z)^2/(1-z) + std::complex Logz3minusplus( + std::complex N); // Log(z)^3 /((1-z)(1+z)) + std::complex Logplusplus( + std::complex N); // Log(1+z)/(1+z) + std::complex Li2zLogplusplus( + std::complex N); // Li2(z)Log(1+z)/(1+z) + std::complex Logz2Logplusplus( + std::complex N); // Log(1+z)Log(z)^2/(1+z) + std::complex LogzLogplus( + std::complex N); // Log(z)Log(1+z) + std::complex Logz2Logplus( + std::complex N); // Log(z)^2Log(1+z) + std::complex LogzLogplus2( + std::complex N); // Log(z)Log(1+z)^2 + std::complex LogzLogplus2plus( + std::complex N); // Log(z)Log(1+z)^2/(1+z) - std::complex - Logplus3plus(std::complex N); // Log(1+z)^3/(1+z) - std::complex - Li3zregminusplus(std::complex N); // (Li3-Zeta[3])/(1-z)/(1+z) - std::complex - Li3zoverplusplus(std::complex N); // (Li3(z/(1+z))/(1+z) - std::complex Logplus3(std::complex N); // Log[1+z]^3 - std::complex Li2z2(std::complex N); // Li2(z^2) - std::complex - Li2z2Logz(std::complex N); // Li2(z^2) Log(z) - std::complex Li3z2(std::complex N); // Li3(z^2) - std::complex - Li3overplus(std::complex N); // Li3(1/(1+z)) - std::complex Li2minus(std::complex N); // Li2(1-z) + std::complex Logplus3plus( + std::complex N); // Log(1+z)^3/(1+z) + std::complex Li3zregminusplus( + std::complex N); // (Li3-Zeta[3])/(1-z)/(1+z) + std::complex Li3zoverplusplus( + std::complex N); // (Li3(z/(1+z))/(1+z) + std::complex Logplus3( + std::complex N); // Log[1+z]^3 + std::complex Li2z2(std::complex N); // Li2(z^2) + std::complex Li2z2Logz( + std::complex N); // Li2(z^2) Log(z) + std::complex Li3z2(std::complex N); // Li3(z^2) + std::complex Li3overplus( + std::complex N); // Li3(1/(1+z)) + std::complex Li2minus(std::complex N); // Li2(1-z) -private: + private: HSum H; long double zeta2; long double zeta3; diff --git a/include/SmallptExp.h b/include/SmallptExp.h index 0ab9f8d..b1bf341 100644 --- a/include/SmallptExp.h +++ b/include/SmallptExp.h @@ -18,23 +18,23 @@ #pragma once -#include #include -#include #include #include #include + +#include +#include #include #include #include "./AnomalousDim.h" #include "./ComplexDefs.h" #include "./MellinFunc.h" - #include "higgs-fo/params.h" class SmallptExp { -public: + public: SmallptExp(int order, int channel, void *params); virtual ~SmallptExp(); @@ -42,7 +42,7 @@ class SmallptExp { std::complex SmallptExpExpr(std::complex N, long double pt); -private: + private: AnomDimensions AD; int NC, NF, CA, ORD, CHANNEL; diff --git a/include/ThresExp.h b/include/ThresExp.h index c2ba72f..a98fffa 100644 --- a/include/ThresExp.h +++ b/include/ThresExp.h @@ -18,22 +18,22 @@ #pragma once -#include #include -#include #include #include #include + +#include +#include #include #include #include "./ComplexDefs.h" #include "./MellinFunc.h" - #include "higgs-fo/params.h" class ThresExp { -public: + public: ThresExp(int order, int channel, void *params); virtual ~ThresExp(); @@ -62,7 +62,7 @@ class ThresExp { long double GOgqqH(long double); long double GOqqgH(long double); -private: + private: int NC, NF, CA, ORD, CHANNEL; long double LF, LR, LQ; long double MH2, MUF2, MUR2, SROOT; diff --git a/meson.build b/meson.build index 38fe664..ccf0a61 100644 --- a/meson.build +++ b/meson.build @@ -31,7 +31,7 @@ subdir('include') subdir('src') # create executable for cross section as a function of pt -higgsexe = executable( +higgsexep = executable( 'higgs-pt', 'dSigmadpt.cpp', include_directories : inc, dependencies : [glib_dep, higgsfo_lib, yaml_lib], @@ -40,7 +40,7 @@ higgsexe = executable( ) # create executable for cross section as a function of N -higgsexe = executable( +higgsexen = executable( 'higgs-n', 'dSigmadptN.cpp', include_directories : inc, dependencies : [glib_dep, higgsfo_lib, yaml_lib], diff --git a/runcards/Mellin-HpT-as-N.yaml b/runcards/Mellin-HpT-as-N.yaml index 50b4a49..b38c6fd 100644 --- a/runcards/Mellin-HpT-as-N.yaml +++ b/runcards/Mellin-HpT-as-N.yaml @@ -14,7 +14,7 @@ scheme : 2 # 1<==>(NLO) order : 1 -# pt bins +# Mellin N bins Nmin : 1.5 Nmax : 10 Nbin : 0.1 @@ -43,5 +43,5 @@ sroot : 13000 # square root of the center of mass energy in GeV # 3<==>qqb channel : 0 -# mellin moment +# pt value pt : 5 diff --git a/src/AnomalousDim.cpp b/src/AnomalousDim.cpp index ebfdee4..d8eb443 100644 --- a/src/AnomalousDim.cpp +++ b/src/AnomalousDim.cpp @@ -1,4 +1,5 @@ #include "AnomalousDim.h" + #include // Construct the Anomalous dimension Computation @@ -17,10 +18,10 @@ AnomDimensions::~AnomDimensions() {} // Functions that computes the Anomalous Dimensions void AnomDimensions::ComputeGamma(std::complex N, int order) { if (order >= 0) { - gg0 = gammagg0(N) * M_PIl; // \gamma_gg^0 - gq0 = gammagS0(N) * M_PIl; // \gamma_gq^0 - qg0 = gammaSg0(N) * M_PIl; // \gamma_qg^0 - qq0 = gammansplus0(N) * M_PIl; // \gamma_qq^0 + gg0 = gammagg0(N) * M_PIl; // \gamma_gg^0 + gq0 = gammagS0(N) * M_PIl; // \gamma_gq^0 + qg0 = gammaSg0(N) * M_PIl; // \gamma_qg^0 + qq0 = gammansplus0(N) * M_PIl; // \gamma_qq^0 // Define the Eigenvalues of the Singlet Matrix plus0 = 0.5 * @@ -124,8 +125,8 @@ void AnomDimensions::sums(std::complex N) { } // LO anomalous dimensions -std::complex -AnomDimensions::gammagg0(std::complex N) { +std::complex AnomDimensions::gammagg0( + std::complex N) { return (-1. / (4. * M_PIl) * (CA * (4. * (HAP.HS(1, N - 2.) - 2. * HAP.HS(1, N - 1.) - 2. * HAP.HS(1, N + 1.) + HAP.HS(1, N + 2.) + @@ -134,24 +135,24 @@ AnomDimensions::gammagg0(std::complex N) { 2. / 3. * NF)); } -std::complex -AnomDimensions::gammagS0(std::complex N) { +std::complex AnomDimensions::gammagS0( + std::complex N) { return -1. / (4. * M_PIl) * (2. * CF * (2. * HAP.HS(1, N - 2.) - 4. * HAP.HS(1, N - 1.) - HAP.HS(1, N + 1.) + 3. * HAP.HS(1, N))); } -std::complex -AnomDimensions::gammaSg0(std::complex N) { +std::complex AnomDimensions::gammaSg0( + std::complex N) { return -1. / (4. * M_PIl) * (2. * NF * (HAP.HS(1, N - 1.) + 4. * HAP.HS(1, N + 1.) - 2. * HAP.HS(1, N + 2) - 3. * HAP.HS(1, N))); } -std::complex -AnomDimensions::gammansplus0(std::complex N) { +std::complex AnomDimensions::gammansplus0( + std::complex N) { return -1. / (4. * M_PIl) * CF * (2. * (HAP.HS(1, N - 1.) + HAP.HS(1, N + 1.)) - 3.); } @@ -231,42 +232,42 @@ void AnomDimensions::DEF1(std::complex N) { } // NLO anomalous dimensions -std::complex -AnomDimensions::gammagg1(std::complex N) { +std::complex AnomDimensions::gammagg1( + std::complex N) { return (CA * CA * PGGA + 0.5 * NF * (CA * PGGB + CF * PGGC)) * 4. / pow(4 * M_PIl, 2); } -std::complex -AnomDimensions::gammagS1(std::complex N) { +std::complex AnomDimensions::gammagS1( + std::complex N) { return (CF * CF * PGQA + CF * CA * PGQB + 0.5 * NF * CF * PGQC) * 4. / pow(4 * M_PIl, 2); } -std::complex -AnomDimensions::gammaSg1(std::complex N) { +std::complex AnomDimensions::gammaSg1( + std::complex N) { return 0.5 * NF * (CA * PQGA + CF * PQGB) * 4. / pow(4 * M_PIl, 2); } -std::complex -AnomDimensions::gammaps1(std::complex N) { +std::complex AnomDimensions::gammaps1( + std::complex N) { return (0.5 * NF * CF * PPSA * 4. / pow(4 * M_PIl, 2)); } -std::complex -AnomDimensions::gammansplus1(std::complex N) { +std::complex AnomDimensions::gammansplus1( + std::complex N) { return CF * ((CF - CA / 2.) * PNPA + CA * PNSB + 0.5 * (NF)*PNSC) / pow(4 * M_PIl, 2); } -std::complex -AnomDimensions::gammansminus1(std::complex N) { +std::complex AnomDimensions::gammansminus1( + std::complex N) { return CF * ((CF - CA / 2.) * PNMA + CA * PNSB + 0.5 * NF * PNSC) / (pow(4 * M_PIl, 2)); } -std::complex -AnomDimensions::gammGGnlo(std::complex N) { +std::complex AnomDimensions::gammGGnlo( + std::complex N) { // A. Vogt https://arxiv.org/pdf/hep-ph/0404111.pdf std::complex GG1 = 4. * CA * NF * @@ -506,37 +507,37 @@ void AnomDimensions::DEF2(std::complex N) { } // NNLO anomalous dimensions -std::complex -AnomDimensions::gammagg2(std::complex N) { +std::complex AnomDimensions::gammagg2( + std::complex N) { return P2GGN / pow(4 * M_PIl, 3); } -std::complex -AnomDimensions::gammagS2(std::complex N) { +std::complex AnomDimensions::gammagS2( + std::complex N) { return P2GQN / pow(4 * M_PIl, 3); } -std::complex -AnomDimensions::gammaSg2(std::complex N) { +std::complex AnomDimensions::gammaSg2( + std::complex N) { return P2QGN / pow(4 * M_PIl, 3); } -std::complex -AnomDimensions::gammaps2(std::complex N) { +std::complex AnomDimensions::gammaps2( + std::complex N) { return P2PSN / pow(4 * M_PIl, 3); } -std::complex -AnomDimensions::gammansplus2(std::complex N) { +std::complex AnomDimensions::gammansplus2( + std::complex N) { return P2PLSN / pow(4 * M_PIl, 3); } -std::complex -AnomDimensions::gammansminus2(std::complex N) { +std::complex AnomDimensions::gammansminus2( + std::complex N) { return P2MINN / pow(4 * M_PIl, 3); } -std::complex -AnomDimensions::gammansval2(std::complex N) { +std::complex AnomDimensions::gammansval2( + std::complex N) { return P2VALN / pow(4 * M_PIl, 3); } diff --git a/src/CombinedRes.cpp b/src/CombinedRes.cpp index fc896fe..b66ba10 100644 --- a/src/CombinedRes.cpp +++ b/src/CombinedRes.cpp @@ -16,10 +16,11 @@ * ===================================================================================== */ -#include +#include "../include/CombinedRes.h" + #include -#include "../include/CombinedRes.h" +#include CombinedRes::CombinedRes(int order, int channel, std::string pdfname, void *params) { @@ -47,13 +48,13 @@ std::complex CombinedRes::Matching(std::complex N, long double pt, int scheme) { // TODO: re-check definition MH2 vs. Qs2 long double xp = std::pow(pt, 2) / MH2; - if (scheme == 0) // small-pt only + if (scheme == 0) // small-pt only { return (0.); - } else if (scheme == 1) // threshold only + } else if (scheme == 1) // threshold only { return (1.); - } else // combined + } else // combined { long double k = 2.; long double m = 9.75; @@ -62,9 +63,8 @@ std::complex CombinedRes::Matching(std::complex N, } } -std::complex -CombinedRes::CombinedResExpr(std::complex N, long double pt, - int scheme) { +std::complex CombinedRes::CombinedResExpr( + std::complex N, long double pt, int scheme) { double pp = static_cast(pt); // take only real part. Does not work for complex double nn = static_cast(N.real()); diff --git a/src/ComplexDefs.cpp b/src/ComplexDefs.cpp index bb70d81..2bd9635 100644 --- a/src/ComplexDefs.cpp +++ b/src/ComplexDefs.cpp @@ -197,9 +197,9 @@ bool isfinite(const std::complex &z) { } // Inner Product -std::complex -operator*(const std::vector> &c1, - const std::vector> &c2) { +std::complex operator*( + const std::vector> &c1, + const std::vector> &c2) { std::complex result; std::complex zero(0., 0.); result = std::inner_product(c1.begin(), c1.end(), c2.begin(), zero); @@ -244,8 +244,7 @@ std::complex log_c_angle(std::complex z, std::complex expm1(const std::complex &z) { const long double x = std::real(z), y = std::imag(z); - if ((std::abs(x) >= 1.0) || (std::abs(y) >= 1.0)) - return (std::exp(z) - 1.0); + if ((std::abs(x) >= 1.0) || (std::abs(y) >= 1.0)) return (std::exp(z) - 1.0); const long double expm1_x = std::expm1(x); const long double exp_x = 1.0 + expm1_x; @@ -261,8 +260,7 @@ std::complex log1p(const std::complex &z) { const long double xp1 = 1.0 + x, abs_x = std::abs(x), abs_y = std::abs(y); - if ((abs_x >= 1.0) || (abs_y >= 1.0)) - return std::log(1.0 + z); + if ((abs_x >= 1.0) || (abs_y >= 1.0)) return std::log(1.0 + z); const long double y_over_xp1 = y / xp1; @@ -295,8 +293,7 @@ std::complex LogGamma(std::complex z) { 0.8441822398385274E-4, -0.26190838401581408E-4, 0.3689918265953162E-5}; std::complex sum = c[0]; - for (int i = 1; i < 15; i++) - sum += c[i] / (zm1 + i); + for (int i = 1; i < 15; i++) sum += c[i] / (zm1 + i); const std::complex log_Gamma_z = log_sqrt_2Pi + std::log(sum) + z_m_0p5 * std::log(z_pg_m0p5) - @@ -341,8 +338,7 @@ std::complex Gamma_inv(const std::complex &z) { 0.844182239838527E-4, -0.26190838401581408E-4, 0.3689918265953162E-5}; std::complex sum = c[0]; - for (int i = 1; i < 15; i++) - sum += c[i] / (zm1 + i); + for (int i = 1; i < 15; i++) sum += c[i] / (zm1 + i); const std::complex Gamma_inv_z = exp(z_pg_m0p5 - z_m_0p5 * log(z_pg_m0p5) - log_sqrt_2Pi) / sum; @@ -491,9 +487,8 @@ std::complex CBesselK(std::complex nu, } // Useful functions to evaluate Hypergeometric2F1 (see AEAE for more details) -std::complex -Gamma_ratio_diff_small_eps(const std::complex &z, - const std::complex &eps) { +std::complex Gamma_ratio_diff_small_eps( + const std::complex &z, const std::complex &eps) { const long double g = 4.7421875; if (inf_norm(eps) > 0.1) std::cout << "One must have |eps|oo < 0.1 in Gamma_ratio_diff_small_eps.", @@ -560,9 +555,8 @@ Gamma_ratio_diff_small_eps(const std::complex &z, } } -std::complex -Gamma_inv_diff_eps(const std::complex &z, - const std::complex &eps) { +std::complex Gamma_inv_diff_eps( + const std::complex &z, const std::complex &eps) { const std::complex eps_pz = z + eps; const double x = std::real(z), eps_px = std::real(eps_pz); const int n = static_cast(std::rint(x)); @@ -590,28 +584,25 @@ Gamma_inv_diff_eps(const std::complex &z, } } else if (is_z_negative_integer && is_eps_pz_negative_integer) { long double fact = -1.0; - for (int k = -1; k >= n; k--) - fact *= k; + for (int k = -1; k >= n; k--) fact *= k; return fact; } else return Gamma_ratio_diff_small_eps(z, eps) * Gamma_inv(eps_pz); } -std::complex -A_sum_init(const int m, const std::complex &eps, - const std::complex &Gamma_inv_one_meps) { +std::complex A_sum_init( + const int m, const std::complex &eps, + const std::complex &Gamma_inv_one_meps) { const std::complex one_meps = 1.0 - eps; if (one_meps - m != 1 - m) { std::complex Gamma_inv_one_meps_mm = Gamma_inv_one_meps; - for (int i = 1; i <= m; i++) - Gamma_inv_one_meps_mm *= one_meps - i; + for (int i = 1; i <= m; i++) Gamma_inv_one_meps_mm *= one_meps - i; return Gamma_inv_one_meps_mm / eps; } else { long double fact = 1.0; - for (int n = 2; n < m; n++) - fact *= n; + for (int n = 2; n < m; n++) fact *= n; return (m % 2 == 0) ? (fact) : (-fact); } @@ -634,16 +625,15 @@ std::complex log_A_sum_init(const int m, } } -std::complex -B_sum_init_PS_one(const std::complex &a, - const std::complex &b, - const std::complex &c, - const std::complex &Gamma_c, - const std::complex &Gamma_inv_one_meps, - const std::complex &Gamma_inv_eps_pa_pm, - const std::complex &Gamma_inv_eps_pb_pm, - const std::complex &one_minus_z, const int m, - const std::complex &eps) { +std::complex B_sum_init_PS_one( + const std::complex &a, const std::complex &b, + const std::complex &c, + const std::complex &Gamma_c, + const std::complex &Gamma_inv_one_meps, + const std::complex &Gamma_inv_eps_pa_pm, + const std::complex &Gamma_inv_eps_pb_pm, + const std::complex &one_minus_z, const int m, + const std::complex &eps) { const long double inf_norm_eps = inf_norm(eps); const long double phase = (m % 2 == 0) ? (1) : (-1); const std::complex a_pm = a + m; @@ -653,15 +643,13 @@ B_sum_init_PS_one(const std::complex &a, const std::complex Pi_eps_pm = M_PIl * (eps + m); std::complex Gamma_inv_one_meps_mm = Gamma_inv_one_meps; - for (int i = 1; i <= m; i++) - Gamma_inv_one_meps_mm *= one_meps - i; + for (int i = 1; i <= m; i++) Gamma_inv_one_meps_mm *= one_meps - i; if (inf_norm_eps > 0.1) { const std::complex Gamma_inv_eps_pm_p1 = phase * std::sin(Pi_eps) / (Pi_eps_pm * Gamma_inv_one_meps_mm); std::complex prod1 = Gamma_inv_one_meps * Gamma_inv_eps_pa_pm * Gamma_inv_eps_pb_pm; - for (int n = 0; n < m; n++) - prod1 *= (a + n) * (b + n) / (n + 1.0); + for (int n = 0; n < m; n++) prod1 *= (a + n) * (b + n) / (n + 1.0); const std::complex prod2 = Gamma_inv(a) * Gamma_inv(b) * Gamma_inv_eps_pm_p1 * std::pow(one_minus_z, eps); @@ -704,8 +692,7 @@ B_sum_init_PS_one(const std::complex &a, phase * sin(Pi_eps) / (Pi_eps_pm * Gamma_inv_one_meps_mm); std::complex prod1 = Gamma_inv_one_meps * Gamma_inv_eps_pa_pm * Gamma_inv_eps_pb_pm; - for (int n = 0; n < m; n++) - prod1 *= (a + n) * (b + n) / (n + 1.0); + for (int n = 0; n < m; n++) prod1 *= (a + n) * (b + n) / (n + 1.0); const std::complex prod2 = Gamma_inv(a) * Gamma_inv(b) * Gamma_inv_eps_pm_p1 * pow(one_minus_z, eps); @@ -716,15 +703,14 @@ B_sum_init_PS_one(const std::complex &a, } } -std::complex -B_sum_init_PS_infinity(const std::complex &a, - const std::complex &c, - const std::complex &Gamma_c, - const std::complex &Gamma_inv_cma, - const std::complex &Gamma_inv_one_meps, - const std::complex &Gamma_inv_eps_pa_pm, - const std::complex &z, const int m, - const std::complex &eps) { +std::complex B_sum_init_PS_infinity( + const std::complex &a, const std::complex &c, + const std::complex &Gamma_c, + const std::complex &Gamma_inv_cma, + const std::complex &Gamma_inv_one_meps, + const std::complex &Gamma_inv_eps_pa_pm, + const std::complex &z, const int m, + const std::complex &eps) { const double inf_norm_eps = inf_norm(eps); const double phase = (m % 2 == 0) ? (1) : (-1); const std::complex cma = c - a, a_mc_p1 = 1.0 - c + a; @@ -738,8 +724,7 @@ B_sum_init_PS_infinity(const std::complex &a, const std::complex Pi_eps_pm = M_PIl * (eps + m); std::complex Gamma_inv_one_meps_mm = Gamma_inv_one_meps; - for (int i = 1; i <= m; i++) - Gamma_inv_one_meps_mm *= one_meps - i; + for (int i = 1; i <= m; i++) Gamma_inv_one_meps_mm *= one_meps - i; if (inf_norm_eps > 0.1) { const std::complex Gamma_inv_eps_pm_p1 = @@ -781,8 +766,7 @@ B_sum_init_PS_infinity(const std::complex &a, Gamma_inv_mp1 /= n + 1.0; if (n != n0) { - if (is_n0_here) - prod_eps_pa_mc_p1_n0 *= eps_pa_mc_p1_pn; + if (is_n0_here) prod_eps_pa_mc_p1_n0 *= eps_pa_mc_p1_pn; sum += (is_eps_non_zero) ? (log1p(eps / a_mc_p1_pn)) : (1.0 / a_mc_p1_pn); } @@ -926,11 +910,10 @@ std::complex hyp_PS_zero(const std::complex &a, return sum; } -std::complex -hyp_PS_one(const std::complex &a, - const std::complex &b, - const std::complex &c, - const std::complex &one_minus_z) { +std::complex hyp_PS_one( + const std::complex &a, const std::complex &b, + const std::complex &c, + const std::complex &one_minus_z) { const int m = static_cast(std::rint(std::real(c - a - b))); const int phase = (m % 2 == 0) ? (1) : (-1); const int m_m1 = m - 1, m_p1 = m + 1; @@ -975,8 +958,7 @@ hyp_PS_one(const std::complex &a, A_sum += A_term, prod_B *= ratio; } - if (m > 0) - prod_B *= (a + m - 1.0) * (b + m - 1.0) / m; + if (m > 0) prod_B *= (a + m - 1.0) * (b + m - 1.0) / m; std::complex B_extra_term = prod_B * Gamma_prod * Gamma_inv_one_meps; @@ -1070,8 +1052,7 @@ std::complex hyp_PS_infinity(const std::complex &a, prod_B *= ratio; } - if (m > 0) - prod_B *= (a + m - 1.0) * (a_mc_p1 + m - 1.0) / m; + if (m > 0) prod_B *= (a + m - 1.0) * (a_mc_p1 + m - 1.0) / m; std::complex B_extra_term = prod_B * Gamma_prod * Gamma_inv_one_meps; @@ -1218,8 +1199,7 @@ std::complex Hyp2F1(std::complex a, } const long double abs_zm1 = std::abs(zm1); - if (abs_zm1 < 1E-5) - return hyp_PS_one(a, b, c, -zm1); + if (abs_zm1 < 1E-5) return hyp_PS_one(a, b, c, -zm1); const long double abs_z = abs(z), abs_z_inv = 1.0 / abs_z, abs_z_over_zm1 = abs_z / abs_zm1; @@ -1235,8 +1215,7 @@ std::complex Hyp2F1(std::complex a, for (unsigned int i = 0; i < 5; i++) { const long double R = R_tab[i]; - if (abs_z <= R) - return hyp_PS_zero(a, b, c, z); + if (abs_z <= R) return hyp_PS_zero(a, b, c, z); if (is_cmb_small && (abs_z_over_zm1 <= R)) return std::pow(-zm1, -a) * hyp_PS_zero(a, c - b, c, z / zm1); } @@ -1244,13 +1223,11 @@ std::complex Hyp2F1(std::complex a, for (unsigned int i = 0; i < 5; i++) { const long double R = R_tab[i]; - if (abs_z_inv <= R) - return hyp_PS_infinity(a, b, c, z); + if (abs_z_inv <= R) return hyp_PS_infinity(a, b, c, z); if (is_cmb_small && (abs_zm1_over_z <= R)) return std::pow(-zm1, -a) * hyp_PS_infinity(a, c - b, c, z / zm1); - if (are_abc_small && (abs_zm1 <= R)) - return hyp_PS_one(a, b, c, -zm1); + if (are_abc_small && (abs_zm1 <= R)) return hyp_PS_one(a, b, c, -zm1); if (are_a_cmb_c_small && (abs_zm1_inv <= R)) return std::pow(-zm1, -a) * hyp_PS_one(a, c - b, c, -1.0 / zm1); } diff --git a/src/HarmonicSum.cpp b/src/HarmonicSum.cpp index 713c44c..b5b38fc 100644 --- a/src/HarmonicSum.cpp +++ b/src/HarmonicSum.cpp @@ -2,9 +2,9 @@ // Table of interpolated function HSum::HSum(bool verbose, bool testinterfun, bool testharmsums) - : _verbose(verbose), _testinterpolatedfunction(testinterfun), + : _verbose(verbose), + _testinterpolatedfunction(testinterfun), _testharmonicsums(testharmsums) { - (_verbose) ? std::cout << "HSUM is a class to evaluate Harmonic Sums up to four order" @@ -40,241 +40,241 @@ HSum::HSum(bool verbose, bool testinterfun, bool testharmsums) std::complex HSum::HS(int i, std::complex N) { switch (i) { - case (-4): - return H_m4(N); - case (-3): - return H_m3(N); - case (-2): - return H_m2(N); - case (-1): - return H_m1(N); - case (1): - return H_1(N); - case (2): - return H_2(N); - case (3): - return H_3(N); - case (4): - return H_4(N); - } - cout << "ERROR: not valid Harmonic Sums index (Implemented Weight up to four " - "yet)" - << endl; - return (0.); -} - -std::complex HSum::HS(int i, int j, std::complex N) { - switch (i) { - case (-3): { - switch (j) { - case (-1): - return H_m3_m1(N); - case (1): - return H_m3_1(N); - } - } break; - case (-2): { - switch (j) { - case (-2): - return H_m2_m2(N); - case (-1): - return H_m2_m1(N); - case (1): - return H_m2_1(N); - case (2): - return H_m2_2(N); - } - } break; - case (-1): { - switch (j) { - case (-3): - return H_m1_m3(N); - case (-2): - return H_m1_m2(N); - case (-1): - return H_m1_m1(N); - case (1): - return H_m1_1(N); - case (2): - return H_m1_2(N); - case (3): - return H_m1_3(N); - } - } break; - case (1): { - switch (j) { + case (-4): + return H_m4(N); case (-3): - return H_1_m3(N); + return H_m3(N); case (-2): - return H_1_m2(N); + return H_m2(N); case (-1): - return H_1_m1(N); + return H_m1(N); case (1): - return H_1_1(N); + return H_1(N); case (2): - return H_1_2(N); + return H_2(N); case (3): - return H_1_3(N); - } - } break; - case (2): { - switch (j) { - case (-2): - return H_2_m2(N); - case (-1): - return H_2_m1(N); - case (1): - return H_2_1(N); - case (2): - return H_2_2(N); - } - } break; - case (3): { - switch (j) { - case (-1): - return H_3_m1(N); - case (1): - return H_3_1(N); - } - } break; + return H_3(N); + case (4): + return H_4(N); } - cout << "ERROR: not valid Harmonic Sums index (Implemented Weight up to 4th " - "order.)" + cout << "ERROR: not valid Harmonic Sums index (Implemented Weight up to four " + "yet)" << endl; return (0.); } -std::complex HSum::HS(int i, int j, int k, - std::complex N) { +std::complex HSum::HS(int i, int j, std::complex N) { switch (i) { - case (-2): { - switch (j) { - case (-1): { - switch (k) { - case (-1): - return H_m2_m1_m1(N); - case (1): - return H_m2_m1_1(N); - } - } break; - case (1): { - switch (k) { - case (-1): - return H_m2_1_m1(N); - case (1): - return H_m2_1_1(N); + case (-3): { + switch (j) { + case (-1): + return H_m3_m1(N); + case (1): + return H_m3_1(N); } } break; - } - } break; - case (-1): { - switch (j) { case (-2): { - switch (k) { - case (-1): - return H_m1_m2_m1(N); - case (1): - return H_m1_m2_1(N); + switch (j) { + case (-2): + return H_m2_m2(N); + case (-1): + return H_m2_m1(N); + case (1): + return H_m2_1(N); + case (2): + return H_m2_2(N); } } break; case (-1): { - switch (k) { - case (-2): - return H_m1_m1_m2(N); - case (-1): - return H_m1_m1_m1(N); - case (1): - return H_m1_m1_1(N); - case (2): - return H_m1_m1_2(N); + switch (j) { + case (-3): + return H_m1_m3(N); + case (-2): + return H_m1_m2(N); + case (-1): + return H_m1_m1(N); + case (1): + return H_m1_1(N); + case (2): + return H_m1_2(N); + case (3): + return H_m1_3(N); } } break; case (1): { - switch (k) { - case (-2): - return H_m1_1_m2(N); - case (-1): - return H_m1_1_m1(N); - case (1): - return H_m1_1_1(N); - case (2): - return H_m1_1_2(N); + switch (j) { + case (-3): + return H_1_m3(N); + case (-2): + return H_1_m2(N); + case (-1): + return H_1_m1(N); + case (1): + return H_1_1(N); + case (2): + return H_1_2(N); + case (3): + return H_1_3(N); } } break; case (2): { - switch (k) { - case (-1): - return H_m1_2_m1(N); - case (1): - return H_m1_2_1(N); + switch (j) { + case (-2): + return H_2_m2(N); + case (-1): + return H_2_m1(N); + case (1): + return H_2_1(N); + case (2): + return H_2_2(N); } } break; - } - } break; - case (1): { - switch (j) { + case (3): { + switch (j) { + case (-1): + return H_3_m1(N); + case (1): + return H_3_1(N); + } + } break; + } + cout << "ERROR: not valid Harmonic Sums index (Implemented Weight up to 4th " + "order.)" + << endl; + return (0.); +} + +std::complex HSum::HS(int i, int j, int k, + std::complex N) { + switch (i) { case (-2): { - switch (k) { - case (-1): - return H_1_m2_m1(N); - case (1): - return H_1_m2_1(N); + switch (j) { + case (-1): { + switch (k) { + case (-1): + return H_m2_m1_m1(N); + case (1): + return H_m2_m1_1(N); + } + } break; + case (1): { + switch (k) { + case (-1): + return H_m2_1_m1(N); + case (1): + return H_m2_1_1(N); + } + } break; } } break; case (-1): { - switch (k) { - case (-2): - return H_1_m1_m2(N); - case (-1): - return H_1_m1_m1(N); - case (1): - return H_1_m1_1(N); - case (2): - return H_1_m1_2(N); + switch (j) { + case (-2): { + switch (k) { + case (-1): + return H_m1_m2_m1(N); + case (1): + return H_m1_m2_1(N); + } + } break; + case (-1): { + switch (k) { + case (-2): + return H_m1_m1_m2(N); + case (-1): + return H_m1_m1_m1(N); + case (1): + return H_m1_m1_1(N); + case (2): + return H_m1_m1_2(N); + } + } break; + case (1): { + switch (k) { + case (-2): + return H_m1_1_m2(N); + case (-1): + return H_m1_1_m1(N); + case (1): + return H_m1_1_1(N); + case (2): + return H_m1_1_2(N); + } + } break; + case (2): { + switch (k) { + case (-1): + return H_m1_2_m1(N); + case (1): + return H_m1_2_1(N); + } + } break; } } break; case (1): { - switch (k) { - case (-2): - return H_1_1_m2(N); - case (-1): - return H_1_1_m1(N); - case (1): - return H_1_1_1(N); - case (2): - return H_1_1_2(N); + switch (j) { + case (-2): { + switch (k) { + case (-1): + return H_1_m2_m1(N); + case (1): + return H_1_m2_1(N); + } + } break; + case (-1): { + switch (k) { + case (-2): + return H_1_m1_m2(N); + case (-1): + return H_1_m1_m1(N); + case (1): + return H_1_m1_1(N); + case (2): + return H_1_m1_2(N); + } + } break; + case (1): { + switch (k) { + case (-2): + return H_1_1_m2(N); + case (-1): + return H_1_1_m1(N); + case (1): + return H_1_1_1(N); + case (2): + return H_1_1_2(N); + } + } break; + case (2): { + switch (k) { + case (-1): + return H_1_2_m1(N); + case (1): + return H_1_2_1(N); + } + } break; } } break; case (2): { - switch (k) { - case (-1): - return H_1_2_m1(N); - case (1): - return H_1_2_1(N); - } - } break; - } - } break; - case (2): { - switch (j) { - case (-1): { - switch (k) { - case (-1): - return H_2_m1_m1(N); - case (1): - return H_2_m1_1(N); - } - } break; - case (1): { - switch (k) { - case (-1): - return H_2_1_m1(N); - case (1): - return H_2_1_1(N); + switch (j) { + case (-1): { + switch (k) { + case (-1): + return H_2_m1_m1(N); + case (1): + return H_2_m1_1(N); + } + } break; + case (1): { + switch (k) { + case (-1): + return H_2_1_m1(N); + case (1): + return H_2_1_1(N); + } + } break; } } break; - } - } break; } cout << "ERROR: not valid Harmonic Sums index (Implemented Weight up to four " "yet)" @@ -285,94 +285,94 @@ std::complex HSum::HS(int i, int j, int k, std::complex HSum::HS(int i, int j, int k, int m, std::complex N) { switch (i) { - case (-1): { - switch (j) { case (-1): { - switch (k) { - case (-1): { - switch (m) { - case (-1): - return H_m1_m1_m1_m1(N); - case (1): - return H_m1_m1_m1_1(N); - } - } break; - case (1): { - switch (m) { - case (-1): - return H_m1_m1_1_m1(N); - case (1): - return H_m1_m1_1_1(N); - } - } break; - } - } break; - case (1): { - switch (k) { - case (-1): { - switch (m) { - case (-1): - return H_m1_1_m1_m1(N); - case (1): - return H_m1_1_m1_1(N); - } - } break; - case (1): { - switch (m) { - case (-1): - return H_m1_1_1_m1(N); - case (1): - return H_m1_1_1_1(N); - } - } break; - } - } break; - } - } break; - case (1): { - switch (j) { - case (-1): { - switch (k) { - case (-1): { - switch (m) { - case (-1): - return H_1_m1_m1_m1(N); - case (1): - return H_1_m1_m1_1(N); - } - } break; - case (1): { - switch (m) { - case (-1): - return H_1_m1_1_m1(N); - case (1): - return H_1_m1_1_1(N); - } - } break; + switch (j) { + case (-1): { + switch (k) { + case (-1): { + switch (m) { + case (-1): + return H_m1_m1_m1_m1(N); + case (1): + return H_m1_m1_m1_1(N); + } + } break; + case (1): { + switch (m) { + case (-1): + return H_m1_m1_1_m1(N); + case (1): + return H_m1_m1_1_1(N); + } + } break; + } + } break; + case (1): { + switch (k) { + case (-1): { + switch (m) { + case (-1): + return H_m1_1_m1_m1(N); + case (1): + return H_m1_1_m1_1(N); + } + } break; + case (1): { + switch (m) { + case (-1): + return H_m1_1_1_m1(N); + case (1): + return H_m1_1_1_1(N); + } + } break; + } + } break; } } break; case (1): { - switch (k) { - case (-1): { - switch (m) { - case (-1): - return H_1_1_m1_m1(N); - case (1): - return H_1_1_m1_1(N); - } - } break; - case (1): { - switch (m) { - case (-1): - return H_1_1_1_m1(N); - case (1): - return H_1_1_1_1(N); - } - } break; + switch (j) { + case (-1): { + switch (k) { + case (-1): { + switch (m) { + case (-1): + return H_1_m1_m1_m1(N); + case (1): + return H_1_m1_m1_1(N); + } + } break; + case (1): { + switch (m) { + case (-1): + return H_1_m1_1_m1(N); + case (1): + return H_1_m1_1_1(N); + } + } break; + } + } break; + case (1): { + switch (k) { + case (-1): { + switch (m) { + case (-1): + return H_1_1_m1_m1(N); + case (1): + return H_1_1_m1_1(N); + } + } break; + case (1): { + switch (m) { + case (-1): + return H_1_1_1_m1(N); + case (1): + return H_1_1_1_1(N); + } + } break; + } + } break; } } break; - } - } break; } cout << "ERROR: not valid Harmonic Sums index (Implemented Weight up to four " "yet)" diff --git a/src/MellinFunc.cpp b/src/MellinFunc.cpp index 5a5aaf3..891ebd8 100644 --- a/src/MellinFunc.cpp +++ b/src/MellinFunc.cpp @@ -59,8 +59,8 @@ std::complex MellinFunc::Li3mzplus(std::complex N) { 0.5 * zeta2 * H.HS(-2, N - 1.) - 3. / 4. * zeta3 * H.HS(-1, N - 1.) + 1. / 8. * zeta2q - 3. / 4. * zeta3 * log2)); } -std::complex -MellinFunc::Li2zLogminus(std::complex N) { +std::complex MellinFunc::Li2zLogminus( + std::complex N) { return (1. / N * (-2. * zeta3 - zeta2 * H.HS(1, N) + 1. / N * (H.HS(1, N) * H.HS(1, N) + H.HS(2, N)) + H.HS(2, 1, N))); @@ -77,8 +77,8 @@ std::complex MellinFunc::Li2zLogplus(std::complex N) { (2. * N * H.HS(-2, 1, N) - 2. * N * zeta2 * (H.HS(-1, N) + log2) + (zeta2 + 2. * H.HS(-1, 1, N) - log2q) + 5. / 4. * N * zeta3))); } -std::complex -MellinFunc::Li2mzLogplus(std::complex N) { +std::complex MellinFunc::Li2mzLogplus( + std::complex N) { long double n = 0.; return (1. / (N * N) * (log2 * (-0.5 * zeta2 * N + log2) + @@ -95,8 +95,8 @@ std::complex MellinFunc::Li2mzLogz(std::complex N) { pow(-1., n) * (1. / (N * N) * (0.5 * zeta2 + H.HS(-2, N)) + 1. / (N * N * N) * (2. * H.HS(-1, N) + 2. * log2))); } -std::complex -MellinFunc::Li2zLogzplusminus(std::complex N) { +std::complex MellinFunc::Li2zLogzplusminus( + std::complex N) { long double n = 0.; return ( 1. / 960. * @@ -107,8 +107,8 @@ MellinFunc::Li2zLogzplusminus(std::complex N) { 60. * (H.HS(2, N - 1.) * H.HS(2, N - 1.) + H.HS(4, N - 1.)) + 240. * H.HS(3, 1, N - 1.)))); } -std::complex -MellinFunc::Li2mzLogzplus(std::complex N) { +std::complex MellinFunc::Li2mzLogzplus( + std::complex N) { long double n = 0.; return (-pow(-1., n) * (2. * H.HS(3, -1, N - 1.) + H.HS(2, -2, N - 1.) - @@ -117,15 +117,15 @@ MellinFunc::Li2mzLogzplus(std::complex N) { 4. * Li4 + 13. / 8. * zeta2q - 7. / 2. * zeta3 * log2 + zeta2 * log2q - 1. / 6. * log2q * log2q)); } -std::complex -MellinFunc::LogzLogminus2minus(std::complex N) { +std::complex MellinFunc::LogzLogminus2minus( + std::complex N) { return (H.HS(2, N - 1.) * H.HS(2, N - 1.) - zeta2 * H.HS(2, N - 1.) + 2. * H.HS(4, N - 1.) + H.HS(1, N - 1.) * H.HS(1, N - 1.) * (H.HS(2, N - 1.) - zeta2) + 2. * H.HS(1, N - 1.) * (H.HS(3, N - 1.) - zeta3) - 4. / 5. * zeta2q); } -std::complex -MellinFunc::Li2mzLogplusplus(std::complex N) { +std::complex MellinFunc::Li2mzLogplusplus( + std::complex N) { long double n = 0.; return (-pow(-1., n) * (H.HS(1, 2, -1, N - 1.) + 2. * H.HS(2, 1, -1, N - 1.) + @@ -139,15 +139,15 @@ MellinFunc::Li2mzLogplusplus(std::complex N) { 6. / 5. * zeta2q - 21. / 8. * zeta3 * log2 + 0.5 * zeta2 * log2q - 1. / 8. * log2q * log2q)); } -std::complex -MellinFunc::Logz2Logminusminus(std::complex N) { +std::complex MellinFunc::Logz2Logminusminus( + std::complex N) { return (-2. * zeta3 * H.HS(1, N - 1.) + 2. * H.HS(1, N - 1.) * H.HS(3, N - 1.) - 2. * zeta2 * H.HS(2, N - 1.) + H.HS(2, N - 1.) * H.HS(2, N - 1.) + 3. * H.HS(4, N - 1.) - zeta2q / 5.); } -std::complex -MellinFunc::Logz3minusplus(std::complex N) { +std::complex MellinFunc::Logz3minusplus( + std::complex N) { long double n = 0.; return (pow(-1., n) * (21. / 20. * zeta2q + 3. * H.HS(-4, N - 1.)) - 6. / 5. * zeta2q + 3. * H.HS(4, N - 1.)); @@ -157,8 +157,8 @@ std::complex MellinFunc::Logplusplus(std::complex N) { return (pow(-1., n) * (H.HS(1, -1, N - 1.) - 0.5 * log2q + (-H.HS(-1, N - 1.) + H.HS(1, N - 1.)) * log2)); } -std::complex -MellinFunc::Li2zLogplusplus(std::complex N) { +std::complex MellinFunc::Li2zLogplusplus( + std::complex N) { long double n = 0.; return (pow(-1., n) * (3. / 40. * zeta2q - H.HS(-2, -1, -1, N - 1.) - @@ -171,8 +171,8 @@ MellinFunc::Li2zLogplusplus(std::complex N) { 4. * (H.HS(-1, N - 1.) - H.HS(1, N - 1.)) * log2) - 5. / 8. * H.HS(1, N - 1.) * zeta3)); } -std::complex -MellinFunc::Logz2Logplusplus(std::complex N) { +std::complex MellinFunc::Logz2Logplusplus( + std::complex N) { long double n = 0.; return (pow(-1., n) * (2. * H.HS(3, -1, N - 1.) + 2. * H.HS(1, -3, N - 1.) + @@ -187,16 +187,16 @@ std::complex MellinFunc::LogzLogplus(std::complex N) { (-12. * log2 + pow(-1., n) * (6. * N * zeta2 + 12. * N * H.HS(-2, N) + 12. * H.HS(-1, N) + 12. * log2))); } -std::complex -MellinFunc::Logz2Logplus(std::complex N) { +std::complex MellinFunc::Logz2Logplus( + std::complex N) { long double n = 0.; return (1. / (N * N * N) * 2. * log2 - pow(-1., n) * (1. / (N * N) * (zeta2 + 2. * H.HS(-2, N)) + 1. / (N * N * N) * (2. * H.HS(-1, N) + 2. * log2) + 1. / N * (2. * H.HS(-3, N) + 3. / 2. * zeta3))); } -std::complex -MellinFunc::LogzLogplus2(std::complex N) { +std::complex MellinFunc::LogzLogplus2( + std::complex N) { long double n = 0.; return (-1. / (12. * N * N) * (12. * log2q + @@ -208,8 +208,8 @@ MellinFunc::LogzLogplus2(std::complex N) { 2. * N * H.HS(2, N) + log2) - N * zeta3)))); } -std::complex -MellinFunc::LogzLogplus2plus(std::complex N) { +std::complex MellinFunc::LogzLogplus2plus( + std::complex N) { long double n = 0.; return (2. * pow(-1., n) * (H.HS(1, 1, -2, N - 1.) + H.HS(1, 2, -1, N - 1.) + @@ -283,19 +283,19 @@ std::complex MellinFunc::Logz2(std::complex N) { std::complex MellinFunc::Logz3(std::complex N) { return (-6. / (N * N * N * N)); } -std::complex -MellinFunc::LogzLogminus(std::complex N) { +std::complex MellinFunc::LogzLogminus( + std::complex N) { return (-1. / N * zeta2 + 1. / N * H.HS(2, N) + 1. / (N * N) * H.HS(1, N)); } -std::complex -MellinFunc::LogzLogminus2(std::complex N) { +std::complex MellinFunc::LogzLogminus2( + std::complex N) { return (1. / (N * N) * (-H.HS(1, N) * H.HS(1, N) + 2. * N * H.HS(1, N) * (zeta2 - H.HS(2, N)) - H.HS(2, N) - 2. * N * H.HS(3, N) + 2. * N * zeta3)); } -std::complex -MellinFunc::Logz2Logminus(std::complex N) { +std::complex MellinFunc::Logz2Logminus( + std::complex N) { return (-2. / (N * N * N) * H.HS(1, N) + 1. / (N * N) * (2. * zeta2 - 2. * H.HS(2, N)) - 2. / N * (H.HS(3, N) - zeta3)); @@ -311,16 +311,16 @@ std::complex MellinFunc::Logminus3(std::complex N) { (H.HS(1, N) * H.HS(1, N) * H.HS(1, N) + 3. * H.HS(1, N) * H.HS(2, N) + 2. * H.HS(3, N))); } -std::complex -MellinFunc::LogzLogminusminus(std::complex N) { +std::complex MellinFunc::LogzLogminusminus( + std::complex N) { return (zeta3 + zeta2 * H.HS(1, N - 1.) - H.HS(1, N - 1.) * H.HS(2, N - 1.) - H.HS(3, N - 1.)); } std::complex MellinFunc::Logzminus(std::complex N) { return (-zeta2 + H.HS(2, N - 1.)); } -std::complex -MellinFunc::Logzminusplus(std::complex N) { +std::complex MellinFunc::Logzminusplus( + std::complex N) { long double n = 0.; return (pow(-1., n) * (0.25 * zeta2 + 0.5 * H.HS(-2, N - 1.)) + 0.5 * H.HS(2, N - 1.) - 0.5 * zeta2); @@ -329,30 +329,30 @@ std::complex MellinFunc::plus(std::complex N) { long double n = 0.; return (-pow(-1., n) * (H.HS(-1, N - 1.) + log2)); } -std::complex -MellinFunc::Li2minusminus(std::complex N) { +std::complex MellinFunc::Li2minusminus( + std::complex N) { return (-zeta2 * H.HS(1, N - 1.) + H.HS(1, 2, N - 1.) + zeta3); } -std::complex -MellinFunc::S12zregminus(std::complex N) { +std::complex MellinFunc::S12zregminus( + std::complex N) { return (-6. / 5. * zeta2q + H.HS(2, 1, 1, N - 1.)); } -std::complex -MellinFunc::Li2zregminus(std::complex N) { +std::complex MellinFunc::Li2zregminus( + std::complex N) { return (H.HS(2, 1, N - 1.) - 2. * zeta3); } -std::complex -MellinFunc::Li3zregminus(std::complex N) { +std::complex MellinFunc::Li3zregminus( + std::complex N) { return (-0.5 * zeta2q + zeta2 * H.HS(2, N - 1.) - H.HS(3, 1, N - 1.)); } -std::complex -MellinFunc::Li2zregLogminusminus(std::complex N) { +std::complex MellinFunc::Li2zregLogminusminus( + std::complex N) { return (6. / 5. * zeta2q - H.HS(1, 2, 1, N - 1.) - 2. * H.HS(2, 1, 1, N - 1.) + 2. * H.HS(1, N - 1.) * zeta3); } -std::complex -MellinFunc::Logplus3plus(std::complex N) { +std::complex MellinFunc::Logplus3plus( + std::complex N) { long double n = 0.; return ( 1. / 4. * pow(-1., n) * @@ -365,8 +365,8 @@ MellinFunc::Logplus3plus(std::complex N) { 12. * H.HS(1, -1, N - 1.) - 4. * H.HS(1, N - 1.) * log2 + log2q + 4. * H.HS(-1, N - 1.) * log2)))); } -std::complex -MellinFunc::Li3zregminusplus(std::complex N) { +std::complex MellinFunc::Li3zregminusplus( + std::complex N) { long double n = 0.; return (1. / 720. * (-5. * (pow(M_PIl, 4) - 12. * pow(M_PIl, 2) * H.HS(2, N - 1.) + @@ -376,8 +376,8 @@ MellinFunc::Li3zregminusplus(std::complex N) { 60. * H.HS(-3, 1, N - 1.) + 5. * M_PIl * M_PIl * log2q - 5. * (log2q * log2q + 24. * Li4 + 21. * log2 * zeta3)))); } -std::complex -MellinFunc::Li3zoverplusplus(std::complex N) { +std::complex MellinFunc::Li3zoverplusplus( + std::complex N) { long double n = 0.; return (1. / 288. * pow(-1., n) * (pow(M_PIl, 4) + 288. * H.HS(3, -1, N - 1.) - diff --git a/src/SmallptExp.cpp b/src/SmallptExp.cpp index d9a6c4e..ab41a9e 100644 --- a/src/SmallptExp.cpp +++ b/src/SmallptExp.cpp @@ -16,11 +16,12 @@ * ===================================================================================== */ -#include -#include +#include "../include/SmallptExp.h" + #include -#include "../include/SmallptExp.h" +#include +#include SmallptExp::SmallptExp(int order, int channel, void *params) : AD(params) { PhysParams param = *reinterpret_cast(params); @@ -45,7 +46,7 @@ SmallptExp::SmallptExp(int order, int channel, void *params) : AD(params) { LR = std::log(MH2 / MUR2); aass = static_cast( - param.alphas / M_PIl); // TODO; long double-check PI normalization + param.alphas / M_PIl); // TODO; long double-check PI normalization SROOT = static_cast(param.sroot); SIGMA0 = static_cast(param.sigma0); @@ -112,15 +113,15 @@ std::complex SmallptExp::C2GG(std::complex N) { // Sigma pt-enhanced functions // //------------------------------------------------------------------------------------------// -std::complex -SmallptExp::SmallptExpExpr(std::complex N, long double pt) { +std::complex SmallptExp::SmallptExpExpr( + std::complex N, long double pt) { // Init. Anomalous Dimensions; // Notice that in order to compare the AD here and // in the notebooks, the AD here have to be shifted by // (+1) i.e computed at (N+1). - AD.ComputeGamma(N + 1., 1); // Init. Anomalous Dimensions + AD.ComputeGamma(N + 1., 1); // Init. Anomalous Dimensions - int pc = 2; // TODO: long double-check this + int pc = 2; // TODO: long double-check this std::complex zero(0., 0.); std::complex ones(1., 0.); std::complex result(0., 0.); @@ -132,81 +133,88 @@ SmallptExp::SmallptExpExpr(std::complex N, long double pt) { long double xp = std::pow(pt, 2) / MH2; switch (ORD) { - case (0): // order as^1 - { - if ((CHANNEL == 0) || (CHANNEL == 4)) // gg-channel or ALL + case (0): // order as^1 { - // pt-enhanced Sigma terms - std::complex Sigma12gg = -Apt1g / 2.; - std::complex Sigma11gg = -Bpt1g - 2 * AD.gg0 - Apt1g * LQ; - - // constant-terms when pt->0 - std::complex HH1GG = - h1gg + 2. * C1GG(N) + 2. * (LF - LQ) * AD.gg0 - - LQ * (Bpt1g + (Bpt1g * LQ) / 2.) - pc * LR * Beta0; - - result += aass * (Sigma12gg * LC2(xp) + Sigma11gg * LC1(xp) + HH1GG); - } - if ((CHANNEL == 1) || (CHANNEL == 4)) { - result += zero; - } - if ((CHANNEL == 2) || (CHANNEL == 4)) { - result += zero; - } - } break; - case (1): // order as^2 - { - if ((CHANNEL == 0) || (CHANNEL == 4)) // gg-channel or ALL + if ((CHANNEL == 0) || (CHANNEL == 4)) // gg-channel or ALL + { + // pt-enhanced Sigma terms + std::complex Sigma12gg = -Apt1g / 2.; + std::complex Sigma11gg = -Bpt1g - 2 * AD.gg0 - Apt1g * LQ; + + // constant-terms when pt->0 + std::complex HH1GG = + h1gg + 2. * C1GG(N) + 2. * (LF - LQ) * AD.gg0 - + LQ * (Bpt1g + (Bpt1g * LQ) / 2.) - pc * LR * Beta0; + + result += aass * (Sigma12gg * LC2(xp) + Sigma11gg * LC1(xp) + HH1GG); + } + if ((CHANNEL == 1) || (CHANNEL == 4)) { + result += zero; + } + if ((CHANNEL == 2) || (CHANNEL == 4)) { + result += zero; + } + } break; + case (1): // order as^2 { - // pt-enhanced Sigma terms - std::complex Sigma24gg = std::pow(Apt1g, 2) / 8.; - std::complex Sigma23gg = - -Apt1g * (Beta0 / 3. + 1. / 2. * (-Bpt1g - Apt1g * LQ - 2. * AD.gg0)); - std::complex Sigma22gg = - (-Apt2g - - (-Bpt1g - 2. * AD.gg0 - Apt1g * LQ) * (Bpt1g - Beta0 + Apt1g * LQ)) / - 2. - - (Apt1g * (2 * C1GG(N) + h1gg + 2. * AD.gg0 * (LF - LQ) - - LQ * (Bpt1g + (Bpt1g * LQ) / 2.) - (Beta0 * (-LQ + LR)) - - Beta0 * LR * pc)) / - 2. + - (-2. * AD.gg0 * (-Bpt1g - 2. * AD.gg0 - Apt1g * LQ) + - 2. * AD.gq0 * AD.qg0) / - 2.; - std::complex Sigma21gg = - -Bpt2g + 2. * Beta0 * C1GG(N) - 2. * AD.gg1 - Apt2g * LQ + - Beta0 * (-Bpt1g - AD.gg0 - Apt1g * LQ) * (LQ - LR) - - (Bpt1g + 2. * AD.gg0 + Apt1g * LQ) * - (2. * C1GG(N) + h1gg + 2. * AD.gg0 * (LF - LQ) - - LQ * (Bpt1g + (Bpt1g * LQ) / 2.) - Beta0 * LR * pc) - - 2. * (C1GQ(N) + AD.gq0 * (LF - LQ)) * AD.qg0; - - // constant terms when pt->0 - // TODO: Implement/check expression of h2gg & C2GG!!! - /* long double h2gg = 0.; */ - /* std::complex HH2GG = std::pow(C1GG(N),2)+2.*C2GG(N) \ */ - /* +2.*C1GG(N)*h1gg+h2gg+2.*AD.gg1*LF+Beta0*AD.gg0*std::pow(LF,2) \ */ - /* +(Apt1g*Beta0*std::pow(LQ,3))/6.+(std::pow(LQ,2)*(Apt2g+Beta0* \ */ - /* (-Bpt1g-2.*AD.gg0-Apt1g*LQ)))/2.-(Beta1*LR+(std::pow(Beta0,2) \ */ - /* *std::pow(LR,2))/2.)*pc-Beta0*LR*(2.*C1GG(N)+h1gg+2.*AD.gg0* \ */ - /* (LF-LQ)-LQ*(Bpt1g+(Bpt1g*LQ)/2.)-Beta0*LR*pc)+AD.gg0*(LF-LQ)* \ */ - /* (4.*C1GG(N)+2.*h1gg+2.*AD.gg0*(LF-LQ)-LQ*(Bpt1g+(Bpt1g*LQ)/2.) \ */ - /* -Beta0*LR*pc)*(LQ*(Bpt1g+(Apt1g*LQ)/2.)+Beta0*LR*pc)-LQ*(Bpt2g \ */ - /* +2.*AD.gg1+Apt2g*LQ-2.*C1GG(N)*Beta0); */ - std::complex HH2GG = zero; - - // NLO contributions - result += aass * aass * - (Sigma24gg * LC4(xp) + Sigma23gg * LC3(xp) + - Sigma22gg * LC2(xp) + Sigma21gg * LC1(xp) + HH2GG); - } - if ((CHANNEL == 1) || (CHANNEL == 4)) { - result += zero; - } - if ((CHANNEL == 2) || (CHANNEL == 4)) { - result += zero; - } - } break; + if ((CHANNEL == 0) || (CHANNEL == 4)) // gg-channel or ALL + { + // pt-enhanced Sigma terms + std::complex Sigma24gg = std::pow(Apt1g, 2) / 8.; + std::complex Sigma23gg = + -Apt1g * + (Beta0 / 3. + 1. / 2. * (-Bpt1g - Apt1g * LQ - 2. * AD.gg0)); + std::complex Sigma22gg = + (-Apt2g - (-Bpt1g - 2. * AD.gg0 - Apt1g * LQ) * + (Bpt1g - Beta0 + Apt1g * LQ)) / + 2. - + (Apt1g * (2 * C1GG(N) + h1gg + 2. * AD.gg0 * (LF - LQ) - + LQ * (Bpt1g + (Bpt1g * LQ) / 2.) - (Beta0 * (-LQ + LR)) - + Beta0 * LR * pc)) / + 2. + + (-2. * AD.gg0 * (-Bpt1g - 2. * AD.gg0 - Apt1g * LQ) + + 2. * AD.gq0 * AD.qg0) / + 2.; + std::complex Sigma21gg = + -Bpt2g + 2. * Beta0 * C1GG(N) - 2. * AD.gg1 - Apt2g * LQ + + Beta0 * (-Bpt1g - AD.gg0 - Apt1g * LQ) * (LQ - LR) - + (Bpt1g + 2. * AD.gg0 + Apt1g * LQ) * + (2. * C1GG(N) + h1gg + 2. * AD.gg0 * (LF - LQ) - + LQ * (Bpt1g + (Bpt1g * LQ) / 2.) - Beta0 * LR * pc) - + 2. * (C1GQ(N) + AD.gq0 * (LF - LQ)) * AD.qg0; + + // constant terms when pt->0 + // TODO: Implement/check expression of h2gg & C2GG!!! + /* long double h2gg = 0.; */ + /* std::complex HH2GG = std::pow(C1GG(N),2)+2.*C2GG(N) \ */ + /* +2.*C1GG(N)*h1gg+h2gg+2.*AD.gg1*LF+Beta0*AD.gg0*std::pow(LF,2) \ + */ + /* +(Apt1g*Beta0*std::pow(LQ,3))/6.+(std::pow(LQ,2)*(Apt2g+Beta0* \ + */ + /* (-Bpt1g-2.*AD.gg0-Apt1g*LQ)))/2.-(Beta1*LR+(std::pow(Beta0,2) \ + */ + /* *std::pow(LR,2))/2.)*pc-Beta0*LR*(2.*C1GG(N)+h1gg+2.*AD.gg0* \ */ + /* (LF-LQ)-LQ*(Bpt1g+(Bpt1g*LQ)/2.)-Beta0*LR*pc)+AD.gg0*(LF-LQ)* \ + */ + /* (4.*C1GG(N)+2.*h1gg+2.*AD.gg0*(LF-LQ)-LQ*(Bpt1g+(Bpt1g*LQ)/2.) \ + */ + /* -Beta0*LR*pc)*(LQ*(Bpt1g+(Apt1g*LQ)/2.)+Beta0*LR*pc)-LQ*(Bpt2g \ + */ + /* +2.*AD.gg1+Apt2g*LQ-2.*C1GG(N)*Beta0); */ + std::complex HH2GG = zero; + + // NLO contributions + result += aass * aass * + (Sigma24gg * LC4(xp) + Sigma23gg * LC3(xp) + + Sigma22gg * LC2(xp) + Sigma21gg * LC1(xp) + HH2GG); + } + if ((CHANNEL == 1) || (CHANNEL == 4)) { + result += zero; + } + if ((CHANNEL == 2) || (CHANNEL == 4)) { + result += zero; + } + } break; } return 2. * pt / MH2 * SIGMA0 * result; diff --git a/src/ThresExp.cpp b/src/ThresExp.cpp index dc53340..0b23f4a 100644 --- a/src/ThresExp.cpp +++ b/src/ThresExp.cpp @@ -16,11 +16,12 @@ * ===================================================================================== */ -#include +#include "../include/ThresExp.h" + #include -#include -#include "../include/ThresExp.h" +#include +#include ThresExp::ThresExp(int order, int channel, void *params) { PhysParams param = *reinterpret_cast(params); @@ -41,28 +42,29 @@ ThresExp::ThresExp(int order, int channel, void *params) { LR = std::log(MH2 / MUR2); aass = static_cast( - param.alphas); // TODO; long double-check PI normalization + param.alphas); // TODO; long double-check PI normalization SROOT = static_cast(param.sroot); SIGMA0 = static_cast(param.sigma0); // Compute Beta Functions - Beta0 = (11.*CA-2.*NF)/(12.*M_PIl); - Beta1 = ((17.*CA*CA-5.*CA*NF-3.*CF*NF)*2./3.)/ - (16.*M_PIl*M_PIl); - Beta2 = ((2857./54.*CA*CA*CA+(CF*CF-205./18.*CF*CA-1415. \ - /54.*CA*CA)*NF+(11./9.*CF+79./54.*CA)*NF*NF))/ \ + Beta0 = (11. * CA - 2. * NF) / (12. * M_PIl); + Beta1 = ((17. * CA * CA - 5. * CA * NF - 3. * CF * NF) * 2. / 3.) / + (16. * M_PIl * M_PIl); + Beta2 = ((2857. / 54. * CA * CA * CA + + (CF * CF - 205. / 18. * CF * CA - 1415. / 54. * CA * CA) * NF + + (11. / 9. * CF + 79. / 54. * CA) * NF * NF)) / std::pow(4. * M_PIl, 3); // One & too loop Cusp Anomalous Dimensions // In the following, the 1/PI is contained in the definitions // of the cusp Anomalous dimensions not in the normalization of as Bth1g = -Beta0; - Ath1g = CA/M_PIl; - Ath1q = CF/M_PIl; - Bth1q = -3./4.*CF/M_PIl; - Ath2g = (CA/2.*(CA*(67./18.-zeta2)-5./9.*NF))/ + Ath1g = CA / M_PIl; + Ath1q = CF / M_PIl; + Bth1q = -3. / 4. * CF / M_PIl; + Ath2g = (CA / 2. * (CA * (67. / 18. - zeta2) - 5. / 9. * NF)) / std::pow(M_PIl, 2); - Ath2q = (CF/2.*(CA*(67./18.-zeta2)-5./9.*NF))/ + Ath2q = (CF / 2. * (CA * (67. / 18. - zeta2) - 5. / 9. * NF)) / std::pow(M_PIl, 2); } @@ -80,18 +82,25 @@ std::complex ThresExp::LOgggH(std::complex NN, std::complex CLOgggH; std::complex N = NN; std::complex half(0.5, 0.); - std::complex xprad(std::pow(std::sqrt(1.+xp)-std::sqrt(xp),4),0.); - - CLOgggH = 2.*aass*CA/std::sqrt(M_PIl)*1./xp* - std::exp(LogGamma(N)-LogGamma(N+0.5))*(Hyp2F1(half,N,N+0.5,xprad) \ - -2.*(1.+xp)/(std::pow(std::sqrt(1.+xp)+std::sqrt(xp),2.))*N/(N+0.5) \ - *Hyp2F1(half,N+1.,N+1.5,xprad)+((1.+xp)*(3.+xp))/ \ - (std::pow(std::sqrt(1.+xp)+std::sqrt(xp),4.))*N*(N+1.)/ - ((N+0.5)*(N+1.5))*Hyp2F1(half,N+2.,N+2.5,xprad)-2.*(1.+xp) \ - /(std::pow(std::sqrt(1.+xp)+std::sqrt(xp),6.))*N* \ - (N+1.)*(N+2.)/((N+0.5)*(N+1.5)*(N+2.5))*Hyp2F1(half,N+3.,N+3.5,xprad) \ - +1./(std::pow(std::sqrt(1.+xp)+std::sqrt(xp),8))*N*(N+1.)*(N+2.) \ - *(N+3.)/((N+0.5)*(N+1.5)*(N+2.5)*(N+3.5))*Hyp2F1(half,N+4.,N+4.5,xprad)); + std::complex xprad( + std::pow(std::sqrt(1. + xp) - std::sqrt(xp), 4), 0.); + + CLOgggH = + 2. * aass * CA / std::sqrt(M_PIl) * 1. / xp * + std::exp(LogGamma(N) - LogGamma(N + 0.5)) * + (Hyp2F1(half, N, N + 0.5, xprad) - + 2. * (1. + xp) / (std::pow(std::sqrt(1. + xp) + std::sqrt(xp), 2.)) * N / + (N + 0.5) * Hyp2F1(half, N + 1., N + 1.5, xprad) + + ((1. + xp) * (3. + xp)) / + (std::pow(std::sqrt(1. + xp) + std::sqrt(xp), 4.)) * N * (N + 1.) / + ((N + 0.5) * (N + 1.5)) * Hyp2F1(half, N + 2., N + 2.5, xprad) - + 2. * (1. + xp) / (std::pow(std::sqrt(1. + xp) + std::sqrt(xp), 6.)) * N * + (N + 1.) * (N + 2.) / ((N + 0.5) * (N + 1.5) * (N + 2.5)) * + Hyp2F1(half, N + 3., N + 3.5, xprad) + + 1. / (std::pow(std::sqrt(1. + xp) + std::sqrt(xp), 8)) * N * (N + 1.) * + (N + 2.) * (N + 3.) / + ((N + 0.5) * (N + 1.5) * (N + 2.5) * (N + 3.5)) * + Hyp2F1(half, N + 4., N + 4.5, xprad)); return CLOgggH; } @@ -103,16 +112,21 @@ std::complex ThresExp::LOgqqH(std::complex NN, // TODO: verify correspondence with small-pt std::complex N = NN; std::complex half(0.5, 0.); - std::complex xprad(std::pow(std::sqrt(1.+xp)-std::sqrt(xp),4.),0.); - - CLOgqqH = aass*CF/std::sqrt(M_PIl)*1./xp* - std::exp(LogGamma(N)-LogGamma(N+0.5))*(Hyp2F1(half,N,N+0.5,xprad) \ - -(4.+3.*xp)/(std::pow(std::sqrt(1.+xp)+std::sqrt(xp),2.))*N/(N+0.5) \ - *Hyp2F1(half, N+1.,N+1.5,xprad)+3.*(1.+xp)/(std::pow(std::sqrt(1.+xp) - +std::sqrt(xp),4.))*N*(N+1.)/((N+0.5)*(N+1.5))* \ - Hyp2F1(half,N+2.,N+2.5,xprad)-1./(std::pow(std::sqrt(1.+xp)+ \ - std::sqrt(xp),6.))*N*(N+1.)*(N+2.)/((N+0.5)*(N+1.5)*(N+2.5))* \ - Hyp2F1(half,N+3.,N+3.5,xprad)); + std::complex xprad( + std::pow(std::sqrt(1. + xp) - std::sqrt(xp), 4.), 0.); + + CLOgqqH = + aass * CF / std::sqrt(M_PIl) * 1. / xp * + std::exp(LogGamma(N) - LogGamma(N + 0.5)) * + (Hyp2F1(half, N, N + 0.5, xprad) - + (4. + 3. * xp) / (std::pow(std::sqrt(1. + xp) + std::sqrt(xp), 2.)) * N / + (N + 0.5) * Hyp2F1(half, N + 1., N + 1.5, xprad) + + 3. * (1. + xp) / (std::pow(std::sqrt(1. + xp) + std::sqrt(xp), 4.)) * N * + (N + 1.) / ((N + 0.5) * (N + 1.5)) * + Hyp2F1(half, N + 2., N + 2.5, xprad) - + 1. / (std::pow(std::sqrt(1. + xp) + std::sqrt(xp), 6.)) * N * (N + 1.) * + (N + 2.) / ((N + 0.5) * (N + 1.5) * (N + 2.5)) * + Hyp2F1(half, N + 3., N + 3.5, xprad)); return CLOgqqH; } @@ -124,15 +138,18 @@ std::complex ThresExp::LOqqgH(std::complex NN, // TODO: verify correspondence with small-pt std::complex N = NN; std::complex half(0.5, 0.); - std::complex xprad(std::pow(std::sqrt(1.+xp)-std::sqrt(xp),4.),0.); - - CLOqqgH = 2.*aass*CF*CF/std::sqrt(M_PIl)*1./ \ - (std::pow(std::sqrt(1.+xp)+std::sqrt(xp),2.))* \ - std::exp(LogGamma(N)-LogGamma(N+0.5))*(Hyp2F1(half,N,N+0.5,xprad) \ - -2.*(1.+xp)/(std::pow(std::sqrt(1.+xp)+std::sqrt(xp),2.))*N/(N+0.5) - *Hyp2F1(half,N+1.,N+1.5,xprad)+1./(std::pow(std::sqrt(1.+xp) \ - +std::sqrt(xp),4.))*N*(N+1.)/((N+0.5)*(N+1.5))* \ - Hyp2F1(half,N+2.,N+2.5,xprad)); + std::complex xprad( + std::pow(std::sqrt(1. + xp) - std::sqrt(xp), 4.), 0.); + + CLOqqgH = + 2. * aass * CF * CF / std::sqrt(M_PIl) * 1. / + (std::pow(std::sqrt(1. + xp) + std::sqrt(xp), 2.)) * + std::exp(LogGamma(N) - LogGamma(N + 0.5)) * + (Hyp2F1(half, N, N + 0.5, xprad) - + 2. * (1. + xp) / (std::pow(std::sqrt(1. + xp) + std::sqrt(xp), 2.)) * N / + (N + 0.5) * Hyp2F1(half, N + 1., N + 1.5, xprad) + + 1. / (std::pow(std::sqrt(1. + xp) + std::sqrt(xp), 4.)) * N * (N + 1.) / + ((N + 0.5) * (N + 1.5)) * Hyp2F1(half, N + 2., N + 2.5, xprad)); return CLOqqgH; } @@ -145,12 +162,13 @@ std::complex ThresExp::LOqqgH(std::complex NN, long double ThresExp::Sigma22ggg(long double xp) { return 3. * Ath1g / 8.; } long double ThresExp::Sigma21ggg(long double xp) { - long double xQp = std::sqrt(1.+xp)+std::sqrt(xp); - return -Bth1g/2.-Ath1g/2.*std::log(std::sqrt(xp)/xQp)-2.*Ath1g*LF; + long double xQp = std::sqrt(1. + xp) + std::sqrt(xp); + return -Bth1g / 2. - Ath1g / 2. * std::log(std::sqrt(xp) / xQp) - + 2. * Ath1g * LF; } long double ThresExp::Sigma20ggg(long double xp) { - return -3./2.*Ath1g*zeta2; + return -3. / 2. * Ath1g * zeta2; } //==========================================================================================// @@ -158,55 +176,72 @@ long double ThresExp::Sigma20ggg(long double xp) { //------------------------------------------------------------------------------------------// // gg->g long double ThresExp::GOgggH(long double xp) { - const long double g1gg = 67./36.*CA-5./18.*NF+CA*gsl_sf_zeta(2.) \ - -Beta0*M_PIl*std::log(xp/(1.+xp))-1./8.*CA* \ - std::pow(std::log(xp/(1.+xp)),2)+2.*CA* \ - gsl_sf_dilog(1.-std::sqrt(xp/(1.+xp)))+CA* \ - std::log(1.-std::sqrt(xp/(1.+xp)))*std::log(xp/(1.+xp)) \ - -0.5*CA*std::log(1.+std::sqrt(xp/(1.+xp)))*std::log(xp/(1.+xp)) \ - +0.5*CA*std::pow(std::log(1.+std::sqrt(xp/(1.+xp))),2)+2.*Beta0 \ - *M_PIl*std::pow(std::log(1.+std::sqrt(xp/(1.+xp))),2)+ \ - CA*gsl_sf_dilog((2.*std::sqrt(xp))/(std::sqrt(1.+xp)+std::sqrt(xp))) \ - -((CA-NF)*(std::sqrt(xp)*std::sqrt(1.+xp)*(1.+xp)-2.*xp-xp*xp))/ \ - (6.*(1.+8.*xp+9.*xp*xp)); - return g1gg/M_PIl; + const long double g1gg = + 67. / 36. * CA - 5. / 18. * NF + CA * gsl_sf_zeta(2.) - + Beta0 * M_PIl * std::log(xp / (1. + xp)) - + 1. / 8. * CA * std::pow(std::log(xp / (1. + xp)), 2) + + 2. * CA * gsl_sf_dilog(1. - std::sqrt(xp / (1. + xp))) + + CA * std::log(1. - std::sqrt(xp / (1. + xp))) * std::log(xp / (1. + xp)) - + 0.5 * CA * std::log(1. + std::sqrt(xp / (1. + xp))) * + std::log(xp / (1. + xp)) + + 0.5 * CA * std::pow(std::log(1. + std::sqrt(xp / (1. + xp))), 2) + + 2. * Beta0 * M_PIl * + std::pow(std::log(1. + std::sqrt(xp / (1. + xp))), 2) + + CA * gsl_sf_dilog((2. * std::sqrt(xp)) / + (std::sqrt(1. + xp) + std::sqrt(xp))) - + ((CA - NF) * + (std::sqrt(xp) * std::sqrt(1. + xp) * (1. + xp) - 2. * xp - xp * xp)) / + (6. * (1. + 8. * xp + 9. * xp * xp)); + return g1gg / M_PIl; } // gq->q long double ThresExp::GOgqqH(long double xp) { - const long double g1gq = -7./4.*CF+134./36.*CA-20./36.*NF \ - -8.*CF*gsl_sf_zeta(2.)+12.*CA*gsl_sf_zeta(2.)-4.*Beta0 \ - *M_PIl*std::log(xp/(1.+xp))+3./2.*CF*std::log(xp/(1.+xp)) \ - -0.5*CA*std::pow(std::log(xp/(1.+xp)),2)+4.*(CF+CA)* \ - gsl_sf_dilog(1.-std::sqrt(xp/(1.+xp)))+(2.*(CA-CF)* \ - (1.+3.*xp+3.*std::sqrt(xp*(1.+xp))))/(2.*std::sqrt(xp* \ - (1.+xp))+1.+3.*xp)+8.*Beta0*M_PIl* - std::log(1.+std::sqrt(xp/(1.+xp)))-3.*CF* \ - std::log(1.+std::sqrt(xp/(1.+xp)))+2.*CF* \ - std::log(1.-std::sqrt(xp/(1.+xp)))*std::log(xp/(1.+xp)) \ - +2.*CA*std::log(1.-std::sqrt(xp/(1.+xp)))*std::log(xp/(1.+xp)) \ - -2.*CF*std::log(1.+std::sqrt(xp/(1.+xp)))*std::log(xp/(1.+xp)) \ - -2.*CF*std::pow(std::log(1.+std::sqrt(xp/(1.+xp))),2)+4.*CF \ - *gsl_sf_dilog(2.*std::sqrt(xp)/(std::sqrt(1.+xp)+std::sqrt(xp))); - return g1gq/M_PIl; + const long double g1gq = + -7. / 4. * CF + 134. / 36. * CA - 20. / 36. * NF - + 8. * CF * gsl_sf_zeta(2.) + 12. * CA * gsl_sf_zeta(2.) - + 4. * Beta0 * M_PIl * std::log(xp / (1. + xp)) + + 3. / 2. * CF * std::log(xp / (1. + xp)) - + 0.5 * CA * std::pow(std::log(xp / (1. + xp)), 2) + + 4. * (CF + CA) * gsl_sf_dilog(1. - std::sqrt(xp / (1. + xp))) + + (2. * (CA - CF) * (1. + 3. * xp + 3. * std::sqrt(xp * (1. + xp)))) / + (2. * std::sqrt(xp * (1. + xp)) + 1. + 3. * xp) + + 8. * Beta0 * M_PIl * std::log(1. + std::sqrt(xp / (1. + xp))) - + 3. * CF * std::log(1. + std::sqrt(xp / (1. + xp))) + + 2. * CF * std::log(1. - std::sqrt(xp / (1. + xp))) * + std::log(xp / (1. + xp)) + + 2. * CA * std::log(1. - std::sqrt(xp / (1. + xp))) * + std::log(xp / (1. + xp)) - + 2. * CF * std::log(1. + std::sqrt(xp / (1. + xp))) * + std::log(xp / (1. + xp)) - + 2. * CF * std::pow(std::log(1. + std::sqrt(xp / (1. + xp))), 2) + + 4. * CF * + gsl_sf_dilog(2. * std::sqrt(xp) / + (std::sqrt(1. + xp) + std::sqrt(xp))); + return g1gq / M_PIl; } // qq->g long double ThresExp::GOqqgH(long double xp) { - const long double g1qq = -9./2.*CF+79./12.*CA-5./6.*NF \ - +12.*CF*gsl_sf_zeta(2.)-10.*CA*gsl_sf_zeta(2.) \ - -(CF-CA)*std::sqrt(1.+xp)/std::sqrt(xp)+4.*CF \ - *gsl_sf_dilog(1.-std::sqrt(xp/(1.+xp)))-3./4.*CF* \ - std::log(xp/(1.+xp))-Beta0*M_PIl*std::log(xp/(1.+xp)) \ - +0.25*CA*std::pow(std::log(xp/(1.+xp)),2)-0.5*CF* \ - std::pow(std::log(xp/(1.+xp)),2)+2.*CF* \ - std::log(1.-std::sqrt(xp/(1.+xp)))*std::log(xp/(1.+xp)) \ - +1.5*CF*std::log(1.+std::sqrt(xp/(1.+xp)))+2.*Beta0*M_PIl \ - *std::log(1.+std::sqrt(xp/(1.+xp)))+CA \ - *std::pow(std::log(1.+std::sqrt(xp/(1.+xp))),2)-CA* \ - std::log(1.+std::sqrt(xp/(1.+xp)))*std::log(xp/(1.+xp))+2.*CA \ - *gsl_sf_dilog(2.*std::sqrt(xp)/(std::sqrt(1.+xp)+std::sqrt(xp))); - return g1qq/M_PIl; + const long double g1qq = + -9. / 2. * CF + 79. / 12. * CA - 5. / 6. * NF + + 12. * CF * gsl_sf_zeta(2.) - 10. * CA * gsl_sf_zeta(2.) - + (CF - CA) * std::sqrt(1. + xp) / std::sqrt(xp) + + 4. * CF * gsl_sf_dilog(1. - std::sqrt(xp / (1. + xp))) - + 3. / 4. * CF * std::log(xp / (1. + xp)) - + Beta0 * M_PIl * std::log(xp / (1. + xp)) + + 0.25 * CA * std::pow(std::log(xp / (1. + xp)), 2) - + 0.5 * CF * std::pow(std::log(xp / (1. + xp)), 2) + + 2. * CF * std::log(1. - std::sqrt(xp / (1. + xp))) * + std::log(xp / (1. + xp)) + + 1.5 * CF * std::log(1. + std::sqrt(xp / (1. + xp))) + + 2. * Beta0 * M_PIl * std::log(1. + std::sqrt(xp / (1. + xp))) + + CA * std::pow(std::log(1. + std::sqrt(xp / (1. + xp))), 2) - + CA * std::log(1. + std::sqrt(xp / (1. + xp))) * std::log(xp / (1. + xp)) + + 2. * CA * + gsl_sf_dilog(2. * std::sqrt(xp) / + (std::sqrt(1. + xp) + std::sqrt(xp))); + return g1qq / M_PIl; } //==========================================================================================// @@ -217,43 +252,43 @@ long double ThresExp::GOqqgH(long double xp) { std::complex ThresExp::ThresExpExpr(std::complex N, long double pt) { // TODO: re-check definition MH2 vs. Qs2 - long double xp = std::pow(pt,2)/MH2; - std::complex zero(0.,0.); - std::complex result(0.,0.); + long double xp = std::pow(pt, 2) / MH2; + std::complex zero(0., 0.); + std::complex result(0., 0.); - std::complex Nbar = N*std::exp(EulerGamma); - std::complex LNbar = 2.*aass*Beta0*std::log(Nbar); + std::complex Nbar = N * std::exp(EulerGamma); + std::complex LNbar = 2. * aass * Beta0 * std::log(Nbar); switch (ORD) { - // TODO: Match CASES with HpT-MON - case (0): // order as^1 - { - if ((CHANNEL==0) || (CHANNEL==5)) { - result += LOgggH(N, xp); - } - if ((CHANNEL==1) || (CHANNEL==5)) { - result += LOgqqH(N, xp); - } - if ((CHANNEL==2) || (CHANNEL==5)) { - result += LOqqgH(N, xp); - } - } break; - case (1): // order as^2 - { - if ((CHANNEL==0) || (CHANNEL==5)) { - std::complex SIGMAGG = - Sigma22ggg(xp)*std::pow(LNbar,2) + - Sigma21ggg(xp)*LNbar; /* +Sigma20ggg(xp); */ - result += aass*LOgggH(N, xp)*(GOgggH(xp)+SIGMAGG); - } - if ((CHANNEL==1) || (CHANNEL==5)) { - return zero; - } - if ((CHANNEL==2) || (CHANNEL==5)) { - return zero; - } - } break; + // TODO: Match CASES with HpT-MON + case (0): // order as^1 + { + if ((CHANNEL == 0) || (CHANNEL == 5)) { + result += LOgggH(N, xp); + } + if ((CHANNEL == 1) || (CHANNEL == 5)) { + result += LOgqqH(N, xp); + } + if ((CHANNEL == 2) || (CHANNEL == 5)) { + result += LOqqgH(N, xp); + } + } break; + case (1): // order as^2 + { + if ((CHANNEL == 0) || (CHANNEL == 5)) { + std::complex SIGMAGG = + Sigma22ggg(xp) * std::pow(LNbar, 2) + + Sigma21ggg(xp) * LNbar; /* +Sigma20ggg(xp); */ + result += aass * LOgggH(N, xp) * (GOgggH(xp) + SIGMAGG); + } + if ((CHANNEL == 1) || (CHANNEL == 5)) { + return zero; + } + if ((CHANNEL == 2) || (CHANNEL == 5)) { + return zero; + } + } break; } - return 2.*pt/MH2*SIGMA0*result; + return 2. * pt / MH2 * SIGMA0 * result; }