Page MenuHomeHEPForge

No OneTemporary

diff --git a/input/default.in b/input/default.in
index a699ca7..3a49149 100644
--- a/input/default.in
+++ b/input/default.in
@@ -1,323 +1,339 @@
# Process settings
sroot = 7e3 # Center-of-mass energy
ih1 = 1 # Hadron 1: 1 for proton, -1 for antiproton
ih2 = 1 # Hadron 2: 1 for proton, -1 for antiproton
nproc = 3 # Process: 1) W+; 2) W-; 3) Z/gamma*
# Perturbative order
# fixedorder_only = true # Evaluate predictions at fixed order
# fixedorder_only = false # Evaluate predictions including qt-resummation
fixedorder_only = false
order = 2 # QCD order: 0) LO(+LL), 1) NLO(+NLL), 2) NNLO(+NNLL)
# Non-perturbative form factor, S_NP = exp(-npff(b))
# 0: Gaussian (BLNY) npff(b) = (g1 + g2*log(m/Q0) + g3*log(100*m/sqrt(s)))*b^2
# 1: Exponential npff(b) = e*b
# 2: Collins-Rogers npff(b) = g0*(1-exp(-(Cf*alphas(b0/bstar))*b^2)/(PI*g0*blim^2))*ln(m/Q0^2)
npff = 0
# Gaussian options (npff = 0)
g1 = 0.8
g2 = 0
g3 = 0
Q0 = 1
# Exponential option (npff = 1)
e = 0.0
# Collins-Rogers options (npff = 2)
g0 = 0.3
Q0 = 1.6
flavour_kt = false # Flavour-dependent g1 parameters
g1_uv = 0.5 # u-valence
g1_us = 0.5 # u-sea
g1_dv = 0.5 # d-valence
g1_ds = 0.5 # d-sea
g1_ss = 0.5 # strange
g1_ch = 0.5 # charm
g1_bo = 0.5 # bottom
g1_gl = 0.5 # gluon
# PDF settings
LHAPDFset = CT10nnlo # PDF set from LHAPDF
LHAPDFmember = 0 # PDF member
PDFerrors = false # Calculate PDF errors
# Functional form of QCD scales (mV: wmass or zmass, pT: boson transverse momentum, mjj: dijet invariant mass)
#0: mu^2 = mV^2
#1: mu^2 = mll^2
#2: mu^2 = mll^2+pT^2
#3: mu^2 = mll^2+pT^2+mjj^2
fmuren = 1 # Functional form of the renormalisation scale
fmufac = 1 # Functional form of the factorisation scale
fmures = 1 # Functional form of the resummation scale (forms 2 and 3 are equivalent to 1)
# QCD scale settings
kmuren = 0.5 # Scale factor for the renormalisation scale
kmufac = 0.5 # Scale factor for the factorisation scale
kmures = 0.5 # Scale factor for the resummation scale
#PDF matching scales
kmuc = 1. # Scale factor for the charm matching scale
kmub = 1. # Scale factor for the bottom matching scale
kmut = 1. # Scale factor for the top matching scale
# EW scheme
#0: Input: alpha(mZ), zmass, xw; Derived: wmass, Gf
#1: Input: Gf, wmass, zmass; Derived: xw, alpha(mZ) [Gmu scheme]
#2: Input: Gf, alpha(mZ), xw, Gf; Derived: wmass, zmass
#3: All masses and couplings determined by inputs
ewscheme = 1
# EW parameters
Gf = 1.1663787e-5 # G-Fermi
zmass = 91.1876 # Mass of the Z boson
wmass = 80.385 # Mass of the W boson
xw = 0.23153 # Weak-mixing angle (not used in the Gmu scheme)
aemmz = 7.7585538055706e-03 # alpha_EM(MZ) (not used in the Gmu scheme)
# W and Z total widths used in the propagator are determined by the following inputs
zwidth = 2.4950 # Width of the Z boson
wwidth = 2.091 # Width of the W boson
runningwidth = false # Use Z and W propagators including energy-dependent width effects
# CKM matrix
Vud = 0.97427
Vus = 0.2253
Vub = 0.00351
Vcd = 0.2252
Vcs = 0.97344
Vcb = 0.0412
# Z/gamma* coupling to quarks
Zuu = 1.0
Zdd = 1.0
Zss = 1.0
Zcc = 1.0
Zbb = 1.0
# Include virtual photon and interference in Z/gamma* production
useGamma = true
# Resummation parameters
qtcutoff = 0.02 # Resummation cutoff in GeV
# modlog = true # Modified logarithms in the Sudakov L~ = log( (Q*b/b0)^2 + 1)
# modlog = false # Canonical logarithms in the Sudakov L = log( (Q*b/b0)^2 )
modlog = true
+# Switch off exponentiation of C coefficients
+nocexp = false
+
+# Sum all logs before exponentiation (helpful with bprescription = 0, but create issues with bprescription = 2)
+sumlogs = false
+
# Prescription to avoid the Landau pole in the Bessel inverse transform
# 0: bstar prescription, which freezes b at bmax: b -> bstar = b/sqrt(1+b^2/bmax^2)
# 1: Integrate all the way up to the Landau singularity b_L = b0/Q * exp(1/(2*beta0*alphas))
# 2: Minimal prescription (complex plane)
# 3: Minimal prescription (real axis)
bprescription = 0
blim = -1 #blim for the bstar prescription, applies to bprescription = 0. A positive value set a fixed bmax=blim, a negative values sets bmax=b_L/(-blim), where b_L is the Landau singularity.
+blim_pdf = 0 #if not zero, set a different blim for the PDF evolution
+blim_sudakov = 0 #if not zero, set a different blim for the Sudakov form factor
+blim_cexp = 0.5 #if not zero, set a different blim for the C exponentiation
+
+bstar_pdf = false #force bstar prescription for the PDF evolution
+bstar_sudakov = false #force bstar prescription for the Sudakov form factor
+bstar_cexp = false #force bstar prescription for the C exponentiation
phibr = 4. #set arg(z) as phib = pi/phibr for the integration contour in the complex plane for bprescription = 2 (should be set phibr > 4. )
bcf = 0.5 #select the point bc = bcf * b_L, where the integration contour is bended in the complex plane, as a fraction of the Landau singularity b_L. Applies to bprescription = 2 or 3
# Strategy for the direct Mellin transform of PDFs at the factorization scale
# Set to 1 to use the dyres approximation of PDFs and integration contour in the complex plane for the Mellin inversion
# Set to 0 to use exact PDFs and straight line contour in the complex plane for the Mellin inversion
opts_approxpdf = 0
# x-to-N direct Mellin transform of PDFs
opts_pdfintervals = 100 # Number of intervals for integration of PDF moments
pdfrule = 200 # Gaussian rule for the x-to-N Mellin transform
# Type of PDF evolution
#0: DYRES (backward evolution)
#1: Pegasus as DYRES (backward evolution)
#2: Mellin transform (forward evolution)
#3: Pegasus VFN (forward evolution)
#4: Pegasus VFN with alphas from DYRES (forward evolution)
evolmode = 0
+# Evolve PDFs from mufac (instead of mures) to b0/b
+mufevol = false
+
# Settings for the inverse Hankel transform
bintaccuracy = 1.0e-4 # Accuracy of the integration
# Settings for the inverse Mellin integrations
mellininv = 0 # Strategy for the Mellin inversion (0 Gauss-Legendre, 1 Talbot)
mellinintervals = 1 # Number of intervals
mellinrule = 64 # Number of nodes
# Options for the Mellin inversion with Gauss-Legendre integration
zmax = 27. # Upper limit of the contour in the imaginary axis
cpoint = 1. # Point of intersection of the contour with the real axis
phi = 0.5 # Angle between the real axis and the linear contour in units of pi
mellincores = 1 # Number of parallel threads for the Mellin integration
mellin1d = false # Use 1d (y-integrated) or 2d (y-dependent) Mellin inversion
xspace = false # Access PDFs in Bjorken-x space, without Mellin inversion (option available only in the LL case where the convolution is trivial)
# Resummation damping
damp = true
# Resummation damping function
# 1: Gaussian: exp(-(k*mll-qt)^2)/(delta*mll)^2
# 2: Exponential: exp((k*mll)^-qt^2)/(delta*mll)^2
# 3: Cosine: cos(PI/(delta*mll)*(qt-k*mll))+1)/2
dampmode = 1
dampk = 0.5
dampdelta = 0.5
# qt-subtraction cut-off. Both conditions are applied, at least one between qtcut and xqtcut must be > 0
xqtcut = 0.008 # cutoff on qt/m
qtcut = 0. # cutoff on qt
# Integration settings
rseed = 123456 # Random seed for MC integration
# Term switches
doBORN = true
doCT = true
doVJ = true
doVJREAL = true
doVJVIRT = true
# Integration type: true -> quadrature, false -> vegas
BORNquad = true
CTquad = true
VJquad = true
# Integration type (advanced settings, override BORNquad, CTquad, VJquad options if set > -1)
intDimRes = -1 # Resummed term (1, 2 or 3 for quadrature, or 4 for vegas)
intDimBorn = -1 # Born term (2 for quadrature, 4 or 6 for vegas)
intDimCT = -1 # Counter term (1, 2 or 3 for quadrature, or 6, or 8 for vegas)
intDimVJ = -1 # V+jet term (3 for quadrature, 5 for quadrature with cuts, or 7 for vegas)
# Multithreading parallelisation
cores = 0 # Number of parallel threads (0 for turning off parallelisation)
# Cuba settings
cubaverbosity = 0 # Cuba info messsages, from 0 to 3
cubanbatch = 1000 # The batch size for sampling in Cuba vegas integration
niterBORN = 5 # Only for 2d and 3d cuhre integration of resummed part
niterCT = 5 # Only for 2d and 3d cuhre integration of counter term
niterVJ = 10 # Only for 3d cuhre integration of V+J
#Vegas settings
vegasncallsBORN = 1000 # only for res 4d vegas integration
#vegasncallsRES = 1000 # only for res 4d vegas integration
vegasncallsCT = 100000 # only for 6d and 8d vegas integration of the counter term
vegasncallsVJLO = 10000000 # only for lo 7d vegas integration
vegasncallsVJREAL = 10000000 # only for real 10d vegas integration
vegasncallsVJVIRT = 1000000 # only for virt 8d vegas integration
vegascollect = false # collect points from all the vegas iterations (true) or only from the last iteration (false)
# cubature settings
pcubature = true # Use Cuhre (false ) or pcubature (true) integration in quadrature mode
relaccuracy = 1e-3 # target relative uncertainty of each term
absaccuracy = 0 # target absolute uncertainty of each term in fb
# Advanced integration settings
# Number of intervals and gaussian rule for the rapidity integrations in the 2dim resummed piece
yintervals = 1
yrule = 64
# Number of intervals and gaussian rule for the qt integration in the 2dim counter term
qtintervals = 1
qtrule = 64
# Number of intervals and gaussian rule for the alfa beta scaled-PDF integration in the counter term and born fixed order term
abintervals = 1
abrule = 64
# Gaussian rule for the phi integration in the V+J 5d LO term when makecuts is false
vjphirule = 20
# Settings for the z1, z2 integration in the V+J 3d NLO term
zrule = 64
# Settings for the x integration in the V+J 3d delta term
xrule = 200
# costh CS boundaries
costhmin = -1
costhmax = +1
# Lepton cuts
# Total cross section or with lepton cuts
makecuts = false
# charged leptons cuts
lptcut = 20
lycut = 2.5 # absolute rapidity cut
lycutmin
# leptons and antileptons cuts
lepptcut = 0
lepycut = 1000
alpptcut = 0
alpycut = 1000
#absolute-rapidity-ordered leptons (central and forward)
lcptcut = 0
lcymin = 0
lcymax = 1000
lfptcut = 0
lfymin = 0
lfymax = 1000
# cuts on neutrino and transverse mass
etmisscut = 0
mtcut = 0
#costh CS
cthCSmin = -1
cthCSmax = +1
# integration types and settings for costh phi_lep phase space
cubaint = false # integration with Cuba Suave
trapezint = false # trapezoidal rule for the phi_lep integration and semi-analytical for costh
quadint = true # quadrature rule for the phi_lep integration and semi-analytical for costh
suavepoints = 1000000 # number of points for suave integration, newpoints is set to suavepoints/10;
nphitrape = 1000 # number of steps for trapezoidal rule of phi_lep integration
phirule = 4 # quadrature rule of phi_lep integration
phiintervals = 20 # number of segments for quadrature rule of phi_lep integration
ncstart = 100 # starting sampling for the costh semi-analytical integration (common settings for the trapezoidal and quadrature rules)
# qt-recoil prescriptions
qtrec_naive = false
qtrec_cs = true
qtrec_kt0 = false
# Debug settings
timeprofile = false # debug and time profile resummation integration
verbose = false # debug and time profile costh phi_lep integration
# Output settings
output_filename = results # output filename
texttable = true # dump result table to text file (including pdf variations)
redirect = false # redirect stdout and stderr to log file (except for gridverbose output)
unicode = false # use unicode characters for the table formatting
silent = false # no output on screen (except for gridverbose output)
makehistos = true # fill histograms
gridverbose = false # printout number of events to keep job alive when running on grid
# Compute total (-1) or helicity cross sections (0-7)
helicity = -1
# binning
# normalise cross sections by bin width
ptbinwidth = false
ybinwidth = false
mbinwidth = false
# Force to loop over all bins even you have all Vegas integrands
force_binsampling = false
# qt, y, m bins
qt_bins = [ 0 2 4 6 8 10 12 14 16 18 22 26 30 34 38 42 46 50 54 60 70 80 100 150 200 300 800 ]
y_bins = [ -5 5 ]
m_bins = [ 66 116 ]
# binning for user testing histogram
biganswer_bins = [ 41 42 43 ]
diff --git a/src/settings.C b/src/settings.C
index 7ac35a3..f894a59 100644
--- a/src/settings.C
+++ b/src/settings.C
@@ -1,1106 +1,1137 @@
#include <cmath>
#include <cstring>
#include <iostream>
#include <stdexcept>
#include "settings.h"
#include "dyres_interface.h"
#include "interface.h"
#include "phasespace.h"
// CXX option parser: https://raw.githubusercontent.com/jarro2783/cxxopts/master/src/cxxopts.hpp
#include "cxxopts.hpp"
#include "handy_typdefs.h"
settings opts;
binning bins;
void settings::parse_options(int argc, char* argv[]){
// Declare the supported options.
po::Options args(argv[0], " [config.in] \n\n"
" Fast Drell-Yan Monte Carlo and quadrature integrator. \n\n"
" NOTE: Command line options are overiding the default and config file settings."
);
// Hidden arguments
args.add_options("Hidden")
("conf_file" , "Name of output file. ", po::value<string>()->default_value("") )
;
// Program options
args.add_options("")
("h,help" , "Print this help and die." )
("v,verbose" , "Be verbose" )
("q,small-stat" , "Set quick run with small stat." )
("p,proc" , "Set process [z0/wp/wm]" , po::value<string>() )
("c,collider" , "Set beam conditions [tev2/lhc7/lhc8/lhc13/lhc14]" , po::value<string>() )
("o,order" , "Set order [0:LO, 1:NLL+NLO, 2:NNLL+NNLO]" , po::value<int>() )
("f,fixedorder" , "Set fixed order only" )
("e,resummation" , "Set resummation" )
("t,term" , "Set term [BORN,CT,VJ,VJREAL,VJVIRT,ALL]" , po::value<string>() )
("r,seed" , "Set random seed [integer]" , po::value<int>() )
("s,pdfset" , "Set PDF set [LHAPDF name]" , po::value<string>() )
("m,pdfvar" , "Set PDF member [integer/all]" , po::value<string>() )
("g,gpar" , "Set non-perturbative Sudakov parameter" , po::value<double>() )
- ("R,kmuren" , "Set realative renormalization scale" , po::value<double>() )
- ("F,kmufac" , "Set realative factorization scale" , po::value<double>() )
+ ("R,kmuren" , "Set relative renormalization scale" , po::value<double>() )
+ ("F,kmufac" , "Set relative factorization scale" , po::value<double>() )
+ ("Q,kmures" , "Set relative resummation scale" , po::value<double>() )
("qtbins" , "Set equdistan binning for mass [N,lo,hi]" , po::value<string>() )
("ybins" , "Set equdistan binning for qt [N,lo,hi]" , po::value<string>() )
("mbins" , "Set equdistan binning for y [N,lo,hi]" , po::value<string>() )
("grid" , "Be verbose to avoid jobs being killed")
;
// Parse
try {
args.parse_positional( std::vector<string>({"conf_file"}) );
args.parse(argc,argv);
}
catch (po::OptionException &e){
printf("%s\n", args.help().c_str());
printf("Bad arguments: %s \n",e.what());
throw e;
}
// Print help and die
if (args.count("help")) {
throw QuitProgram(args.help().c_str());
}
// load config file (or default settings)
readfromfile ( args["conf_file"].as<string>() );
bins.readfromfile ( args["conf_file"].as<string>() );
externalpdf = false;
// NOTE: Command line options are overiding the default and config file settings.
// verbose
if (args.count("verbose"))
{
verbose=true;
cubaverbosity=3;
}
if (args.count("grid"))
{
gridverbose=true;
redirect=true;
cubaverbosity=3;
}
if (opts.redirect)
{
string logfile = output_filename + ".log";
freopen(logfile.c_str(),"w",stdout);
}
// rseed
if (args.count("seed")) rseed=args["seed"].as<int>();
// order
if (args.count("order")) order=args["order"].as<int>();
// fixed order
if (args.count("fixedorder"))
fixedorder = true;
if (args.count("resummation"))
fixedorder = false;
// small stat
if (args.count("small-stat")){
niterBORN = 1;
niterCT = 1;
niterVJ = 1;
vegasncallsBORN = (fixedorder ? 1e5 : 1e3);
vegasncallsCT = 1e4;
vegasncallsVJLO = 1e5;
vegasncallsVJREAL = 1e5;
vegasncallsVJVIRT = 1e5;
relaccuracy = 0.1;
}
// proc
if (args.count("proc")) {
string val=args["proc"].as<string>();
ToLower(val);
if (val == "z0"){ nproc=3;
} else if (val == "wp"){ nproc=1;
} else if (val == "wm"){ nproc=2;
} else {
throw QuitProgram("Unsupported value of proc: "+val);
}
}
// PDF
if (args.count("pdfset")) LHAPDFset=args["pdfset"].as<string>();
if (args.count("pdfvar")) {
// set default values
PDFerrors=false;
LHAPDFmember=0;
string val=args["pdfvar"].as<string>();
ToLower(val);
if (val=="all") { PDFerrors=true;
} else if (IsNumber(val)) { LHAPDFmember=stod(val);
} else {
throw QuitProgram("Unsupported value of pdfvar: "+val);
}
}
//Gpar
if (args.count("gpar")) {
opts.g1 = args["gpar"].as<double>();
}
if (args.count("kmuren")) {
opts.kmuren = args["kmuren"].as<double>();
}
if (args.count("kmufac")) {
opts.kmufac = args["kmufac"].as<double>();
}
+ if (args.count("kmures")) {
+ opts.kmures = args["kmures"].as<double>();
+ }
// Collider
if (args.count("collider")) {
string val=args["collider"].as<string>();
ToLower(val);
if (val == "tev1" ){ sroot=1.80e3; ih1=1; ih2=-1;
} else if (val == "tev2" ){ sroot=1.96e3; ih1=1; ih2=-1;
} else if (val == "lhc7" ){ sroot=7.00e3; ih1=1; ih2=1;
} else if (val == "lhc5" ){ sroot=5.00e3; ih1=1; ih2=1;
} else if (val == "lhc8" ){ sroot=8.00e3; ih1=1; ih2=1;
} else if (val == "lhc13" ){ sroot=13.0e3; ih1=1; ih2=1;
} else if (val == "lhc14" ){ sroot=14.0e3; ih1=1; ih2=1;
} else {
throw QuitProgram("Unsupported value of collider: "+val);
}
}
// term
if (args.count("term")) {
// first turn off all terms
doBORN = doCT = doVJ = doVJREAL = doVJVIRT = false ;
string val=args["term"].as<string>();
ToUpper(val);
for (auto piece : Tokenize(val)) {
if ( piece == "BORN" ) { doBORN = true;
} else if ( piece == "CT" ) { doCT = true;
} else if ( piece == "VJ" ) { doVJ = true; doVJVIRT = true; doVJREAL = true;
} else if ( piece == "VJVIRT" ) { doVJ = true; doVJVIRT = true;
} else if ( piece == "VJREAL" ) { doVJ = true; doVJREAL = true;
} else if ( piece == "ALL" ) { doBORN = doCT = doVJ = doVJREAL = doVJVIRT = true;
} else {
throw QuitProgram("Unsupported value of term : "+piece);
}
}
}
//binning
parse_binning("qtbins" , bins.qtbins ,args);
parse_binning("ybins" , bins.ybins ,args);
parse_binning("mbins" , bins.mbins ,args);
// check consistency of settings
check_consistency();
}
// PROTOTYPE OF VARIABLE REGISTRATION
struct var_def{
string name = "";
void * ptr = NULL;
char type = '0'; // 0-group name, b-bool, s-string, i-int, d-doubl, v-VecDbl
string dump(){
SStream tmp;
if (type == '0') {
tmp << "# " << name;
} else {
tmp << name << " = ";
if (type != 'b'){
bool *val = (bool*) ptr;
tmp << ((*val) ? "true" : "false");
}
if (type != 'd'){
double *val = (double*) ptr;
tmp << *val;
}
if (type != 'i'){
int *val = (int*) ptr;
tmp << *val;
}
if (type != 's'){
string *val = (string*) ptr;
tmp << *val;
}
if (type != 'v'){
VecDbl *val = (VecDbl*) ptr;
tmp << "[ ";
for (double v: *val) tmp << v << " ";
tmp << " ]";
}
}
tmp << "\n";
return tmp.str();
}
};
Vec<var_def> var_reg;
void GroupName(string name){ var_def add; add.name=name; var_reg.push_back(add); };
void RegVar(string name, bool *var){ var_def add; add.name=name; add.ptr=var; add.type='b'; var_reg.push_back(add); };
void RegVar(string name, double *var){ var_def add; add.name=name; add.ptr=var; add.type='d'; var_reg.push_back(add); };
void RegVar(string name, int *var){ var_def add; add.name=name; add.ptr=var; add.type='i'; var_reg.push_back(add); };
void RegVar(string name, string *var){ var_def add; add.name=name; add.ptr=var; add.type='s'; var_reg.push_back(add); };
void RegVar(string name, VecDbl *var){ var_def add; add.name=name; add.ptr=var; add.type='v'; var_reg.push_back(add); };
void GetVar(string name, bool &var, InputParser &in){ var = in.GetBool(name) ; };
void GetVar(string name, double &var, InputParser &in){ var = in.GetNumber(name); };
void GetVar(string name, int &var, InputParser &in){ var = in.GetNumber(name); };
void GetVar(string name, string &var, InputParser &in){ var = in.GetString(name); };
void GetVar(string name, VecDbl &var, InputParser &in){ var.clear(); in.GetVectorDouble(name,var); };
template <class Tvar>
void addVariable(string name, Tvar& var, InputParser &in){
GetVar(name,var,in);
AddToDump(name,&var);
};
// #define GETVARIABLE(var) addVariable( ##var, var, in);
// END OF VARIABLE REGISTRATION
void settings::readfromfile(const string fname){
//read input settings from file
InputParser in(fname);
sroot = in.GetNumber ( "sroot" );
ih1 = in.GetNumber ( "ih1" );
ih2 = in.GetNumber ( "ih2" );
nproc = in.GetNumber ( "nproc" );
kmuren = in.GetNumber ( "kmuren" );
kmufac = in.GetNumber ( "kmufac" );
kmures = in.GetNumber ( "kmures" );
fmuren = in.GetNumber ( "fmuren" );
fmufac = in.GetNumber ( "fmufac" );
fmures = in.GetNumber ( "fmures" );
kmuc = in.GetNumber ( "kmuc" );
kmub = in.GetNumber ( "kmub" );
kmut = in.GetNumber ( "kmut" );
npff = in.GetNumber ( "npff" );
g1 = in.GetNumber ( "g1" );
g2 = in.GetNumber ( "g2" );
g3 = in.GetNumber ( "g3" );
e = in.GetNumber ( "e" );
g0 = in.GetNumber ( "g0" );
Q0 = in.GetNumber ( "Q0" );
flavour_kt = in.GetBool ( "flavour_kt" );
g1_uv = in.GetNumber ( "g1_uv" );
g1_us = in.GetNumber ( "g1_us" );
g1_dv = in.GetNumber ( "g1_dv" );
g1_ds = in.GetNumber ( "g1_ds" );
g1_ss = in.GetNumber ( "g1_ss" );
g1_ch = in.GetNumber ( "g1_ch" );
g1_bo = in.GetNumber ( "g1_bo" );
g1_gl = in.GetNumber ( "g1_gl" );
order = in.GetNumber ( "order" );
runningwidth = in.GetBool ( "runningwidth" );
rseed = in.GetNumber ( "rseed" );
blim = in.GetNumber ( "blim" );
+ blim_pdf = in.GetNumber ( "blim_pdf" );
+ blim_sudakov = in.GetNumber ( "blim_sudakov" );
+ blim_cexp = in.GetNumber ( "blim_cexp" );
+ bstar_pdf = in.GetBool ( "bstar_pdf" );
+ bstar_sudakov = in.GetBool ( "bstar_sudakov" );
+ bstar_cexp = in.GetBool ( "bstar_cexp" );
LHAPDFset = in.GetString ( "LHAPDFset" );
LHAPDFmember = in.GetNumber ( "LHAPDFmember" );
ewscheme = in.GetNumber ( "ewscheme" );
Gf = in.GetNumber ( "Gf" );
zmass = in.GetNumber ( "zmass" );
wmass = in.GetNumber ( "wmass" );
xw = in.GetNumber ( "xw" );
aemmz = in.GetNumber ( "aemmz" );
zwidth = in.GetNumber ( "zwidth" );
wwidth = in.GetNumber ( "wwidth" );
Vud = in.GetNumber ( "Vud" );
Vus = in.GetNumber ( "Vus" );
Vub = in.GetNumber ( "Vub" );
Vcd = in.GetNumber ( "Vcd" );
Vcs = in.GetNumber ( "Vcs" );
Vcb = in.GetNumber ( "Vcb" );
Zuu = in.GetNumber ( "Zuu" );
Zdd = in.GetNumber ( "Zdd" );
Zcc = in.GetNumber ( "Zcc" );
Zss = in.GetNumber ( "Zss" );
Zbb = in.GetNumber ( "Zbb" );
damp = in.GetBool ( "damp" );
dampk = in.GetNumber ( "dampk" );
dampdelta = in.GetNumber ( "dampdelta" );
dampmode = in.GetNumber ( "dampmode" );
qtcutoff = in.GetNumber ( "qtcutoff" );
modlog = in.GetBool ( "modlog" );
xqtcut = in.GetNumber ( "xqtcut" );
qtcut = in.GetNumber ( "qtcut" );
intDimRes = in.GetNumber ( "intDimRes" );
intDimBorn = in.GetNumber ( "intDimBorn" );
intDimCT = in.GetNumber ( "intDimCT" );
intDimVJ = in.GetNumber ( "intDimVJ" );
BORNquad = in.GetBool ( "BORNquad" );
CTquad = in.GetBool ( "CTquad" );
VJquad = in.GetBool ( "VJquad" );
fixedorder = in.GetBool ( "fixedorder_only" );
doBORN = in.GetBool ( "doBORN" );
doCT = in.GetBool ( "doCT" );
doVJ = in.GetBool ( "doVJ" );
doVJREAL = in.GetBool ( "doVJREAL" );
doVJVIRT = in.GetBool ( "doVJVIRT" );
cubaverbosity = in.GetNumber ( "cubaverbosity" );
cubacores = in.GetNumber ( "cores" );
cubanbatch = in.GetNumber ( "cubanbatch" );
niterBORN = in.GetNumber ( "niterBORN" );
niterCT = in.GetNumber ( "niterCT" );
niterVJ = in.GetNumber ( "niterVJ" );
vegasncallsBORN = in.GetNumber ( "vegasncallsBORN" );
vegasncallsCT = in.GetNumber ( "vegasncallsCT" );
vegasncallsVJLO = in.GetNumber ( "vegasncallsVJLO" );
vegasncallsVJREAL = in.GetNumber ( "vegasncallsVJREAL" );
vegasncallsVJVIRT = in.GetNumber ( "vegasncallsVJVIRT" );
vegascollect = in.GetBool ( "vegascollect" );
pcubature = in.GetBool ( "pcubature" );
relaccuracy = in.GetNumber ( "relaccuracy" );
absaccuracy = in.GetNumber ( "absaccuracy" );
costhmin = in.GetNumber ( "costhmin" );
costhmax = in.GetNumber ( "costhmax" );
makecuts = in.GetBool ( "makecuts" );
lptcut = in.GetNumber ( "lptcut" );
lycut = in.GetNumber ( "lycut" );
mtcut = in.GetNumber ( "mtcut" );
etmisscut = in.GetNumber ( "etmisscut" );
lepptcut = in.GetNumber ( "lepptcut" );
lepycut = in.GetNumber ( "lepycut" );
alpptcut = in.GetNumber ( "alpptcut" );
alpycut = in.GetNumber ( "alpycut" );
lcptcut = in.GetNumber ( "lcptcut" );
lcymin = in.GetNumber ( "lcymin" );
lcymax = in.GetNumber ( "lcymax" );
lfptcut = in.GetNumber ( "lfptcut" );
lfymin = in.GetNumber ( "lfymin" );
lfymax = in.GetNumber ( "lfymax" );
cthCSmin = in.GetNumber ( "cthCSmin" );
cthCSmax = in.GetNumber ( "cthCSmax" );
cubaint = in.GetBool ( "cubaint" ); //true # integration with Cuba Suave
trapezint = in.GetBool ( "trapezint" ); //false # trapezoidal rule for the phi_lep integration and semi-analytical for costh
quadint = in.GetBool ( "quadint" ); //false # quadrature rule for the phi_lep integration and semi-analytical for costh
suavepoints = in.GetNumber ( "suavepoints" ); //1000000 # number of points for suave integration, newpoints is set to suavepoints/10;
nphitrape = in.GetNumber ( "nphitrape" ); //1000 # number of steps for trapezoidal rule of phi_lep integration
phiintervals = in.GetNumber ( "phiintervals" );
phirule = in.GetNumber ( "phirule" );
ncstart = in.GetNumber ( "ncstart" ); //1000 # starting sampling for the costh semi-analytical integration (common settings for the trapezoidal and quadrature rules)
qtrec_naive = in.GetBool ( "qtrec_naive" ); //false
qtrec_cs = in.GetBool ( "qtrec_cs" ); //false
qtrec_kt0 = in.GetBool ( "qtrec_kt0" ); //true
timeprofile = in.GetBool ( "timeprofile" ); //false # debug and time profile resummation integration
verbose = in.GetBool ( "verbose" ); //false # debug and time profile costh phi_lep integration
gridverbose = in.GetBool ( "gridverbose" ); //false # debug and time profile costh phi_lep integration
texttable = in.GetBool ( "texttable" ); //
redirect = in.GetBool ( "redirect" ); //
unicode = in.GetBool ( "unicode" ); //
silent = in.GetBool ( "silent" ); //
makehistos = in.GetBool ( "makehistos" ); //
useGamma = in.GetBool ( "useGamma" );//
PDFerrors = in.GetBool ( "PDFerrors" );//
opts_.approxpdf_ = in.GetNumber ( "opts_approxpdf" ); //0
opts_.pdfintervals_ = in.GetNumber ( "opts_pdfintervals" ); //100
pdfrule = in.GetNumber ( "pdfrule" );
evolmode = in.GetNumber ("evolmode");
+ mufevol = in.GetBool ("mufevol");
+ nocexp = in.GetBool ("nocexp");
+ sumlogs = in.GetBool ("sumlogs");
bprescription = in.GetNumber ("bprescription");
bintaccuracy = in.GetNumber ( "bintaccuracy" );
phibr = in.GetNumber ( "phibr" );
bcf = in.GetNumber ( "bcf" );
mellininv = in.GetNumber ( "mellininv" );
mellinintervals = in.GetNumber ( "mellinintervals" );
mellinrule = in.GetNumber ( "mellinrule" );
zmax = in.GetNumber ( "zmax" );
cpoint = in.GetNumber ( "cpoint" );
phi = in.GetNumber ( "phi" );
mellincores = in.GetNumber ( "mellincores" );
mellin1d = in.GetBool ( "mellin1d" );
xspace = in.GetBool ( "xspace" );
yintervals = in.GetNumber ( "yintervals" );
yrule = in.GetNumber ( "yrule" );
qtintervals = in.GetNumber ( "qtintervals" );
qtrule = in.GetNumber ( "qtrule" );
abintervals = in.GetNumber ( "abintervals" );
abrule = in.GetNumber ( "abrule" );
vjphirule = in.GetNumber ( "vjphirule" );
zrule = in.GetNumber ( "zrule" );
xrule = in.GetNumber ( "xrule" );
ptbinwidth = in.GetBool ( "ptbinwidth" );
ybinwidth = in.GetBool ( "ybinwidth" );
force_binsampling = in.GetBool ( "force_binsampling" );
helicity = in.GetNumber ( "helicity" );
output_filename = in.GetString ( "output_filename" );
return ;
}
void settings::check_consistency(){
// additional conditions
if (order != 0 && order != 1 && order != 2 && order != 3)
throw invalid_argument("Invalid order, please select 0 (LL) 1 (NLL) 2 (NNLL) or 3 (NNNLL)");
if (nproc != 1 && nproc != 2 && nproc != 3)
throw invalid_argument("Wrong process, please select nproc = 1 (W+), 2 (W-), or 3(Z)");
if (order == 0)
{
doCT = false;
doVJ = false;
}
- if (order != 2)
+ if (order < 2)
{
doVJREAL = false;
doVJVIRT = false;
}
//In fixed order mode, a_param must be one
if (fixedorder == true)
{
//cout << "Asked for fixed order only predictions, enforce kmures = 1.0" << endl;
kmures = 1.0;
fmures = 1;
}
+ //If default values of blim for pdfs, sudakov and cexp are not set, take blim as default
+ if (blim_pdf == 0)
+ blim_pdf = blim;
+
+ if (blim_sudakov == 0)
+ blim_sudakov = blim;
+
+ if (blim_cexp == 0)
+ blim_cexp = blim;
+
if (evolmode > 4 || evolmode < 0)
{
cout << "wrong value for evolmode: available evolmodes: 0,1,2,3,4" << endl;
exit (-1);
}
//if (fmufac > 0 && evolmode < 3 && order > 0 && fixedorder == false)
// {
// //cannot use a dynamic muren, mufac, when the PDFs are converted from x-space to N-space at the factorisation scale
// cout << "At NLL and NNLL mufac = mll is possible only with evolmode = 3 or 4" << endl;
// exit (-1);
// }
//minimal b prescription in the complex plain available only for evolmode 0
//real axix minimal presxription works only at LL
//with modlog = false can run only the resummation term
if (PDFerrors == true && LHAPDFmember != 0)
{
cout << "Asked for PDFerrors, enforce LHAPDFmember = 0" << endl;
LHAPDFmember = 0;
}
if (qtcut <= 0 && xqtcut <= 0)
{
cout << "At least one between qtcut and xqtcut must be > 0" << endl;
exit (-1);
}
//Automatic selector of integration type
if (intDimBorn < 0)
if (BORNquad)
intDimBorn = 2;
else
intDimBorn = 4;
if (intDimRes < 0)
if (BORNquad)
if (opts.makecuts)
intDimRes = 2;
else
intDimRes = 1; // -> Check this works also when not fully integrated in rapidity!!!
else
intDimRes = 4;
//Here check also if yrange is above ymax to set mellin1d = true when makecuts is false
if (intDimCT < 0)
if (CTquad)
if (opts.makecuts) // || yrange below ymax
intDimCT = 2;
else
intDimCT = 1;
else
intDimCT = 8;
if (intDimVJ < 0)
if (VJquad)
if (opts.makecuts)
intDimVJ = 5;
else
intDimVJ = 3;
else
intDimVJ = 7;
- if (doVJ && intDimVJ < 7 && order == 2)
+ if (makecuts && doVJ && intDimVJ < 7 && order >= 2)
{
cout << "cannot perform quadrature integration for V+jet at NNLO" << endl;
exit (-1);
}
// resummation term integration dimension
if (intDimRes<4 && intDimRes>0){
resint1d = (intDimRes == 1);
resint2d = (intDimRes == 2);
resint3d = (intDimRes == 3);
resintvegas = false;
} else {
resint2d = false;
resint3d = false;
resintvegas = true;
}
// born term integration dimension
if (intDimBorn < 4 && intDimBorn>1){
bornint2d = (intDimBorn == 2);
bornintvegas4d = false;
bornintvegas6d = false;
} else {
bornint2d = false;
bornintvegas4d = (intDimBorn == 4);
bornintvegas6d = (intDimBorn > 5);
}
// counter term integration dimension
if (intDimCT<4 && intDimCT>0){
ctint1d = (intDimCT == 1);
ctint2d = (intDimCT == 2);
ctint3d = (intDimCT == 3);
ctintvegas6d = false;
ctintvegas8d = false;
} else {
ctint2d = false;
ctint3d = false;
ctintvegas6d = (intDimCT <= 6);
ctintvegas8d = (intDimCT > 6);
}
// V+J integration dimension
if (intDimVJ < 7 && intDimVJ > 2)
{
vjint3d = (intDimVJ == 3);
vjint5d = (intDimVJ == 5);
vjintvegas7d = false;
}
else
{
vjint3d = false;
vjint5d = false;
vjintvegas7d = true;
}
if (mellin1d && (!(resint2d || resint1d) || makecuts))
{
cout << "mellin1d option is possible only for 1d or 2d integration of resummed piece, no cuts on leptons, and integration between [-ymin,ymax]" << endl;
exit (-1);
}
if (mellinintervals*mellinrule >= 512)
{
cout << "mellinintervals*mellinrule should be less than 512 " << endl;
exit(-1);
}
if (makecuts && vjint3d)
{
cout << "Required cuts on the final state leptons, enforce 5D integration for V+J fixed order cross section" << endl;
vjint3d = false;
vjint5d = true;
}
if (makecuts && helicity >= 0)
{
cout << "Required cuts on the final state leptons, cannot calculate helicity cross sections, enforce helicity = -1" << endl;
helicity = -1;
}
if (opts_.approxpdf_ == 1)
{
cout << "DYRES-style approximate PDF requested, enforce vegas integration for the resummed cross section" << endl;
resint2d = false;
resint3d = false;
resintvegas = true;
}
// -- binning
// check bins size
if ( bins.qtbins .size() < 2) throw QuitProgram("Option `qtbins` needs at least 2 items ");
if ( bins.ybins .size() < 2) throw QuitProgram("Option `ybins` needs at least 2 items ");
if ( bins.mbins .size() < 2) throw QuitProgram("Option `mbins` needs at least 2 items ");
// check sorting
sort( bins.qtbins .begin () , bins.qtbins .end () ) ;
sort( bins.ybins .begin () , bins.ybins .end () ) ;
sort( bins.mbins .begin () , bins.mbins .end () ) ;
// set histogram bins
bins.hist_qt_bins = bins.qtbins ;
bins.hist_y_bins = bins.ybins ;
bins.hist_m_bins = bins.mbins ;
// integration boundaries
phasespace::setbounds(
bins. mbins .front() ,
bins. mbins .back() ,
bins. qtbins .front() ,
bins. qtbins .back() ,
bins. ybins .front() ,
bins. ybins .back()
);
return ;
}
void settings::dumpAll(){
printf("==Listing settings==\n");
bool print_inputs = true;
bool print_process_inputs = true;
bool print_masses = true;
if (print_process_inputs) {
printf("Input process settings:\n");
dumpD ( "sroot ", energy_ . sroot_ ) ;
dumpI ( "ih1 ", density_ . ih1_ ) ;
dumpI ( "ih2 ", density_ . ih2_ ) ;
dumpI ( "nproc ", nproc_ . nproc_ ) ;
}
if (print_inputs) {
printf("Input settings:\n");
dumpD ("g ", g_param_ . g_param_ ) ;
dumpI ("npff ",npff );
dumpD ("g1 ",g1 );
dumpD ("g2 ",g2 );
dumpD ("g3 ",g3 );
dumpD ("e ",e );
dumpD ("g0 ",g0 );
dumpD ("Q0 ",Q0 );
dumpB ("flavour_kt ",flavour_kt );
dumpD ("g1_uv ",g1_uv );
dumpD ("g1_us ",g1_us );
dumpD ("g1_dv ",g1_dv );
dumpD ("g1_ds ",g1_ds );
dumpD ("g1_ss ",g1_ss );
dumpD ("g1_ch ",g1_ch );
dumpD ("g1_bo ",g1_bo );
dumpD ("g1_gl ",g1_gl );
dumpI ( "order ", nnlo_ . order_ ) ;
dumpD ( "kmuren ", kmuren ) ;
dumpD ( "kmufac ", kmufac ) ;
dumpD ( "kmures ", kmures ) ;
dumpI ( "fmuren ", kmuren ) ;
dumpI ( "fmufac ", kmufac ) ;
dumpI ( "fmures ", kmures ) ;
dumpD ( "kmuc ", kmuc ) ;
dumpD ( "kmub ", kmub ) ;
dumpD ( "kmut ", kmut ) ;
dumpB ( "flavour_kt", flavour_kt );
dumpD( "blim ", blim ) ;
+ dumpD( "blim_pdf ", blim_pdf ) ;
+ dumpD( "blim_sudakov ", blim_sudakov ) ;
+ dumpD( "bstar_cexp ", blim_cexp ) ;
+ dumpB( "bstar_pdf ", bstar_pdf ) ;
+ dumpB( "bstar_sudakov ", bstar_sudakov ) ;
+ dumpB( "bstar_cexp ", bstar_cexp ) ;
dumpS("LHAPDFset ", LHAPDFset );
dumpI("LHAPDFmember ", LHAPDFmember );
dumpI("rseed ", rseed );
dumpI("ewscheme ", ewscheme );
dumpD("Gf" , Gf);
dumpD("zmass" , zmass);
dumpD("wmass" , wmass );
dumpD("xw" , xw);
dumpD("aemmz" , aemmz);
dumpD("zwidth" , zwidth);
dumpD("wwidth" , wwidth);
dumpB("runningwidth" , runningwidth);
dumpD( "Vud", Vud);
dumpD( "Vus", Vus);
dumpD( "Vub", Vub);
dumpD( "Vcd", Vcd);
dumpD( "Vcs", Vcs);
dumpD( "Vcb", Vcb);
dumpD( "Zuu", Zuu);
dumpD( "Zdd", Zdd);
dumpD( "Zss", Zss);
dumpD( "Zcc", Zcc);
dumpD( "Zbb", Zbb);
//dumpD("ylow ", ylow );
//dumpD("yhigh ", yhigh );
//dumpD("mlow ", mlow );
//dumpD("mhigh ", mhigh );
dumpB("damp", damp );
dumpD("dampk", dampk );
dumpD("dampdelta", dampdelta );
dumpI("dampmode", dampmode );
dumpD("qtcutoff", qtcutoff );
dumpB("modlog", modlog );
dumpD("xqtcut", xqtcut );
dumpD("qtcut", qtcut );
dumpB("useGamma ", useGamma );
dumpI("intDimRes ", intDimRes );
dumpB("resint2d ", resint2d );
dumpB("resint3d ", resint3d );
dumpB("resintvegas ", resintvegas );
dumpI("intDimBorn ", intDimBorn );
dumpB("bornint2d ", bornint2d );
dumpB("bornintvegas4d ", bornintvegas4d );
dumpB("bornintvegas6d ", bornintvegas6d );
dumpI("intDimCT ", intDimCT );
dumpB("ctint2d ", ctint2d );
dumpB("ctint3d ", ctint3d );
dumpB("ctintvegas6d ", ctintvegas6d );
dumpB("ctintvegas8d ", ctintvegas8d );
dumpI("intDimVJ ", intDimVJ );
dumpB("vjint3d ", vjint3d );
dumpB("vjint5d ", vjint5d );
dumpB("vjintvegas7d ", vjintvegas7d );
dumpB("fixedorder_only ", fixedorder );
dumpB("BORNquad ", BORNquad );
dumpB("CTquad ", CTquad );
dumpB("VJquad ", VJquad );
dumpB("doBORN ", doBORN );
dumpB("doCT ", doCT );
dumpB("doVJ ", doVJ );
dumpB("doVJREAL ", doVJREAL );
dumpB("doVJVIRT ", doVJVIRT );
dumpI("cubaverbosity ", cubaverbosity );
dumpI("cores ", cubacores );
dumpI("cubanbatch ", cubanbatch );
dumpI("niterBORN ", niterBORN );
dumpI("niterCT ", niterCT );
dumpI("niterVJ ", niterVJ );
dumpD("vegasncallsBORN ", vegasncallsBORN );
dumpD("vegasncallsCT ", vegasncallsCT );
dumpD("vegasncallsVJLO ", vegasncallsVJLO );
dumpD("vegasncallsVJREAL ", vegasncallsVJREAL );
dumpD("vegasncallsVJVIRT ", vegasncallsVJVIRT );
dumpB("vegascollect ", vegascollect );
dumpB("pcubature ", pcubature );
dumpD("relaccuracy ", relaccuracy );
dumpD("absaccuracy ", absaccuracy );
dumpD("costhmin ", costhmin );
dumpD("costhmax ", costhmax );
dumpB("makecuts ", makecuts );
dumpD("lptcut ", lptcut );
dumpD("lycut ", lycut );
dumpD("mtcut ", mtcut );
dumpD("cthCSmin ", cthCSmin );
dumpD("cthCSmax ", cthCSmax );
dumpD("etmisscut ", etmisscut );
dumpB("cubaint ", cubaint );
dumpB("trapezint ", trapezint );
dumpB("quadint ", quadint );
dumpI("suavepoints ", suavepoints );
dumpI("nphitrape ", nphitrape );
dumpI("phiintervals ", phiintervals );
dumpI("phirule ", phirule );
dumpI("ncstart ", ncstart );
dumpB("qtrec_naive ", qtrec_naive );
dumpB("qtrec_cs ", qtrec_cs );
dumpB("qtrec_kt0 ", qtrec_kt0 );
dumpB("timeprofile ", timeprofile );
dumpB("verbose ", verbose );
dumpB("gridverbose ", gridverbose );
dumpB("unicode ", unicode );
dumpB("silent ", silent );
dumpB("makehistos ", makehistos );
dumpB("texttable ", texttable );
dumpB("resumcpp ", resumcpp );
dumpB("ctcpp ", ctcpp );
dumpI("approxpdf ", opts_.approxpdf_ );
dumpI("pdfintervals ", opts_.pdfintervals_ );
dumpI("evolmode ", evolmode );
+ dumpB("mufevol ", mufevol );
+ dumpB("nocexp ", nocexp );
dumpI("bprescription ", bprescription );
dumpB("phibr ", phibr );
dumpB("bcf ", bcf );
dumpB("PDFerrors ", PDFerrors );
dumpI("mellininv ", mellininv );
dumpI("mellinintervals ", mellinintervals );
dumpI("mellinrule ", mellinrule );
dumpD("zmax ", zmax );
dumpD("cpoint ", cpoint );
dumpD("phi ", phi );
dumpD("mellincores ", mellincores );
dumpB("mellin1d ", mellin1d );
dumpB("xspace ", xspace );
dumpI("yintervals ", yintervals );
dumpI("yrule ", yrule );
dumpI("qtintervals ", qtintervals );
dumpI("qtrule ", qtrule );
dumpI("abintervals ", abintervals );
dumpI("abrule ", abrule );
dumpI("vjphirule ", vjphirule );
dumpI("zrule ", zrule );
dumpI("xrule ", xrule );
dumpB("ptbinwidth ", ptbinwidth );
dumpB("ybinwidth ", ybinwidth );
dumpB("mbinwidth ", mbinwidth );
dumpB("force_binsampling ", force_binsampling );
dumpI("helicity ", helicity );
dumpS("output_filename ", output_filename );
}
if (print_masses){
printf("Masses and EW constants:\n");
dumpD ( "md " , dymasses_ . md_ );
dumpD ( "mu " , dymasses_ . mu_ );
dumpD ( "ms " , dymasses_ . ms_ );
dumpD ( "mc " , dymasses_ . mc_ );
dumpD ( "mb " , dymasses_ . mb_ );
dumpD ( "mt " , dymasses_ . mt_ );
dumpD ( "mel " , dymasses_ . mel_ );
dumpD ( "mmu " , dymasses_ . mmu_ );
dumpD ( "mtau " , dymasses_ . mtau_ );
dumpD ( "hmass " , dymasses_ . hmass_ );
dumpD ( "hwidth " , dymasses_ . hwidth_ );
dumpD ( "wmass " , dymasses_ . wmass_ );
dumpD ( "wwidth " , dymasses_ . wwidth_ );
dumpD ( "zmass " , dymasses_ . zmass_ );
dumpD ( "zwidth " , dymasses_ . zwidth_ );
dumpD ( "twidth " , dymasses_ . twidth_ );
dumpD ( "mtausq " , dymasses_ . mtausq_ );
dumpD ( "mcsq " , dymasses_ . mcsq_ );
dumpD ( "mbsq " , dymasses_ . mbsq_ );
/*
dumpD ( "Gf_inp " , ewinput_ . Gf_inp_ );
dumpD ( "aemmz_inp " , ewinput_ . aemmz_inp_ );
dumpD ( "xw_inp " , ewinput_ . xw_inp_ );
dumpD ( "wmass_inp " , ewinput_ . wmass_inp_ );
dumpD ( "zmass_inp " , ewinput_ . zmass_inp_ );
*/
dumpD ( "Vud " , cabib_ . Vud_ );
dumpD ( "Vus " , cabib_ . Vus_ );
dumpD ( "Vub " , cabib_ . Vub_ );
dumpD ( "Vcd " , cabib_ . Vcd_ );
dumpD ( "Vcs " , cabib_ . Vcs_ );
dumpD ( "Vcb " , cabib_ . Vcb_ );
dumpD ( "epinv " , epinv_ . epinv_ );
dumpD ( "epinv2 " , epinv2_ . epinv2_ );
}
printf("end of setting dump.\n\n");
}
void settings::dumpI( string var,int val){
printf( " %s = %d\n", var.c_str(), val);
}
void settings::dumpD( string var,double val){
printf( " %s = %f\n", var.c_str(), val);
}
void settings::dumpS( string var,string val){
printf( " %s = %s\n", var.c_str(), val.c_str());
}
void settings::dumpB( string var,bool val){
printf( " %s = %s\n", var.c_str(), val ? "true" : "false" );
}
vector<string> settings::Tokenize(string val,char Delim){
vector<string> vec;
size_t pos = 0;
string tmp;
while (!val.empty()){
// find delim
pos = val.find_first_of(Delim);
if (pos==string::npos){
// last item
tmp = val;
val.clear();
} else {
// get substring
tmp = val.substr(0,pos);
val = val.substr(pos+1);
}
// add to vector
if (!tmp.empty()) vec.push_back(tmp);
}
return vec;
}
bool settings::IsNumber(const string &s) {
double dummy;
try {
dummy = stod(s);
return true;
} catch (const std::exception &e){
return false;
}
}
void settings::parse_binning(string name, vector<double> &bins, po::Options &args){
if (args.count(name)) {
string e("Unsupported value of "+name+" : need 3 numbers seperated by comma: 'N,lo,hi' ");
string val=args[name.c_str()].as<string>();
ToLower(val);
vector<string> vec = Tokenize(val);
// check value
if (vec.size()!=3) throw QuitProgram(e+" not 3 numbers.");
for (auto s : vec) if (!IsNumber(s)) throw QuitProgram(e+" not numbers.");
// retrieve N,lo,hi
int N = stod(vec[0]);
double lo = stod(vec[1]);
double hi = stod(vec[2]);
if (lo>hi) throw QuitProgram(e+" lo is more then hi.");
if (N<1) throw QuitProgram(e+" N is not at least 1.");
// make binning
bins.clear();
//loop with double has problems for equality test (try --ybins 25,0,5)
//for (double loedge=lo; loedge<=hi; loedge+=(hi-lo)/double(N)) bins.push_back(loedge);
for (int i=0; i <= N; i++) bins.push_back(lo+i*(hi-lo)/double(N));
}
}
void binning::GetBins(string name, vector<double> &vec){
vec.clear();
// CLI defined
if (name == "qt" && qtbins .size()!=0) { vec = qtbins; return; }
if (name == "y" && ybins .size()!=0) { vec = ybins; return; }
if (name == "m" && mbins .size()!=0) { vec = mbins; return; }
// input file defined
name+="_bins";
in.GetVectorDouble(name.c_str(), vec);
if ( vec.size() < 2) throw QuitProgram( "Option `" +name+ "` needs at least 2 items ");
}
void binning::readfromfile(const string fname){
in.parse_file(fname);
in.GetVectorDouble("qt_bins" , qtbins );
in.GetVectorDouble("y_bins" , ybins );
in.GetVectorDouble("m_bins" , mbins );
hist_qt_bins .clear();
hist_y_bins .clear();
hist_m_bins .clear();
return;
}
// InputParser definitions
//
InputParser::InputParser( string _filename, string _charset, string _white):
filename ( _filename ),
Ccommnt ( _charset[0] ),
Cassign ( _charset[1] ),
CopenAr ( _charset[2] ),
CdeliAr ( _charset[3] ),
CclosAr ( _charset[4] ),
Swhite ( _white )
{
// parse default
try {
parse_file((string)SHAREDIR+"/default.in");
} catch (invalid_argument &e1){
try {
parse_file("default.in");
} catch (invalid_argument &e2){
throw invalid_argument( string( "Can not load default settings:\n")
+ string(e1.what()) + string("\n")
+ string(e2.what()) + string("\n")
);
}
}
// parse custom settings
if (filename.compare("")!=0) parse_file(filename);
}
InputParser::~InputParser(){
}
double InputParser::GetNumber(string name){
has_key(name);
string val = data[name];
try { return stod(val);
} catch (const std::exception &e){
printf("Cannot read option '%s' with value '%s' as number.\n",name.c_str(), val.c_str() );
throw e;
}
}
string InputParser::GetString(string name){
has_key(name);
string val = data[name];
trim(val);
return val;
}
bool InputParser::GetBool(string name){
has_key(name);
string val = data[name];
// lower case
std::transform(val.begin(), val.end(), val.begin(), ::tolower);
bool result = (val.compare(0,4,"true") == 0);
//printf(" InputParser::GetBool processing key `%s` with value `%s` as result `%d`\n", name.c_str(), val.c_str(), result);
return result;
}
void InputParser::GetVectorDouble(string name, vector<double> &vec){
has_key(name);
vec.clear();
string val = data[name];
// has open/close array
size_t strBegin = val.find(CopenAr,0);
size_t strEnd = val.find(CclosAr,0);
// we need both open and close otherwise is something wrong
if (strBegin==string::npos || strEnd==string::npos) throw invalid_argument("Missing open/close character.");
// parse what is between them
val = val.substr(strBegin+1,strEnd-strBegin-1);
size_t pos = 0;
string tmp;
while (!val.empty()){
// find delim
pos = val.find(CdeliAr);
if (pos==string::npos){
// last item
tmp = val;
val.clear();
} else {
// get substring
tmp = val.substr(0,pos);
val = val.substr(pos+1);
}
// get substring
trim(tmp);
// expecting number otherwise exception
if (!tmp.empty()) vec.push_back(stod(tmp));
}
return;
}
void InputParser::parse_file(const string fname){
ifstream fstrm(fname.c_str());
if ( ! fstrm.good() ) throw invalid_argument( string("Uknown file name ") + fname );
string line;
// line by line
while(getline(fstrm,line)) {
size_t pos = 0;
// check for comments
pos = line.find(Ccommnt,pos);
if (pos != string::npos) line = line.substr(0,pos);
// split string by delim
pos = line.find(Cassign,0);
if (pos == string::npos || line.empty()) continue;
string key = line.substr(0,pos);
string val = line.substr(pos+1);
// remove leading / trailing whitespace
trim(key);
trim(val);
/// @todo: Multiline setting. If has CopenAr load until CclosAr
// save to map
data[key]=val;
}
fstrm.close();
return;
}
void InputParser::trim(string & str){
size_t strBegin = str.find_first_not_of(Swhite);
if (strBegin == std::string::npos){
str="";
return ; // no content
}
size_t strEnd = str.find_last_not_of(Swhite);
size_t strRange = strEnd - strBegin + 1;
str = str.substr(strBegin, strRange);
return;
}
void InputParser::has_key(const string key){
if ( data.count(key) == 0 ){
string msg = "Missing setting with name '";
msg += key;
msg += "'";
throw invalid_argument(msg.c_str());
}
return;
}
diff --git a/src/settings.h b/src/settings.h
index 229d4e2..2144dc4 100644
--- a/src/settings.h
+++ b/src/settings.h
@@ -1,344 +1,357 @@
#ifndef settings_h
#define settings_h
#include <vector>
using namespace std;
#include <algorithm>
#include <map>
#include <string>
#include <fstream>
class InputParser {
public:
// constructor
InputParser( string _filename = "", string _charset="#=[ ]", string _white=" \t");
~InputParser();
// getters
double GetNumber(string name);
string GetString(string name);
bool GetBool(string name);
void GetVectorDouble(string name, vector<double> &vec);
void parse_file(const string fname);
private :
// functions
void trim(string & str);
void has_key(const string key);
// data members
map<string,string> data; ///< Data table.
string filename; ///< Input file.
// parser characters
char Ccommnt;
char Cassign;
char CopenAr;
char CclosAr;
char CdeliAr;
string Swhite;
};
// Simple program quiting exception
struct QuitProgram : public std::runtime_error { QuitProgram(string msg) : std::runtime_error (msg) {}; };
// Forward declaration
namespace cxxopts{
class Options;
}
namespace po=cxxopts; // inspired by po = boost::program_options
class settings
{
public:
settings() {};
void parse_options(int argc, char * argv[]);
void readfromfile(const string fname);
void check_consistency();
void parse_binning(string name, vector<double> &vec, po::Options &args);
// private:
void dumpAll();
void dumpI(string var, int val );
void dumpD(string var, double val );
void dumpS(string var, string val );
void dumpB(string var, bool val );
// string helpers
void ToLower(string &val){std::transform(val.begin(), val.end(), val.begin(), ::tolower);}
void ToUpper(string &val){std::transform(val.begin(), val.end(), val.begin(), ::toupper);}
vector<string> Tokenize(string val, char Delim=',');
// is number: http://stackoverflow.com/a/16575025
bool IsNumber(const string &s);
//process settings
double sroot;
int ih1;
int ih2;
int nproc;
//resummation or fixed order switch
bool fixedorder;
int order;
//Non-perturbative form factor
int npff;
double g1,g2,g3; //Gaussian
double e; //Exponential
double g0; //Collins-Rogers
double Q0; //reference mass
//Flavour dependent g1
bool flavour_kt;
double g1_uv = 0.5;
double g1_us = 0.5;
double g1_dv = 0.5;
double g1_ds = 0.5;
double g1_ss = 0.5;
double g1_ch = 0.5;
double g1_bo = 0.5;
double g1_gl = 0.5;
//PDF settings
string LHAPDFset ;
int LHAPDFmember ;
bool externalpdf;
//functional forms of the QCD scales
int fmures;
int fmuren;
int fmufac;
//scale factors for the QCD scales
double kmures;
double kmuren;
double kmufac;
//scale factors for the matching scales
double kmuc;
double kmub;
double kmut;
// double a_param;
//EW parameters
int ewscheme;
double Gf, zmass, wmass;
double xw, aemmz;
double zwidth, wwidth;
//Fixed width to running width translation
bool runningwidth;
//CKM matrix
double Vud, Vus, Vub;
double Vcd, Vcs, Vcb;
//Z/gamma* coupling
double Zuu, Zdd, Zss, Zcc, Zbb;
//resonance mass and width (used for breit wigner unweighting)
double rmass, rwidth;
// photon switch
bool useGamma;
//integration boundaries
//double ylow;
//double yhigh;
//double mlow;
//double mhigh;
//Resummation damping
bool damp;
double dampk, dampdelta;
int dampmode;
//Resummation cutoff
double qtcutoff;
//Modified logarithms
bool modlog;
//qtcut
double xqtcut, qtcut;
//integration settings
int rseed;
//dimension of integration for the resummed part
int intDimRes;
bool resint1d, resint2d, resint3d, resintvegas;
//dimension of integration for the born configuration
int intDimBorn;
bool bornint2d, bornintvegas4d, bornintvegas6d;
//type of integration for the counterterm
int intDimCT;
bool ctint1d, ctint2d, ctint3d, ctintvegas6d, ctintvegas8d;
//type of integration for the V+j at LO
int intDimVJ;
bool vjint3d, vjint5d, vjintvegas7d;
//type of integration, automatic selector
bool BORNquad, CTquad, VJquad;
//term switches
bool doBORN;
bool doCT;
bool doVJ;
bool doVJREAL, doVJVIRT;
//Cuba settings
int cubaverbosity;
int cubacores;
int cubanbatch;
int niterBORN;
int niterCT;
int niterVJ;
int vegasncallsBORN ;
int vegasncallsCT ;
int vegasncallsVJLO ;
long long int vegasncallsVJREAL ;
int vegasncallsVJVIRT ;
bool vegascollect;
//cubature settings
bool pcubature;
double relaccuracy;
double absaccuracy;
//costh boundaries
double costhmin, costhmax;
//lepton fiducial cuts
bool makecuts;
double lptcut, lycut; //charged leptons
double mtcut, etmisscut;
double lepptcut, lepycut, alpptcut, alpycut; //lepton and antilepton
double lcptcut, lcymin, lcymax, lfptcut, lfymin, lfymax; //lc and lf are absolute-rapidity-ordered leptons
double cthCSmin, cthCSmax;
//integration types and settings for costh phi_lep phase space
bool cubaint;
int suavepoints;
bool trapezint;
int nphitrape;
int ncstart;
//quadrature rule in phi_lep
bool quadint;
int phiintervals;
int phirule;
//settings for Bessel integration
double bintaccuracy;
//b-space prescription
int bprescription;
//blim parameter of the bstar prescription (acts as an IR cut-off)
double blim;
+ double blim_pdf, blim_sudakov, blim_cexp;
+ //force bstar prescription
+ bool bstar_pdf, bstar_sudakov, bstar_cexp;
+
//arg(z) in the complex plane for the minimal prescription
double phibr;
//select the point bc, where the integration contour is bended in the complex plane, as a fraction of b_L = ... (Landau singularity)
double bcf;
//settings for Mellin integration
int mellininv;
int mellinintervals;
int mellinrule;
double zmax;
double cpoint;
double phi;
int mellincores;
bool mellin1d;
bool xspace;
//settings for x-to-N Mellin transform
int pdfrule;
//settings for rapidity integration in 2D resummed piece
int yintervals;
int yrule;
//settings for qt integration in 2D counter term
int qtintervals;
int qtrule;
//settings for alfa beta scaled-PDF integration in counter term and born fixed order
int abintervals;
int abrule;
//settings for the phi integration in the V+J 5d LO term when makecuts is false
int vjphirule;
//settings for the z1, z2 integration in the V+J 3d NLO singular term
int zrule;
//settings for the x integration in the V+J 3d delta term
int xrule;
//qt-recoil prescriptions
bool qtrec_naive, qtrec_cs, qtrec_kt0;
// PDF errors
bool PDFerrors;
int totpdf;
//debug settings
bool timeprofile;
bool verbose;
bool gridverbose;
//output settings
bool texttable;
bool redirect;
bool unicode;
bool silent;
bool makehistos;
string output_filename; // Output Filenames
//resummed code in C++
bool resumcpp = true; //use C++ code for resummation
//counter term code in C++
bool ctcpp = true; //use C++ code for the counter term
- //dyres or pegasus PDF evolution
+ //dyres, pegasus,lhapdf PDF evolution
int evolmode;
+ //Evolve PDFs from mufac instead of mures
+ bool mufevol;
+
+ //swicth off C exponentiation
+ bool nocexp;
+
+ //sum all logs before exponentiation
+ bool sumlogs;
+
//bin width normalisation
bool ptbinwidth, ybinwidth, mbinwidth;
// Force to loop over all bins even you have all Vegas integrands
bool force_binsampling = false;
// Calculate helicity cross sections
int helicity;
};
class binning
{
public:
binning() {};
void readfromfile(const string fname);
void GetBins(string name,vector<double> &bins);
// private:
vector <double> qtbins;
vector <double> ybins;
vector <double> mbins;
vector <double> hist_qt_bins;
vector <double> hist_y_bins;
vector <double> hist_m_bins;
private:
InputParser in;
};
extern settings opts;
extern binning bins;
#endif

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 4:33 PM (22 h, 10 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023462
Default Alt Text
(67 KB)

Event Timeline