update plot-muons to fit energy scale and resolution difference

2 files changed, 73 insertions, 28 deletions

diff --git a/utils/chi2 b/utils/chi2
index 9fd4151..f969498 100755
--- a/utils/chi2
+++ b/utils/chi2
@@ -48,7 +48,7 @@ def printoptions(*args, **kwargs):
 # FIXME: These are just placeholders! Should get real number from stopping
 # muons.
 ENERGY_SCALE_MEAN = {'e': 1.0, 'u': 1.0, 'eu': 1.0}
-ENERGY_SCALE_UNCERTAINTY = {'e':0.1, 'u': 0.1, 'eu': 0.1}
+ENERGY_SCALE_UNCERTAINTY = {'e':0.1, 'u': 1.1, 'eu': 0.1}
 ENERGY_RESOLUTION_MEAN = {'e': 0.0, 'u': 0.0, 'eu': 0.0}
 ENERGY_RESOLUTION_UNCERTAINTY = {'e':0.1, 'u': 0.1, 'eu': 0.1}
 
diff --git a/utils/plot-muons b/utils/plot-muons
index d14f8b6..5217299 100755
--- a/utils/plot-muons
+++ b/utils/plot-muons
@@ -27,6 +27,8 @@ import numpy as np
 from scipy.stats import iqr, poisson
 from scipy.stats import iqr, norm, beta
 from sddm.stats import *
+import nlopt
+from sddm.dc import estimate_errors
 
 particle_id = {20: 'e', 22: 'u'}
 
@@ -42,6 +44,28 @@ def print_particle_probs(data):
         particle_id_str = ''.join([particle_id[int(''.join(x))] for x in grouper(str(id),2)])
         print("P(%s) = %.2f +/- %.2f" % (particle_id_str,mode[i]*100,std[i]*100))
 
+def fit_straight_line(y, yerr):
+    def nll(x, grad=None):
+        nll = -norm.logpdf(x,y,yerr).sum()
+        return nll
+
+    x0 = np.array([np.mean(y)])
+    opt = nlopt.opt(nlopt.LN_SBPLX, len(x0))
+    opt.set_min_objective(nll)
+    low = np.array([-1e9])
+    high = np.array([1e9])
+    opt.set_lower_bounds(low)
+    opt.set_upper_bounds(high)
+    opt.set_ftol_abs(1e-10)
+    opt.set_initial_step([0.01])
+
+    xopt = opt.optimize(x0)
+    nll_xopt = nll(xopt)
+
+    stepsizes = estimate_errors(nll,xopt,low,high)
+
+    return xopt[0], stepsizes[0]
+
 if __name__ == '__main__':
     import argparse
     import numpy as np
@@ -131,15 +155,15 @@ if __name__ == '__main__':
     fig = plt.figure()
     plt.subplot(2,1,1)
     plt.hist(muons.ke.values, bins=np.logspace(3,7,100), histtype='step', color='C0', label="Data")
-    norm = len(muons.ke.values)/len(muons_mc.ke.values)
-    plt.hist(muons_mc.ke.values, weights=np.tile(norm,len(muons_mc.ke.values)), bins=np.logspace(3,7,100), histtype='step', color='C1', label="Monte Carlo")
+    scale = len(muons.ke.values)/len(muons_mc.ke.values)
+    plt.hist(muons_mc.ke.values, weights=np.tile(scale,len(muons_mc.ke.values)), bins=np.logspace(3,7,100), histtype='step', color='C1', label="Monte Carlo")
     plt.legend()
     plt.xlabel("Energy (MeV)")
     plt.gca().set_xscale("log")
     plt.subplot(2,1,2)
     plt.hist(np.cos(muons.theta.values), bins=np.linspace(-1,1,100), histtype='step', color='C0', label="Data")
-    norm = len(muons.theta.values)/len(muons_mc.theta.values)
-    plt.hist(np.cos(muons_mc.theta.values), weights=np.tile(norm,len(muons_mc.ke.values)), bins=np.linspace(-1,1,100), histtype='step', color='C1', label="Monte Carlo")
+    scale = len(muons.theta.values)/len(muons_mc.theta.values)
+    plt.hist(np.cos(muons_mc.theta.values), weights=np.tile(scale,len(muons_mc.ke.values)), bins=np.linspace(-1,1,100), histtype='step', color='C1', label="Monte Carlo")
     plt.legend()
     despine(fig,trim=True)
     plt.xlabel(r"$\cos(\theta)$")
@@ -172,15 +196,15 @@ if __name__ == '__main__':
     fig = plt.figure()
     plt.subplot(2,1,1)
     plt.hist(stopping_muons.ke.values, bins=np.logspace(3,7,100), histtype='step', color='C0', label="Data")
-    norm = len(stopping_muons.ke.values)/len(stopping_muons_mc.ke.values)
-    plt.hist(stopping_muons_mc.ke.values, weights=np.tile(norm,len(stopping_muons_mc.ke.values)), bins=np.logspace(3,7,100), histtype='step', color='C1', label="Monte Carlo")
+    scale = len(stopping_muons.ke.values)/len(stopping_muons_mc.ke.values)
+    plt.hist(stopping_muons_mc.ke.values, weights=np.tile(scale,len(stopping_muons_mc.ke.values)), bins=np.logspace(3,7,100), histtype='step', color='C1', label="Monte Carlo")
     plt.legend()
     plt.xlabel("Energy (MeV)")
     plt.gca().set_xscale("log")
     plt.subplot(2,1,2)
     plt.hist(np.cos(stopping_muons.theta.values), bins=np.linspace(-1,1,100), histtype='step', color='C0', label="Data")
-    norm = len(stopping_muons.theta.values)/len(stopping_muons_mc.theta.values)
-    plt.hist(np.cos(stopping_muons_mc.theta.values), weights=np.tile(norm,len(stopping_muons_mc.theta.values)), bins=np.linspace(-1,1,100), histtype='step', color='C1', label="Monte Carlo")
+    scale = len(stopping_muons.theta.values)/len(stopping_muons_mc.theta.values)
+    plt.hist(np.cos(stopping_muons_mc.theta.values), weights=np.tile(scale,len(stopping_muons_mc.theta.values)), bins=np.linspace(-1,1,100), histtype='step', color='C1', label="Monte Carlo")
     plt.legend()
     despine(fig,trim=True)
     plt.xlabel(r"$\cos(\theta)$")
@@ -236,27 +260,48 @@ if __name__ == '__main__':
     dT = stopping_muons.groupby(pd_bins)['dT'].agg(['mean','sem','std',std_err,median,median_err,iqr_std,iqr_std_err])
     dT_mc = stopping_muons_mc.groupby(pd_bins_mc)['dT'].agg(['mean','sem','std',std_err,median,median_err,iqr_std,iqr_std_err])
 
-    fig = plt.figure()
-    plt.errorbar(T, dT['median']*100/T, yerr=dT['median_err']*100/T, color='C0', label="Data")
-    plt.errorbar(T, dT_mc['median']*100/T, yerr=dT_mc['median_err']*100/T, color='C1', label="Monte Carlo")
-    plt.legend()
-    despine(fig,trim=True)
-    plt.xlabel("Kinetic Energy (MeV)")
-    plt.ylabel(r"Energy bias (\%)")
+    y = (dT['median']*100/T-dT_mc['median']*100/T).values
+    yerr = np.sqrt((dT['median_err']*100/T)**2+(dT_mc['median_err']*100/T)**2).values
+    mean, std = fit_straight_line(y,yerr)
+
+    print("Data energy bias = %.2f +/- %.2f" % (mean, std))
+
+    fig, (a0, a1) = plt.subplots(2, 1, gridspec_kw={'height_ratios': [3, 1]})
+    a0.errorbar(T, dT['median']*100/T, yerr=dT['median_err']*100/T, color='C0', label="Data")
+    a0.errorbar(T, dT_mc['median']*100/T, yerr=dT_mc['median_err']*100/T, color='C1', label="Monte Carlo")
+    despine(ax=a0,trim=True)
+    a0.set_ylabel(r"Energy bias (\%)")
+    a0.legend()
+    a1.errorbar(T, dT['median']*100/T-dT_mc['median']*100/T, yerr=np.sqrt((dT['median_err']*100/T)**2+(dT_mc['median_err']*100/T)**2), color='C0')
+    a1.axhline(mean,T[0],T[-1],ls='--',color='r')
+    a1.set_ylim(0,25)
+    despine(ax=a1,trim=True)
+    a1.set_xlabel("Kinetic Energy (MeV)")
+    a1.set_ylabel(r"Difference (\%)")
     plt.tight_layout()
     if args.save:
         plt.savefig("stopping_muon_energy_bias.pdf")
         plt.savefig("stopping_muon_energy_bias.eps")
     else:
-        plt.title("Stopping Muon Energy Bias")
-
-    fig = plt.figure()
-    plt.errorbar(T, dT['iqr_std']*100/T, yerr=dT['iqr_std_err']*100/T, color='C0', label="Data")
-    plt.errorbar(T, dT_mc['iqr_std']*100/T, yerr=dT_mc['iqr_std_err']*100/T, color='C1', label="Monte Carlo")
-    plt.legend()
-    despine(fig,trim=True)
-    plt.xlabel("Kinetic Energy (MeV)")
-    plt.ylabel(r"Energy resolution (\%)")
+        plt.suptitle("Stopping Muon Energy Bias")
+
+    y = (dT['iqr_std']*100/T-dT_mc['iqr_std']*100/T).values
+    yerr = np.sqrt((dT['iqr_std_err']*100/T)**2+(dT_mc['iqr_std_err']*100/T)**2).values
+    mean, std = fit_straight_line(y,yerr)
+
+    print("Data energy resolution = %.2f +/- %.2f" % (mean, std))
+
+    fig, (a0, a1) = plt.subplots(2, 1, gridspec_kw={'height_ratios': [3, 1]})
+    a0.errorbar(T, dT['iqr_std']*100/T, yerr=dT['iqr_std_err']*100/T, color='C0', label="Data")
+    a0.errorbar(T, dT_mc['iqr_std']*100/T, yerr=dT_mc['iqr_std_err']*100/T, color='C1', label="Monte Carlo")
+    a0.set_ylabel(r"Energy resolution (\%)")
+    despine(ax=a0,trim=True)
+    a0.legend()
+    a1.errorbar(T, dT['iqr_std']*100/T-dT_mc['iqr_std']*100/T, yerr=np.sqrt((dT['iqr_std_err']*100/T)**2+(dT_mc['iqr_std_err']*100/T)**2), color='C0')
+    a1.axhline(mean,T[0],T[-1],ls='--',color='r')
+    despine(ax=a1,trim=True)
+    a1.set_xlabel("Kinetic Energy (MeV)")
+    a1.set_ylabel(r"Difference (\%)")
     plt.tight_layout()
     if args.save:
         plt.savefig("stopping_muon_energy_resolution.pdf")
@@ -295,10 +340,10 @@ if __name__ == '__main__':
     hist = np.histogram(michel_low_nhit.ke.values,bins=bins)[0]
     hist_mc = np.histogram(michel_low_nhit_mc.ke.values,bins=bins)[0]
     if hist_mc.sum() > 0:
-        norm = hist.sum()/hist_mc.sum()
+        scale = hist.sum()/hist_mc.sum()
     else:
-        norm = 1.0
-    p = get_multinomial_prob(hist,hist_mc,norm)
+        scale = 1.0
+    p = get_multinomial_prob(hist,hist_mc,scale)
     plt.text(0.95,0.95,"p = %.2f" % p,horizontalalignment='right',verticalalignment='top',transform=plt.gca().transAxes)
     despine(fig,trim=True)
     plt.xlabel("Energy (MeV)")