diff options
Diffstat (limited to 'utils/plot-fit-results')
| -rwxr-xr-x | utils/plot-fit-results | 332 |
1 files changed, 44 insertions, 288 deletions
diff --git a/utils/plot-fit-results b/utils/plot-fit-results index 18139ff..e154b21 100755 --- a/utils/plot-fit-results +++ b/utils/plot-fit-results @@ -15,265 +15,62 @@ # this program. If not, see <https://www.gnu.org/licenses/>. from __future__ import print_function, division -import numpy as np -from scipy.stats import iqr, norm, beta -from matplotlib.lines import Line2D -import itertools - -IDP_E_MINUS = 20 -IDP_MU_MINUS = 22 - -SNOMAN_MASS = { - 20: 0.511, - 21: 0.511, - 22: 105.658, - 23: 105.658 -} - -def plot_hist(x, label=None): - # determine the bin width using the Freedman Diaconis rule - # see https://en.wikipedia.org/wiki/Freedman%E2%80%93Diaconis_rule - h = 2*iqr(x)/len(x)**(1/3) - n = max(int((np.max(x)-np.min(x))/h),10) - bins = np.linspace(np.min(x),np.max(x),n) - plt.hist(x, bins=bins, histtype='step', label=label) - -def plot_legend(n): - plt.figure(n) - ax = plt.gca() - handles, labels = ax.get_legend_handles_labels() - new_handles = [Line2D([],[],c=h.get_edgecolor()) for h in handles] - plt.legend(handles=new_handles,labels=labels) - -def get_stats(x): - """ - Returns a tuple (mean, error mean, std, error std) for the values in x. - - The formula for the standard error on the standard deviation comes from - https://stats.stackexchange.com/questions/156518. - """ - mean = np.mean(x) - std = np.std(x) - n = len(x) - u4 = np.mean((x-mean)**4) - error = np.sqrt((u4-(n-3)*std**4/(n-1))/n)/(2*std) - return mean, std/np.sqrt(n), std, error - -def iqr_std_err(x): - """ - Returns the approximate standard deviation assuming the central part of the - distribution is gaussian. - """ - x = x.dropna() - n = len(x) - if n == 0: - return np.nan - # see https://stats.stackexchange.com/questions/110902/error-on-interquartile-range - std = iqr(x.values)/1.3489795 - return 1.573*std/np.sqrt(n) - -def iqr_std(x): - """ - Returns the approximate standard deviation assuming the central part of the - distribution is gaussian. - """ - x = x.dropna() - n = len(x) - if n == 0: - return np.nan - return iqr(x.values)/1.3489795 - -def quantile_error(x,q): - """ - Returns the standard error for the qth quantile of `x`. The error is - computed using the Maritz-Jarrett method described here: - https://www.itl.nist.gov/div898/software/dataplot/refman2/auxillar/quantse.htm. - """ - x = np.sort(x) - n = len(x) - m = int(q*n+0.5) - A = m - 1 - B = n - m - i = np.arange(1,len(x)+1) - w = beta.cdf(i/n,A,B) - beta.cdf((i-1)/n,A,B) - return np.sqrt(np.sum(w*x**2)-np.sum(w*x)**2) - -def q90_err(x): - """ - Returns the error on the 90th percentile for all the non NaN values in a - Series `x`. - """ - x = x.dropna() - n = len(x) - if n == 0: - return np.nan - return quantile_error(x.values,0.9) - -def q90(x): - """ - Returns the 90th percentile for all the non NaN values in a Series `x`. - """ - x = x.dropna() - n = len(x) - if n == 0: - return np.nan - return np.percentile(x.values,90.0) - -def median(x): - """ - Returns the median for all the non NaN values in a Series `x`. - """ - x = x.dropna() - n = len(x) - if n == 0: - return np.nan - return np.median(x.values) - -def median_err(x): - """ - Returns the approximate error on the median for all the non NaN values in a - Series `x`. The error on the median is approximated assuming the central - part of the distribution is gaussian. - """ - x = x.dropna() - n = len(x) - if n == 0: - return np.nan - # First we estimate the standard deviation using the interquartile range. - # Here we are essentially assuming the central part of the distribution is - # gaussian. - std = iqr(x.values)/1.3489795 - median = np.median(x.values) - # Now we estimate the error on the median for a gaussian - # See https://stats.stackexchange.com/questions/45124/central-limit-theorem-for-sample-medians. - return 1/(2*np.sqrt(n)*norm.pdf(median,median,std)) - -def std_err(x): - x = x.dropna() - mean = np.mean(x) - std = np.std(x) - n = len(x) - if n == 0: - return np.nan - elif n == 1: - return 0.0 - u4 = np.mean((x-mean)**4) - error = np.sqrt((u4-(n-3)*std**4/(n-1))/n)/(2*std) - return error - -# Taken from https://raw.githubusercontent.com/mwaskom/seaborn/c73055b2a9d9830c6fbbace07127c370389d04dd/seaborn/utils.py -def despine(fig=None, ax=None, top=True, right=True, left=False, - bottom=False, offset=None, trim=False): - """Remove the top and right spines from plot(s). - - fig : matplotlib figure, optional - Figure to despine all axes of, default uses current figure. - ax : matplotlib axes, optional - Specific axes object to despine. - top, right, left, bottom : boolean, optional - If True, remove that spine. - offset : int or dict, optional - Absolute distance, in points, spines should be moved away - from the axes (negative values move spines inward). A single value - applies to all spines; a dict can be used to set offset values per - side. - trim : bool, optional - If True, limit spines to the smallest and largest major tick - on each non-despined axis. - - Returns - ------- - None - - """ - # Get references to the axes we want - if fig is None and ax is None: - axes = plt.gcf().axes - elif fig is not None: - axes = fig.axes - elif ax is not None: - axes = [ax] - - for ax_i in axes: - for side in ["top", "right", "left", "bottom"]: - # Toggle the spine objects - is_visible = not locals()[side] - ax_i.spines[side].set_visible(is_visible) - if offset is not None and is_visible: - try: - val = offset.get(side, 0) - except AttributeError: - val = offset - _set_spine_position(ax_i.spines[side], ('outward', val)) - - # Potentially move the ticks - if left and not right: - maj_on = any( - t.tick1line.get_visible() - for t in ax_i.yaxis.majorTicks - ) - min_on = any( - t.tick1line.get_visible() - for t in ax_i.yaxis.minorTicks - ) - ax_i.yaxis.set_ticks_position("right") - for t in ax_i.yaxis.majorTicks: - t.tick2line.set_visible(maj_on) - for t in ax_i.yaxis.minorTicks: - t.tick2line.set_visible(min_on) - - if bottom and not top: - maj_on = any( - t.tick1line.get_visible() - for t in ax_i.xaxis.majorTicks - ) - min_on = any( - t.tick1line.get_visible() - for t in ax_i.xaxis.minorTicks - ) - ax_i.xaxis.set_ticks_position("top") - for t in ax_i.xaxis.majorTicks: - t.tick2line.set_visible(maj_on) - for t in ax_i.xaxis.minorTicks: - t.tick2line.set_visible(min_on) - - if trim: - # clip off the parts of the spines that extend past major ticks - xticks = ax_i.get_xticks() - if xticks.size: - firsttick = np.compress(xticks >= min(ax_i.get_xlim()), - xticks)[0] - lasttick = np.compress(xticks <= max(ax_i.get_xlim()), - xticks)[-1] - ax_i.spines['bottom'].set_bounds(firsttick, lasttick) - ax_i.spines['top'].set_bounds(firsttick, lasttick) - newticks = xticks.compress(xticks <= lasttick) - newticks = newticks.compress(newticks >= firsttick) - ax_i.set_xticks(newticks) - - yticks = ax_i.get_yticks() - if yticks.size: - firsttick = np.compress(yticks >= min(ax_i.get_ylim()), - yticks)[0] - lasttick = np.compress(yticks <= max(ax_i.get_ylim()), - yticks)[-1] - ax_i.spines['left'].set_bounds(firsttick, lasttick) - ax_i.spines['right'].set_bounds(firsttick, lasttick) - newticks = yticks.compress(yticks <= lasttick) - newticks = newticks.compress(newticks >= firsttick) - ax_i.set_yticks(newticks) if __name__ == '__main__': import argparse import numpy as np import h5py import pandas as pd + import itertools + from sddm import IDP_E_MINUS, IDP_MU_MINUS, SNOMAN_MASS + from sddm.plot import plot_hist, plot_legend, get_stats, despine, iqr_std_err, iqr_std, quantile_error, q90_err, q90, median, median_err, std_err parser = argparse.ArgumentParser("plot fit results") parser.add_argument("filenames", nargs='+', help="input files") parser.add_argument("--save", action="store_true", default=False, help="save plots") args = parser.parse_args() + if args.save: + # default \textwidth for a fullpage article in Latex is 16.50764 cm. + # You can figure this out by compiling the following TeX document: + # + # \documentclass{article} + # \usepackage{fullpage} + # \usepackage{layouts} + # \begin{document} + # textwidth in cm: \printinunitsof{cm}\prntlen{\textwidth} + # \end{document} + + width = 16.50764 + width /= 2.54 # cm -> inches + # According to this page: + # http://www-personal.umich.edu/~jpboyd/eng403_chap2_tuftegospel.pdf, + # Tufte suggests an aspect ratio of 1.5 - 1.6. + height = width/1.5 + FIGSIZE = (width,height) + + import matplotlib.pyplot as plt + + font = {'family':'serif', 'serif': ['computer modern roman']} + plt.rc('font',**font) + + plt.rc('text', usetex=True) + else: + # on retina screens, the default plots are way too small + # by using Qt5 and setting QT_AUTO_SCREEN_SCALE_FACTOR=1 + # Qt5 will scale everything using the dpi in ~/.Xresources + import matplotlib + matplotlib.use("Qt5Agg") + + import matplotlib.pyplot as plt + + # Default figure size. Currently set to my monitor width and height so that + # things are properly formatted + FIGSIZE = (13.78,7.48) + + # Make the defalt font bigger + plt.rc('font', size=22) + # Read in all the data. # # Note: We have to add the filename as a column to each dataset since this @@ -351,47 +148,6 @@ if __name__ == '__main__': markers = itertools.cycle(('o', 'v')) - if args.save: - # default \textwidth for a fullpage article in Latex is 16.50764 cm. - # You can figure this out by compiling the following TeX document: - # - # \documentclass{article} - # \usepackage{fullpage} - # \usepackage{layouts} - # \begin{document} - # textwidth in cm: \printinunitsof{cm}\prntlen{\textwidth} - # \end{document} - - width = 16.50764 - width /= 2.54 # cm -> inches - # According to this page: - # http://www-personal.umich.edu/~jpboyd/eng403_chap2_tuftegospel.pdf, - # Tufte suggests an aspect ratio of 1.5 - 1.6. - height = width/1.5 - FIGSIZE = (width,height) - - import matplotlib.pyplot as plt - - font = {'family':'serif', 'serif': ['computer modern roman']} - plt.rc('font',**font) - - plt.rc('text', usetex=True) - else: - # on retina screens, the default plots are way too small - # by using Qt5 and setting QT_AUTO_SCREEN_SCALE_FACTOR=1 - # Qt5 will scale everything using the dpi in ~/.Xresources - import matplotlib - matplotlib.use("Qt5Agg") - - import matplotlib.pyplot as plt - - # Default figure size. Currently set to my monitor width and height so that - # things are properly formatted - FIGSIZE = (13.78,7.48) - - # Make the defalt font bigger - plt.rc('font', size=22) - fig3, ax3 = plt.subplots(3,1,figsize=FIGSIZE,num=3,sharex=True) fig4, ax4 = plt.subplots(3,1,figsize=FIGSIZE,num=4,sharex=True) |