move sddm.tex to doc/

4 files changed, 12 insertions, 431 deletions

diff --git a/Makefile b/Makefile
index 2082cdb..24a8580 100644
--- a/Makefile
+++ b/Makefile
@@ -7,16 +7,7 @@ default: all
 	cd src && "$(MAKE)" $@
 	cd utils && "$(MAKE)" $@
 
-all: sddm.pdf
+doc:
+	cd doc && "$(MAKE)"
 
-all: sddm.pdf
-
-calculate_limits: calculate_limits.c
-
-sddm.pdf: sddm.tex
-	pdflatex $^
-
-clean:
-	rm calculate_limits
-
-.PHONY: clean sddm.pdf
+.PHONY: doc
diff --git a/calculate_limits.py b/calculate_limits.py
deleted file mode 100755
index f297a69..0000000
--- a/calculate_limits.py
+++ /dev/null
@@ -1,419 +0,0 @@
-#!/usr/bin/env python3
-import numpy as np
-from scipy.integrate import quad
-from scipy.stats import expon
-from scipy.special import spherical_jn, erf
-from numpy import pi
-from functools import lru_cache
-
-# on retina screens, the default plots are way too small
-# by using Qt5 and setting QT_AUTO_SCREEN_SCALE_FACTOR=1
-# Qt5 will scale everything using the dpi in ~/.Xresources
-import matplotlib
-matplotlib.use("Qt5Agg")
-
-# speed of light (exact)
-SPEED_OF_LIGHT = 299792458 # m/s
-
-# depth of the SNO detector (m)
-# currently just converted 6800 feet -> meters
-h = 2072.0
-
-# radius of earth (m)
-# from the "Earth Fact Sheet" at https://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html
-# we use the volumetric mean radius here
-R = 6.371e6
-
-# Approximate dark matter velocity in m/s. The true distribution is expected to
-# be a Maxwell Boltzmann distribution which is modulated annually by the
-# earth's rotation around the sun, but we just assume a single constant
-# velocity here. From Lewin and Smith Appendix B
-DM_VELOCITY = 230e3
-
-# Approximate earth velocity in the galactic rest frame (?)
-# from Lewin and Smith Equation 3.6
-EARTH_VELOCITY = 244e3
-
-# Approximate dark matter density in GeV/m^3. From Tom Caldwell's thesis.
-# from Lewin and Smith page 91
-DM_DENSITY = 0.4e6
-
-# Number density of scatterers in the Earth.
-#
-# FIXME: Currently just set to the number density of atoms in water. Need to
-# update this for rock, and in fact this will change near the detector since
-# there is water outside the AV.
-DENSITY_WATER = 1e3 # In kg/m^3
-
-# mean density of norite rock which is the rock surrounding SNO
-# probably conservative
-NORITE_DENSITY = 3e3 # In kg/m^3
-
-# atomic masses for various elements
-# from https://www.angelo.edu/faculty/kboudrea/periodic/structure_mass.htm
-element_mass = {
-  'H':1.0,
-  'C':12.011,
-  'O':15.9994,
-  'Na':22.98977,
-  'Mg':24.305,
-  'Al':26.98154,
-  'Si':28.0855,
-  'K':39.0983,
-  'Ca':40.08,
-  'Mn':54.9380,
-  'Fe':55.847,
-  'Ti':47.90
-}
-
-# composition and mean density of norite rock which is the rock surrounding SNO
-# the composition is from Table 3.2 in the SNOLAB User's handbook
-water = {'composition':
- {'H':20,
-  'O':80},
-  'density':1e3}
-
-# composition and mean density of the mantle
-# from www.knowledgedoor.com/2/elements_handbook/element_abundances_in_the_earth_s_mantle.html
-# density from hyperphysics.phy-astr.gsu.edu/hbase/Geophys/earthstruct.html
-mantle = {'composition':
- {'O':44.33,
-  'Mg':22.17,
-  'Si':21.22,
-  'Fe':6.3},
-  'density':4.4e3}
-
-# composition and mean density of norite rock which is the rock surrounding SNO
-# the composition is from Table 3.2 in the SNOLAB User's handbook
-norite = {'composition':
- {'H':0.15,
-  'C':0.04,
-  'O':46.0,
-  'Na':2.2,
-  'Mg':3.3,
-  'Al':9.0,
-  'Si':26.2,
-  'K': 1.2,
-  'Ca':5.2,
-  'Mn':0.1,
-  'Fe':6.2,
-  'Ti':0.5},
-  'density':3e3}
-
-# Fiducial volume (m)
-FIDUCIAL_RADIUS = 5
-
-# Fiducial volume (m^3)
-FIDUCIAL_VOLUME = 4*pi*FIDUCIAL_RADIUS**3/3
-
-# proton mass from the PDG (2018)
-PROTON_MASS = 0.938 # GeV
-
-# proton mass from the PDG (2018)
-ATOMIC_MASS_UNIT = 0.931 # GeV
-
-# mass of Xenon in atomic units
-XENON_MASS = 131.293
-
-# mass of Neon in atomic units
-NEON_MASS = 20.18
-
-# mass of argon in atomic units
-ARGON_MASS = 39.948
-
-# mass of germanium in atomic units
-GERMANIUM_MASS = 72.64
-
-# mass of tungsten in atomic units
-TUNGSTEN_MASS = 183.84
-
-# mass of oxygen in atomic units
-OXYGEN_MASS = 15.999
-
-# mass of silicon in atomic units
-SILICON_MASS = 28.0855
-
-# mass of iron in atomic units
-IRON_MASS = 55.845
-
-# mass of magnesium in atomic units
-MAGNESIUM_MASS = 24.305
-
-# galactic escape velocity (m/s)
-# from Tom Caldwell's thesis page 25
-ESCAPE_VELOCITY = 244e3
-
-# conversion constant from PDG
-HBARC = 197.326978812e-3 # GeV fm
-
-# Avogadros number (kg^-1)
-N0 = 6.02214085774e26
-
-def get_probability(r, l):
-    """
-    Returns the probability of a dark photon decaying in the SNO detector from
-    a dark matter decay distributed uniformly in the Earth. Assumes that the
-    depth of SNO is much larger than the dimensions of the SNO detector.
-    """
-    if r <= h:
-        theta_min = 0
-    elif r <= 2*R - h:
-        theta_min = pi - np.arccos((h**2 + r**2 - 2*R*h)/(2*r*(R-h)))
-    else:
-        return 0
-    return (1 + np.cos(theta_min))*np.exp(-r/l)/(2*l)
-
-@lru_cache()
-def get_probability2(l):
-    p, err = quad(get_probability, 0, min(2*R-h,10*l), args=(l,), epsabs=0, epsrel=1e-7, limit=1000)
-
-    return p
-
-def get_event_rate(m, cs0, l, A):
-    """
-    Returns the event rate of leptons produced from self-destructing dark
-    matter in the SNO detector for a given dark matter mass m, a cross section
-    cs, and a mediator decay length l.
-    """
-    # For now we assume the event rate is constant throughout the earth, so we
-    # are implicitly assuming that the cross section is pretty small.
-    flux = DM_VELOCITY*DM_DENSITY/m
-
-    #p, err = quad(get_probability, 0, min(2*R-h,10*l), args=(l,), epsabs=0, epsrel=1e-7, limit=1000)
-    p = get_probability2(l)
-
-    cs = get_cross_section(cs0, m, A)
-
-    # FIXME: factor of 2 because the DM particle decays into two mediators?
-    return p*cs*(N0/A)*flux*FIDUCIAL_VOLUME
-
-def get_event_rate_sno(m, cs0, l, composition):
-    """
-    Returns the event rate of leptons produced from self-destructing dark
-    matter in the SNO detector for a given dark matter mass m, a cross section
-    cs, and a mediator decay length l.
-    """
-    rate = 0.0
-
-    for element, mass_fraction in composition['composition'].items():
-        rate += mass_fraction/100*get_event_rate(m,cs0,l,element_mass[element])
-
-    return rate*composition['density']
-
-def get_nuclear_form_factor(A, e):
-    """
-    Returns the nuclear form factor for a WIMP-nucleus interaction.
-
-    From Tom Caldwell's thesis page 24.
-
-    Also used Mark Pepin's thesis page 50
-    """
-    # mass of xenon nucleus
-    mn = A*ATOMIC_MASS_UNIT
-
-    # calculate approximate size of radius
-    s = 0.9 # fm
-    a = 0.52 # fm
-    c = 1.23*A**(1/3) - 0.60 # fm
-    r1 = np.sqrt(c**2 + (7/3)*pi**2*a**2 - 5*s**2)
-
-    q = np.sqrt(2*mn*e)
-
-    if q*r1/HBARC < 1e-10:
-        return 1.0
-
-    # Helm form factor
-    # from Mark Pepin's thesis page 50
-    f = 3*spherical_jn(1,q*r1/HBARC)*np.exp(-(q*s/HBARC)**2/2)/(q*r1/HBARC)
-
-    return f
-
-def get_cross_section(cs0, m, A):
-    """
-    Returns the WIMP cross section from the target-independent WIMP-nucleon
-    cross section cs0 at zero momentum transfer.
-
-    From Tom Caldwell's thesis page 21.
-    """
-    # mass of xenon nucleus
-    mn = A*ATOMIC_MASS_UNIT
-
-    # reduced mass of the nucleus and the WIMP
-    mr = (m*mn)/(m + mn)
-
-    # reduced mass of the proton and the WIMP
-    mp = (m*PROTON_MASS)/(m + PROTON_MASS)
-
-    return cs0*A**2*mr**2/mp**2
-
-def get_differential_event_rate_xenon(e, m, cs0, A):
-    """
-    Returns the event rate of WIMP scattering in the Xenon 100T detector.
-    """
-    # mass of nucleus
-    mn = A*ATOMIC_MASS_UNIT
-
-    # reduced mass of the nucleus and the WIMP
-    mr = (m*mn)/(m + mn)
-
-    cs = get_cross_section(cs0, m, A)
-
-    v0 = DM_VELOCITY
-
-    vesc = ESCAPE_VELOCITY
-
-    # earth's velocity through the galaxy
-    ve = EARTH_VELOCITY
-
-    # minimum wimp velocity needed to produce a recoil of energy e
-    vmin = np.sqrt(mn*e/2)*(mn+m)*SPEED_OF_LIGHT/(mn*m)
-
-    f = get_nuclear_form_factor(A, e)
-
-    x = vmin/v0
-    y = ve/v0
-    z = vesc/v0
-
-    # Equation 3.49 in Mark Pepin's thesis
-    k0 = (pi*v0**2)**(3/2)
-
-    # Equation 3.49 in Mark Pepin's thesis
-    k1 = k0*(erf(z)-(2/np.sqrt(pi))*z*np.exp(-z**2))
-
-    # From Mark Pepin's CDMS thesis page 59
-    if x <= z - y:
-        I = (k0/k1)*(DM_DENSITY/(2*ve))*(erf(x+y) - erf(x-y) - (4/np.sqrt(pi))*y*np.exp(-z**2))
-    elif x <= y + z:
-        I = (k0/k1)*(DM_DENSITY/(2*ve))*(erf(z) - erf(x-y) - (2/np.sqrt(pi))*(y+z-x)*np.exp(-z**2))
-    else:
-        return 0
-
-    return (N0/A)*(mn/(2*m*mr**2))*cs*f**2*I*SPEED_OF_LIGHT**2
-
-def get_event_rate_xenon(m, cs, A, threshold):
-    """
-    Returns the event rate of WIMP scattering in a dark matter detector using
-    an element with atomic mass A. Rate is in Hz kg^-1.
-    """
-    # mass of nucleus
-    mn = A*ATOMIC_MASS_UNIT
-
-    # reduced mass of the nucleus and the WIMP
-    mr = (m*mn)/(m + mn)
-
-    vesc = ESCAPE_VELOCITY
-
-    # earth's velocity through the galaxy
-    ve = EARTH_VELOCITY
-
-    # max recoil they looked for was 40 keV
-    emax = 2*mr**2*((ve+vesc)/SPEED_OF_LIGHT)**2/mn
-
-    rate, err = quad(get_differential_event_rate_xenon, threshold, emax, args=(m,cs,A), epsabs=0, epsrel=1e-2, limit=1000)
-
-    return rate
-
-if __name__ == '__main__':
-    import matplotlib.pyplot as plt
-
-    ls = np.logspace(-1,8,1000)
-
-    cs = 1e-50 # cm^2
-
-    # FIXME: should use water density for L < 1 m, silicon density for L ~ 1
-    # km, and iron for L >> 1 km
-    rate = np.array([get_event_rate_sno(1, cs*1e-4, l, water) for l in ls])
-
-    colors = plt.rcParams['axes.prop_cycle'].by_key()['color']
-
-    plt.figure()
-    plt.subplot(111)
-    plt.plot(ls, rate/np.max(rate),color=colors[0])
-    plt.xlabel("Mediator Decay Length (m)")
-    plt.ylabel("Event Rate (arbitrary units)")
-    plt.axvline(x=FIDUCIAL_RADIUS, color=colors[1], ls='--', label="SNO radius")
-    plt.axvline(x=h, ls='--', color=colors[2], label="Depth of SNO")
-    plt.axvline(x=R, ls='--', color=colors[3], label="Earth radius")
-    plt.gca().set_xscale('log')
-    plt.gca().set_xlim((ls[0], ls[-1]))
-    plt.title("Event Rate of Self-Destructing Dark Matter in SNO")
-    plt.legend()
-
-    # threshold is ~5 keVr
-    xenon_100t_threshold = 5e-6
-
-    # threshold is ~1 keVr
-    cdms_threshold = 1e-6
-
-    # threshold is ~100 eV
-    # FIXME: is this correct?
-    cresst_threshold = 1e-7
-
-    ms = np.logspace(-2,3,200)
-    cs0s = np.logspace(-50,-40,200)
-
-    mm, cs0cs0 = np.meshgrid(ms, cs0s)
-
-    rate1 = np.empty(mm.shape)
-    rate2 = np.empty(mm.shape)
-    rate3 = np.empty(mm.shape)
-    rate4 = np.empty(mm.shape)
-    rate5 = np.empty(mm.shape)
-    rate6 = np.empty(mm.shape)
-    for i in range(mm.shape[0]):
-        print("\r%i/%i" % (i+1,mm.shape[0]),end='')
-        for j in range(mm.shape[1]):
-            rate1[i,j] = get_event_rate_xenon(mm[i,j], cs0cs0[i,j]*1e-4, XENON_MASS, xenon_100t_threshold)
-            rate2[i,j] = get_event_rate_sno(mm[i,j], cs0cs0[i,j]*1e-4, 1.0, water)
-            rate3[i,j] = get_event_rate_xenon(mm[i,j], cs0cs0[i,j]*1e-4, GERMANIUM_MASS, cdms_threshold)
-            rate4[i,j] = get_event_rate_xenon(mm[i,j], cs0cs0[i,j]*1e-4, TUNGSTEN_MASS, cresst_threshold)
-            rate5[i,j] = get_event_rate_sno(mm[i,j], cs0cs0[i,j]*1e-4, 1e3, norite)
-            rate6[i,j] = get_event_rate_sno(mm[i,j], cs0cs0[i,j]*1e-4, 1e6, mantle)
-    print()
-
-    # Fiducial volume of the Xenon1T detector is 1042 +/- 12 kg
-    # from arxiv:1705.06655
-    xenon_100t_fiducial_volume = 1042 # kg
-
-    # Livetime of the Xenon1T results
-    # from arxiv:1705.06655
-    xenon_100t_livetime = 34.2 # days
-
-    plt.figure()
-    plt.subplot(111)
-    plt.gca().set_xscale('log')
-    plt.gca().set_yscale('log')
-    CS1 = plt.contour(mm,cs0cs0,rate1*3600*24*xenon_100t_fiducial_volume*xenon_100t_livetime,[10.0], colors=[colors[0]])
-    CS2 = plt.contour(mm,cs0cs0,rate2*3600*24*668.8,[10.0],colors=[colors[1]])
-    CS3 = plt.contour(mm,cs0cs0,rate3*3600*24*70.10,[10.0],colors=[colors[2]])
-    # FIXME: I used 2.39 kg day because that's what CRESST-3 reports in their paper
-    # but! I only use Tungsten here so do I need to multiply by the mass fraction of tungsten?
-    CS4 = plt.contour(mm,cs0cs0,rate4*3600*24*2.39,[10.0],colors=[colors[3]])
-    CS5 = plt.contour(mm,cs0cs0,rate5*3600*24*668.8,[10.0],colors=[colors[1]], linestyles=['dashed'])
-    CS6 = plt.contour(mm,cs0cs0,rate6*3600*24*668.8,[10.0],colors=[colors[1]], linestyles=['dotted'])
-    plt.clabel(CS1, inline=1, fmt="XENON1T", fontsize=10, use_clabeltext=True)
-    plt.clabel(CS2, inline=1, fmt=r"SNO ($\mathrm{L}_V$ = 1 m)", fontsize=10)
-    plt.clabel(CS3, inline=1, fmt="CDMSLite", fontsize=10)
-    plt.clabel(CS4, inline=1, fmt="CRESST-3", fontsize=10)
-    plt.clabel(CS5, inline=1, fmt=r"SNO ($\mathrm{L}_V$ = 1 km)", fontsize=10)
-    plt.clabel(CS6, inline=1, fmt=r"SNO ($\mathrm{L}_V$ = 1000 km)", fontsize=10)
-    plt.xlabel(r"$m_\chi$ (GeV)")
-    plt.ylabel(r"WIMP-nucleon scattering cross section ($\mathrm{cm}^2$)")
-    plt.title("Self-Destructing Dark Matter Limits")
-
-    x = np.linspace(0,300e-6,1000)
-    
-    # reproducing Figure 2.1 in Tom Caldwell's thesis
-    plt.figure()
-    plt.plot(x*1e6, list(map(lambda x: get_nuclear_form_factor(XENON_MASS, x)**2,x)), label="Xe")
-    plt.plot(x*1e6, list(map(lambda x: get_nuclear_form_factor(NEON_MASS, x)**2,x)), label="Ne")
-    plt.plot(x*1e6, list(map(lambda x: get_nuclear_form_factor(ARGON_MASS, x)**2,x)), label="Ar")
-    plt.plot(x*1e6, list(map(lambda x: get_nuclear_form_factor(GERMANIUM_MASS, x)**2,x)), label="Ge")
-    plt.xlabel("Recoil Energy (keV)")
-    plt.ylabel(r"$F^2(E)$")
-    plt.gca().set_yscale('log')
-    plt.legend()
-    plt.gca().set_ylim((1e-5,1))
-    plt.gca().set_xlim((0,300))
-    plt.title("Nuclear Form Factors for various elements")
-    plt.show()
diff --git a/doc/Makefile b/doc/Makefile
new file mode 100644
index 0000000..10cf3fe
--- /dev/null
+++ b/doc/Makefile
@@ -0,0 +1,9 @@
+all: sddm.pdf
+
+sddm.pdf: sddm.tex
+	pdflatex $^
+
+clean:
+	rm sddm.pdf
+
+.PHONY: clean sddm.pdf
diff --git a/sddm.tex b/doc/sddm.tex
index 1f3fa96..1f3fa96 100644
--- a/sddm.tex
+++ b/doc/sddm.tex