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

a_pion_p,a_pion_n,a_eta_p,a_eta_n,b_p,b_n=0.2166300422469397, -0.21663004224694038, 0.07607523965672636, 0.07605541831688267, -0.41945356821526036, 0.013814848551653033
#constants of proportionality between the non-relativistic couplings a_pion, a_meson and b defined in Eq. (2.10) of 1707.06998 and the effective coupling of the relativistic Hamiltonian for the operator Q^(7)_3 defined in the same paper. Valid at leading order, calculated with directDM, arXiv:1708.02678

a_pion_0=a_pion_p+a_pion_n
a_pion_1=a_pion_p-a_pion_n
a_eta_0=a_eta_p+a_eta_n
a_eta_1=a_eta_p-a_eta_n
b_0=b_p+b_n
b_1=b_p-b_n

m_pion=134.9766*1e-3 # GeV
m_eta=547.862*1e-3 # GeV

c73_model=WD.eft_hamiltonian('c73',{10: lambda c73: c73*np.array([b_0,b_1]), (10,'pion'): lambda c73,q: c73*np.array([a_pion_0,a_pion_1])*q**2/(m_pion**2+q**2), (10,'eta'): lambda c73,q: c73*np.array([a_eta_0,a_eta_1])*q**2/(m_eta**2+q**2)})
#effective hamiltonian including meson-pole effect

vmin,delta_eta0=WD.streamed_halo_function()

mchi_vec=np.logspace(-1,3,100)

mchi_vec,c73_squared=WD.mchi_vs_exclusion(WD.PICO60_2019,c73_model,vmin,delta_eta0,mchi_vec=mchi_vec)
#returns two arrays of m_chi and lambda_tilde sixth

c73_lambda_tilde=c73_squared**(-1/6)
#to get lambda_tilde
mchi_vec_new=mchi_vec
c73_lambda_tilde_new=c73_lambda_tilde
pl.clf()
pl.plot(mchi_vec_new,c73_lambda_tilde_new,label='90% C.L.')
pl.xlabel('$m_{\chi}$  (GeV)',fontsize=10)
pl.ylabel('$\~\Lambda$  (GeV)',fontsize=10)
pl.xscale('log')
pl.yscale('log')

pl.legend(title='$Q^{(7)}_3$',loc=4,fontsize=15,title_fontsize=15,frameon=False)

pl.show()