update contamination closure test script

- print out mean acceptance fraction and autocorrelation time
- print standard normal distribution pdf with pull plot histograms

1 files changed, 13 insertions, 2 deletions

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))