diff --git a/invisible_cities/icaros/correction_functions.py b/invisible_cities/icaros/correction_functions.py
new file mode 100644
index 000000000..0e234376c
--- /dev/null
+++ b/invisible_cities/icaros/correction_functions.py
@@ -0,0 +1,51 @@
+import pandas as pd
+import numpy  as np
+from scipy.interpolate import griddata
+
+
+
+def normalization(krmap, method, x_low, x_high, y_low, y_high):
+
+    mu_values = krmap.mu
+
+    if method == 'max':
+        E_reference_max = mu_values.dropna().max()
+        return E_reference_max
+
+    if method == 'mean chamber':
+        E_reference_chamber = mu_values.dropna().mean()
+        return E_reference_chamber
+
+    if method == 'mean anode':
+        mu_values_anode = krmap[krmap.k == 0].mu
+        E_reference_anode = mu_values_anode.dropna().mean()
+        return E_reference_anode
+
+    mask_region = (krmap['x'] <= x_high) & (krmap['x'] >= x_low) & (krmap['y'] <= y_high) & (krmap['y'] >= y_low)
+
+    if method == 'mean region anode':
+        region = krmap[krmap.k == 0][mask_region]
+        E_reference_slice_anode = region.mu.dropna().mean()
+        return E_reference_slice_anode
+
+    if method == 'mean region chamber':
+        region = krmap[mask_region]
+        E_reference_region = region.mu.dropna().mean()
+        return E_reference_region
+
+
+def apply_3Dmap(krmap, method,x_low, x_high, y_low, y_high, dt, x, y, E, keV = False):
+
+    map_points = krmap['dt x y'.split()].values
+    norm = normalization(krmap, method, x_low, x_high, y_low, y_high)
+
+    data_points = np.stack([dt, x, y], axis = 1)
+    E_interpolated_data = griddata(map_points, krmap.mu.values, data_points, method = 'nearest')
+
+    correction_factor = norm/E_interpolated_data
+    Ec = E * correction_factor
+
+    if keV:
+        Ec = Ec * (41.55 / norm)
+
+    return Ec
diff --git a/invisible_cities/icaros/correction_functions_test.py b/invisible_cities/icaros/correction_functions_test.py
new file mode 100644
index 000000000..c6b72896c
--- /dev/null
+++ b/invisible_cities/icaros/correction_functions_test.py
@@ -0,0 +1,43 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import pandas as pd
+from correction_functions import normalization
+
+
+def test_method_norm():
+
+    k_vals = np.arange(10)
+    i_vals = np.arange(12)
+    j_vals = np.arange(12)
+
+    k, i, j = np.meshgrid(k_vals, i_vals, j_vals, indexing='ij')
+    k = k.ravel()
+    i = i.ravel()
+    j = j.ravel()
+
+    x = i * 25
+    y = j * 25
+    dt = k * 45
+
+    mu = (k*i + 800)*j
+
+    map_test = pd.DataFrame({
+        'k': k,
+        'i': i,
+        'j': j,
+        'x': x,
+        'y': y,
+        'dt': dt,
+        'mu': mu
+    })
+
+    region_mask = (
+        (map_test.x <= 100) & (map_test.x >= -100) &
+        (map_test.y <= 100) & (map_test.y >= -100)
+    )
+
+    assert normalization(map_test, 'max', None, None, None, None) == 9889
+    assert normalization(map_test, 'mean chamber', None, None, None, None) == 4536.125
+    assert normalization(map_test, 'mean anode', None, None, None, None) == 4400
+    assert normalization(map_test, 'mean region anode', -100, 100, -100, 100) == map_test.loc[(map_test.k == 0) & region_mask,'mu'].mean()
+    assert normalization(map_test, 'mean region chamber', -100, 100, -100, 100) == map_test.loc[region_mask,'mu'].mean()
diff --git a/invisible_cities/icaros/krmap_functions.py b/invisible_cities/icaros/krmap_functions.py
new file mode 100644
index 000000000..0c614a1aa
--- /dev/null
+++ b/invisible_cities/icaros/krmap_functions.py
@@ -0,0 +1,150 @@
+import pandas as pd
+import numpy  as np
+
+from invisible_cities.core.fit_functions import expo, gauss, fit
+from invisible_cities.io.dst_io import load_dst, load_dsts
+from invisible_cities.core.core_functions import in_range, shift_to_bin_centers
+
+import tables            as tb
+from   scipy             import stats
+from   scipy             import optimize
+from   scipy.optimize import curve_fit
+
+import itertools
+
+"""
+
+This script contains functions to compute NaN_maps, 3D_maps from kdsts and to merge them to create an 'uniform' map.
+
+"""
+
+def create_NaN_map(xy_range, dt_range, xy_nbins, dt_nbins):
+    """
+    - Uses every possible indices combination
+    - Creates a DataFrame filled with NaNs for each (i, j, k) combination
+    """
+
+    xy_bins = np.linspace(xy_range[0], xy_range[1], xy_nbins + 1)
+    dt_bins = np.linspace(dt_range[0], dt_range[1], dt_nbins + 1)
+
+    #shift to bin centers invisible_cities.core.core_functions
+
+    i_range = np.arange(0, len(xy_bins)-1)
+    j_range = np.arange(0, len(xy_bins)-1)
+    k_range = np.arange(0, len(dt_bins)-1)
+
+    combinations = pd.DataFrame(itertools.product(k_range, i_range, j_range),
+                                columns=['k', 'i', 'j'])
+    Nan_map = pd.DataFrame(
+        np.nan,
+        index=range(len(combinations)),
+        columns=['nevents', 'mu', 'sigma', 'mu_error', 'sigma_error']
+    )
+    Nan_map = pd.concat([combinations, Nan_map], axis=1)
+
+    return Nan_map
+
+
+def med_fun(df):
+    sigma = (0.04/2.35)*df.S2e.median()
+
+    return pd.DataFrame({
+         'nevents': len(df),
+         'median': df.S2e.median(),
+         'sigma': sigma,
+         'median_error': sigma/np.sqrt(len(df)),
+         'sigma_error': np.nan}, index = [0]) #error de std?
+
+
+
+def fit_function(df, bins, size =10):
+    if len(df)<size:
+        results = med_fun(df)
+
+    counts, bin_edges = np.histogram(df.S2e, bins)
+    bin_centers = shift_to_bin_centers(bin_edges)
+
+    try:
+        f = fit(gauss, bin_centers, counts, seed = (len(df), 8000, 400), maxfev = 2500)
+    except:
+
+        results = pd.DataFrame({'nevents': len(df), 'mu': np.nan, 'sigma': np.nan, 'mu_error': np.nan, 'sigma_error': np.nan}, index = [0])
+
+        return results
+    results = pd.DataFrame({'nevents': len(df), 'mu': f.values[1], 'sigma': f.values[2], 'mu_error': f.errors[1], 'sigma_error': f.errors[2]}, index = [0])
+
+    return results
+
+
+def map_3D_fits(df, xy_range, dt_range, xy_nbins, dt_nbins):
+
+    df = df.loc[:,['X', 'Y', 'DT', 'S2e']]
+
+    mask = in_range(df.X, *xy_range) & in_range(df.Y, *xy_range) & in_range(df.DT, *dt_range)
+
+    df = df[mask]
+
+    x = df.X
+    y = df.Y
+    dt = df.DT
+
+    xy_bins = np.linspace(*xy_range, xy_nbins + 1)
+    dt_bins = np.linspace(*dt_range, dt_nbins + 1)
+
+    map_df = df.assign( i = np.digitize(x.values, bins = xy_bins) - 1,
+                        j = np.digitize(y.values, bins = xy_bins) - 1,
+                        k = np.digitize(dt.values, bins = dt_bins) - 1)
+
+
+    result = map_df.groupby(['k', 'i', 'j']).apply(fit_function, bins = np.linspace(3000, 9000, 101))
+
+    result = result.reset_index()
+
+    if 'level_3' in result.columns:
+        result = result.drop(columns = 'level_3')
+
+    return result
+
+
+def merge_maps(NaN_map, map_3D):
+
+    full_map = pd.concat([map_3D, NaN_map], ignore_index = True)
+    full_map = full_map.groupby(['k', 'i','j']).first().reset_index()
+
+    return full_map
+
+
+def include_coordinates(map, xy_range, dt_range, xy_nbins, dt_nbins):
+
+    xy_bins = np.linspace(*xy_range, xy_nbins + 1)
+    dt_bins = np.linspace(*dt_range, dt_nbins + 1)
+
+    x = shift_to_bin_centers(xy_bins)
+    y = shift_to_bin_centers(xy_bins)
+    dt = shift_to_bin_centers(dt_bins)
+
+    combinaciones_bins = pd.DataFrame(
+    itertools.product(dt, x, y),
+    columns=['dt', 'x', 'y'])
+
+
+    map = pd.concat([combinaciones_bins.reset_index(drop=True), map.reset_index(drop=True)], axis = 1)
+
+    return map
+
+
+def compute_3D_map(df, xy_range, dt_range, xy_nbins, dt_nbins):
+    NaN_map = create_NaN_map(xy_range, dt_range, xy_nbins, dt_nbins)
+    map_3D_fit = map_3D_fits(df, xy_range = xy_range, dt_range = dt_range, xy_nbins = xy_nbins, dt_nbins = dt_nbins)
+
+    map = merge_maps(NaN_map, map_3D_fit)
+    full_map = include_coordinates(map, xy_range, dt_range, xy_nbins, dt_nbins)
+
+    return full_map
+
+#def compute_metadata():
+#    ['rmax', 'zmax', 'bin_size_z', 'bin_size_x', 'bin_size_y', 'method',
+#       'zbins', 'xbins', 'ybins', 'nbins_z', 'nbins_x', 'nbins_y', 'map_shape',
+#       'map_extent']
+
+#def save_map(map, metadata)
diff --git a/invisible_cities/icaros/krmap_functions_test.py b/invisible_cities/icaros/krmap_functions_test.py
new file mode 100644
index 000000000..29456fa6a
--- /dev/null
+++ b/invisible_cities/icaros/krmap_functions_test.py
@@ -0,0 +1,88 @@
+import numpy as np
+import pandas as pd
+from invisible_cities.core.core_functions import in_range, shift_to_bin_centers
+import itertools
+from krmap_functions import create_NaN_map, med_fun, fit_function, map_3D_fits, merge_maps, include_coordinates
+from pytest import raises
+from invisible_cities.core.core_functions import fix_random_seed
+
+
+def test_empty_medfun():
+    empty_dst = pd.DataFrame(columns = ['DT', 'x', 'y', 'S2e'])
+    result = med_fun(empty_dst)
+
+    assert result.shape == (1, 5)
+    assert result['nevents'].sum() == 0
+    for col in ['median', 'sigma', 'median_error', 'sigma_error']:
+        assert pd.isna(result[col].iloc[0])
+
+
+
+def test_empty_medfun_2():
+    empty_dst = pd.DataFrame(columns = ['DT', 'x', 'y', 'S2e'])
+
+    #with raises(ValueError):
+    med_fun(empty_dst)
+
+
+def test_fit_fun():
+    empty_dst = pd.DataFrame(columns = ['DT', 'x', 'y', 'S2e'])
+    result_fit = fit_function(empty_dst, 1)
+    result_fun = med_fun(empty_dst)
+
+    assert np.allclose(result_fit.mu.values, result_fun['median'].values, equal_nan=True)
+
+
+
+def test_Nevents():
+    d = {'DT': np.empty(6, dtype=int), 'x' : np.empty(6), 'y' : np.empty(6),'S2e' : np.ones(6)}
+    df_test = pd.DataFrame(data = d, index = range(0, 6))
+    N = len(df_test)
+    result_med = med_fun(df_test)
+    assert result_med['nevents'].iloc[0] == N
+    assert fit_function(df_test, bins = 6)['nevents'].iloc[0] == N
+
+
+def test_merge_maps():
+    Nan_map_test = create_NaN_map(xy_range = (-500, 500), dt_range = (0, 1400), xy_nbins = 100, dt_nbins = 10)
+    d = {'nevents': np.empty(6, dtype = int), 'mu' : np.empty(6), 'sigma' : np.empty(6), 'mu_error' : np.empty(6), 'sigma_error' : np.empty(6)}
+    map_3D_test = pd.DataFrame(data = d, index = range(0, 6))
+    assert merge_maps(Nan_map_test, map_3D_test).shape == Nan_map_test.shape
+
+test_merge_maps()
+
+
+def test_medfun_even():
+    d = {'DT': np.empty(6, dtype=int), 'x' : np.empty(6), 'y' : np.empty(6),'S2e' : [8000, 7500, 8300, 7900, 9000, 8100]}
+    df_test = pd.DataFrame(data = d, index = range(0,6))
+    S2e = df_test['S2e']
+    result_med_fun = med_fun(df_test)
+    result_med_fun = result_med_fun['median']
+    result_med_data = 8050
+
+    assert (result_med_fun == result_med_data).all()
+
+def test_medfun_odd():
+    d = {'DT': np.empty(7, dtype=int), 'x' : np.empty(7), 'y' : np.empty(7),'S2e' : [8000, 7500, 8300, 7900, 9000, 8100, 9100]}
+    df_test = pd.DataFrame(data = d, index = range(0,7))
+    S2e = df_test['S2e']
+    result_med_fun = med_fun(df_test)
+    result_med_fun = result_med_fun['median']
+    result_med_data = 8100
+    assert (result_med_fun == result_med_data).all()
+
+
+
+def test_fit_fun2():
+
+    with fix_random_seed(42):
+        S2e_df = pd.DataFrame({
+            'S2e': np.random.normal(loc = 8000, scale = 10.0, size = 10000)
+        })
+
+    results = fit_function(S2e_df, bins = 50)
+
+    assert np.isclose(results['mu'][0], S2e_df.S2e.mean(), atol = 1)
+    assert np.isclose(results['sigma'][0], S2e_df.S2e.std(), atol = 0.5)
+    assert np.isclose(results['mu_error'][0], S2e_df.S2e.std()/np.sqrt(len(S2e_df)), atol = 0.1)
+    #error de sigma?
diff --git a/invisible_cities/icaros/lifetime_vdrift_functions.py b/invisible_cities/icaros/lifetime_vdrift_functions.py
new file mode 100644
index 000000000..c91d087ea
--- /dev/null
+++ b/invisible_cities/icaros/lifetime_vdrift_functions.py
@@ -0,0 +1,137 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy import stats
+from matplotlib.patches import Circle, Rectangle
+from typing import Tuple
+from invisible_cities.core.stat_functions import poisson_sigma
+import invisible_cities.core.fit_functions as fit
+
+
+
+def select_lifetime_region(df, bins, x0, y0, r, half_size, circle=True):
+
+    df = df.dropna(subset=['X','Y','S2e', 'DT'])
+
+    xmin, xmax = df.X.min(), df.X.max()
+    ymin, ymax = df.Y.min(), df.Y.max()
+
+    mean, ebins, _  = stats.binned_statistic_dd(
+        (df.X, df.Y), df.S2e,
+        bins=[np.linspace(xmin, xmax, bins),
+              np.linspace(ymin, ymax, bins)],
+        statistic="mean"
+    )
+
+    x_centers = shift_to_bin_centers(ebins[0])
+    y_centers = shift_to_bin_centers(ebins[1])
+
+    Xc,Yc = np.meshgrid(x_centers, y_centers)
+
+
+    if circle == False:
+
+        x1, y1 = x0-half_size, y0-half_size
+        mask = ((Xc >= x1) & (Xc <= x1+2*half_size) & (Yc >= y1) & (Yc <= y1+2*half_size))
+        mean_values = mean.T[mask]
+
+
+        mask_df = ((df.X >= x1) & (df.X <= x1+2*half_size) & (df.Y >= y1) &
+                   (df.Y <= y1+2*half_size))
+
+        df_region = df[mask_df]
+
+    else:
+        mask = (Xc - x0)**2 + (Yc - y0)**2 <= r**2
+        mean_values = mean.T[mask]
+
+        mask_df = (df.X - x0)**2 + (df.Y - y0)**2 <= r**2
+        df_region = df[mask_df]
+
+
+    return df_region, mean_values
+
+
+
+def LT_fit(DT, s2e, p0):
+
+    def exponential(t, e, tau):
+        return e*np.exp(-t/tau)
+
+    popt, pcov = curve_fit(exponential, DT, s2e, p0 = p0)
+    magnitudes = (popt[0], popt[1])
+    uncertainties = (np.sqrt(pcov[0,0]), np.sqrt(pcov[1,1]))
+    return magnitudes, uncertainties
+
+
+
+def sigmoid(x          : np.array,
+            scale      : float,
+            inflection : float,
+            slope      : float,
+            offset     : float)->np.array:
+    '''
+    Sigmoid function, it computes the sigmoid of the input array x using the specified
+    parameters for scaling, inflection point, slope, and offset.
+    Parameters
+    ----------
+    x : np.array
+        The input array.
+    scale : float
+        The scaling factor determining the maximum value of the sigmoid function.
+    inflection : float
+        The x-value of the sigmoid's inflection point (where the function value is half of the scale).
+    slope : float
+        The slope parameter that controls the steepness of the sigmoid curve.
+    offset : float
+        The vertical offset added to the sigmoid function.
+    Returns
+    -------
+    np.array
+        Array of computed sigmoid values for x array.
+    '''
+    sigmoid = scale / (1 + np.exp(-slope * (x - inflection))) + offset
+    return sigmoid
+
+
+def compute_drift_v(dtdata   : np.array,
+                    nbins    : int,
+                    dtrange  : Tuple[float, float],
+                    seed     : Tuple[float, float, float, float]):
+    '''
+    Computes the drift velocity for a given distribution
+    using the sigmoid function to get the cathode edge.
+    Parameters
+    ----------
+    dtdata: array_like
+        Values of DT coordinate.
+    nbins: int (optional)
+        The number of bins in the z coordinate for the binned fit
+    dtrange: length-2 tuple (optional)
+        Fix the range in DT.
+    seed: length-4 tuple (optional)
+        Seed for the fit.
+    detector: string (optional)
+        Used to get the cathode position from DB.
+    Returns
+    -------
+    dv: float
+        Drift velocity.
+    dvu: float
+        Drift velocity uncertainty.
+    '''
+    y, x = np.histogram(dtdata, nbins, dtrange)
+    x    = shift_to_bin_centers(x)
+    if seed is None: seed = np.max(y), np.mean(dtrange), 0.5, np.min(y) # CHANGE: dtrange should be established from db
+    # At the moment there is not NEXT-100 DB so this won't work for that geometry
+    z_cathode = 1187.37
+    print(seed)
+    try:
+        f   = fit.fit(sigmoid, x, y, seed, sigma = poisson_sigma(y), fit_range = dtrange)
+        par = f.values
+        err = f.errors
+        dv  = (z_cathode + 10/2)/par[1]
+        dvu = dv/par[1] * err[1]
+    except RuntimeError:
+        print("WARNING: Sigmoid fit for dv computation fails. NaN value will be set in its place.")
+        dv, dvu = np.nan, np.nan
+    return dv, dvu
diff --git a/invisible_cities/icaros/selection_functions.py b/invisible_cities/icaros/selection_functions.py
new file mode 100644
index 000000000..226d1ca01
--- /dev/null
+++ b/invisible_cities/icaros/selection_functions.py
@@ -0,0 +1,148 @@
+import numpy             as np
+import pandas            as pd
+from   scipy.optimize     import curve_fit
+
+import matplotlib.pyplot as plt
+from invisible_cities.core.core_functions import in_range, shift_to_bin_centers
+from scipy import stats
+from apply_3Dmap import apply_3Dmap, normalization
+import itertools
+
+
+
+def eff_of_sel(sel, name = "", do_print = True):
+    nsel, ntot = np.sum(sel), len(sel)
+    eff = nsel/float(ntot)
+    if do_print:
+        seff = "{:.4f}".format(eff)
+        print(f"{name} efficiency {seff}, {nsel}/{ntot}.")
+    return eff
+
+
+dtrms2_low = lambda dt: -0.7 + 0.030 * (dt-20) # Gonzalo's
+dtrms2_upp = lambda dt: 2.6 + 0.036 * (dt-20) # Gonzalo'2
+dtrms2_cen = lambda dt:  1.0 + 0.033 * (dt-20)
+def dist_to_bandcenter(df): return df.Zrms**2 - dtrms2_cen(df.DT)
+
+
+def load_files(path_data, n):
+    filenames = []
+    for i in range(len(n)):
+        path_run =  f"/trigger1/ldc{i}"
+        path = path_data + path_run
+        file = [os.path.join(path, f) for f in os.listdir(path) if f.endswith(".h5")]
+        filenames.append(file)
+    filenames = list(itertools.chain.from_iterable(filenames))
+    return filenames
+
+
+
+def eff_of_selection(df_before, df_after, name = ""):
+    events_before = df_before.event.nunique()
+    events_after = df_after.event.nunique()
+
+    eff = events_after/events_before
+
+    print(f"{name} efficiency {eff}, {events_after}/{events_before}.")
+
+    return eff
+
+
+def apply_correctionmap(kdst, map3D, norm_method, keV = True):
+    corrected_energy = apply_3Dmap(map3D, norm_method, kdst.DT, kdst.X, kdst.Y, kdst.S2e, keV = keV)
+    kdst['Ec'] = corrected_energy.values
+
+    return kdst
+
+
+def select_diffusion_band(kdst, dtrms2_low, dtrms2_upp):
+
+    sel_DTband = in_range(kdst.Zrms**2, dtrms2_low(kdst.DT), dtrms2_upp(kdst.DT))
+    df_DTband = kdst[sel_DTband]
+
+    eff_DT = eff_of_selection(kdst, df_DTband, 'band selection')
+
+    return df_DTband, eff_DT
+
+
+def select_Xrays(kdst, low_xrays, high_xrays):
+
+    sel_xrays = (kdst.Ec >= low_xrays) & (kdst.Ec <= high_xrays)
+    df_Xrays = kdst[sel_xrays]
+
+    eff_Xrays = eff_of_selection(kdst, df_Xrays, 'remove xrays')
+
+    return df_Xrays, eff_Xrays
+
+
+def select_1S1_1S2(kdst):
+    a = kdst.groupby('event')['time'].count()
+    events_1S1_1S2 = a[a == 1].index.values
+    df_1S1_1S2 = kdst[kdst.event.isin(events_1S1_1S2)]
+
+    eff_1S1_1S2 = eff_of_selection(kdst, df_1S1_1S2, '1S1 & 1S2')
+
+    return df_1S1_1S2, eff_1S1_1S2
+
+
+def select_S2t(kdst, low_S2t, high_S2t):
+    range_S2t = (low_S2t, high_S2t)
+    sel_S2t = in_range(kdst.S2t, *range_S2t)
+
+    df_S2t = kdst[sel_S2t]
+
+    eff_S2t = eff_of_selection(kdst, df_S2t, "S2 in trigger time")
+
+    return df_S2t, eff_S2t
+
+
+def select_Rmax(kdst, R_max):
+    df_Rmax = kdst[kdst.R <= R_max]
+
+    eff_Rmax = eff_of_selection(kdst, df_Rmax, f'events with R less than {R_max}')
+
+    return df_Rmax, eff_Rmax
+
+
+def select_DTrange(kdst, low_DT, high_DT):
+    df_DTrange = kdst[(kdst.DT >= low_DT) & (kdst.DT <= high_DT)]
+
+    eff_DTrange = eff_of_selection(kdst, df_DTrange, f'events in DT range [{low_DT}, {high_DT}]')
+
+    return df_DTrange, eff_DTrange
+
+
+def select_nsipm(kdst, low_nsipm, high_nsipm):
+    sel_nsipm = (kdst.Nsipm >= low_nsipm) & (kdst.Nsipm <= high_nsipm)
+    df_nsipm = kdst[sel_nsipm]
+
+    eff_nsipm = eff_of_selection(kdst, df_nsipm, f'events in NSipm range [{low_nsipm}, {high_nsipm}]')
+
+    return df_nsipm, eff_nsipm
+
+
+def apply_selections(kdst, run_number,  dtrms2_low, dtrms2_upp, low_xrays, high_xrays, low_S2t, high_S2t, R_max, low_DT, high_DT, low_nsipm, high_nsipm):
+
+    df, eff_DTband = select_diffusion_band(kdst, dtrms2_low, dtrms2_upp)
+
+    df, eff_Xrays = select_Xrays(df, low_xrays, high_xrays)
+
+    df, eff_1S1_1S2 = select_1S1_1S2(df)
+
+    df, eff_S2t = select_S2t(df, low_S2t, high_S2t)
+
+    df, eff_Rmax = select_Rmax(df, R_max)
+
+    df, eff_DTrange = select_Dtrange(df, low_DT, high_DT)
+
+    df_final, eff_nsipm = select_nsipm(df, low_nsipm, high_nsipm)
+
+    total_efficiency = eff_of_selection(kdst, df_final, f'total events after all selections')
+
+    d = {'eff diffusion band': [eff_DT], 'eff X rays': [eff_Xrays], 'eff 1S1 & 1S2': [eff_1S1_1S2],
+         'eff S2 trigger time': [eff_S2t], 'eff Rmax': [eff_Rmax], 'eff range DT': [eff_DTrange],
+         'eff number of SiPMS': [eff_nsipm], 'total efficiency': [total_eff]}
+
+    df_efficiencies = pd.DataFrame(data = d)
+
+    return df_final, df_efficiencies