Up to date as of 09/18/24

darklim/elf/_elf.py

@@ -1,13 +1,15 @@
 from IPython.utils import io
 import numpy as np
 import sys
-sys.path.insert(0, '/global/cfs/cdirs/lz/users/vvelan/Test/DarkELF/')
+sys.path.insert(0, '/home/michael/DarkELF/')
 from darkelf import darkelf
 from darklim import constants
 
 __all__ = [
     "get_dRdE_lambda_Al2O3_electron",
     "get_dRdE_lambda_GaAs_electron",
+    "get_dRdE_lambda_Si_electron",
+    "get_dRdE_lambda_Si_phonon",
 ]
 
 def get_dRdE_lambda_Al2O3_electron(mX_eV=1e8, mediator='massless', sigmae=1e-31, kcut=0, method='grid', withscreening=True, suppress_darkelf_output=False, gain=1.):
@@ -208,3 +210,100 @@ def get_dRdE_lambda_GaAs_phonon(mX_eV=1e8, mediator='massless', sigman=1e-31, da
 
     return fun
 
+def get_dRdE_lambda_Si_electron(mX_eV=1e8, mediator='massless', sigmae=1e-31, kcut=0, method='grid', withscreening=True, suppress_darkelf_output=False, gain=1.):
+    """
+    Function to get an anonymous lambda function, which calculates dRdE
+    for DM-electron scattering in Si given only deposited energy.
+
+    Parameters
+    ----------
+    mX_eV : float
+        Dark matter mass in eV
+    mediator : str
+        Dark photon mediator mass. Must be "massive" (infinity) or
+        "massless" (zero).
+    sigmae : float
+        DM-electron scattering cross section in cm^2
+    kcut : float
+        Maximum k value in the integration, in eV. If kcut=0 (default), the
+        integration is cut off at the highest k-value of the grid at hand.
+    method : str
+        Must be "grid" or "Lindhard". Choice to use interpolated grid of
+        epsilon, or Lindhard analytic epsilon
+    withscreening : bool
+        Whether to include the 1/|epsilon|^2 factor in the scattering rate
+    suppress_darkelf_output : bool
+        Whether to suppress the (useful but long) output that DarkELF gives
+        when loading a material's properties.
+
+    Returns
+    -------
+    fun : lambda function
+        A function to calculate dRdE in DRU given E 
+
+    """
+
+    # Set up DarkELF GaAs object
+    if suppress_darkelf_output:
+        print('WARNING: You are suppressing DarkELF output')
+        with io.capture_output() as captured:
+            Si = darkelf(target='Si', filename="Si_mermin.dat")
+    else:
+        Si = darkelf(target='Si', filename="Si_mermin.dat")
+
+    # Create anonymous function to get rate with only deposited energy
+    # Note DarkELF expects recoil energies and WIMP masses in eV, and returns rates in counts/kg/yr/eV
+    # But DarkLim expects recoil energies in keV, WIMP masses in GeV, and rates in counts/kg/day/keV (DRU)
+    Si.update_params(mX=mX_eV, mediator=mediator)
+    fun = lambda keV : np.heaviside(keV * 1000 / gain - constants.bandgap_Si_eV, 1) * \
+            Si.dRdomega_electron(keV * 1000 / gain, method=method, sigmae=sigmae, kcut=kcut, withscreening=withscreening) * \
+            (1000 / 365.25) / gain
+
+    return fun
+
+
+
+
+def get_dRdE_lambda_Si_phonon(mX_eV=1e8, mediator='massless', sigman=1e-31, dark_photon=False, suppress_darkelf_output=False, gain=1.):
+    """
+    Function to get an anonymous lambda function, which calculates dRdE
+    for DM-nuclear scattering via phonons in Al2O3 given only deposited energy.
+
+    Parameters
+    ----------
+    mX_eV : float
+        Dark matter mass in eV
+    mediator : str
+        Dark photon mediator mass. Must be "massive" (infinity) or
+        "massless" (zero).
+    sigman : float
+        DM-nucleon scattering cross section in cm^2
+    dark_photon : bool
+        Whether to treat this as a dark photon
+    suppress_darkelf_output : bool
+        Whether to suppress the (useful but long) output that DarkELF gives
+        when loading a material's properties.
+
+    Returns
+    -------
+    fun : lambda function
+        A function to calculate dRdE in DRU given E 
+
+    """
+
+    # Set up DarkELF GaAs object
+    if suppress_darkelf_output:
+        print('WARNING: You are suppressing DarkELF output')
+        with io.capture_output() as captured:
+            Si = darkelf(target='Si', filename="Si_mermin.dat", phonon_filename='Si_epsphonon_o.dat')
+    else:
+        Si = darkelf(target='Si', filename="Si_mermin.dat", phonon_filename='Si_epsphonon_o.dat')
+
+    # Create anonymous function to get rate with only deposited energy
+    # Note DarkELF expects recoil energies and WIMP masses in eV, and returns rates in counts/kg/yr/eV
+    # But DarkLim expects recoil energies in keV, WIMP masses in GeV, and rates in counts/kg/day/keV (DRU)
+    Si.update_params(mX=mX_eV, mediator=mediator)
+    fun = lambda keV : Si._dR_domega_multiphonons_no_single(keV * 1000 / gain, sigman=sigman, dark_photon=dark_photon) * \
+            (1000 / 365.25) / gain
+
+    return fun

darklim/limit/_limit.py

@@ -12,6 +12,7 @@
 
 from darklim import constants, feldman_cousins
 from darklim.limit import _upper
+import darklim.elf._elf as elf
 import mendeleev
 
 import matplotlib.pyplot as plt
@@ -23,6 +24,7 @@
     "drde",
     "drde_max_q",
     "gauss_smear",
+    "drde_darkelf",
     "optimuminterval",
     "fc_limits",
     "get_fc_ul",
@@ -389,11 +391,40 @@ def gauss_smear(x, f, res, nres=1e5, gauss_width=10):
 
     return s(x)
 
+def drde_darkelf(m_dm,elf_model=None, elf_params=None, elf_target=None,sigma0=1e-41):
+
+    drdefunction = None
+
+    if elf_model == 'electron' and elf_target == 'Si':
+        elf_mediator = elf_params['mediator'] if 'mediator' in elf_params else 'massless'
+        elf_kcut = elf_params['kcut'] if 'kcut' in elf_params else 0
+        elf_method = elf_params['method'] if 'method' in elf_params else 'grid'
+        elf_screening = elf_params['withscreening'] if 'withscreening' in elf_params else True
+        elf_suppress = elf_params['suppress_darkelf_output'] if 'suppress_darkelf_output' in elf_params else False
+
+        drdefunction = \
+            [elf.get_dRdE_lambda_Si_electron(mX_eV=m*1e9, sigmae=sigma0, mediator=elf_mediator,
+                                                kcut=elf_kcut, method=elf_method, withscreening=elf_screening,
+                                                suppress_darkelf_output=elf_suppress, gain=1)
+            for m in m_dm]
+
+    elif elf_model == 'phonon' and elf_target == 'Si':
+        elf_mediator = elf_params['mediator'] if 'mediator' in elf_params else 'massless'
+        elf_suppress = elf_params['suppress_darkelf_output'] if 'suppress_darkelf_output' in elf_params else False
+        elf_darkphoton = elf_params['dark_photon'] if 'dark_photon' in elf_params else False
+
+        drdefunction = \
+            [elf.get_dRdE_lambda_Si_phonon(mX_eV=m*1e9, sigman=sigma0, mediator=elf_mediator,
+                                                dark_photon=elf_darkphoton,
+                                                suppress_darkelf_output=elf_suppress, gain=1)
+            for m in m_dm]
+
+    return drdefunction
 
 
 def optimuminterval(eventenergies, effenergies, effs, masslist, exposure,
                     tm="Si", cl=0.9, res=None, gauss_width=10, verbose=False,
-                    drdefunction=None, hard_threshold=0.0, sigma0=1e-41):
+                    drdefunction=None, hard_threshold=0.0, sigma0=1e-41,iself=False):
     """
     Function for running Steve Yellin's Optimum Interval code on an inputted spectrum and efficiency curve.
 
@@ -438,12 +469,17 @@ def optimuminterval(eventenergies, effenergies, effs, masslist, exposure,
         `masslist` and the cross section sigma=10^-41 cm^2. The experiment efficiency must be taken
         into account. The energy unit is keV, the rate unit is 1/keV/kg/day.
         By default (or if None is provided) the standard Lewin&Smith signal model is used with gaussian
-        smearing of width `res`, truncated at `gauss_width` standard deviations.
+        smearing of width `res`, truncated at `gauss_width` standard deviations. 
+        If iself is True, then this function does not need to be a callable function, and can instead be
+        a nested list of rates for given masses. Assumes guassian smearing has already been done. 
     hard_threshold : float, optional
         The energy value (keV) below which the efficiency is zero.
         This argument is not required in a case of smooth efficiency curve, however it must be provided 
         in a case of step-function-like efficiency.
-
+    iself : bool, optional
+        If True, then the algoroithm will allow you to insert a non-callable function for drdefunction. 
+        This allows for calculation of a drde spectrum produced by DarkElf to be created as a callable 
+        function and then executed, producing a list of rates. 
     Returns
     -------
     sigma : ndarray
@@ -472,7 +508,7 @@ def optimuminterval(eventenergies, effenergies, effs, masslist, exposure,
     elow = max(hard_threshold, min(effenergies))
     ehigh = max(effenergies)
 
-    en_interp = np.logspace(np.log10(elow), np.log10(ehigh), int(1e5))
+    en_interp = np.logspace(np.log10(elow), np.log10(ehigh), int(1e4))
 
     #sigma0 = 1e-41
 
@@ -485,30 +521,43 @@ def optimuminterval(eventenergies, effenergies, effs, masslist, exposure,
     for ii, mass in enumerate(masslist):
         if verbose:
             print(f"On mass {ii+1} of {len(masslist)}.")
-
-        if drdefunction is None:
+
+        if iself is False:    
+
+            if drdefunction is None:
+                exp = effs * exposure
+
+                curr_exp = interpolate.interp1d(
+                    effenergies, exp, kind="linear", bounds_error=False, fill_value=(0, exp[-1]),
+                )
+
+                init_rate = drde(
+                    en_interp, mass, sigma0, tm=tm,
+                )
+                if res is not None:
+                    init_rate = gauss_smear(en_interp, init_rate, res, gauss_width=gauss_width)
+                rate = init_rate * curr_exp(en_interp)
+            else:
+                try:
+                    rate = drdefunction[ii](en_interp) * exposure
+                except ValueError:
+                    rate = np.array([drdefunction[ii](en) for en in en_interp]) * exposure
+
+        else:
+            print("iself")
             exp = effs * exposure
-
             curr_exp = interpolate.interp1d(
                 effenergies, exp, kind="linear", bounds_error=False, fill_value=(0, exp[-1]),
             )
-
-            init_rate = drde(
-                en_interp, mass, sigma0, tm=tm,
-            )
-            if res is not None:
-                init_rate = gauss_smear(en_interp, init_rate, res, gauss_width=gauss_width)
-            rate = init_rate * curr_exp(en_interp)
-        else:
-            try:
-                rate = drdefunction[ii](en_interp) * exposure
-            except ValueError:
-                rate = np.array([drdefunction[ii](en) for en in en_interp]) * exposure
-
+            print(drdefunction[ii])
+            rate = np.array(drdefunction[ii] * curr_exp(en_interp))
+            print(rate)
 
         integ_rate = integrate.cumtrapz(rate, x=en_interp, initial=0)
+        print(integ_rate)
 
         tot_rate = integ_rate[-1]
+        print(tot_rate)
 
         x_val_fcn = interpolate.interp1d(
             en_interp,