v0.2.4 : Calculate 2nd Order Energy Correction

This calculation was done using a 1M event sample generated by shooting 4GeV electrons directly into the ECal in order to avoid any upstream effects. The fit was done with the code shown below which also produced a plot to visually confirm a decent mean was found. The value of 3940.5 was used.

Code

import uproot
import numpy as np
from scipy.stats import norm
import scipy
%matplotlib inline
import matplotlib.pyplot as plt
import mplhep
plt.style.use(mplhep.style.ROOT)

t = uproot.open('data/v3.2.0-gamma-14-geb86ea1/plain/type_events_sim_plain_geometry_v14_events_1000000_run_1.root:LDMX_Events')
tre = t['EcalVeto_valid/summedDet_'].array(library='np')

vals, edges, p = plt.hist(tre,bins='auto', range=(0,8000), density=True,
                         histtype='step',linewidth=2, label='1M Plain 4GeV')
centers = (edges[1:]+edges[:-1])/2

def two_sided_normal(x, mu, low_sig, high_sig) :
    low_side = norm.pdf(x[x<mu], mu, low_sig)
    high_side = norm.pdf(x[x>=mu], mu, high_sig)
    return np.concatenate((low_side,high_side))

# filter to reasonable values, specifically dropping high energy events
# from wacky hits late in the detector and low energy events probably
# due to PN interactions
high_cut = 6000
low_cut = 2500
tre_f = tre[(tre < high_cut)&(tre > low_cut)]

mean = tre_f.mean()
stdd = np.sqrt(((tre_f - mean)**2).mean())

bin_selection = (centers < high_cut) & (centers > low_cut)

optim, cov = scipy.optimize.curve_fit(two_sided_normal, centers[bin_selection], vals[bin_selection], 
                                      sigma=np.sqrt(vals[bin_selection]), p0=[mean,stdd,stdd])

plt.plot(centers, norm.pdf(centers, mean, stdd), 
         label=f'$\mu = {mean:.1f}$MeV $\sigma = {stdd:.1f}$MeV')
plt.plot(centers, two_sided_normal(centers, *optim), 
         label=f'$\mu = {optim[0]:.1f}$MeV\n$\sigma_l = {optim[1]:.1f}$MeV\n$\sigma_h = {optim[2]:.1f}$MeV')
plt.yscale('log')
plt.ylim(1e-9,1e-2)
plt.legend(loc='upper left')
plt.xlabel('Total Rec Energy [MeV]')
plt.ylabel('Event Fraction')
plt.show()

Full Changelog: v0.2.3...v0.2.4