From 5ce4c5c0efc3f5820434cf94a73251426ca40c6c Mon Sep 17 00:00:00 2001
From: Klaus Rabbertz <klaus.rabbertz@cern.ch>
Date: Mon, 7 Feb 2022 06:00:07 +0100
Subject: [PATCH] A few fixes for stat. uncertainties with NLOjet++

---
 tools/fnlo-add-tables.pl            | 17 ++++++++++-------
 v2.5/toolkit/src/fnlo-tk-example.cc |  8 ++++----
 v2.5/toolkit/src/fnlo-tk-rootout.cc | 12 ++++++------
 v2.5/toolkit/src/fnlo-tk-statunc.cc | 15 +++++++++++----
 v2.5/toolkit/src/fnlo-tk-yodaout.cc | 16 +++++++++++-----
 5 files changed, 42 insertions(+), 26 deletions(-)

diff --git a/tools/fnlo-add-tables.pl b/tools/fnlo-add-tables.pl
index 638eb73e..a3658c5d 100755
--- a/tools/fnlo-add-tables.pl
+++ b/tools/fnlo-add-tables.pl
@@ -1,7 +1,7 @@
 #!/usr/bin/env perl
 #
 # fastNLO table addition script:
-#    This perl script creates either a total sum table (scen.tab) or
+#    This perl script creates either a total sum table (scen.tab.gz) or
 #    sum tables for statistical uncertainty evaluation:
 #    - For each single LO table with *one* 'dummy' NLO table
 #    - For each single NLO table with *all* LO tables or
@@ -11,6 +11,7 @@
 # created by K. Rabbertz: 05.03.2006
 # last modified:
 # adapted by K. Rabbertz from fastadd.pl: 01.08.2014
+# adapted by K. Rabbertz to use gzip'ped tables: 05.02.2022
 #
 #-----------------------------------------------------------------------
 #
@@ -56,6 +57,8 @@ unless ( @ARGV == 1 ) {
     die "fnlo-add-tables.pl: Error! Need one scenario name!\n";
 }
 my $scen   = shift;
+# Remove potential trailing slash
+$scen =~ s/\/$//;
 my $debug  = $opt_d;
 my $vers   = $opt_v;
 
@@ -205,7 +208,7 @@ if ( $opt_s && ! @alltabs ) {
             my $scmd1 = $cmd;
             $scmd1 .= " $lodir/$lotab $nlodir/$nlotabs[0]";
             $scmd1 .= " $lotab";
-            $scmd1 =~ s/\.${tabext}$/\.tab/;
+            $scmd1 =~ s/\.${tabext}$/\.tab\.gz/;
             print "fnlo-add-tables.pl: Creating sum table for $lotab ...\n";
             if ( $debug ) {print "fnlo-add-tables.pl: Running command $scmd1\n";}
             system("$scmd1 >> ${scen}_addst.log");
@@ -213,9 +216,9 @@ if ( $opt_s && ! @alltabs ) {
         $scmd2 .= " $lodir/$lotab";
     }
     my $losum = $lotabs[0];
-    $losum =~ s/\d\d\d\d\.tab/sum\.tab/;
+    $losum =~ s/\d+\.tab/sum\.tab/;
     $scmd2 .= " $losum";
-    $scmd2 =~ s/\.${tabext}$/\.tab/;
+    $scmd2 =~ s/\.${tabext}$/\.tab\.gz/;
     if ( -f $losum && ! -z $losum ) {
         print "fnlo-add-tables.pl: Non-empty all LO sum table $losum exists already, skipped!\n";
     } else {
@@ -232,18 +235,18 @@ if ( $opt_s && ! @alltabs ) {
             $scmd1 .= " $losum";
             $scmd1 .= " $nlodir/$nlotab";
             $scmd1 .= " $nlotab";
-            $scmd1 =~ s/\.${tabext}$/\.tab/;
+            $scmd1 =~ s/\.${tabext}$/\.tab\.gz/;
             print "fnlo-add-tables.pl: Creating sum table for $nlotab ...\n";
             if ( $debug ) {print "fnlo-add-tables.pl: Running command $scmd1\n";}
             system("$scmd1 >> ${scen}_addst.log");
         }
     }
-    $losum =~ s/\.${tabext}$/\.tab/;
+    $losum =~ s/\.${tabext}$/\.tab\.gz/;
     print "\nfnlo-add-tables.pl: Removing temporary LO sum table $losum ...\n";
     unlink $losum;
 # Normal mode: All LO with all NLO/NNLO tables
 } else {
-    my $scentab = ${scen}.".tab";
+    my $scentab = ${scen}.".tab.gz";
     if ( -f $scentab && ! -z $scentab ) {
         print "fnlo-add-tables.pl: Non-empty sum table for $scentab exists already, skipped!\n";
     } else {
diff --git a/v2.5/toolkit/src/fnlo-tk-example.cc b/v2.5/toolkit/src/fnlo-tk-example.cc
index dd37b9eb..ffa91444 100644
--- a/v2.5/toolkit/src/fnlo-tk-example.cc
+++ b/v2.5/toolkit/src/fnlo-tk-example.cc
@@ -158,8 +158,8 @@ int main(int argc, char** argv) {
    //   XsUnc  = fnlo.GetPDFUncertainty(ePDFUnc);
 
    cout << _CSEPSC << endl;
-   cout << " # Relative Scale Uncertainties (6P)" << endl;
-   cout << " # bin      cross section           lower uncertainty       upper uncertainty" << endl;
+   cout << " # Relative scale uncertainties (6P)" << endl;
+   cout << " # bin      cross_section           lower_uncertainty       upper_uncertainty" << endl;
    cout << _SSEPSC << endl;
    for ( unsigned int iobs=0;iobs<XsUnc.xs.size();iobs++ ) {
       printf("%5.i      %#18.11E      %#18.11E      %#18.11E\n",iobs+1,XsUnc.xs[iobs],XsUnc.dxsl[iobs],XsUnc.dxsu[iobs]);
@@ -184,8 +184,8 @@ int main(int argc, char** argv) {
    vector<double> central  = fnlo.CalcPDFUncertaintyCentral(PDFUnc);
 
    cout << _CSEPSC << endl;
-   cout << " # Relative and Absolute PDF Uncertainties via LHAPDF6" << endl;
-   cout << " # bin      cross section      lower rel. uncertainty   upper rel. uncertainty   lower abs. uncertainty   upper abs. uncertainty" << endl;
+   cout << " # Relative and absolute PDF uncertainties via LHAPDF6" << endl;
+   cout << " # bin      cross_section      lower_rel._uncertainty   upper_rel._uncertainty   lower_abs._uncertainty   upper_abs._uncertainty" << endl;
    cout << _SSEPSC << endl;
    for ( unsigned int iobs=0;iobs<central.size();iobs++ ) {
       printf("%5.i      %#18.11E      %#18.11E      %#18.11E      %#18.11E      %#18.11E\n",iobs+1,central[iobs],errdn[iobs],errup[iobs],errdnabs[iobs],errupabs[iobs]);
diff --git a/v2.5/toolkit/src/fnlo-tk-rootout.cc b/v2.5/toolkit/src/fnlo-tk-rootout.cc
index b29aad49..06e70cfc 100644
--- a/v2.5/toolkit/src/fnlo-tk-rootout.cc
+++ b/v2.5/toolkit/src/fnlo-tk-rootout.cc
@@ -476,7 +476,7 @@ int main(int argc, char** argv) {
                vector < double > xstmp = fnlo->GetCrossSection(lNorm);
                XsUnc = fnlo->GetPDFUncertainty(PDFUncStyles[iPDF], lNorm);
                if ( PDFUncStyles[iPDF] == kPDFNone ) XsUnc.xs = xstmp;
-               snprintf(buffer, sizeof(buffer), " # Relative PDF Uncertainties (%s %s %s)",sOrder.c_str(),PDFFiles[iPDF].c_str(),chunc.c_str());
+               snprintf(buffer, sizeof(buffer), " # Relative PDF uncertainties (%s %s %s)",sOrder.c_str(),PDFFiles[iPDF].c_str(),chunc.c_str());
                snprintf(titlel, sizeof(titlel), "-dsigma_%s/sigma",PDFFiles[iPDF].c_str());
                snprintf(titleu, sizeof(titleu), "+dsigma_%s/sigma",PDFFiles[iPDF].c_str());
             }
@@ -487,7 +487,7 @@ int main(int argc, char** argv) {
                }
                eAddUnc = kAddStat;
                XsUnc = fnlo->GetAddUncertainty(eAddUnc, lNorm);
-               snprintf(buffer, sizeof(buffer), " # Relative Statistical Uncertainties (%s %s)",sOrder.c_str(),PDFFiles[iPDF].c_str());
+               snprintf(buffer, sizeof(buffer), " # Relative statistical uncertainties (%s %s)",sOrder.c_str(),PDFFiles[iPDF].c_str());
                snprintf(titlel, sizeof(titlel), "-dsigma_stat/sigma");
                snprintf(titleu, sizeof(titleu), "+dsigma_stat/sigma");
             }
@@ -495,7 +495,7 @@ int main(int argc, char** argv) {
             else if ( iUnc==2 ) {
                eScaleUnc = kSymmetricTwoPoint;
                XsUnc = fnlo->GetScaleUncertainty(eScaleUnc, lNorm);
-               snprintf(buffer, sizeof(buffer), " # 2P Relative Scale Uncertainties (%s %s)",sOrder.c_str(),PDFFiles[iPDF].c_str());
+               snprintf(buffer, sizeof(buffer), " # 2P Relative scale uncertainties (%s %s)",sOrder.c_str(),PDFFiles[iPDF].c_str());
                snprintf(titlel, sizeof(titlel), "-dsigma_2P/sigma");
                snprintf(titleu, sizeof(titleu), "+dsigma_2P/sigma");
             }
@@ -503,7 +503,7 @@ int main(int argc, char** argv) {
             else if ( iUnc==3 ) {
                eScaleUnc = kAsymmetricSixPoint;
                XsUnc = fnlo->GetScaleUncertainty(eScaleUnc, lNorm);
-               snprintf(buffer, sizeof(buffer), " # 6P Relative Scale Uncertainties (%s %s)",sOrder.c_str(),PDFFiles[iPDF].c_str());
+               snprintf(buffer, sizeof(buffer), " # 6P Relative scale uncertainties (%s %s)",sOrder.c_str(),PDFFiles[iPDF].c_str());
                snprintf(titlel, sizeof(titlel), "-dsigma_6P/sigma");
                snprintf(titleu, sizeof(titleu), "+dsigma_6P/sigma");
             }
@@ -517,11 +517,11 @@ int main(int argc, char** argv) {
                yell <<  buffer << endl;
                yell << _SSEPSC << endl;
                shout << "bin      cross_section           lower_uncertainty       upper_uncertainty" << endl;
-               yell << _SSEPSC << endl;
+               yell << _TSEPSC << endl;
                for ( unsigned int iobs=0;iobs<XsUnc.xs.size();iobs++ ) {
                   printf("%5.i      %#18.11E      %#18.11E      %#18.11E\n",iobs+1,XsUnc.xs[iobs],XsUnc.dxsl[iobs],XsUnc.dxsu[iobs]);
                }
-               yell << _SSEPSC << endl;
+               yell << _TSEPSC << endl;
             }
 
             //! --- Initialize dimension bin counter
diff --git a/v2.5/toolkit/src/fnlo-tk-statunc.cc b/v2.5/toolkit/src/fnlo-tk-statunc.cc
index 62d1eecc..d4a7d264 100644
--- a/v2.5/toolkit/src/fnlo-tk-statunc.cc
+++ b/v2.5/toolkit/src/fnlo-tk-statunc.cc
@@ -434,9 +434,9 @@ int main(int argc, char** argv) {
    shout << "fnlo-tk-statunc: Evaluating uncertainties" << endl;
    yell  << _CSEPSC << endl;
    yell  << _DSEPSC << endl;
-   shout << "Relative Statistical Uncertainties" << endl;
+   shout << "Relative statistical uncertainties (ST)" << endl;
    yell  << _SSEPSC << endl;
-   shout << "bin      cross section           lower uncertainty       upper uncertainty" << endl;
+   shout << "bin      cross_section           lower_uncertainty       upper_uncertainty" << endl;
    yell  << _TSEPSC << endl;
 
    vector < double > dxsu;
@@ -460,6 +460,13 @@ int main(int argc, char** argv) {
       warn["fnlo-tk-statunc"] << "No Rivet ID found in fastNLO Table, stopping here without YODA output!" << endl;
       exit(0);
    }
+
+   size_t i0 = RivetId.find("/");
+   if (i0 == std::string::npos) {
+      warn["fnlo-tk-yodaout"] << "Found malformed Rivet ID without \"/\": " << RivetId << " Stopping here without YODA output!" <<  endl;
+      exit(0);
+   }
+
    if ( NDim == 2 ) {
       for (size_t i = fnlo.GetRivetId().find("/"); i < fnlo.GetRivetId().size(); i++) {
          if (isupper(fnlo.GetRivetId()[i])) {                                      // Find capital letter
@@ -469,8 +476,8 @@ int main(int argc, char** argv) {
          }
       }
       if (capital_pos == 0) {
-         error["fnlo-tk-statunc"] << "Rivet ID found in fastNLO table does not indicate the 2-dim. histogram counter, aborted." << endl;
-         exit(1);
+         warn["fnlo-tk-statunc"] << "Rivet ID found in fastNLO table does not indicate the 2-dim. histogram counter, stopping here without YODA output!" << endl;
+         exit(0);
       }
    }
 
diff --git a/v2.5/toolkit/src/fnlo-tk-yodaout.cc b/v2.5/toolkit/src/fnlo-tk-yodaout.cc
index 669cfa64..e3efb205 100644
--- a/v2.5/toolkit/src/fnlo-tk-yodaout.cc
+++ b/v2.5/toolkit/src/fnlo-tk-yodaout.cc
@@ -484,7 +484,7 @@ int main(int argc, char** argv) {
       yell <<  buffer << endl;
       yell << _SSEPSC << endl;
       shout << "bin      cross_section           lower_uncertainty       upper_uncertainty" << endl;
-      yell << _SSEPSC << endl;
+      yell << _TSEPSC << endl;
    }
    vector < double > dxsu;
    vector < double > dxsl;
@@ -500,21 +500,27 @@ int main(int argc, char** argv) {
       }
    }
    if ( XsUnc.xs.size() ) {
-      yell << _SSEPSC << endl;
+      yell << _TSEPSC << endl;
    }
 
    //! --- Get RivetID
    //!     For 2+-dimensions determine running number in Rivet plot name by spotting the capital letter in "RIVET_ID=" in the fnlo table
    //!     For inverted order, the Id starts with "-" after "/".
-   size_t capital_pos  = 0;
    int    invert_order = 1;
+   size_t capital_pos  = 0;
    string RivetId = fnlo->GetRivetId();
    if (RivetId.empty()) {
-      warn["fnlo-tk-yodaout"] << "No Rivet ID found in fastNLO Table, no YODA formatted output possible, exiting!" << endl;
+      warn["fnlo-tk-yodaout"] << "No Rivet ID found in fastNLO Table, stopping here without YODA output!" << endl;
       exit(0);
    }
+
+   size_t i0 = RivetId.find("/");
+   if (i0 == std::string::npos) {
+      warn["fnlo-tk-yodaout"] << "Found malformed Rivet ID without \"/\": " << RivetId << " Stopping here without YODA output!" <<  endl;
+      exit(0);
+   }
+
    if ( NDim > 1 ) {
-      size_t i0 = RivetId.find("/");
       /// Check for inversion of histogram order
       if ( RivetId.substr(i0,2) == "/-" ) {
          invert_order = -1;
-- 
GitLab