Optimum interval updated to consider only the energies simulated

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, 1e5)
+        en_interp = np.geomspace(elow, ehigh, int(1e5))
 
     event_inds = (eventenergies > elow) & (eventenergies < ehigh)
 

darklim/sensitivity/_sens_est.py

@@ -1,6 +1,7 @@
 import numpy as np
 import matplotlib.pyplot as plt
 from scipy import stats, special
+from scipy.interpolate import interp1d
 
 import mendeleev
 from darklim import constants
@@ -11,7 +12,10 @@
 import darklim.elf._elf as elf
 import darklim.detector._detector as detector
 import time
-np.random.seed(int(time.time()))
+
+E_LOW_GLOBAL_KEV = 1e-6
+E_HIGH_GLOBAL_KEV = 100.
+NPTS_GLOBAL = int(1e5)
 
 __all__ = [
     "calculate_substrate_mass",
@@ -157,7 +161,7 @@ class SensEst(object):
 
     """
 
-    def __init__(self, m_det, time_elapsed, eff=1, tm="Si", gain=1):#, signal_name='SI-NR'):
+    def __init__(self, m_det, time_elapsed, eff=1, tm="Si", gain=1, seed=None):#, signal_name='SI-NR'):
         """
         Initialization of the SensEst class.
 
@@ -185,6 +189,9 @@ def __init__(self, m_det, time_elapsed, eff=1, tm="Si", gain=1):#, signal_name='
         self._backgrounds = []
         self._background_labels = []
 
+        if seed is not None:
+            np.random.seed(seed)
+
         #self.signal_name = signal_name
 
         #self.signal = None
@@ -294,8 +301,8 @@ def reset_sim(self):
         self._backgrounds = []
 
 
-    def run_sim(self, threshold, e_high=1., e_low=1e-6, m_dms=np.geomspace(0.01, 2, num=5),
-                nexp=1, npts=1000, plot_bkgd=False, res=None, verbose=False, sigma0=1e-41,
+    def run_sim(self, threshold, e_high=E_HIGH_GLOBAL_KEV, e_low=E_LOW_GLOBAL_KEV, m_dms=np.geomspace(0.01, 2, num=5),
+                nexp=1, npts=NPTS_GLOBAL, plot_bkgd=False, res=None, verbose=False, sigma0=1e-41,
                 elf_model=None, elf_params=None, elf_target=None,
                 gaas_params=None, return_only_drde=False):
 
@@ -437,30 +444,41 @@ def run_sim(self, threshold, e_high=1., e_low=1e-6, m_dms=np.geomspace(0.01, 2,
 
         sigs = []
 
-        en_interp = np.geomspace(e_low, e_high, num=npts) # keV
-        rate_interp = []
+        en_interp_wide = np.geomspace(max(e_low, threshold), e_high, num=npts)
+        rate_interp_wide = [None for _ in range(len(m_dms))]
 
         for ii in range(len(m_dms)):
 
             t_start = time.time()
             try:
-                rate_temp = drdefunction[ii](en_interp) * self.exposure
+                rate_temp = drdefunction[ii](en_interp_wide) * self.exposure
             except ValueError:
-                rate_temp = np.array([drdefunction[ii](en) for en in en_interp]) * self.exposure
+                rate_temp = np.array([drdefunction[ii](en) for en in en_interp_wide]) * self.exposure
 
-            rate_interp.append(rate_temp)
+            rate_interp_wide[ii] = rate_temp
 
             print(f'Finished mass {ii}. Took {(time.time()-t_start)/60:.2f} minutes.')
 
         for jj in range(nexp):
 
             evts_sim = self._generate_background(
-                en_interp, plot_bkgd=(plot_bkgd and jj==0))
+                en_interp_wide, plot_bkgd=(plot_bkgd and jj==0))
             if jj == 0:
-                print(f'In the first pseudoexperiment, we simulated {len(evts_sim)} events')
+                print(f'In the first pseudoexperiment for mass {m_dms[0]} GeV, we simulated {len(evts_sim)} events')
+
+            # Combine original en_interp with event energies and sort them
+            combined_energies = np.unique(np.concatenate((en_interp_wide, evts_sim)))
+            min_event, max_event = min(evts_sim), max(evts_sim)
+            en_interp = combined_energies[(combined_energies >= min_event) & (combined_energies <= max_event)]
+
+            # Define interpolation function based on en_interp and rate_interp
+            rate_interp = [None for _ in range(len(m_dms))]
+            for ii in range(len(m_dms)):
+                interp_func = interp1d(en_interp_wide, rate_interp_wide[ii], kind='linear', bounds_error=True)
+                rate_interp[ii] = np.copy(interp_func(en_interp))
 
             sig_temp, _, _ = optimuminterval(
-                evts_sim[evts_sim >= threshold], # evt energies
+                evts_sim, # evt energies
                 en_interp, # efficiency curve energies
                 np.ones_like(en_interp), # efficiency curve values
                 m_dms, # mass list
@@ -471,7 +489,7 @@ def run_sim(self, threshold, e_high=1., e_low=1e-6, m_dms=np.geomspace(0.01, 2,
                 gauss_width=10, # if smearing, number of sigma to go out to
                 verbose=(verbose*(jj==0)), # print outs
 #                drdefunction=drdefunction, # lambda function for dRdE(E)
-#                hard_threshold=0., # hard threshold for energies
+                hard_threshold=threshold, # hard threshold for energies
                 sigma0=sigma0, # Starting guess for sigma
                 en_interp=en_interp,
                 rate_interp=rate_interp,

examples/sapphire_oi_scan.py

@@ -58,56 +58,21 @@ def plot_dm_rates(m_dms,dm_rates,raw_dm_rates,sigma0,savename=None):
 def process_mass(mass, args):
     # All the code that processes the mass value goes here, extracted from the original loop.
 
-    SE = darklim.sensitivity.SensEst(args.target_mass_kg, args.t_days, tm=args.target, eff=1., gain=1.)
+    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)
+    #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)
 
     per_device_threshold_keV = args.nsigma * args.baseline_res_eV * 1e-3
     threshold_keV = args.coincidence * per_device_threshold_keV
 
-    # First, figure out what the maximum energy from this dRdE is
-    e_high_keV = 100. # keV
-    e_low_keV = 1e-6 # keV
-    drdefunction = SE.run_sim(
-        threshold_keV,
-        e_high=e_high_keV,
-        e_low=e_low_keV,
-        m_dms=[mass],
-        sigma0=args.sigma0,
-        elf_model=args.elf_model,
-        elf_target=args.target,
-        elf_params=args.elf_params,
-        return_only_drde=True,
-#            gaas_params=None
-        )
-    drdefunction = drdefunction[0]
-
-    e_high_guesses = np.geomspace(e_low_keV, e_high_keV, 3000)
-    skip = False
-    try:
-        drdefunction_guesses = drdefunction(e_high_guesses)
-    except ValueError:
-        drdefunction_guesses = np.array([drdefunction(en) for en in e_high_guesses])
-    indices = np.where(drdefunction_guesses > 0)
-    if len(indices[0]) == 0:
-        e_high_keV = threshold_keV * 1.1
-    else:
-        j = int(indices[0][-1])
-        e_high_keV = e_high_guesses[j] * 1.1
-        if e_high_keV < threshold_keV:
-            e_high_keV = threshold_keV * 1.1
-
-    if skip:
-        sigma = np.inf
-    else:
-        _, sigma = SE.run_sim(
+    _, sigma = SE.run_sim(
             threshold_keV,
-            e_high=e_high_keV,
-            e_low=1e-6,
+            #e_high=e_high_keV,
+            #e_low=1e-6,
             m_dms=[mass],
             nexp=args.nexp,
-            npts=100000,
+            #npts=100000,
             plot_bkgd=False,
             res=None,
             verbose=True,
@@ -117,7 +82,7 @@ def process_mass(mass, args):
             elf_params=args.elf_params,
             return_only_drde=False,
 #            gaas_params=None
-        )
+    )
 
     print(f'Done mass = {mass}, sigma = {sigma}')