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()