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;
}