<jupyter_start><jupyter_text>Comparison for WIMpy_NREFT This is an example file to go with the `WIMpy_NREFT` code, available [here](https://github.com/bradkav/WIMpy_NREFT). The main functions are in `DMUtils.py`, which is loaded like this:<jupyter_code>import DMUtils as DMU
import DMUtils_old as DMUold

#We'll also import some useful libraries
import numpy as np
import matplotlib.pyplot as pl

from tqdm import tqdm<jupyter_output><empty_output><jupyter_text>Loading nuclear form factors We load the NREFT form factor coefficients for the nucleus we're interested in (say Iodine for now) as follows. The form factor coefficients (`FF_I` in this case) are required to calculate the recoil rate.<jupyter_code>A_I = 19 #Nucleon number
Z_I = 9 #Atomic number
root = "FormFactors" #Folder where the form factor files are stored
FF_I = DMU.LoadFormFactors(root, A_I, Z_I)<jupyter_output>Loading Form Factor for (A, Z) = (19, 9)...<jupyter_text>Differential recoil rate Before 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 numbers of arrays starts from zero in python):<jupyter_code>cp[0] = 1.0
cn[0] = 1.0
print("Couplings to protons [GeV^{-2}]:", cp)
print("Couplings to neutrons [GeV^{-2}]:", cn)<jupyter_output>('Couplings to protons [GeV^{-2}]:', array([ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.]))
('Couplings to neutrons [GeV^{-2}]:', array([ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.]))<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:<jupyter_code>m_x = 100 #DM mass in GeV

E_list = np.linspace(0, 100, 60)
R_list_SI = E_list*0.0
R_list_SI_old = E_list*0.0

for j in tqdm(range(int(100))):
    R_list_SI = (4*np.pi/(2*(0.5)+1))*DMU.dRdE_NREFT(E_list, A_I, m_x, cp, cn, "Fluorine") 

for i, E in enumerate(E_list):
    R_list_SI_old[i] = DMUold.dRdE_NREFT_old(E, A_I, m_x, cp, cn, FF_I)<jupyter_output>100%|██████████| 100/100 [00:00<00:00, 774.23it/s]<jupyter_text>Now let's plot the result, which will be in units of keV$^{-1}$ kg$^{-1}$ day$^{-1}$:<jupyter_code>#print R_list_SI_old

pl.figure()
pl.plot(E_list, R_list_SI)
pl.plot(E_list, R_list_SI_old, 'g--', lw=2.0)
pl.xlabel(r'$E_R$ [keV]')
pl.ylabel(r'$\mathrm{d}R/\mathrm{d}E_R$ [keV$^{-1}$ kg$^{-1}$ day$^{-1}$]')
pl.title('Standard SI interactions')
pl.show()<jupyter_output><empty_output><jupyter_text>Now 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}]:', array([ 0.        , -0.37659593,  1.71680264,  1.56974372,  0.37376856,
        1.63652149, -0.68606561,  0.83867336, -1.98127741,  1.12393638,
       -0.8585    ]))
('Couplings to neutrons [GeV^{-2}]:', array([ 0.        ,  0.03986585, -0.19475956,  0.87252854, -0.54769899,
       -0.72719728,  1.14118449,  1.41142361,  0.69142929, -1.4057233 ,
        0.18440249]))<jupyter_text>Then let's calculate and plot again:<jupyter_code>E_list = np.linspace(0, 100, 1000)
R_list_random = E_list*0.0
R_list_random_old = E_list*0.0
for i, E in enumerate(E_list):
    R_list_random[i] = (4*np.pi/(2*(0.5)+1))*DMU.dRdE_NREFT(E, A_I, m_x, cp_random, cn_random,"Fluorine") 
    R_list_random_old[i] = DMUold.dRdE_NREFT_old(E, A_I, m_x, cp_random, cn_random, FF_I) 
pl.figure()
pl.plot(E_list, R_list_random)
pl.plot(E_list, R_list_random_old,'g--', lw=2.0)
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.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, A_I, m_x, cp_random, cn_random, FF_I, eff=efficiency))<jupyter_output>Number of signal events per day per kg of target mass:
1.95125821474e-05
{'fraction_non_alphanumeric': 0.14689482791770814, 'fraction_numerical': 0.06517977312055374, 'mean_word_length': 4.938356164383562, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 32, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '388459', 'license': 'MIT', 'n_tokens_mistral': 2271, 'n_tokens_neox': 2021, 'n_words': 575, 'path': 'Examples/.ipynb_checkpoints/NREFT_comparison-checkpoint.ipynb', 'repo_name': 'bradkav/WIMpy_NREFT', 'size': '5201'}