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