Updates on zemrude

invisible_cities/cities/zemrude.py

@@ -18,6 +18,7 @@
 
 from .. icaros  .correction_functions    import apply_correctionmap_inplace_kdst
 from .. icaros  .selection_functions     import apply_selections
+from .. icaros  .selection_functions     import discard_nan_values
 from .. icaros  .krmap_functions         import compute_3D_map
 from .. icaros  .krmap_functions         import gaussian_fit_ready
 from .. icaros  .krmap_functions         import get_median
@@ -36,14 +37,18 @@
 
 
 def concatenated_dsts_from_files(path: List[str], group: str, node:str)-> Iterator[Dict[str,Union[pd.DataFrame, int, np.ndarray]]]:
-    df = load_dsts(path, group, node)
+    df = load_dsts(path, group, node, ignore_errors = True)
     with tb.open_file(path[0], 'r') as h5in:
         run_number = get_run_number(h5in)
 
     yield dict(dst = df,
                run_number = run_number
                )
 
+def select_no_nan():
+    def discard_nans(df):
+        return discard_nan_values(df)
+    return discard_nans
 
 def apply_map(pre_map, norm_method, xy_params, col_name, unit):
     pre_map = pd.read_hdf(pre_map)
@@ -97,9 +102,9 @@ def save(efficiencies, krmap, metadata, t_evol):
     return save
 
 
-def do_control_plots(plots_out,ebins1, ns1bins, s1hbins, s1wbins, ebins2, ns2bins, s2hbins, s2qbins, qmaxbins, s2wbins, dtrms2_low, dtrms2_upp, drms2_cen, dtbins2, bins, dtrs2_bins, statistic, x0, y0, shape, shape_size, xy_range_plot):
+def do_control_plots(plots_out,ebins1, ns1bins, s1hbins, s1wbins, ebins2, ns2bins, s2hbins, s2qbins, qmaxbins, s2wbins, dtrms2_low, dtrms2_upp, drms2_cen, dtbins2, high_S2e, low_S2e, high_DT, low_DT, bins, dtrs2_bins, statistic, x0, y0, shape, shape_size, xy_range_plot):
     def control_plots(df, df_corr, efficiencies, run_number):
-        return make_control_plots(df, df_corr, efficiencies, run_number, plots_out, ebins1, ns1bins, s1hbins, s1wbins, ebins2, ns2bins, s2hbins, s2qbins, qmaxbins, s2wbins, dtrms2_low, dtrms2_upp, drms2_cen,dtbins2, bins, dtrs2_bins, statistic, x0, y0, shape, shape_size, xy_range_plot)
+        return make_control_plots(df, df_corr, efficiencies, run_number, plots_out, ebins1, ns1bins, s1hbins, s1wbins, ebins2, ns2bins, s2hbins, s2qbins, qmaxbins, s2wbins, dtrms2_low, dtrms2_upp, drms2_cen,dtbins2,high_S2e, low_S2e, high_DT, low_DT, bins, dtrs2_bins, statistic, x0, y0, shape, shape_size, xy_range_plot)
     return control_plots
 
 
@@ -158,6 +163,8 @@ def zemrude(files_in           : OneOrManyFiles
             , qmaxbins         : np.ndarray
             , s2wbins          : np.ndarray
             , dtbins2          : np.ndarray
+            , high_S2e         : float
+            , low_S2e          : float
             , bins             : int
             , dtr2_bins        : tuple
             , statistic        : str
@@ -166,12 +173,16 @@ def zemrude(files_in           : OneOrManyFiles
             , xy_params        : dict = None
             ):
 
+    apply_nan_cut          = fl.map( select_no_nan()
+                                     ,item  = 'dst'
+                                     )
+
     apply_preliminary_map  = fl.map( apply_map(pre_map,
                                               norm_method,
                                               xy_params,
                                               'Ec',
                                               unit = keV)
-                                  , item = 'dst')
+                                  , item  = 'dst')
 
 
     apply_selections = fl.map( select_dst(dtrms2_low,
@@ -248,6 +259,10 @@ def zemrude(files_in           : OneOrManyFiles
                                                   dtrms2_upp,
                                                   dtrms2_cen,
                                                   dtbins2,
+                                                  high_S2e,
+                                                  low_S2e,
+                                                  high_DT,
+                                                  low_DT,
                                                   bins,
                                                   dtr2_bins,
                                                   statistic,
@@ -269,7 +284,8 @@ def zemrude(files_in           : OneOrManyFiles
     with tb.open_file(file_out, "w", filters=tbl.filters(compression)):
         pass
     fl.push( source = concatenated_dsts_from_files(files_in, "DST", "Events")
-            ,pipe   = fl.pipe(apply_preliminary_map,
+            ,pipe   = fl.pipe( apply_nan_cut,
+                               apply_preliminary_map,
                                apply_selections,
                                compute_3D_map,
                                compute_metadata,

invisible_cities/config/zemrude.conf

@@ -60,6 +60,10 @@ s2qbins       = np.linspace(0, 1.5e3, 100)
 qmaxbins      = np.linspace(0, 300, 100)
 s2wbins       = np.linspace(0, 60, 100)
 dtbins2       = np.linspace(0, 1400, 51)
+high_S2e      = 1e4
+low_S2e       = 4e3
+high_DT       = 1350
+low_DT        = 20
 bins          = 100
 dtr2_bins     = (20, 20)
 statistic     = 'mean'

invisible_cities/config/zemrude_mc.conf

@@ -0,0 +1,70 @@
+files_in    = '/Users/maria/NEXT/NEXT-100/data/MC/kr/fullsim/202608/sophronia/*.h5'
+file_out    = 'MC_3Dmap_fullsim_202608_50bins.h5'
+compression = 'ZLIB4'
+event_range = all
+
+# run number 0 is for MC
+run_number = 0
+
+# Detector database used
+detector_db = 'next100'
+
+pre_map     = '/Users/maria/IC/invisible_cities/database/test_data/preliminary_map_15502.h5'
+norm_method = maximum
+
+# ------ Selection parameters ------
+dtrms2_low = lambda dt: -0.7 * np.ones_like(dt) #less than zrms2.min()
+dtrms2_upp = lambda dt: 999998000002.0 * np.ones_like(dt) #greater than zrms2.max()
+dtrms2_cen = lambda dt: 0.0 * (dt-20)
+low_xrays  = 0
+high_xrays = 1e6
+low_S2t    = 0
+high_S2t   = 1e12
+R_max      = 500
+low_DT     = 20
+high_DT    = 1350
+low_nsipm  = 0
+high_nsipm = 1000
+
+# ------ Create map parameters ------
+xy_range     = (-500, 500)
+dt_range     = (20, 1350)
+xy_nbins     = 100
+dt_nbins     = 10
+S2e_range    = (1000, 20000)
+fit_function = gaussian
+min_events   = 40
+nbins        = 100
+
+# ------ Time evolution parameters ------
+slice_hours = 1e12
+x0          = 0
+y0          = 0
+shape       = circle
+shape_size  = 100
+dtbins_dv   = np.linspace(1200, 1400)
+s1_DTrange  = (1000, 1350)
+bins_Ec     = np.linspace(20, 60, 101)
+error       = False
+
+# ------ Make control plots parameters ------
+plots_out     = f'plots_out_{run_number}'
+ebins1        = np.linspace(0, 50, 101)
+ns1bins       = np.linspace(0, 10, 10)
+s1hbins       = np.linspace(0, 20, 100)
+s1wbins       = np.linspace(50, 1000, 20)
+ebins2        = np.arange(4e3, 1.4e4, 51)
+ns2bins       = np.linspace(0, 20, 20)
+s2hbins       = np.linspace(0, 3e3, 100)
+s2qbins       = np.linspace(0, 1.5e3, 100)
+qmaxbins      = np.linspace(0, 300, 100)
+s2wbins       = np.linspace(0, 60, 100)
+dtbins2       = np.linspace(0, 1400, 51)
+high_S2e      = 1.4e4
+low_S2e       = 4e3
+high_DT       = 1350
+low_DT        = 20
+bins          = 100
+dtr2_bins     = (20, 20)
+statistic     = 'mean'
+xy_range_plot = np.linspace(-500, 500, 100)

invisible_cities/icaros/control_plots_functions.py

@@ -330,26 +330,26 @@ def monitor_dtime(df          : pd.DataFrame,
     axs[0,0].plot(df1.DT, dtrms2_low(df1.DT), ".r", ms=2);
     axs[0,0].plot(df1.DT, dtrms2_upp(df1.DT), ".r", ms=2);
     axs[0,0].plot(df1.DT, dtrms2_cen(df1.DT), '.g', ms = 2);
-    axs[0,0].set_xlabel("Drift time ($\mu$s)"); axs[0,0].set_ylabel("DT$_{rms}^2$ ($\mu$s)"); axs[0,0].set_xlim(0, 1300)
+    axs[0,0].set_xlabel(r"Drift time ($\mu$s)"); axs[0,0].set_ylabel(r"DT$_{rms}^2$ ($\mu$s)"); axs[0,0].set_xlim(0, 1300)
     axs[0, 0].set_title('Before selection')
 
     axs[0,1].hist2d(df2.DT, df2.Zrms**2, (dtbins, dtrms2bins));
     axs[0,1].plot(df2.DT, dtrms2_low(df2.DT), ".r", ms=2);
     axs[0,1].plot(df2.DT, dtrms2_upp(df2.DT), ".r", ms=2);
     axs[0,1].plot(df2.DT, dtrms2_cen(df2.DT), '.g', ms = 2);
-    axs[0,1].set_xlabel("Drift time ($\mu$s)"); axs[0,1].set_ylabel("DT$_{rms}^2$ ($\mu$s)"); axs[0,1].set_xlim(0, 1300)
+    axs[0,1].set_xlabel(r"Drift time ($\mu$s)"); axs[0,1].set_ylabel(r"DT$_{rms}^2$ ($\mu$s)"); axs[0,1].set_xlim(0, 1300)
     axs[0,1].set_title('After selection')
 
     axs[1,0].hist(df1_.DT, dtbins, histtype = 'step', color = 'mediumpurple',lw = 2, label = 'before selection');
     axs[1,0].hist(df2.DT, dtbins, histtype = 'step', color = 'black', lw = 2, label = 'after selection');
     axs[1,0].legend();
-    axs[1,0].set_xlabel("Drift time ($\mu$s)");
+    axs[1,0].set_xlabel(r"Drift time ($\mu$s)");
     axs[1,0].grid(True)
 
     axs[1,1].hist(df1_.Zrms**2, 100, (0, 40), histtype = 'step',color = 'mediumpurple',lw = 2, label = 'before selection');
     axs[1,1].hist(df2.Zrms**2, 100, (0, 40), histtype = 'step', color = 'black',lw = 2, label = 'after selection');
     axs[1,1].legend();
-    axs[1,1].set_xlabel("DT$_{rms}^2$ ($\mu$s)");
+    axs[1,1].set_xlabel(r"DT$_{rms}^2$ ($\mu$s)");
     axs[1,1].grid(True)
 
     fig.tight_layout();
@@ -370,26 +370,30 @@ def monitor_lifetime(df        : pd.DataFrame,
 
 
 def monitor_kr_distribution(df        : pd.DataFrame,
+                            high_S2e  : float,
+                            low_S2e   : float,
+                            high_DT   : float,
+                            low_DT    : float,
                             bins      : int,
                             dtr2_bins : tuple):
     """
     Plots the square radial distribution and a 2D distribution of the
     square radius as a function of drift time.
     """
-    sel = in_range(df.S2e, 7.5e3, 9.5e3) & in_range(df.DT, 20, 1350)
+    sel = in_range(df.S2e, low_S2e, high_S2e) & in_range(df.DT, low_DT, high_DT)
 
     DT = (df.DT[sel]).dropna()
     R2 = (df.X[sel]**2 + df.Y[sel]**2).dropna()
 
     fig, axs = plt.subplots(1, 2, figsize = (21, 7))
 
     axs[0].hist(R2, bins, histtype = 'step', color = 'mediumpurple', lw = 2);
-    axs[0].set_xlabel("R$^2$ (mm$^2$)"); freq();
+    axs[0].set_xlabel(r"R$^2$ (mm$^2$)"); freq();
     axs[0].grid(True)
 
     axs[1].hist2d(DT, R2, dtr2_bins);
-    axs[1].set_xlabel("DT ($\mu$s)");
-    axs[1].set_ylabel("R$^2$ (mm$^2$)");
+    axs[1].set_xlabel(r"DT ($\mu$s)");
+    axs[1].set_ylabel(r"R$^2$ (mm$^2$)");
 
 
 def hist2D(df        : pd.DataFrame,
@@ -465,14 +469,14 @@ def plot_Ec(Ec   : pd.core.series.Series,
              label = f'mean Ec: {mean_Ec:.2f}keV\n'
              f'median Ec: {median_Ec:.2f}keV\n'
              f'std Ec: {stdEc:.2f}keV\n'
-             f'umean Ec: {umeanEc:.2f}keV'
+             f'umean Ec: {umeanEc:.4f}keV'
              )
 
     axs.hist(Ec_2, 100, range = (25, 60), histtype = 'step', color = 'mediumpurple',lw = 2,
              label = f'mean Ec_2: {mean_Ec2:.2f}keV\n'
              f'median Ec_2: {median_Ec2:.2f}keV\n'
              f'std Ec_2: {stdEc2:.2f}keV\n'
-             f'umean Ec_2: {umeanEc2:.2f}keV'
+             f'umean Ec_2: {umeanEc2:.4f}keV'
              )
     axs.set_xlabel('Ec (keV)'); freq();
     axs.grid()
@@ -516,7 +520,7 @@ def plot_lifetime_fit(df         : pd.DataFrame,
                     f'u_lifetime : {uncertainties[1]:.2f}'
              )
     axs.errorbar(dt, e, yerr = se, fmt = '.');
-    axs.set_ylim(6000, 10000);
+    #axs.set_ylim(6000, 14000);
     axs.legend();
 
 
@@ -592,7 +596,7 @@ def plot_efficiencies(efficiencies : pd.DataFrame):
     for i in range(ncuts):
         plt.scatter(efficiencies.columns[i], efficiencies_plot[i],
                     marker = 'x', color = colors[i], lw = 2,
-                    s = 100, label = f'{efficiencies.columns[i]}: {efficiencies_plot[i]:.2f}')
+                    s = 100, label = f'{efficiencies.columns[i]}: {efficiencies_plot[i]:.3f}')
     plt.xticks(rotation=20, ha='right')
     plt.grid()
     plt.legend()
@@ -722,6 +726,10 @@ def make_control_plots(df               : pd.DataFrame,
                        dtrms2_upp       : Callable,
                        dtrms2_cen       : Callable,
                        dtbins2          : np.array,
+                       high_S2e         : float,
+                       low_S2e          : float,
+                       high_DT          : float,
+                       low_DT           : float,
                        bins             : int,
                        dtr2_bins        : tuple,
                        statistic        : str,
@@ -749,7 +757,7 @@ def make_control_plots(df               : pd.DataFrame,
     plt.gcf().savefig(f"{plots_out}/monitor_lifetime_run{run_number}.png")
     plt.close()
 
-    monitor_kr_distribution(df, bins, dtr2_bins)
+    monitor_kr_distribution(df, high_S2e, low_S2e, high_DT, low_DT,  bins, dtr2_bins)
     plt.gcf().savefig(f"{plots_out}/monitor_kr_distribution_run{run_number}.png")
     plt.close()
 

invisible_cities/icaros/krmap_functions.py

@@ -725,3 +725,16 @@ def save_map(name          : str,
         df_writer(file, efficiencies, group_name = 'data', table_name = 'selection_efficiencies')
         df_writer(file, krmap, group_name = 'krmap', table_name = 'krmap')
         df_writer(file, t_evol, group_name = 't_evol', table_name = 't_evol')
+
+def merge_multiple_maps(maps : pd.DataFrame) -> pd.DataFrame:
+    """
+    Merge multiple maps into a single map.
+    Group maps by their (k i j) indices. The output 'mu' (new_mu) is computed as the weighted mean of the input 'mu' values using 'nevents' as weights.
+    """
+    def merge_bins(df):
+        EPSILON = np.finfo(np.float64).eps
+        new_mu  = np.sum(df.mu * df.nevents) / (df.nevents.sum() + EPSILON)
+        return pd.DataFrame( dict( mu = new_mu
+                             , nevents  = df.nevents.sum())
+                       , index=[0])
+    return maps.groupby('k i j'.split()).apply(merge_bins).reset_index()

invisible_cities/icaros/selection_functions.py

@@ -32,6 +32,26 @@ def eff_of_selection(df_before  : pd.DataFrame,
 
     return eff
 
+def discard_nan_values(kdst : pd.DataFrame) -> pd.DataFrame :
+
+    """
+    Discard every kdst row that has any NaN value in DT, X or Y.
+    Parameters
+    ----------
+    kdst : pd.DataFrame
+       Dataframe (kdst sophronia output) to clean from NaN values.
+    Returns
+    -------
+    kdst_nan : pd.DataFrame
+       Dataframe with no NaNs
+    """
+
+    kdst     = kdst[~(kdst.DT.isna())]
+    kdst     = kdst[~(kdst.X.isna())]
+    kdst_nan = kdst[~(kdst.Y.isna())]
+
+    return kdst_nan
+
 
 def select_var_inrange(kdst      : pd.DataFrame,
                        col_name  : str,