<jupyter_start><jupyter_text>Examples for WIMpy_NREFT Example notebook to go with the `WIMpy_NREFT` code, available [here](https://github.com/bradkav/WIMpy_NREFT). This notebook let's you calculate and plot recoil spectra $\mathrm{d}R/\mathrm{d}E_R$ for a range of interactions, relevant for the direct detection of Dark Matter. In particular, we focus on the non-relativistic effective field theory (NREFT) of direct detection. Don't know how to use a Jupyter notebook? [Try this link](http://jupyter-notebook.readthedocs.io/en/stable/examples/Notebook/Running%20Code.html), or just hit `Shift+Enter` over and over again to run the cells... ReferencesFor a review of the general 'direct detection' formalism, see [arXiv:1002.1912](https://arxiv.org/abs/1002.1912).For more information about the NREFT, see e.g. [arXiv:1203.3542](https://arxiv.org/abs/1203.3542) and [arXiv:1307.5955](https://arxiv.org/abs/1307.5955). Contents * [Getting started with WIMpy](Getting-started-with-WIMpy) * [Standard spin-independent interactions](Spin-independent-interactions) * [Comparing spin-indepedent and spin-dependent interactions](Comparing-$\mathcal{O}_1$-and-$\mathcal{O}_4$) * [Non-relativistic effective field theory (NREFT) spectra](Random-NREFT-spectra) * [Specific interactions (millicharge, magnetic dipole, etc.)](Specific-interactions) * [Using different nuclear targets](Different-targets)  ------------  Getting started with WIMpy The main functions are in `DMUtils.py`, which is loaded like this:<jupyter_code>from WIMpy import DMUtils as DMU<jupyter_output><empty_output><jupyter_text>We'll also load a bunch of useful libraries (numpy, matplotlib, etc.):<jupyter_code>%matplotlib inline
#We'll also import some useful libraries
import numpy as np
import matplotlib.pyplot as pl

import matplotlib as mpl

font = {'family' : 'sans-serif',
        'size'   : 16}

mpl.rcParams['xtick.major.size'] = 5
mpl.rcParams['xtick.major.width'] = 1
mpl.rcParams['xtick.minor.size'] = 3
mpl.rcParams['xtick.minor.width'] = 1
mpl.rcParams['ytick.major.size'] = 5
mpl.rcParams['ytick.major.width'] = 1
mpl.rcParams['ytick.minor.size'] = 3
mpl.rcParams['ytick.minor.width'] = 1
mpl.rc('font', **font)

mpl.rcParams['xtick.direction'] = 'in'
mpl.rcParams['ytick.direction'] = 'in'
mpl.rcParams['xtick.top'] = True
mpl.rcParams['ytick.right'] = True

from tqdm import tqdm
from scipy.integrate import quad<jupyter_output><empty_output><jupyter_text>--------------- Spin-independent interactions  Differential recoil rateBefore we calculate the differential recoil rate $\mathrm{d}R/\mathrm{d}E_R$, we have to specify the couplings. We do this by defining two vectors of couplings: 1 for the couplings to protons and 1 for the couplings to neutrons. Each vector should have 11 elements (one entry for each of the 11 NREFT operators - note that $\mathcal{O}_2$ is set to zero by default, but we keep 11 entries so that the numbering matches):<jupyter_code>#Initialise the empty vectors
cp = np.zeros(11)
cn = np.zeros(11)<jupyter_output><empty_output><jupyter_text>Let's start with the standard spin-independent interaction (so we want to set the first element $\mathcal{O}_1$, noting that the numbering of arrays starts from zero in python):<jupyter_code>cp[0] = 1e-9
cn[0] = 1e-9
print("Couplings to protons [GeV^{-2}]:", cp)
print("Couplings to neutrons [GeV^{-2}]:", cn)<jupyter_output>Couplings to protons [GeV^{-2}]: [1.e-09 0.e+00 0.e+00 0.e+00 0.e+00 0.e+00 0.e+00 0.e+00 0.e+00 0.e+00
 0.e+00]
Couplings to neutrons [GeV^{-2}]: [1.e-09 0.e+00 0.e+00 0.e+00 0.e+00 0.e+00 0.e+00 0.e+00 0.e+00 0.e+00
 0.e+00]<jupyter_text>**Note that the couplings in the NREFT are *dimensionful* and should have units of GeV$^{-2}$**. Now let's calculate the differential spectrum for a DM particle of mass 100 GeV (using the NREFT framework and the standard formalism for spin-independent interactions - to make sure everything matches):<jupyter_code>m_x = 100 #DM mass in GeV

#Number of protons and neutrons in Xenon
A_Xe = 131.
N_p_Xe = 54.
N_n_Xe = 131 - N_p_Xe

#Convert the coupling to a standard SI cross section
sig = DMU.coupling_to_xsec(cp[0], m_x)

E_list = np.logspace(-3, 2,1000)

R_SI = DMU.dRdE_standard(E_list, N_p_Xe, N_n_Xe, m_x, sig)
R_SI_NREFT = DMU.dRdE_NREFT(E_list, m_x, cp, cn, target="Xe131")<jupyter_output><empty_output><jupyter_text>Now let's plot the result, which will be in units of keV$^{-1}$ kg$^{-1}$ day$^{-1}$:<jupyter_code>pl.figure(figsize=(6,4))

pl.plot(E_list, R_SI, label="Standard SI")
pl.plot(E_list, R_SI_NREFT, '--', lw=2.0, label=r'$\mathcal{O}_1$')
pl.xlabel(r'$E_R$ [keV]')
pl.ylabel(r'$\mathrm{d}R/\mathrm{d}E_R$ [keV$^{-1}$ kg$^{-1}$ day$^{-1}$]')

pl.show()
E_list = np.logspace(-1, 3, 1000)
dRdE_list = DMU.dRdE_standard(E_list, 131.0, 0.0, 100.0, 1e-46, vlag=230.0, sigmav=156.0,vesc=544.0)


np.savetxt("../tests/dRdE_Xenon.txt", list(zip(E_list, dRdE_list)))

pl.figure()

pl.loglog(E_list, dRdE_list)

pl.show()<jupyter_output><empty_output><jupyter_text>------------------ Comparing $\mathcal{O}_1$ and $\mathcal{O}_4$Now we take a look at the spectra for $\mathcal{O}_1$ and $\mathcal{O}_4$ (standard spin-independent and spin-dependent interactions):<jupyter_code>#Spin-independent coupling
cp_SI = np.zeros(11)
cp_SI[0] = 2.5e-5

cn_SI = 1.0*cp_SI

#Spin-dependent coupling
cp_SD = np.zeros(11)
cp_SD[3] = 1.0

cn_SD = np.zeros(11)
cn_SD[3] = 2.5e-2

m_x = 100 #DM mass in GeV

E_list = np.logspace(-3, 2,1000)

pl.figure()

pl.plot(E_list, 1e-3*DMU.dRdE_NREFT(E_list, m_x, cp_SI, cn_SI, "Xe131"),'k--' , lw=1.5,label ='SI')
pl.plot(E_list, 1e-3*DMU.dRdE_NREFT(E_list, m_x, cp_SD, 0.0*cp_SD, "Xe131"),'r-' , lw=1.5,label ='SD (p)')
pl.plot(E_list, 1e-3*DMU.dRdE_NREFT(E_list, m_x, 0.0*cn_SD, cn_SD, "Xe131"),'b-' , lw=1.5,label ='SD (n)')
pl.xlabel(r'$E_R$ [keV]')
#pl.ylabel(r'$\mathrm{d}R/\mathrm{d}E_R$ [keV$^{-1}$ kg$^{-1}$ day$^{-1}$]')
#pl.title('Random NREFT interactions')
pl.ylabel(r'$\mathrm{d}R/\mathrm{d}E_R$ [arb. units]')


pl.legend(loc='best',fancybox = True)
pl.savefig("../plots/SIvsSD.pdf", bbox_inches="tight")
pl.show()<jupyter_output><empty_output><jupyter_text>------------ Random NREFT spectraNow let's try something more exotic. Let's just take some random values for all the couplings to protons and neutrons (setting the couplings to$\mathcal{O}_1$ to zero):<jupyter_code>cp_random = np.random.randn(11)
cn_random = np.random.randn(11)
cp_random[0] = 0
cn_random[0] = 0
print("Couplings to protons [GeV^{-2}]:", cp_random)
print("Couplings to neutrons [GeV^{-2}]:", cn_random)<jupyter_output>Couplings to protons [GeV^{-2}]: [ 0.          0.08471479  0.72814372 -0.43902207  1.32343557  0.92118775
  0.0719406   0.47921031 -1.16790334 -1.40743246 -1.77060201]
Couplings to neutrons [GeV^{-2}]: [ 0.          0.48246571 -0.12894153 -2.06375125  0.18499707  0.68374067
 -0.10273287 -0.01755086  0.4826994   1.1853268  -0.43329687]<jupyter_text>Now let's plot a bunch of different random spectra:<jupyter_code>m_x = 100 #DM mass in GeV

E_list = np.logspace(-3, 2,100)
R_random = np.zeros((len(E_list), 10))

pl.figure(figsize=(7,5))

for i in range(10):
    cp_random = np.random.randn(11)
    cn_random = np.random.randn(11)
    cp_random[0] = 0
    cn_random[0] = 0
    
    R_random[:,i] = DMU.dRdE_NREFT(E_list, m_x, cp_random, cn_random, "Xe131") 
    pl.loglog(E_list, R_random)
    
pl.xlabel(r'$E_R$ [keV]')
pl.ylabel(r'$\mathrm{d}R/\mathrm{d}E_R$ [keV$^{-1}$ kg$^{-1}$ day$^{-1}$]')
pl.show()<jupyter_output><empty_output><jupyter_text>Event Rates You can also calculate the total number of events, using `Nevents_NREFT`. You have to specify a minimum and maximum energy, and you can also specify an efficiency function (i.e. the efficiency of the detector as a function of recoil energy) which multiplies the recoil rate:<jupyter_code>Emin = 5 #keV
Emax = 50 #keV

#Just some made-up efficiency function
def efficiency(E_R):
    return 0.5/(1 + np.exp(-E_R))
print("Number of signal events per day per kg of target mass:")
print(DMU.Nevents_NREFT(Emin, Emax,  m_x, cp_random, cn_random, eff=efficiency,target="Xe131"))<jupyter_output>Number of signal events per day per kg of target mass:
5521468766.542366<jupyter_text>----------------- Specific interactionsLet's look at a few specific interactions. In particular, we look at **millicharged DM**, with millicharge $\epsilon$:$$ \mathcal{L} \supset \epsilon \overline{\chi}\gamma_\mu \chi\,A^\mu\,,$$Dark Matter which has a **magnetic dipole moment** $\mu_\chi$:$$ \mathcal{L} \supset\frac{\mu_\chi}{2}\overline{\chi}\sigma^{\mu \nu} \chi \,F_{\mu\nu}\,,$$and Dark Matter which has an **anapole moment** $c_A$:$$ \mathcal{L} \supset c_A \overline{\chi}\gamma^\mu\gamma^5\chi\,\partial^\nu F_{\mu\nu}\,.$$<jupyter_code>#Number of protons and neutrons in Xenon
A_Xe = 131.
N_p_Xe = 54.
N_n_Xe = 131 - N_p_Xe

E_list = np.logspace(-3, 2,100)

mu_x = 1e-8 # in units of the Bohr Magneton
eps = 1e-10 #in units of the electron charge
c_A = 1e-5 #in units of GeV^-2

#DM mass in GeV
m_x = 100.0

#SI cross section
sig = 1e-46

pl.figure(figsize=(7,5))
pl.loglog(E_list, DMU.dRdE_standard(E_list, N_p_Xe, N_n_Xe, m_x, sig), lw=1.5, label='Spin-independent')
pl.loglog(E_list, DMU.dRdE_millicharge(E_list, m_x, eps, "Xe131"), lw=1.5, label='Millicharge DM')
pl.loglog(E_list, DMU.dRdE_magnetic(E_list, m_x, mu_x, "Xe131"), lw=1.5, label='Magnetic DM')
pl.loglog(E_list, DMU.dRdE_anapole(E_list, m_x, c_A, "Xe131"), lw=1.5, label='Anapole DM')

pl.xlabel(r'$E_R$ [keV]')
pl.ylabel(r'$\mathrm{d}R/\mathrm{d}E_R$ [keV$^{-1}$ kg$^{-1}$ day$^{-1}$]')
pl.legend(loc = 'best')
pl.show()<jupyter_output><empty_output><jupyter_text>---------------- Different targetsYou can specify which nucleus you want to calculate the scattering rate for using the `target` keyword. So far, we have coded in the following targets:<jupyter_code>print(DMU.target_list)<jupyter_output>['Xe128' 'Xe129' 'Xe130' 'Xe131' 'Xe132' 'Xe134' 'Xe136' 'Ar40' 'C12'
 'Ge70' 'Ge72' 'Ge73' 'Ge74' 'Ge76' 'Fluorine' 'Iodine' 'Na23']<jupyter_text>So for example, to see the number of signal events in Carbon 12 for a given set of couplings, you just do:<jupyter_code>c_SI = np.zeros(11)
c_SI[0] = 1.0

print("Number of signal events per day per kg of target mass (Carbon-12):")
print(DMU.Nevents_NREFT(E_min=1.0, E_max=100.0,  m_x=100.0, cp=c_SI, cn=c_SI, eff=None,target="C12"))
Np = 37
Nn = 0

from scipy.integrate import quad
quad(lambda x: x*DMU.dRdE_standard(x, Np, Nn, 50, 5e-46), 1e-1, 1e3)[0]/quad(lambda x: DMU.dRdE_standard(x, Np, Nn, 50, 5e-46), 1e-1, 1e3)[0]<jupyter_output><empty_output><jupyter_text>---------------------------------- Relativistic kinematicsThis is not part of the full module yet - it's only just being tested:<jupyter_code>import sys
sys.path.append("../WIMpy/")
import DMUtils_relativistic as DMU_rel
m_x = 0.1
m_A = 40

ER_list = np.logspace(-3, 9)

pl.figure(figsize=(7,5))

pl.loglog(1e-3*ER_list, DMU_rel.vmin_rel(ER_list, m_x, m_A), label='Relativistic')
pl.loglog(1e-3*ER_list, DMU.vmin(ER_list, m_x, m_A), linestyle='--', label='Non-relativistic')

pl.legend(fontsize=14, loc='best')

pl.xlabel(r'$E_R$ [eV]')
pl.ylabel(r'$v_\mathrm{min}$ [km/s]')

pl.axhline(3e5, color='grey', linestyle=':')
pl.title(r'$m_\chi = 0.1\,\mathrm{GeV}$; $A = 40$', fontsize=16)

pl.savefig("../plots/RelativisticKinematics.pdf", bbox_inches='tight')

pl.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17101220261610042, 'fraction_numerical': 0.059608462821525765, 'mean_word_length': 4.9364252214695155, 'pattern_counts': {'":': 0, '<': 53, '<?xml version=': 0, '>': 53, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '192105', 'license': 'MIT', 'n_tokens_mistral': 5103, 'n_tokens_neox': 4546, 'n_words': 1104, 'path': 'Examples/NREFT_example.ipynb', 'repo_name': 'bradkav/WIMpy_NREFT', 'size': '11391'}