From 8e82dc389eb5f67ccd0240e9b0bc165a83b27564 Mon Sep 17 00:00:00 2001 From: tlatorre Date: Tue, 2 Jun 2020 10:51:58 -0500 Subject: update contamination closure test script - print out mean acceptance fraction and autocorrelation time - print standard normal distribution pdf with pull plot histograms --- utils/dc-closure-test | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) (limited to 'utils/dc-closure-test') diff --git a/utils/dc-closure-test b/utils/dc-closure-test index 43d26d3..5b9f77d 100755 --- a/utils/dc-closure-test +++ b/utils/dc-closure-test @@ -18,7 +18,7 @@ from __future__ import print_function, division import numpy as np from scipy.stats import iqr import nlopt -from scipy.stats import poisson +from scipy.stats import poisson, norm import contextlib import sys from math import exp @@ -400,6 +400,13 @@ def fit(data, sacrifice, steps): with np.errstate(invalid='ignore'): sampler.run_mcmc(pos, steps) + print("Mean acceptance fraction: {0:.3f}".format(np.mean(sampler.acceptance_fraction))) + + try: + print("autocorrelation time: ", sampler.get_autocorr_time(quiet=True)) + except Exception as e: + print(e) + samples = sampler.chain.reshape((-1,len(x0))) return samples @@ -539,11 +546,15 @@ if __name__ == '__main__': std = np.std(contamination) contamination_pull[bg].append((mean - nbg[bg])/std) + bins = np.linspace(-10,10,101) + bincenters = (bins[1:] + bins[:-1])/2 + fig = plt.figure() axes = [] for i, bg in enumerate(['signal','muon','noise','neck','flasher','breakdown']): axes.append(plt.subplot(3,2,i+1)) - plt.hist(contamination_pull[bg],bins=100,histtype='step') + plt.hist(contamination_pull[bg],bins=bins,histtype='step',normed=True) + plt.plot(bincenters,norm.pdf(bincenters)) plt.title(bg.capitalize()) for ax in axes: ax.set_xlim((-10,10)) -- cgit