From 8d3a8b22f8059532bd0f6eec067ad6f36d16c7b4 Mon Sep 17 00:00:00 2001
From: Klaus Rabbertz <klaus.rabbertz@cern.ch>
Date: Fri, 31 May 2024 18:05:16 +0200
Subject: [PATCH] Implement 30 scale variations for ratio
---
v2.5/toolkit/fastnlotoolkit/fastNLOReader.cc | 63 ++++++++++++++-----
.../include/fastnlotk/fastNLOConstants.h.in | 5 +-
.../include/fastnlotk/fastNLOReader.h | 2 +-
v2.5/toolkit/src/fnlo-tk-yodaout.cc | 11 +++-
4 files changed, 60 insertions(+), 21 deletions(-)
diff --git a/v2.5/toolkit/fastnlotoolkit/fastNLOReader.cc b/v2.5/toolkit/fastnlotoolkit/fastNLOReader.cc
index a04cc405..816e515b 100644
--- a/v2.5/toolkit/fastnlotoolkit/fastNLOReader.cc
+++ b/v2.5/toolkit/fastnlotoolkit/fastNLOReader.cc
@@ -1028,7 +1028,7 @@ vector < double > fastNLOReader::GetUncertainty(bool lNorm) {
}
//______________________________________________________________________________
-vector < double > fastNLOReader::GetNormCrossSection() {
+vector < double > fastNLOReader::GetNormCrossSection(bool lNormScale, double xmurd, double xmufd) {
// Check whether normalization is defined
if (INormFlag == 0) {
logger.error["GetNormCrossSection"]<<"Normalization not defined for this scenario, aborting!"<<endl;
@@ -1038,6 +1038,13 @@ vector < double > fastNLOReader::GetNormCrossSection() {
if (XSection.empty()) CalcCrossSection();
vector < double > XSectionNorm = XSection;
+ // Recalculate with modified scale factors for normalisation, if requested
+ if (lNormScale) {
+ SetScaleFactorsMuRMuF(xmurd,xmufd);
+ CalcCrossSection();
+ }
+ vector < double > XSectionDen = XSection;
+
// // Second table to be loaded?
// if ( INormFlag < 0 ) {
// string denomtable = GetDenomTable();
@@ -1082,7 +1089,7 @@ vector < double > fastNLOReader::GetNormCrossSection() {
}
for (int in = idivlo; in <= idivup; in++) {
double bwidth = GetObsBinUpBound(in,iDim) - GetObsBinLoBound(in,iDim);
- xsnorm += XSection[in]*bwidth;
+ xsnorm += XSectionDen[in]*bwidth;
twidth += bwidth;
}
}
@@ -3402,22 +3409,35 @@ double fastNLOReader::RescaleCrossSectionUnits(double binsize, int xunits) {
// Scale uncertainty
//______________________________________________________________________________
XsUncertainty fastNLOReader::GetXsUncertainty(const EScaleUncertaintyStyle eScaleUnc, bool lNorm, double sclfac) {
- // Get 2-, 6- or, 31-point scale uncertainty around sclfac * central scale (31 only for normalised x sections/ratios)
- const double xmur0[7] = {1.0, 0.5, 2.0, 0.5, 1.0, 1.0, 2.0};
- const double xmuf0[7] = {1.0, 0.5, 2.0, 1.0, 0.5, 2.0, 1.0};
- // const double xmur0[7] = {1.0, 0.5, 2.0, 0.5, 1.0, 1.0, 2.0, 0.5, 2.0, 0.5, 1.0, 1.0, 2.0};
- // const double xmuf0[7] = {1.0, 0.5, 2.0, 1.0, 0.5, 2.0, 1.0, 0.5, 2.0, 1.0, 0.5, 2.0, 1.0};
- // const double xmurn[7] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.5, 2.0, 0.5, 1.0, 1.0, 2.0};
- // const double xmufn[7] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.5, 2.0, 1.0, 0.5, 2.0, 1.0};
+ // Get 2-, 6- or, 30-point scale uncertainty around sclfac * central scale (30 only for normalised x sections/ratios)
+ // const double xmur0[7] = {1.0, 0.5, 2.0, 0.5, 1.0, 1.0, 2.0};
+ // const double xmuf0[7] = {1.0, 0.5, 2.0, 1.0, 0.5, 2.0, 1.0};
+ const double xmurn0[31] = {1.0, 0.5, 2.0, 0.5, 1.0, 1.0, 2.0,
+ 0.5, 0.5, 1.0, 1.0, 0.5, 0.5, 0.5, 1.0, 1.0, 1.0, 1.0, 0.5,
+ 2.0, 2.0, 1.0, 1.0, 2.0, 2.0, 2.0, 1.0, 1.0, 1.0, 1.0, 2.0};
+ const double xmufn0[31] = {1.0, 0.5, 2.0, 1.0, 0.5, 2.0, 1.0,
+ 0.5, 1.0, 1.0, 0.5, 0.5, 0.5, 1.0, 0.5, 1.0, 1.0, 0.5, 1.0,
+ 2.0, 1.0, 1.0, 2.0, 2.0, 2.0, 1.0, 2.0, 1.0, 1.0, 2.0, 1.0};
+ const double xmurd0[31] = {1.0, 0.5, 2.0, 0.5, 1.0, 1.0, 2.0,
+ 1.0, 1.0, 0.5, 0.5, 0.5, 1.0, 0.5, 0.5, 1.0, 0.5, 1.0, 1.0,
+ 1.0, 1.0, 2.0, 2.0, 2.0, 1.0, 2.0, 2.0, 1.0, 2.0, 1.0, 1.0};
+ const double xmufd0[31] = {1.0, 0.5, 2.0, 1.0, 0.5, 2.0, 1.0,
+ 1.0, 0.5, 0.5, 1.0, 1.0, 0.5, 0.5, 0.5, 0.5, 1.0, 1.0, 1.0,
+ 1.0, 2.0, 2.0, 1.0, 1.0, 2.0, 2.0, 2.0, 2.0, 1.0, 1.0, 1.0};
+
if ( sclfac < 0.1 || sclfac > 10. ) {
logger.error["GetScaleUncertainty"]<<"ERROR! Illegal value for sclfac, exiting. sclfac = "<< sclfac <<endl;
exit(1);
}
- double xmurs[7];
- double xmufs[7];
- for (int i=0; i<7; i++) {
- xmurs[i] = sclfac*xmur0[i];
- xmufs[i] = sclfac*xmuf0[i];
+ double xmurn[31];
+ double xmufn[31];
+ double xmurd[31];
+ double xmufd[31];
+ for (int i=0; i<31; i++) {
+ xmurn[i] = sclfac*xmurn0[i];
+ xmufn[i] = sclfac*xmufn0[i];
+ xmurd[i] = sclfac*xmurd0[i];
+ xmufd[i] = sclfac*xmufd0[i];
}
XsUncertainty XsUnc;
@@ -3427,6 +3447,8 @@ XsUncertainty fastNLOReader::GetXsUncertainty(const EScaleUncertaintyStyle eScal
npoint = 2;
} else if (eScaleUnc == kAsymmetricSixPoint) {
npoint = 6;
+ } else if (eScaleUnc == kAsymmetricRatio) {
+ npoint = 30;
}
logger.debug["GetScaleUncertainty"]<<"npoint = "<<npoint<<endl;
@@ -3436,18 +3458,25 @@ XsUncertainty fastNLOReader::GetXsUncertainty(const EScaleUncertaintyStyle eScal
logger.info["GetScaleUncertainty"]<<"Symmetric 2-point scale variations selected,"<<endl;
} else if (npoint == 6) {
logger.info["GetScaleUncertainty"]<<"Asymmetric 6-point scale variations selected,"<<endl;
+ } else if (npoint == 30 && lNorm) {
+ logger.info["GetScaleUncertainty"]<<"Asymmetric 30-point scale variations for ratios selected,"<<endl;
} else {
logger.error["GetScaleUncertainty"]<<"ERROR! No usual scale variation scheme selected, exiting."<<endl;
logger.error["GetScaleUncertainty"]<<"npoint = "<<npoint<<endl;
+ logger.error["GetScaleUncertainty"]<<"lNorm = "<<lNorm<<endl;
exit(1);
}
vector < double > MyXSection;
//! Cross section and absolute uncertainties
for (unsigned int iscl = 0; iscl <= npoint; iscl++) {
- SetScaleFactorsMuRMuF(xmurs[iscl],xmufs[iscl]);
+ SetScaleFactorsMuRMuF(xmurn[iscl],xmufn[iscl]);
CalcCrossSection();
- MyXSection = GetCrossSection(lNorm);
+ bool lNormScale = false;
+ if (lNorm && eScaleUnc == kAsymmetricRatio) {
+ lNormScale = true;
+ }
+ MyXSection = GetNormCrossSection(lNormScale,xmurd[iscl],xmufd[iscl]);
for (unsigned int iobs = 0; iobs < NObsBin; iobs++) {
if (iscl == 0) {
XsUnc.xs.push_back(MyXSection[iobs]);
@@ -3473,7 +3502,7 @@ XsUncertainty fastNLOReader::GetXsUncertainty(const EScaleUncertaintyStyle eScal
}
logger.info["GetScaleUncertainty"]<<"Setting scale factors back to default of unity."<<endl;
- SetScaleFactorsMuRMuF(xmurs[0],xmufs[0]);
+ SetScaleFactorsMuRMuF(xmurn[0],xmufn[0]);
return XsUnc;
}
diff --git a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOConstants.h.in b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOConstants.h.in
index 5a4ddb7b..8f225cfe 100644
--- a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOConstants.h.in
+++ b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOConstants.h.in
@@ -108,8 +108,9 @@ namespace fastNLO {
enum EScaleUncertaintyStyle {
kScaleNone = 0, // no scale uncertainty, only central scale (mu_r,mu_f) = (1,1) evaluated
kSymmetricTwoPoint = 1, // symmetric (mu_r,mu_f) scale variations by factors (1/2,1/2), (2,2)
- kAsymmetricSixPoint = 2 // asymmetric (mu_r,mu_f) scale variations by factors (1/2,1/2), (2,2) plus
+ kAsymmetricSixPoint = 2, // asymmetric (mu_r,mu_f) scale variations by factors (1/2,1/2), (2,2) plus
// (1/2,1), (1,1/2), (1,2), (2,1)
+ kAsymmetricRatio = 3 // 30(+1) combinations for mu_r, mu_f in numerator and denominator
};
constexpr std::string_view ScaleUncertaintyStyle_to_string(EScaleUncertaintyStyle estyle) {
@@ -117,6 +118,7 @@ namespace fastNLO {
case kScaleNone: return "NN";
case kSymmetricTwoPoint: return "2P";
case kAsymmetricSixPoint: return "6P";
+ case kAsymmetricRatio: return "30";
default: return "??";
}
}
@@ -125,6 +127,7 @@ namespace fastNLO {
if (sstyle == "NN") return kScaleNone;
if (sstyle == "2P") return kSymmetricTwoPoint;
if (sstyle == "6P") return kAsymmetricSixPoint;
+ if (sstyle == "30") return kAsymmetricRatio;
return{kScaleNone};
}
diff --git a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOReader.h b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOReader.h
index d76909a0..23ac03c4 100644
--- a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOReader.h
+++ b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOReader.h
@@ -89,7 +89,7 @@ public:
std::vector < double > GetCrossSection(bool lNorm = false); //!< Return vector with all cross section values, normalize on request
std::vector < double > GetUncertainty(bool lNorm = false); //!< Return vector with additional uncertainty of cross section values, normalise on request (NOT YET IMPLEMENTED)
- std::vector < double > GetNormCrossSection(); //!< Return vector with all normalized cross section values
+ std::vector < double > GetNormCrossSection(bool lNormScale = false, double xmurd = 1.0, double xmufd = 1.0); //!< Return vector with all normalized cross section values
std::vector < std::map< double, double > > GetCrossSection_vs_x1(); //! Cross section vs. x1 ( XSection_vsX1[bin][<x,xs>] )
std::vector < std::map< double, double > > GetCrossSection_vs_x2(); //! Cross section vs. x2 ( XSection_vsX1[bin][<x,xs>] )
diff --git a/v2.5/toolkit/src/fnlo-tk-yodaout.cc b/v2.5/toolkit/src/fnlo-tk-yodaout.cc
index 73fc42e5..6b6d299d 100644
--- a/v2.5/toolkit/src/fnlo-tk-yodaout.cc
+++ b/v2.5/toolkit/src/fnlo-tk-yodaout.cc
@@ -121,6 +121,7 @@ int main(int argc, char** argv) {
man << " Alternatives: NN (none, but correct MC sampling average value --> NNPDF PDFs)" << endl;
man << " 2P (symmetric 2-point scale factor variation)" << endl;
man << " 6P (asymmetric 6-point scale factor variation)" << endl;
+ man << " 30 (asymmetric 30-point scale factor variation for ratios)" << endl;
man << " HS (symmetric Hessian PDF uncertainty --> ABM, (G)JR PDFs)" << endl;
man << " HA (asymmetric Hessian PDF uncertainty)" << endl;
man << " HP (pairwise asymmetric Hessian PDF uncertainty --> CTEQ|MSTW PDFs)" << endl;
@@ -189,6 +190,9 @@ int main(int argc, char** argv) {
} else if ( chunc == "6P" ) {
eScaleUnc = kAsymmetricSixPoint;
shout["fnlo-tk-yodaout"] << "Showing uncertainty from asymmetric 6-point scale factor variation." << endl;
+ } else if ( chunc == "30" ) {
+ eScaleUnc = kAsymmetricRatio;
+ shout["fnlo-tk-yodaout"] << "Showing uncertainty from asymmetric 30-point scale factor variation for ratios." << endl;
} else if ( chunc == "HS" ) {
ePDFUnc = kHessianSymmetric;
shout["fnlo-tk-yodaout"] << "Showing symmetric Hessian PDF uncertainty (--> ABM, (G)JR PDFs)." << endl;
@@ -442,9 +446,12 @@ int main(int argc, char** argv) {
if ( fnlo->IsNorm() ) {
lNorm = true;
} else {
- error["fnlo-read"] << "Normalization requested but not defined for this table, aborted!" << endl;
+ error["fnlo-tk-yodaout"] << "Normalization requested but not defined for this table, aborted!" << endl;
exit(1);
}
+ } else if ( eScaleUnc == kAsymmetricRatio) {
+ error["fnlo-tk-yodaout"] << "Scale uncertainty for ratio only works with normalisation on, aborted!" << endl;
+ exit(1);
}
//! --- Determine dimensioning of observable bins in table
@@ -492,7 +499,7 @@ int main(int argc, char** argv) {
XsUncertainty XsUnc;
string LineName;
string UncName;
- if ( chunc == "2P" || chunc == "6P" ) {
+ if ( chunc == "2P" || chunc == "6P" || chunc == "30") {
XsUnc = fnlo->GetXsUncertainty(eScaleUnc, lNorm, sclfac);
UncName = " # Relative scale uncertainties";
LineName += "_dxscl";
--
GitLab