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


pl.clf()
vmin,delta_eta1=WD.streamed_halo_function(yearly_modulation=True)
# delta_eta1 for yearly modulation calculated setting yearly_modulation=True
#

# define Wilson coefficient c9 in terms of the reference cross section
# sigma_ref=c9^2*mu_chi_n^2/pi, with mu_chi_n the WIMP-nucleon reduced mass, # and of r=c9^n/c9^p (WIMP-neutron over WIMP-proton ratio)
# sigma_ref in cm^2 and c9 in GeV^-2.

def c9_tau(mchi,cross_section,r):
    hbarc2=0.389e-27
    mn=0.931
    mu=mchi*mn/(mchi+mn)
    cp=(np.pi*cross_section/(mu**2*hbarc2))**(1./2.)
    return cp*np.array([1.+r,1.-r])

# instantiates the effective Hamiltonian containing only O9.
c9=WD.eft_hamiltonian('c9',{9: c9_tau})

# best-fit values from Table 1 of rXiv:1804.07528.
cross_section=8.29e-33
mchi=9.3
r=4.36

# load the response functions (j_chi=0.5 is the spin default value) 
WD.load_response_functions(WD.DAMA_LIBRA_2019,c9,verbose=False)

e,sm=WD.wimp_dd_rate(WD.DAMA_LIBRA_2019, c9, vmin, delta_eta1,cross_section=cross_section, mchi=mchi , r=r)

pl.plot(e,sm,'-k',linewidth=3)

# modulation amplitudes measured by DAMA-LIBRA (from Fig. 11 of NUCL. PHYS. AT. ENERGY 19 (2018) 307-325)
experimental_s1=np.array([0.0232    , 0.0164    , 0.0178    , 0.019     , 0.0178,0.0109, 0.011, 0.004, 0.0065, 0.0066,0.0009,0.0016,0.0007,0.0016, 0.00039922])

errors_on_s1=np.array([0.0052, 0.0043, 0.0028    , 0.0029    , 0.0028,0.0025    , 0.0022    , 0.002     , 0.002, 0.0019,0.0018    , 0.0018    , 0.0018    , 0.0018    , 0.00046397])

errors_on_energy=np.array([0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25,0.25, 0.25, 0.25, 4.  ])
pl.errorbar(e,experimental_s1,yerr=errors_on_s1,xerr=errors_on_energy,fmt='none',ecolor='red')

pl.xlabel('E$^{\prime}$ (keV)')
pl.ylabel('$S^{1}_{[E_1^{\prime},E_2^{\prime}]}$ (events/kg/day/keV)')

pl.show()