From 7feb0e21c822478db5db2f023e74047e0e9918c8 Mon Sep 17 00:00:00 2001
From: Klaus Rabbertz <klaus.rabbertz@cern.ch>
Date: Fri, 12 Jul 2024 18:14:33 +0200
Subject: [PATCH] Combined old subprocess plotting stuff into one new py3
script always producing one stack plot. Use -h for help
---
tools/plotting/fnlo-py-subprocs.py | 445 ++++++++++
tools/plotting/sp_contrib_flexible.py | 822 ------------------
tools/plotting/yb_ystar/compare_spLO_xstot.py | 377 --------
.../plotting/yb_ystar/compare_spNLO_xstot.py | 377 --------
tools/plotting/yb_ystar/kfactor_ybystar.py | 294 -------
.../yb_ystar/subproc_contrib_ybystar.py | 381 --------
6 files changed, 445 insertions(+), 2251 deletions(-)
create mode 100755 tools/plotting/fnlo-py-subprocs.py
delete mode 100755 tools/plotting/sp_contrib_flexible.py
delete mode 100755 tools/plotting/yb_ystar/compare_spLO_xstot.py
delete mode 100755 tools/plotting/yb_ystar/compare_spNLO_xstot.py
delete mode 100755 tools/plotting/yb_ystar/kfactor_ybystar.py
delete mode 100755 tools/plotting/yb_ystar/subproc_contrib_ybystar.py
diff --git a/tools/plotting/fnlo-py-subprocs.py b/tools/plotting/fnlo-py-subprocs.py
new file mode 100755
index 00000000..6f8127df
--- /dev/null
+++ b/tools/plotting/fnlo-py-subprocs.py
@@ -0,0 +1,445 @@
+#!/usr/bin/env python3
+# -*- coding:utf-8 -*-
+#
+########################################################################
+#
+# Plot the subprocess decomposition
+#
+# Created by B. Schillinger, 01.06.2018
+# Modified by K. Rabbertz, 15.04.2019
+# Simplified, standardised, and prepared for python3 by K. Rabbertz, 06/07 2024
+#
+########################################################################
+#
+# fastNLO for direct evaluation of interpolation grids
+# ATTENTION: fastNLO python extension is required for Python 3!
+from fastnlo import SetGlobalVerbosity
+from fastnlo import fastNLOLHAPDF
+import fastnlo
+import matplotlib.pyplot as plt
+import argparse
+# import glob
+# import datetime
+import os
+import re
+# import shutil
+# import sys
+# import time
+import numpy as np
+import matplotlib as mpl
+import matplotlib.gridspec as gridspec
+from matplotlib.ticker import (FormatStrFormatter, LogFormatter, NullFormatter, ScalarFormatter, AutoMinorLocator, MultipleLocator)
+from matplotlib.pyplot import cm
+# import warnings
+# warnings.filterwarnings("error")
+
+
+# We do not want any interactive plotting! Figures are saved to files instead.
+# This also avoids the ANNOYANCE of frequently missing Tkinter/tkinter (python2/3) GUI backends!
+# To produce scalable graphics for publication use eps, pdf, or svg as file format.
+# For this to work we try the Cairo backend, which can do all of these plus the raster format png.
+# If this is not usable, we fall back to the Agg backend capable only of png for nice web plots.
+# ngbackends = mpl.rcsetup.non_interactive_bk
+# print('Non GUI backends are: ', ngbackends)
+# 1st try cairo
+backend = 'cairo'
+usecairo = True
+try:
+ import cairocffi as cairo
+except ImportError:
+ try:
+ import cairo
+ except ImportError:
+ usecairo = False
+ else:
+ if cairo.version_info < (1, 11, 0):
+ # Introduced create_for_data for Py3.
+ usecairo = False
+if usecairo:
+ mpl.use('cairo')
+else:
+ backend = 'agg'
+ useagg = True
+ try:
+ mpl.use(backend, force=True)
+ print('[fnlo-py-subprocs]: Warning! Could not import cairo backend :-( Using agg instead for raster plots only!')
+ except:
+ useagg = False
+ print('[fnlo-py-subprocs]: Can not use agg backend :-(')
+ raise ImportError('[fnlo-py-subprocs]: Neither cairo nor agg backend found :-( Cannot produce any plots. Good bye!')
+ mpl.use('agg')
+
+
+# Redefine ScalarFormatter
+class ScalarFormatterForceFormat(ScalarFormatter):
+ # Override function that finds format to use.
+ def _set_format(self, vmin, vmax):
+ self.format = "%1.2f" # Give format here
+
+
+# Action class to allow comma-separated (or empty) list in options
+class SplitArgs(argparse.Action):
+ def __call__(self, parser, namespace, values, option_string=None):
+ if values:
+ setattr(namespace, self.dest, values[0].split(','))
+ else:
+ setattr(namespace, self.dest, [''])
+
+
+# Some global definitions
+_fntrans = str.maketrans({'[': '', ']': '', '(': '', ')': '', ',': ''}) # Filename 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'}
+
+# Some process dictionaries
+_process_color = {"gg": "red", "gq": "orange", "gaq": "yellow",
+ "qq": "green", "qaq": "blue", "aqaq": "purple"}
+_process_names = {
+ "gg": "Gluon-Gluon",
+ "gq": "Gluon-Quark",
+ "gaq": "Gluon-Antiquark",
+ "qq": "Quark-Quark",
+ "qaq": "Quark-Antiquark",
+ "aqaq": "Antiquark-Antiquark",
+}
+_process_texnames = {
+ "gg": r"$\mathrm{gg}$",
+ "gq": r"$\mathrm{gq}$",
+ "gaq": r"$\mathrm{g\overline{q}}$",
+ "qq": r"$\mathrm{qq}$",
+ "qaq": r"$\mathrm{q\overline{q}}$",
+ "qaq": r"$\mathrm{q\overline{q}}$",
+ "aqaq": r"$\mathrm{\overline{q}\overline{q}}$",
+}
+
+#####################################################################################
+
+
+def main():
+
+ # Define arguments & options
+ parser = argparse.ArgumentParser(epilog='', formatter_class=argparse.ArgumentDefaultsHelpFormatter)
+ # Positional arguments
+ parser.add_argument('tabfile', type=argparse.FileType('r'),
+ help='fastNLO table to be evaluated including path and extension .tab or tab.gz.')
+
+ # Optional arguments
+ parser.add_argument('-s', '--subprocs', required=False, nargs=1, type=str, action=SplitArgs,
+ default=['gg', 'gq', 'gaq', 'qq', 'qaq', 'aqaq'],
+ choices=['gg', 'gq', 'gaq', 'qq', 'qaq', 'aqaq'],
+ help='Subprocesses to be compared.\n'
+ 'Options: gg, gq, gaq, qq, qaq, aqaq (TODO: q=q, q!q, uu, uau, auau, etc.)\n'
+ 'See SelectProcesses() in fastNLOReader.cc for more information.')
+
+ parser.add_argument('-o', '--order', required=False, nargs='?', type=str, default='NNLO', choices=['LO', 'NLO', 'NNLO'],
+ help='Order to which the subprocess contributions are calculated.')
+
+ parser.add_argument('-n', '--norm-order', required=False, nargs='?', type=str, default='SAME', choices=['LO', 'NLO', 'NNLO', 'SAME'],
+ help='Normalisation order to which the subprocesses contributions are normalised.')
+
+ parser.add_argument('-p', '--pdfset', default='CT18NNLO',
+ help='PDFset to use with fastNLO table.')
+
+ parser.add_argument('-m', '--member', default=0, type=int,
+ help='Member of PDFset to use, default is 0.')
+
+ parser.add_argument('-v', '--verbosity', type=int, default=1, choices=[-1000, 0, 1, 2, 1000],
+ help='Set verbosity of fastNLO. \n'
+ 'Choices are -1000 (DEBUG), 0 (INFO), 1 (WARNING), 2 (ERROR), 1000 (SILENT)')
+ # Options for plotting details
+ parser.add_argument('--filename', default=None, type=str,
+ help='Replace default filename derived from table name and orders by given string.')
+ parser.add_argument('--format', required=False, nargs=1, type=str, action=SplitArgs,
+ help='Comma-separated list of plot formats to use: eps, pdf, png, svg. If nothing is chosen, png is used.')
+ parser.add_argument('--title', default=None, type=str,
+ help='Replace table name as default title by given string.')
+ parser.add_argument('--legtitle', default=None, type=str,
+ help='Replace order & PDF set as default title by given string.')
+ parser.add_argument('--xlabel', default=None, type=str,
+ help='Replace default x axis label from table by given string, e.g. \'$p_\mathrm{T}$ [GeV]\'')
+ parser.add_argument('--ylabel', default='Subprocess decomposition', type=str,
+ help='Replace y axis label by given string.')
+ parser.add_argument('--xmin', default=None, type=float,
+ help='Replace default x axis minimum from table by given number.')
+ parser.add_argument('--xmax', default=None, type=float,
+ help='Replace default x axis maximum from table by given number.')
+ parser.add_argument('--ymin', default=None, type=float,
+ help='Replace default y axis minimum from table by given number.')
+ parser.add_argument('--ymax', default=None, type=float,
+ help='Replace default y axis maximum from table by given number.')
+
+ # Print header
+ print("\n###########################################################################################")
+ print("# fnlo_py_subprocs:")
+ print("# Plot the subprocess decomposition")
+ print("###########################################################################################\n")
+
+ # Parse arguments
+ args = vars(parser.parse_args())
+
+ # Table and table name
+ table = args['tabfile'].name
+ tabbase = os.path.basename(table)
+ # Eliminate extensions tab and gz
+ tabname = re.sub('.gz$', '', tabbase)
+ tabname = re.sub('.tab$', '', tabname)
+ print('[fnlo-py-subprocs]: Going to read table {}.'.format(tabbase))
+
+ # PDF set
+ pdfset_ = os.path.basename(args['pdfset'])
+ pdfname = os.path.splitext(pdfset_)[0]
+ print('[fnlo-py-subprocs]: Using PDF set {}.'.format(pdfname))
+
+ # Subprocesses to investigate
+ print('[fnlo-py-subprocs]: Subprocesses that will be investigated: ')
+ print('[fnlo-py-subprocs]: ', args['subprocs'])
+
+ # Set subprocess and normalisation orders
+ order = args['order']
+ norm_order = args['norm_order']
+ if norm_order == 'SAME':
+ norm_order = order
+ print('[fnlo-py-subprocs]: The chosen subprocess order is: {}'.format(order))
+ print('[fnlo-py-subprocs]: The chosen normalisation order is: {}'.format(norm_order))
+
+ # Plot formats to use
+ formats = args['format']
+ if formats is None:
+ formats = ['png']
+ for fmt in formats:
+ if fmt not in _formats.keys():
+ print('[fnlo-py-subprocs]: Illegal format specified, aborted!')
+ print('[fnlo-py-subprocs]: Format list:', args['format'])
+ exit(1)
+
+ # Verbosity
+ ldebug = False
+ if args['verbosity'] < 1:
+ ldebug = True
+
+ ###################### Start EVALUATION with fastNLO library ###################################################
+ print('[fnlo-py-subprocs]: Starting table evaluation.')
+ SetGlobalVerbosity(args['verbosity'])
+ fnlo = fastnlo.fastNLOLHAPDF(table, args['pdfset'], args['member'])
+
+ # Get info from table
+ # Get dimensionality of table:
+ ndim = fnlo.GetNumDiffBin()
+ if ndim > 1:
+ print('[fnlo-py-subprocs]: Error! Script only tested for use with 1-dim. grids, but ndim = {}. Aborted!'.format(ndim))
+
+ # Determine maximum of subsequent orders available in table
+ max_iorder = -1
+ o_existence = []
+ for i in range(_text_to_order['NNLO']+1):
+ o_existence.append(fnlo.SetContributionON(fastnlo.kFixedOrder, i, True))
+ if o_existence[i]:
+ max_iorder = max_iorder + 1
+ else:
+ break
+ print('[fnlo-py-subprocs]: Maximum subsequent order found in table is: {}, corresponding to {}.'.format(max_iorder, _order_to_text[max_iorder]))
+ if max_iorder < _text_to_order[order] or max_iorder < _text_to_order[norm_order]:
+ print('[fnlo-py-subprocs]: Error! Maximum subsequent order available smaller than requested order. Aborted!')
+ print('[fnlo-py-subprocs]: Requested: {}, found in table: {} vs. {}'.format(max_iorder, _text_to_order[order], _text_to_order[norm_order]))
+ exit(1)
+
+ # Get labels of all dimensions in table
+ labels = fnlo.GetDimLabels()
+ # Get x label from first dimension:
+ xlabel = fnlo.GetDimLabel(0)
+ if args['xlabel'] is not None:
+ xlabel = args['xlabel']
+ print("[fnlo-py-subprocs]: x-axis label: {}".format(xlabel))
+ # y label
+ ylabel = args['ylabel']
+ print("[fnlo-py-subprocs]: y-axis label: {}".format(ylabel))
+ # Plot title
+ title = tabname
+ if args['title'] is not None:
+ title = args['title']
+ print("[fnlo-py-subprocs]: Plot title: {}".format(title))
+
+ # Get bin bounds in x
+ bin_bounds = np.array(fnlo.GetObsBinsBounds(0))
+ xlow = bin_bounds[0][0]
+ xhig = bin_bounds[-1][1]
+
+ # Set x axis limits
+ if args['xmin'] is not None:
+ xlow = args['xmin']
+ if args['xmax'] is not None:
+ xhig = args['xmax']
+ print("[fnlo-py-subprocs]: x-axis minimum: {}".format(xlow))
+ print("[fnlo-py-subprocs]: x-axis maximum: {}".format(xhig))
+
+ # Set y axis limits (ylow, yhig might be None, i.e. unset!)
+ ylow = args['ymin']
+ yhig = args['ymax']
+
+ # Evaluate cross section for normalisation
+ print("[fnlo-py-subprocs]: Normalisation order: {}".format(norm_order))
+ for j in range(0, max_iorder+1):
+ fnlo.SetContributionON(fastnlo.kFixedOrder, j, j <= _text_to_order[norm_order])
+ fnlo.CalcCrossSection()
+ xs_all_oneplot = np.array(fnlo.GetCrossSection())
+ if ldebug:
+ print("[fnlo-py-subprocs]: The normalisation cross section is:")
+ print("[fnlo-py-subprocs]: ", xs_all_oneplot)
+
+ print("[fnlo-py-subprocs]: Subprocess Order: {}".format(order))
+ for l in range(0, max_iorder+1):
+ fnlo.SetContributionON(fastnlo.kFixedOrder, l, l <= _text_to_order[order])
+ xs_subproc_list = []
+ for i in range(0, len(args['subprocs'])):
+ subproc_name = args['subprocs'][i]
+ print('[fnlo-py-subprocs]: Selected subprocess: {}'.format(subproc_name))
+ fnlo.SelectProcesses(subproc_name, True)
+ fnlo.CalcCrossSection()
+ xs_subproc = np.array(fnlo.GetCrossSection())
+ xs_subproc_list.append(xs_subproc)
+ xs_sub_oneplot = np.array(xs_subproc_list)
+
+ # Calculate fraction of subprocess at order on total cross sectio at norm_order
+ fractions_sub_oneplot = []
+ for k in range(0, len(args['subprocs'])):
+ fractions_sub_oneplot.append(np.divide(xs_sub_oneplot[k, :], xs_all_oneplot[:]))
+ fractions_oneplot = np.array(fractions_sub_oneplot)
+
+ # Checking whether contribution is positive or negative
+ pos_contr = np.array(fractions_oneplot)
+ neg_contr = np.array(fractions_oneplot)
+
+ # Go through 'fractions' = subprocess after subprocess
+ for i in range(0, len(args['subprocs'])):
+ for j in range(0, len(bin_bounds)): # Check the different bins
+ if (fractions_oneplot[i, j] < 0):
+ pos_contr[i, j] = 0.0
+ else:
+ neg_contr[i, j] = 0.0
+ if ldebug:
+ print('[fnlo-py-subprocs]: Positive contributions of {} in {} on total xs in {}:'.format(args['subprocs'][i], order, norm_order))
+ print('[fnlo-py-subprocs]: ', pos_contr[i, :], '\n')
+ print('[fnlo-py-subprocs]: Negative contributions of {} in {} on total xs in {}:'.format(args['subprocs'][i], order, norm_order))
+ print('[fnlo-py-subprocs]: ', neg_contr[i, :], '\n')
+
+ if ldebug:
+ print('[fnlo-py-subprocs]: Summary of positive and negative contributions of the subprocesses:\n')
+ print('[fnlo-py-subprocs]: Positive:', pos_contr, '\n')
+ print('[fnlo-py-subprocs]: Negative:', neg_contr, '\n')
+
+ # Create stackplot that compares the contribution of the subprocesses to the total cross section
+ plt.close()
+
+ fig = plt.figure(figsize=(8, 7))
+ ax0 = fig.add_subplot(111)
+ # Number of bins via bin_bounds variable
+ # For tracking the current 'height' of stackplot, in order to staple and not overlap the values for each subprocess
+ y0_pos = np.zeros([len(bin_bounds)])
+ # Same as above for the negative contributions, plotted below x-axis
+ y0_neg = np.zeros([len(bin_bounds)])
+ patches = [] # Needed later for legend
+
+ # Go through fractions, subproc by subproc
+ for i in range(0, len(args['subprocs'])):
+ cl = _process_color[args['subprocs'][i]]
+ lb = _process_texnames[args['subprocs'][i]]
+
+ # Positive contributions
+ ax0.fill_between(bin_bounds.flatten(), steppify_bin(y0_pos), steppify_bin(y0_pos+pos_contr[i, :]), facecolor=cl, alpha=0.75)
+ # Calculate new height of stackplot
+ y0_pos += pos_contr[i, :]
+
+ # Negative contributions
+ ax0.fill_between(bin_bounds.flatten(), steppify_bin(y0_neg), steppify_bin(y0_neg+neg_contr[i, :]), facecolor=cl, hatch='X', alpha=0.75)
+ # Calculate new height below x-axis
+ y0_neg += neg_contr[i, :]
+
+ # Patch for subprocess i (for the legend)
+ patches.append(mpl.patches.Rectangle((0, 0), 0, 0, color=cl, label=lb, alpha=0.75))
+
+ # Derive y-axis limits from content, if not set otherwise
+ if ylow is None:
+ ylow = np.amin(y0_neg)
+ if yhig is None:
+ yhig = np.amax(y0_pos)
+ print("[fnlo-py-subprocs]: y-axis minimum: {}".format(ylow))
+ print("[fnlo-py-subprocs]: y-axis maximum: {}".format(yhig))
+
+ # Set axes limits
+ ax0.set_xlim(xlow, xhig)
+ ax0.set_ylim(ylow, yhig)
+ ax0.set_xscale('log', nonpositive='clip')
+ # ax0.autoscale(enable=True, axis='x', tight=True)
+
+ for i in range(len(args['subprocs'])-1, -1, -1):
+ ax0.add_patch(patches[i])
+
+ # Dotted line at xs_sub/xs = 1 = 100%
+ plt.axhline(y=1, xmin=0, xmax=1, color='k', linestyle='dotted')
+ # Dotted line at xs_sub/xs = 0 = 0%
+ plt.axhline(y=0, xmin=0, xmax=1, color='k', linestyle='-')
+
+ plt.title(label=title, fontsize=16, loc='center')
+ ax0.set_xlabel(xlabel, fontsize=16, loc='right')
+ ax0.set_ylabel(ylabel, fontsize=16, loc='top')
+ plt.xticks(fontsize=16)
+ plt.yticks(fontsize=16)
+ # TODO: ticks & labels
+ # ax0.get_xaxis().set_major_locator(mpl.ticker.LogLocator(subs=(1.0, 2.0, 4.0, 7.0)))
+ # ax0.get_xaxis().set_major_formatter(mpl.ticker.LogFormatter())
+ legtitle = '{} / {}'.format(order, norm_order)
+ if order == norm_order:
+ legtitle = '{}'.format(order)
+ legtitle = legtitle+'\n'+args['pdfset']
+ if args['legtitle'] is not None:
+ legtitle = args['legtitle']
+ plt.legend(title=legtitle, loc='best') # reverse=True requires matplotlib 3.7
+
+ plt.tight_layout()
+ fig.tight_layout()
+
+ # Naming of the plots
+ figname = tabname
+ if args['filename'] is not None:
+ figname = args['filename']
+ figname = figname + '_' + order + '-vs-' + norm_order
+
+ # Store in all requested formats
+ for fmt in formats:
+ filename = '%s.%s' % (figname, fmt)
+ # ban characters defined in _fntrans for filenames (parentheses and commas)
+ filename = filename.translate(_fntrans)
+ fig.savefig(filename, bbox_inches='tight')
+ print('[fnlo-py-subprocs]: Plot saved as:', filename)
+ plt.close(fig)
+
+
+##################################################
+# function to make better uncertainty plot (shaded)
+# could maybe be replaced by .flatten() ??
+
+
+def steppify_bin(arr, isx=False):
+ """
+ Produce stepped array of arr, needed for example for stepped fill_betweens.
+ Pass all x bin edges to produce stepped x arr and all y bincontents to produce
+ stepped bincontents representation
+ steppify_bin([1,2,3], True)
+ -> [1,2,2,3]
+ steppify_bin([5,6])
+ -> [5,5,6,6]
+ """
+ if isx:
+ # newarr = np.array(zip(arr[:-1], arr[1:])).ravel()
+ newarr = np.array(list(zip(arr[:-1], arr[1:]))).flatten()
+ else:
+ # newarr = np.array(zip(arr, arr)).ravel()
+ newarr = np.array(list(zip(arr, arr))).flatten()
+ return newarr
+
+
+#################################################
+
+if __name__ == '__main__':
+ main()
diff --git a/tools/plotting/sp_contrib_flexible.py b/tools/plotting/sp_contrib_flexible.py
deleted file mode 100755
index f126d255..00000000
--- a/tools/plotting/sp_contrib_flexible.py
+++ /dev/null
@@ -1,822 +0,0 @@
-#! /usr/bin/env python
-
-###############################################################
-#Script for comparing the contribution of several subprocesses
-#to the total cross section.
-#Either (per default) as stackplot or for single process (-e).
-#
-#Choice of specific orders for subprocess (-o) and normalisation
-#(-n) is possible here.
-#If no specific order is chosen --> all 5 plots are produced.
-#
-#General script, not customized for any certain tables.
-###############################################################
-
-
-import argparse
-import shutil
-import numpy as np
-import os
-import time, datetime
-import matplotlib
-matplotlib.use('Agg') #
-import matplotlib.pyplot as plt
-import matplotlib.gridspec as gridspec
-from matplotlib.pyplot import cm
-import sys
-
-import matplotlib.style
-matplotlib.style.use('classic') #to produce plots that look like they used to in matplotlib 1.x
-
-import fastnlo
-
-
-
-def main():
-
- parser = argparse.ArgumentParser()
-
- #make it possible to use different options when running
- parser.add_argument('-t','--table', default='table.tab', required=True,
- help='FastNLO tables that shall be evaluated.')
- parser.add_argument('-p', '--pdfset', default='CT14nlo',
- help='PDFset(s) to evaluate fastNLO tables.')
- parser.add_argument('-m', '--member', default=0, type=int,
- help='Member of PDFset, default is 0.')
-
- parser.add_argument('-s', '--subproc', default=['all'], type=str, nargs='*',
- help='Subprocesses that shall be compared. \n'
- 'Options: gg, gq, gaq, qq, q=q, q!q, uu, auu, auau, etc. \n'
- 'See SelectProcesses() in fastNLOReader.cc for more information.')
-
- parser.add_argument('-e', '--extraplot', default=False, const=True, type=bool, nargs='?',
- help='Produces extra single-process-plots if option -e is chosen. \n'
- 'Otherwise stackplot (containing all chosen processes) per default.')
-
- parser.add_argument('-f', '--filename', type=str, #default='stackplot.png',
- help='Optional name for output file. \n'
- 'If nothing is chosen: string containing tablename and orders.')
-
- parser.add_argument('-x', '--xaxis', default='x quantity', type=str,
- help='Set label for x-axis. Default is \'x quantity\'.')
-
- #parser.add_argument('-c', '--constyaxis', default=False, const=True, type=bool, nargs='?',
- # help='Set constant limits for y-axis from -0.4 to 1.2. \n'
- # 'Otherwise (per default) y-axis will be flexibly adjusted to the data.')
-
- parser.add_argument('-c', '--constyaxis', default=None, type=float, nargs='*', #nargs=2
- help='Set constant limits for y-axis by giving lower and upper limit. \n'
- 'Otherwise (per default) -c sets y-axis from -0.2 to 1.2. \n'
- 'Without -c, y-axis will be flexibly adjusted to the data.')
-
-
- #flexible choice of order for the subprocess contribution (-o) and normalisation (-n)
- #(NLO to NLO+LO for instance)
- parser.add_argument('-o', '--order', type=int,
- help='Order in which the (single) subprocess contributions to XS are calculated. \n'
- 'Type 0 for LO, 1 for NLO, 2 for NNLO, etc. '
- 'If nothing is chosen order=normorder. '
- 'If normorder is not chosen either, all five plots are produced.')
-
- parser.add_argument('-n', '--normorder', type=int,
- help='Normalisation order in which the XS with all subprocesses included is calculated. \n'
- 'Type 0 for LO, 1 for NLO, 2 for NNLO, etc. '
- 'If nothing is chosen normorder=order. '
- 'If order is not chosen either, all five plots are produced.')
-
-
-################# TAKING CARE OF THE INPUT #################################
- #parse arguments
- args = vars(parser.parse_args())
- namesp = parser.parse_args() #later used for checking if -c is set or not
-
- #table name
- table_ = os.path.basename(args['table'])
- tablename = os.path.splitext(table_)[0]
- tablename = os.path.splitext(tablename)[0] #because yb tables have .tab.gz ending (getting rid of both here)
- print('\n')
- print('Table: ', tablename, '\n')
-
- #pdfset name
- pdfset_ = os.path.basename(args['pdfset'])
- pdfname = os.path.splitext(pdfset_)[0]
- print('PDF Set: ', pdfname, '\n')
-
- #Which subprocesses are evaluated?
- print('Subprocesses that are investigated: ')
- print(args['subproc'])
-
-
- #handling of subprocess order and normalisation order:
-
- #for labeling purposes
- print('subprocess order: ', args['order'])
- print('normalisation order: ', args['normorder'])
-
- #subprocess order (either chosen or per default all combinations)
- sp_order_all = ['LO', 'NLO', 'NNLO']
- if (args['order'] is not None) and (args['order'] in range(0,3)):
- sp_order = sp_order_all[args['order']]
- elif (args['order'] is None) and (args['normorder'] is not None):
- sp_order = sp_order_all[args['normorder']]
- print("Set subprocess order to chosen normalisation order.")
- else:
- print("No certain (or valid) subprocess order given.")
- sp_order = sp_order_all #in case of no specification --> sp_order is array, will produce all 5 plots
-
- #normalisation order (as above)
- norm_order_all = ['LO', 'NLO', 'NNLO']
- if (args['normorder'] is not None) and (args['normorder'] in range(0,3)):
- norm_order = norm_order_all[args['normorder']]
- elif (args['normorder'] is None) and (args['order'] is not None):
- norm_order = norm_order_all[args['order']]
- print("Set normalisation order to chosen subprocess order.")
- else:
- print("No certain (or valid) normalisation order given.")
- norm_order = norm_order_all #watch out! here norm_order becomes an array
-
-
-
-############################ EVALUATING ####################################
- fnlo = fastnlo.fastNLOLHAPDF(args['table'], args['pdfset'], args['member'])
-
- #dictionary with settings for the different orders (when switching them on/off)
- b_order = { "LO": [True, False, False], "NLO": [True, True, False], "NNLO": [True, True, True] }
-
- '''
- if (norm_order < sp_order):
- print("Choice of orders not reasonable.") #depends on purpose of the plot. (maybe useful for kfactors if args['subproc']='all')
- else:
- #usual plotting
- '''
-
- #labeling of the x-axis
- #get dimensionality of table:
- ndim = fnlo.GetNumDiffBin()
- print("\n", "Dimensions: ", ndim)
- #if (ndim > 1):
-
- #get labels of all the dimensions:
- labels = fnlo.GetDimLabels()
- print("Labels: ", labels)
-
- #get label of first dimension:
- xlabel = fnlo.GetDimLabel(0)
- print("x-label: ", xlabel)
-
-
-
-
- print("\n")
- #check whether all 5 plots are needed or just one certain --> evaluate accordingly
- if isinstance(norm_order, list) and isinstance(sp_order, list): #actually one condition to check should be enough in both cases
- #produce all 5 plots (this also happens if both input orders were invalid...)
- print('Start table evaluation for creating 5 plots. \n')
-
- #cross section with every process included
- xs_all_list = []
- for n in norm_order:
- print("\n", "n: ", n)
- for j in range(0,3):
- fnlo.SetContributionON(fastnlo.kFixedOrder, j, b_order[n][j]);
- #for order, setting in sorted(b_order.iteritems()):
- # index = 0
- # for j in setting:
- # fnlo.SetContributionON(fastnlo.kFixedOrder, index, j)
- # index+=1
-
- fnlo.CalcCrossSection()
- xs_all_list.append(fnlo.GetCrossSection())
- xs_all = np.array(xs_all_list)
- print('\n')
- print('Cross section with all subprocesses xs_all: \n')
- print(xs_all, '\n \n')
-
- #cross section of single processes
- xs_sub_list = []
- for s in sp_order: #go through lo, nlo, nnlo
- print("\n", "Subprocess order s: ", s)
- xs_subproc_list = [] #will be emptied after each iteration
- for l in range(0,3): #go through the settings for this order s
- fnlo.SetContributionON(fastnlo.kFixedOrder, l, b_order[s][l])
-
- for i in range(0, len(args['subproc'])):
- subproc_name = args['subproc'][i]
- print('Selected subprocess: %s' %subproc_name)
- fnlo.SelectProcesses(subproc_name, True) #optional second argument: symmetric==True per default
- fnlo.CalcCrossSection()
- xs_subproc = np.array(fnlo.GetCrossSection())
- xs_subproc_list.append(xs_subproc) #append xs of chosen subprocess in that certain order
- print('XS for subprocess %s: \n' %args['subproc'][i])
- print(xs_subproc, '\n \n') #should be one line in xs_sub[k,i,:]
- #all selected processes in lo or nlo or nnlo
- xs_sub_list.append(np.array(xs_subproc_list))
- xs_sub = np.array(xs_sub_list)
- print("XS of the selected processes in different pTZ-regions, in LO, NLO and NNLO: ")
- print("From outmost level to innermost: order, subprocess, pTZ-bin. \n")
- print(xs_sub, '\n \n')
-
-
- #### Investigating the fractions ####
- #Now, check all five combinations of lo, nlo, nnlo that shall be plotted
- #How much do the selected subprocesses contribute to the total xs?
-
- #Three plots where subprocess is of same order as normalisation is
- fractions_eq_order = []
- for i in range(0,3):
- fractions_eq_order.append(np.divide(xs_sub[i,:],xs_all[i]))
- fractions_eq = np.array(fractions_eq_order)
- print('fractions, equal order of subprocess and normalisation: \n', fractions_eq, '\n')
-
- '''
- #output for checking what has been calculated
- for n in range(0, 3):
- for i in range(0, len(fractions_eq[0,:])):
- print('fraction of %s on total XS, both in %s: ' %(args['subproc'][i], sp_order_all[n]))
- print(fractions_eq[n,i,:], '\n \n')
- '''
-
- #make the two additional plots for (LO to NNLO) and (NLO to NNLO)
- fractions_diff_order = []
- for s in range(0,2): #s is current subprocess order (LO, NLO)
- fractions_diff_order.append(np.divide(xs_sub[s,:,:], xs_all[2,:]))
- fractions_diff = np.array(fractions_diff_order)
- print('fractions, different order for subprocess/normalisation: LO/NNLO, NLO/NNLO: ')
- print(fractions_diff, '\n')
-
- #combine these arrays to 'fractions'-array (5 outmost indices, one for each plot)
- fractions = np.append(fractions_eq, fractions_diff, axis=0)
- print("------------------------------------------------------------------")
- print('Fractions for all the plots to be produced, first index indicates:')
- print('LO/LO, NLO/NLO, NNLO/NNLO, LO/NNLO, NLO/NNLO')
- print(fractions, '\n')
- print("------------------------------------------------------------------ \n")
-
-
- #checking whether contribution is positive or negative
- #(process all 5 "lines" of fractions-array)
- pos_contr = np.array(fractions)
- neg_contr = np.array(fractions)
- bin_bounds = np.array(fnlo.GetObsBinsBounds(0))
-
- for p in range(0, len(fractions[:,0,0])): #should be range(0,6), as there are 5 plots
- for i in range(0, len(args['subproc'])): #go through 'fractions' (second index line by line) = subprocess after subprocess
- for j in range(0, len(bin_bounds)): #check the different pTZ bins
- if (fractions[p,i,j] < 0):
- pos_contr[p,i,j] = 0.0
- else:
- neg_contr[p,i,j] = 0.0
- #print('Positive contributions of %s on total XS in plot %s: \n' %(args['subproc'][i], p))
- #print(pos_contr[p,i,:], '\n')
- #print('Negative contributions of %s on total XS in plot %s: \n' %(args['subproc'][i], p))
- #print(neg_contr[p,i,:], '\n \n')
-
- print('Summary of positive and negative contributions of the subprocesses for the 5 plots: \n')
- print('pos_contr: \n', pos_contr, '\n')
- print('neg_contr: \n', neg_contr, '\n')
-
- ############################ Five-Plots-PLOTTING ############################################
- #Bins for x-axis (according to first pdf (0))
- bin_bounds = np.array(fnlo.GetObsBinsBounds(0))
- #print(bin_bounds.T) #bin_bounds.T[0] = lower bounds, bin_bounds.T[1] upper bounds
- print('bin bounds:', bin_bounds.flatten())
-
- #for naming the plots afterwards and give correct label to y-axis:
- plots_order = np.array([["LO", "LO"], ["NLO", "NLO"], ["NNLO", "NNLO"], ["LO", "NNLO"], ["NLO", "NNLO"]], dtype=str)
- print('The five plots will contain: ')
- print(plots_order, '\n')
-
- #Check whether comparison-plot or single-plot is required:
- print('Single plots for each process?: ', args['extraplot'])
- if (args['extraplot']==False): #stackplot that compares the contribution of the subprocesses to the total XS
- plt.close()
-
- for p in range(0, 5): #p is index that says which plot 0 to 5 is produced
- #plt.close('all')
- fig0 = plt.figure(figsize=(8,7))
- ax0 = fig0.add_subplot(111)
- #number of bins via bin_bounds variable
- y0_pos = np.zeros([len(bin_bounds)]) #for tracking the current 'height' of stackplot, in order to staple and not overlap the values for each subprocess
- y0_neg = np.zeros([len(bin_bounds)]) #same as above for the negative contributions, plotted below x-axis
- print('- - - -')
- print('Start plotting of plot %s \n' %p)
- color = iter(cm.rainbow(np.linspace(0,1,len(args['subproc']))))
- patches = [] #later needed for legend
-
- for i in range(0, len(args['subproc'])): #go through fractions, subproc by subproc
- c = next(color)
- #positive contributions
- ax0.fill_between(bin_bounds.flatten(), steppify_bin(y0_pos), steppify_bin(y0_pos+pos_contr[p, i]), #label=labeling[args['subproc'][i]],
- facecolor=c, alpha=0.4)
- y0_pos += pos_contr[p, i] #calculate new height of stackplot
-
- #negative contributions
- ax0.fill_between(bin_bounds.flatten(), steppify_bin(y0_neg), steppify_bin(y0_neg+neg_contr[p, i]), #label=labeling[args['subproc'][i]],
- facecolor=c, hatch='X', alpha=0.4)
- y0_neg += neg_contr[p, i] #new 'height' below x-axis
-
- #patch for subprocess i (for the legend) ##labeling happens here
- patches.append(matplotlib.patches.Rectangle((0, 0), 0, 0, color=c,
- label=args['subproc'][i], alpha=0.4))
- ax0.add_patch(patches[i])
-
- #get limits for axes:
- xlim = ax0.get_xlim()
- ylim = ax0.get_ylim()
-
- print('xlim:', xlim)
- print('ylim:', ylim, '\n')
-
- xmin_ = xlim[0]
- xmax_ = xlim[1]
- #ymin = ylim[0]
- #ymax = ylim[1]
- ymin_ = np.amin(y0_neg)+0.3*np.amin(y0_neg)
- ymax_ = np.amax(y0_pos)+0.1*np.amax(y0_pos)
-
- #settings for the whole stackplot
- ax0.set_xscale('log', nonposx='clip')
-
-
- '''
- if args['constyaxis']==True:
- ax0.axis([0, xmax, -0.20, 1.20])
- elif args['constyaxis']==False:
- ax0.set_xlim(xmin_, xmax_) ##
- ax0.set_ylim(ymin_, ymax_) ##set axis limit according to data
- '''
-
- #print("namesp", namesp, "\n")
-
- if (args['constyaxis'] is not None): #take limits that user has chosen
- if (len(args['constyaxis'])==2):
- ylow = args['constyaxis'][0]
- yhigh = args['constyaxis'][1]
- if ylow < yhigh:
- ax0.axis([0, xmax_, ylow, yhigh])
- else:
- print(" # ERROR! when trying to set y-limits.")
- print("Check choice of y-limits. First value should be lower than second.")
- print("Syntax: -c <lower ylimit> <higher ylimit>.")
- break
- elif (len(args['constyaxis'])==0): #set constant y-axis
- ax0.axis([0, xmax_, -0.20, 1.20])
- else:
- print(" # ERROR! when trying to set y-limits.")
- print("Something went wrong. Check input after -c.")
- print("Possible: no value, two values or don't use -c at all for flexible axis. \n"
- break
- elif (args['constyaxis'] is None): #use flexible y-axis
- ax0.set_xlim(xmin_, xmax_) ##
- ax0.set_ylim(ymin_, ymax_) ##set axis limit according to data
-
-
- ax0.autoscale(enable=True, axis='x', tight=True)
- plt.axhline(y=1, xmin=0, xmax=1, color='k', linestyle='dotted') #dotted line at xs_sub/xs=1=100%
- plt.axhline(y=0, xmin=0, xmax=1, color='k', linestyle='-') #line at xs_sub/xs_tot=0
-
- plt.title(x=0.5, y=1.06, s='%s' %(table_), fontsize=22)
-
- ax0.set_xlabel('%s' %xlabel, fontsize=16, horizontalalignment='right')
- ax0.xaxis.set_label_coords(1.0, -0.04)
- ax0.set_ylabel(r'$\frac{\mathrm{XS_{sp, %s}}}{\mathrm{XS_{tot, %s}}}$' %(plots_order[p,0], plots_order[p,1]), rotation=0, fontsize=24)
- ax0.yaxis.set_label_coords(-0.16, 0.9)
-
- ax0.text(0.96, 0.03, args['pdfset']+', %s to %s' %(plots_order[p,0], plots_order[p,1]), alpha=0.6, transform=ax0.transAxes, fontsize=14, verticalalignment='bottom', horizontalalignment='right')
- plt.legend(bbox_to_anchor=(1.02, 1), loc='upper left') #'upper left' refers to the legend box
-
- plt.tight_layout()
- fig0.tight_layout()
-
- #naming of the plots
- if args['filename'] is not None:
- stackplotname = '%s_%s.png' %(args['filename'], plots_order[p,0]+'_'+plots_order[p,1])
- else:
- stackplotname = '%s.stackplot_%s.png' %(tablename, plots_order[p,0]+'_'+plots_order[p,1])
- fig0.savefig(stackplotname, bbox_inches='tight')
-
- print('Stackplot for subprocess comparison saved as: %s' %(stackplotname), '\n \n')
-
-
- else:
- print("----------------------------------------------------------")
- print("Create individual plot for each single subprocess \n")
-
- #individual plot for each subproc (NumberOfPlots = 5*NumberOfSubprocesses)
- #create folder for saving the plots, will be overwritten in case it already exists (avoid countless folders while testing)
- current_dir = os.getcwd()
-
- if args['filename'] is not None:
- prename = args['filename']
- else: prename = tablename
-
- final_dir = os.path.join(current_dir, r'%s_subproc' %(prename))
-
- if os.path.exists(final_dir):
- shutil.rmtree(final_dir) #remove final_dir if it already exists
- os.makedirs(final_dir)
-
- #loop where single plot for each of the selected subprocesses is created
-
- #plot the subprocess contributions for each chosen process individually
- #do this for all the five combinations: LO/LO, NLO/NLO, NNLO/NNLO, LO/NNLO, NLO/NNLO
- for p in range(0, 5): #p is index that says which plot 0 to 5 is produced
- #plt.close('all')
- color = iter(cm.rainbow(np.linspace(0,1,len(args['subproc']))))
- fig0 = plt.figure(figsize=(8,7))
- ax0 = fig0.add_subplot(111)
-
- for i in range(0, len(args['subproc'])):
- c = next(color)
- plt.close()
- fig0 = plt.figure(figsize=(8, 7))
- ax0 = fig0.add_subplot(111)
-
- #plot the fractions
- y0 = np.zeros([len(bin_bounds)])
-
- ax0.fill_between(bin_bounds.flatten(), steppify_bin(y0), steppify_bin(pos_contr[p,i]),
- facecolor=c, alpha=0.4)
- ax0.fill_between(bin_bounds.flatten(), steppify_bin(y0), steppify_bin(neg_contr[p,i]),
- facecolor=c, hatch='X', alpha=0.4)
-
- #get limits for axes:
- xlim = ax0.get_xlim()
- ylim = ax0.get_ylim()
-
-
- ax0.set_xscale('log', nonposx='clip')
-
- xmin_ = xlim[0]
- xmax_ = xlim[1]
- ymin_ = ylim[0]+0.3*ylim[0]
- ymax_ = ylim[1]+0.1*ylim[1]
-
- if (args['constyaxis'] is not None): #take limits that user has chosen
- if (len(args['constyaxis'])==2):
- ylow = args['constyaxis'][0]
- yhigh = args['constyaxis'][1]
- if ylow < yhigh:
- ax0.axis([0, xmax_, ylow, yhigh])
- else:
- print(" # ERROR! when trying to set y-limits.")
- print("Check choice of y-limits. First value should be lower than second.")
- print("Syntax: -c <lower ylimit> <higher ylimit>.")
- break
- elif (len(args['constyaxis'])==0): #set constant y-axis
- ax0.axis([0, xmax_, -0.20, 1.20])
- else:
- print(" # ERROR! when trying to set y-limits.")
- print("Something went wrong. Check input after -c.")
- print("Possible: no value, two values or don't use -c at all for flexible axis. \n")
- break
- elif (args['constyaxis'] is None): #use flexible y-axis
- ax0.set_xlim(xmin_, xmax_) ##
- ax0.set_ylim(ymin_, ymax_) ##set axis limit according to data
-
-
- #ax0.axis([bin_bounds.flatten()[0], xlim[1], ylim[0]-0.1*ylim[0], ylim[1]+0.1*ylim[1]]) #leave y-limits fixed for better comparison between the subprocess-plots
- ax0.autoscale(enable=True, axis='x', tight=True)
- plt.axhline(y=1, xmin=0, xmax=1, color='k', linestyle='dotted')
- plt.axhline(y=0, xmin=0, xmax=1, color='k', linestyle='-')
-
-
- titlename = args['subproc'][i]
-
- plt.title(x=0.5, y=1.01, s=titlename, fontsize=20)
- ax0.set_xlabel('%s' %xlabel, fontsize=16, horizontalalignment='right')
- ax0.xaxis.set_label_coords(1.0, -0.04)
- ax0.set_ylabel(r'$\frac{\mathrm{XS_{sp, %s}}}{\mathrm{XS_{tot, %s}}}$' %(plots_order[p,0], plots_order[p,1]), rotation=0, fontsize=24)
- ax0.yaxis.set_label_coords(-0.16, 0.9)
- ax0.text(0.96, 0.03, args['pdfset']+', %s to %s' %(plots_order[p,0], plots_order[p,1]), transform=ax0.transAxes, alpha=0.6, fontsize=14, verticalalignment='bottom', horizontalalignment='right')
- ax0.text(0.96, 0.08, '%s' %table_, transform=ax0.transAxes, alpha=0.6, fontsize=14, verticalalignment='bottom', horizontalalignment='right')
-
- #plt.legend() #location parameter loc='upper right'
-
- plt.tight_layout()
- fig0.tight_layout()
- plotname = '%s_%s_singlefrac_%s_%s.png' %(tablename, args['subproc'][i], plots_order[p,0], plots_order[p,1])
-
- #save the plots
- fig0.savefig(os.path.join(final_dir,plotname))
- print('Plot saved in %s' %(final_dir))
- print('fraction plot saved as: %s' %(plotname), "\n \n")
-
-
-
- elif isinstance(norm_order, str) and isinstance(sp_order, str):
- #produce only demanded plot for chosen orders
- print('Start table evaluation for creating 1 plot. \n')
-
- print("Normalisation Order: ", norm_order)
- for j in range(0,3):
- fnlo.SetContributionON(fastnlo.kFixedOrder, j, b_order[norm_order][j])
- fnlo.CalcCrossSection()
- xs_all_oneplot = np.array(fnlo.GetCrossSection())
- print("\n","Cross Section with all subprocesses included in %s: \n" %norm_order)
- print(xs_all_oneplot, "\n \n")
-
- print("Subprocess Order: ", sp_order)
- for l in range(0,3):
- fnlo.SetContributionON(fastnlo.kFixedOrder, l, b_order[sp_order][l])
- xs_subproc_list = []
- for i in range(0, len(args['subproc'])): #go through subprocesses
- subproc_name = args['subproc'][i]
- print('Selected subprocess: %s' %subproc_name)
- fnlo.SelectProcesses(subproc_name, True) #optional second argument: symmetric==True per default
- fnlo.CalcCrossSection()
- xs_subproc = np.array(fnlo.GetCrossSection())
- xs_subproc_list.append(xs_subproc)
- print("XS for subprocess %s in %s: \n" %(subproc_name, sp_order)
- print(xs_subproc, '\n \n')
- xs_sub_oneplot = np.array(xs_subproc_list)
- print("XS of the selected processes in different pTZ-regions in %s: " %sp_order)
- print(xs_sub_oneplot, "\n \n")
-
- #calculate fraction of subprocess in sp_order on total XS in norm_order
- fractions_sub_oneplot = []
- for k in range(0, len(args['subproc'])):
- fractions_sub_oneplot.append(np.divide(xs_sub_oneplot[k,:], xs_all_oneplot[:]))
- fractions_oneplot = np.array(fractions_sub_oneplot)
- print("--------------------------------------------------------------")
- print("Fractions for subprocesses in %s to total XS in %s" %(sp_order, norm_order))
- print(fractions_oneplot, "\n")
- print("--------------------------------------------------------------")
-
- #checking whether contribution is positive or negative
- pos_contr = np.array(fractions_oneplot)
- neg_contr = np.array(fractions_oneplot)
- bin_bounds = np.array(fnlo.GetObsBinsBounds(0))
-
- for i in range(0, len(args['subproc'])): #go through 'fractions' = subprocess after subprocess
- for j in range(0, len(bin_bounds)): #check the different pTZ bins
- if (fractions_oneplot[i,j] < 0):
- pos_contr[i,j] = 0.0
- else:
- neg_contr[i,j] = 0.0
- print('Positive contributions of %s in %s on total XS in %s: \n' %(args['subproc'][i], sp_order, norm_order))
- print(pos_contr[i,:], '\n')
- print('Negative contributions of %s in %s on total XS in %s: \n' %(args['subproc'][i], sp_order, norm_order))
- print(neg_contr[i,:], '\n \n')
-
- print('Summary of positive and negative contributions of the subprocesses: \n')
- print("positive: \n", pos_contr, '\n')
- print("negative: \n", neg_contr, '\n')
-
-
- ######################## One-Plot-PLOTTING #########################################
- #Bins for x-axis (according to first pdf (0))
- bin_bounds = np.array(fnlo.GetObsBinsBounds(0))
- print('bin bounds:', bin_bounds.flatten())
-
- #Check whether comparison-plot or single-plot is required:
- print('Single plots for each process?: ', args['extraplot'])
- if (args['extraplot']==False): #stackplot that compares the contribution of the subprocesses to the total XS
- plt.close()
-
- fig0 = plt.figure(figsize=(8,7))
- ax0 = fig0.add_subplot(111)
- #number of bins via bin_bounds variable
- y0_pos = np.zeros([len(bin_bounds)]) #for tracking the current 'height' of stackplot, in order to staple and not overlap the values for each subprocess
- y0_neg = np.zeros([len(bin_bounds)]) #same as above for the negative contributions, plotted below x-axis
-
- color = iter(cm.rainbow(np.linspace(0,1,len(args['subproc']))))
- patches = [] #later needed for legend
-
- for i in range(0, len(args['subproc'])): #go through fractions, subproc by subproc
- c = next(color)
- #positive contributions
- ax0.fill_between(bin_bounds.flatten(), steppify_bin(y0_pos), steppify_bin(y0_pos+pos_contr[i,:]), #label=labeling[args['subproc'][i]],
- facecolor=c, alpha=0.4)
- y0_pos += pos_contr[i, :] #calculate new height of stackplot
-
- #negative contributions
- ax0.fill_between(bin_bounds.flatten(), steppify_bin(y0_neg), steppify_bin(y0_neg+neg_contr[i, :]), #label=labeling[args['subproc'][i]],
- facecolor=c, hatch='X', alpha=0.4)
- y0_neg += neg_contr[i, :] #new 'height' below x-axis
-
- #patch for subprocess i (for the legend) ##labeling happens here
- patches.append(matplotlib.patches.Rectangle((0, 0), 0, 0, color=c,
- label=args['subproc'][i], alpha=0.4))
- ax0.add_patch(patches[i])
-
- #get limits for axes:
- xlim = ax0.get_xlim()
- ylim = ax0.get_ylim()
-
- print('xlim:', xlim)
- print('ylim:', ylim, '\n')
-
- xmin_ = xlim[0]
- xmax_ = xlim[1]
- ymin_ = np.amin(y0_neg)+0.3*np.amin(y0_neg)
- ymax_ = np.amax(y0_pos)+0.1*np.amax(y0_pos)
-
-
- if (args['constyaxis'] is not None): #take limits that user has chosen
- if (len(args['constyaxis'])==2):
- ylow = args['constyaxis'][0]
- yhigh = args['constyaxis'][1]
- if ylow < yhigh:
- ax0.axis([0, xmax_, ylow, yhigh])
- else:
- print(" # ERROR! when trying to set y-limits.")
- print("Check choice of y-limits. First value should be lower than second.")
- print("Syntax: -c <lower ylimit> <higher ylimit>.")
- break
- elif (len(args['constyaxis'])==0): #set constant y-axis
- ax0.axis([0, xmax_, -0.20, 1.20])
- else:
- print(" # ERROR! when trying to set y-limits.")
- print("Something went wrong. Check input after -c.")
- print("Possible: no value, two values or don't use -c at all for flexible axis. \n")
- break
- elif (args['constyaxis'] is None): #use flexible y-axis
- ax0.set_xlim(xmin_, xmax_) ##
- ax0.set_ylim(ymin_, ymax_) ##set axis limit according to data
-
- ax0.set_xscale('log', nonposx='clip')
- ax0.autoscale(enable=True, axis='x', tight=True)
-
- plt.axhline(y=1, xmin=0, xmax=1, color='k', linestyle='dotted') #dotted line at xs_sub/xs=1=100%
- plt.axhline(y=0, xmin=0, xmax=1, color='k', linestyle='-') #line at xs_sub/xs_tot=0
-
- plt.title(x=0.5, y=1.06, s='%s' %(table_), fontsize=22)
-
- #ax0.set_xlabel('x quantity', fontsize=20)
- ax0.set_xlabel('%s' %xlabel, fontsize=16, horizontalalignment='right')
- ax0.xaxis.set_label_coords(1.0, -0.04)
- ax0.set_ylabel(r'$\frac{\mathrm{XS_{sp, %s}}}{\mathrm{XS_{tot, %s}}}$' %(sp_order, norm_order), rotation=0, fontsize=24)
- ax0.yaxis.set_label_coords(-0.16, 0.9)
-
- ax0.text(0.96, 0.03, args['pdfset']+', %s to %s' %(sp_order, norm_order), alpha=0.6, transform=ax0.transAxes, fontsize=14, verticalalignment='bottom', horizontalalignment='right')
- plt.legend(bbox_to_anchor=(1.02, 1), loc='upper left') #'upper left' refers to the legend box
-
- plt.tight_layout()
- fig0.tight_layout()
-
- #naming of the plots
- if args['filename'] is not None:
- stackplotname = '%s_%s.png' %(args['filename'], sp_order+'_'+norm_order)
- else:
- stackplotname = '%s.stackplot_%s.png' %(tablename, sp_order+'_'+norm_order)
- fig0.savefig(stackplotname, bbox_inches='tight')
-
- print('Stackplot for subprocess comparison saved as: %s' %(stackplotname))
-
-
-
- else:
- print("\n")
- print("----------------------------------------------------------")
- print("Create individual plot for each single subprocess \n")
-
- #individual plot for each subproc
- #create folder for saving the plots, will be overwritten in case it already exists (avoid countless folders while testing)
- current_dir = os.getcwd()
-
- if args['filename'] is not None:
- prename = args['filename']
- else: prename = tablename
-
- final_dir = os.path.join(current_dir, r'%s_sp_%s_%s' %(prename, sp_order, norm_order))
-
- if os.path.exists(final_dir):
- shutil.rmtree(final_dir) #remove final_dir if it already exists
- os.makedirs(final_dir)
-
- #loop where single plot for each of the selected subprocesses is created
- #start the plot
- color = iter(cm.rainbow(np.linspace(0,1,len(args['subproc']))))
-
- for i in range(0, len(args['subproc'])):
- c = next(color)
- plt.close()
- fig0 = plt.figure(figsize=(8, 7))
- ax0 = fig0.add_subplot(111)
-
- #plot the fractions
- y0 = np.zeros([len(bin_bounds)])
-
- ax0.fill_between(bin_bounds.flatten(), steppify_bin(y0), steppify_bin(pos_contr[i]),
- facecolor=c, alpha=0.4)
- ax0.fill_between(bin_bounds.flatten(), steppify_bin(y0), steppify_bin(neg_contr[i]),
- facecolor=c, hatch='X', alpha=0.4)
-
- #get limits for axes:
- xlim = ax0.get_xlim()
- ylim = ax0.get_ylim()
-
- ax0.set_xscale('log', nonposx='clip')
-
- xmin_ = xlim[0]
- xmax_ = xlim[1]
- ymin_ = ylim[0]+0.3*ylim[0]
- ymax_ = ylim[1]+0.1*ylim[1]
-
- if (args['constyaxis'] is not None): #take limits that user has chosen
- if (len(args['constyaxis'])==2):
- ylow = args['constyaxis'][0]
- yhigh = args['constyaxis'][1]
- if ylow < yhigh:
- ax0.axis([0, xmax_, ylow, yhigh])
- else:
- print(" # ERROR! when trying to set y-limits.")
- print("Check choice of y-limits. First value should be lower than second.")
- print("Syntax: -c <lower ylimit> <higher ylimit>.")
- break
- elif (len(args['constyaxis'])==0): #set constant y-axis
- ax0.axis([0, xmax_, -0.20, 1.20])
- else:
- print(" # ERROR! when trying to set y-limits.")
- print("Something went wrong. Check input after -c.")
- print("Possible: no value, two values or don't use -c at all for flexible axis. \n")
- break
- elif (args['constyaxis'] is None): #use flexible y-axis
- ax0.set_xlim(xmin_, xmax_) ##
- ax0.set_ylim(ymin_, ymax_) ##set axis limit according to data
-
- ax0.autoscale(enable=True, axis='x', tight=True)
- plt.axhline(y=1, xmin=0, xmax=1, color='k', linestyle='dotted')
- plt.axhline(y=0, xmin=0, xmax=1, color='k', linestyle='-')
-
-
- titlename = args['subproc'][i]
-
- plt.title(x=0.5, y=1.01, s=titlename, fontsize=20)
- ax0.set_xlabel('%s' %xlabel, fontsize=16, horizontalalignment='right')
- ax0.xaxis.set_label_coords(1.0, -0.04)
- ax0.set_ylabel(r'$\frac{\mathrm{XS_{sp, %s}}}{\mathrm{XS_{tot, %s}}}$' %(sp_order, norm_order), rotation=0, fontsize=24)
- ax0.yaxis.set_label_coords(-0.16, 0.9)
- ax0.text(0.96, 0.03, args['pdfset']+', %s to %s' %(sp_order, norm_order), transform=ax0.transAxes, alpha=0.6, fontsize=14, verticalalignment='bottom', horizontalalignment='right')
- ax0.text(0.96, 0.08, '%s' %table_, transform=ax0.transAxes, alpha=0.6, fontsize=14, verticalalignment='bottom', horizontalalignment='right')
-
- #plt.legend() #location parameter loc='upper right'
-
- plt.tight_layout()
- fig0.tight_layout()
- plotname = '%s_%s_singlefrac_%s_%s.png' %(tablename, args['subproc'][i], sp_order, norm_order)
-
- #save the plots
- fig0.savefig(os.path.join(final_dir,plotname))
- print('Plot saved in %s' %(final_dir))
- print('fraction plot saved as: %s' %(plotname), "\n \n")
-
-
-
- else:
- print("Something went wrong. Check choice of order(s). \n")
-
-
-
-
-
-
-
-""" #other alternative possibility to switch orders on/off
- for order, setting in sorted(b_order.iteritems()): #go through tuples in b_order, evaluate
- print("\n", "order: ", order)
- print("setting: ", setting)
- index=0
- for j in setting:
- print("index", index)
- print("j", j)
- fnlo.SetContributionON(fastnlo.kFixedOrder,index,j);
- index+=1
-"""
-
-
-
-
-
-
-
-
-
-
-
-
-##################################################
-##function to make better uncertainty plot (shaded)
-##could maybe be replaced by .flatten() ??
-def steppify_bin(arr, isx=False):
- """
- Produce stepped array of arr, needed for example for stepped fill_betweens.
- Pass all x bin edges to produce stepped x arr and all y bincontents to produce
- stepped bincontents representation
- steppify_bin([1,2,3], True)
- -> [1,2,2,3]
- steppify_bin([5,6])
- -> [5,5,6,6]
- """
- if isx:
- newarr = np.array(zip(arr[:-1], arr[1:])).ravel()
- else:
- newarr = np.array(zip(arr, arr)).ravel()
- return newarr
-
-
-#################################################
-
-if __name__ == '__main__':
- main()
diff --git a/tools/plotting/yb_ystar/compare_spLO_xstot.py b/tools/plotting/yb_ystar/compare_spLO_xstot.py
deleted file mode 100755
index 391809fb..00000000
--- a/tools/plotting/yb_ystar/compare_spLO_xstot.py
+++ /dev/null
@@ -1,377 +0,0 @@
-#! /usr/bin/env python
-
-####################################################
-#
-#Here, the contributions of subprocesses in LO are
-#compared to the total cross section in LO+NLO.
-#
-####################################################
-
-import argparse
-import numpy as np
-import os
-import time, datetime
-import matplotlib
-import matplotlib.pyplot as plt
-import matplotlib.gridspec as gridspec
-
-import fastnlo
-
-
-
-def main():
-
- parser = argparse.ArgumentParser()
-
- #make it possible to use different options when running
- parser.add_argument('-t','--table', default='table.tab', required=True,
- help='FastNLO tables that shall be evaluated.')
- parser.add_argument('-p', '--pdfset', default='CT14nlo',
- help='PDFset(s) to evaluate fastNLO tables.')
- parser.add_argument('-m', '--member', default=0, type=int,
- help='Member of PDFset, default is 0.')
-
- parser.add_argument('-s', '--subproc', default=['gg'], type=str, nargs='*',
- help='Subprocesses that shall be compared. \n'
- 'Options: gg, gq, qiqi, qiai, qiqj, qiaj')
- parser.add_argument('-c', '--compare', default=False, const=True, type=bool, nargs='?',
- help='Stackplot if option -c is chosen. Otherwise single-process-plots.')
- #help='Stackplot if True. Single-Process-Plot if False.')
-
-
-
- #Dictionary for how the subprocesses are called
- processes = {"gg":"Gluon-Gluon",
- "gq":"Gluon-Quark/Antiquark",
- "qiqi":"Quark-Quark and $\mathrm{\overline{q_i} \; \overline{q_i}}$ (same flavor)",
- "qiai":"Quark-Antiquark (same flavor)",
- "qiqj":"Quark-Quark and $\mathrm{\overline{q_i} \; \overline{q_j}}}$ (different flavors)",
- "qiaj":"Quark-Antiquark (different flavors)"}
-
- #shorter for labeling the stackplot
- labeling = {"gg":r"$\mathrm{gg}$",
- "gq":r"$\mathrm{gq} \ & \ \mathrm{g\overline{q}}$",
- "qiqi":r"$\mathrm{q_i q_i} \ & \ \mathrm{\overline{q_i} \; \overline{q_i}}$",
- "qiai":r"$\mathrm{q_i} \; \mathrm{\overline{q_i}}$",
- "qiqj":r"$\mathrm{q_i q_j} \ & \ \mathrm{\overline{q_i} \; \overline{q_j}}$",
- "qiaj":r"$\mathrm{q_i} \; \mathrm{\overline{q_j}}$"}
-
-
- #parse arguments
- args = vars(parser.parse_args())
-
- #table name
- table_ = os.path.basename(args['table'])
- tablename = os.path.splitext(table_)[0]
- tablename = os.path.splitext(tablename)[0] #because yb tables have .tab.gz ending (getting rid of both here)
- print '\n'
- print 'Table: ', tablename, '\n'
-
- #pdfset name
- pdfset_ = os.path.basename(args['pdfset'])
- pdfname = os.path.splitext(pdfset_)[0]
- print 'PDF Set: ', pdfname, '\n'
-
-
- #checks for investigated phase space region --> mainly for plot-labeling
-
- #range of y_boost
- def yb():
- if 'yb0' in args['table']:
- yb_range = r'$\mathrm{0.0 \leq \/ y_b < 0.5}$'
- info_yb = r'0 <= yb < 0.5'
- elif 'yb1' in args['table']:
- yb_range = r'$\mathrm{0.5 \leq \/ y_b < 1.0}$'
- info_yb = r'0.5 <= yb < 1.0'
- elif 'yb2' in args['table']:
- yb_range = r'$\mathrm{1.0 \leq \/ y_b < 1.5}$'
- info_yb = r'1.0 <= yb < 1.5'
- elif 'yb3' in args['table']:
- yb_range = r'$\mathrm{1.5 \leq \/ y_b < 2.0}$'
- info_yb = r'1.5 <= yb < 2.0'
- elif 'yb4' in args['table']:
- yb_range = r'$\mathrm{2.0 \leq \/ y_b < 2.5}$'
- info_yb = r'2.0 <= yb < 2.5'
- else:
- print 'No info about range of yb.'
- yb_range = ''
- info_yb = ''
- return yb_range, info_yb #yb_range is in correct format for labeling
-
-
- #range of ystar
- def ystar():
- if 'ystar0' in args['table']:
- ystar_range = r'$\mathrm{0.0 \leq \/ y^{\ast} < 0.5}$'
- info_ystar = r'0 <= ystar < 0.5'
- elif 'ystar1' in args['table']:
- ystar_range = r'$\mathrm{0.5 \leq \/ y^{\ast} < 1.0}$'
- info_ystar = r'0.5 <= ystar < 1.0'
- elif 'ystar2' in args['table']:
- ystar_range = r'$\mathrm{1.0 \leq \/ y^{\ast} < 1.5}$'
- info_ystar = r'1.0 <= ystar < 1.5'
- elif 'ystar3' in args['table']:
- ystar_range = r'$\mathrm{1.5 \leq \/ y^{\ast} < 2.0}$'
- info_ystar = r'1.5 <= ystar < 2.0'
- elif 'ystar4' in args['table']:
- ystar_range = r'$\mathrm{2.0 \leq \/ y^{\ast} < 2.5}$'
- info_ystar = r'2.0 <= ystar < 2.5'
- else:
- print 'No info about range of ystar.'
- ystar_range = ''
- info_ystar = ''
- return ystar_range, info_ystar
-
-
-
- #show information about phase space region
- yb_range = yb()[0]
- ystar_range = ystar()[0]
-
- info_yb = yb()[1]
- info_ystar = ystar()[1]
-
- print 'Phase Space Region:'
- print '%s' %info_yb
- print '%s' %info_ystar
-
-
-
- #Which subprocesses are evaluated?
- print 'Subprocesses that are investigated: \n'
- print args['subproc']
-
-
- ######## EVALUATING ########
- fnlo = fastnlo.fastNLOLHAPDF(args['table'], args['pdfset'], args['member'])
- fnlo.SetContributionON(fastnlo.kFixedOrder,0,True);
- fnlo.SetContributionON(fastnlo.kFixedOrder,1,True);
-
-
- #cross section with every process included (for LO+NLO)
- fnlo.CalcCrossSection()
- xs_all = np.array(fnlo.GetCrossSection())
- print 'Cross section with all subprocesses xs_all: \n'
- print xs_all, '\n \n'
-
- #cross section of single processes in LO (only LO-contribution considered, no NLO!)
- fnlo.SetContributionON(fastnlo.kFixedOrder, 0, True); #switch LO on
- fnlo.SetContributionON(fastnlo.kFixedOrder, 1, False); #switch NLO off
- xs_subproc_list = []
- for i in range(0, len(args['subproc'])):
- subproc_name = args['subproc'][i]
- print 'Selected subprocess: %s' %subproc_name
- fnlo.SelectProcesses(subproc_name)
- fnlo.CalcCrossSection()
- xs_subproc = np.array(fnlo.GetCrossSection())
- xs_subproc_list.append(xs_subproc)
- print 'XS for subprocess %s: \n' %args['subproc'][i]
- print xs_subproc, '\n \n'
- xs_sub = np.array(xs_subproc_list) #Array with xs for selected processes in certain pTZ-bins
- print 'XS of the selected processes in different pTZ-regions: \n'
- print xs_sub, '\n \n'
-
- #how much do the selected subprocesses contribute to the total xs?
- fractions = np.divide(xs_sub,xs_all)
- print 'fractions: \n', fractions, '\n'
-
- for i in range(0, len(fractions)):
- print 'fraction of %s in LO on total XS: ' %args['subproc'][i]
- print fractions[i,:], '\n \n'
-
-
- #checking whether contribution is positive or negative
- pos_contr = np.array(fractions)
- neg_contr = np.array(fractions)
- bin_bounds = np.array(fnlo.GetObsBinsBounds(0))
-
- for i in range(0, len(args['subproc'])): #go through 'fractions' line by line = subprocess after subprocess
- for j in range(0, len(bin_bounds)): #check the different pTZ bins
- if (fractions[i,j] < 0):
- pos_contr[i,j] = 0.0
- else:
- neg_contr[i,j] = 0.0
- print 'Positive LO contributions of %s on total XS: \n' %args['subproc'][i]
- print pos_contr[i,:], '\n'
- print 'Negative LO contributions of %s on total XS: \n' %args['subproc'][i]
- print neg_contr[i,:], '\n \n'
-
- print 'Summary of positive and negative LO contributions of the subprocesses: \n'
- print pos_contr, '\n'
- print neg_contr, '\n'
-
-
-
- ######## PLOTTING ########
- #Bins for x-axis (according to first pdf (0))
- bin_bounds = np.array(fnlo.GetObsBinsBounds(0))
- #print bin_bounds.T #bin_bounds.T[0] = lower bounds, bin_bounds.T[1] upper bounds
- #print bin_bounds.flatten()
-
-
- #dictionary for using specific color per process
- process_color = {"gg":"g",
- "gq":"r",
- "qiqi":"y",
- "qiai":"c",
- "qiqj":"m",
- "qiaj":"b"}
-
- #Check whether comparison-plot or single-plot is required:
- print args['compare']
- if (args['compare']==True): #stackplot that compares the contribution of the subprocesses to the total XS
- plt.close()
- fig0 = plt.figure(figsize=(8,7))
- ax0 = fig0.add_subplot(111)
-
- patches = [] #later needed for legend
- #number of bins via bin_bounds variable
- y0_pos = np.zeros([len(bin_bounds)]) #for tracking the current 'height' of stackplot, in order to staple and not overlap the values for each subprocess
- y0_neg = np.zeros([len(bin_bounds)]) #same as above for the negative contributions, plotted below x-axis
-
- for i in range(0, len(args['subproc'])): #go through fractions, subproc by subproc
-
- #positive contributions
- ax0.fill_between(bin_bounds.flatten(), steppify_bin(y0_pos), steppify_bin(y0_pos+pos_contr[i]), #label=labeling[args['subproc'][i]],
- facecolor=process_color[args['subproc'][i]], alpha=0.4)
- y0_pos += pos_contr[i] #calculate new height of stackplot
-
- #negative contributions
- ax0.fill_between(bin_bounds.flatten(), steppify_bin(y0_neg), steppify_bin(y0_neg+neg_contr[i]), #label=labeling[args['subproc'][i]],
- facecolor=process_color[args['subproc'][i]], hatch='X', alpha=0.4)
- y0_neg += neg_contr[i] #new 'height' below x-axis
-
- #patch for subprocess i (for the legend) ##labeling happens here
- patches.append(matplotlib.patches.Rectangle((0, 0), 0, 0, color=process_color[args['subproc'][i]],
- label=labeling[args['subproc'][i]], alpha=0.4))
- ax0.add_patch(patches[i])
-
- #settings for the whole stackplot
- ax0.set_xscale('log', nonposx='clip')
- ax0.axis([30, 1000, -0.20, 1.20])
- plt.axhline(y=1, xmin=0, xmax=1, color='k', linestyle='dotted') #dotted line at xs_sub/xs=1=100%
- plt.axhline(y=0, xmin=0, xmax=1, color='k', linestyle='-') #line at xs_sub/xs_tot=0
-
- #use phasespace region as title
- title_phasespace = yb_range + ", " + ystar_range
- plt.title(x=0.5, y=1.01, s='%s' %(title_phasespace), fontsize=22)
-
- ax0.set_xlabel('$\mathrm{p_{T,Z}} \ \mathrm{[GeV]}$', fontsize=20)
- ax0.xaxis.set_label_coords(0.9, -0.05)
- ax0.set_ylabel(r'$\frac{\mathrm{XS_{sp,LO}}}{\mathrm{XS_{tot,LO+NLO}}}$', rotation=0, fontsize=22)
- ax0.yaxis.set_label_coords(-0.16, 0.9)
- ax0.text(0.96, 0.03, args['pdfset'], transform=ax0.transAxes, fontsize=14, verticalalignment='bottom', horizontalalignment='right')
- plt.legend()
-
- plt.tight_layout()
- fig0.tight_layout()
-
- #make string with the subprocess names
- processnames = ""
- for i in range(0, len(args['subproc'])):
- processnames += ("_"+args['subproc'][i])
- stackplotname = 'compare_spLO_xstot%s_%s_%s.png' %(processnames, tablename, args['pdfset'])
- fig0.savefig(stackplotname)
-
- print 'Stackplot for subprocess comparison saved as: %s' %(stackplotname)
-
-
-
-
- else: #individual plot for each subproc
- #create folder for saving the plots
- current_dir = os.getcwd()
-
- final_dir = os.path.join(current_dir, r'single_spLO_plots_%s' %tablename)
- now = datetime.datetime.now()
- nowstr = 'single_spLO_xstot_plots_%s_' %(tablename)
- nowstr += now.strftime('%Y-%m-%d_%H-%M-%S')
- final_dir_date = os.path.join(current_dir, nowstr)
-
- #check if directory already exists, if so: include datetime to name
- if not os.path.exists(final_dir):
- os.makedirs(final_dir)
- else:
- os.makedirs(final_dir_date)
-
- #loop where single plot for each of the selected subprocesses is created
- #start the plot
- for i in range(0, len(args['subproc'])):
- plt.close()
- fig0 = plt.figure(figsize=(8, 7))
- ax0 = fig0.add_subplot(111)
-
- #plot the fractions
- y0 = np.zeros([len(bin_bounds)])
- ax0.fill_between(bin_bounds.flatten(), steppify_bin(y0), steppify_bin(fractions[i]),
- facecolor=process_color[args['subproc'][i]], alpha=0.4) #improve labeling?
-
-
-
- #info about phase space region included as text (see below)
- ax0.set_xscale('log', nonposx='clip')
- ax0.axis([30, 1000, -0.4, 1.2])
- plt.axhline(y=1, xmin=0, xmax=1, color='k', linestyle='dotted')
- plt.axhline(y=0, xmin=0, xmax=1, color='k', linestyle='-')
-
- titlename = processes[args['subproc'][i]]
-
- plt.title(x=0.5, y=1.01, s=titlename, fontsize=20)
- ax0.set_xlabel('$\mathrm{p_{T,Z}} \ \mathrm{[GeV]}$', fontsize=20)
- ax0.xaxis.set_label_coords(0.9, -0.05)
- ax0.set_ylabel(r'$\frac{\mathrm{XS_{sp,LO}}}{\mathrm{XS_{tot,LO+NLO}}}$', rotation=0, fontsize=22)
- ax0.yaxis.set_label_coords(-0.16, 0.9)
- ax0.text(0.96, 0.03, args['pdfset'], transform=ax0.transAxes, fontsize=14, verticalalignment='bottom', horizontalalignment='right')
- ax0.text(0.96, 0.94, yb_range, transform=ax0.transAxes, fontsize=16, verticalalignment='bottom', horizontalalignment='right')
- ax0.text(0.96, 0.88, ystar_range, transform=ax0.transAxes, fontsize=16, verticalalignment='bottom', horizontalalignment='right')
- plt.legend() #location parameter loc='upper right'
-
- plt.tight_layout()
- fig0.tight_layout()
- plotname = '%s_LOfrac_xstot_%s_%s.png' %(args['subproc'][i], tablename, args['pdfset'])
-
-
- #save the plots
- if os.path.exists(final_dir_date):
- fig0.savefig(os.path.join(final_dir_date,plotname))
- print 'Plot saved in %s' %(final_dir_date)
- else:
- fig0.savefig(os.path.join(final_dir,plotname))
- print 'Plot saved in %s' %(final_dir)
-
- print '\n'
- print 'fraction plot saved as: %s' %(plotname)
-
-
-
-
-
-##################################################
-##function to make better uncertainty plot (shaded)
-##could maybe be replaced by .flatten() ??
-def steppify_bin(arr, isx=False):
- """
- Produce stepped array of arr, needed for example for stepped fill_betweens.
- Pass all x bin edges to produce stepped x arr and all y bincontents to produce
- stepped bincontents representation
- steppify_bin([1,2,3], True)
- -> [1,2,2,3]
- steppify_bin([5,6])
- -> [5,5,6,6]
- """
- if isx:
- newarr = np.array(zip(arr[:-1], arr[1:])).ravel()
- else:
- newarr = np.array(zip(arr, arr)).ravel()
- return newarr
-
-
-#################################################
-
-if __name__ == '__main__':
- main()
-
-
-
-
-
diff --git a/tools/plotting/yb_ystar/compare_spNLO_xstot.py b/tools/plotting/yb_ystar/compare_spNLO_xstot.py
deleted file mode 100755
index 57ac4668..00000000
--- a/tools/plotting/yb_ystar/compare_spNLO_xstot.py
+++ /dev/null
@@ -1,377 +0,0 @@
-#! /usr/bin/env python
-
-##################################
-#Here, the contribution of subprocesses in NLO are
-#compared to the total cross section in LO+NLO.
-#Six available subprocess categories.
-##################################
-
-import argparse
-import numpy as np
-import os
-import time, datetime
-import matplotlib
-import matplotlib.pyplot as plt
-import matplotlib.gridspec as gridspec
-
-import fastnlo
-
-
-
-def main():
-
- parser = argparse.ArgumentParser()
-
- #make it possible to use different options when running
- parser.add_argument('-t','--table', default='table.tab', required=True,
- help='FastNLO tables that shall be evaluated.')
- parser.add_argument('-p', '--pdfset', default='CT14nlo',
- help='PDFset(s) to evaluate fastNLO tables.')
- parser.add_argument('-m', '--member', default=0, type=int,
- help='Member of PDFset, default is 0.')
-
- parser.add_argument('-s', '--subproc', default=['gg'], type=str, nargs='*',
- help='Subprocesses that shall be compared. \n'
- 'Options: gg, gq, qiqi, qiai, qiqj, qiaj')
- parser.add_argument('-c', '--compare', default=False, const=True, type=bool, nargs='?',
- help='Stackplot if option -c is chosen. Otherwise single-process-plots.')
- #help='Stackplot if True. Single-Process-Plot if False.')
-
-
-
- #Dictionary for how the subprocesses are called
- processes = {"gg":"Gluon-Gluon",
- "gq":"Gluon-Quark/Antiquark",
- "qiqi":"Quark-Quark and $\mathrm{\overline{q_i} \; \overline{q_i}}$ (same flavor)",
- "qiai":"Quark-Antiquark (same flavor)",
- "qiqj":"Quark-Quark and $\mathrm{\overline{q_i} \; \overline{q_j}}}$ (different flavors)",
- "qiaj":"Quark-Antiquark (different flavors)"}
-
- #shorter for labeling the stackplot
- labeling = {"gg":r"$\mathrm{gg}$",
- "gq":r"$\mathrm{gq} \ & \ \mathrm{g\overline{q}}$",
- "qiqi":r"$\mathrm{q_i q_i} \ & \ \mathrm{\overline{q_i} \; \overline{q_i}}$",
- "qiai":r"$\mathrm{q_i} \; \mathrm{\overline{q_i}}$",
- "qiqj":r"$\mathrm{q_i q_j} \ & \ \mathrm{\overline{q_i} \; \overline{q_j}}$",
- "qiaj":r"$\mathrm{q_i} \; \mathrm{\overline{q_j}}$"}
-
-
- #parse arguments
- args = vars(parser.parse_args())
-
- #table name
- table_ = os.path.basename(args['table'])
- tablename = os.path.splitext(table_)[0]
- tablename = os.path.splitext(tablename)[0] #because yb tables have .tab.gz ending (getting rid of both here)
- print '\n'
- print 'Table: ', tablename, '\n'
-
- #pdfset name
- pdfset_ = os.path.basename(args['pdfset'])
- pdfname = os.path.splitext(pdfset_)[0]
- print 'PDF Set: ', pdfname, '\n'
-
-
- #checks for investigated phase space region --> mainly for plot-labeling
-
- #range of y_boost
- def yb():
- if 'yb0' in args['table']:
- yb_range = r'$\mathrm{0.0 \leq \/ y_b < 0.5}$'
- info_yb = r'0 <= yb < 0.5'
- elif 'yb1' in args['table']:
- yb_range = r'$\mathrm{0.5 \leq \/ y_b < 1.0}$'
- info_yb = r'0.5 <= yb < 1.0'
- elif 'yb2' in args['table']:
- yb_range = r'$\mathrm{1.0 \leq \/ y_b < 1.5}$'
- info_yb = r'1.0 <= yb < 1.5'
- elif 'yb3' in args['table']:
- yb_range = r'$\mathrm{1.5 \leq \/ y_b < 2.0}$'
- info_yb = r'1.5 <= yb < 2.0'
- elif 'yb4' in args['table']:
- yb_range = r'$\mathrm{2.0 \leq \/ y_b < 2.5}$'
- info_yb = r'2.0 <= yb < 2.5'
- else:
- print 'No info about range of yb.'
- yb_range = ''
- info_yb = ''
- return yb_range, info_yb #yb_range is in correct format for labeling
-
-
- #range of ystar
- def ystar():
- if 'ystar0' in args['table']:
- ystar_range = r'$\mathrm{0.0 \leq \/ y^{\ast} < 0.5}$'
- info_ystar = r'0 <= ystar < 0.5'
- elif 'ystar1' in args['table']:
- ystar_range = r'$\mathrm{0.5 \leq \/ y^{\ast} < 1.0}$'
- info_ystar = r'0.5 <= ystar < 1.0'
- elif 'ystar2' in args['table']:
- ystar_range = r'$\mathrm{1.0 \leq \/ y^{\ast} < 1.5}$'
- info_ystar = r'1.0 <= ystar < 1.5'
- elif 'ystar3' in args['table']:
- ystar_range = r'$\mathrm{1.5 \leq \/ y^{\ast} < 2.0}$'
- info_ystar = r'1.5 <= ystar < 2.0'
- elif 'ystar4' in args['table']:
- ystar_range = r'$\mathrm{2.0 \leq \/ y^{\ast} < 2.5}$'
- info_ystar = r'2.0 <= ystar < 2.5'
- else:
- print 'No info about range of ystar.'
- ystar_range = ''
- info_ystar = ''
- return ystar_range, info_ystar
-
-
-
- #show information about phase space region
- yb_range = yb()[0]
- ystar_range = ystar()[0]
-
- info_yb = yb()[1]
- info_ystar = ystar()[1]
-
- print 'Phase Space Region:'
- print '%s' %info_yb
- print '%s' %info_ystar
-
-
-
- #Which subprocesses are evaluated?
- print 'Subprocesses that are investigated: \n'
- print args['subproc']
-
-
- ######## EVALUATING ########
- fnlo = fastnlo.fastNLOLHAPDF(args['table'], args['pdfset'], args['member'])
- fnlo.SetContributionON(fastnlo.kFixedOrder,0,True);
- fnlo.SetContributionON(fastnlo.kFixedOrder,1,True);
-
-
- #cross section with every process included (for LO+NLO)
- fnlo.CalcCrossSection()
- xs_all = np.array(fnlo.GetCrossSection())
- print 'Cross section with all subprocesses xs_all: \n'
- print xs_all, '\n \n'
-
- #cross section of single processes in NLO (no LO-contribution considered!)
- fnlo.SetContributionON(fastnlo.kFixedOrder, 1, True); #switch NLO on
- fnlo.SetContributionON(fastnlo.kFixedOrder, 0, False); #switch LO off
- xs_subproc_list = []
- for i in range(0, len(args['subproc'])):
- subproc_name = args['subproc'][i]
- print 'Selected subprocess: %s' %subproc_name
- fnlo.SelectProcesses(subproc_name)
- fnlo.CalcCrossSection()
- xs_subproc = np.array(fnlo.GetCrossSection())
- xs_subproc_list.append(xs_subproc)
- print 'XS for subprocess %s: \n' %args['subproc'][i]
- print xs_subproc, '\n \n'
- xs_sub = np.array(xs_subproc_list) #Array with xs for selected processes in certain pTZ-bins
- print 'XS of the selected processes in different pTZ-regions: \n'
- print xs_sub, '\n \n'
-
- #how much do the selected subprocesses contribute to the total xs?
- fractions = np.divide(xs_sub,xs_all)
- print 'fractions: \n', fractions, '\n'
-
- for i in range(0, len(fractions)):
- print 'fraction of %s in NLO on total XS: ' %args['subproc'][i]
- print fractions[i,:], '\n \n'
-
-
- #checking whether contribution is positive or negative
- pos_contr = np.array(fractions)
- neg_contr = np.array(fractions)
- bin_bounds = np.array(fnlo.GetObsBinsBounds(0))
-
- for i in range(0, len(args['subproc'])): #go through 'fractions' line by line = subprocess after subprocess
- for j in range(0, len(bin_bounds)): #check the different pTZ bins
- if (fractions[i,j] < 0):
- pos_contr[i,j] = 0.0
- else:
- neg_contr[i,j] = 0.0
- print 'Positive NLO contributions of %s on total XS: \n' %args['subproc'][i]
- print pos_contr[i,:], '\n'
- print 'Negative NLO contributions of %s on total XS: \n' %args['subproc'][i]
- print neg_contr[i,:], '\n \n'
-
- print 'Summary of positive and negative NLO contributions of the subprocesses: \n'
- print pos_contr, '\n'
- print neg_contr, '\n'
-
-
-
- ######## PLOTTING ########
- #Bins for x-axis (according to first pdf (0))
- bin_bounds = np.array(fnlo.GetObsBinsBounds(0))
- #print bin_bounds.T #bin_bounds.T[0] = lower bounds, bin_bounds.T[1] upper bounds
- #print bin_bounds.flatten()
-
-
- #dictionary for using specific color per process
- process_color = {"gg":"g",
- "gq":"r",
- "qiqi":"y",
- "qiai":"c",
- "qiqj":"m",
- "qiaj":"b"}
-
- #Check whether comparison-plot or single-plot is required:
- print args['compare']
- if (args['compare']==True): #stackplot that compares the contribution of the subprocesses to the total XS
- plt.close()
- fig0 = plt.figure(figsize=(8,7))
- ax0 = fig0.add_subplot(111)
-
- patches = [] #later needed for legend
- #number of bins via bin_bounds variable
- y0_pos = np.zeros([len(bin_bounds)]) #for tracking the current 'height' of stackplot, in order to staple and not overlap the values for each subprocess
- y0_neg = np.zeros([len(bin_bounds)]) #same as above for the negative contributions, plotted below x-axis
-
- for i in range(0, len(args['subproc'])): #go through fractions, subproc by subproc
-
- #positive contributions
- ax0.fill_between(bin_bounds.flatten(), steppify_bin(y0_pos), steppify_bin(y0_pos+pos_contr[i]), #label=labeling[args['subproc'][i]],
- facecolor=process_color[args['subproc'][i]], alpha=0.4)
- y0_pos += pos_contr[i] #calculate new height of stackplot
-
- #negative contributions
- ax0.fill_between(bin_bounds.flatten(), steppify_bin(y0_neg), steppify_bin(y0_neg+neg_contr[i]), #label=labeling[args['subproc'][i]],
- facecolor=process_color[args['subproc'][i]], hatch='X', alpha=0.4)
- y0_neg += neg_contr[i] #new 'height' below x-axis
-
- #patch for subprocess i (for the legend) ##labeling happens here
- patches.append(matplotlib.patches.Rectangle((0, 0), 0, 0, color=process_color[args['subproc'][i]],
- label=labeling[args['subproc'][i]], alpha=0.4))
- ax0.add_patch(patches[i])
-
- #settings for the whole stackplot
- ax0.set_xscale('log', nonposx='clip')
- ax0.axis([30, 1000, -0.20, 1.20])
- plt.axhline(y=1, xmin=0, xmax=1, color='k', linestyle='dotted') #dotted line at xs_sub/xs=1=100%
- plt.axhline(y=0, xmin=0, xmax=1, color='k', linestyle='-') #line at xs_sub/xs_tot=0
-
- #use phasespace region as title
- title_phasespace = yb_range + ", " + ystar_range
- plt.title(x=0.5, y=1.01, s='%s' %(title_phasespace), fontsize=22)
-
- ax0.set_xlabel('$\mathrm{p_{T,Z}} \ \mathrm{[GeV]}$', fontsize=20)
- ax0.xaxis.set_label_coords(0.9, -0.05)
- ax0.set_ylabel(r'$\frac{\mathrm{XS_{sp,NLO}}}{\mathrm{XS_{tot,LO+NLO}}}$', rotation=0, fontsize=22)
- ax0.yaxis.set_label_coords(-0.16, 0.9)
- ax0.text(0.96, 0.03, args['pdfset'], transform=ax0.transAxes, fontsize=14, verticalalignment='bottom', horizontalalignment='right')
- plt.legend()
-
- plt.tight_layout()
- fig0.tight_layout()
-
- #make string with the subprocess names
- processnames = ""
- for i in range(0, len(args['subproc'])):
- processnames += ("_"+args['subproc'][i])
- stackplotname = 'compare_spNLO_xstot%s_%s_%s.png' %(processnames, tablename, args['pdfset'])
- fig0.savefig(stackplotname)
-
- print 'Stackplot for subprocess comparison saved as: %s' %(stackplotname)
-
-
-
-
- else: #individual plot for each subproc
- #create folder for saving the plots
- current_dir = os.getcwd()
-
- final_dir = os.path.join(current_dir, r'single_spNLO_plots_%s' %tablename)
- now = datetime.datetime.now()
- nowstr = 'single_spNLO_xstot_plots_%s_' %(tablename)
- nowstr += now.strftime('%Y-%m-%d_%H-%M-%S')
- final_dir_date = os.path.join(current_dir, nowstr)
-
- #check if directory already exists, if so: include datetime to name
- if not os.path.exists(final_dir):
- os.makedirs(final_dir)
- else:
- os.makedirs(final_dir_date)
-
- #loop where single plot for each of the selected subprocesses is created
- #start the plot
- for i in range(0, len(args['subproc'])):
- plt.close()
- fig0 = plt.figure(figsize=(8, 7))
- ax0 = fig0.add_subplot(111)
-
- #plot the fractions
- y0 = np.zeros([len(bin_bounds)])
- ax0.fill_between(bin_bounds.flatten(), steppify_bin(y0), steppify_bin(fractions[i]),
- facecolor=process_color[args['subproc'][i]], alpha=0.4) #improve labeling?
-
-
-
- #info about phase space region included as text (see below)
- ax0.set_xscale('log', nonposx='clip')
- ax0.axis([30, 1000, -0.4, 1.2])
- plt.axhline(y=1, xmin=0, xmax=1, color='k', linestyle='dotted')
- plt.axhline(y=0, xmin=0, xmax=1, color='k', linestyle='-')
-
- titlename = processes[args['subproc'][i]]
-
- plt.title(x=0.5, y=1.01, s=titlename, fontsize=20)
- ax0.set_xlabel('$\mathrm{p_{T,Z}} \ \mathrm{[GeV]}$', fontsize=20)
- ax0.xaxis.set_label_coords(0.9, -0.05)
- ax0.set_ylabel(r'$\frac{\mathrm{XS_{sp,NLO}}}{\mathrm{XS_{tot,LO+NLO}}}$', rotation=0, fontsize=22)
- ax0.yaxis.set_label_coords(-0.16, 0.9)
- #ax0.set_ylabel('$\mathrm{XS_{sp,NLO}} / \mathrm{XS_{tot,LO+NLO}}$', fontsize=20)
- ax0.text(0.96, 0.03, args['pdfset'], transform=ax0.transAxes, fontsize=14, verticalalignment='bottom', horizontalalignment='right')
- ax0.text(0.96, 0.94, yb_range, transform=ax0.transAxes, fontsize=16, verticalalignment='bottom', horizontalalignment='right')
- ax0.text(0.96, 0.88, ystar_range, transform=ax0.transAxes, fontsize=16, verticalalignment='bottom', horizontalalignment='right')
- plt.legend() #location parameter loc='upper right'
-
- plt.tight_layout()
- fig0.tight_layout()
- plotname = '%s_NLOfrac_xstot_%s_%s.png' %(args['subproc'][i], tablename, args['pdfset'])
-
-
- #save the plots
- if os.path.exists(final_dir_date):
- fig0.savefig(os.path.join(final_dir_date,plotname))
- print 'Plot saved in %s' %(final_dir_date)
- else:
- fig0.savefig(os.path.join(final_dir,plotname))
- print 'Plot saved in %s' %(final_dir)
-
- print '\n'
- print 'fraction plot saved as: %s' %(plotname)
-
-
-
-
-
-##################################################
-##function to make better uncertainty plot (shaded)
-##could maybe be replaced by .flatten() ??
-def steppify_bin(arr, isx=False):
- """
- Produce stepped array of arr, needed for example for stepped fill_betweens.
- Pass all x bin edges to produce stepped x arr and all y bincontents to produce
- stepped bincontents representation
- steppify_bin([1,2,3], True)
- -> [1,2,2,3]
- steppify_bin([5,6])
- -> [5,5,6,6]
- """
- if isx:
- newarr = np.array(zip(arr[:-1], arr[1:])).ravel()
- else:
- newarr = np.array(zip(arr, arr)).ravel()
- return newarr
-
-
-#################################################
-
-if __name__ == '__main__':
- main()
-
-
-
-
-
diff --git a/tools/plotting/yb_ystar/kfactor_ybystar.py b/tools/plotting/yb_ystar/kfactor_ybystar.py
deleted file mode 100755
index 3d04b301..00000000
--- a/tools/plotting/yb_ystar/kfactor_ybystar.py
+++ /dev/null
@@ -1,294 +0,0 @@
-#! /usr/bin/env python
-
-######################################################################
-#Looking at the k-factor
-#(total cross section in NLO to total xs in LO)
-#
-#Comparison of multiple tables that have been chosen.
-#No distinction between different subprocesses (always include ALL).
-######################################################################
-
-import argparse
-import numpy as np
-import os
-import time, datetime
-import matplotlib
-import matplotlib.pyplot as plt
-import matplotlib.gridspec as gridspec
-from matplotlib.pyplot import cm
-
-import fastnlo
-
-
-def main():
-
- parser = argparse.ArgumentParser()
-
- #make it possible to use different options when running
- parser.add_argument('-t','--table', default=['table.tab'], required=True, type=str, nargs='*',
- help='FastNLO tables that shall be evaluated.')
- parser.add_argument('-p', '--pdfset', default='CT14nlo',
- help='PDFset(s) to evaluate fastNLO tables.')
- parser.add_argument('-m', '--member', default=0, type=int,
- help='Member of PDFset, default is 0.')
-
-
-
- #parse arguments
- args = vars(parser.parse_args())
-
- #table name
- table_names = []
- for i in range(0, len(args['table'])):
- table_ = os.path.basename(args['table'][i])
- tablename = os.path.splitext(table_)[0]
- tablename = os.path.splitext(tablename)[0] #because yb tables have .tab.gz ending (getting rid of both here)
- table_names.append(table_)
- tables = np.array(table_names)
- print '\n'
- print 'Tables: ', tables, '\n'
-
- #pdfset name
- pdfset_ = os.path.basename(args['pdfset'])
- pdfname = os.path.splitext(pdfset_)[0]
- print 'PDF Set: ', pdfname, '\n'
-
-
- #checks for investigated phase space region --> mainly for plot-labeling
- #range of y_boost
- def yb(table):
- if 'yb0' in table:
- yb_range = r'$\mathrm{0.0 \leq \/ y_b < 0.5}$'
- info_yb = r'0 <= yb < 0.5'
- elif 'yb1' in table:
- yb_range = r'$\mathrm{0.5 \leq \/ y_b < 1.0}$'
- info_yb = r'0.5 <= yb < 1.0'
- elif 'yb2' in table:
- yb_range = r'$\mathrm{1.0 \leq \/ y_b < 1.5}$'
- info_yb = r'1.0 <= yb < 1.5'
- elif 'yb3' in table:
- yb_range = r'$\mathrm{1.5 \leq \/ y_b < 2.0}$'
- info_yb = r'1.5 <= yb < 2.0'
- elif 'yb4' in table:
- yb_range = r'$\mathrm{2.0 \leq \/ y_b < 2.5}$'
- info_yb = r'2.0 <= yb < 2.5'
- else:
- print 'No info about range of yb.'
- yb_range = ''
- info_yb = ''
- return yb_range, info_yb #yb_range is in correct format for labeling
-
-
- #range of ystar
- def ystar(table):
- if 'ystar0' in table:
- ystar_range = r'$\mathrm{0.0 \leq \/ y^{\ast} < 0.5}$'
- info_ystar = r'0 <= ystar < 0.5'
- elif 'ystar1' in table:
- ystar_range = r'$\mathrm{0.5 \leq \/ y^{\ast} < 1.0}$'
- info_ystar = r'0.5 <= ystar < 1.0'
- elif 'ystar2' in table:
- ystar_range = r'$\mathrm{1.0 \leq \/ y^{\ast} < 1.5}$'
- info_ystar = r'1.0 <= ystar < 1.5'
- elif 'ystar3' in table:
- ystar_range = r'$\mathrm{1.5 \leq \/ y^{\ast} < 2.0}$'
- info_ystar = r'1.5 <= ystar < 2.0'
- elif 'ystar4' in table:
- ystar_range = r'$\mathrm{2.0 \leq \/ y^{\ast} < 2.5}$'
- info_ystar = r'2.0 <= ystar < 2.5'
- else:
- print 'No info about range of ystar.'
- ystar_range = ''
- info_ystar = ''
- return ystar_range, info_ystar
-
-
-
-
- ######## EVALUATING ########
-
- fractions_list = []
- binbounds_list = [[] for i in range(len(args['table']))]
- #binbounds_array = np.array(
- for i in range(0, len(args['table'])):
- tab = tables[i]
- print 'Selected table: %s' %tab
- #show information about phase space region
- yb_range = yb(tab)[0]
- ystar_range = ystar(tab)[0]
- info_yb = yb(tab)[1]
- info_ystar = ystar(tab)[1]
-
- print 'Phase Space Region:'
- print '%s' %info_yb
- print '%s' %info_ystar, '\n'
-
- #Evaluate selected table:
- fnlo = fastnlo.fastNLOLHAPDF(args['table'][i], args['pdfset'], args['member'])
-
- ####leading order cross section##################################################
- fnlo.SetContributionON(fastnlo.kFixedOrder,0,True); #Switch LO on
- fnlo.SetContributionON(fastnlo.kFixedOrder,1,False); #Switch NLO off
-
- #cross section with every process included (only Leading Order)
- fnlo.CalcCrossSection()
- xs_lo = np.array(fnlo.GetCrossSection())
- print '\n'
- print 'Cross section in Leading Order with all subprocesses xs_lo: \n'
- print xs_lo, '\n \n'
-
- ###next-to-leading order cross section (LO+NLO):#################################
- fnlo.SetContributionON(fastnlo.kFixedOrder, 0, True); #switch LO on!
- fnlo.SetContributionON(fastnlo.kFixedOrder, 1, True); #switch NLO on
-
- #cross section with every process included (Next-to-Leading Order)
- fnlo.CalcCrossSection()
- xs_nlo = np.array(fnlo.GetCrossSection())
- print '\n'
- print 'Cross section in Next-to-Leading Order with all subprocesses xs_nlo: \n'
- print xs_nlo, '\n \n'
-
- #XS in NLO compared to XS in LO
- #alternative could be: first calculate all the nlo and lo XS for each table, then divide
- #here it is the other way round: first k-factors for each table, then put them into one list
- fraction = np.divide(xs_nlo,xs_lo) #nlo/lo for one table (in pTZ bins)
- fractions_list.append(fraction) #collect k-factor for all the tables in this list
- obs_binbounds = np.array(fnlo.GetObsBinsBounds(0)) #must be done for each table, as binbounds may differ
- #obs_binbounds = fnlo.GetObsBinBounds(0) #leave it as a list
- binbounds_list[i]= obs_binbounds
- fractions = np.array(fractions_list)
- binbounds = np.array(binbounds_list)
- print 'k-factors for the selected tables: \n'
- print 'fractions: \n', fractions, '\n'
-
- #print 'binbounds:', binbounds
-
- #NOTE: fractions and binbounds sometimes are arrays of lists that differ in length
-
-
-
- ######## PLOTTING ########
- #taking care of the directory
- current_dir = os.getcwd()
-
- final_dir = os.path.join(current_dir, r'kfactor_plots_%s' %pdfname)
- now = datetime.datetime.now()
- nowstr = 'kfactor_plots_%s' %(pdfname)
- nowstr += now.strftime('_%Y-%m-%d_%H-%M-%S')
- final_dir_date = os.path.join(current_dir, nowstr)
-
- #check if final directory already exists, if so: include datetime to name
- if not os.path.exists(final_dir):
- os.makedirs(final_dir)
- else:
- os.makedirs(final_dir_date)
-
-
- plt.close()
- fig0 = plt.figure(figsize=(8,7))
- ax0 = fig0.add_subplot(111)
-
-
- color = iter(cm.rainbow(np.linspace(0,1,len(args['table']))))
- #print 'color:', color
-
-
- ymin = -0.20
- ymax = 1.20
- xmin = 30
- xmax = 1000
- for i in range(0, len(args['table'])): #go through fractions, line by line = table by table
- bin_bounds = binbounds[i][:]
-# print 'x',bin_bounds
- fraction = fractions[i][:]
-# print 'y',fraction
-
- #plot the k-factors in each pTZ bin for table i:
- c = next(color) #use the iterator for choosing the linecolor fractions[i,:]
- ax0.plot(bin_bounds.flatten(), steppify_bin(fraction),
- label=yb(args['table'][i])[0] + ", " + ystar(args['table'][i])[0],
- color=c, alpha=1.0)
-
- #check current limits of y-axis:
- ylim = ax0.get_ylim()
-# if (ylim[0]-0.05 < ymin):
-# ymin = ylim[0]-0.05
-# if (ylim[1]+0.05 > ymax):
-# ymax = ylim[1]+0.05
-
- #check current limits of x-axis:
- xlim = ax0.get_xlim()
-# if (xlim[0]-0.05 < xmin):
-# xmin = ylim[0]-0.05
-# if (xlim[1]+0.05 > xmax):
-# xmax = xlim[1]+0.05
-
-
- #settings for the whole plot
- ax0.set_xscale('log', nonposx='clip')
-# ax0.axis([xmin, xmax, ymin, ymax])
- ax0.axis([xlim[0], xlim[1], ylim[0]-0.05, ylim[1]+0.05]) #flexible axes
- plt.axhline(y=1, xmin=0, xmax=1, color='k', linestyle='dotted')
- plt.axhline(y=0, xmin=0, xmax=1, color='k', linestyle='dotted')
-
-
- #plt.title(x=0.5, y=1.01, s='k-factor', fontsize=22)
- ax0.set_xlabel('$\mathrm{p_{T,Z}} \ \mathrm{[GeV]}$', fontsize=18)
- ax0.xaxis.set_label_coords(0.9, -0.05)
- ax0.set_ylabel(r'$\frac{\mathrm{XS_{tot,NLO}}}{\mathrm{XS_{tot,LO}}}$', rotation=0, fontsize=24)
- ax0.yaxis.set_label_coords(-0.14, 0.9)
-
- ax0.text(0.96, 0.03, args['pdfset'], transform=ax0.transAxes, fontsize=14, verticalalignment='bottom', horizontalalignment='right')
- #plt.legend(loc="best", fancybox=True, framealpha=0.5)
- plt.legend(bbox_to_anchor=(1.02,1), loc='upper left') #'upper left' refers to the legend box
-
- plt.tight_layout()
- fig0.tight_layout()
-
- figname = 'kfactors_%s.png' %(pdfname)
- #save the plot
- if os.path.exists(final_dir_date):
- fig0.savefig(os.path.join(final_dir_date,figname), bbox_inches='tight')
- print 'Plot saved in %s' %(final_dir_date)
- else:
- fig0.savefig(os.path.join(final_dir,figname), bbox_inches='tight')
- print 'Plot saved in %s' %(final_dir)
- print '\n'
-
- print 'k-factor plot saved as: %s' %(figname)
-
-
-
-
-
-
-
-##################################################
-##function to make better uncertainty plot (shaded)
-##could maybe be replaced by .flatten() ??
-def steppify_bin(arr, isx=False):
- """
- Produce stepped array of arr, needed for example for stepped fill_betweens.
- Pass all x bin edges to produce stepped x arr and all y bincontents to produce
- stepped bincontents representation
- steppify_bin([1,2,3], True)
- -> [1,2,2,3]
- steppify_bin([5,6])
- -> [5,5,6,6]
- """
- if isx:
- newarr = np.array(zip(arr[:-1], arr[1:])).ravel()
- else:
- newarr = np.array(zip(arr, arr)).ravel()
- return newarr
-
-
-#################################################
-
-if __name__ == '__main__':
- main()
-
-
-
-
-
diff --git a/tools/plotting/yb_ystar/subproc_contrib_ybystar.py b/tools/plotting/yb_ystar/subproc_contrib_ybystar.py
deleted file mode 100755
index 7b5e33ac..00000000
--- a/tools/plotting/yb_ystar/subproc_contrib_ybystar.py
+++ /dev/null
@@ -1,381 +0,0 @@
-#! /usr/bin/env python
-
-###############################################################
-#Script for comparing the contribution of several subprocesses
-#to the total cross section.
-#Either as stackplot ('-c') or (per default) for single process.
-#
-#Customized for yb-ystar-tables.
-#Six available subprocess categories.
-###############################################################
-
-
-import argparse
-import numpy as np
-import os
-import time, datetime
-import matplotlib
-import matplotlib.pyplot as plt
-import matplotlib.gridspec as gridspec
-
-import fastnlo
-
-
-
-def main():
-
- parser = argparse.ArgumentParser()
-
- #make it possible to use different options when running
- parser.add_argument('-t','--table', default='table.tab', required=True,
- help='FastNLO tables that shall be evaluated.')
- parser.add_argument('-p', '--pdfset', default='CT14nlo',
- help='PDFset(s) to evaluate fastNLO tables.')
- parser.add_argument('-m', '--member', default=0, type=int,
- help='Member of PDFset, default is 0.')
-
- parser.add_argument('-s', '--subproc', default=['gg'], type=str, nargs='*',
- help='Subprocesses that shall be compared. \n'
- 'Options: gg, gq, qiqi, qiai, qiqj, qiaj')
- parser.add_argument('-c', '--compare', default=False, const=True, type=bool, nargs='?',
- help='Stackplot if option -c is chosen. Otherwise single-process-plots.')
- #help='Stackplot if True. Single-Process-Plot if False.')
-
-
-
-
- #Dictionary for how the subprocesses are called
- processes = {"gg":"Gluon-Gluon",
- "gq":"Gluon-Quark/Antiquark",
- "qiqi":"Quark-Quark and $\mathrm{\overline{q_i} \; \overline{q_i}}$ (same flavor)",
- "qiai":"Quark-Antiquark (same flavor)",
- "qiqj":"Quark-Quark and $\mathrm{\overline{q_i} \; \overline{q_j}}}$ (different flavors)",
- "qiaj":"Quark-Antiquark (different flavors)"}
-
- #shorter for labeling the stackplot
- labeling = {"gg":r"$\mathrm{gg}$",
- "gq":r"$\mathrm{gq} \ & \ \mathrm{g\overline{q}}$",
- "qiqi":r"$\mathrm{q_i q_i} \ & \ \mathrm{\overline{q_i} \; \overline{q_i}}$",
- "qiai":r"$\mathrm{q_i} \; \mathrm{\overline{q_i}}$",
- "qiqj":r"$\mathrm{q_i q_j} \ & \ \mathrm{\overline{q_i} \; \overline{q_j}}$",
- "qiaj":r"$\mathrm{q_i} \; \mathrm{\overline{q_j}}$"}
-
-
- #parse arguments
- args = vars(parser.parse_args())
-
- #table name
- table_ = os.path.basename(args['table'])
- tablename = os.path.splitext(table_)[0]
- tablename = os.path.splitext(tablename)[0] #because yb tables have .tab.gz ending (getting rid of both here)
- print '\n'
- print 'Table: ', tablename, '\n'
-
- #pdfset name
- pdfset_ = os.path.basename(args['pdfset'])
- pdfname = os.path.splitext(pdfset_)[0]
- print 'PDF Set: ', pdfname, '\n'
-
- #checks for investigated phase space region --> mainly for plot-labeling
- #range of y_boost
- def yb():
- if 'yb0' in args['table']:
- yb_range = r'$\mathrm{0.0 \leq \/ y_b < 0.5}$'
- info_yb = r'0 <= yb < 0.5'
- elif 'yb1' in args['table']:
- yb_range = r'$\mathrm{0.5 \leq \/ y_b < 1.0}$'
- info_yb = r'0.5 <= yb < 1.0'
- elif 'yb2' in args['table']:
- yb_range = r'$\mathrm{1.0 \leq \/ y_b < 1.5}$'
- info_yb = r'1.0 <= yb < 1.5'
- elif 'yb3' in args['table']:
- yb_range = r'$\mathrm{1.5 \leq \/ y_b < 2.0}$'
- info_yb = r'1.5 <= yb < 2.0'
- elif 'yb4' in args['table']:
- yb_range = r'$\mathrm{2.0 \leq \/ y_b < 2.5}$'
- info_yb = r'2.0 <= yb < 2.5'
- else:
- print 'No info about range of yb.'
- yb_range = ''
- info_yb = ''
- return yb_range, info_yb #yb_range is in correct format for labeling
-
-
- #range of ystar
- def ystar():
- if 'ystar0' in args['table']:
- ystar_range = r'$\mathrm{0.0 \leq \/ y^{\ast} < 0.5}$'
- info_ystar = r'0 <= ystar < 0.5'
- elif 'ystar1' in args['table']:
- ystar_range = r'$\mathrm{0.5 \leq \/ y^{\ast} < 1.0}$'
- info_ystar = r'0.5 <= ystar < 1.0'
- elif 'ystar2' in args['table']:
- ystar_range = r'$\mathrm{1.0 \leq \/ y^{\ast} < 1.5}$'
- info_ystar = r'1.0 <= ystar < 1.5'
- elif 'ystar3' in args['table']:
- ystar_range = r'$\mathrm{1.5 \leq \/ y^{\ast} < 2.0}$'
- info_ystar = r'1.5 <= ystar < 2.0'
- elif 'ystar4' in args['table']:
- ystar_range = r'$\mathrm{2.0 \leq \/ y^{\ast} < 2.5}$'
- info_ystar = r'2.0 <= ystar < 2.5'
- else:
- print 'No info about range of ystar.'
- ystar_range = ''
- info_ystar = ''
- return ystar_range, info_ystar
-
-
-
- #show information about phase space region
- #
- yb_range = yb()[0]
- ystar_range = ystar()[0]
-
- info_yb = yb()[1]
- info_ystar = ystar()[1]
-
- print 'Phase Space Region:'
- print '%s' %info_yb
- print '%s' %info_ystar, '\n'
-
-
-
-
- #Which subprocesses are evaluated?
- print 'Subprocesses that are investigated: '
- print args['subproc']
-
-
- ######## EVALUATING ########
- fnlo = fastnlo.fastNLOLHAPDF(args['table'], args['pdfset'], args['member'])
- fnlo.SetContributionON(fastnlo.kFixedOrder,1,True);
-
-
- #cross section with every process included
- fnlo.CalcCrossSection()
- xs_all = np.array(fnlo.GetCrossSection())
- print '\n'
- print 'Cross section with all subprocesses xs_all: \n'
- print xs_all, '\n \n'
-
- #cross section of single processes
- xs_subproc_list = []
- for i in range(0, len(args['subproc'])):
- subproc_name = args['subproc'][i]
- print 'Selected subprocess: %s' %subproc_name
- fnlo.SelectProcesses(subproc_name)
- fnlo.CalcCrossSection()
- xs_subproc = np.array(fnlo.GetCrossSection())
- xs_subproc_list.append(xs_subproc)
- print 'XS for subprocess %s: \n' %args['subproc'][i]
- print xs_subproc, '\n \n'
- xs_sub = np.array(xs_subproc_list) #Array with xs for selected processes in certain pTZ-bins
- print 'XS of the selected processes in different pTZ-regions: \n'
- print xs_sub, '\n \n'
-
- #how much do the selected subprocesses contribute to the total xs?
- fractions = np.divide(xs_sub,xs_all)
- print 'fractions: \n', fractions, '\n'
-
- for i in range(0, len(fractions)):
- print 'fraction of %s on total XS: ' %args['subproc'][i]
- print fractions[i,:], '\n \n'
-
-
- #checking whether contribution is positive or negative
- pos_contr = np.array(fractions)
- neg_contr = np.array(fractions)
- bin_bounds = np.array(fnlo.GetObsBinsBounds(0))
-
- for i in range(0, len(args['subproc'])): #go through 'fractions' line by line = subprocess after subprocess
- for j in range(0, len(bin_bounds)): #check the different pTZ bins
- if (fractions[i,j] < 0):
- pos_contr[i,j] = 0.0
- else:
- neg_contr[i,j] = 0.0
- print 'Positive contributions of %s on total XS: \n' %args['subproc'][i]
- print pos_contr[i,:], '\n'
- print 'Negative contributions of %s on total XS: \n' %args['subproc'][i]
- print neg_contr[i,:], '\n \n'
-
- print 'Summary of positive and negative contributions of the subprocesses: \n'
- print pos_contr, '\n'
- print neg_contr, '\n'
-
-
-
- ######## PLOTTING ########
- #Bins for x-axis (according to first pdf (0))
- bin_bounds = np.array(fnlo.GetObsBinsBounds(0))
- #print bin_bounds.T #bin_bounds.T[0] = lower bounds, bin_bounds.T[1] upper bounds
- #print bin_bounds.flatten()
-
-
- #dictionary for using specific color per process
- process_color = {"gg":"g",
- "gq":"r",
- "qiqi":"y",
- "qiai":"c",
- "qiqj":"m",
- "qiaj":"b"}
-
- #Check whether comparison-plot or single-plot is required:
- print args['compare']
- if (args['compare']==True): #stackplot that compares the contribution of the subprocesses to the total XS
- plt.close()
- fig0 = plt.figure(figsize=(8,7))
- ax0 = fig0.add_subplot(111)
-
- patches = [] #later needed for legend
- #number of bins via bin_bounds variable
- y0_pos = np.zeros([len(bin_bounds)]) #for tracking the current 'height' of stackplot, in order to staple and not overlap the values for each subprocess
- y0_neg = np.zeros([len(bin_bounds)]) #same as above for the negative contributions, plotted below x-axis
-
- for i in range(0, len(args['subproc'])): #go through fractions, subproc by subproc
-
- #positive contributions
- ax0.fill_between(bin_bounds.flatten(), steppify_bin(y0_pos), steppify_bin(y0_pos+pos_contr[i]), #label=labeling[args['subproc'][i]],
- facecolor=process_color[args['subproc'][i]], alpha=0.4)
- y0_pos += pos_contr[i] #calculate new height of stackplot
-
- #negative contributions
- ax0.fill_between(bin_bounds.flatten(), steppify_bin(y0_neg), steppify_bin(y0_neg+neg_contr[i]), #label=labeling[args['subproc'][i]],
- facecolor=process_color[args['subproc'][i]], hatch='X', alpha=0.4)
- y0_neg += neg_contr[i] #new 'height' below x-axis
-
- #patch for subprocess i (for the legend) ##labeling happens here
- patches.append(matplotlib.patches.Rectangle((0, 0), 0, 0, color=process_color[args['subproc'][i]],
- label=labeling[args['subproc'][i]], alpha=0.4))
- ax0.add_patch(patches[i])
-
- #settings for the whole stackplot
- ax0.set_xscale('log', nonposx='clip')
- ax0.axis([30, 1000, -0.20, 1.20])
- plt.axhline(y=1, xmin=0, xmax=1, color='k', linestyle='dotted') #dotted line at xs_sub/xs=1=100%
- plt.axhline(y=0, xmin=0, xmax=1, color='k', linestyle='-') #line at xs_sub/xs_tot=0
-
- #use phasespace region as title
- title_phasespace = yb_range + ", " + ystar_range
-
- plt.title(x=0.5, y=1.01, s='%s' %(title_phasespace), fontsize=22)
-
- ax0.set_xlabel('$\mathrm{p_{T,Z}} \ \mathrm{[GeV]}$', fontsize=20)
- ax0.xaxis.set_label_coords(0.9, -0.05)
- ax0.set_ylabel(r'$\frac{\mathrm{XS_{subproc}}}{\mathrm{XS_{total}}}$', rotation=0, fontsize=24)
- ax0.yaxis.set_label_coords(-0.16, 0.9)
-
- ax0.text(0.96, 0.03, args['pdfset']+', NLO', transform=ax0.transAxes, fontsize=14, verticalalignment='bottom', horizontalalignment='right')
- plt.legend()
-
- plt.tight_layout()
- fig0.tight_layout()
-
- #make string with the subprocess names
- processnames = ""
- for i in range(0, len(args['subproc'])):
- processnames += ("_"+args['subproc'][i])
- stackplotname = 'compare%s_%s_%s.png' %(processnames, tablename, args['pdfset'])
- fig0.savefig(stackplotname)
-
- print 'Stackplot for subprocess comparison saved as: %s' %(stackplotname)
-
-
-
-
- else: #individual plot for each subproc
- #create folder for saving the plots
- current_dir = os.getcwd()
-
- final_dir = os.path.join(current_dir, r'single_subproc_plots_%s' %tablename)
- now = datetime.datetime.now()
- nowstr = 'single_subproc_plots_%s_' %(tablename)
- nowstr += now.strftime('%Y-%m-%d_%H-%M-%S')
- final_dir_date = os.path.join(current_dir, nowstr)
-
- #check if directory already exists, if so: include datetime to name
- if not os.path.exists(final_dir):
- os.makedirs(final_dir)
- else:
- os.makedirs(final_dir_date)
-
- #loop where single plot for each of the selected subprocesses is created
- #start the plot
- for i in range(0, len(args['subproc'])):
- plt.close()
- fig0 = plt.figure(figsize=(8, 7))
- ax0 = fig0.add_subplot(111)
-
- #plot the fractions
- y0 = np.zeros([len(bin_bounds)])
- ax0.fill_between(bin_bounds.flatten(), steppify_bin(y0), steppify_bin(fractions[i]),
- facecolor=process_color[args['subproc'][i]], alpha=0.4) #improve labeling?
-
-
- #info about phase space region included as text (see below)
- ax0.set_xscale('log', nonposx='clip')
- ax0.axis([30, 1000, -0.4, 1.2])
- plt.axhline(y=1, xmin=0, xmax=1, color='k', linestyle='dotted')
- plt.axhline(y=0, xmin=0, xmax=1, color='k', linestyle='-')
-
- titlename = processes[args['subproc'][i]]
-
- plt.title(x=0.5, y=1.01, s=titlename, fontsize=20)
- ax0.set_xlabel('$\mathrm{p_{T,Z}} \ \mathrm{[GeV]}$', fontsize=20)
- #ax0.xaxis.set_label_coords(0.9, -0.05)
- ax0.set_ylabel('$\mathrm{XS_{subproc}} / \mathrm{XS_{total}}$', fontsize=20)
- ax0.text(0.96, 0.03, args['pdfset'], transform=ax0.transAxes, fontsize=14, verticalalignment='bottom', horizontalalignment='right')
- ax0.text(-0.12, 0.98, 'NLO', transform=ax0.transAxes, fontsize=16, rotation=90, verticalalignment='top', horizontalalignment='left') #manually added 'NLO' to y-label
- ax0.text(0.96, 0.94, yb_range, transform=ax0.transAxes, fontsize=16, verticalalignment='bottom', horizontalalignment='right')
- ax0.text(0.96, 0.88, ystar_range, transform=ax0.transAxes, fontsize=16, verticalalignment='bottom', horizontalalignment='right')
-
- plt.legend() #location parameter loc='upper right'
-
- plt.tight_layout()
- fig0.tight_layout()
- plotname = '%s_frac_%s_%s.png' %(args['subproc'][i], tablename, args['pdfset'])
-
-
- #save the plots
- if os.path.exists(final_dir_date):
- fig0.savefig(os.path.join(final_dir_date,plotname))
- print 'Plot saved in %s' %(final_dir_date)
- else:
- fig0.savefig(os.path.join(final_dir,plotname))
- print 'Plot saved in %s' %(final_dir)
-
- print '\n'
- print 'fraction plot saved as: %s' %(plotname)
-
-
-
-
-
-##################################################
-##function to make better uncertainty plot (shaded)
-##could maybe be replaced by .flatten() ??
-def steppify_bin(arr, isx=False):
- """
- Produce stepped array of arr, needed for example for stepped fill_betweens.
- Pass all x bin edges to produce stepped x arr and all y bincontents to produce
- stepped bincontents representation
- steppify_bin([1,2,3], True)
- -> [1,2,2,3]
- steppify_bin([5,6])
- -> [5,5,6,6]
- """
- if isx:
- newarr = np.array(zip(arr[:-1], arr[1:])).ravel()
- else:
- newarr = np.array(zip(arr, arr)).ravel()
- return newarr
-
-
-#################################################
-
-if __name__ == '__main__':
- main()
-
-
-
-
-
--
GitLab