import WimPyDD as WD
import numpy as np
import matplotlib.pyplot as pl

pl.clf()
vmin, delta_eta0=WD.streamed_halo_function()


# defines c_1^tau=[c_1^0, c_1^1] in terms of the WIMP-nucleon cross section sigma_N=(c_N)^2*mu_chi_N^2/pi, with c_N=c_1^p=c_1^n, c_1^0=c_1^p+c_1^n,
# c_1^1=c_1^p-c_1^n (for r=c_1^n/c_1^p=1)
def c_tau_SI(sigma_N,mchi,cn_over_cp=1):
    hbarc2=0.389e-27 #(hbar*c)^2 in GeV^2 * cm^2
    mn=0.931
    mu=mchi*mn/(mchi+mn)
    return np.sqrt(np.pi*sigma_N/hbarc2)/mu*np.array([1+cn_over_cp,1-cn_over_cp])

SI=WD.eft_hamiltonian('Spin-independent',{1: c_tau_SI})

# uses built-in XENON_1T_2018 experiment 
mchi,sigma=WD.mchi_vs_exclusion(WD.XENON_1T_2018, SI,vmin, delta_eta0)
pl.plot(mchi,sigma,':r',label='WimPyDD')


# published exclusion plot, from arXiv:1805.12562
mchi_xenon1t=np.array([  6.38084285,   6.82568656,   7.38401692,   8.16949332,
         9.14061906,  10.1129544 ,  11.44291272,  13.39150131,
        16.76448605,  19.61927374,  23.2195425 ,  27.48048481,
        31.80087168,  36.3894599 ,  45.04623899,  54.52370581,
        74.67422632, 102.2718466 , 160.27952577, 262.73153013,
       376.36545636, 539.14715404, 705.96063944])

sigma_xenon1t=np.array([1.34339933e-44, 7.44380301e-45, 3.66524124e-45, 1.80472177e-45,
       9.15247311e-46, 5.37983840e-46, 2.89426612e-46, 1.65176674e-46,
       8.37677640e-47, 6.23550734e-47, 4.78065253e-47, 4.24820170e-47,
       4.24820170e-47, 4.37547938e-47, 4.92388263e-47, 5.54102033e-47,
       7.22727132e-47, 9.42668455e-47, 1.42510267e-46, 2.28546386e-46,
       3.25702066e-46, 4.64158883e-46, 5.87801607e-46])

pl.plot(mchi_xenon1t,sigma_xenon1t,'-k',label='Published')
pl.title('SI, XENON1T exclusion plot')
pl.xscale('log')
pl.yscale('log')
pl.xlabel('$m_\chi$ (GeV)')
pl.ylabel('$\sigma_p$ (cm$^2$)')
pl.legend()
pl.show()