merging with vetri

darklim/constants/_constants.py

@@ -13,4 +13,5 @@
 Al2O3_density = 3.98
 GaAs_density = 5.32
 GaAs_light_fraction = 0.60
-GaAs_average_phonon_energy_eV = 0.022
+GaAs_average_phonon_energy_eV = 0.022
+GaAs_average_phonon_energy_eV = 0.022

darklim/detector/__init__.py

@@ -0,0 +1 @@
+from ._detector import *

darklim/detector/_detector.py

@@ -0,0 +1,103 @@
+import numpy as np
+import math
+from darklim import constants
+import time
+from scipy import integrate, interpolate
+from darklim import elf
+import darklim.sensitivity._sens_est as sens_est
+import matplotlib.pyplot as plt
+import matplotlib as mpl
+
+def get_deposited_energy_gaas(E_recoil_eV, pce, lce_per_channel, res, n_coincidence_light, threshold_eV, coincidence_window_us, phonon_tau_us, n_samples=1):
+
+    # Get the light signal in each channel
+    E_light_eV = constants.GaAs_light_fraction * E_recoil_eV
+
+    n_photons_generated_average = np.floor(E_light_eV / constants.bandgap_GaAs_eV)
+    if n_photons_generated_average == 0:
+        if n_samples == 1:
+            return 0.
+        else:
+            return np.full(n_samples, 0.)
+
+    n_photons_detected_ch1 = np.random.binomial(n_photons_generated_average, lce_per_channel, n_samples)
+    n_photons_detected_ch2 = np.random.binomial(n_photons_generated_average, lce_per_channel, n_samples)
+
+    E_ch1_eV = n_photons_detected_ch1 * constants.bandgap_GaAs_eV * np.random.normal(1, res, n_samples)
+    E_ch2_eV = n_photons_detected_ch2 * constants.bandgap_GaAs_eV * np.random.normal(1, res, n_samples)
+
+    # Get the heat signal observed within coincidence window
+    E_heat_eV = (1 - constants.GaAs_light_fraction) * E_recoil_eV
+
+    n_phonons_generated_average = int(E_heat_eV / constants.GaAs_average_phonon_energy_eV)
+    phonon_arrival_times_us = np.random.exponential(phonon_tau_us, (n_phonons_generated_average, n_samples))
+    n_phonons_detected = np.sum(phonon_arrival_times_us < coincidence_window_us, axis=0)
+
+    E_ch0_eV = n_phonons_detected * pce * constants.GaAs_average_phonon_energy_eV * np.random.normal(1, res, n_samples)
+
+    if n_coincidence_light == 1:
+        E_det_eV = (E_ch0_eV + E_ch1_eV + E_ch2_eV) * (E_ch0_eV > threshold_eV) * ((E_ch1_eV > threshold_eV) + (E_ch2_eV > threshold_eV))
+    elif n_coincidence_light == 2:
+        E_det_eV = (E_ch0_eV + E_ch1_eV + E_ch2_eV) * (E_ch0_eV > threshold_eV) * (E_ch1_eV > threshold_eV) * (E_ch2_eV > threshold_eV)
+
+    if n_samples == 1:
+        return E_det_eV[0]
+    else:
+        return E_det_eV
+
+
+
+def convert_dRdE_dep_to_obs_gaas(E_dep_keV, dRdE_dep_DRU, pce=0.40, lce_per_channel=0.10, res=0.10, n_coincidence_light=1,
+    calorimeter_threshold_eV=0.37, coincidence_window_us=100., phonon_tau_us=100., E_min_keV=None, E_max_keV=None, n_samples=int(1e6)):
+
+    # Reduce data to the appropriate energy range
+    if E_min_keV is None or E_min_keV < E_dep_keV[0]:
+        E_min_keV = E_dep_keV[0]
+    if E_max_keV is None or E_max_keV > E_dep_keV[-1]:
+        E_max_keV = E_dep_keV[-1]
+
+    E_pdf = E_dep_keV[(E_dep_keV >= E_min_keV) * (E_dep_keV <= E_max_keV)]
+    dRdE_pdf = dRdE_dep_DRU[(E_dep_keV >= E_min_keV) * (E_dep_keV <= E_max_keV)]
+
+    # Draw samples from the distribution
+    cdf = integrate.cumtrapz(dRdE_pdf, x=E_pdf, initial=0.0)
+    cdf /= cdf[-1]
+
+    inv_cdf = interpolate.interp1d(cdf, E_pdf)
+
+    samples = np.random.rand(n_samples)
+
+    energies_sim_keV = inv_cdf(samples)
+    energies_obs_keV = np.zeros_like(energies_sim_keV)
+    energies_obs_keV = np.copy(energies_sim_keV)
+
+    for i, E in enumerate(energies_sim_keV):
+        energies_obs_keV[i] = get_deposited_energy_gaas(E * 1000, pce, lce_per_channel, res, n_coincidence_light, calorimeter_threshold_eV,
+            coincidence_window_us, phonon_tau_us) / 1000
+
+    # Perhaps no energy is ever observed
+    if sum(energies_obs_keV > 0) == 0:
+        return E_pdf, np.zeros_like(dRdE_pdf), np.array([])
+
+    # Convert to E vs dRdE that we can later interpolate from
+    # Normalize to the number of events that are detected
+    bins = np.geomspace(min(energies_obs_keV[energies_obs_keV > 0]) * 0.95, max(energies_obs_keV) * 1.05, 10000)
+    counts, bin_edges = np.histogram(energies_obs_keV, bins)
+    counts = counts * 1.0 / np.diff(bin_edges)
+    bin_centers = 0.5 * (bin_edges[1:] + bin_edges[:-1])
+
+    integral_original = sum(0.5 * (dRdE_pdf[1:] + dRdE_pdf[:-1]) * np.diff(E_pdf))
+    fraction_surviving = sum(energies_obs_keV > 0) / len(energies_obs_keV)
+    integral_desired = integral_original * fraction_surviving
+    integral_observed = sum(counts * np.diff(bin_edges))
+
+    E_obs_keV = np.copy(bin_centers)
+    dRdE_obs_DRU = counts * integral_desired / integral_observed
+
+    E_obs_keV = E_obs_keV[dRdE_obs_DRU > 0]
+    dRdE_obs_DRU = dRdE_obs_DRU[dRdE_obs_DRU > 0]
+
+    # Return arrays and the list of energies
+    return E_obs_keV, dRdE_obs_DRU, energies_obs_keV
+
+

darklim/elf/__init__.py

@@ -0,0 +1 @@
+from ._elf import *

darklim/elf/_elf.py

@@ -0,0 +1,210 @@
+from IPython.utils import io
+import numpy as np
+import sys
+sys.path.insert(0, '/global/cfs/cdirs/lz/users/vvelan/Test/DarkELF/')
+from darkelf import darkelf
+from darklim import constants
+
+__all__ = [
+    "get_dRdE_lambda_Al2O3_electron",
+    "get_dRdE_lambda_GaAs_electron",
+]
+
+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.):
+    """
+    Function to get an anonymous lambda function, which calculates dRdE
+    for DM-electron scattering 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).
+    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:
+            sapphire = darkelf(target='Al2O3', filename="Al2O3_mermin.dat")
+    else:
+        sapphire = darkelf(target='Al2O3', filename="Al2O3_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)
+    sapphire.update_params(mX=mX_eV, mediator=mediator)
+    fun = lambda keV : np.heaviside(keV * 1000 / gain - constants.bandgap_Al2O3_eV, 1) * \
+            sapphire.dRdomega_electron(keV * 1000 / gain, method=method, sigmae=sigmae, kcut=kcut, withscreening=withscreening) * \
+            (1000 / 365.25) / gain
+
+    return fun
+
+
+
+
+def get_dRdE_lambda_Al2O3_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:
+            sapphire = darkelf(target='Al2O3', filename="Al2O3_mermin.dat", phonon_filename='Al2O3_epsphonon_o.dat')
+    else:
+        sapphire = darkelf(target='Al2O3', filename="Al2O3_mermin.dat", phonon_filename='Al2O3_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)
+    sapphire.update_params(mX=mX_eV, mediator=mediator)
+    fun = lambda keV : sapphire._dR_domega_multiphonons_no_single(keV * 1000 / gain, sigman=sigman, dark_photon=dark_photon) * \
+            (1000 / 365.25) / gain
+
+    return fun
+
+
+
+def get_dRdE_lambda_GaAs_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 GaAs 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:
+            gaas = darkelf(target='GaAs', filename="GaAs_mermin.dat")
+    else:
+        gaas = darkelf(target='GaAs', filename="GaAs_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)
+    gaas.update_params(mX=mX_eV, mediator=mediator)
+    fun = lambda keV : np.heaviside(keV * 1000 / gain - constants.bandgap_GaAs_eV, 1) * \
+            gaas.dRdomega_electron(keV * 1000 / gain, method=method, sigmae=sigmae, kcut=kcut, withscreening=withscreening) * \
+            (1000 / 365.25) / gain
+
+    return fun
+
+
+
+
+def get_dRdE_lambda_GaAs_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 GaAs 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:
+            gaas = darkelf(target='GaAs', filename="GaAs_mermin.dat", phonon_filename='GaAs_epsphonon_data10K.dat')
+    else:
+        gaas = darkelf(target='GaAs', filename="GaAs_mermin.dat", phonon_filename='GaAs_epsphonon_data10K.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)
+    gaas.update_params(mX=mX_eV, mediator=mediator)
+    fun = lambda keV : gaas._dR_domega_multiphonons_no_single(keV * 1000 / gain, sigman=sigman, dark_photon=dark_photon) * \
+            (1000 / 365.25) / gain
+
+    return fun
+

darklim/limit/_limit.py

@@ -16,6 +16,7 @@
 
 import matplotlib.pyplot as plt
 
+
 __all__ = [
     "upper",
     "helmfactor",
@@ -499,7 +500,11 @@ def optimuminterval(eventenergies, effenergies, effs, masslist, exposure,
                 init_rate = gauss_smear(en_interp, init_rate, res, gauss_width=gauss_width)
             rate = init_rate * curr_exp(en_interp)
         else:
-            rate = drdefunction[ii](en_interp) * exposure
+            try:
+                rate = drdefunction[ii](en_interp) * exposure
+            except ValueError:
+                rate = np.array([drdefunction[ii](en) for en in en_interp]) * exposure
+
 
         integ_rate = integrate.cumtrapz(rate, x=en_interp, initial=0)
 
@@ -610,7 +615,8 @@ def fc_limits(known_bkg_func, eventenergies, effenergies, effs, masslist, exposu
             ax.set_xscale('log')
             ax.set_yscale('log')
             ax.legend()
-            outdir = '/global/cfs/cdirs/lz/users/haselsco/TESSERACT_Limits/DarkLim/examples/'
+            outdir = '/global/cfs/cdirs/lz/users/vvelan/Test/DarkLim/examples/'
+
             plt.savefig(outdir+'testplot_{:0.3f}GeV.png'.format(mass),dpi=300, facecolor='white',bbox_inches='tight')
 
         sigma[ii] = (sigma0 / tot_rate) * ul
@@ -640,7 +646,8 @@ def get_fc_ul(known_bkg_func, eventenergies, threshold, ehigh, exposure, verbose
 
 def get_signal_rate(effenergies, effs, masslist, exposure,
                     tm="Si", 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,
+                    savedir=''):
     """
     return the signal rate for each of the masses in masslist for the 
     given reference cross section, exposure, and efficiency curve.
@@ -665,7 +672,10 @@ def get_signal_rate(effenergies, effs, masslist, exposure,
         if drdefunction is None:
             init_rate = drde(en_interp, mass, sigma0, tm=tm)
         else:
-            init_rate = drdefunction[ii](en_interp,mass) # note here.. mass also an input
+            try:
+                init_rate = drdefunction[ii](en_interp)
+            except ValueError:
+                init_rate = np.array([drdefunction[ii](en) for en in en_interp])
 
         if res is not None:
             init_rate = gauss_smear(en_interp, init_rate, res, gauss_width=gauss_width)
@@ -691,7 +701,7 @@ def get_signal_rate(effenergies, effs, masslist, exposure,
             ax.set_yscale('log')
             ax.legend(loc='lower left',frameon=False)
             ax.set_title('m={:0.3f}GeV,\n rate over threshold={:0.3e} evts'.format(mass,signal_rates[ii]))
-            outdir = '/global/cfs/cdirs/lz/users/haselsco/TESSERACT_Limits/DarkLim/examples/'
-            plt.savefig(outdir+'testplot_{:0.3f}GeV.png'.format(mass),facecolor='white',bbox_inches='tight')
-                    
-    return signal_rates, raw_signal_rates
+            outdir = '/global/cfs/cdirs/lz/users/vvelan/Test/DarkLim/examples/'
+            plt.savefig(outdir+savedir+'/testplot_{:0.3f}GeV.png'.format(mass),facecolor='white',bbox_inches='tight')
+
+    return signal_rates, raw_signal_rates