diff --git a/v2.5/generators/nlojet++/interface/hadron/InclusiveNJetDiffs.cc b/v2.5/generators/nlojet++/interface/hadron/InclusiveNJetDiffs.cc
index e1eb191a98e7f7d5568bc54bd53639379afe0ce1..108c4e1fee6f99208e5ff2e34af8ceda3e2168c6 100644
--- a/v2.5/generators/nlojet++/interface/hadron/InclusiveNJetDiffs.cc
+++ b/v2.5/generators/nlojet++/interface/hadron/InclusiveNJetDiffs.cc
@@ -142,7 +142,7 @@ private:
    double obs2[3];
    vector<string> ScaleLabel; // Scale labels (Scale1: must be defined; Scale2: only for flex-scale tables)
    // enum to switch between implemented scale definitions (max. of 2 simultaneously)
-   enum Scales { PTMAX, PTJETMIN, PTJETAVE, PTJETMAX };
+   enum Scales { PTMAX, PTJETMIN, PTJETAVE, PTJETMAX, HTP, HTPHALF };
    Scales mudef[2];
    double mu[2];
    int jetalgo;               // Define 1st fastjet jet algorithm (no default, must be defined)
@@ -181,6 +181,11 @@ struct fNLOSorter {
    bool operator() (const lorentzvector<double> &a, const lorentzvector<double> &b) {return (a.perp() > b.perp());};
 };
 
+// sign function returning -1, 0, 1
+template <typename T> int sgn(T val) {
+    return (T(0) < val) - (val < T(0));
+}
+
 // --- fastNLO user: check and get steering parameters once and store into static vars
 static std::map < std::string, bool > SteeringPars;
 
@@ -247,34 +252,38 @@ void UserHHC::phys_output(const std::basic_string<char>& __file_name, unsigned l
       }
    }
    // scale descriptions
-   string label;
    SteeringPars["ScaleDescriptionScale1"] = ftable->TestParameterInSteering("ScaleDescriptionScale1");
+   ScaleLabel.resize(2);
    if ( SteeringPars["ScaleDescriptionScale1"] ) {
-      ftable->GetParameterFromSteering("ScaleDescriptionScale1",label);
-      ScaleLabel.push_back(label);
+      ftable->GetParameterFromSteering("ScaleDescriptionScale1",ScaleLabel[0]);
    } else {
       say::error["ScenarioCode"] << "No description of scale 1, aborted!" << endl;
       exit(1);
    }
    SteeringPars["ScaleDescriptionScale2"] = ftable->TestParameterInSteering("ScaleDescriptionScale2");
+   ScaleLabel[1] = "pT_max_[GeV]"; // default
    if ( SteeringPars["ScaleDescriptionScale2"] ) {
-      ftable->GetParameterFromSteering("ScaleDescriptionScale2",label);
-      ScaleLabel.push_back(label);
+      ftable->GetParameterFromSteering("ScaleDescriptionScale2",ScaleLabel[1]);
+   } else {
+      ScaleLabel[1] = "-";
+      say::warn["ScenarioCode"] << "No description of scale 2, flexible-scale tables not possible!" << endl;
    }
    // scale descriptions define the scales
-   lptmax = true;
+   lptmax = false;
    for ( unsigned int i = 0; i < ScaleLabel.size(); i++ ) {
       if ( ScaleLabel[i] == "pT_max_[GeV]" ) {
          mudef[i] = PTMAX;
+         lptmax = true;
       } else if ( ScaleLabel[i] == "pT_jet_min_[GeV]" ) {
          mudef[i] = PTJETMIN;
-         lptmax = false;
       } else if ( ScaleLabel[i] == "pT_jet_ave_[GeV]" ) {
          mudef[i] = PTJETAVE;
-         lptmax = false;
       } else if ( ScaleLabel[i] == "pT_jet_max_[GeV]" ) {
          mudef[i] = PTJETMAX;
-         lptmax = false;
+      } else if ( ScaleLabel[i] == "HT_part_[GeV]" ) {
+         mudef[i] = HTP;
+      } else if ( ScaleLabel[i] == "HT_part/2_[GeV]" ) {
+         mudef[i] = HTPHALF;
       } else {
          say::error["ScenarioCode"] << "Unknown scale, i.e. scale description, aborted!" << endl;
          say::error["ScenarioCode"] << "ScaleLabel[" << i << "] = " << ScaleLabel[i] << endl;
@@ -283,7 +292,7 @@ void UserHHC::phys_output(const std::basic_string<char>& __file_name, unsigned l
       }
    }
 
-   // definition of default jet algorithm and jet phase space limits (no defaults)
+   // definition of jet algorithm and jet phase space limits (no defaults)
    //
    // --- fastNLO user: set the jet algorithm and size via steering file
    // fastjet clustering jet algos: 0 = kT, 1 = CA, 2 = anti-kT
@@ -471,6 +480,17 @@ void UserHHC::userfunc(const event_hhc& p, const amplitude_hhc& amp) {
    //     Usually, pT and E are in GeV, but this may be changed.
    //     ATTENTION: Scales must always be in GeV!
 
+   // derive partonic HT scale
+   double htp = 0.;
+   unsigned int np = p.upper();
+   for (unsigned int i = 1; i <= np; i++) {
+      htp += p[i].perp();
+      // --- debug output
+      if ( say::debug.GetSpeak() ) {
+         say::debug["ScenarioCode"] << "partonic HT: parton # i, pt, htp: " << i << ", " << p[i].perp() << ", " << htp << endl;
+      }
+   }
+
    // apply the jet algorithm(s) to partonic 4-vector array p of NLOJet++
    pj  = jetclusfj(p);
    pj2 = pj;
@@ -733,6 +753,15 @@ void UserHHC::userfunc(const event_hhc& p, const amplitude_hhc& amp) {
                   //    cout << "MAX iPair.first = " << iPair.first << ", sclmin = " << sclmin[iPair.first] <<", sclave = " << sclave[iPair.first]/acount[iPair.first] << ", sclmax = " << sclmax[iPair.first] << endl;}
                   mu[i] = sclmax[iPair.first];
                   break;
+               case HTP :
+                  // partonic sum pT
+                  mu[i] = htp;
+                  say::debug["ScenarioCode"]  << "HTP scale values: " << i <<  " : mu[i] = " << mu[i] << endl;
+                  break;
+               case HTPHALF :
+                  // half partonic sum pT
+                  mu[i] = htp/2.;
+                  say::debug["ScenarioCode"]  << "HTPHALF scale values: " << i <<  " : mu[i] = " << mu[i] << endl;
                default :
                   say::error["ScenarioCode"] << "Scale not yet implemented, aborted!" << endl;
                   say::error["ScenarioCode"] << "ScaleLabel[" << i << "] = " << ScaleLabel[i] << endl;
diff --git a/v2.5/generators/nlojet++/interface/hadron/InclusiveNJetEvents.cc b/v2.5/generators/nlojet++/interface/hadron/InclusiveNJetEvents.cc
index 725d49c64773c553d1553df3608f7cd21adb3c93..607f1424b4e5a788ae81136f5b599272d0ff9bca 100644
--- a/v2.5/generators/nlojet++/interface/hadron/InclusiveNJetEvents.cc
+++ b/v2.5/generators/nlojet++/interface/hadron/InclusiveNJetEvents.cc
@@ -55,6 +55,8 @@ void inputfunc(unsigned int&, unsigned int&, unsigned int&);
 void psinput(phasespace_hhc *, double&);
 // --- fastNLO v2.2: interface to NLOJet++: user class
 user_base_hhc * userfunc();
+// --- dphi
+double dphi(double phi2, double phi1);
 
 //----- array of the symbols symbols -----
 extern "C"{
@@ -132,13 +134,13 @@ private:
    int NDim;                  // Dimensionality of distributions (no default, must be defined)
    vector<string> DimLabel;   // Dimension labels (no default, must be defined)
    // enum to switch between implemented observables (max. of 3 simultaneously)
-   enum Obs { YMAX, YSTAR, Y1ABSRAP, Y2SIGN, MJJGEV, MJJTEV, PT12GEV, CHIJJ, HTHALFGEV, PTMAXGEV, DPHI12 };
+   enum Obs { YMAX, YSTAR, Y1ABSRAP, Y2SIGN, MJJGEV, MJJTEV, PT12GEV, CHIJJ, HTGEV, HTHALFGEV, PTMAXGEV, DPHI12 };
    Obs obsdef[3];
    double obs[3];
    vector<string> ScaleLabel; // Scale labels (Scale1: no default, must be defined; Scale2: default is "pT_max_[GeV]")
    // enum to switch between implemented scale definitions (max. of 2 simultaneously)
    // (Njet > 1!)
-   enum Scales { PTMAX, PT12AVE, PT123AVE, MJJHALF, PTMAXEXPYSTAR, EXPYSTAR, HTHALF };
+   enum Scales { PTMAX, PT12AVE, PT123AVE, MJJHALF, PTMAXEXPYSTAR, EXPYSTAR, HT, HTHALF, HTP, HTPHALF };
    Scales mudef[2];
    double mu[2];
    int jetalgo;               // Define fastjet jet algorithm (no default, must be defined)
@@ -249,6 +251,8 @@ void UserHHC::phys_output(const std::basic_string<char>& __file_name, unsigned l
          obsdef[i] = PT12GEV;
       } else if ( DimLabel[i] == "Chi" ) {
          obsdef[i] = CHIJJ;
+      } else if ( DimLabel[i] == "HT_[GeV]" ) {
+         obsdef[i] = HTGEV;
       } else if ( DimLabel[i] == "HT/2_[GeV]" ) {
          obsdef[i] = HTHALFGEV;
       } else if ( DimLabel[i] == "pT_max_[GeV]" ) {
@@ -280,7 +284,7 @@ void UserHHC::phys_output(const std::basic_string<char>& __file_name, unsigned l
       say::warn["ScenarioCode"] << "No description of scale 2, flexible-scale tables not possible!" << endl;
    }
    // scale descriptions define the scales
-   for ( unsigned int i = 0; i < 2; i++ ) {
+   for ( unsigned int i = 0; i < ScaleLabel.size(); i++ ) {
       if ( ScaleLabel[i] == "pT_max_[GeV]" ) {
          mudef[i] = PTMAX;
       } else if ( ScaleLabel[i] == "<pT_1,2>_[GeV]" ) {
@@ -293,8 +297,14 @@ void UserHHC::phys_output(const std::basic_string<char>& __file_name, unsigned l
          mudef[i] = PTMAXEXPYSTAR;
       } else if ( ScaleLabel[i] == "exp(0.3*y_star)" ) {
          mudef[i] = EXPYSTAR;
+      } else if ( ScaleLabel[i] == "HT_[GeV]" ) {
+         mudef[i] = HT;
       } else if ( ScaleLabel[i] == "HT/2_[GeV]" ) {
          mudef[i] = HTHALF;
+      } else if ( ScaleLabel[i] == "HT_part_[GeV]" ) {
+         mudef[i] = HTP;
+      } else if ( ScaleLabel[i] == "HT_part/2_[GeV]" ) {
+         mudef[i] = HTPHALF;
       } else {
          say::error["ScenarioCode"] << "Unknown scale, i.e. scale description, aborted!" << endl;
          say::error["ScenarioCode"] << "ScaleLabel[" << i << "] = " << ScaleLabel[i] << endl;
@@ -494,6 +504,17 @@ void UserHHC::userfunc(const event_hhc& p, const amplitude_hhc& amp) {
    //     Usually, pT and E are in GeV, but this may be changed.
    //     ATTENTION: Scales must always be in GeV!
 
+   // derive partonic HT scale
+   double htp = 0.;
+   unsigned int np = p.upper();
+   for (unsigned int i = 1; i <= np; i++) {
+      htp += p[i].perp();
+      // --- debug output
+      if ( say::debug.GetSpeak() ) {
+         say::debug["ScenarioCode"] << "partonic HT: parton # i, pt, htp: " << i << ", " << p[i].perp() << ", " << htp << endl;
+      }
+   }
+
    // apply the jet algorithm to partonic 4-vector array p of NLOJet++
    pj = jetclusfj(p);
    unsigned int nj = pj.upper();
@@ -574,11 +595,11 @@ void UserHHC::userfunc(const event_hhc& p, const amplitude_hhc& amp) {
    // no. of central jets, pT of third central jet, and HT of central jets (i.e. inside ycjjmax)
    int Ncjet   = 0;
    //   double pT3c = 0;
-   double HT2  = 0.;
+   double ht  = 0.;
    for (unsigned int k = 1; k <= njet; k++) {
       if ( abs(pj[k].rapidity()) < ycjjmax ) {
          Ncjet++;
-         HT2 += pj[k].perp()/2.;
+         ht += pj[k].perp();
       }
       //      if ( Ncjet == 3 ) {pT3c = pj[k].perp();}
    }
@@ -618,9 +639,13 @@ void UserHHC::userfunc(const event_hhc& p, const amplitude_hhc& amp) {
          // dijet chi
          obs[i] = exp(abs(y1-y2));
          break;
+      case HTGEV :
+         // jet pT sum
+         obs[i] = ht;
+         break;
       case HTHALFGEV :
          // half of jet pT sum
-         obs[i] = HT2;
+         obs[i] = ht/2.;
          break;
       case PTMAXGEV :
          // leading jet pT
@@ -701,9 +726,23 @@ void UserHHC::userfunc(const event_hhc& p, const amplitude_hhc& amp) {
             // ATLAS definition, this works only as second scale choice!
             mu[i] = exp(0.3*ystar);
             break;
+         case HT :
+            // jet pT sum (must be in GeV!)
+            mu[i] = ht;
+            break;
          case HTHALF :
             // half of jet pT sum (must be in GeV!)
-            mu[i] = HT2;
+            mu[i] = ht/2.;
+            break;
+         case HTP :
+            // partonic sum pT
+            mu[i] = htp;
+            say::debug["ScenarioCode"]  << "HTP scale values: " << i <<  " : mu[i] = " << mu[i] << endl;
+            break;
+         case HTPHALF :
+            // half partonic sum pT
+            mu[i] = htp/2.;
+            say::debug["ScenarioCode"]  << "HTPHALF scale values: " << i <<  " : mu[i] = " << mu[i] << endl;
             break;
          default :
             say::error["ScenarioCode"] << "Scale not yet implemented, aborted!" << endl;
@@ -744,14 +783,12 @@ void UserHHC::userfunc(const event_hhc& p, const amplitude_hhc& amp) {
    } else {
       // --- event rejected
       if ( say::debug.GetSpeak() ) {
-         say::debug["ScenarioCode"]  << "----------------- Event/jet rejected! ------------------" << endl;
+         say::debug["ScenarioCode"]  << "----------------- Event rejected! ------------------" << endl;
          say::debug["ScenarioCode"]  << "====================  End of event  ====================" << endl;
       }
    }
 }
 
-
-
 //------ DON'T TOUCH THIS PART! ------
 
 // --- fastNLO v2.2: interface to NLOJet++: read steering file, set LO of selected process
@@ -833,6 +870,19 @@ void InitfNLO(const std::basic_string<char>& __file_name) {
    // --- set fastNLO filename according to NLOJet++ command line arguments
    string tabFilename = __file_name.c_str();
    tabFilename += ".tab";
+
+#ifdef HAVE_ZLIB
+   bool lgzip = true;
+#else
+   bool lgzip = false;
+#endif
+
+   if ( SteeringPars["OutputCompression"] ) {
+      ftable->GetParameterFromSteering("OutputCompression",lgzip);
+   }
+   //   string filename = fScenConsts.OutputFilename;
+   //   if ( lgzip ) tabFilename += ".gz";
+   tabFilename += ".gz";
    ftable->SetFilename(tabFilename);
 }
 
@@ -860,3 +910,14 @@ void UserHHC::end_of_event() {
       }
    }
 }
+
+// --- define function for azimuthal angular distance in [-Pi,Pi]
+double dphi(double phi2, double phi1) {
+   double delta_phi = phi2-phi1;
+   if (delta_phi > M_PI) {
+      delta_phi -= 2*M_PI;
+   } else if (delta_phi < -M_PI) {
+      delta_phi += 2*M_PI;
+   }
+   return delta_phi;
+}
diff --git a/v2.5/generators/nlojet++/interface/hadron/InclusiveNJetPairs.cc b/v2.5/generators/nlojet++/interface/hadron/InclusiveNJetPairs.cc
index c04c1330809b2cdb58e6f2c1ae0561a9a364cb40..9fe3caff430909f3d0510da2b7fe3cb558cf1bf4 100644
--- a/v2.5/generators/nlojet++/interface/hadron/InclusiveNJetPairs.cc
+++ b/v2.5/generators/nlojet++/interface/hadron/InclusiveNJetPairs.cc
@@ -1,5 +1,5 @@
 //
-// fastNLO v2.2 creator code for inclusive N jet pairs scenarios
+// fastNLO v2.2 creator code for inclusive N jet-pairs scenarios
 //
 // ============== fastNLO user: ========================================
 // To create your own scenario, it is recommended to take this
@@ -139,7 +139,8 @@ private:
    double obs[3];
    vector<string> ScaleLabel; // Scale labels (Scale1: no default, must be defined; Scale2: default is "pT_max_[GeV]")
    // enum to switch between implemented scale definitions (max. of 2 simultaneously)
-   enum Scales { PTMAX, PTJET };
+   // (Njet > 1!)
+   enum Scales { PTMAX, PTJET, HTP, HTPHALF };
    Scales mudef[2];
    double mu[2];
    int jetalgo;               // Define fastjet jet algorithm (no default, must be defined)
@@ -154,7 +155,7 @@ private:
    double djlkmax;            // Maximal distance of neighbouring jet (default is +DBL_MAX)
    bool ldphi;                // Switch to use jet distance dphi/true or dR/false (default is true)
    bool lunique;              // Switch between unique entry per jet pair (true) or multiple ones (false) (default is false)
-   int Njetmin;               // Minimal number of jets in phase space (default is 1)
+   int Njetmin;               // Minimal number of jets in phase space (default is 2)
    double obsmin[3];          // Minimum in observable in nth dimension (default derived from binning)
    double obsmax[3];          // Maximum in observable in nth dimension (default derived from binning)
 };
@@ -177,6 +178,11 @@ struct fNLOSorter {
    bool operator() (const lorentzvector<double> &a, const lorentzvector<double> &b) {return (a.perp() > b.perp());};
 };
 
+// sign function returning -1, 0, 1
+template <typename T> int sgn(T val) {
+    return (T(0) < val) - (val < T(0));
+}
+
 // --- fastNLO user: check and get steering parameters once and store into static vars
 static std::map < std::string, bool > SteeringPars;
 
@@ -260,11 +266,15 @@ void UserHHC::phys_output(const std::basic_string<char>& __file_name, unsigned l
       say::warn["ScenarioCode"] << "No description of scale 2, flexible-scale tables not possible!" << endl;
    }
    // scale descriptions define the scales
-   for ( unsigned int i = 0; i < 2; i++ ) {
+   for ( unsigned int i = 0; i < ScaleLabel.size(); i++ ) {
       if ( ScaleLabel[i] == "pT_max_[GeV]" ) {
          mudef[i] = PTMAX;
       } else if ( ScaleLabel[i] == "pT_jet_[GeV]" ) {
          mudef[i] = PTJET;
+      } else if ( ScaleLabel[i] == "HT_part_[GeV]" ) {
+         mudef[i] = HTP;
+      } else if ( ScaleLabel[i] == "HT_part/2_[GeV]" ) {
+         mudef[i] = HTPHALF;
       } else {
          say::error["ScenarioCode"] << "Unknown scale, i.e. scale description, aborted!" << endl;
          say::error["ScenarioCode"] << "ScaleLabel[" << i << "] = " << ScaleLabel[i] << endl;
@@ -351,14 +361,14 @@ void UserHHC::phys_output(const std::basic_string<char>& __file_name, unsigned l
       say::error["ScenarioCode"] << "If you really want to mix, the code needs to be adapted." << endl;
       exit(1);
    }
-   // minimal number of jets required to be within preselected jet phase space (for inclusive jets this must be one!)
+   // minimal number of jets required to be within preselected jet phase space (for jet pairs this must be two!)
    SteeringPars["Njetmin"] = ftable->TestParameterInSteering("Njetmin");
-   Njetmin = 1;
+   Njetmin = 2;
    if ( SteeringPars["Njetmin"] ) {
       ftable->GetParameterFromSteering("Njetmin",Njetmin);
    }
-   if ( Njetmin < 1 ) {
-      say::error["ScenarioCode"] << "This is a 1+-jet scenario. At least one jet must be present, aborted!" << endl;
+   if ( Njetmin < 2 ) {
+      say::error["ScenarioCode"] << "This is a 2+-jet scenario. At least two jets must be present, aborted!" << endl;
       say::error["ScenarioCode"] << "Please correct the Njetmin requirement. Njetmin = " << Njetmin << endl;
       exit(1);
    }
@@ -457,6 +467,17 @@ void UserHHC::userfunc(const event_hhc& p, const amplitude_hhc& amp) {
    //     Usually, pT and E are in GeV, but this may be changed.
    //     ATTENTION: Scales must always be in GeV!
 
+   // derive partonic HT scale
+   double htp = 0.;
+   unsigned int np = p.upper();
+   for (unsigned int i = 1; i <= np; i++) {
+      htp += p[i].perp();
+      // --- debug output
+      if ( say::debug.GetSpeak() ) {
+         say::debug["ScenarioCode"] << "partonic HT: parton # i, pt, htp: " << i << ", " << p[i].perp() << ", " << htp << endl;
+      }
+   }
+
    // apply the jet algorithm to partonic 4-vector array p of NLOJet++
    pj = jetclusfj(p);
    unsigned int nj = pj.upper();
@@ -559,6 +580,16 @@ void UserHHC::userfunc(const event_hhc& p, const amplitude_hhc& amp) {
             // jet pT
             mu[i] = pj[k].perp();
             break;
+         case HTP :
+            // partonic sum pT
+            mu[i] = htp;
+            say::debug["ScenarioCode"]  << "HTP scale values: " << i <<  " : mu[i] = " << mu[i] << endl;
+            break;
+         case HTPHALF :
+            // half partonic sum pT
+            mu[i] = htp/2.;
+            say::debug["ScenarioCode"]  << "HTPHALF scale values: " << i <<  " : mu[i] = " << mu[i] << endl;
+            break;
          default :
             say::error["ScenarioCode"] << "Scale not yet implemented, aborted!" << endl;
             say::error["ScenarioCode"] << "ScaleLabel[" << i << "] = " << ScaleLabel[i] << endl;
@@ -640,7 +671,7 @@ void UserHHC::userfunc(const event_hhc& p, const amplitude_hhc& amp) {
                ftable->FillAllSubprocesses(contribsfix,scen);
             }
          } else {
-            // --- jet-pair observable rejected
+            // --- jet pair rejected
             if ( say::debug.GetSpeak() ) {
                say::debug["ScenarioCode"]  << "----------------- Jet-pair observable rejected! ------------------" << endl;
             }
diff --git a/v2.5/generators/nlojet++/interface/hadron/InclusiveNJets.cc b/v2.5/generators/nlojet++/interface/hadron/InclusiveNJets.cc
index 6f361e2d5d3c8092db05e794ed6601b352ce98c0..2db3cc9d67d09b655a58a16ee97906744692b2da 100644
--- a/v2.5/generators/nlojet++/interface/hadron/InclusiveNJets.cc
+++ b/v2.5/generators/nlojet++/interface/hadron/InclusiveNJets.cc
@@ -55,6 +55,8 @@ void inputfunc(unsigned int&, unsigned int&, unsigned int&);
 void psinput(phasespace_hhc *, double&);
 // --- fastNLO v2.2: interface to NLOJet++: user class
 user_base_hhc * userfunc();
+// --- dphi
+double dphi(double phi2, double phi1);
 
 //----- array of the symbols symbols -----
 extern "C"{
@@ -137,7 +139,7 @@ private:
    double obs[3];
    vector<string> ScaleLabel; // Scale labels (Scale1: no default, must be defined; Scale2: default is "pT_max_[GeV]")
    // enum to switch between implemented scale definitions (max. of 2 simultaneously)
-   enum Scales { PTMAX, PTJET };
+   enum Scales { PTMAX, PTJET, HTP, HTPHALF };
    Scales mudef[2];
    double mu[2];
    int jetalgo;               // Define fastjet jet algorithm (no default, must be defined)
@@ -170,6 +172,11 @@ struct fNLOSorter {
    bool operator() (const lorentzvector<double> &a, const lorentzvector<double> &b) {return (a.perp() > b.perp());};
 };
 
+// sign function returning -1, 0, 1
+template <typename T> int sgn(T val) {
+    return (T(0) < val) - (val < T(0));
+}
+
 // --- fastNLO user: check and get steering parameters once and store into static vars
 static std::map < std::string, bool > SteeringPars;
 
@@ -253,11 +260,15 @@ void UserHHC::phys_output(const std::basic_string<char>& __file_name, unsigned l
       say::warn["ScenarioCode"] << "No description of scale 2, flexible-scale tables not possible!" << endl;
    }
    // scale descriptions define the scales
-   for ( unsigned int i = 0; i < 2; i++ ) {
+   for ( unsigned int i = 0; i < ScaleLabel.size(); i++ ) {
       if ( ScaleLabel[i] == "pT_max_[GeV]" ) {
          mudef[i] = PTMAX;
       } else if ( ScaleLabel[i] == "pT_jet_[GeV]" ) {
          mudef[i] = PTJET;
+      } else if ( ScaleLabel[i] == "HT_part_[GeV]" ) {
+         mudef[i] = HTP;
+      } else if ( ScaleLabel[i] == "HT_part/2_[GeV]" ) {
+         mudef[i] = HTPHALF;
       } else {
          say::error["ScenarioCode"] << "Unknown scale, i.e. scale description, aborted!" << endl;
          say::error["ScenarioCode"] << "ScaleLabel[" << i << "] = " << ScaleLabel[i] << endl;
@@ -415,6 +426,17 @@ void UserHHC::userfunc(const event_hhc& p, const amplitude_hhc& amp) {
    //     Usually, pT and E are in GeV, but this may be changed.
    //     ATTENTION: Scales must always be in GeV!
 
+   // derive partonic HT scale
+   double htp = 0.;
+   unsigned int np = p.upper();
+   for (unsigned int i = 1; i <= np; i++) {
+      htp += p[i].perp();
+      // --- debug output
+      if ( say::debug.GetSpeak() ) {
+         say::debug["ScenarioCode"] << "partonic HT: parton # i, pt, htp: " << i << ", " << p[i].perp() << ", " << htp << endl;
+      }
+   }
+
    // apply the jet algorithm to partonic 4-vector array p of NLOJet++
    pj = jetclusfj(p);
    unsigned int nj = pj.upper();
@@ -477,16 +499,6 @@ void UserHHC::userfunc(const event_hhc& p, const amplitude_hhc& amp) {
    double ptmax = pj[1].perp();
    for (unsigned int k = 1; k <= njet; k++) {
 
-      // derive some jet quantities
-      //      double pt = pj[k].perp();
-      // double yeta;
-      // if ( ! lpseudo ) {
-      //    yeta = abs(pj[k].rapidity());
-      // } else {
-      //    yeta = abs(pj[k].prapidity());
-      // }
-      // double phi = atan2(pj[k].Y(), pj[k].X());
-
       // --- calculate observable of nth dimension
       for ( int i = 0; i<NDim; i++ ) {
          switch(obsdef[i]) {
@@ -545,6 +557,16 @@ void UserHHC::userfunc(const event_hhc& p, const amplitude_hhc& amp) {
                // jet pT
                mu[i] = pj[k].perp();
                break;
+         case HTP :
+            // partonic sum pT
+            mu[i] = htp;
+            say::debug["ScenarioCode"]  << "HTP scale values: " << i <<  " : mu[i] = " << mu[i] << endl;
+            break;
+         case HTPHALF :
+            // half partonic sum pT
+            mu[i] = htp/2.;
+            say::debug["ScenarioCode"]  << "HTPHALF scale values: " << i <<  " : mu[i] = " << mu[i] << endl;
+            break;
             default :
                say::error["ScenarioCode"] << "Scale not yet implemented, aborted!" << endl;
                say::error["ScenarioCode"] << "ScaleLabel[" << i << "] = " << ScaleLabel[i] << endl;
@@ -591,8 +613,6 @@ void UserHHC::userfunc(const event_hhc& p, const amplitude_hhc& amp) {
    }
 }
 
-
-
 //------ DON'T TOUCH THIS PART! ------
 
 // --- fastNLO v2.2: interface to NLOJet++: read steering file, set LO of selected process
@@ -714,3 +734,14 @@ void UserHHC::end_of_event() {
       }
    }
 }
+
+// --- define function for azimuthal angular distance in [-Pi,Pi]
+double dphi(double phi2, double phi1) {
+   double delta_phi = phi2-phi1;
+   if (delta_phi > M_PI) {
+      delta_phi -= 2*M_PI;
+   } else if (delta_phi < -M_PI) {
+      delta_phi += 2*M_PI;
+   }
+   return delta_phi;
+}
diff --git a/v2.5/generators/nlojet++/interface/hadron/InclusiveNJets_new.cc b/v2.5/generators/nlojet++/interface/hadron/InclusiveNJets_new.cc
index e242f1cd02cfb3c2e2dda6b96938b9017e36eb81..f12c9a2562abbf4fd1f76ed785ae878f4e19e272 100644
--- a/v2.5/generators/nlojet++/interface/hadron/InclusiveNJets_new.cc
+++ b/v2.5/generators/nlojet++/interface/hadron/InclusiveNJets_new.cc
@@ -59,446 +59,503 @@ using namespace std;
 // --- fastNLO v2.2: define user class to be used with NLOJet++
 class UserHHC : public FastNLOUserHHC {
 public:
-    // --- fastNLO user: evaluate steering file and define physics output (called once before first event)
-    virtual void read_steering();
-    // --- fastNLO user: analyze parton event (called once for each event)
-    virtual void userfunc(const event_hhc&, const amplitude_hhc&);
+   // --- fastNLO user: evaluate steering file and define physics output (called once before first event)
+   virtual void read_steering();
+   // --- fastNLO user: analyze parton event (called once for each event)
+   virtual void userfunc(const event_hhc&, const amplitude_hhc&);
 
 private:
-    // --- fastNLO user: define the jet algorithm (for the choice of included header file above)
-    fastjet_jets jetclusfj;
-
-    // --- define the jet structure
-    bounded_vector<lorentzvector<double> > pj;
-
-    // --- fastNLO steering
-    bool lFlexibleScaleTable; // Fill fixed- or flexible-scale table
-    int NDim; // Dimensionality of distributions
-    vector<string> DimLabel; // Dimension labels
-    // enum to switch between implemented observables (max. of 3 simultaneously)
-    enum Obs { PTJETGEV,
-               YJET,
-               PHIJET };
-    Obs obsdef[3];
-    double obs[3];
-    vector<string> ScaleLabel; // Scale labels
-    // enum to switch between implemented scale definitions (max. of 2 simultaneously)
-    // (Njet > 1!)
-    enum Scales { PTMAX,
-                  PTJET };
-    Scales mudef[2];
-    double mu[2];
-    int jetalgo; // Define fastjet jet algorithm
-    double jetsize; // Define jet size R
-    double overlapthreshold; // Define overlap threshold (for some jet algorithms)
-    double ptjmin; // Minimal jet pT (should be >= minimum of 1 GeV specified in interface to fastjet)
-    double yetajmin; // Minimal jet (pseudo-)rapidity
-    double yetajmax; // Maximal jet (pseudo-)rapidity
-    bool lpseudo; // Switch to use either jet rapidity y or jet eta
-    int Njetmin; // Minimal number of jets in phase space
-    double obsmin[3]; // Minimum in observable in nth dimension (default derived from binning)
-    double obsmax[3]; // Maximum in observable in nth dimension (default derived from binning)
+   // --- fastNLO user: define the jet algorithm (for the choice of included header file above)
+   fastjet_jets jetclusfj;
+
+   // --- define the jet structure
+   bounded_vector<lorentzvector<double> > pj;
+
+   // --- fastNLO steering
+   bool lFlexibleScaleTable;  // Fill fixed- or flexible-scale table (default is fixed-scale)
+   int NDim;                  // Dimensionality of distributions (no default, must be defined)
+   vector<string> DimLabel;   // Dimension labels (no default, must be defined)
+   // enum to switch between implemented observables (max. of 3 simultaneously)
+   enum Obs { PTJETGEV, YJET, ETAJET, PHIJET };
+   Obs obsdef[3];
+   double obs[3];
+   vector<string> ScaleLabel; // Scale labels (Scale1: no default, must be defined; Scale2: default is "pT_max_[GeV]")
+   // enum to switch between implemented scale definitions (max. of 2 simultaneously)
+   enum Scales { PTMAX, PTJET, HTP, HTPHALF };
+   Scales mudef[2];
+   double mu[2];
+   int jetalgo;               // Define fastjet jet algorithm (no default, must be defined)
+   double jetsize;            // Define jet size R (no default, must be defined)
+   double overlapthreshold;   // Define overlap threshold (default is 0.5)
+   double ptjmin;             // Minimal jet pT (no default, must be defined; should be >= minimum of 1 GeV specified in interface to fastjet)
+   double yetajmin;           // Minimal jet (pseudo-)rapidity (no default, must be defined)
+   double yetajmax;           // Maximal jet (pseudo-)rapidity (no default, must be defined)
+   bool lpseudo;              // Switch to use either jet rapidity y or jet eta
+   int Njetmin;               // Minimal number of jets in phase space (default is 1)
+   double obsmin[3];          // Minimum in observable in nth dimension (default derived from binning)
+   double obsmax[3];          // Maximum in observable in nth dimension (default derived from binning)
 };
 
 // --- fastNLO user: modify the jet selection in UserHHC::userfunc (default = cutting in |y| min, |y| max and pt min)
 //                   (the return value must be true for jets to be UNselected)
 struct fNLOSelector {
-    fNLOSelector(double ymin, double ymax, double ptmin, bool pseudo = false)
-        : _ymin(ymin)
-        , _ymax(ymax)
-        , _ptmin(ptmin)
-        , _pseudo(pseudo){};
-    double _ymin, _ymax, _ptmin;
-    bool _pseudo;
-    bool operator()(const lorentzvector<double>& a)
-    {
-        if (!_pseudo)
-            return !(_ymin <= abs(a.rapidity()) && abs(a.rapidity()) < _ymax && _ptmin <= a.perp());
-        else
-            return !(_ymin <= abs(a.prapidity()) && abs(a.prapidity()) < _ymax && _ptmin <= a.perp());
-    };
+   fNLOSelector(double ymin, double ymax, double ptmin, bool pseudo=false):
+      _ymin (ymin), _ymax (ymax), _ptmin (ptmin), _pseudo (pseudo){};
+   double _ymin, _ymax, _ptmin;
+   bool _pseudo;
+   bool operator() (const lorentzvector<double> &a) {
+      if (!_pseudo) return ! (_ymin <= abs(a.rapidity())  && abs(a.rapidity())  < _ymax && _ptmin <= a.perp());
+      else          return ! (_ymin <= abs(a.prapidity()) && abs(a.prapidity()) < _ymax && _ptmin <= a.perp());
+   };
 };
 
 FastNLOUserHHC* FastNLOUserHHC::instance = new UserHHC;
 
 // --- fastNLO user: modify the jet sorting in UserHHC::userfunc (default = descending in jet pt)
 struct fNLOSorter {
-    bool operator()(const lorentzvector<double>& a, const lorentzvector<double>& b) { return (a.perp() > b.perp()); };
+   bool operator() (const lorentzvector<double> &a, const lorentzvector<double> &b) {return (a.perp() > b.perp());};
 };
 
+// sign function returning -1, 0, 1
+template <typename T> int sgn(T val) {
+    return (T(0) < val) - (val < T(0));
+}
+
 // --- fastNLO user: check and get steering parameters once and store into static vars
-static std::map<std::string, bool> SteeringPars;
+static std::map < std::string, bool > SteeringPars;
 
 // --- fastNLO user: class UserHHC: evaluate steering file and define physics output (called once before first event)
 void UserHHC::read_steering()
 {
-    // --- fastNLO user:
-    //     Here is your playground where you can evaluate the steering file settings.
-    //     ATTENTION: Some settings are mandatory for the correct functioning!
-
-    // get general steering parameters needed here from steering file
-    say::debug["UserHHC::phys_output"] << "Evaluating steering parameters ..." << endl;
-    // fixed- or flexible-scale table
-    SteeringPars["FlexibleScaleTable"] = ftable->TestParameterInSteering("FlexibleScaleTable");
-    lFlexibleScaleTable = false; // default
-    if (SteeringPars["FlexibleScaleTable"]) {
-        ftable->GetParameterFromSteering("FlexibleScaleTable", lFlexibleScaleTable);
-    }
-    // dimensionality
-    SteeringPars["DifferentialDimension"] = ftable->TestParameterInSteering("DifferentialDimension");
-    if (SteeringPars["DifferentialDimension"]) {
-        ftable->GetParameterFromSteering("DifferentialDimension", NDim);
-    } else {
-        say::error["ScenarioCode"] << "Dimensioning of binning not set, aborted!" << endl;
-        exit(1);
-    }
-    if (NDim < 1 || 3 < NDim) {
-        say::error["ScenarioCode"] << "Only 1- to 3-dimensional binning implemented, aborted!" << endl;
-        say::error["ScenarioCode"] << "Please implement the requested " << NDim << "-dimensional binning." << endl;
-        exit(1);
-    }
-    // dimension labels
-    SteeringPars["DimensionLabels"] = ftable->TestParameterInSteering("DimensionLabels");
-    DimLabel.resize(NDim);
-    if (SteeringPars["DimensionLabels"]) {
-        ftable->GetParameterFromSteering("DimensionLabels", DimLabel);
-    } else {
-        say::error["ScenarioCode"] << "Dimension labels not set, aborted!" << endl;
-        exit(1);
-    }
-    // define the observables according to the dimension labels
-    for (int i = 0; i < NDim; i++) {
-        if (DimLabel[i] == "|y|") {
-            obsdef[i] = YJET;
-        } else if (DimLabel[i] == "pT_[GeV]") {
-            obsdef[i] = PTJETGEV;
-        } else if (DimLabel[i] == "phi") {
-            obsdef[i] = PHIJET;
-        } else {
-            say::error["ScenarioCode"] << "Unknown observable, i.e. dimension label, aborted!" << endl;
-            say::error["ScenarioCode"] << "DimLabel[" << i << "] = " << DimLabel[i] << endl;
-            say::error["ScenarioCode"] << "Please complement this scenario to include the requested observable." << endl;
-            exit(1);
-        }
-    }
-    // scale descriptions
-    SteeringPars["ScaleDescriptionScale1"] = ftable->TestParameterInSteering("ScaleDescriptionScale1");
-    ScaleLabel.resize(2);
-    if (SteeringPars["ScaleDescriptionScale1"]) {
-        ftable->GetParameterFromSteering("ScaleDescriptionScale1", ScaleLabel[0]);
-    } else {
-        say::error["ScenarioCode"] << "No description of scale 1, aborted!" << endl;
-        exit(1);
-    }
-    SteeringPars["ScaleDescriptionScale2"] = ftable->TestParameterInSteering("ScaleDescriptionScale2");
-    ScaleLabel[1] = "pT_max_[GeV]"; // default
-    if (SteeringPars["ScaleDescriptionScale2"]) {
-        ftable->GetParameterFromSteering("ScaleDescriptionScale2", ScaleLabel[1]);
-    } else {
-        ScaleLabel[1] = "-";
-        say::warn["ScenarioCode"] << "No description of scale 2, flexible-scale tables not possible!" << endl;
-    }
-    // scale descriptions define the scales
-    for (unsigned int i = 0; i < 2; i++) {
-        if (ScaleLabel[i] == "pT_max_[GeV]") {
-            mudef[i] = PTMAX;
-        } else if (ScaleLabel[i] == "pT_jet_[GeV]") {
-            mudef[i] = PTJET;
-        } else {
-            say::error["ScenarioCode"] << "Unknown scale, i.e. scale description, aborted!" << endl;
-            say::error["ScenarioCode"] << "ScaleLabel[" << i << "] = " << ScaleLabel[i] << endl;
-            say::error["ScenarioCode"] << "Please complement this scenario to include the requested scale." << endl;
-            exit(1);
-        }
-    }
-
-    // definition of jet algorithm and jet phase space limits (no defaults)
-    //
-    // --- fastNLO user: set the jet algorithm and size via steering file
-    // fastjet clustering jet algos: 0 = kT, 1 = CA, 2 = anti-kT
-    // fastjet cone jet algos: 10 = SISCone, 11 = CDFMidPointCone, 12 = D0RunIICone
-    SteeringPars["JetAlgo"] = ftable->TestParameterInSteering("JetAlgo");
-    if (SteeringPars["JetAlgo"]) {
-        ftable->GetParameterFromSteering("JetAlgo", jetalgo);
-        if (jetalgo < 0 || (2 < jetalgo && jetalgo < 10) || 12 < jetalgo) {
-            say::error["ScenarioCode"] << "Unknown jet algorithm " << jetalgo << ", aborted!" << endl;
-            exit(1);
-        }
-    } else {
-        say::error["ScenarioCode"] << "No jet algorithm selected, aborted!" << endl;
-        exit(1);
-    }
-    SteeringPars["Rjet"] = ftable->TestParameterInSteering("Rjet");
-    if (SteeringPars["Rjet"]) {
-        ftable->GetParameterFromSteering("Rjet", jetsize);
-    } else {
-        say::error["ScenarioCode"] << "Jet size R not defined, aborted!" << endl;
-        exit(1);
-    }
-    SteeringPars["OvThr"] = ftable->TestParameterInSteering("OvThr");
-    overlapthreshold = 0.5; // default
-    if (SteeringPars["OvThr"]) {
-        ftable->GetParameterFromSteering("OvThr", overlapthreshold);
-    } else if (jetalgo > 9) {
-        say::error["ScenarioCode"] << "Overlap threshold not defined for jet algorithm " << jetalgo << ", aborted!" << endl;
-        exit(1);
-    }
-    // --- fastNLO user: declare and initialize overall jet phase space cuts via steering file
-    // overall lowest pT for jets to be considered
-    SteeringPars["ptjmin"] = ftable->TestParameterInSteering("ptjmin");
-    if (SteeringPars["ptjmin"]) {
-        ftable->GetParameterFromSteering("ptjmin", ptjmin);
-    } else {
-        say::error["ScenarioCode"] << "Minimal jet pT (ptjmin) not defined, aborted!" << endl;
-        exit(1);
-    }
-    // overall highest pT for jets not implemented, since uncritical with respect to CPU time consumption
-    // overall smallest |(pseudo-)rapidity| for jets to be considered, use either y or eta but not both
-    SteeringPars["yjmin"] = ftable->TestParameterInSteering("yjmin");
-    SteeringPars["etajmin"] = ftable->TestParameterInSteering("etajmin");
-    if (SteeringPars["yjmin"] && !SteeringPars["etajmin"]) {
-        ftable->GetParameterFromSteering("yjmin", yetajmin);
-    } else if (!SteeringPars["yjmin"] && SteeringPars["etajmin"]) {
-        ftable->GetParameterFromSteering("etajmin", yetajmin);
-    } else {
-        say::error["ScenarioCode"] << "Minimal jet (pseudo)rapidity (yjmin or etajmin) not uniquely defined, aborted!" << endl;
-        exit(1);
-    }
-    // overall largest |(pseudo-)rapidity| for jets to be considered, use either y or eta but not both
-    SteeringPars["yjmax"] = ftable->TestParameterInSteering("yjmax");
-    SteeringPars["etajmax"] = ftable->TestParameterInSteering("etajmax");
-    if (SteeringPars["yjmax"] && !SteeringPars["etajmax"]) {
-        ftable->GetParameterFromSteering("yjmax", yetajmax);
-    } else if (!SteeringPars["yjmax"] && SteeringPars["etajmax"]) {
-        ftable->GetParameterFromSteering("etajmax", yetajmax);
-    } else {
-        say::error["ScenarioCode"] << "Maximal jet (pseudo)rapidity (yjmax or etajmax) not uniquely defined, aborted!" << endl;
-        exit(1);
-    }
-    // define logical for decision on cuts in (pseudo-)rapidity, no mixing allowed here
-    if (SteeringPars["yjmin"] && SteeringPars["yjmax"]) {
-        lpseudo = false;
-    } else if (SteeringPars["etajmin"] && SteeringPars["etajmax"]) {
-        lpseudo = true;
-    } else {
-        say::error["ScenarioCode"] << "Phase space cuts mixed in (pseudo-)rapidity, aborted!" << endl;
-        say::error["ScenarioCode"] << "Booleans for cut selections are"
-                                   << " yjmin " << SteeringPars["yjmin"] << ", yjmax " << SteeringPars["yjmax"] << ", etajmin " << SteeringPars["etajmin"] << ", etajmax " << SteeringPars["etajmax"] << endl;
-        say::error["ScenarioCode"] << "If you really want to mix, the code needs to be adapted." << endl;
-        exit(1);
-    }
-    // minimal number of jets required to be within preselected jet phase space (for inclusive jets this must be one!)
-    SteeringPars["Njetmin"] = ftable->TestParameterInSteering("Njetmin");
-    Njetmin = 1;
-    if (SteeringPars["Njetmin"]) {
-        ftable->GetParameterFromSteering("Njetmin", Njetmin);
-    }
-
-    // --- fastNLO user: declare and initialize Njet phase space cuts and definitions via steering file
-    // overall minimum & maximum for 1st observable, e.g. maximal absolute rapidity |y_max|
-    SteeringPars["obs0min"] = ftable->TestParameterInSteering("obs0min");
-    obsmin[0] = ftable->GetObsBinsLoBoundsMin(0); // by default derived from binning in obs0
-    if (SteeringPars["obs0min"]) {
-        ftable->GetParameterFromSteering("obs0min", obsmin[0]);
-    }
-    SteeringPars["obs0max"] = ftable->TestParameterInSteering("obs0max");
-    obsmax[0] = ftable->GetObsBinsUpBoundsMax(0); // by default derived from binning in obs0
-    if (SteeringPars["obs0max"]) {
-        ftable->GetParameterFromSteering("obs0max", obsmax[0]);
-    }
-    // overall minimum & maximum for 2nd observable, e.g. dijet mass mjj
-    obsmin[1] = -DBL_MAX;
-    obsmax[1] = +DBL_MAX;
-    if (NDim > 1) {
-        SteeringPars["obs1min"] = ftable->TestParameterInSteering("obs1min");
-        obsmin[1] = ftable->GetObsBinsLoBoundsMin(1); // by default derived from binning in obs1
-        if (SteeringPars["obs1min"]) {
-            ftable->GetParameterFromSteering("obs1min", obsmin[1]);
-        }
-        SteeringPars["obs1max"] = ftable->TestParameterInSteering("obs1max");
-        obsmax[1] = ftable->GetObsBinsUpBoundsMax(1); // by default derived from binning in obs1
-        if (SteeringPars["obs1max"]) {
-            ftable->GetParameterFromSteering("obs1max", obsmax[1]);
-        }
-    }
-    // overall minimum & maximum for 3rd observable
-    obsmin[2] = -DBL_MAX;
-    obsmax[2] = +DBL_MAX;
-    if (NDim > 2) {
-        SteeringPars["obs2min"] = ftable->TestParameterInSteering("obs2min");
-        obsmin[2] = ftable->GetObsBinsLoBoundsMin(2); // by default derived from binning in obs2
-        if (SteeringPars["obs2min"]) {
-            ftable->GetParameterFromSteering("obs2min", obsmin[2]);
-        }
-        SteeringPars["obs2max"] = ftable->TestParameterInSteering("obs2max");
-        obsmax[2] = ftable->GetObsBinsUpBoundsMax(2); // by default derived from binning in obs2
-        if (SteeringPars["obs2max"]) {
-            ftable->GetParameterFromSteering("obs2max", obsmax[2]);
-        }
-    }
-    jetclusfj.setup(static_cast<fastjet_jets::JetAlgorithm>(jetalgo), jetsize, overlapthreshold);
+   // --- fastNLO user:
+   //     Here is your playground where you can evaluate the steering file settings.
+   //     ATTENTION: Some settings are mandatory for the correct functioning!
+
+   // get general steering parameters needed here from steering file
+   say::debug["UserHHC::phys_output"] << "Evaluating steering parameters ..." << endl;
+   // fixed- or flexible-scale table
+   SteeringPars["FlexibleScaleTable"] = ftable->TestParameterInSteering("FlexibleScaleTable");
+   lFlexibleScaleTable = false; // default
+   if ( SteeringPars["FlexibleScaleTable"] ) {
+      ftable->GetParameterFromSteering("FlexibleScaleTable",lFlexibleScaleTable);
+   }
+   // dimensionality
+   SteeringPars["DifferentialDimension"] = ftable->TestParameterInSteering("DifferentialDimension");
+   if ( SteeringPars["DifferentialDimension"] ) {
+      ftable->GetParameterFromSteering("DifferentialDimension",NDim);
+   } else {
+      say::error["ScenarioCode"] << "Dimensioning of binning not set, aborted!" << endl;
+      exit(1);
+   }
+   if ( NDim < 1 || 3 < NDim ) {
+      say::error["ScenarioCode"] << "Only 1- to 3-dimensional binning implemented, aborted!" << endl;
+      say::error["ScenarioCode"] << "Please implement the requested " << NDim << "-dimensional binning." << endl;
+      exit(1);
+   }
+   // dimension labels
+   SteeringPars["DimensionLabels"] = ftable->TestParameterInSteering("DimensionLabels");
+   DimLabel.resize(NDim);
+   if ( SteeringPars["DimensionLabels"] ) {
+      ftable->GetParameterFromSteering("DimensionLabels",DimLabel);
+   } else {
+      say::error["ScenarioCode"] << "Dimension labels not set, aborted!" << endl;
+      exit(1);
+   }
+   // define the observables according to the dimension labels
+   for ( int i = 0; i<NDim; i++ ) {
+      if ( DimLabel[i]        == "|y|" ) {
+         obsdef[i] = YJET;
+      } else if ( DimLabel[i] == "|eta|" ) {
+         obsdef[i] = ETAJET;
+      } else if ( DimLabel[i] == "pT_[GeV]" ) {
+         obsdef[i] = PTJETGEV;
+      } else if ( DimLabel[i] == "phi" ) {
+         obsdef[i] = PHIJET;
+      } else {
+         say::error["ScenarioCode"] << "Unknown observable, i.e. dimension label, aborted!" << endl;
+         say::error["ScenarioCode"] << "DimLabel[" << i << "] = " << DimLabel[i] << endl;
+         say::error["ScenarioCode"] << "Please complement this scenario to include the requested observable." << endl;
+         exit(1);
+      }
+   }
+   // scale descriptions
+   SteeringPars["ScaleDescriptionScale1"] = ftable->TestParameterInSteering("ScaleDescriptionScale1");
+   ScaleLabel.resize(2);
+   if ( SteeringPars["ScaleDescriptionScale1"] ) {
+      ftable->GetParameterFromSteering("ScaleDescriptionScale1",ScaleLabel[0]);
+   } else {
+      say::error["ScenarioCode"] << "No description of scale 1, aborted!" << endl;
+      exit(1);
+   }
+   SteeringPars["ScaleDescriptionScale2"] = ftable->TestParameterInSteering("ScaleDescriptionScale2");
+   ScaleLabel[1] = "pT_max_[GeV]"; // default
+   if ( SteeringPars["ScaleDescriptionScale2"] ) {
+      ftable->GetParameterFromSteering("ScaleDescriptionScale2",ScaleLabel[1]);
+   } else {
+      ScaleLabel[1] = "-";
+      say::warn["ScenarioCode"] << "No description of scale 2, flexible-scale tables not possible!" << endl;
+   }
+   // scale descriptions define the scales
+   for ( unsigned int i = 0; i < ScaleLabel.size(); i++ ) {
+      if ( ScaleLabel[i] == "pT_max_[GeV]" ) {
+         mudef[i] = PTMAX;
+      } else if ( ScaleLabel[i] == "pT_jet_[GeV]" ) {
+         mudef[i] = PTJET;
+      } else if ( ScaleLabel[i] == "HT_part_[GeV]" ) {
+         mudef[i] = HTP;
+      } else if ( ScaleLabel[i] == "HT_part/2_[GeV]" ) {
+         mudef[i] = HTPHALF;
+      } else {
+         say::error["ScenarioCode"] << "Unknown scale, i.e. scale description, aborted!" << endl;
+         say::error["ScenarioCode"] << "ScaleLabel[" << i << "] = " << ScaleLabel[i] << endl;
+         say::error["ScenarioCode"] << "Please complement this scenario to include the requested scale." << endl;
+         exit(1);
+      }
+   }
+
+   // definition of jet algorithm and jet phase space limits (no defaults)
+   //
+   // --- fastNLO user: set the jet algorithm and size via steering file
+   // fastjet clustering jet algos: 0 = kT, 1 = CA, 2 = anti-kT
+   // fastjet cone jet algos: 10 = SISCone, 11 = CDFMidPointCone, 12 = D0RunIICone
+   SteeringPars["JetAlgo"] = ftable->TestParameterInSteering("JetAlgo");
+   if ( SteeringPars["JetAlgo"] ) {
+      ftable->GetParameterFromSteering("JetAlgo",jetalgo);
+      if ( jetalgo < 0 || (2 < jetalgo && jetalgo < 10) || 12 < jetalgo ) {
+         say::error["ScenarioCode"] << "Unknown jet algorithm " << jetalgo << ", aborted!" << endl;
+         exit(1);
+      }
+   } else {
+      say::error["ScenarioCode"] << "No jet algorithm selected, aborted!" << endl;
+      exit(1);
+   }
+   SteeringPars["Rjet"] = ftable->TestParameterInSteering("Rjet");
+   if ( SteeringPars["Rjet"] ) {
+      ftable->GetParameterFromSteering("Rjet",jetsize);
+   } else {
+      say::error["ScenarioCode"] << "Jet size R not defined, aborted!" << endl;
+      exit(1);
+   }
+   SteeringPars["OvThr"] = ftable->TestParameterInSteering("OvThr");
+   overlapthreshold = 0.5; // default
+   if ( SteeringPars["OvThr"] ) {
+      ftable->GetParameterFromSteering("OvThr",overlapthreshold);
+   } else if ( jetalgo > 9 ) {
+      say::error["ScenarioCode"] << "Overlap threshold not defined for jet algorithm " << jetalgo << ", aborted!" << endl;
+      exit(1);
+   }
+   // --- fastNLO user: declare and initialize overall jet phase space cuts via steering file
+   // overall lowest pT for jets to be considered
+   SteeringPars["ptjmin"] = ftable->TestParameterInSteering("ptjmin");
+   if ( SteeringPars["ptjmin"] ) {
+      ftable->GetParameterFromSteering("ptjmin",ptjmin);
+   } else {
+      say::error["ScenarioCode"] << "Minimal jet pT (ptjmin) not defined, aborted!" << endl;
+      exit(1);
+   }
+   // overall highest pT for jets not implemented, since uncritical with respect to CPU time consumption
+   // overall smallest |(pseudo-)rapidity| for jets to be considered, use either y or eta but not both
+   SteeringPars["yjmin"]   = ftable->TestParameterInSteering("yjmin");
+   SteeringPars["etajmin"] = ftable->TestParameterInSteering("etajmin");
+   if ( SteeringPars["yjmin"] && !SteeringPars["etajmin"] ) {
+      ftable->GetParameterFromSteering("yjmin",yetajmin);
+   } else if ( !SteeringPars["yjmin"] && SteeringPars["etajmin"] ) {
+      ftable->GetParameterFromSteering("etajmin",yetajmin);
+   } else {
+      say::error["ScenarioCode"] << "Minimal jet (pseudo)rapidity (yjmin or etajmin) not uniquely defined, aborted!" << endl;
+      exit(1);
+   }
+   // overall largest |(pseudo-)rapidity| for jets to be considered, use either y or eta but not both
+   SteeringPars["yjmax"]   = ftable->TestParameterInSteering("yjmax");
+   SteeringPars["etajmax"] = ftable->TestParameterInSteering("etajmax");
+   if ( SteeringPars["yjmax"] && !SteeringPars["etajmax"] ) {
+      ftable->GetParameterFromSteering("yjmax",yetajmax);
+   } else if ( !SteeringPars["yjmax"] && SteeringPars["etajmax"] ) {
+      ftable->GetParameterFromSteering("etajmax",yetajmax);
+   } else {
+      say::error["ScenarioCode"] << "Maximal jet (pseudo)rapidity (yjmax or etajmax) not uniquely defined, aborted!" << endl;
+      exit(1);
+   }
+   // define logical for decision on cuts in (pseudo-)rapidity, no mixing allowed here
+   if ( SteeringPars["yjmin"] && SteeringPars["yjmax"] ) {
+      lpseudo = false;
+   } else if ( SteeringPars["etajmin"] && SteeringPars["etajmax"] ) {
+      lpseudo = true;
+   } else {
+      say::error["ScenarioCode"] << "Phase space cuts mixed in (pseudo-)rapidity, aborted!" << endl;
+      say::error["ScenarioCode"] << "Booleans for cut selections are" <<
+         " yjmin "    << SteeringPars["yjmin"] <<
+         ", yjmax "   << SteeringPars["yjmax"] <<
+         ", etajmin " << SteeringPars["etajmin"] <<
+         ", etajmax " << SteeringPars["etajmax"] << endl;
+      say::error["ScenarioCode"] << "If you really want to mix, the code needs to be adapted." << endl;
+      exit(1);
+   }
+   // minimal number of jets required to be within preselected jet phase space (for inclusive jets this must be one!)
+   SteeringPars["Njetmin"] = ftable->TestParameterInSteering("Njetmin");
+   Njetmin = 1;
+   if ( SteeringPars["Njetmin"] ) {
+      ftable->GetParameterFromSteering("Njetmin",Njetmin);
+   }
+   if ( Njetmin < 1 ) {
+      say::error["ScenarioCode"] << "This is a 1+-jet scenario. At least one jet must be present, aborted!" << endl;
+      say::error["ScenarioCode"] << "Please correct the Njetmin requirement. Njetmin = " << Njetmin << endl;
+      exit(1);
+   }
+
+   // --- fastNLO user: declare and initialize Njet phase space cuts and definitions via steering file
+   // overall minimum & maximum for 1st observable, e.g. maximal absolute rapidity |y_max|
+   SteeringPars["obs0min"] = ftable->TestParameterInSteering("obs0min");
+   obsmin[0] = ftable->GetObsBinsLoBoundsMin(0); // by default derived from binning in obs0
+   if ( SteeringPars["obs0min"] ) {
+      ftable->GetParameterFromSteering("obs0min",obsmin[0]);
+   }
+   SteeringPars["obs0max"] = ftable->TestParameterInSteering("obs0max");
+   obsmax[0] = ftable->GetObsBinsUpBoundsMax(0); // by default derived from binning in obs0
+   if ( SteeringPars["obs0max"] ) {
+      ftable->GetParameterFromSteering("obs0max",obsmax[0]);
+   }
+   // overall minimum & maximum for 2nd observable, e.g. dijet mass mjj
+   obsmin[1] = -DBL_MAX;
+   obsmax[1] = +DBL_MAX;
+   if (NDim > 1) {
+      SteeringPars["obs1min"] = ftable->TestParameterInSteering("obs1min");
+      obsmin[1] = ftable->GetObsBinsLoBoundsMin(1); // by default derived from binning in obs1
+      if ( SteeringPars["obs1min"] ) {
+         ftable->GetParameterFromSteering("obs1min",obsmin[1]);
+      }
+      SteeringPars["obs1max"] = ftable->TestParameterInSteering("obs1max");
+      obsmax[1] = ftable->GetObsBinsUpBoundsMax(1); // by default derived from binning in obs1
+      if ( SteeringPars["obs1max"] ) {
+         ftable->GetParameterFromSteering("obs1max",obsmax[1]);
+      }
+   }
+   // overall minimum & maximum for 3rd observable
+   obsmin[2] = -DBL_MAX;
+   obsmax[2] = +DBL_MAX;
+   if (NDim > 2) {
+      SteeringPars["obs2min"] = ftable->TestParameterInSteering("obs2min");
+      obsmin[2] = ftable->GetObsBinsLoBoundsMin(2); // by default derived from binning in obs2
+      if ( SteeringPars["obs2min"] ) {
+         ftable->GetParameterFromSteering("obs2min",obsmin[2]);
+      }
+      SteeringPars["obs2max"] = ftable->TestParameterInSteering("obs2max");
+      obsmax[2] = ftable->GetObsBinsUpBoundsMax(2); // by default derived from binning in obs2
+      if ( SteeringPars["obs2max"] ) {
+         ftable->GetParameterFromSteering("obs2max",obsmax[2]);
+      }
+   }
+   jetclusfj.setup(static_cast<fastjet_jets::JetAlgorithm>(jetalgo), jetsize, overlapthreshold);
 }
 
 // --- fastNLO v2.2: class UserHHC: analyze parton event (called once for each event)
-void UserHHC::userfunc(const event_hhc& p, const amplitude_hhc& amp)
-{
-    say::debug["UserHHC::userfunc"] << "---------- UserHHC::userfunc called ----------" << endl;
-
-    // --- fastNLO user:
-    //     Here is your playground where you compute your observables and
-    //     scales for each jet or event.
-    //     The bin number ("obsbin") gets passed to fastNLO's table filling code.
-    //     Usually, pT and E are in GeV, but this may be changed.
-    //     ATTENTION: Scales must always be in GeV!
-
-    // apply the jet algorithm to partonic 4-vector array p of NLOJet++
-    pj = jetclusfj(p);
-    unsigned int nj = pj.upper();
-
-    // --- check on minimal and maximal no. of jets
-    // ATTENTION: In principle, without cuts, there should always be two.
-    //            For efficiency reasons though, our interface to the fastjet algorithms
-    //            requires a minimal jet pT of 1 GeV. If this is a problem, the ptmin value
-    //            in fj-jets.cc needs to be changed.
-    // There should never be more than four jets in NLOJet++
-    if (nj < 1) {
-        say::debug["ScenarioCode"] << "This event from NLOJet++ has no jets with pT > 1 GeV. Skipped!" << endl;
-        return;
-    } else if (nj > 4) {
-        say::error["ScenarioCode"] << "This event from NLOJet++ has more than four jets, which should never happen. Aborted!" << endl;
-        exit(1);
-    }
-
-    // --- give some debug output before selection and sorting
-    if (say::info.GetSpeak()) {
-        for (unsigned int i = 1; i <= nj; i++) {
-            double pti = pj[i].perp();
-            double yi = pj[i].rapidity();
-            double etai = pj[i].prapidity();
-            say::info["ScenarioCode"] << "before cuts: jet # i, pt, y, eta: " << i << ", " << pti << ", " << yi << ", " << etai << endl;
-        }
-    }
-
-    // --- select jets in y (lpseudo = false) or eta (lpseudo = true) and ptjmin
-    //     note: failing jets are not deleted but moved to the end of the jet array pj!
-    fNLOSelector SelJets(yetajmin, yetajmax, ptjmin, lpseudo);
-
-    // --- count number of selected jets left at this stage
-    size_t njet = std::remove_if(pj.begin(), pj.end(), SelJets) - pj.begin();
-    if ((int)njet < Njetmin)
-        return; // Nothing to be done
-
-    // --- sort selected n jets at beginning of jet array pj, by default decreasing in pt
-    static fNLOSorter SortJets;
-    std::sort(pj.begin(), pj.begin() + njet, SortJets);
-
-    // --- give some debug output after selection and sorting
-    if (say::info.GetSpeak()) {
-        say::info["ScenarioCode"] << "# jets before and after phase space cuts: nj, njet = " << nj << ", " << njet << endl;
-        if (!lpseudo) {
-            say::info["ScenarioCode"] << "phase space cuts: yjmin, yjmax, ptjmin: " << yetajmin << ", " << yetajmax << ", " << ptjmin << endl;
-        } else {
-            say::info["ScenarioCode"] << "phase space cuts: etajmin, etajmax, ptjmin: " << yetajmin << ", " << yetajmax << ", " << ptjmin << endl;
-        }
-        for (unsigned int i = 1; i <= njet; i++) {
-            double pti = pj[i].perp();
-            double yi = pj[i].rapidity();
-            double etai = pj[i].prapidity();
-            say::info["ScenarioCode"] << "after cuts: jet # i, pt, y, eta: " << i << ", " << pti << ", " << yi << ", " << etai << endl;
-        }
-    }
-
-    // ---- fastNLO v2.2
-    // Analyze inclusive jets in jet loop
-    // set one possible scale choice (Attention: Only correct if jets sorted descending in pT)
-    double ptmax = pj[1].perp();
-    for (unsigned int k = 1; k <= njet; k++) {
-
-        // derive some jet quantities
-        double pt = pj[k].perp();
-        double yeta;
-        if (!lpseudo) {
-            yeta = abs(pj[k].rapidity());
-        } else {
-            yeta = abs(pj[k].prapidity());
-        }
-        double phi = atan2(pj[k].Y(), pj[k].X());
-
-        // --- calculate observable of nth dimension
-        for (int i = 0; i < NDim; i++) {
-            switch (obsdef[i]) {
-            case PTJETGEV:
-                // jet pT
-                obs[i] = pt;
-                break;
-            case YJET:
-                // jet rapidity
-                obs[i] = yeta;
-                break;
-            case PHIJET:
-                // jet azimuthal angle
-                obs[i] = phi;
-                break;
-            default:
-                say::error["ScenarioCode"] << "Observable not yet implemented, aborted!" << endl;
-                say::error["ScenarioCode"] << "DimLabel[" << i << "] = " << DimLabel[i] << endl;
-                say::error["ScenarioCode"] << "Please complement this scenario to include the requested observable." << endl;
-                exit(1);
-            }
-        }
-
-        // --- Further Njet phase space cuts?
-        if (obsmin[0] <= obs[0] && obs[0] < obsmax[0] && (NDim < 2 || (obsmin[1] <= obs[1] && obs[1] < obsmax[1])) && (NDim < 3 || (obsmin[2] <= obs[2] && obs[2] < obsmax[2]))) {
-
-            // --- set the renormalization and factorization scales
-            // --- calculate the requested scales
-            for (unsigned int i = 0; i < 2; i++) {
-                switch (mudef[i]) {
-                case PTMAX:
-                    // maximal jet pT
-                    mu[i] = ptmax;
-                    break;
-                case PTJET:
-                    // jet pT
-                    mu[i] = pt;
-                    break;
-                default:
-                    say::error["ScenarioCode"] << "Scale not yet implemented, aborted!" << endl;
-                    say::error["ScenarioCode"] << "ScaleLabel[" << i << "] = " << ScaleLabel[i] << endl;
-                    say::error["ScenarioCode"] << "Please complement this scenario to include the requested scale." << endl;
-                    exit(1);
-                }
-            }
-            static vector<double> scalevars;
-            if (!lFlexibleScaleTable)
-                scalevars = ftable->GetScaleVariations();
-
-            // get matrix elements
-            static vector<fnloEvent> contribsflex;
-            static vector<vector<fnloEvent> > contribsfix;
-            if (lFlexibleScaleTable) {
-                contribsflex = UsefulNlojetTools::GetFlexibleScaleNlojetContribHHC(p, amp);
-            } else {
-                contribsfix = UsefulNlojetTools::GetFixedScaleNlojetContribHHC(p, amp, mu[0], scalevars);
-            }
-
-            // set and fill scenario specific quantites
-            fnloScenario scen;
-            if (NDim < 1 || NDim > 3) {
-                say::error["ScenarioCode"] << "Less than 1D(?!) or more than 3D binning not implemented here, aborted!" << endl;
-                say::error["ScenarioCode"] << "DifferentialDimension NDim = " << NDim << endl;
-            }
-            for (int i = 0; i < NDim; i++) {
-                scen.SetObservableDimI(obs[i], i);
+void UserHHC::userfunc(const event_hhc& p, const amplitude_hhc& amp) {
+   if ( say::debug.GetSpeak() ) {
+      say::debug["UserHHC::userfunc"] << "---------- UserHHC::userfunc called ----------" << endl;
+      say::debug["ScenarioCode"]  << "==================== Start of event ====================" << endl;
+   }
+
+   // --- fastNLO user:
+   //     Here is your playground where you compute your observables and
+   //     scales for each jet or event.
+   //     The bin number ("obsbin") gets passed to fastNLO's table filling code.
+   //     Usually, pT and E are in GeV, but this may be changed.
+   //     ATTENTION: Scales must always be in GeV!
+
+   // derive partonic HT scale
+   double htp = 0.;
+   unsigned int np = p.upper();
+   for (unsigned int i = 1; i <= np; i++) {
+      htp += p[i].perp();
+      // --- debug output
+      if ( say::debug.GetSpeak() ) {
+         say::debug["ScenarioCode"] << "partonic HT: parton # i, pt, htp: " << i << ", " << p[i].perp() << ", " << htp << endl;
+      }
+   }
+
+   // apply the jet algorithm to partonic 4-vector array p of NLOJet++
+   pj = jetclusfj(p);
+   unsigned int nj = pj.upper();
+
+   // --- check on minimal and maximal no. of jets
+   // ATTENTION: In principle, without cuts, there should always be two.
+   //            For efficiency reasons though, our interface to the fastjet algorithms
+   //            requires a minimal jet pT of 1 GeV. If this is a problem, the ptmin value
+   //            in fj-jets.cc needs to be changed.
+   // There should never be more than four jets in NLOJet++
+   if (nj < 1) {
+      say::debug["ScenarioCode"] << "This event from NLOJet++ has no jets with pT > 1 GeV. Skipped!" << endl;
+      return;
+   } else if (nj > 4) {
+      say::error["ScenarioCode"] << "This event from NLOJet++ has more than four jets, which should never happen. Aborted!" << endl;
+      exit(1);
+   }
+
+   // --- give some debug output before selection and sorting
+   if ( say::debug.GetSpeak() ) {
+      for (unsigned int i=1; i<=nj; i++) {
+         double pti  = pj[i].perp();
+         double yi   = pj[i].rapidity();
+         double etai = pj[i].prapidity();
+         say::debug["ScenarioCode"] << "before cuts: jet # i, pt, y, eta: " << i << ", " << pti << ", " << yi << ", " << etai << endl;
+      }
+   }
+
+   // --- select jets in y (lpseudo = false) or eta (lpseudo = true) and ptjmin
+   //     note: failing jets are not deleted but moved to the end of the jet array pj!
+   fNLOSelector SelJets (yetajmin,yetajmax,ptjmin,lpseudo);
+
+   // --- count number of selected jets left at this stage
+   size_t njet = std::remove_if(pj.begin(), pj.end(), SelJets) - pj.begin();
+   if ( (int)njet < Njetmin ) return; // Nothing to be done
+
+   // --- sort selected n jets at beginning of jet array pj, by default decreasing in pt
+   static fNLOSorter SortJets;
+   std::sort(pj.begin(), pj.begin() + njet, SortJets);
+
+   // --- give some debug output after selection and sorting
+   if ( say::debug.GetSpeak() ) {
+      say::debug["ScenarioCode"] << "# jets before and after phase space cuts: nj, njet = " << nj << ", " << njet << endl;
+      if ( ! lpseudo ) {
+         say::debug["ScenarioCode"] << "phase space cuts: yjmin, yjmax, ptjmin: " << yetajmin << ", " << yetajmax << ", " << ptjmin << endl;
+      } else {
+         say::debug["ScenarioCode"] << "phase space cuts: etajmin, etajmax, ptjmin: " << yetajmin << ", " << yetajmax << ", " << ptjmin << endl;
+      }
+      for (unsigned int i=1; i<=njet; i++) {
+         double pti  = pj[i].perp();
+         double yi   = pj[i].rapidity();
+         double etai = pj[i].prapidity();
+         say::debug["ScenarioCode"] << "after cuts: jet # i, pt, y, eta: " << i << ", " << pti << ", " << yi << ", " << etai << endl;
+      }
+   }
+
+   // ---- fastNLO v2.2
+   // Analyze inclusive jets in jet loop
+   // set one possible scale choice (Attention: Only correct if jets sorted descending in pT)
+   double ptmax = pj[1].perp();
+   for (unsigned int k = 1; k <= njet; k++) {
+
+      // --- calculate observable of nth dimension
+      for ( int i = 0; i<NDim; i++ ) {
+         switch(obsdef[i]) {
+         case PTJETGEV :
+            // jet pT
+            obs[i] = pj[k].perp();
+            break;
+         case YJET :
+            // jet rapidity
+            obs[i] = abs(pj[k].rapidity());
+            break;
+         case ETAJET :
+            // jet pseudorapidity
+            obs[i] = abs(pj[k].prapidity());
+            break;
+         case PHIJET :
+            // jet azimuthal angle
+            obs[i] = atan2(pj[k].Y(), pj[k].X());
+            break;
+         default :
+            say::error["ScenarioCode"] << "Observable not yet implemented, aborted!" << endl;
+            say::error["ScenarioCode"] << "DimLabel[" << i << "] = " << DimLabel[i] << endl;
+            say::error["ScenarioCode"] << "Please complement this scenario to include the requested observable." << endl;
+            exit(1);
+         }
+      }
+
+      // --- give some debug output before final selection
+      if ( say::debug.GetSpeak() ) {
+         say::debug["ScenarioCode"]  << "---------------- Before final selection ----------------" << endl;
+         for ( int i = 0; i<NDim; i++ ) {
+            say::debug["ScenarioCode"]  << "Obs. min/max values: " << i <<  " : obsmin = " << obsmin[i] << ", obs = " << obs[i] << ", obsmax = " << obsmax[i] << endl;
+         }
+      }
+
+      // --- Further Njet phase space cuts?
+      if ( obsmin[0] <= obs[0] && obs[0] < obsmax[0] &&
+           (NDim < 2 || (obsmin[1] <= obs[1] && obs[1] < obsmax[1])) &&
+           (NDim < 3 || (obsmin[2] <= obs[2] && obs[2] < obsmax[2])) ) {
+
+         // --- jet accepted
+         if ( say::debug.GetSpeak() ) {
+            say::debug["ScenarioCode"]  << "----------------- Event/jet accepted! ------------------" << endl;
+            say::debug["ScenarioCode"]  << "====================  End of event  ====================" << endl;
+         }
+
+         // --- set the renormalization and factorization scales
+         // --- calculate the requested scales
+         for ( unsigned int i = 0; i < 2; i++ ) {
+            switch(mudef[i]) {
+            case PTMAX :
+               // maximal jet pT
+               mu[i] = ptmax;
+               break;
+            case PTJET :
+               // jet pT
+               mu[i] = pj[k].perp();
+               break;
+         case HTP :
+            // partonic sum pT
+            mu[i] = htp;
+            say::debug["ScenarioCode"]  << "HTP scale values: " << i <<  " : mu[i] = " << mu[i] << endl;
+            break;
+         case HTPHALF :
+            // half partonic sum pT
+            mu[i] = htp/2.;
+            say::debug["ScenarioCode"]  << "HTPHALF scale values: " << i <<  " : mu[i] = " << mu[i] << endl;
+            break;
+            default :
+               say::error["ScenarioCode"] << "Scale not yet implemented, aborted!" << endl;
+               say::error["ScenarioCode"] << "ScaleLabel[" << i << "] = " << ScaleLabel[i] << endl;
+               say::error["ScenarioCode"] << "Please complement this scenario to include the requested scale." << endl;
+               exit(1);
             }
+         }
+         static vector<double> scalevars;
+         if ( ! lFlexibleScaleTable ) scalevars = ftable->GetScaleVariations();
+
+         // get matrix elements
+         static vector<fnloEvent> contribsflex;
+         static vector< vector<fnloEvent> > contribsfix;
+         if (lFlexibleScaleTable) {
+            contribsflex = UsefulNlojetTools::GetFlexibleScaleNlojetContribHHC(p,amp);
+         } else {
+            contribsfix  = UsefulNlojetTools::GetFixedScaleNlojetContribHHC(p,amp,mu[0],scalevars);
+         }
+
+         // set and fill scenario specific quantites
+         fnloScenario scen;
+         if ( NDim < 1 || NDim > 3 ) {
+            say::error["ScenarioCode"] << "Less than 1D(?!) or more than 3D binning not implemented here, aborted!" << endl;
+            say::error["ScenarioCode"] << "DifferentialDimension NDim = " << NDim << endl;
+         }
+         for ( int i = 0; i<NDim; i++ ) {
+            scen.SetObservableDimI( obs[i], i );
+         }
+
+         scen.SetObsScale1( mu[0] );   // must be consistent with 'mu' from contribs
+         if (lFlexibleScaleTable) {
+            scen.SetObsScale2( mu[1] );
+            ftable->FillAllSubprocesses(contribsflex,scen);
+         } else {
+            ftable->FillAllSubprocesses(contribsfix,scen);
+         }
+      } else {
+         // --- jet rejected
+         if ( say::debug.GetSpeak() ) {
+            say::debug["ScenarioCode"]  << "----------------- Event/jet rejected! ------------------" << endl;
+            say::debug["ScenarioCode"]  << "====================  End of event  ====================" << endl;
+         }
+      }
+   }
+}
 
-            scen.SetObsScale1(mu[0]); // must be consistent with 'mu' from contribs
-            if (lFlexibleScaleTable) {
-                scen.SetObsScale2(mu[1]);
-                ftable->FillAllSubprocesses(contribsflex, scen);
-            } else {
-                ftable->FillAllSubprocesses(contribsfix, scen);
-            }
-        }
-    }
+// --- define function for azimuthal angular distance in [-Pi,Pi]
+double dphi(double phi2, double phi1) {
+   double delta_phi = phi2-phi1;
+   if (delta_phi > M_PI) {
+      delta_phi -= 2*M_PI;
+   } else if (delta_phi < -M_PI) {
+      delta_phi += 2*M_PI;
+   }
+   return delta_phi;
 }