Calculate limit on real data with DarkELF or NR

darklim/constants/_constants.py

@@ -9,6 +9,7 @@
 
 bandgap_GaAs_eV = 1.42
 bandgap_Al2O3_eV = 8.8
+bandgap_Si_eV = 1.12
 Al2O3_density = 3.98
 GaAs_density = 5.32
 GaAs_light_fraction = 0.60

darklim/elf/_elf.py

@@ -8,6 +8,7 @@
 __all__ = [
     "get_dRdE_lambda_Al2O3_electron",
     "get_dRdE_lambda_GaAs_electron",
+    "get_dRdE_lambda_Si_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.):
@@ -59,6 +60,8 @@ def get_dRdE_lambda_Al2O3_electron(mX_eV=1e8, mediator='massless', sigmae=1e-31,
             sapphire.dRdomega_electron(keV * 1000 / gain, method=method, sigmae=sigmae, kcut=kcut, withscreening=withscreening) * \
             (1000 / 365.25) / gain
 
+    print(f'The mass is {mX_eV/1e6} MeV/c2')
+
     return fun
 
 
@@ -208,3 +211,57 @@ 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
+
+
+

darklim/limit/_limit.py

@@ -473,7 +473,7 @@ def optimuminterval(eventenergies, effenergies, effs, masslist, exposure,
     ehigh = max(effenergies)
 
     if en_interp is None:
-        en_interp = np.geomspace(elow, ehigh, int(1e5))
+        en_interp = np.geomspace(elow, ehigh, int(1e4))
 
     event_inds = (eventenergies > elow) & (eventenergies < ehigh)
 

examples/calculate_oi_limit.py

@@ -0,0 +1,106 @@
+import os
+import sys
+import numpy as np
+import matplotlib.pyplot as plt
+import scipy.stats as stats
+
+import darklim
+from darklim import constants
+import darklim.elf as elf
+
+import scanparser
+
+from multihist import Hist1d
+import time
+import datetime
+
+import multiprocessing as mp
+
+
+##################################################################
+
+
+def process_mass(mass, args, shared_energies_eV):
+    # All the code that processes the mass value goes here, extracted from the original loop.
+
+    eventenergies_keV = sorted(np.random.choice(shared_energies_eV, size=30, replace=False) * 1e-3)
+    print('The chosen energies are:', eventenergies_keV)
+
+    effenergies_keV = np.linspace(min(eventenergies_keV), max(eventenergies_keV), 10)
+    effs = np.full(len(effenergies_keV), 0.7955)
+
+    elf_mediator = 'massive'
+    elf_kcut = 0
+    elf_method = 'grid'
+    elf_screening = True
+    elf_suppress = False
+
+    drdefunction = \
+        [elf.get_dRdE_lambda_Si_electron(mX_eV=m*1e9, sigmae=args.sigma0, mediator=elf_mediator,
+                                            kcut=elf_kcut, method=elf_method, withscreening=elf_screening,
+                                            suppress_darkelf_output=elf_suppress, gain=1.)
+         for m in mass]
+
+    E_deposited_keV_arr = np.geomspace(min(eventenergies_keV), max(eventenergies_keV), 1000)
+    for j, m in enumerate(mass):
+        #try:
+        #    dRdE_DRU_arr = drdefunction[j](E_deposited_keV_arr)
+        #except ValueError:
+        dRdE_DRU_arr = np.array([drdefunction[j](en) for en in E_deposited_keV_arr])
+
+        check = sum(dRdE_DRU_arr > 0)
+        if check == 0:
+            print(f'No observed rate for mass {j}, {m} GeV')
+
+        drdefunction[j] = lambda E: np.interp(E, E_deposited_keV_arr, dRdE_DRU_arr, left=0., right=0.)
+
+    sigma, _, _ = darklim.limit.optimuminterval(eventenergies_keV, effenergies_keV, effs, mass, args.exposure_kgd,
+                    tm=args.target, cl=0.9, res=args.baseline_res_eV*1e-3, gauss_width=3, verbose=True,
+                    drdefunction=drdefunction, hard_threshold=args.baseline_res_eV*5*1e-3, sigma0=1e-41,
+                    en_interp=None, rate_interp=None)
+
+    print(f'Done mass = {mass}, sigma = {sigma}')
+
+    return mass, sigma
+
+
+
+def sapphire_scan():
+
+    # Read command-line arguments
+    args = scanparser.get_scan_parameters()
+
+    # Write input parameters to a text file
+    scanparser.write_info(args)
+
+    # Main parallel execution block
+    shared_energies_eV = np.load('spectra/BigFins_shared_0719.npy')
+    #process_mass(args.masses_GeV, args, shared_energies_eV)
+
+    if True:
+        with mp.Pool(processes=min(args.max_cpus, mp.cpu_count())) as pool:
+            results = pool.starmap(process_mass, [([mass], args, shared_energies_eV) for mass in args.masses_GeV])
+
+        # save results to txt file
+        sigma = np.zeros_like(args.masses_GeV)
+        for i, result in enumerate(results):
+            sigma[i] = result[1][0]
+
+        outname = args.results_dir + 'limit.txt'
+        tot = np.column_stack( (args.masses_GeV, sigma) )
+        np.savetxt(outname, tot, fmt=['%.5e','%0.5e'], delimiter=' ')
+
+    return
+
+
+# ------------------------------------------------------
+# ------------------------------------------------------
+
+
+if __name__ == "__main__":
+
+    t_start = time.time()
+    sapphire_scan()
+    t_end = time.time()
+    print(f'Full scan took {(t_end - t_start)/60:.2f} minutes.')
+

examples/defaults.py

@@ -11,7 +11,7 @@ def __init__(self):
         self.nexp = 100
         self.t_days = 5 / 60 / 24.
 
-        self.target = 'Al2O3'
+        self.target = 'Si'
         self.volume_cm3 = 1.
 
         self.n_sensors = 1
@@ -21,7 +21,7 @@ def __init__(self):
 
         self.baseline_res_eV = 0.1
 
-        self.masses_GeV = [3e-4, 1e3, 50]
+        self.masses_GeV = [1e-4, 1e2, 35]
         self.sigma0 = 1e-35
 
         self.elf_params_NR = {}

examples/sapphire_oi_scan.py

@@ -61,7 +61,7 @@ def process_mass(mass, args):
     SE = darklim.sensitivity.SensEst(args.target_mass_kg, args.t_days, tm=args.target, eff=1., gain=1., seed=(int(time.time() + mass*1e6)))
     SE.reset_sim()
     #SE.add_flat_bkgd(1) # flat background of 1 DRU
-    SE.add_nfold_lee_bkgd(m=args.n_sensors, n=args.coincidence, w=args.window_s, e0=0.41e-3, R=28.5)
+    SE.add_nfold_lee_bkgd(m=args.n_sensors, n=args.coincidence, w=args.window_s)#, e0=0.41e-3, R=28.5)
 
     per_device_threshold_keV = args.nsigma * args.baseline_res_eV * 1e-3
     threshold_keV = args.coincidence * per_device_threshold_keV

examples/scanparser.py

@@ -88,6 +88,9 @@ def get_scan_parameters():
     parser.add_argument('--results_dir', type=str, default=df.results_dir,
                         help='Output directory')
 
+    parser.add_argument('--fake', type=str, default=0,
+                        help='Ignore this')
+
     parser.add_argument('--max_cpus', type=int, default=df.max_cpus,
                         help='Maximum number of CPU cores to use')