diff --git a/tools/plotting/fnlo-py-print.py b/tools/plotting/fnlo-py-print.py index 2d4cc6719713cbd36ca56ef650e8ecb6901ca0ba..fbbddb7151b303fc21a67f94571811533f625f58 100755 --- a/tools/plotting/fnlo-py-print.py +++ b/tools/plotting/fnlo-py-print.py @@ -1,4 +1,4 @@ -#!/usr/bin/env python3 +#!/usr/bin/env -S python3 -u #-*- coding:utf-8 -*- # ######################################################################## @@ -34,22 +34,25 @@ class SplitArgs(argparse.Action): # Some global definitions -_fntrans = str.maketrans({'[': '', ']': '', '(': '', ')': '', ',': ''}) # Filename translation table -_sntrans = str.maketrans({'[': '', ']': '', '(': '', ')': '', ',': '', '/': '÷'}) # Scalename translation table -_formats = {'eps': 0, 'pdf': 1, 'png': 2, 'svg': 3} +#_fntrans = str.maketrans({'[': '', ']': '', '(': '', ')': '', ',': ''}) # Filename translation table +#_sntrans = str.maketrans({'[': '', ']': '', '(': '', ')': '', ',': '', '/': '÷'}) # Scalename translation table +#_formats = {'eps': 0, 'pdf': 1, 'png': 2, 'svg': 3} _text_to_order = {'LO': 0, 'NLO': 1, 'NNLO': 2} _order_to_text = {0: 'LO', 1: 'NLO', 2: 'NNLO'} -_order_color = {'LO': 'g', 'NLO': 'b', 'NNLO': 'r'} -_order_symbol = {'LO': '^', 'NLO': 's', 'NNLO': 'o'} -_colors = ['tab:orange', 'tab:green', 'tab:purple', 'tab:blue', 'tab:brown'] -_symbols = ['s', 'X', 'o', '^', 'v'] -_hatches = ['', '//', '\\', '|', '-'] -_scale_to_text = {0: 'kScale1', 1: 'kScale2', 2: 'kQuadraticSum', 3: 'kQuadraticMean', 4: 'kQuadraticSumOver4', - 5: 'kLinearMean', 6: 'kLinearSum', 7: 'kScaleMax', 8: 'kScaleMin', 9: 'kProd', - 10: 'kS2plusS1half', 11: 'kPow4Sum', 12: 'kWgtAvg', 13: 'kS2plusS1fourth', 14: 'kExpProd2', 15: 'kExtern'} -_pdfbasenames = ['ABM11', 'ABMP15', 'ABMP16', 'CJ12', 'CJ15', 'CT10', 'CT14', 'HERAPDF20', 'JR14', - 'MMHT2014', 'MSTW2008', 'NNPDF23', 'NNPDF30', 'NNPDF31'] -_debug = False +#_order_color = {'LO': 'g', 'NLO': 'b', 'NNLO': 'r'} +#_order_symbol = {'LO': '^', 'NLO': 's', 'NNLO': 'o'} +#_colors = ['tab:orange', 'tab:green', 'tab:purple', 'tab:blue', 'tab:brown'] +#_symbols = ['s', 'X', 'o', '^', 'v'] +#_hatches = ['', '//', '\\', '|', '-'] +#_scale_to_text = {0: 'kScale1', 1: 'kScale2', 2: 'kQuadraticSum', 3: 'kQuadraticMean', 4: 'kQuadraticSumOver4', +# 5: 'kLinearMean', 6: 'kLinearSum', 7: 'kScaleMax', 8: 'kScaleMin', 9: 'kProd', +# 10: 'kS2plusS1half', 11: 'kPow4Sum', 12: 'kWgtAvg', 13: 'kS2plusS1fourth', 14: 'kExpProd2', 15: 'kExtern'} +#_pdfbasenames = ['ABM11', 'ABMP15', 'ABMP16', 'CJ12', 'CJ15', 'CT10', 'CT14', 'HERAPDF20', 'JR14', +# 'MMHT2014', 'MSTW2008', 'NNPDF23', 'NNPDF30', 'NNPDF31'] +#_debug = False +# TODO: How can this be done using fastNLO directly? +_verb_to_enum = {'SILENT': 1000, 'ERROR': 2, 'WARNING': 1, 'INFO': 0, 'MANUAL': -1000, 'DEBUG': -1000} + ######################################################################## @@ -77,47 +80,64 @@ def main(): parser.add_argument('-v', '--verbosity', default='WARNING', type=str, help="Set fastNLO output verbosity to one of SILENT, ERROR, WARNING, INFO, MANUAL, DEBUG") - # Print header - print("\n###########################################################################################") - print("# fnlo-py-print:") - print("# Read fastNLO grids and print cross section & uncertainty") - print("###########################################################################################\n") - # Parse arguments args = vars(parser.parse_args()) + # Verbosity + verb = args['verbosity'] + eVerb = _verb_to_enum[verb] + # fVerbosity = toVerbosity()[VerbosityLevel] + SetGlobalVerbosity(eVerb) + + # Print header + if eVerb < 1: + print(" ##################################################################################") + print(" # [fnlo-py-print] Program to read fastNLO tables and write out") + print(" # [fnlo-py-print] QCD cross sections with uncertainties") + print(" #---------------------------------------------------------------------------------") + print(" # [fnlo-py-print] For an explanation of command line arguments type:") + print(" # [fnlo-py-print] ./fnlo-py-print.py -h") + print(" ##################################################################################") + print("") + # List of table names files = args['table'] - print('\n') - print('[fnlo-py-print]: Analysing table list:') + if eVerb < 1000: + print(" ##################################################################################") + print(" # [fnlo-py-print] Iterating over table list:") for file in files: - print('[fnlo-py-print]: ', file) + if eVerb < 1000: + print(' # [fnlo-py-print] ', file) # PDF set name pdfset = os.path.basename(args['pdfset']) - print('[fnlo-py-print]: Using PDF set', pdfset) + if eVerb < 1000: + print(' # [fnlo-py-print] Using PDF set:', pdfset) # Orders to be shown iorders = [] iordmin = _text_to_order['LO'] iordmax = _text_to_order['NNLO'] if args['order'] is None: - print('[fnlo-py-print]: Evaluate table up to highest available order.') + if eVerb < 1000: + print(' # [fnlo-py-print] Evaluate tables up to highest available order.') else: for ord in args['order']: if ord in _text_to_order: iorders.append(_text_to_order[ord]) else: - print('[fnlo-py-print]: Illegal order specified, aborted!') - print('[fnlo-py-print]: Order list:', args['order']) + print(' # [fnlo-py-print] Illegal order specified, aborted!') + print(' # [fnlo-py-print] Order list:', args['order']) exit(1) iordmin = min(iorders) iordmax = max(iorders) - print('[fnlo-py-print]: Evaluate table up to order(s)', args['order']) + if eVerb < 1000: + print(' # [fnlo-py-print] Evaluate tables up to order(s):', args['order']) + print(" ##################################################################################") # Type of uncertainty (None, scale variation [2P,6P], PDF [L6], statistical [ST]) unc_type = None - unc_style = None + unc_style = fastnlo.kScaleNone unc_label = '' if args['uncertainty']: unc_type = args['uncertainty'] @@ -131,19 +151,16 @@ def main(): unc_style = fastnlo.kLHAPDF6 unc_label = 'PDF uncertainty (LHAPDF)' elif unc_type == 'ST': - unc_style = fastnlo.kAddStat + unc_style = fastnlo.kStatInt unc_label = 'numerical uncertainty (theory)' else: - print('[fnlo-py-print]: Illegal uncertainty specified, aborted!') - print('[fnlo-py-print]: Uncertainty:', args['uncertainty']) + print(' # [fnlo-py-print] Illegal uncertainty specified, aborted!') + print(' # [fnlo-py-print] Uncertainty:', args['uncertainty']) exit(11) # Scale choice scale_choice = args['scale'] - # Verbosity - verb = args['verbosity'] - # Loop over table list for table in files: # Table name @@ -156,32 +173,35 @@ def main(): elif tablename.endswith('.tab'): tablename = tablename.replace('.tab', '', 1) else: - print('[fnlo-py-print]: Error! Wrong extension for table: ', table) + print(' # [fnlo-py-print] Error! Wrong extension for table: ', table) exit(1) - print('[fnlo-py-print]: Analysing table: ', table) + if eVerb < 1000: + print("") + print(" ##################################################################################") + print(' # [fnlo-py-print] Evaluating table: ', table) + print(" ##################################################################################") + print("") ###################### Start EVALUATION with fastNLO library ################################################### - # SetGlobalVerbosity(0) # Does not work since changed to default in the following call # Take the general information (bin_bounds, labels, order_existence, etc.) from given pdfset. - print("HHHHHHHHHHHHHHHHH: {}".format(verb)) - fnlo = fastnlo.fastNLOLHAPDF(table, args['pdfset'], args['member'], args['verbosity']) + fnlo = fastnlo.fastNLOLHAPDF(table, args['pdfset'], args['member']) # Get labeling for the x-axis # Dimensionality of the table: ndim = fnlo.GetNumDiffBin() - if verb: - print('\n') - print('[fnlo-py-print]: Table Dimensions: ', ndim) + if eVerb < 0: + print('') + print(' # [fnlo-py-print] Table Dimensions: ', ndim) # Labels of all the dimensions: labels = fnlo.GetDimLabels() - if verb: - print('[fnlo-py-print]: Labels:', labels) + if eVerb < 0: + print(' # [fnlo-py-print] Labels:', labels) # x label of first dimension from table: xlabel = fnlo.GetDimLabel(0) - if verb: - print('[fnlo-py-print]: x-label:', xlabel) + if eVerb < 0: + print(' # [fnlo-py-print] x-label:', xlabel) # Generic y label ylabel = '$\sigma \pm \Delta\sigma(\mu_R,\mu_F)$' @@ -189,22 +209,21 @@ def main(): # Creating x-axis bin_bounds = np.array(fnlo.GetObsBinsBounds(0)) # bin_bounds = np.array(fnlo.GetDim1BinBounds(0)) Mods needed for 2D tables ... - if verb: - print('[fnlo-py-print]: bin_bounds.T: \n', bin_bounds.T, '\n') - print('[fnlo-py-print]: bin_bounds.flatten()', - bin_bounds.flatten(), '\n') + if eVerb < 0: + print(' # [fnlo-py-print] bin_bounds.T: \n', bin_bounds.T, '\n') + print(' # [fnlo-py-print] bin_bounds.flatten()', bin_bounds.flatten(), '\n') x_axis = (bin_bounds.T[0]+bin_bounds.T[1]) / 2. # this is a list of bin centers xmin = 0.95*min(bin_bounds.ravel()) xmax = 1.05*max(bin_bounds.ravel()) - if verb: - print('[fnlo-py-print]: xmin=%s, xmax=%s. \n' % (xmin, xmax)) + if eVerb < 0: + print(' # [fnlo-py-print] xmin=%s, xmax=%s. \n' % (xmin, xmax)) # Preparing x-errors (via bin_bounds) --> x_errors[0, :] are initially negative (positive via -1*), x_errors[1, :] positive x_errors = np.array( [-1*(bin_bounds.T[0]-x_axis), bin_bounds.T[1]-x_axis]) - if verb: - print('[fnlo-py-print]: \n x_errors: ', x_errors, '\n') + if eVerb < 0: + print(' # [fnlo-py-print] \n x_errors: ', x_errors, '\n') # Check existence of orders in table lflex = fnlo.GetIsFlexibleScaleTable() @@ -219,10 +238,8 @@ def main(): max_order = i if not lflex: if scale_choice != 0: - print( - '[fnlo-py-print]: Invalid choice of scale = ', scale_choice, ' Aborted!') - print( - '[fnlo-py-print]: For fixed-scale tables only the default=0 is allowed.') + print(' # [fnlo-py-print] Invalid choice of scale = ', scale_choice, ' Aborted!') + print(' # [fnlo-py-print] For fixed-scale tables only the default=0 is allowed.') exit(1) else: scale_name = fnlo.GetScaleDescription(i, 0) @@ -236,8 +253,8 @@ def main(): '_'+scl0+'_'+scl1 if cnt_order == i-1: cnt_order += 1 - if verb: - print('[fnlo-py-print]: Table has continuous orders up to', + if eVerb < 0: + print(' # [fnlo-py-print] Table has continuous orders up to', cnt_order, 'and a maximal order of', max_order) # If previously undefined set iordmax to maximum found in table @@ -245,8 +262,8 @@ def main(): iordmax = max_order if iordmax > cnt_order: - print('[fnlo-py-print]: Invalid choice of orders. Aborted!') - print('[fnlo-py-print]: Highest order requested is', + print(' # [fnlo-py-print] Invalid choice of orders. Aborted!') + print(' # [fnlo-py-print] Highest order requested is', _order_to_text[iordmax], 'but continuous orders are available only up to', _order_to_text[cnt_order]) exit(1) @@ -257,23 +274,23 @@ def main(): else: order_list = args['order'] - print('[fnlo-py-print]: List of requested orders:', order_list) + if eVerb < 1: + print(' # [fnlo-py-print] List of requested orders:', order_list) # For flexible-scale tables set scale to user choice (default is 0) if lflex: - print( - '[fnlo-py-print]: Setting requested scale choice for flexible-scale table:', scale_choice) + if eVerb < 0: + print(' # [fnlo-py-print] Setting requested scale choice for flexible-scale table:', scale_choice) fnlo.SetMuRFunctionalForm(scale_choice) fnlo.SetMuFFunctionalForm(scale_choice) else: if scale_choice == 0: - print( - '[fnlo-py-print]: Evaluating fixed-scale table. Scale choice must be', scale_choice) + if eVerb < 0: + print(' # [fnlo-py-print] Evaluating fixed-scale table. Scale choice must be', scale_choice) else: - print( - '[fnlo-py-print]: No scale choice possible for fixed-scale table. Aborted!') - print('[fnlo-py-print]: scale_choice = ', scale_choice) + print(' # [fnlo-py-print] No scale choice possible for fixed-scale table. Aborted!') + print(' # [fnlo-py-print] scale_choice = ', scale_choice) exit(1) # Now evaluate fastNLO table having a look at the requested uncertainty @@ -286,53 +303,37 @@ def main(): fnlo.SetContributionON(fastnlo.kFixedOrder, j, True) else: fnlo.SetContributionON(fastnlo.kFixedOrder, j, False) - if verb: - print('[fnlo-py-print]: \n') - print('[fnlo-py-print]: Calculate XS for order: %s' % n, '\n') - print( - '[fnlo-py-print]: ---- ---- ---- ---- ---- ---- ---- ----') + if eVerb < 0: + print(' # [fnlo-py-print] \n') + print(' # [fnlo-py-print] Calculate XS for order: %s' % n, '\n') + print(' # [fnlo-py-print] ---- ---- ---- ---- ---- ---- ---- ----') fnlo.CalcCrossSection() xs_list.append(fnlo.GetCrossSection()) ### Get requested uncertainty ### - if verb: - print('[fnlo-py-print]: Used scale factor MuF: ', fnlo.GetScaleFactorMuF()) - print('[fnlo-py-print]: Used scale factor MuR: ', fnlo.GetScaleFactorMuR(), '\n') - print('[fnlo-py-print]: Calculate uncertainty for this central scale.\n') # RELATIVE uncertainty for chosen uncertainty type and style # Calculate this for all accessible orders just in case + if eVerb < 0: + print(' # [fnlo-py-print] Used scale factor MuF: ', fnlo.GetScaleFactorMuF()) + print(' # [fnlo-py-print] Used scale factor MuR: ', fnlo.GetScaleFactorMuR(), '\n') + print(' # [fnlo-py-print] Calculate uncertainty for this central scale.\n') + elif eVerb < 1000: + print(" ##################################################################################") + print("") + # Start of fastNLO print out lnorm = False iprint = 1 if unc_type == '2P' or unc_type == '6P': - #unc_name = " # Relative scale uncertainties (" + unc_type + ")" - #fnlo.PrintScaleUncertaintyVec(unc_style,unc_name) - rel_unc_item = np.array(fnlo.GetScaleUncertaintyVec(unc_style,lnorm,iprint)) + rel_unc_item = np.array(fnlo.GetXsUncertaintyVec(unc_style,lnorm,iprint)) elif unc_type == 'L6': - #unc_name = " # Relative PDF Uncertainties (" + n + " " + pdfset + " " + unc_type + ")" - #fnlo.PrintPDFUncertaintyVec(unc_style,unc_name) - rel_unc_item = np.array(fnlo.GetPDFUncertaintyVec(unc_style,lnorm,iprint)) + rel_unc_item = np.array(fnlo.GetXsUncertaintyVec(unc_style,lnorm,iprint)) elif unc_type == 'ST': - #unc_name = " # Relative statistical uncertainties (" + unc_type + ")"; - #fnlo.PrintAddUncertaintyVec(unc_style,unc_name) - rel_unc_item = np.array(fnlo.GetAddUncertaintyVec(unc_style,lnorm,iprint)) + rel_unc_item = np.array(fnlo.GetXsUncertaintyVec(unc_style,lnorm,iprint)) else: - #unc_name = " # Without uncertainties"; - #fnlo.PrintScaleUncertaintyVec('NN',unc_name) - rel_unc_item = np.array(fnlo.GetAddUncertaintyVec(fastnlo.kAddNone,lnorm,iprint)) - + rel_unc_item = np.array(fnlo.GetXsUncertaintyVec(fastnlo.kScaleNone,lnorm,iprint)) - exit(0) rel_unc_list.append(rel_unc_item) - if verb: - print('[fnlo-py-print]: \n') - print('[fnlo-py-print]: Relative uncertainty in %s: \n' % n) - # 3 entries: central value (xs), unc_low, unc_high - print(rel_unc_item, '\n') - print( - '---------------------------------------------------------------------------------------') - print( - '---------------------------------------------------------------------------------------') xs_all = np.array(xs_list) rel_unc = np.array(rel_unc_list) @@ -344,44 +345,10 @@ def main(): # rel_unc[0,2,:] means rel. uncertainty downwards (in LO) ######### - if verb: - print('[fnlo-py-print]: Cross section xs_all uses ', order_list) - - if verb: - print (xs_all, '\n \n') - - # ABSOLUTE scale uncertainty - # length of axis 0 in xs_all equals number of orders - num_orders = np.size(xs_all, 0) - abs_scale_unc = np.empty([num_orders, 2, len(x_axis)]) - for k in range(0, len(xs_all)): - abs_scale_unc[k, 0, :] = rel_unc[k, 2, :] * \ - xs_all[k] # absolute uncertainties downwards (low) - abs_scale_unc[k, 1, :] = rel_unc[k, 1, :] * \ - xs_all[k] # absolute uncertainties upwards (high) - - if verb: - print( - '[fnlo-py-print]: Absolute Scale uncertainties downwards, upwards (order by order): \n') - print(abs_scale_unc, '\n') - - ############################## Do the printing #################################################### - - print(' ##################################################################################') - print(' # [fnlo-py-print]: Evaluating uncertainties') - print(' ##################################################################################') - print(' #=================================================================================') - print(' # {}'.format('BLUBBER')) - print(' #---------------------------------------------------------------------------------') - print(' # bin cross_section lower_uncertainty upper_uncertainty') - print(' #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -') - - - stop_time = timeit.default_timer() timediff = stop_time-start_time - print('fnlo-py-print: Elapsed time: %s sec = %s min' % - (timediff, round(timediff/60, 2))) + if eVerb < 1: + print(' # [fnlo-py-print] Elapsed time: %s sec = %s min' % (timediff, round(timediff/60, 2))) if __name__ == '__main__': diff --git a/v2.5/toolkit/fastnlotoolkit/fastNLOLHAPDF.cc b/v2.5/toolkit/fastnlotoolkit/fastNLOLHAPDF.cc index 2ff3e9052243256e555b209d87cd477c5ef3809c..4e9329bc73dea4e2a4fa7040d4102237c6d39660 100644 --- a/v2.5/toolkit/fastnlotoolkit/fastNLOLHAPDF.cc +++ b/v2.5/toolkit/fastnlotoolkit/fastNLOLHAPDF.cc @@ -59,7 +59,7 @@ fastNLOLHAPDF::~fastNLOLHAPDF() { //______________________________________________________________________________ -fastNLOLHAPDF::fastNLOLHAPDF(string name, string LHAPDFFile, int PDFMember, std::string verbosity) : fastNLOReader(name, verbosity), fnPDFs(0) , fiPDFMember(0) , fchksum(0.) { +fastNLOLHAPDF::fastNLOLHAPDF(string name, string LHAPDFFile, int PDFMember) : fastNLOReader(name), fnPDFs(0) , fiPDFMember(0) , fchksum(0.) { PDFSet = NULL; PDF = NULL; SetLHAPDFFilename(LHAPDFFile); @@ -71,7 +71,7 @@ fastNLOLHAPDF::fastNLOLHAPDF(string name, string LHAPDFFile, int PDFMember, std: //______________________________________________________________________________ -fastNLOLHAPDF::fastNLOLHAPDF(const fastNLOTable& table, string LHAPDFFile, int PDFMember, std::string verbosity) : fastNLOReader(table, verbosity), fnPDFs(0) , fiPDFMember(0) , fchksum(0.) { +fastNLOLHAPDF::fastNLOLHAPDF(const fastNLOTable& table, string LHAPDFFile, int PDFMember) : fastNLOReader(table), fnPDFs(0) , fiPDFMember(0) , fchksum(0.) { PDFSet = NULL; PDF = NULL; SetLHAPDFFilename(LHAPDFFile); @@ -346,8 +346,100 @@ vector<double> fastNLOLHAPDF::CalcPDFUncertaintySymm(const vector<LHAPDF::PDFUn } +// +// Evaluation of uncertainties +// +// alpha_s uncertainty //______________________________________________________________________________ -XsUncertainty fastNLOLHAPDF::GetPDFUncertainty(const fastNLO::EPDFUncertaintyStyle ePDFUnc, bool lNorm) { +XsUncertainty fastNLOLHAPDF::GetXsUncertainty(const fastNLO::EAsUncertaintyStyle eAsUnc, bool lNorm) { + // Get a_s(M_Z) uncertainty, da_s(M_Z) as of PDG2016 + const double dasmz[2] = {-0.0011, 0.0011}; + double asmz = GetAlphasMz(); + XsUncertainty XsUnc; + + unsigned int NObsBin = GetNObsBin(); + + logger.info["GetAsUncertainty"]<<"Current a_s(M_Z) = a_s("<<PDG_MZ<<") = "<<asmz<<endl; + logger.info["GetAsUncertainty"]<<"da_s(M_Z) = + "<<dasmz[1]<<" - "<<-dasmz[0]<<endl; + if ( eAsUnc == fastNLO::kAsNone ) { + logger.info["GetAsUncertainty"]<<"Only default value selected, uncertainties will be zero."<<endl; + } else if ( eAsUnc == fastNLO::kAsGRV ) { + logger.info["GetAsUncertainty"]<<"GRV evolution used for a_s(M_Z) uncertainty."<<endl; + } else { + logger.error["GetAsUncertainty"]<<"ERROR! Unknown a_s(M_Z) uncertainty type selected, exiting."<<endl; + logger.error["GetAsUncertainty"]<<"type = "<<eAsUnc<<endl; + exit(1); + } + + vector < double > MyAsMz; + MyAsMz.push_back(asmz); + MyAsMz.push_back(asmz+dasmz[0]); + MyAsMz.push_back(asmz+dasmz[1]); + vector < double > MyXSection; + //! Cross section and absolute uncertainties + for ( unsigned int ias = 0; ias < MyAsMz.size(); ias++ ) { + SetAlphasMz(MyAsMz[ias]); + CalcCrossSection(); + MyXSection = GetCrossSection(lNorm); + for ( unsigned int iobs = 0; iobs < NObsBin; iobs++ ) { + if ( ias == 0 ) { + XsUnc.xs.push_back(MyXSection[iobs]); + XsUnc.dxsu.push_back(0); + XsUnc.dxsl.push_back(0); + } else { + XsUnc.dxsu[iobs] = max(XsUnc.dxsu[iobs],MyXSection[iobs]-XsUnc.xs[iobs]); + XsUnc.dxsl[iobs] = min(XsUnc.dxsl[iobs],MyXSection[iobs]-XsUnc.xs[iobs]); + } + } + } + + //! Divide by cross section != 0 to give relative uncertainties + for ( unsigned int iobs = 0; iobs < NObsBin; iobs++ ) { + if ( fabs(XsUnc.xs[iobs]) > DBL_MIN ) { + XsUnc.dxsu[iobs] = +fabs(XsUnc.dxsu[iobs] / XsUnc.xs[iobs]); + XsUnc.dxsl[iobs] = -fabs(XsUnc.dxsl[iobs] / XsUnc.xs[iobs]); + } else { + XsUnc.dxsu[iobs] = 0.; + XsUnc.dxsl[iobs] = 0.; + } + logger.debug["GetAsUncertainty"]<<"iobs = " << iobs << ", dxsl = " << XsUnc.dxsl[iobs] << ", dxsu = " << XsUnc.dxsu[iobs] <<endl; + } + logger.info["GetAsUncertainty"]<<"Setting a_s(M_Z) back to initial value of "<<asmz<<endl; + SetAlphasMz(asmz); + + return XsUnc; +} + + +//______________________________________________________________________________ + +std::vector< std::vector<double> > fastNLOLHAPDF::GetXsUncertaintyVec(const fastNLO::EAsUncertaintyStyle eAsUnc, bool lNorm, int iprint) { + XsUncertainty xsUnc = fastNLOLHAPDF::GetXsUncertainty(eAsUnc, lNorm); + if (iprint > 0) { + string style{AsUncertaintyStyle_to_string(eAsUnc)}; + string UncName = " # Relative a_s(M_Z) uncertainties (" + style + ")"; + fastNLOTools::PrintXSUncertainty(xsUnc, UncName); + } + std::vector<std::vector<double> > xsUncVec; + xsUncVec.resize(3); + xsUncVec[0] = xsUnc.xs; + xsUncVec[1] = xsUnc.dxsu; + xsUncVec[2] = xsUnc.dxsl; + return xsUncVec; +} + + +//______________________________________________________________________________ + +void fastNLOLHAPDF::PrintXsUncertaintyVec(fastNLO::EAsUncertaintyStyle eAsUnc, std::string UncName, bool lNorm) { + XsUncertainty xsUnc = GetXsUncertainty(eAsUnc, lNorm); + fastNLOTools::PrintXSUncertainty(xsUnc, UncName); +} + + +// PDF uncertainty +//______________________________________________________________________________ +XsUncertainty fastNLOLHAPDF::GetXsUncertainty(const fastNLO::EPDFUncertaintyStyle ePDFUnc, bool lNorm) { XsUncertainty XsUnc; unsigned int NObsBin = GetNObsBin(); unsigned int nMem = GetNPDFMaxMember(); @@ -530,8 +622,8 @@ XsUncertainty fastNLOLHAPDF::GetPDFUncertainty(const fastNLO::EPDFUncertaintySty //______________________________________________________________________________ -std::vector< std::vector<double> > fastNLOLHAPDF::GetPDFUncertaintyVec(const fastNLO::EPDFUncertaintyStyle ePDFUnc, bool lNorm, int iprint) { - XsUncertainty xsUnc = GetPDFUncertainty(ePDFUnc, lNorm); +std::vector< std::vector<double> > fastNLOLHAPDF::GetXsUncertaintyVec(const fastNLO::EPDFUncertaintyStyle ePDFUnc, bool lNorm, int iprint) { + XsUncertainty xsUnc = GetXsUncertainty(ePDFUnc, lNorm); if (iprint > 0) { string style{PDFUncertaintyStyle_to_string(ePDFUnc)}; string PDFFile = GetLHAPDFFilename(); @@ -549,94 +641,37 @@ std::vector< std::vector<double> > fastNLOLHAPDF::GetPDFUncertaintyVec(const fas //______________________________________________________________________________ -void fastNLOLHAPDF::PrintPDFUncertaintyVec(fastNLO::EPDFUncertaintyStyle ePDFUnc, std::string UncName, bool lNorm) { - XsUncertainty xsUnc = GetPDFUncertainty(ePDFUnc, lNorm); +void fastNLOLHAPDF::PrintXsUncertaintyVec(fastNLO::EPDFUncertaintyStyle ePDFUnc, std::string UncName, bool lNorm) { + XsUncertainty xsUnc = GetXsUncertainty(ePDFUnc, lNorm); fastNLOTools::PrintXSUncertainty(xsUnc, UncName); } +// Scale uncertainty (from fastNLOReader) //______________________________________________________________________________ -XsUncertainty fastNLOLHAPDF::GetAsUncertainty(const fastNLO::EAsUncertaintyStyle eAsUnc, bool lNorm) { - // Get a_s(M_Z) uncertainty, da_s(M_Z) as of PDG2016 - const double dasmz[2] = {-0.0011, 0.0011}; - double asmz = GetAlphasMz(); - XsUncertainty XsUnc; - - unsigned int NObsBin = GetNObsBin(); - - logger.info["GetAsUncertainty"]<<"Current a_s(M_Z) = a_s("<<PDG_MZ<<") = "<<asmz<<endl; - logger.info["GetAsUncertainty"]<<"da_s(M_Z) = + "<<dasmz[1]<<" - "<<-dasmz[0]<<endl; - if ( eAsUnc == fastNLO::kAsNone ) { - logger.info["GetAsUncertainty"]<<"Only default value selected, uncertainties will be zero."<<endl; - } else if ( eAsUnc == fastNLO::kAsGRV ) { - logger.info["GetAsUncertainty"]<<"GRV evolution used for a_s(M_Z) uncertainty."<<endl; - } else { - logger.error["GetAsUncertainty"]<<"ERROR! Unknown a_s(M_Z) uncertainty type selected, exiting."<<endl; - logger.error["GetAsUncertainty"]<<"type = "<<eAsUnc<<endl; - exit(1); - } - - vector < double > MyAsMz; - MyAsMz.push_back(asmz); - MyAsMz.push_back(asmz+dasmz[0]); - MyAsMz.push_back(asmz+dasmz[1]); - vector < double > MyXSection; - //! Cross section and absolute uncertainties - for ( unsigned int ias = 0; ias < MyAsMz.size(); ias++ ) { - SetAlphasMz(MyAsMz[ias]); - CalcCrossSection(); - MyXSection = GetCrossSection(lNorm); - for ( unsigned int iobs = 0; iobs < NObsBin; iobs++ ) { - if ( ias == 0 ) { - XsUnc.xs.push_back(MyXSection[iobs]); - XsUnc.dxsu.push_back(0); - XsUnc.dxsl.push_back(0); - } else { - XsUnc.dxsu[iobs] = max(XsUnc.dxsu[iobs],MyXSection[iobs]-XsUnc.xs[iobs]); - XsUnc.dxsl[iobs] = min(XsUnc.dxsl[iobs],MyXSection[iobs]-XsUnc.xs[iobs]); - } - } - } - - //! Divide by cross section != 0 to give relative uncertainties - for ( unsigned int iobs = 0; iobs < NObsBin; iobs++ ) { - if ( fabs(XsUnc.xs[iobs]) > DBL_MIN ) { - XsUnc.dxsu[iobs] = +fabs(XsUnc.dxsu[iobs] / XsUnc.xs[iobs]); - XsUnc.dxsl[iobs] = -fabs(XsUnc.dxsl[iobs] / XsUnc.xs[iobs]); - } else { - XsUnc.dxsu[iobs] = 0.; - XsUnc.dxsl[iobs] = 0.; - } - logger.debug["GetAsUncertainty"]<<"iobs = " << iobs << ", dxsl = " << XsUnc.dxsl[iobs] << ", dxsu = " << XsUnc.dxsu[iobs] <<endl; - } - logger.info["GetAsUncertainty"]<<"Setting a_s(M_Z) back to initial value of "<<asmz<<endl; - SetAlphasMz(asmz); - - return XsUnc; +XsUncertainty fastNLOLHAPDF::GetXsUncertainty(const fastNLO::EScaleUncertaintyStyle eScaleUnc, bool lNorm, double sclfac) { + return fastNLOReader::GetXsUncertainty(eScaleUnc, lNorm, sclfac); } - - //______________________________________________________________________________ - -std::vector< std::vector<double> > fastNLOLHAPDF::GetAsUncertaintyVec(const fastNLO::EAsUncertaintyStyle eAsUnc, bool lNorm, int iprint) { - XsUncertainty xsUnc = fastNLOLHAPDF::GetAsUncertainty(eAsUnc, lNorm); - if (iprint > 0) { - string style{AsUncertaintyStyle_to_string(eAsUnc)}; - string UncName = " # Relative a_s(M_Z) uncertainties (" + style + ")"; - fastNLOTools::PrintXSUncertainty(xsUnc, UncName); - } - std::vector<std::vector<double> > xsUncVec; - xsUncVec.resize(3); - xsUncVec[0] = xsUnc.xs; - xsUncVec[1] = xsUnc.dxsu; - xsUncVec[2] = xsUnc.dxsl; - return xsUncVec; +std::vector< std::vector<double> > fastNLOLHAPDF::GetXsUncertaintyVec(const fastNLO::EScaleUncertaintyStyle eScaleUnc, bool lNorm, int iprint, double sclfac) { + return fastNLOReader::GetXsUncertaintyVec(eScaleUnc, lNorm, iprint, sclfac); +} +//______________________________________________________________________________ +void fastNLOLHAPDF::PrintXsUncertaintyVec(fastNLO::EScaleUncertaintyStyle eScaleUnc, std::string UncName, bool lNorm, double sclfac) { + fastNLOReader::PrintXsUncertaintyVec(eScaleUnc, UncName, lNorm, sclfac); } +// Statistical uncertainty (from fastNLOReader) //______________________________________________________________________________ - -void fastNLOLHAPDF::PrintAsUncertaintyVec(fastNLO::EAsUncertaintyStyle eAsUnc, std::string UncName, bool lNorm) { - XsUncertainty xsUnc = GetAsUncertainty(eAsUnc, lNorm); - fastNLOTools::PrintXSUncertainty(xsUnc, UncName); +XsUncertainty fastNLOLHAPDF::GetXsUncertainty(const fastNLO::EStatUncertaintyStyle eStatUnc, bool lNorm) { + return fastNLOReader::GetXsUncertainty(eStatUnc, lNorm); +} +//______________________________________________________________________________ +std::vector< std::vector<double> > fastNLOLHAPDF::GetXsUncertaintyVec(const fastNLO::EStatUncertaintyStyle eStatUnc, bool lNorm, int iprint) { + return fastNLOReader::GetXsUncertaintyVec(eStatUnc, lNorm, iprint); +} +//______________________________________________________________________________ +void fastNLOLHAPDF::PrintXsUncertaintyVec(fastNLO::EStatUncertaintyStyle eStatUnc, std::string UncName, bool lNorm) { + fastNLOReader::PrintXsUncertaintyVec(eStatUnc, UncName, lNorm); } diff --git a/v2.5/toolkit/fastnlotoolkit/fastNLOReader.cc b/v2.5/toolkit/fastnlotoolkit/fastNLOReader.cc index e08f051c4490ee44d1f9b56da3c99bd80fb07e32..a04cc405c232908c6dd557ffe8111ac1b6c98961 100644 --- a/v2.5/toolkit/fastnlotoolkit/fastNLOReader.cc +++ b/v2.5/toolkit/fastnlotoolkit/fastNLOReader.cc @@ -457,10 +457,8 @@ fastNLOReader::fastNLOReader() : fastNLOTable() { //______________________________________________________________________________ -fastNLOReader::fastNLOReader(string filename, string verbosity) : fastNLOTable(filename, verbosity) { - cout << "AAAAAAAAAAAAAAAAAAA " << GetGlobalVerbosity() << endl; - say::SetGlobalVerbosity(say::toVerbosity()[verbosity]); - cout << "BBBBBBBBBBBBBBBBBBB " << GetGlobalVerbosity() << endl; +fastNLOReader::fastNLOReader(string filename) : fastNLOTable(filename) { + // say::SetGlobalVerbosity(say::toVerbosity()[verbosity]); logger.SetClassName("fastNLOReader"); logger.debug["fastNLOReader"]<<"New fastNLOReader reading filename="<<filename<<endl; fUnits = fastNLO::kPublicationUnits; @@ -475,8 +473,8 @@ fastNLOReader::fastNLOReader(string filename, string verbosity) : fastNLOTable(f //______________________________________________________________________________ -fastNLOReader::fastNLOReader(const fastNLOTable& table, string verbosity) : fastNLOTable(table, verbosity) { - say::SetGlobalVerbosity(say::toVerbosity()[verbosity]); +fastNLOReader::fastNLOReader(const fastNLOTable& table) : fastNLOTable(table) { + // say::SetGlobalVerbosity(say::toVerbosity()[verbosity]); logger.SetClassName("fastNLOReader"); fUnits = fastNLO::kPublicationUnits; fMuRFunc = fastNLO::kScale1; @@ -493,7 +491,7 @@ fastNLOReader::~fastNLOReader(void) { } //______________________________________________________________________________ -fastNLOReader::fastNLOReader(const fastNLOReader& other, string verbosity) : +fastNLOReader::fastNLOReader(const fastNLOReader& other) : fastNLOTable(other), ffilename(other.ffilename), fScalevar(other.fScalevar), fScaleFacMuR(other.fScaleFacMuR), fUnits(other.fUnits), fPDFSuccess(other.fPDFSuccess), fPDFCached(other.fPDFCached), @@ -503,7 +501,7 @@ fastNLOReader::fastNLOReader(const fastNLOReader& other, string verbosity) : XSectionRef_s1(other.XSectionRef_s1), XSectionRef_s2(other.XSectionRef_s2) { //! Copy constructor - say::SetGlobalVerbosity(say::toVerbosity()[verbosity]); + // say::SetGlobalVerbosity(say::toVerbosity()[verbosity]); OrderCoefficients(); // initialize pointers to fCoeff's } @@ -3398,8 +3396,12 @@ double fastNLOReader::RescaleCrossSectionUnits(double binsize, int xunits) { } +// +// Evaluation of uncertainties +// +// Scale uncertainty //______________________________________________________________________________ -XsUncertainty fastNLOReader::GetScaleUncertainty(const EScaleUncertaintyStyle eScaleUnc, bool lNorm, double sclfac) { +XsUncertainty fastNLOReader::GetXsUncertainty(const EScaleUncertaintyStyle eScaleUnc, bool lNorm, double sclfac) { // Get 2-, 6- or, 31-point scale uncertainty around sclfac * central scale (31 only for normalised x sections/ratios) const double xmur0[7] = {1.0, 0.5, 2.0, 0.5, 1.0, 1.0, 2.0}; const double xmuf0[7] = {1.0, 0.5, 2.0, 1.0, 0.5, 2.0, 1.0}; @@ -3479,8 +3481,8 @@ XsUncertainty fastNLOReader::GetScaleUncertainty(const EScaleUncertaintyStyle eS //______________________________________________________________________________ -std::vector< std::vector<double> > fastNLOReader::GetScaleUncertaintyVec(const EScaleUncertaintyStyle eScaleUnc, bool lNorm, int iprint, double sclfac) { - XsUncertainty xsUnc = fastNLOReader::GetScaleUncertainty(eScaleUnc, lNorm, sclfac); +std::vector< std::vector<double> > fastNLOReader::GetXsUncertaintyVec(const EScaleUncertaintyStyle eScaleUnc, bool lNorm, int iprint, double sclfac) { + XsUncertainty xsUnc = fastNLOReader::GetXsUncertainty(eScaleUnc, lNorm, sclfac); if (iprint > 0) { string style{ScaleUncertaintyStyle_to_string(eScaleUnc)}; string UncName = " # Relative scale uncertainties (" + style + ")"; @@ -3497,16 +3499,15 @@ std::vector< std::vector<double> > fastNLOReader::GetScaleUncertaintyVec(const E //______________________________________________________________________________ -void fastNLOReader::PrintScaleUncertaintyVec(fastNLO::EScaleUncertaintyStyle eScaleUnc, std::string UncName, bool lNorm, double sclfac) { - XsUncertainty xsUnc = GetScaleUncertainty(eScaleUnc, lNorm, sclfac); +void fastNLOReader::PrintXsUncertaintyVec(fastNLO::EScaleUncertaintyStyle eScaleUnc, std::string UncName, bool lNorm, double sclfac) { + XsUncertainty xsUnc = GetXsUncertainty(eScaleUnc, lNorm, sclfac); fastNLOTools::PrintXSUncertainty(xsUnc, UncName); } -// Added to include CoeffInfoBlocks -// +// Statistical uncertainty (from CoeffInfoBlocks) //______________________________________________________________________________ -XsUncertainty fastNLOReader::GetAddUncertainty(const EAddUncertaintyStyle eAddUnc, bool lNorm) { +XsUncertainty fastNLOReader::GetXsUncertainty(const EStatUncertaintyStyle eStatUnc, bool lNorm) { // XsUncertainty XsUnc; vector < double > MyXSection; @@ -3519,23 +3520,23 @@ XsUncertainty fastNLOReader::GetAddUncertainty(const EAddUncertaintyStyle eAddUn MydXSection = GetUncertainty(lNorm); //! Fill return struct - if (eAddUnc == kAddNone) { - logger.info["GetAddUncertainty"]<<"No additional uncertainty selected, uncertainties will be zero."<<endl; + if (eStatUnc == kStatNone) { + logger.info["GetStatUncertainty"]<<"No statistical uncertainty selected, uncertainties will be zero."<<endl; for (unsigned int iobs = 0; iobs < NObsBin; iobs++) { XsUnc.xs.push_back(MyXSection[iobs]); XsUnc.dxsu.push_back(0); XsUnc.dxsl.push_back(0); } - } else if (eAddUnc == kAddStat) { - logger.info["GetAddUncertainty"]<<"Statistical/numerical uncertainties selected."<<endl; + } else if (eStatUnc == kStatInt) { + logger.info["GetStatUncertainty"]<<"Statistical integration uncertainty selected."<<endl; for (unsigned int iobs = 0; iobs < NObsBin; iobs++) { XsUnc.xs.push_back(MyXSection[iobs]); XsUnc.dxsu.push_back(MydXSection[iobs]); XsUnc.dxsl.push_back(-MydXSection[iobs]); } } else { - logger.error["GetAddUncertainty"]<<"ERROR! No valid additional uncertainty style selected, exiting."<<endl; - logger.error["GetAddUncertainty"]<<"Style enum = "<<eAddUnc<<endl; + logger.error["GetStatUncertainty"]<<"ERROR! No valid statistical uncertainty style selected, exiting."<<endl; + logger.error["GetStatUncertainty"]<<"Style enum = "<<eStatUnc<<endl; exit(1); } @@ -3548,7 +3549,7 @@ XsUncertainty fastNLOReader::GetAddUncertainty(const EAddUncertaintyStyle eAddUn XsUnc.dxsu[iobs] = 0.; XsUnc.dxsl[iobs] = 0.; } - logger.debug["GetAddUncertainty"]<<"iobs = " << iobs << ", dxsl = " << XsUnc.dxsl[iobs] << ", dxsu = " << XsUnc.dxsu[iobs] <<endl; + logger.debug["GetStatUncertainty"]<<"iobs = " << iobs << ", dxsl = " << XsUnc.dxsl[iobs] << ", dxsu = " << XsUnc.dxsu[iobs] <<endl; } return XsUnc; @@ -3556,10 +3557,10 @@ XsUncertainty fastNLOReader::GetAddUncertainty(const EAddUncertaintyStyle eAddUn //______________________________________________________________________________ -std::vector< std::vector<double> > fastNLOReader::GetAddUncertaintyVec(const EAddUncertaintyStyle eAddUnc, bool lNorm, int iprint) { - XsUncertainty xsUnc = fastNLOReader::GetAddUncertainty(eAddUnc, lNorm); +std::vector< std::vector<double> > fastNLOReader::GetXsUncertaintyVec(const EStatUncertaintyStyle eStatUnc, bool lNorm, int iprint) { + XsUncertainty xsUnc = fastNLOReader::GetXsUncertainty(eStatUnc, lNorm); if (iprint > 0) { - string style{AddUncertaintyStyle_to_string(eAddUnc)}; + string style{StatUncertaintyStyle_to_string(eStatUnc)}; string UncName = " # Relative statistical uncertainties (" + style + ")"; fastNLOTools::PrintXSUncertainty(xsUnc, UncName); } @@ -3574,7 +3575,7 @@ std::vector< std::vector<double> > fastNLOReader::GetAddUncertaintyVec(const EAd //______________________________________________________________________________ -void fastNLOReader::PrintAddUncertaintyVec(fastNLO::EAddUncertaintyStyle eAddUnc, std::string UncName, bool lNorm) { - XsUncertainty xsUnc = GetAddUncertainty(eAddUnc, lNorm); +void fastNLOReader::PrintXsUncertaintyVec(fastNLO::EStatUncertaintyStyle eStatUnc, std::string UncName, bool lNorm) { + XsUncertainty xsUnc = GetXsUncertainty(eStatUnc, lNorm); fastNLOTools::PrintXSUncertainty(xsUnc, UncName); } diff --git a/v2.5/toolkit/fastnlotoolkit/fastNLOTable.cc b/v2.5/toolkit/fastnlotoolkit/fastNLOTable.cc index 03e5ec921fba941d9a0f07be4e65e259d7a513a9..c8850c94008fdeb21553aae821e54e5abaa2f786 100644 --- a/v2.5/toolkit/fastnlotoolkit/fastNLOTable.cc +++ b/v2.5/toolkit/fastnlotoolkit/fastNLOTable.cc @@ -36,10 +36,8 @@ fastNLOTable::fastNLOTable() // ___________________________________________________________________________________________________ -fastNLOTable::fastNLOTable(string name, string verbosity) : ffilename(name), fPrecision(8), logger("fastNLOTable") { +fastNLOTable::fastNLOTable(string name) : ffilename(name), fPrecision(8), logger("fastNLOTable") { //logger.SetClassName("fastNLOTable"); - say::SetGlobalVerbosity(say::toVerbosity()[verbosity]); - cout << "DDDDDDDDDDDDDDDDDDDDDDDD " << say::GetGlobalVerbosity() << endl; if (!fWelcomeOnce) PrintWelcomeMessage(); ReadTable(); } @@ -52,7 +50,7 @@ fastNLOTable::~fastNLOTable(){ } // ___________________________________________________________________________________________________ -fastNLOTable::fastNLOTable(const fastNLOTable& other, string verbosity) +fastNLOTable::fastNLOTable(const fastNLOTable& other) : ffilename(other.ffilename), fPrecision(other.fPrecision), ITabVersionRead(other.ITabVersionRead), ITabVersionWrite(other.ITabVersionRead), @@ -67,7 +65,7 @@ fastNLOTable::fastNLOTable(const fastNLOTable& other, string verbosity) IDivUpPointer(other.IDivUpPointer) { //! Copy constructor - say::SetGlobalVerbosity(say::toVerbosity()[verbosity]); + // say::SetGlobalVerbosity(say::toVerbosity()[verbosity]); logger.SetClassName("fastNLOTable"); for (size_t i = 0; i < other.fCoeff.size(); ++i) { fCoeff[i] = other.fCoeff[i]->Clone(); @@ -1958,7 +1956,7 @@ void fastNLOTable::PrintContributionSummary(int iprint) const { // logger.debug["PrintContributionSummary"] << "Printing flag iprint = " << iprint << endl; char buffer[1024]; - cout << endl; + // cout << endl; cout << fastNLO::_CSEPSC << endl; logger.shout << "Overview on contribution types and numbers contained in table: " << ffilename << endl; cout << fastNLO::_SSEPSC << endl; @@ -2596,7 +2594,6 @@ void fastNLOTable::PrintWelcomeMessage() { char years[100] = FNLO_YEARS; speaker &infosep = logger.infosep; - infosep << "" << endl; infosep << fastNLO::_CSEPSC << endl; infosep << " #" << endl; infosep << " # " << fnlo << "_" << subproject << endl; diff --git a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOConstants.h.in b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOConstants.h.in index ba7623cbef12ec173570bab9eb98865095ee2ce6..5a4ddb7b8879c1a1e14b9fbe099e213f3bbca0c3 100644 --- a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOConstants.h.in +++ b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOConstants.h.in @@ -105,14 +105,6 @@ namespace fastNLO { kPublicationUnits = 1 // calculate the cross section in units as given in the according publication }; - enum EUncertaintyType { // TODO: KR to be further developed - kNoUncertainty = 0, // no uncertainty, only central cross section evaluated - kScaleUncertainty = 1, // scale uncertainty (as proxy for missing higher orders) - kPDFUncertainty = 2, // PDF uncertainty - kNumUncertainty = 3, // Statistical numerical uncertainty of integration - kNPUncertainty = 4 // Nonperturbative uncertainty - }; - enum EScaleUncertaintyStyle { kScaleNone = 0, // no scale uncertainty, only central scale (mu_r,mu_f) = (1,1) evaluated kSymmetricTwoPoint = 1, // symmetric (mu_r,mu_f) scale variations by factors (1/2,1/2), (2,2) @@ -136,23 +128,23 @@ namespace fastNLO { return{kScaleNone}; } - enum EAddUncertaintyStyle { - kAddNone = 0, // no additional uncertainty - kAddStat = 1 // statistical/numerical uncertainty + enum EStatUncertaintyStyle { + kStatNone = 0, // no statistical uncertainty + kStatInt = 1 // statistical uncertainty from MC integration }; - constexpr std::string_view AddUncertaintyStyle_to_string(EAddUncertaintyStyle estyle) { + constexpr std::string_view StatUncertaintyStyle_to_string(EStatUncertaintyStyle estyle) { switch (estyle) { - case kAddNone: return "NN"; - case kAddStat: return "ST"; - default: return "??"; + case kStatNone: return "NN"; + case kStatInt: return "ST"; + default: return "??"; } } - constexpr std::optional<EAddUncertaintyStyle> AddUncertaintyStyle_to_enum(std::string_view sstyle) { - if (sstyle == "NN") return kAddNone; - if (sstyle == "ST") return kAddStat; - return{kAddNone}; + constexpr std::optional<EStatUncertaintyStyle> StatUncertaintyStyle_to_enum(std::string_view sstyle) { + if (sstyle == "NN") return kStatNone; + if (sstyle == "ST") return kStatInt; + return{kStatNone}; } enum EPDFUncertaintyStyle { diff --git a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOLHAPDF.h b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOLHAPDF.h index 9f9a33b83b755640d054ab23f0ddc53a925d26bb..b6adee55fa7be149d5da4fdfb37989dc5bdfd1a3 100644 --- a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOLHAPDF.h +++ b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOLHAPDF.h @@ -32,8 +32,8 @@ public: fastNLOLHAPDF(std::string name); fastNLOLHAPDF(const fastNLOTable&); ~fastNLOLHAPDF(); - fastNLOLHAPDF(std::string name, std::string LHAPDFfile, int PDFSet = 0, std::string verbosity = "INFO"); - fastNLOLHAPDF(const fastNLOTable&, std::string LHAPDFfile, int PDFSet = 0, std::string verbosity = "INFO"); + fastNLOLHAPDF(std::string name, std::string LHAPDFfile, int PDFSet = 0); + fastNLOLHAPDF(const fastNLOTable&, std::string LHAPDFfile, int PDFSet = 0); // Initializer. Necessary for some alternative evolutions. virtual void InitEvolveAlphas(); @@ -60,14 +60,24 @@ public: double GetAlphasMz() const; - //! Return struct with vectors/vectors containing the cross section values and the selected uncertainty - XsUncertainty GetAsUncertainty( const fastNLO::EAsUncertaintyStyle eAsUnc, bool lNorm = false); - std::vector< std::vector<double> > GetAsUncertaintyVec( const fastNLO::EAsUncertaintyStyle eAsUnc, bool lNorm = false, int iprint = 0); - void PrintAsUncertaintyVec(fastNLO::EAsUncertaintyStyle, std::string UncName, bool lNorm = false); + //! Return struct with vectors (for C++) or vector of vectors (for Python) containing the cross section values and the selected uncertainty + //! Enum of Uncertaintstyle decides on method to call + // Use implementations in fastNLOReader for these + XsUncertainty GetXsUncertainty(const fastNLO::EStatUncertaintyStyle eStatUnc, bool lNorm = false); + std::vector< std::vector<double> > GetXsUncertaintyVec(const fastNLO::EStatUncertaintyStyle eStatUnc, bool lNorm = false, int iprint = 0); + void PrintXsUncertaintyVec(fastNLO::EStatUncertaintyStyle, std::string UncName, bool lNorm = false); // - XsUncertainty GetPDFUncertainty( const fastNLO::EPDFUncertaintyStyle ePDFUnc, bool lNorm = false); - std::vector<std::vector<double> > GetPDFUncertaintyVec(fastNLO::EPDFUncertaintyStyle, bool lNorm = false, int iprint = 0); - void PrintPDFUncertaintyVec(fastNLO::EPDFUncertaintyStyle, std::string UncName, bool lNorm = false); + XsUncertainty GetXsUncertainty(const fastNLO::EScaleUncertaintyStyle eScaleUnc, bool lNorm = false, double sclfac = 1.); + std::vector< std::vector<double> > GetXsUncertaintyVec(const fastNLO::EScaleUncertaintyStyle eScaleUnc, bool lNorm = false, int iprint = 0, double sclfac = 1.); + void PrintXsUncertaintyVec(fastNLO::EScaleUncertaintyStyle, std::string UncName, bool lNorm = false, double sclfac =1.); + // Specific implementations in fastNLOLHAPDF + XsUncertainty GetXsUncertainty(const fastNLO::EAsUncertaintyStyle eAsUnc, bool lNorm = false); + std::vector< std::vector<double> > GetXsUncertaintyVec(const fastNLO::EAsUncertaintyStyle eAsUnc, bool lNorm = false, int iprint = 0); + void PrintXsUncertaintyVec(fastNLO::EAsUncertaintyStyle, std::string UncName, bool lNorm = false); + // + XsUncertainty GetXsUncertainty(const fastNLO::EPDFUncertaintyStyle ePDFUnc, bool lNorm = false); + std::vector<std::vector<double> > GetXsUncertaintyVec(const fastNLO::EPDFUncertaintyStyle, bool lNorm = false, int iprint = 0); + void PrintXsUncertaintyVec(fastNLO::EPDFUncertaintyStyle, std::string UncName, bool lNorm = false); std::vector<LHAPDF::PDFUncertainty> GetPDFUncertaintyLHAPDF(double cl=100*erf(1/sqrt(2)), bool alternative=false); //!< return PDF uncertainty, formulae taken from LHAPDF6 std::vector<double> CalcPDFUncertaintyMinus(const std::vector<LHAPDF::PDFUncertainty>& ) const; //!<get vector<double> for PDF-minus uncertainty. Uncertainties are POSITIVE! diff --git a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOReader.h b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOReader.h index 3f57d4a3db5d12d21c3a46f59773940685d3a68c..d76909a02350dd5be45f72521813f58cb0490a8c 100644 --- a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOReader.h +++ b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOReader.h @@ -24,9 +24,12 @@ class fastNLOReader : public fastNLOTable , public fastNLOPDFLinearCombinations public: typedef double(*mu_func)(double,double); - fastNLOReader(std::string filename, std::string verbosity = "INFO"); - fastNLOReader(const fastNLOTable&, std::string verbosity = "INFO"); - fastNLOReader(const fastNLOReader&, std::string verbosity = "INFO"); + fastNLOReader(std::string filename); + fastNLOReader(const fastNLOTable&); + fastNLOReader(const fastNLOReader&); + // fastNLOReader(std::string filename, std::string verbosity = "INFO"); + // fastNLOReader(const fastNLOTable&, std::string verbosity = "INFO"); + // fastNLOReader(const fastNLOReader&, std::string verbosity = "INFO"); virtual ~fastNLOReader(); void SetFilename(std::string filename) ; void InitScalevariation(); @@ -94,14 +97,15 @@ public: std::vector < double > GetQScales(); //!< Order (power of alpha_s) rel. to LO: 0 --> LO, 1 --> NLO std::vector < std::vector < double > > GetCrossSection2Dim(); - //! Return struct with vectors/vectors containing the cross section values and the selected uncertainty - XsUncertainty GetAddUncertainty( const fastNLO::EAddUncertaintyStyle eAddUnc, bool lNorm = false); - std::vector< std::vector<double> > GetAddUncertaintyVec( const fastNLO::EAddUncertaintyStyle eAddUnc, bool lNorm = false, int iprint = 0); - void PrintAddUncertaintyVec(fastNLO::EAddUncertaintyStyle, std::string UncName, bool lNorm = false); + //! Return struct with vectors (for C++) or vector of vectors (for Python) containing the cross section values and the selected uncertainty + //! Enum of Uncertaintstyle decides on method to call + XsUncertainty GetXsUncertainty(const fastNLO::EScaleUncertaintyStyle eScaleUnc, bool lNorm = false, double sclfac = 1.); + std::vector< std::vector<double> > GetXsUncertaintyVec(const fastNLO::EScaleUncertaintyStyle eScaleUnc, bool lNorm = false, int iprint = 0, double sclfac = 1.); + void PrintXsUncertaintyVec(fastNLO::EScaleUncertaintyStyle, std::string UncName, bool lNorm = false, double sclfac =1.); // - XsUncertainty GetScaleUncertainty( const fastNLO::EScaleUncertaintyStyle eScaleUnc, bool lNorm = false, double sclfac = 1.); - std::vector< std::vector<double> > GetScaleUncertaintyVec( const fastNLO::EScaleUncertaintyStyle eScaleUnc, bool lNorm = false, int iprint = 0, double sclfac = 1.); - void PrintScaleUncertaintyVec(fastNLO::EScaleUncertaintyStyle, std::string UncName, bool lNorm = false, double sclfac =1.); + XsUncertainty GetXsUncertainty(const fastNLO::EStatUncertaintyStyle eStatUnc, bool lNorm = false); + std::vector< std::vector<double> > GetXsUncertaintyVec(const fastNLO::EStatUncertaintyStyle eStatUnc, bool lNorm = false, int iprint = 0); + void PrintXsUncertaintyVec(fastNLO::EStatUncertaintyStyle, std::string UncName, bool lNorm = false); // ---- Getters for fastNLOReader member variables ---- // fastNLO::EScaleFunctionalForm GetMuRFunctionalForm() const { return fMuRFunc; }; diff --git a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOTable.h b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOTable.h index 590fff5c74bb29eb276ed008f1cd46b0858b1c99..de0638f08aefcede200f33ddd0a68cc7f281c523 100644 --- a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOTable.h +++ b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOTable.h @@ -21,9 +21,11 @@ class fastNLOTable { public: fastNLOTable(); - fastNLOTable(std::string filename, std::string verbosity = "INFO"); + fastNLOTable(std::string filename); + // fastNLOTable(std::string filename, std::string verbosity = "INFO"); virtual ~fastNLOTable(); - fastNLOTable(const fastNLOTable&, std::string verbosity = "INFO"); + // fastNLOTable(const fastNLOTable&, std::string verbosity = "INFO"); + fastNLOTable(const fastNLOTable&); virtual void ReadTable(); virtual void WriteTable(); diff --git a/v2.5/toolkit/src/fnlo-tk-cppread.cc b/v2.5/toolkit/src/fnlo-tk-cppread.cc index 1207e263e91bec9757feb969eb572252f0a560ff..f30861dab0fb461048fb66d1c59b808c883e75ff 100644 --- a/v2.5/toolkit/src/fnlo-tk-cppread.cc +++ b/v2.5/toolkit/src/fnlo-tk-cppread.cc @@ -1,8 +1,8 @@ ///******************************************************************** /// /// fastNLO_toolkit: fnlo-tk-cppread -/// Program to read fastNLO tables and derive -/// QCD cross sections using PDFs e.g. from LHAPDF +/// Program to read fastNLO tables and derive QCD cross sections +/// using PDFs e.g. from LHAPDF /// /// For more explanations type: /// ./fnlo-tk-cppread -h @@ -23,7 +23,6 @@ #include <iomanip> #include <iostream> #include <string> -#include <string> #include <unordered_map> #include <utility> #include <vector> @@ -55,14 +54,27 @@ int main(int argc, char** argv) { //! namespaces using namespace std; - using namespace say; //! namespace for 'speaker.h'-verbosity levels - using namespace fastNLO; //! namespace for fastNLO constants + using namespace say; //! namespace for 'speaker.h'-verbosity levels + using namespace fastNLO; //! namespace for fastNLO constants //! --- Set initial verbosity level - SetGlobalVerbosity(INFO); + // Could set verbosity level directly using enum here, but we'll need fVerbosity anyway for cmd line parsing + // SetGlobalVerbosity(WARNING); + string VerbosityLevel = "WARNING"; + say::Verbosity fVerbosity = toVerbosity()[VerbosityLevel]; + SetGlobalVerbosity(fVerbosity); + + //--- Set verbosity level of table evaluation from command line argument if existing + if (argc > 13) { + VerbosityLevel = (const char*) argv[13]; + //! --- Reset verbosity level from here on + fVerbosity = toVerbosity()[VerbosityLevel]; + SetGlobalVerbosity(fVerbosity); + } //! --- Parse command line - // TODO: Use getopt or getopt_long standard! + // TODO: Use getopt_long or argparse, e.g. https://github.com/p-ranav/argparse + // or https://github.com/morrisfranken/argparse ? char buffer[1024]; // TODO: Use PATH_MAX instead? string tablename; if (argc <= 1) { @@ -71,7 +83,7 @@ int main(int argc, char** argv) { shout["fnlo-tk-cppread"] << "fastNLO Cross-Section Calculator"<<endl; yell << _SSEPSC << endl; error["fnlo-tk-cppread"] << "No fastNLO table specified!" << endl; - shout["fnlo-tk-cppread"] << "For more explanations type:" << endl; + shout["fnlo-tk-cppread"] << "For an explanation of command line arguments type:" << endl; shout["fnlo-tk-cppread"] << "./fnlo-tk-cppread -h" << endl; shout["fnlo-tk-cppread"] << "For version number printout type:" << endl; shout["fnlo-tk-cppread"] << "./fnlo-tk-cppread -v" << endl; @@ -82,24 +94,27 @@ int main(int argc, char** argv) { if (tablename == "-v") { fastNLOTools::PrintFastnloVersion(); return 0; + } else if (tablename == "-h") { + // Avoid suppressing usage printouts + SetGlobalVerbosity(INFO); } //! --- Print program purpose yell << _CSEPSC << endl; info["fnlo-tk-cppread"] << "Program to read fastNLO tables and derive" << endl; info["fnlo-tk-cppread"] << "QCD cross sections using PDFs e.g. from LHAPDF" << endl; - yell << _SSEPSC << endl; - info["fnlo-tk-cppread"] << "For more explanations type:" << endl; + infosep << _SSEPSC << endl; + info["fnlo-tk-cppread"] << "For an explanation of command line arguments type:" << endl; info["fnlo-tk-cppread"] << "./fnlo-tk-cppread -h" << endl; info["fnlo-tk-cppread"] << "For version number printout type:" << endl; info["fnlo-tk-cppread"] << "./fnlo-tk-cppread -v" << endl; - yell << _CSEPSC << endl; - yell << "" << endl; + infosep << _CSEPSC << endl; + infosep << "" << endl; + infosep << _CSEPSC << endl; //! --- Usage info if (tablename == "-h") { - yell << _CSEPSC << endl; info["fnlo-tk-cppread"] << "fastNLO Cross-Section Calculator"<<endl; - yell << _SSEPSC << endl; - yell << " #" << endl; + infosep << _SSEPSC << endl; + infosep << " #" << endl; info["fnlo-tk-cppread"] << "This program evaluates a fastNLO table for a set of specified options and" << endl; info["fnlo-tk-cppread"] << "prints out a table with detailed binning and cross-section information" << endl; info["fnlo-tk-cppread"] << "for each observable bin." << endl; @@ -145,22 +160,18 @@ int main(int argc, char** argv) { man << "[newEcms]: Reweight grid from default Ecms to new value newEcms, def. = no reweighting" << endl; man << " ATTENTION! Missing phase space when increasing Ecms; still loss in precision when decreasing Ecms!" << endl; man << "[Verbosity]: Set verbosity level of table evaluation [DEBUG,INFO,WARNING,ERROR], def. = WARNING" << endl; - yell << " #" << endl; + infosep << " #" << endl; man << "Use \"_\" to skip changing a default argument." << endl; - yell << " #" << endl; - yell << _CSEPSC << endl; + infosep << " #" << endl; + infosep << _CSEPSC << endl; return 0; } else { - yell << _CSEPSC << endl; - shout["fnlo-tk-cppread"] << "fastNLO Cross-Section Calculator"<<endl; - yell << _SSEPSC << endl; shout["fnlo-tk-cppread"] << "Evaluating table: " << tablename << endl; } } - //--- PDF set - string chtmp = "X"; - string PDFFile = "X"; + //--- PDF choice + string PDFFile; if (argc > 2) { PDFFile = (const char*) argv[2]; } @@ -183,11 +194,12 @@ int main(int argc, char** argv) { const double fixmuflow[] = { 1.0, 2.718281828459045, 4.481689070338065, 4.481689070338065, 2.718281828459045, 12.18249396070347, 2.718281828459045 }; const double fixmurhigh[] = { 1.0, 90.0171313005, 54.5981500331, 148.4131591026, 90.0171313005, 54.5981500331, 90.0171313005 }; const double fixmufhigh[] = { 1.0, 90.0171313005, 54.5981500331, 148.4131591026, 54.5981500331, 90.0171313005, 148.4131591026 }; + string chtmp; if (argc > 3) { chtmp = (const char*) argv[3]; } if (argc <= 3 || chtmp == "_") { - shout["fnlo-tk-cppread"] << "No request given for PDF or scale variations, using primary scale only." << endl; + shout["fnlo-tk-cppread"] << "No request for PDF or scale variations, using primary scale only." << endl; } else { nvars = atoi(argv[3]); if (nvars < -nfixmax) { @@ -224,19 +236,19 @@ int main(int argc, char** argv) { } //--- alpha_s evolution code - string AsEvolCode = "GRV"; + string AsEvolCode; if (argc > 4) { AsEvolCode = (const char*) argv[4]; } if (argc <= 4 || AsEvolCode == "_") { AsEvolCode = "GRV"; - shout["fnlo-tk-cppread"] << "No request given for alpha_s evolution code, using GRV default." << endl; + shout["fnlo-tk-cppread"] << "No request for alpha_s evolution code, using GRV default." << endl; } else { shout["fnlo-tk-cppread"] << "Using alpha_s evolution code: " << AsEvolCode << endl; } //--- Normalization - string chnorm = "no"; + string chnorm; if (argc > 5) { chnorm = (const char*) argv[5]; } @@ -247,7 +259,7 @@ int main(int argc, char** argv) { } //--- Scale choice (flex-scale tables only; ignored for fix-scale tables) - string chflex = "kScale1"; + string chflex; if (argc > 6) { chflex = (const char*) argv[6]; } @@ -264,8 +276,7 @@ int main(int argc, char** argv) { chtmp = (const char*) argv[7]; } if (argc <= 7 || chtmp == "_") { - shout["fnlo-tk-cppread"] << "No request given for central scale factor," << endl; - shout << " using unmodified central scale." << endl; + shout["fnlo-tk-cppread"] << "No request for central scale factor, using default of unity." << endl; } else { sclfac = atof(argv[7]); if ( sclfac < 0.1 || sclfac > 10. ) { @@ -289,8 +300,7 @@ int main(int argc, char** argv) { chtmp = (const char*) argv[8]; } if (argc <= 8 || chtmp == "_") { - shout["fnlo-tk-cppread"] << "No request given for no. of flavours," << endl; - shout << " using default value for LHAPDF or Nf = " << Nf << "." << endl; + shout["fnlo-tk-cppread"] << "No request for no. of flavours, using default of LHAPDF or Nf = " << Nf << "." << endl; } else { Nf = atoi(argv[8]); if ( AsEvolCode == "LHAPDF" ) { @@ -312,8 +322,7 @@ int main(int argc, char** argv) { chtmp = (const char*) argv[9]; } if (argc <= 9 || chtmp == "_") { - shout["fnlo-tk-cppread"] << "No request given for no. of loops," << endl; - shout << " using default value for LHAPDF or NLoop = " << NLoop << "." << endl; + shout["fnlo-tk-cppread"] << "No request for no. of loops, using default of LHAPDF or NLoop = " << NLoop << "." << endl; } else { NLoop = atoi(argv[9]); if ( AsEvolCode == "LHAPDF" ) { @@ -335,8 +344,7 @@ int main(int argc, char** argv) { chtmp = (const char*) argv[10]; } if (argc <= 10 || chtmp == "_") { - shout["fnlo-tk-cppread"] << "No request given for value of alpha_s(M_Z)," << endl; - shout << " using default value for LHAPDF or asMz = " << asMz << "." << endl; + shout["fnlo-tk-cppread"] << "No request for value of alpha_s(M_Z), using default of LHAPDF or asMz = " << asMz << "." << endl; } else { asMz = atof(argv[10]); if ( AsEvolCode == "LHAPDF" ) { @@ -358,8 +366,7 @@ int main(int argc, char** argv) { chtmp = (const char*) argv[11]; } if (argc <= 11 || chtmp == "_") { - shout["fnlo-tk-cppread"] << "No request given for value of M_Z," << endl; - shout << " using default value for LHAPDF or Mz = " << Mz << "." << endl; + shout["fnlo-tk-cppread"] << "No request for value of M_Z, using default of LHAPDF or Mz = " << Mz << "." << endl; } else { Mz = atof(argv[11]); if ( AsEvolCode == "LHAPDF" ) { @@ -381,8 +388,7 @@ int main(int argc, char** argv) { chtmp = (const char*) argv[12]; } if (argc <= 12 || chtmp == "_") { - shout["fnlo-tk-cppread"] << "No request given for Ecms," << endl; - shout << " will use default value as defined in grid." << endl; + shout["fnlo-tk-cppread"] << "No request for Ecms, using default as defined in grid." << endl; } else { newEcms = atof(argv[12]); if ( newEcms < 10. || newEcms > 100000. ) { @@ -394,15 +400,8 @@ int main(int argc, char** argv) { } } - //--- Set verbosity level of table evaluation - string VerbosityLevel = "WARNING"; - if (argc > 13) { - VerbosityLevel = (const char*) argv[13]; - } - if (argc <= 13 || VerbosityLevel == "_") { - VerbosityLevel = "WARNING"; - shout["fnlo-tk-cppread"] << "No request given for verbosity level, using WARNING default." << endl; - } else { + //--- Report about changed verbosity level at the beginning + if (!VerbosityLevel.empty()) { shout["fnlo-tk-cppread"] << "Using verbosity level: " << VerbosityLevel << endl; } @@ -411,12 +410,9 @@ int main(int argc, char** argv) { error["fnlo-tk-cppread"] << "Too many arguments, aborting!" << endl; exit(1); } - yell << _CSEPSC << endl; + yell << _CSEPSC << endl << endl; //--- End of parsing arguments - //! --- Reset verbosity level from here on - SetGlobalVerbosity(toVerbosity()[VerbosityLevel]); - // ************************** fastNLO and example documentation starts here **************************** // --- fastNLO user: Hello! // If you use fastNLO for the first time, please read through the @@ -849,9 +845,12 @@ int main(int argc, char** argv) { //! ONLY if compiled --with-qcdnum support! fnlo = new fastNLOQCDNUMAS(tablename,PDFFile,0); #else - printf("fnlo-tk-cppread: ERROR! The alpha_s evolution code %s was selected!\n",AsEvolCode.c_str()); - printf(" But the fastNLO Toolkit was compiled without the optional support for this!\n"); - printf(" Please choose another alpha_s evolution code or recompile with %s support.\n",AsEvolCode.c_str()); + snprintf(buffer, sizeof(buffer), "ERROR! The alpha_s evolution code %s was selected!\n",AsEvolCode.c_str()); + error["fnlo-tk-cppread"] << buffer << endl; + snprintf(buffer, sizeof(buffer), "But the fastNLO Toolkit was compiled without the optional support for this!\n"); + error["fnlo-tk-cppread"] << buffer << endl; + snprintf(buffer, sizeof(buffer), "Please choose another alpha_s evolution code or recompile with %s support.\n",AsEvolCode.c_str()); + error["fnlo-tk-cppread"] << buffer << endl; exit(1); #endif } else if (AsEvolCode == "HOPPET") { @@ -859,22 +858,30 @@ int main(int argc, char** argv) { //! ONLY if compiled --with-hoppet support! fnlo = new fastNLOHoppetAs(tablename,PDFFile,0); #else - printf("fnlo-tk-cppread: ERROR! The alpha_s evolution code %s was selected!\n",AsEvolCode.c_str()); - printf(" But the fastNLO Toolkit was compiled without the optional support for this!\n"); - printf(" Please choose another alpha_s evolution code or recompile with %s support.\n",AsEvolCode.c_str()); + snprintf(buffer, sizeof(buffer), "ERROR! The alpha_s evolution code %s was selected!\n",AsEvolCode.c_str()); + error["fnlo-tk-cppread"] << buffer << endl; + snprintf(buffer, sizeof(buffer), "But the fastNLO Toolkit was compiled without the optional support for this!\n"); + error["fnlo-tk-cppread"] << buffer << endl; + snprintf(buffer, sizeof(buffer), "Please choose another alpha_s evolution code or recompile with %s support.\n",AsEvolCode.c_str()); + error["fnlo-tk-cppread"] << buffer << endl; exit(1); #endif } else { - printf("fnlo-tk-cppread: ERROR! Unknown alpha_s evolution code %s!\n",AsEvolCode.c_str()); - printf(" If you compiled with optional QCDNUM or HOPPET support, please\n"); - printf(" do not forget to comment in the marked lines in main.cc!\n"); + snprintf(buffer, sizeof(buffer), "fnlo-tk-cppread: ERROR! Unknown alpha_s evolution code %s!\n",AsEvolCode.c_str()); + error["fnlo-tk-cppread"] << buffer << endl; + snprintf(buffer, sizeof(buffer), " If you compiled with optional QCDNUM or HOPPET support, please\n"); + error["fnlo-tk-cppread"] << buffer << endl; + snprintf(buffer, sizeof(buffer), " do not forget to comment in the marked lines in main.cc!\n"); + error["fnlo-tk-cppread"] << buffer << endl; exit(1); } - //! Print some fastNLO table info + //! Print essential table information // TODO: Add print out of scale info, in particular for flex-scale tables - fnlo->PrintContributionSummary(0); - fnlo->Print(0); + if (fVerbosity < toVerbosity()["WARNING"]) { + fnlo->PrintContributionSummary(0); + fnlo->Print(0); + } //! Define the PDF set and member fnlo->SetLHAPDFFilename(PDFFile); @@ -1389,25 +1396,25 @@ int main(int argc, char** argv) { } //! Start print out - yell << _DSEPLC << endl; - shout << "My Cross Sections" << endl; + cout << _DSEPLC << endl; + cout << " # My Cross Sections" << endl; //! PDF members if ( !sclvar ) { - snprintf(buffer, sizeof(buffer), "The PDF member chosen here is: %i",ivar); + snprintf(buffer, sizeof(buffer), " # The PDF member chosen here is: %i",ivar); } //! Scale factor variations else if ( ivar == 0 || !sclfix ) { - snprintf(buffer, sizeof(buffer), "The scale factors xmur, xmuf chosen here are: % #10.3f, % #10.3f",fnlo->GetScaleFactorMuR(),fnlo->GetScaleFactorMuF()); + snprintf(buffer, sizeof(buffer), " # The scale factors xmur, xmuf chosen here are: % #10.3f, % #10.3f",fnlo->GetScaleFactorMuR(),fnlo->GetScaleFactorMuF()); } //! Fixed-scale variations else if ( sclfix) { - snprintf(buffer, sizeof(buffer), "The fixed scales mur, muf chosen here are: % #10.3f, % #10.3f",fixmur[ivar],fixmuf[ivar]); + snprintf(buffer, sizeof(buffer), " # The fixed scales mur, muf chosen here are: % #10.3f, % #10.3f",fixmur[ivar],fixmuf[ivar]); } //! Undefined else { } - shout << buffer << endl; - yell << _SSEPLC << endl; + cout << buffer << endl; + cout << _SSEPLC << endl; //! Get table constants relevant for print out const int NDim = fnlo->GetNumDiffBin(); @@ -1502,8 +1509,8 @@ int main(int argc, char** argv) { if (NDim == 1) { snprintf(buffer, sizeof(buffer), "%s%s [ %-17.17s ] < %-10.10s > %s", header0.c_str(),headdim0.c_str(),DimLabel[0].c_str(),headscl.c_str(),header2.c_str()); - yell << buffer << endl; - yell << _SSEPLC << endl; + cout << buffer << endl; + cout << _SSEPLC << endl; NDimBins[0] = 0; for (int i=0; i<NObsBin; i++) { @@ -1557,8 +1564,8 @@ int main(int argc, char** argv) { } else if (NDim == 2) { snprintf(buffer, sizeof(buffer), "%s%s [ %-17.17s ] %s [ %-17.17s ] < %-10.10s > %s", header0.c_str(),headdim0.c_str(),DimLabel[0].c_str(),headdim2.c_str(),DimLabel[1].c_str(),headscl.c_str(),header2.c_str()); - yell << buffer << endl; - yell << _SSEPLC << endl; + cout << buffer << endl; + cout << _SSEPLC << endl; for (int i=0; i<NObsBin; i++) { for (int j=0; j<NDim; j++) { if (i==0) { @@ -1615,8 +1622,8 @@ int main(int argc, char** argv) { snprintf(buffer, sizeof(buffer), "%s%s [ %-17.17s ] %s [ %-17.17s ] %s [ %-17.17s ] < %-10.10s > %s", header0.c_str(),headdim0.c_str(),DimLabel[0].c_str(),headdim1.c_str(),DimLabel[1].c_str(), headdim2.c_str(),DimLabel[2].c_str(),headscl.c_str(),header2.c_str()); - yell << buffer << endl; - yell << _SSEPLC << endl; + cout << buffer << endl; + cout << _SSEPLC << endl; for (int i=0; i<NObsBin; i++) { if (ilo > -1 && inlo > -1 && ithc2 > -1 && lthcvar && inpc1 > -1 ) { printf(" %5.i % -#10.4g %5.i % -#10.4g % -#10.4g %5.i % -#10.4g % -#10.4g %5.i % -#10.4g % -#10.4g % -#10.4g %#18.11E %#18.11E %#9.5F %#9.5F %#9.5F", diff --git a/v2.5/toolkit/src/fnlo-tk-example.cc b/v2.5/toolkit/src/fnlo-tk-example.cc index 073416df3fddd0cc29fb072bfe8ca5e7ac9b4a73..a10450233ca7d1854dd712de146dcb5d4b18a24b 100644 --- a/v2.5/toolkit/src/fnlo-tk-example.cc +++ b/v2.5/toolkit/src/fnlo-tk-example.cc @@ -149,7 +149,7 @@ int main(int argc, char** argv) { // EPDFUncertaintyStyle ePDFUnc = kHessianCTEQCL68; // Return values are three vectors xs, dxsu, dxsl in struct XsUnc XsUncertainty XsUnc; - XsUnc = fnlo.GetScaleUncertainty(eScaleUnc); + XsUnc = fnlo.GetXsUncertainty(eScaleUnc); // XsUnc = fnlo.GetPDFUncertainty(ePDFUnc); cout << _CSEPSC << endl; diff --git a/v2.5/toolkit/src/fnlo-tk-rootout.cc b/v2.5/toolkit/src/fnlo-tk-rootout.cc index 75f6d3abf263cb639c3e7aed4ffa637004e10ece..4d9060a96c8eae345845aa25b448bf11ab291c2a 100644 --- a/v2.5/toolkit/src/fnlo-tk-rootout.cc +++ b/v2.5/toolkit/src/fnlo-tk-rootout.cc @@ -141,7 +141,7 @@ int main(int argc, char** argv) { //! --- Uncertainty choice EPDFUncertaintyStyle ePDFUnc = kPDFNone; - EAddUncertaintyStyle eAddUnc = kAddNone; + EStatUncertaintyStyle eStatUnc = kStatNone; EScaleUncertaintyStyle eScaleUnc = kScaleNone; string chunc = "none"; if (argc > 3) { @@ -465,7 +465,7 @@ int main(int argc, char** argv) { if ( iUnc==0 ) { // Back up zeroth member result when no proper uncertainty choice was made (only approx. correct for NNPDF) vector < double > xstmp = fnlo->GetCrossSection(lNorm); - XsUnc = fnlo->GetPDFUncertainty(PDFUncStyles[iPDF], lNorm); + XsUnc = fnlo->GetXsUncertainty(PDFUncStyles[iPDF], lNorm); if ( PDFUncStyles[iPDF] == kPDFNone ) XsUnc.xs = xstmp; snprintf(buffer, sizeof(buffer), " # Relative PDF uncertainties (%s %s %s)",sOrder.c_str(),PDFFiles[iPDF].c_str(),chunc.c_str()); snprintf(titlel, sizeof(titlel), "-dsigma_%s/sigma",PDFFiles[iPDF].c_str()); @@ -476,8 +476,8 @@ int main(int argc, char** argv) { if (ITabVersion < 25000) { info["fnlo-tk-rootout"] << "Table version " << ITabVersion << "too small; statistical uncertainties not available." << endl; } - eAddUnc = kAddStat; - XsUnc = fnlo->GetAddUncertainty(eAddUnc, lNorm); + eStatUnc = kStatInt; + XsUnc = fnlo->GetXsUncertainty(eStatUnc, lNorm); snprintf(buffer, sizeof(buffer), " # Relative statistical uncertainties (%s %s)",sOrder.c_str(),PDFFiles[iPDF].c_str()); snprintf(titlel, sizeof(titlel), "-dsigma_stat/sigma"); snprintf(titleu, sizeof(titleu), "+dsigma_stat/sigma"); @@ -485,7 +485,7 @@ int main(int argc, char** argv) { //! 2P scale uncertainties else if ( iUnc==2 ) { eScaleUnc = kSymmetricTwoPoint; - XsUnc = fnlo->GetScaleUncertainty(eScaleUnc, lNorm); + XsUnc = fnlo->GetXsUncertainty(eScaleUnc, lNorm); snprintf(buffer, sizeof(buffer), " # 2P Relative scale uncertainties (%s %s)",sOrder.c_str(),PDFFiles[iPDF].c_str()); snprintf(titlel, sizeof(titlel), "-dsigma_2P/sigma"); snprintf(titleu, sizeof(titleu), "+dsigma_2P/sigma"); @@ -493,7 +493,7 @@ int main(int argc, char** argv) { //! 6P scale uncertainties else if ( iUnc==3 ) { eScaleUnc = kAsymmetricSixPoint; - XsUnc = fnlo->GetScaleUncertainty(eScaleUnc, lNorm); + XsUnc = fnlo->GetXsUncertainty(eScaleUnc, lNorm); snprintf(buffer, sizeof(buffer), " # 6P Relative scale uncertainties (%s %s)",sOrder.c_str(),PDFFiles[iPDF].c_str()); snprintf(titlel, sizeof(titlel), "-dsigma_6P/sigma"); snprintf(titleu, sizeof(titleu), "+dsigma_6P/sigma"); diff --git a/v2.5/toolkit/src/fnlo-tk-yodaout.cc b/v2.5/toolkit/src/fnlo-tk-yodaout.cc index a12c8c540971cd8ff7c1e670bf82fb54dc99d98e..d05b117320dba7e22379762eb18959e2604825ec 100644 --- a/v2.5/toolkit/src/fnlo-tk-yodaout.cc +++ b/v2.5/toolkit/src/fnlo-tk-yodaout.cc @@ -1,8 +1,11 @@ ///******************************************************************** /// /// fastNLO_toolkit: fnlo-tk-yodaout -/// Program to read fastNLO tables and write out -/// QCD cross sections in YODA format for use with Rivet +/// Program to read fastNLO tables and write out QCD cross sections +/// with uncertainty in text and YODA format (for use with Rivet) +/// +/// For more explanations type: +/// ./fnlo-tk-yodaout -h /// /// K. Rabbertz, G. Sieber, S. Tyros /// @@ -21,6 +24,7 @@ #include <string> #include <vector> #include "fastnlotk/fastNLOAlphas.h" +#include "fastnlotk/fastNLOConstants.h" #include "fastnlotk/fastNLOLHAPDF.h" #include "fastnlotk/fastNLOTools.h" #include "fastnlotk/speaker.h" @@ -38,16 +42,18 @@ int main(int argc, char** argv) { using namespace fastNLO; //! namespace for fastNLO constants //! --- Set initial verbosity level - SetGlobalVerbosity(WARNING); - say::Verbosity fvol = toVerbosity()["WARNING"]; + // Could set verbosity level directly using enum here, but we'll need fVerbosity anyway for cmd line parsing + // SetGlobalVerbosity(WARNING); + string VerbosityLevel = "WARNING"; + say::Verbosity fVerbosity = toVerbosity()["WARNING"]; + SetGlobalVerbosity(fVerbosity); //--- Set verbosity level of table evaluation from command line argument if existing - string VerbosityLevel; if (argc > 11) { VerbosityLevel = (const char*) argv[11]; //! --- Reset verbosity level from here on - fvol = toVerbosity()[VerbosityLevel]; - SetGlobalVerbosity(fvol); + fVerbosity = toVerbosity()[VerbosityLevel]; + SetGlobalVerbosity(fVerbosity); } //! --- Parse commmand line @@ -70,14 +76,18 @@ int main(int argc, char** argv) { if (tablename == "-v") { fastNLOTools::PrintFastnloVersion(); return 0; + } else if (tablename == "-h") { + // Avoid suppressing usage printouts + SetGlobalVerbosity(INFO); } //! --- Print program purpose yell << _CSEPSC << endl; info["fnlo-tk-yodaout"] << "Program to read fastNLO tables and write out" << endl; - info["fnlo-tk-yodaout"] << "QCD cross sections in YODA format for use with Rivet" << endl; + info["fnlo-tk-yodaout"] << "QCD cross sections with uncertainties in text format and" << endl; + info["fnlo-tk-yodaout"] << "in YODA format for use with Rivet" << endl; info["fnlo-tk-yodaout"] << "(If compiled without YODA support only text printout is given)" << endl; infosep << _SSEPSC << endl; - info["fnlo-tk-yodaout"] << "For more explanations type:" << endl; + info["fnlo-tk-yodaout"] << "For an explanation of command line arguments type:" << endl; info["fnlo-tk-yodaout"] << "./fnlo-tk-yodaout -h" << endl; info["fnlo-tk-yodaout"] << "For version number printout type:" << endl; info["fnlo-tk-yodaout"] << "./fnlo-tk-yodaout -v" << endl; @@ -87,8 +97,8 @@ int main(int argc, char** argv) { //! --- Usage info if (tablename == "-h") { info["fnlo-tk-yodaout"] << "fastNLO YODA Writer" << endl; - yell << _SSEPSC << endl; - yell << " #" << endl; + infosep << _SSEPSC << endl; + infosep << " #" << endl; info["fnlo-tk-yodaout"] << "This program evaluates a fastNLO table and" << endl; info["fnlo-tk-yodaout"] << "prints out cross sections with statistical (if available), " << endl; info["fnlo-tk-yodaout"] << "scale, or PDF uncertainties in YODA format for use with Rivet." << endl; @@ -135,10 +145,10 @@ int main(int argc, char** argv) { man << "[xrescale]: Rescale x axis values, def. = 1" << endl; man << "[dxsrescale]: Rescale uncertainty values, def. = 1" << endl; man << "[Verbosity]: Set verbosity level of table evaluation [DEBUG,INFO,WARNING,ERROR], def. = WARNING" << endl; - yell << " #" << endl; + infosep << " #" << endl; man << "Use \"_\" to skip changing a default argument." << endl; - yell << " #" << endl; - yell << _CSEPSC << endl; + infosep << " #" << endl; + infosep << _CSEPSC << endl; return 0; } else { shout["fnlo-tk-yodaout"] << "Evaluating table: " << tablename << endl; @@ -162,7 +172,7 @@ int main(int argc, char** argv) { EScaleUncertaintyStyle eScaleUnc = kScaleNone; EPDFUncertaintyStyle ePDFUnc = kPDFNone; EAsUncertaintyStyle eAsUnc = kAsNone; - EAddUncertaintyStyle eAddUnc = kAddNone; + EStatUncertaintyStyle eStatUnc = kStatNone; string chunc; if (argc > 3) { chunc = (const char*) argv[3]; @@ -201,7 +211,7 @@ int main(int argc, char** argv) { eAsUnc = kAsGRV; shout["fnlo-tk-yodaout"] << "Showing a_s(M_Z) uncertainty with GRV evolution." << endl; } else if ( chunc == "ST" ) { - eAddUnc = kAddStat; + eStatUnc = kStatInt; shout["fnlo-tk-yodaout"] << "Showing statistical uncertainty of x section calculation." << endl; } else { error["fnlo-tk-yodaout"] << "Illegal choice of uncertainty, " << chunc << ", aborted!" << endl; @@ -334,20 +344,21 @@ int main(int argc, char** argv) { error["fnlo-tk-yodaout"] << "Too many arguments, aborting!" << endl; exit(1); } - yell << _CSEPSC << endl; + yell << _CSEPSC << endl << endl; //--- End of parsing arguments //! --- fastNLO initialisation, attach table fastNLOTable table = fastNLOTable(tablename); //! Print essential table information - if (fvol < 1) { + if (fVerbosity < toVerbosity()["WARNING"]) { table.PrintContributionSummary(0); table.Print(0); } //! Initialise a fastNLO reader instance //! Note: This also initializes the cross section to the LO/NLO one! + fastNLOLHAPDF* fnlo = NULL; if ( chunc != "AS" ) { fnlo = new fastNLOLHAPDF(table,PDFFile,0); @@ -480,34 +491,30 @@ int main(int argc, char** argv) { fnlo->CalcCrossSection(); //! If required get uncertainties (only for additive perturbative contributions) XsUncertainty XsUnc; - std::vector< std::vector<double> > xsvec; // TEST string LineName; string UncName; if ( chunc == "2P" || chunc == "6P" ) { - XsUnc = fnlo->GetScaleUncertainty(eScaleUnc, lNorm, sclfac); + XsUnc = fnlo->GetXsUncertainty(eScaleUnc, lNorm, sclfac); UncName = " # Relative scale uncertainties (" + chunc + ")"; LineName += "_dxscl"; - xsvec = fnlo->GetScaleUncertaintyVec(eScaleUnc, lNorm, sclfac); } else if ( chunc == "AS" ) { - XsUnc = fnlo->GetAsUncertainty(eAsUnc, lNorm); + XsUnc = fnlo->GetXsUncertainty(eAsUnc, lNorm); UncName = " # Relative a_s(M_Z) uncertainties (" + chunc + ")"; LineName += "_dxa_s"; } else if ( chunc == "ST" ) { - XsUnc = fnlo->GetAddUncertainty(eAddUnc, lNorm); + XsUnc = fnlo->GetXsUncertainty(eStatUnc, lNorm); UncName = " # Relative statistical uncertainties (" + chunc + ")"; LineName += "_dxst"; } else if ( chunc != "none" ) { - XsUnc = fnlo->GetPDFUncertainty(ePDFUnc, lNorm); + XsUnc = fnlo->GetXsUncertainty(ePDFUnc, lNorm); UncName = " # Relative PDF Uncertainties (" + chord + " " + PDFFile + " " + chunc + ")"; LineName += "_dxpdf"; } else { - XsUnc = fnlo->GetScaleUncertainty(kScaleNone, lNorm); + XsUnc = fnlo->GetXsUncertainty(kScaleNone, lNorm); UncName = " # Without uncertainties"; LineName += "_dxnone"; } fastNLOTools::PrintXSUncertainty(XsUnc, UncName); - cout << "VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV" << endl; - fastNLOTools::PrintXSUncertaintyVec(xsvec, UncName); //! --- Get RivetID //! For 2+-dimensions determine running number in Rivet plot name by spotting the capital letter in "RIVET_ID=" in the fnlo table