<jupyter_start><jupyter_text>SummaryThe purpose of this book is to show how the extracted Fano factor, combined with the expected yield with plus multiple-scatter modification, can predict the trend of the EDELWEISS measured ionization yield width. And, how that final predicted width function can be used to extract the systematic deviation due to the finite binning of the EDELWEISS data.   Extracted Fano Factor Matching<jupyter_code>import h5py
filename = 'data/systematic_error_fits.h5'
#remove vars
f = h5py.File(filename,'r')
#for i in f['mcmc/edwdata_sys_error']:
#    print(i)

#save the results for the Edw fit
path='{}/{}/'.format('mcmc','edwdata_sys_error')

Cms = np.asarray(f[path+'Cms'])
slope = np.asarray(f[path+'m'])
a_yield = np.asarray(f[path+'A'])
b_yield = np.asarray(f[path+'B'])
aH = np.asarray(f[path+'aH'])
scale = np.asarray(f[path+'scale'])
samples = np.asarray(f[path+'samples'])
sampsize = np.asarray(f[path+'sampsize'])
xl = np.asarray(f[path+'Er'])
upvec = np.asarray(f[path+'Csig_u'])
dnvec = np.asarray(f[path+'Csig_l'])
Sigtot = np.asarray(f[path+'Sigss'])
Sigss = np.sqrt(Sigtot**2 - (Cms+slope*xl)**2)

print(Cms)
print(samples[0:5,:])
f.close()
import fano_calc as fc
#recall defaults (filename='test.h5', 
#det='GGA3',band='ER',F=0.00001,V=4.0,alpha=(1/10000.0),aH=0.035,Erv=None,sigv=None,erase=False)
E,sig = fc.RWCalc(filename='data/res_calc.h5',band='NR',alpha=(1/18.0),F=999)
#Ebase,sigbase = fc.RWCalc(filename='data/res_calc.h5',band='NR',alpha=(1/18.0),F=000)

import EdwRes as er
#instead I should use Amy's Fit
f = er.get_sig_nuc_func_alt(1.3,1.5,3.1,0.4,2.7,pars={'V':scale*4.0,'eps_eV':3.0,'a':a_yield,'b':b_yield}, \
                            aH=aH,C=0,m=0)
Ebase = np.arange(0.1,200,0.1)
sigbase = f(Ebase)
# import data from Edelweiss
import pandas as pds
res_data = pds.read_csv("data/edelweiss_NRwidth_GGA3_data.txt", skiprows=1, \
                       names=['E_recoil', 'sig_NR', 'E_recoil_err', 'sig_NR_err'], \
                       delim_whitespace=True)

resER_data = pds.read_csv("data/edelweiss_ERwidth_GGA3_data.txt", skiprows=1, \
                         names=['E_recoil', 'sig_ER', 'sig_ER_err'], \
                         delim_whitespace=True)

resER_data = resER_data.sort_values(by='E_recoil')

print (res_data.head(10))
E_recoil = res_data["E_recoil"]
sig_NR = res_data["sig_NR"]
sig_NR_err = res_data['sig_NR_err']
E_recoil_ER = resER_data["E_recoil"]
sig_ER = resER_data["sig_ER"]
sig_ER_err = resER_data['sig_ER_err']

#order those things
idx = np.argsort(E_recoil)
E_recoil = E_recoil[idx]
sig_NR = sig_NR[idx]
sig_NR_err = sig_NR_err[idx]

#make a selectable version of the minimum index
minidx = 2
#set up a 1d plot
fig,axes = plt.subplots(1,1,figsize=(9.0,8.0),sharex=True)
ax1 = axes

for aH_em, Cms_em, m_em, scale_em, A_em, B_em in samples[np.random.randint(len(samples), size=100)]:
    ax1.plot(xl, np.sqrt(Sigss**2+(Cms_em+m_em*xl)**2), color="orange", alpha=0.1)

ax1.plot(xl,upvec,color='r',linestyle="--",linewidth=2, \
         label='1$\sigma$ fluct.')
ax1.plot(xl,dnvec,color='r',linestyle="--",linewidth=2, \
         label='')


ax1.plot(xl,np.sqrt(Sigss**2+(Cms+xl*slope)**2),color='g',linestyle="-",linewidth=3, \
         label='(C$_0$={:01.3}; m={:01.2E})'.format(Cms,slope))

ax1.plot(E,sig,color='b',linestyle=":",linewidth=2, \
        label = 'Extracted Fano')

ax1.errorbar(E_recoil[minidx::],sig_NR[minidx::], yerr=sig_NR_err[minidx::], \
         color='k', marker='o', markersize=4,linestyle='none',label='NR Edw. Measurement', linewidth=2)

ymin = 0.04
ymax = 0.1



ax1.set_yscale('linear')
#ax1.set_yscale('log')
ax1.set_xlim(10, 200) 
ax1.set_ylim(ymin,ymax)
ax1.set_xlabel(r'recoil energy [keV]',**axis_font)
ax1.set_ylabel('ionization yield width',**axis_font)
ax1.grid(True)
ax1.yaxis.grid(True,which='minor',linestyle='--')
ax1.legend(loc=1,prop={'size':22})
#ax1.legend(bbox_to_anchor=(1.04,1),borderaxespad=0,prop={'size':22})

for axis in ['top','bottom','left','right']:
  ax1.spines[axis].set_linewidth(2)

plt.tight_layout()
#plt.savefig('figures/figure.png')
plt.show()<jupyter_output><empty_output><jupyter_text>Binning Systematicenergy bin|interval (keV):-|:-0|\[5,10)1|\[10,20)2|\[20,30)3|\[30,40)4|\[40,50)5|\[50,70)6|\[70,150)<jupyter_code>bins = [5, 10, 20, 30, 40, 50, 70,150]
bins = np.asarray(bins)<jupyter_output><empty_output><jupyter_text>In order to evaluate the systematic we can approximate that the distribution in ionization yield at each energy is normal. That is it follows the probability distribution function:\begin{equation}P(Q|E_r) = \frac{1}{\sqrt{2\pi\sigma_Q(E_r)^2}} \exp{\left[\frac{(Q-\bar{Q})^2}{2\sigma_Q(E_r)^2}\right]},\end{equation}Is valid for a fixed energy, $E_r$. Using the law of conditional probability we can then simply obtain the joint distribution $P(Q,E_r)$ and integrate over the recoil energy bin. We get the following:\begin{equation}P_i(Q) = \frac{\int_{a_i}^{b_i} dE_r P(Q|E_r)P(E_r)}{N},\end{equation}where N is a normalization factor, $P_i(Q)$ is the yield probability distribution over the bin $i$, and $P(E_r)$ is the probability of getting a given recoil energy (approximated throughout this analysis as an exponential with decay coefficient $\alpha$ = (1/18.0) keV$^{-1}$).With this function the _measured_ yield width over a finite bin can be calculated quite simply:\begin{equation}\begin{aligned}\langle Q \rangle &= \int_{-\infty}^{\infty} dQ Q P(Q) \\\langle Q^2 \rangle &= \int_{-\infty}^{\infty} dQ Q^2 P(Q) \\\sigma_{Qi} &= \sqrt{\langle Q^2 \rangle - \langle Q \rangle^2}\end{aligned}\end{equation}<jupyter_code>#let's find the corrections for the best fit points
import binning_sys as bsys


print(bins)
delta = lambda C,m: bsys.getBinSys(C,m,bins=bins)[1]
deltav = np.vectorize(delta)


del110 = deltav(Cms,slope)
del110 = np.array(del110)
print(del110)
#print(np.shape(deltav(Candm[:,0],Candm[:,1])))
#del110 = deltav(Candm[:,0],Candm[:,1])
#print(delta(Candm[0,0],Candm[0,1]))
#let's make a callable out of the predicted single-scatter resolution
import scipy.interpolate as inter

sigQnr_c = inter.InterpolatedUnivariateSpline(Ebase, sigbase, k=3)
sigQnr_c_v = np.vectorize(sigQnr_c)
def log_likelihood(theta, x, y, yerr,sigfunc):
    Cms, m = theta
    model = np.sqrt(sigfunc(x)**2 + (Cms+m*x)**2)
    #sigma2 = yerr**2 + model**2*np.exp(2*log_f)
    sigma2 = yerr**2
    return -0.5*np.sum((y-model)**2/sigma2 + np.log(sigma2))
from scipy.optimize import minimize
np.random.seed(42)
nll = lambda *args: -log_likelihood(*args)
initial = np.array([0.0207, 4.73e-5]) + 0.001*np.random.randn(2)
soln = minimize(nll, initial, args=(E_recoil[minidx::], sig_NR[minidx::], sig_NR_err[minidx::],sigQnr_c_v))
Cms_ml, m_ml = soln.x

#lmfout = lmf.minimize(residual_const, params, args=(E_recoil[2::], sig_NR[2::], sig_NR_err[2::],sigQnr_c_v))
print("Maximum likelihood estimates:")
print("m = {0:.4f}".format(Cms_ml))
print("b = {0:.3e}".format(m_ml))

#do correction to data points
E_edw = E_recoil
sig_edw = sig_NR*del110
sig_edw_err = sig_NR_err*del110
print(bins)
print(del110)

print(E_edw)
print(sig_NR)
print(sig_edw)
print(sig_edw_err)

nll = lambda *args: -log_likelihood(*args)
initial = np.array([0.0207, 4.73e-5]) + 0.001*np.random.randn(2)
soln = minimize(nll, initial, args=(E_edw[minidx::], sig_edw[minidx::], sig_edw_err[minidx::],sigQnr_c_v))
Cms_ml, m_ml = soln.x

#lmfout = lmf.minimize(residual_const, params, args=(E_recoil[2::], sig_NR[2::], sig_NR_err[2::],sigQnr_c_v))
print("Maximum likelihood estimates:")
print("m = {0:.4f}".format(Cms_ml))
print("b = {0:.3e}".format(m_ml))
#set up a 1d plot
fig,axes = plt.subplots(1,1,figsize=(9.0,8.0),sharex=True)
ax1 = axes

for aH_em, Cms_em, m_em, scale_em, A_em, B_em in samples[np.random.randint(len(samples), size=100)]:
    ax1.plot(xl, np.sqrt(Sigss**2+(Cms_em+m_em*xl)**2), color="orange", alpha=0.1)

ax1.plot(xl,upvec,color='r',linestyle="--",linewidth=2, \
         label='1$\sigma$ fluct.')
ax1.plot(xl,dnvec,color='r',linestyle="--",linewidth=2, \
         label='')


ax1.plot(xl,np.sqrt(Sigss**2+(Cms+xl*slope)**2),color='g',linestyle="-",linewidth=3, \
         label='(C$_0$={:01.3}; m={:01.2E})'.format(Cms,slope))

ax1.plot(xl,np.sqrt(Sigss**2+(Cms_ml+xl*m_ml)**2),color='b',linestyle="--",linewidth=3, \
         label='bc (C$_0$={:01.3}; m={:01.2E})'.format(Cms_ml,m_ml))

ax1.plot(E,sig,color='b',linestyle=":",linewidth=2, \
        label = 'Extracted Fano')

ax1.errorbar(E_recoil[minidx::],sig_edw[minidx::], yerr=sig_edw_err[minidx::], \
         color='k', marker='o', markersize=4,linestyle='none',label='NR Edw. Measurement', linewidth=2)

ymin = 0.04
ymax = 0.1



ax1.set_yscale('linear')
#ax1.set_yscale('log')
ax1.set_xlim(10, 200) 
ax1.set_ylim(ymin,ymax)
ax1.set_xlabel(r'recoil energy [keV]',**axis_font)
ax1.set_ylabel('ionization yield width',**axis_font)
ax1.grid(True)
ax1.yaxis.grid(True,which='minor',linestyle='--')
ax1.legend(loc=1,prop={'size':22})
#ax1.legend(bbox_to_anchor=(1.04,1),borderaxespad=0,prop={'size':22})

for axis in ['top','bottom','left','right']:
  ax1.spines[axis].set_linewidth(2)

plt.tight_layout()
#plt.savefig('figures/figure.png')
plt.show()
import h5py
filename = 'data/test_bs_6par.h5'
#remove vars
f = h5py.File(filename,'r')

#save the results for the Edw fit
path='{}/{}/'.format('mcmc','edwdata_binsys')

bins_binsys = np.asarray(f[path+'bins_binsys'])
binsys = np.asarray(f[path+'binsys'])


#print(binsys)
Nsamp = np.shape(binsys)[0]
print(Nsamp)
f.close()
#set up a 1d plot
fig,axes = plt.subplots(1,1,figsize=(9.0,8.0),sharex=True)
ax1 = axes

for di in binsys:
    #print(di)
    soln = minimize(nll, initial, args=(E_edw[minidx::], sig_edw[minidx::]*di[minidx::],  \
                                        sig_edw_err[minidx::]*di[minidx::],sigQnr_c_v))
    Cms_mli, m_mli = soln.x
    #print(Cms_mli,m_mli)
    ax1.plot(xl, np.sqrt(Sigss**2+(Cms_mli+m_mli*xl)**2), color="blue", alpha=0.1)

ax1.plot(xl,upvec,color='r',linestyle="--",linewidth=2, \
         label='1$\sigma$ fluct.')
ax1.plot(xl,dnvec,color='r',linestyle="--",linewidth=2, \
         label='')


ax1.plot(xl,np.sqrt(Sigss**2+(Cms+xl*slope)**2),color='g',linestyle="-",linewidth=3, \
         label='(C$_0$={:01.3}; m={:01.2E})'.format(Cms,slope))

ax1.plot(xl,np.sqrt(Sigss**2+(Cms_ml+xl*m_ml)**2),color='b',linestyle="--",linewidth=3, \
         label='bc (C$_0$={:01.3}; m={:01.2E})'.format(Cms_ml,m_ml))

ax1.plot(E,sig,color='b',linestyle=":",linewidth=2, \
        label = 'Extracted Fano')

ax1.errorbar(E_recoil[minidx::],sig_edw[minidx::], yerr=sig_edw_err[minidx::], \
         color='k', marker='o', markersize=4,linestyle='none',label='NR Edw. Measurement', linewidth=2)

ymin = 0.04
ymax = 0.1



ax1.set_yscale('linear')
#ax1.set_yscale('log')
ax1.set_xlim(10, 200) 
ax1.set_ylim(ymin,ymax)
ax1.set_xlabel(r'recoil energy [keV]',**axis_font)
ax1.set_ylabel('ionization yield width',**axis_font)
ax1.grid(True)
ax1.yaxis.grid(True,which='minor',linestyle='--')
ax1.legend(loc=1,prop={'size':22})
#ax1.legend(bbox_to_anchor=(1.04,1),borderaxespad=0,prop={'size':22})

for axis in ['top','bottom','left','right']:
  ax1.spines[axis].set_linewidth(2)

plt.tight_layout()
#plt.savefig('figures/figure.png')
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.2009233774305247, 'fraction_numerical': 0.041019266625233065, 'mean_word_length': 5.744910179640718, 'pattern_counts': {'":': 3, '<': 13, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '524773', 'license': 'MIT', 'n_tokens_mistral': 5027, 'n_tokens_neox': 4778, 'n_words': 747, 'path': 'analysis_notebooks/binning_systematic_Final.ipynb', 'repo_name': 'villano-lab/nrFano_paper2019', 'size': '11263'}