From 65db8e84b4db3da6230775124dda0948f6058d6f Mon Sep 17 00:00:00 2001 From: carhc Date: Fri, 27 Sep 2024 02:21:09 +0200 Subject: [PATCH 01/48] Add sigmoid function --- invisible_cities/reco/krmap_evolution.py | 34 ++++++++++++++++++++++++ 1 file changed, 34 insertions(+) create mode 100644 invisible_cities/reco/krmap_evolution.py diff --git a/invisible_cities/reco/krmap_evolution.py b/invisible_cities/reco/krmap_evolution.py new file mode 100644 index 0000000000..e7a034c090 --- /dev/null +++ b/invisible_cities/reco/krmap_evolution.py @@ -0,0 +1,34 @@ +import numpy as np + + +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 From 6e91528d2bce68de586a5bc416fd571e86214f28 Mon Sep 17 00:00:00 2001 From: carhc Date: Fri, 27 Sep 2024 02:22:55 +0200 Subject: [PATCH 02/48] Add gauss_seed function --- invisible_cities/reco/krmap_evolution.py | 31 ++++++++++++++++++++++++ 1 file changed, 31 insertions(+) diff --git a/invisible_cities/reco/krmap_evolution.py b/invisible_cities/reco/krmap_evolution.py index e7a034c090..12b09ffdeb 100644 --- a/invisible_cities/reco/krmap_evolution.py +++ b/invisible_cities/reco/krmap_evolution.py @@ -32,3 +32,34 @@ def sigmoid(x : np.array, sigmoid = scale / (1 + np.exp(-slope * (x - inflection))) + offset return sigmoid + + +def gauss_seed(x : np.array, + y : np.array, + sigma_rel : Optional[int] = 0.05): + + ''' + This function estimates the seed for a gaussian fit. + + Parameters + ---------- + x: np.array + Data to fit. + y: int + Number of bins for the histogram. + sigma_rel (Optional): int + Relative error, default 5%. + + Returns + ------- + seed: List + List with the seed estimation. + ''' + + y_max = np.argmax(y) + x_max = x[y_max] + sigma = sigma_rel * x_max + amp = y_max * (2 * np.pi)**0.5 * sigma * np.diff(x)[0] + seed = amp, x_max, sigma + + return seed From e445e30b756f79244106d1febebc78f9156c5f86 Mon Sep 17 00:00:00 2001 From: carhc Date: Fri, 27 Sep 2024 02:25:46 +0200 Subject: [PATCH 03/48] Add resolution function --- invisible_cities/reco/krmap_evolution.py | 30 ++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/invisible_cities/reco/krmap_evolution.py b/invisible_cities/reco/krmap_evolution.py index 12b09ffdeb..0ab858f2e3 100644 --- a/invisible_cities/reco/krmap_evolution.py +++ b/invisible_cities/reco/krmap_evolution.py @@ -63,3 +63,33 @@ def gauss_seed(x : np.array, seed = amp, x_max, sigma return seed + + +def resolution(values : np.array, + errors : np.array): + + ''' + Computes the resolution (FWHM) from the Gaussian parameters. + + Parameters + ---------- + values: np.array + Gaussian parameters: amplitude, center, and sigma. + errors: np.array + Uncertainties for the Gaussian parmeters. + + Returns + ------- + res: float + Resolution. + ures: float + Uncertainty of resolution. + ''' + + amp , mu, sigma = values + u_amp, u_mu, u_sigma = errors + + res = 235.48 * sigma/mu + ures = res * (u_mu**2/mu**2 + u_sigma**2/sigma**2)**0.5 + + return res, ures From a12fcc7f1007314199298525f5f0da65f4f99a3c Mon Sep 17 00:00:00 2001 From: carhc Date: Fri, 27 Sep 2024 02:34:40 +0200 Subject: [PATCH 04/48] Add quick_gauss_fit function --- invisible_cities/reco/krmap_evolution.py | 35 ++++++++++++++++++++++-- 1 file changed, 33 insertions(+), 2 deletions(-) diff --git a/invisible_cities/reco/krmap_evolution.py b/invisible_cities/reco/krmap_evolution.py index 0ab858f2e3..f71656b3cb 100644 --- a/invisible_cities/reco/krmap_evolution.py +++ b/invisible_cities/reco/krmap_evolution.py @@ -1,5 +1,8 @@ import numpy as np +from .. core.fit_functions import fit, gauss +from .. core.core_functions import shift_to_bin_centers + def sigmoid(x : np.array, scale : float, @@ -52,8 +55,8 @@ def gauss_seed(x : np.array, Returns ------- - seed: List - List with the seed estimation. + seed: Tuple + Tuple with the seed estimation. ''' y_max = np.argmax(y) @@ -93,3 +96,31 @@ def resolution(values : np.array, ures = res * (u_mu**2/mu**2 + u_sigma**2/sigma**2)**0.5 return res, ures + + +def quick_gauss_fit(data : np.array, + bins : int): + + ''' + This function histograms input data and then fits it to a Gaussian. + + Parameters + ---------- + data: np.array + Data to fit. + bins: int + Number of bins for the histogram. + + Returns + ------- + fit_output: FitFunction + Object containing the fit results + ''' + + y, x = np.histogram(data, bins) + x = shift_to_bin_centers(x) + seed = gauss_seed(x, y) + + fit_result = fit(gauss, x, y, seed) + + return fit_result From dbc4b41cdb6d0f601b99062c79ca2ff70f03c707 Mon Sep 17 00:00:00 2001 From: carhc Date: Fri, 27 Sep 2024 02:37:14 +0200 Subject: [PATCH 05/48] Add get_number_of_time_bins function --- invisible_cities/reco/krmap_evolution.py | 29 ++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/invisible_cities/reco/krmap_evolution.py b/invisible_cities/reco/krmap_evolution.py index f71656b3cb..a831e5da9e 100644 --- a/invisible_cities/reco/krmap_evolution.py +++ b/invisible_cities/reco/krmap_evolution.py @@ -124,3 +124,32 @@ def quick_gauss_fit(data : np.array, fit_result = fit(gauss, x, y, seed) return fit_result + + +def get_number_of_time_bins(nStimeprofile : int, + tstart : int, + tfinal : int)->int: + + ''' + Computes the number of time bins to use for a given time step + in seconds. + + Parameters + ---------- + nStimeprofile: int + Time step in seconds. + tstart: int + Initial timestamp for the dataset. + tfinal: int + Final timestamp for the dataset. + + Returns + ------- + ntimebins: int + Number of time bins. + ''' + + ntimebins = int(np.floor((tfinal - tstart) / nStimeprofile)) + ntimebins = np.max([ntimebins, 1]) + + return ntimebins From 755d843a0700c8661ba7e3abbc2f54333393ab0a Mon Sep 17 00:00:00 2001 From: carhc Date: Fri, 27 Sep 2024 02:41:40 +0200 Subject: [PATCH 06/48] Add get_time_series_df function --- invisible_cities/reco/krmap_evolution.py | 36 +++++++++++++++++++++++- 1 file changed, 35 insertions(+), 1 deletion(-) diff --git a/invisible_cities/reco/krmap_evolution.py b/invisible_cities/reco/krmap_evolution.py index a831e5da9e..2ab932cf0d 100644 --- a/invisible_cities/reco/krmap_evolution.py +++ b/invisible_cities/reco/krmap_evolution.py @@ -1,7 +1,10 @@ import numpy as np +from typing import List, Tuple, Callable, Optional +from pandas import DataFrame + from .. core.fit_functions import fit, gauss -from .. core.core_functions import shift_to_bin_centers +from .. core.core_functions import in_range, shift_to_bin_centers def sigmoid(x : np.array, @@ -153,3 +156,34 @@ def get_number_of_time_bins(nStimeprofile : int, ntimebins = np.max([ntimebins, 1]) return ntimebins + + +def get_time_series_df(ntimebins : int, + time_range : Tuple[float, float], + dst : DataFrame)->Tuple[np.array, List[np.array]]: + + ''' + Given a dst this function returns a time series (ts) and a list of masks which are used to divide + the event in time intervals. + + Parameters + ---------- + ntimebins : int + Number of time bins + time_range : Tuple + Time range + dst : pd.DataFrame + DataFrame + + Returns + ------- + A Tuple with: + np.array : The time series + List[np.array] : The list of masks to get the events for each time series. + ''' + + modified_right_limit = np.nextafter(time_range[-1], np.inf) + time_bins = np.linspace(time_range[0], modified_right_limit, ntimebins+1) + masks = np.array([in_range(dst['time'].to_numpy(), time_bins[i], time_bins[i + 1]) for i in range(ntimebins)]) + + return shift_to_bin_centers(time_bins), masks From acad63eb418b90c0e5f7ce4c68105aead491cddb Mon Sep 17 00:00:00 2001 From: carhc Date: Fri, 27 Sep 2024 02:47:47 +0200 Subject: [PATCH 07/48] Add compute_drift_v function --- invisible_cities/reco/krmap_evolution.py | 59 +++++++++++++++++++++++- 1 file changed, 58 insertions(+), 1 deletion(-) diff --git a/invisible_cities/reco/krmap_evolution.py b/invisible_cities/reco/krmap_evolution.py index 2ab932cf0d..8533d9e523 100644 --- a/invisible_cities/reco/krmap_evolution.py +++ b/invisible_cities/reco/krmap_evolution.py @@ -1,10 +1,12 @@ import numpy as np -from typing import List, Tuple, Callable, Optional + +from typing import List, Tuple, Optional from pandas import DataFrame from .. core.fit_functions import fit, gauss from .. core.core_functions import in_range, shift_to_bin_centers +from .. database import load_db as DB def sigmoid(x : np.array, @@ -187,3 +189,58 @@ def get_time_series_df(ntimebins : int, masks = np.array([in_range(dst['time'].to_numpy(), time_bins[i], time_bins[i + 1]) for i in range(ntimebins)]) return shift_to_bin_centers(time_bins), masks + + +def compute_drift_v(zdata : np.array, + nbins : int, + zrange : Tuple[float, float], + seed : Tuple[float, float, float, float], + detector : str)->Tuple[float, float]: + + ''' + Computes the drift velocity for a given distribution + using the sigmoid function to get the cathode edge. + + Parameters + ---------- + zdata: array_like + Values of Z coordinate. + nbins: int (optional) + The number of bins in the z coordinate for the binned fit. + zrange: length-2 tuple (optional) + Fix the range in z. + 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(zdata, nbins, zrange) + x = shift_to_bin_centers(x) + + if seed is None: seed = np.max(y), np.mean(zrange), 0.5, np.min(y) + + # At the moment there is not NEXT-100 DB so this won't work for that geometry + z_cathode = DB.DetectorGeo(detector).ZMAX[0] + + try: + f = fit(sigmoid, x, y, seed, sigma = poisson_sigma(y), fit_range = zrange) + + par = f.values + err = f.errors + + dv = z_cathode/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 From 4fbeb39c59ec5852c3de6e3d607de6506c3113dd Mon Sep 17 00:00:00 2001 From: carhc Date: Fri, 27 Sep 2024 02:57:47 +0200 Subject: [PATCH 08/48] Add e0_xy_correction function --- invisible_cities/reco/krmap_evolution.py | 37 ++++++++++++++++++++++-- 1 file changed, 35 insertions(+), 2 deletions(-) diff --git a/invisible_cities/reco/krmap_evolution.py b/invisible_cities/reco/krmap_evolution.py index 8533d9e523..698183d8b7 100644 --- a/invisible_cities/reco/krmap_evolution.py +++ b/invisible_cities/reco/krmap_evolution.py @@ -1,11 +1,16 @@ import numpy as np +import pandas as pd - -from typing import List, Tuple, Optional +from typing import List, Tuple, Optional, Callable from pandas import DataFrame +from .. types.symbols import NormStrategy from .. core.fit_functions import fit, gauss from .. core.core_functions import in_range, shift_to_bin_centers +from .. reco.corrections import get_normalization_factor +from .. reco.corrections import correct_geometry_ +from .. reco.corrections import maps_coefficient_getter +from .. core.stat_functions import poisson_sigma from .. database import load_db as DB @@ -244,3 +249,31 @@ def compute_drift_v(zdata : np.array, dv, dvu = np.nan, np.nan return dv, dvu + + +def e0_xy_correction(map : pd.DataFrame, + norm_strat : NormStrategy)->Callable: + + ''' + Provides the function to compute only the geometrical corrections. + + Parameters + ---------- + map: pd.DataFrame + Map containing the corrections. + norm_strat: + Normalization strategy used when correcting the energy. + + Returns + ------- + A function to compute geometrical corrections given a hit's X,Y position. + ''' + + normalization = get_normalization_factor(map , norm_strat) + get_xy_corr_fun = maps_coefficient_getter (map.mapinfo, map.e0) + + def geo_correction_factor(x : np.array, + y : np.array): + return correct_geometry_(get_xy_corr_fun(x,y))*normalization + + return geo_correction_factor From f1f83e14b03bf88aff01567245623221716e44ec Mon Sep 17 00:00:00 2001 From: carhc Date: Fri, 27 Sep 2024 03:10:24 +0200 Subject: [PATCH 09/48] Add computing_kr_parameters function --- invisible_cities/reco/krmap_evolution.py | 151 +++++++++++++++++++++-- 1 file changed, 139 insertions(+), 12 deletions(-) diff --git a/invisible_cities/reco/krmap_evolution.py b/invisible_cities/reco/krmap_evolution.py index 698183d8b7..4f2ec40dce 100644 --- a/invisible_cities/reco/krmap_evolution.py +++ b/invisible_cities/reco/krmap_evolution.py @@ -1,18 +1,20 @@ import numpy as np import pandas as pd -from typing import List, Tuple, Optional, Callable -from pandas import DataFrame - -from .. types.symbols import NormStrategy -from .. core.fit_functions import fit, gauss -from .. core.core_functions import in_range, shift_to_bin_centers -from .. reco.corrections import get_normalization_factor -from .. reco.corrections import correct_geometry_ -from .. reco.corrections import maps_coefficient_getter -from .. core.stat_functions import poisson_sigma -from .. database import load_db as DB - +from typing import List, Tuple, Optional, Callable +from pandas import DataFrame + +from .. types.symbols import NormStrategy +from .. types.symbols import KrFitFunction # Won't work until previous PR are approved +from .. core.fit_functions import fit, gauss +from .. core.core_functions import in_range, shift_to_bin_centers +from .. reco.corrections import get_normalization_factor +from .. reco.corrections import correct_geometry_ +from .. reco.corrections import maps_coefficient_getter +from .. reco.corrections import apply_all_correction +from .. core.stat_functions import poisson_sigma +from .. database import load_db as DB +from .. reco.icaro_components import get_fit_function_lt def sigmoid(x : np.array, scale : float, @@ -277,3 +279,128 @@ def geo_correction_factor(x : np.array, return correct_geometry_(get_xy_corr_fun(x,y))*normalization return geo_correction_factor + + +def computing_kr_parameters(data : DataFrame, + ts : float, + emaps : pd.DataFrame, + fittype : KrFitFunction, + nbins_dv : int, + zrange_dv : List[float, float], + detector : str, + norm_strat : NormStrategy, + norm_value : float)->DataFrame: # REVISAR NORM_STRAT Y NORM_VALUE + + ''' + Computes some average parameters (e0, lt, drift v, energy + resolution, S1w, S1h, S1e, S2w, S2h, S2e, S2q, Nsipm, 'Xrms, Yrms) + for a given krypton distribution. + + Parameters + ---------- + data: DataFrame + Kdst distribution to analyze. + ts: float + Central time of the distribution. + emaps: correction map + Allows geometrical correction of the energy. + fittype: KrFitFunction + Kind of fit to perform + nbins_dv: int + Number of bins in Z-coordinate to use in the histogram for the + drift velocity calculation. + zrange_dv: List + Range in Z-coordinate to use in the histogram for the drift + velocity calculation. + detector: string + Used to get the cathode position from DB for the drift velocity + computation. + norm_strat: NormStrategy + Normalization strategy to follow. + norm_value: float + Energy value to normalize to. + + Returns + ------- + pars: DataFrame + Each column corresponds to the average value of a different parameter. + ''' + + # Computing E0, LT + + geo_correction_factor = e0_xy_correction(map = emaps, + norm_strat = norm_strat) # PREGUNTAR POR ESTRATEGIA + + fit_func, seed = get_fit_function_lt(fittype) + + x = data.DT, + y = data.S2e.to_numpy()*geo_correction_factor(data.X.to_numpy(), data.Y.to_numpy()) + + fit_output, _, _, _ = fit(func = fit_func, # Misma funcion que en el ajuste del mapa + x = x, + y = y, + seed = seed(x, y), + full_output = False) + + e0, lt = fit_output.values + e0u, ltu = fit_output.errors + + # Computing Drift Velocity + + dv, dvu = compute_drift_v(zdata = data.Z.to_numpy(), + nbins = nbins_dv, + zrange = zrange_dv, + detector = detector) + + # Computing Resolution + + tot_corr_factor = apply_all_correction(maps = emaps, + apply_temp = False, + norm_strat = norm_strat, + norm_value = norm_value) + + nbins = int((len(data.S2e))**0.5) # Binning as a function of nevts. Should we change it? + + f = quick_gauss_fit(data.S2e.to_numpy()*tot_corr_factor(data.X. to_numpy(), + data.Y. to_numpy(), + data.Z. to_numpy(), + data.time.to_numpy()), + bins = nbins) + + par, err = f.values, f.errors + res, err_res = resolution(values = par, + errors = err) + + # Averages of parameters + + parameters = ['S1w', 'S1h', 'S1e', + 'S2w', 'S2h', 'S2e', 'S2q', + 'Nsipm', 'Xrms', 'Yrms'] + + mean_d, var_d = {}, {} + + for parameter in parameters: + + data_value = getattr(data, parameter) + mean_d[parameter] = np.mean(data_value) + var_d [parameter] = (np.var(data_value)/len(data_value))**0.5 + + # Creating parameter evolution table + + evol = DataFrame({'ts' : [ts] , + 'e0' : [e0] , 'e0u' : [e0u] , + 'lt' : [lt] , 'ltu' : [ltu] , + 'dv' : [dv] , 'dvu' : [dvu] , + 'resol': [res] , 'resolu': [err_res] , + 's1w' : [mean_d['S1w']] , 's1wu' : [var_d['S1w']] , + 's1h' : [mean_d['S1h']] , 's1hu' : [var_d['S1h']] , + 's1e' : [mean_d['S1e']] , 's1eu' : [var_d['S1e']] , + 's2w' : [mean_d['S2w']] , 's2wu' : [var_d['S2w']] , + 's2h' : [mean_d['S2h']] , 's2hu' : [var_d['S2h']] , + 's2e' : [mean_d['S2e']] , 's2eu' : [var_d['S2e']] , + 's2q' : [mean_d['S2q']] , 's2qu' : [var_d['S2q']] , + 'Nsipm': [mean_d['Nsipm']], 'Nsipmu': [var_d['Nsipm']], + 'Xrms' : [mean_d['Xrms']] , 'Xrmsu' : [var_d['Xrms']] , + 'Yrms' : [mean_d['Yrms']] , 'Yrmsu' : [var_d['Yrms']]}) + + return evol From 2dbef03a1a27e02852c0ca70fecab87e4d5d478e Mon Sep 17 00:00:00 2001 From: carhc Date: Fri, 27 Sep 2024 03:19:27 +0200 Subject: [PATCH 10/48] Add kr_time_evolution function --- invisible_cities/reco/krmap_evolution.py | 73 ++++++++++++++++++++++++ 1 file changed, 73 insertions(+) diff --git a/invisible_cities/reco/krmap_evolution.py b/invisible_cities/reco/krmap_evolution.py index 4f2ec40dce..6eb2993da8 100644 --- a/invisible_cities/reco/krmap_evolution.py +++ b/invisible_cities/reco/krmap_evolution.py @@ -404,3 +404,76 @@ def computing_kr_parameters(data : DataFrame, 'Yrms' : [mean_d['Yrms']] , 'Yrmsu' : [var_d['Yrms']]}) return evol + + +def kr_time_evolution(ts : np.array[float], + masks_time : List[np.array], + dst : pd.DataFrame, + emaps : pd.DataFrame, + fittype : KrFitFunction, + nbins_dv : int, + zrange_dv : Tuple[float, float], + detector : str, + norm_strat : NormStrategy, + norm_value : float)->DataFrame: + + + ''' + Computes some average parameters (e0, lt, drift v, + S1w, S1h, S1e, S2w, S2h, S2e, S2q, Nsipm, 'Xrms, Yrms) + for a given krypton distribution and for different time slices. + Returns a DataFrame. + + Parameters + ---------- + ts: np.array + Sequence of central times for the different time slices. + masks_time: list of boolean lists + Allows dividing the distribution into time slices. + data: DataFrame + Kdst distribution to analyze. + emaps: correction map + Allows geometrical correction of the energy. + fittype: KrFitFunction + Kind of fit to perform. + nbins_dv: int + Number of bins in Z-coordinate for doing the histogram to compute + the drift velocity. + zrange_dv: int + Range in Z-coordinate for doing the histogram to compute the drift + velocity. + detector: string + Used to get the cathode position from DB for the drift velocity + computation. + norm_strat: NormStrategy + Normalization strategy to follow. + norm_value: float + Energy value to normalize to. + + Returns + ------- + evol_pars: pd.DataFrame + Dataframe containing the parameters evolution. Each column corresponds + to the average value for a given parameter. Each row corresponds to the + parameters for a given time slice. + ''' + + frames = [] + + for index in range(len(masks_time)): + + sel_dst = dst[masks_time[index]] + pars = computing_kr_parameters(data = sel_dst, + ts = ts[index], + emaps = emaps, + fittype = fittype, + nbins_dv = nbins_dv, + zrange_dv = zrange_dv, + detector = detector, + norm_strat = norm_strat, + norm_value = norm_value) + frames.append(pars) + + evol_pars = pd.concat(frames, ignore_index=True) + + return evol_pars From 6df626539fd576cf1d492a6d082e83f2d59f838e Mon Sep 17 00:00:00 2001 From: carhc Date: Fri, 27 Sep 2024 03:24:45 +0200 Subject: [PATCH 11/48] Add cut_effs_evolution function --- invisible_cities/reco/krmap_evolution.py | 61 +++++++++++++++++++++++- 1 file changed, 60 insertions(+), 1 deletion(-) diff --git a/invisible_cities/reco/krmap_evolution.py b/invisible_cities/reco/krmap_evolution.py index 6eb2993da8..618e01cd9b 100644 --- a/invisible_cities/reco/krmap_evolution.py +++ b/invisible_cities/reco/krmap_evolution.py @@ -14,7 +14,8 @@ from .. reco.corrections import apply_all_correction from .. core.stat_functions import poisson_sigma from .. database import load_db as DB -from .. reco.icaro_components import get_fit_function_lt +from .. reco.icaro_components import get_fit_function_lt # Won't work until previous PR are approved + def sigmoid(x : np.array, scale : float, @@ -477,3 +478,61 @@ def kr_time_evolution(ts : np.array[float], evol_pars = pd.concat(frames, ignore_index=True) return evol_pars + + +def cut_effs_evolution(masks_time : List[np.array], + dst : pd.DataFrame, + mask_s1 : np.array, + mask_s2 : np.array, + mask_band : np.array, + evol_table : pd.DataFrame): + + ''' + Computes the efficiencies in time evolution for different time slices. + Returns the input DataFrame updated with S1eff, S2eff, Bandeff. + + Parameters + ---------- + masks_time: list of time masks + Masks which divide the data into time slices. + data: pd.DataFrame + kdst data. + mask_s1: np.array + Mask of S1 cut. + mask_s2: np.array + Mask of S2 cut. + mask_band: np.array + Mask of band cut. + evol_table: pd.DataFrame + Table of Kr evolution parameters. + + Returns + ------- + evol_table_updated: pd.DataFrame + Kr evolution parameters table updated with efficiencies. + ''' + + len_ts = len(masks_time) + + n0 = np.zeros(len_ts) + nS1 = np.zeros(len_ts) + nS2 = np.zeros(len_ts) + nBand = np.zeros(len_ts) + + for index in range(len_ts): + + time_mask = masks_time[index] + nS1mask = time_mask & mask_s1 + nS2mask = nS1mask & mask_s2 + nBandmask = nS2mask & mask_band + + n0 [index] = dst[time_mask].event.nunique() + nS1 [index] = dst[nS1mask] .event.nunique() + nS2 [index] = dst[nS2mask] .event.nunique() + nBand[index] = dst[nBandmask].event.nunique() + + evol_table_updated = evol_table.assign(S1eff = nS1 / n0, + S2eff = nS2 / nS1, + Bandeff = nBand / nS2) + + return evol_table_updated From 075d7c84da7d7a7e115d8225bc2b31c0800e23af Mon Sep 17 00:00:00 2001 From: carhc Date: Fri, 27 Sep 2024 03:35:08 +0200 Subject: [PATCH 12/48] Add add_krevol function --- invisible_cities/reco/krmap_evolution.py | 62 ++++++++++++++++++++++++ 1 file changed, 62 insertions(+) diff --git a/invisible_cities/reco/krmap_evolution.py b/invisible_cities/reco/krmap_evolution.py index 618e01cd9b..49a31acb78 100644 --- a/invisible_cities/reco/krmap_evolution.py +++ b/invisible_cities/reco/krmap_evolution.py @@ -536,3 +536,65 @@ def cut_effs_evolution(masks_time : List[np.array], Bandeff = nBand / nS2) return evol_table_updated + + +def add_krevol(r_fid : float, # Esto sería para meter en la ciudad de ICARO, flow.map(funcion, args, out) etc + nStimeprofile : int, + **map_params): # PREGUNTA: Pongo aquí explícitamente todos los argumentos? O se pueden meter con un diccionario? + + ''' + Adds the time evolution dataframe to the map. + + Parameters + --------- + r_fid: float + Maximum radius for fiducial sample. + nStimeprofile: int + Number of seconds for each time bin. + map_params: dict + Dictionary containing the config file variables. + + Returns + --------- + Function which takes as input map, kr_data, and kr_mask + and returns the time evolution. + ''' + + def add_krevol(map, kdst, mask_s1, # Más de lo mismo con la pregunta anterior + mask_s2, mask_band, fittype, + nbins_dv, zrange_dv, detector): + + fid_sel = (kdst.R < r_fid) & mask_s1 & mask_s2 & mask_band + dstf = kdst[fid_sel] + min_time = dstf.time.min() + max_time = dstf.time.max() + + ntimebins = get_number_of_time_bins(nStimeprofile = nStimeprofile, + tstart = min_time, + tfinal = max_time) + + ts, masks_time = get_time_series_df(ntimebins = ntimebins, + time_range = (min_time, max_time), + dst = kdst) + + masks_timef = [mask[fid_sel] for mask in masks_time] + + evol_table = kr_time_evolution(ts = ts, + masks_time = masks_timef, + dst = dstf, + emaps = map, + fittype = fittype, + nbins_dv = nbins_dv, + zrange_dv = zrange_dv, + detector = detector) + + evol_table_eff = cut_effs_evolution(masks_time = masks_time, + data = kdst, + mask_s1 = mask_s1, + mask_s2 = mask_s2, + mask_band = mask_band, + evol_table = evol_table) + + return evol_table_eff + + return add_krevol From b9bd5cf669d3103fc5657c334e3892e011418f3b Mon Sep 17 00:00:00 2001 From: carhc Date: Tue, 1 Oct 2024 12:21:15 +0200 Subject: [PATCH 13/48] Update comments for clarity --- invisible_cities/reco/krmap_evolution.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/invisible_cities/reco/krmap_evolution.py b/invisible_cities/reco/krmap_evolution.py index 49a31acb78..3bbbbdc4c8 100644 --- a/invisible_cities/reco/krmap_evolution.py +++ b/invisible_cities/reco/krmap_evolution.py @@ -72,6 +72,11 @@ def gauss_seed(x : np.array, Tuple with the seed estimation. ''' + # Looks for the higher Y value and takes its corresponding X value + # as the center of the Gaussian. Based on the "sigma_rel" parameter, + # applies some error to that X value and provides the estimation for the + # amplitude, the center and the sigma. + y_max = np.argmax(y) x_max = x[y_max] sigma = sigma_rel * x_max @@ -85,7 +90,7 @@ def resolution(values : np.array, errors : np.array): ''' - Computes the resolution (FWHM) from the Gaussian parameters. + Computes the resolution (% FWHM) from the Gaussian parameters. Parameters ---------- @@ -174,7 +179,7 @@ def get_time_series_df(ntimebins : int, ''' Given a dst this function returns a time series (ts) and a list of masks which are used to divide - the event in time intervals. + the dst in time intervals. Parameters ---------- From 862ac4f9f584575b30290319cadf68d2e1075b7c Mon Sep 17 00:00:00 2001 From: carhc Date: Wed, 2 Oct 2024 11:01:56 +0200 Subject: [PATCH 14/48] Remove norm_value and norm_strat dependency --- invisible_cities/reco/krmap_evolution.py | 15 +++------------ 1 file changed, 3 insertions(+), 12 deletions(-) diff --git a/invisible_cities/reco/krmap_evolution.py b/invisible_cities/reco/krmap_evolution.py index 3bbbbdc4c8..a253d15948 100644 --- a/invisible_cities/reco/krmap_evolution.py +++ b/invisible_cities/reco/krmap_evolution.py @@ -293,9 +293,7 @@ def computing_kr_parameters(data : DataFrame, fittype : KrFitFunction, nbins_dv : int, zrange_dv : List[float, float], - detector : str, - norm_strat : NormStrategy, - norm_value : float)->DataFrame: # REVISAR NORM_STRAT Y NORM_VALUE + detector : str)->DataFrame: ''' Computes some average parameters (e0, lt, drift v, energy @@ -321,10 +319,6 @@ def computing_kr_parameters(data : DataFrame, detector: string Used to get the cathode position from DB for the drift velocity computation. - norm_strat: NormStrategy - Normalization strategy to follow. - norm_value: float - Energy value to normalize to. Returns ------- @@ -335,7 +329,7 @@ def computing_kr_parameters(data : DataFrame, # Computing E0, LT geo_correction_factor = e0_xy_correction(map = emaps, - norm_strat = norm_strat) # PREGUNTAR POR ESTRATEGIA + norm_strat = NormStrategy.max) fit_func, seed = get_fit_function_lt(fittype) @@ -419,10 +413,7 @@ def kr_time_evolution(ts : np.array[float], fittype : KrFitFunction, nbins_dv : int, zrange_dv : Tuple[float, float], - detector : str, - norm_strat : NormStrategy, - norm_value : float)->DataFrame: - + detector : str) -> pd.DataFrame: ''' Computes some average parameters (e0, lt, drift v, From 3d3dbe1382aa6166d0ec614714af2cb1a0fa264f Mon Sep 17 00:00:00 2001 From: carhc Date: Wed, 2 Oct 2024 11:03:38 +0200 Subject: [PATCH 15/48] Generalize computing_kr_parameters for different fittype --- invisible_cities/reco/krmap_evolution.py | 33 ++++++++++++++++++------ 1 file changed, 25 insertions(+), 8 deletions(-) diff --git a/invisible_cities/reco/krmap_evolution.py b/invisible_cities/reco/krmap_evolution.py index a253d15948..320b82784f 100644 --- a/invisible_cities/reco/krmap_evolution.py +++ b/invisible_cities/reco/krmap_evolution.py @@ -14,7 +14,7 @@ from .. reco.corrections import apply_all_correction from .. core.stat_functions import poisson_sigma from .. database import load_db as DB -from .. reco.icaro_components import get_fit_function_lt # Won't work until previous PR are approved +from .. reco.icaro_components import get_fit_function_lt, select_fit_variables, transform_parameters # Won't work until previous PR are approved def sigmoid(x : np.array, @@ -328,22 +328,39 @@ def computing_kr_parameters(data : DataFrame, # Computing E0, LT + fit_func, seed = get_fit_function_lt(fittype) + geo_correction_factor = e0_xy_correction(map = emaps, norm_strat = NormStrategy.max) - fit_func, seed = get_fit_function_lt(fittype) - x = data.DT, - y = data.S2e.to_numpy()*geo_correction_factor(data.X.to_numpy(), data.Y.to_numpy()) + x, y = select_fit_variables(fittype, data) # Importar de previo PR + + corr = geo_correction_factor(data.X.to_numpy(), data.Y.to_numpy()) + + y_corr = y -np.log(corr) if fittype == KrFitFunction.log_lin else y*corr + + # If we didn't care about using the specific fit function of the map computation, + # this would be reduced to just the following: + # y_corr = y*geo_correction_factor(data.X.to_numpy(), data.Y.to_numpy()) fit_output, _, _, _ = fit(func = fit_func, # Misma funcion que en el ajuste del mapa x = x, - y = y, - seed = seed(x, y), + y = y_corr, + seed = seed(x, y_corr), full_output = False) - e0, lt = fit_output.values - e0u, ltu = fit_output.errors + if fittype == KrFitFunction.log_lin: + + par, err, cov = transform_parameters(fit_output) + + e0, lt = par + e0u, ltu = err + + else: + + e0, lt = fit_output.values + e0u, ltu = fit_output.errors # Computing Drift Velocity From 73e95e022af29130dd86a0615289f015e4a938e7 Mon Sep 17 00:00:00 2001 From: carhc Date: Wed, 2 Oct 2024 11:05:59 +0200 Subject: [PATCH 16/48] Change from Z to DT in compute_drift_v --- invisible_cities/reco/krmap_evolution.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/invisible_cities/reco/krmap_evolution.py b/invisible_cities/reco/krmap_evolution.py index 320b82784f..e4dbb69099 100644 --- a/invisible_cities/reco/krmap_evolution.py +++ b/invisible_cities/reco/krmap_evolution.py @@ -364,9 +364,10 @@ def computing_kr_parameters(data : DataFrame, # Computing Drift Velocity - dv, dvu = compute_drift_v(zdata = data.Z.to_numpy(), + dv, dvu = compute_drift_v(zdata = data.DT.to_numpy(), nbins = nbins_dv, zrange = zrange_dv, + seed = None, detector = detector) # Computing Resolution From be06c3f6e8527d9506e48cf5af0395c90fe8816d Mon Sep 17 00:00:00 2001 From: carhc Date: Wed, 2 Oct 2024 11:09:04 +0200 Subject: [PATCH 17/48] Propagate the removal of norm_strat and norm_value dependency --- invisible_cities/reco/krmap_evolution.py | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/invisible_cities/reco/krmap_evolution.py b/invisible_cities/reco/krmap_evolution.py index e4dbb69099..9e6a18168a 100644 --- a/invisible_cities/reco/krmap_evolution.py +++ b/invisible_cities/reco/krmap_evolution.py @@ -373,9 +373,7 @@ def computing_kr_parameters(data : DataFrame, # Computing Resolution tot_corr_factor = apply_all_correction(maps = emaps, - apply_temp = False, - norm_strat = norm_strat, - norm_value = norm_value) + apply_temp = False) nbins = int((len(data.S2e))**0.5) # Binning as a function of nevts. Should we change it? @@ -421,6 +419,14 @@ def computing_kr_parameters(data : DataFrame, 'Xrms' : [mean_d['Xrms']] , 'Xrmsu' : [var_d['Xrms']] , 'Yrms' : [mean_d['Yrms']] , 'Yrmsu' : [var_d['Yrms']]}) + # PARA METER CUT_EFFS_EVOLUTION: Esta función sólamente se aplica a la dst filtrada, y por lo tanto + # no hay una colección de máscaras. Aquí sólo se mete el timestamp que hace referencia al tiempo del + # rango temporal escogido para el tratamiento de la evolución temporal. Esto va dentro del bucle de + # kr_time_evolution. Si queremos meter la evolución de eficiencias de cortes, yo la metería en esa + # función y no en esta (aunque aquí en un primer momento pueda parecer más lógico) ya que en la otra + # es posible filtrar la dst entera y aquí no según el flow de la función. + + return evol From 58b891a6fb6461cab6d0bf8c27327d99d6f15631 Mon Sep 17 00:00:00 2001 From: carhc Date: Wed, 2 Oct 2024 11:10:48 +0200 Subject: [PATCH 18/48] Modify kr_time_evolution to avoid explicit loops and appends --- invisible_cities/reco/krmap_evolution.py | 58 +++++++++++++++--------- 1 file changed, 36 insertions(+), 22 deletions(-) diff --git a/invisible_cities/reco/krmap_evolution.py b/invisible_cities/reco/krmap_evolution.py index 9e6a18168a..ea51afa4be 100644 --- a/invisible_cities/reco/krmap_evolution.py +++ b/invisible_cities/reco/krmap_evolution.py @@ -441,7 +441,7 @@ def kr_time_evolution(ts : np.array[float], ''' Computes some average parameters (e0, lt, drift v, - S1w, S1h, S1e, S2w, S2h, S2e, S2q, Nsipm, 'Xrms, Yrms) + S1w, S1h, S1e, S2w, S2h, S2e, S2q, Nsipm, Xrms, Yrms) for a given krypton distribution and for different time slices. Returns a DataFrame. @@ -451,7 +451,7 @@ def kr_time_evolution(ts : np.array[float], Sequence of central times for the different time slices. masks_time: list of boolean lists Allows dividing the distribution into time slices. - data: DataFrame + dst: DataFrame Kdst distribution to analyze. emaps: correction map Allows geometrical correction of the energy. @@ -460,42 +460,56 @@ def kr_time_evolution(ts : np.array[float], nbins_dv: int Number of bins in Z-coordinate for doing the histogram to compute the drift velocity. - zrange_dv: int + zrange_dv: Tuple Range in Z-coordinate for doing the histogram to compute the drift velocity. detector: string Used to get the cathode position from DB for the drift velocity computation. - norm_strat: NormStrategy - Normalization strategy to follow. - norm_value: float - Energy value to normalize to. Returns ------- evol_pars: pd.DataFrame - Dataframe containing the parameters evolution. Each column corresponds + DataFrame containing the parameters evolution. Each column corresponds to the average value for a given parameter. Each row corresponds to the parameters for a given time slice. ''' - frames = [] + # Create a DataFrame with the time slices (ts) + mask_time_df = pd.DataFrame({'ts': ts}) + + def compute_parameters(row, masks_time, dst, emaps, fittype, nbins_dv, zrange_dv, detector): + idx = row.name + mask = masks_time[idx] + time = row['ts'] + + # Select the data for the given mask + sel_dst = dst[mask] + + # Compute the krypton parameters using the selected data + pars = computing_kr_parameters(data = sel_dst, + ts = time, + emaps = emaps, + fittype = fittype, + nbins_dv = nbins_dv, + zrange_dv = zrange_dv, + detector = detector) + + + + if pars is None or pars.empty: + return pd.Series(np.nan) # Return nan Series if no parameters found + + if isinstance(pars, pd.DataFrame) and len(pars) == 1: + pars = pars.iloc[0] # Convert single-row DataFrame to Series - for index in range(len(masks_time)): + return pd.Series(pars) - sel_dst = dst[masks_time[index]] - pars = computing_kr_parameters(data = sel_dst, - ts = ts[index], - emaps = emaps, - fittype = fittype, - nbins_dv = nbins_dv, - zrange_dv = zrange_dv, - detector = detector, - norm_strat = norm_strat, - norm_value = norm_value) - frames.append(pars) - evol_pars = pd.concat(frames, ignore_index=True) + # Apply the function to the previous df and compute krypton parameters for each time slice + evol_pars = mask_time_df.apply(lambda row: compute_parameters(row, masks_time, dst, emaps, + fittype, nbins_dv, zrange_dv, detector), + axis=1) return evol_pars From 280c79213f82b9859394795f9883b3ce4502a378 Mon Sep 17 00:00:00 2001 From: carhc Date: Wed, 2 Oct 2024 15:42:10 +0200 Subject: [PATCH 19/48] New ideas in comments (delete this) --- invisible_cities/reco/krmap_evolution.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/invisible_cities/reco/krmap_evolution.py b/invisible_cities/reco/krmap_evolution.py index ea51afa4be..7a81239896 100644 --- a/invisible_cities/reco/krmap_evolution.py +++ b/invisible_cities/reco/krmap_evolution.py @@ -326,6 +326,11 @@ def computing_kr_parameters(data : DataFrame, Each column corresponds to the average value of a different parameter. ''' + #effs1, effs2, effbans = efficiencies(data, masks) + + #data = data[np.logical_and(masks)] + + # Computing E0, LT fit_func, seed = get_fit_function_lt(fittype) @@ -484,10 +489,20 @@ def compute_parameters(row, masks_time, dst, emaps, fittype, nbins_dv, zrange_dv time = row['ts'] # Select the data for the given mask + # nS1mask = time_mask & mask_s1 + # nS2mask = nS1mask & mask_s2 + # nBandmask = nS2mask & mask_band + # ntot = dst[mask].nunique() + # ns1 = dst[mask & masks_s1].nunique() + # ns2 = dst[mask & mask_s1 & mask_s2]. unique() + # nband = dst[mask & mask_s1 & mask_s2 & mask_band].unique() + + sel_dst = dst[mask] # Compute the krypton parameters using the selected data pars = computing_kr_parameters(data = sel_dst, + # masks = masks, ts = time, emaps = emaps, fittype = fittype, From 2ce357108f6f394e860ccfe6afc4f0d9f4da2df0 Mon Sep 17 00:00:00 2001 From: carhc Date: Wed, 9 Oct 2024 10:10:45 +0200 Subject: [PATCH 20/48] Delete this --- invisible_cities/reco/krmap_evolution.py | 290 ++++++++++++++++++----- 1 file changed, 232 insertions(+), 58 deletions(-) diff --git a/invisible_cities/reco/krmap_evolution.py b/invisible_cities/reco/krmap_evolution.py index 7a81239896..83f8bd018b 100644 --- a/invisible_cities/reco/krmap_evolution.py +++ b/invisible_cities/reco/krmap_evolution.py @@ -5,8 +5,9 @@ from pandas import DataFrame from .. types.symbols import NormStrategy -from .. types.symbols import KrFitFunction # Won't work until previous PR are approved -from .. core.fit_functions import fit, gauss +from .. evm.ic_containers import FitFunction +# from .. types.symbols import KrFitFunction # Won't work until previous PR are approved +from .. core.fit_functions import fit, gauss, polynom, expo from .. core.core_functions import in_range, shift_to_bin_centers from .. reco.corrections import get_normalization_factor from .. reco.corrections import correct_geometry_ @@ -14,7 +15,186 @@ from .. reco.corrections import apply_all_correction from .. core.stat_functions import poisson_sigma from .. database import load_db as DB -from .. reco.icaro_components import get_fit_function_lt, select_fit_variables, transform_parameters # Won't work until previous PR are approved +# from .. reco.icaro_components import get_fit_function_lt, select_fit_variables, transform_parameters # Won't work until previous PR are approved + + +from enum import auto +from invisible_cities.evm.ic_containers import FitFunction +from invisible_cities.types.symbols import AutoNameEnumBase + +class KrFitFunction(AutoNameEnumBase): + expo = auto() + linear = auto() + log_lin = auto() + +def lin_seed(x : np.array, + y : np.array): + + ''' + Estimate the seed for a linear fit. + + Parameters + ---------- + x : np.array + Independent variable. + y : np.array + Dependent variable. + + Returns + ------- + seed : tuple + Seed parameters (intercept, slope) for the linear fit. + ''' + + x0, x1 = x.min(), x.max() + y0, y1 = y.min(), y.max() + + if x1 == x0: # If same x value, set slope to 0 and use the mean value of y as interceipt + b = 0 + a = y.mean() + + else: + b = (y1 - y0) / (x1 - x0) + a = y0 - b * x0 + + seed = a, b + + return seed + + +def expo_seed(x : np.array, + y : np.array): + + ''' + Estimate the seed for an exponential fit. + + Parameters + ---------- + x : np.array + Independent variable. + y : np.array + Dependent variable. + + Returns + ------- + seed : tuple + Seed parameters (constant, mean) for the exponential fit. + ''' + + x, y = zip(*sorted(zip(x, y))) + + const = y[0] + + if const <= 0 or y[-1] <= 0: + raise ValueError("y data must be > 0") + + lt = -x[-1] / np.log(y[-1] / const) + + seed = const, lt + return seed + + +def select_fit_variables(fittype : KrFitFunction, + dst : pd.DataFrame): + + ''' + Select the data for fitting based on the specified fit type. + + NOTES: Since x axis (DT) is never altered, maybe we can just + return the y values. However, when we implement the binned fit, + the profile could be done here (?) so it would make sense to + always provide both x and y. We could rename parameters and have + fittype (binned / unbinned) and fitfunction (lin, expo, log-lin...) + + Parameters + ---------- + fittype : KrFitFunction + The type of fit function to prepare data for (e.g., linear, exponential, log-linear). + dst : pd.DataFrame + The DataFrame containing the data to be prepared for fitting. + + Returns + ------- + x_data : pd.Series + The independent variable data prepared for fitting. + y_data : pd.Series + The dependent variable data prepared for fitting. + ''' + + if fittype is KrFitFunction.linear : return dst.DT, dst.S2e + elif fittype is KrFitFunction.expo : return dst.DT, dst.S2e + elif fittype is KrFitFunction.log_lin: return dst.DT, -np.log(dst.S2e) + + +def get_fit_function_lt(fittype : KrFitFunction): + + ''' + Retrieve the fitting function and seed function based on the + specified fittype. + + Parameters + ---------- + fittype : KrFitFunction + The type of fit function to retrieve (e.g., linear, exponential, log-linear). + + Returns + ------- + fit_function : function + The fitting function corresponding to the specified fit type. + seed_function : function + The seed function corresponding to the specified fit type. + ''' + + linear_function = lambda x, y0, slope: polynom(x, y0, slope) + expo_function = lambda x, e0, lt: expo (x, e0, -lt) + + if fittype is KrFitFunction.linear: return linear_function, lin_seed + elif fittype is KrFitFunction.log_lin: return linear_function, lin_seed + elif fittype is KrFitFunction.expo: return expo_function, expo_seed + + +def transform_parameters(fit_output : FitFunction): + + ''' + Transform the parameters obtained from the fitting output into EO and LT. + When using log_lin fit, we need to convert the intermediate variables into + the actual physical magnitudes involved in the process. + + Parameters + ---------- + fit_output : FitFunction + Output from IC's fit containing the parameter values, errors, and + covariance matrix. + + Returns + ------- + par : list + Transformed parameter values. + err : list + Transformed parameter errors. + cov : float + Transformed covariance value. + ''' + + par = fit_output.values + err = fit_output.errors + cov = fit_output.cov[0, 1] + + a, b = par + u_a, u_b = err + + E0 = np.exp(-a) + s_E0 = np.abs(E0 * u_a) + + lt = 1 / b + s_lt = np.abs(lt**2 * u_b) + + cov = E0 * lt**2 * cov # Not sure about this + + par = [ E0, lt] + err = [s_E0, s_lt] + + return par, err, cov def sigmoid(x : np.array, @@ -292,12 +472,12 @@ def computing_kr_parameters(data : DataFrame, emaps : pd.DataFrame, fittype : KrFitFunction, nbins_dv : int, - zrange_dv : List[float, float], + dt_range_dv : Tuple[float, float], detector : str)->DataFrame: ''' Computes some average parameters (e0, lt, drift v, energy - resolution, S1w, S1h, S1e, S2w, S2h, S2e, S2q, Nsipm, 'Xrms, Yrms) + resolution, S1w, S1h, S1e, S2w, S2h, S2e, S2q, Nsipm, Xrms, Yrms) for a given krypton distribution. Parameters @@ -313,7 +493,7 @@ def computing_kr_parameters(data : DataFrame, nbins_dv: int Number of bins in Z-coordinate to use in the histogram for the drift velocity calculation. - zrange_dv: List + dt_range_dv: List Range in Z-coordinate to use in the histogram for the drift velocity calculation. detector: string @@ -326,34 +506,34 @@ def computing_kr_parameters(data : DataFrame, Each column corresponds to the average value of a different parameter. ''' - #effs1, effs2, effbans = efficiencies(data, masks) - - #data = data[np.logical_and(masks)] - - # Computing E0, LT - fit_func, seed = get_fit_function_lt(fittype) - + fit_func, seed = get_fit_function_lt(fittype = fittype) + print('fit_func and seed ok') geo_correction_factor = e0_xy_correction(map = emaps, norm_strat = NormStrategy.max) - + print('geo corr function ok') x, y = select_fit_variables(fittype, data) # Importar de previo PR corr = geo_correction_factor(data.X.to_numpy(), data.Y.to_numpy()) - - y_corr = y -np.log(corr) if fittype == KrFitFunction.log_lin else y*corr - + print('calculating corrections based on map ok') + y_corr = y - np.log(corr) if fittype == KrFitFunction.log_lin else y*corr + print('applying corrections ok') # If we didn't care about using the specific fit function of the map computation, # this would be reduced to just the following: # y_corr = y*geo_correction_factor(data.X.to_numpy(), data.Y.to_numpy()) - + print(x) + print(y) + print(data.X) + print(data.Y) + print(corr) + print(y_corr) fit_output, _, _, _ = fit(func = fit_func, # Misma funcion que en el ajuste del mapa - x = x, - y = y_corr, + x = x.to_numpy(), + y = y_corr.to_numpy(), seed = seed(x, y_corr), - full_output = False) + full_output = True) if fittype == KrFitFunction.log_lin: @@ -371,7 +551,7 @@ def computing_kr_parameters(data : DataFrame, dv, dvu = compute_drift_v(zdata = data.DT.to_numpy(), nbins = nbins_dv, - zrange = zrange_dv, + zrange = dt_range_dv, seed = None, detector = detector) @@ -392,6 +572,17 @@ def computing_kr_parameters(data : DataFrame, res, err_res = resolution(values = par, errors = err) + # Computing Efficiencies + + # n0 = data.event.nunique() + # nS1 = data.query('nS1 == 1').event.nunique() + # nS2 = data.query('nS2 == 1').event.nunique() + # nband = nS2 # Best way to implement band? + + # S1eff = nS1 / n0, + # S2eff = nS2 / nS1 + # B_eff = nband / nS2 + # Averages of parameters parameters = ['S1w', 'S1h', 'S1e', @@ -408,7 +599,8 @@ def computing_kr_parameters(data : DataFrame, # Creating parameter evolution table - evol = DataFrame({'ts' : [ts] , + evol = DataFrame({'ts' : [ts] , # 'S1eff' : [S1eff] , + # 'S2eff': [S2eff] , 'B_eff' : [B_eff] , 'e0' : [e0] , 'e0u' : [e0u] , 'lt' : [lt] , 'ltu' : [ltu] , 'dv' : [dv] , 'dvu' : [dvu] , @@ -435,13 +627,13 @@ def computing_kr_parameters(data : DataFrame, return evol -def kr_time_evolution(ts : np.array[float], +def kr_time_evolution(ts : np.array, masks_time : List[np.array], dst : pd.DataFrame, emaps : pd.DataFrame, fittype : KrFitFunction, nbins_dv : int, - zrange_dv : Tuple[float, float], + dt_range_dv : Tuple[float, float], detector : str) -> pd.DataFrame: ''' @@ -465,7 +657,7 @@ def kr_time_evolution(ts : np.array[float], nbins_dv: int Number of bins in Z-coordinate for doing the histogram to compute the drift velocity. - zrange_dv: Tuple + dt_range_dv: Tuple Range in Z-coordinate for doing the histogram to compute the drift velocity. detector: string @@ -483,33 +675,21 @@ def kr_time_evolution(ts : np.array[float], # Create a DataFrame with the time slices (ts) mask_time_df = pd.DataFrame({'ts': ts}) - def compute_parameters(row, masks_time, dst, emaps, fittype, nbins_dv, zrange_dv, detector): + def compute_parameters(row, masks_time, dst, emaps, fittype, nbins_dv, dt_range_dv, detector): idx = row.name mask = masks_time[idx] time = row['ts'] - - # Select the data for the given mask - # nS1mask = time_mask & mask_s1 - # nS2mask = nS1mask & mask_s2 - # nBandmask = nS2mask & mask_band - # ntot = dst[mask].nunique() - # ns1 = dst[mask & masks_s1].nunique() - # ns2 = dst[mask & mask_s1 & mask_s2]. unique() - # nband = dst[mask & mask_s1 & mask_s2 & mask_band].unique() - - + print(idx) sel_dst = dst[mask] - # Compute the krypton parameters using the selected data - pars = computing_kr_parameters(data = sel_dst, - # masks = masks, - ts = time, - emaps = emaps, - fittype = fittype, - nbins_dv = nbins_dv, - zrange_dv = zrange_dv, - detector = detector) - + # Compute the krypton parameters using the time-binned data + pars = computing_kr_parameters(data = sel_dst, + ts = time, + emaps = emaps, + fittype = fittype, + nbins_dv = nbins_dv, + dt_range_dv = dt_range_dv, + detector = detector) if pars is None or pars.empty: @@ -523,7 +703,7 @@ def compute_parameters(row, masks_time, dst, emaps, fittype, nbins_dv, zrange_dv # Apply the function to the previous df and compute krypton parameters for each time slice evol_pars = mask_time_df.apply(lambda row: compute_parameters(row, masks_time, dst, emaps, - fittype, nbins_dv, zrange_dv, detector), + fittype, nbins_dv, dt_range_dv, detector), axis=1) return evol_pars @@ -611,7 +791,7 @@ def add_krevol(r_fid : float, # Esto sería para meter en la ciudad de def add_krevol(map, kdst, mask_s1, # Más de lo mismo con la pregunta anterior mask_s2, mask_band, fittype, - nbins_dv, zrange_dv, detector): + nbins_dv, dt_range_dv, detector): fid_sel = (kdst.R < r_fid) & mask_s1 & mask_s2 & mask_band dstf = kdst[fid_sel] @@ -634,16 +814,10 @@ def add_krevol(map, kdst, mask_s1, # Más de lo mismo con la pregunta emaps = map, fittype = fittype, nbins_dv = nbins_dv, - zrange_dv = zrange_dv, + dt_range_dv = dt_range_dv, detector = detector) - evol_table_eff = cut_effs_evolution(masks_time = masks_time, - data = kdst, - mask_s1 = mask_s1, - mask_s2 = mask_s2, - mask_band = mask_band, - evol_table = evol_table) - return evol_table_eff + return evol_table return add_krevol From c4af750fe5919af55842bc142d933a840d4f1a5b Mon Sep 17 00:00:00 2001 From: carhc Date: Fri, 11 Oct 2024 11:56:53 +0200 Subject: [PATCH 21/48] Revert "New ideas in comments (delete this)" This reverts commit ddefe930c9d4e3d8994c83a5cee4a565ccf19e09. --- invisible_cities/reco/krmap_evolution.py | 275 ++++------------------- 1 file changed, 43 insertions(+), 232 deletions(-) diff --git a/invisible_cities/reco/krmap_evolution.py b/invisible_cities/reco/krmap_evolution.py index 83f8bd018b..ea51afa4be 100644 --- a/invisible_cities/reco/krmap_evolution.py +++ b/invisible_cities/reco/krmap_evolution.py @@ -5,9 +5,8 @@ from pandas import DataFrame from .. types.symbols import NormStrategy -from .. evm.ic_containers import FitFunction -# from .. types.symbols import KrFitFunction # Won't work until previous PR are approved -from .. core.fit_functions import fit, gauss, polynom, expo +from .. types.symbols import KrFitFunction # Won't work until previous PR are approved +from .. core.fit_functions import fit, gauss from .. core.core_functions import in_range, shift_to_bin_centers from .. reco.corrections import get_normalization_factor from .. reco.corrections import correct_geometry_ @@ -15,186 +14,7 @@ from .. reco.corrections import apply_all_correction from .. core.stat_functions import poisson_sigma from .. database import load_db as DB -# from .. reco.icaro_components import get_fit_function_lt, select_fit_variables, transform_parameters # Won't work until previous PR are approved - - -from enum import auto -from invisible_cities.evm.ic_containers import FitFunction -from invisible_cities.types.symbols import AutoNameEnumBase - -class KrFitFunction(AutoNameEnumBase): - expo = auto() - linear = auto() - log_lin = auto() - -def lin_seed(x : np.array, - y : np.array): - - ''' - Estimate the seed for a linear fit. - - Parameters - ---------- - x : np.array - Independent variable. - y : np.array - Dependent variable. - - Returns - ------- - seed : tuple - Seed parameters (intercept, slope) for the linear fit. - ''' - - x0, x1 = x.min(), x.max() - y0, y1 = y.min(), y.max() - - if x1 == x0: # If same x value, set slope to 0 and use the mean value of y as interceipt - b = 0 - a = y.mean() - - else: - b = (y1 - y0) / (x1 - x0) - a = y0 - b * x0 - - seed = a, b - - return seed - - -def expo_seed(x : np.array, - y : np.array): - - ''' - Estimate the seed for an exponential fit. - - Parameters - ---------- - x : np.array - Independent variable. - y : np.array - Dependent variable. - - Returns - ------- - seed : tuple - Seed parameters (constant, mean) for the exponential fit. - ''' - - x, y = zip(*sorted(zip(x, y))) - - const = y[0] - - if const <= 0 or y[-1] <= 0: - raise ValueError("y data must be > 0") - - lt = -x[-1] / np.log(y[-1] / const) - - seed = const, lt - return seed - - -def select_fit_variables(fittype : KrFitFunction, - dst : pd.DataFrame): - - ''' - Select the data for fitting based on the specified fit type. - - NOTES: Since x axis (DT) is never altered, maybe we can just - return the y values. However, when we implement the binned fit, - the profile could be done here (?) so it would make sense to - always provide both x and y. We could rename parameters and have - fittype (binned / unbinned) and fitfunction (lin, expo, log-lin...) - - Parameters - ---------- - fittype : KrFitFunction - The type of fit function to prepare data for (e.g., linear, exponential, log-linear). - dst : pd.DataFrame - The DataFrame containing the data to be prepared for fitting. - - Returns - ------- - x_data : pd.Series - The independent variable data prepared for fitting. - y_data : pd.Series - The dependent variable data prepared for fitting. - ''' - - if fittype is KrFitFunction.linear : return dst.DT, dst.S2e - elif fittype is KrFitFunction.expo : return dst.DT, dst.S2e - elif fittype is KrFitFunction.log_lin: return dst.DT, -np.log(dst.S2e) - - -def get_fit_function_lt(fittype : KrFitFunction): - - ''' - Retrieve the fitting function and seed function based on the - specified fittype. - - Parameters - ---------- - fittype : KrFitFunction - The type of fit function to retrieve (e.g., linear, exponential, log-linear). - - Returns - ------- - fit_function : function - The fitting function corresponding to the specified fit type. - seed_function : function - The seed function corresponding to the specified fit type. - ''' - - linear_function = lambda x, y0, slope: polynom(x, y0, slope) - expo_function = lambda x, e0, lt: expo (x, e0, -lt) - - if fittype is KrFitFunction.linear: return linear_function, lin_seed - elif fittype is KrFitFunction.log_lin: return linear_function, lin_seed - elif fittype is KrFitFunction.expo: return expo_function, expo_seed - - -def transform_parameters(fit_output : FitFunction): - - ''' - Transform the parameters obtained from the fitting output into EO and LT. - When using log_lin fit, we need to convert the intermediate variables into - the actual physical magnitudes involved in the process. - - Parameters - ---------- - fit_output : FitFunction - Output from IC's fit containing the parameter values, errors, and - covariance matrix. - - Returns - ------- - par : list - Transformed parameter values. - err : list - Transformed parameter errors. - cov : float - Transformed covariance value. - ''' - - par = fit_output.values - err = fit_output.errors - cov = fit_output.cov[0, 1] - - a, b = par - u_a, u_b = err - - E0 = np.exp(-a) - s_E0 = np.abs(E0 * u_a) - - lt = 1 / b - s_lt = np.abs(lt**2 * u_b) - - cov = E0 * lt**2 * cov # Not sure about this - - par = [ E0, lt] - err = [s_E0, s_lt] - - return par, err, cov +from .. reco.icaro_components import get_fit_function_lt, select_fit_variables, transform_parameters # Won't work until previous PR are approved def sigmoid(x : np.array, @@ -472,12 +292,12 @@ def computing_kr_parameters(data : DataFrame, emaps : pd.DataFrame, fittype : KrFitFunction, nbins_dv : int, - dt_range_dv : Tuple[float, float], + zrange_dv : List[float, float], detector : str)->DataFrame: ''' Computes some average parameters (e0, lt, drift v, energy - resolution, S1w, S1h, S1e, S2w, S2h, S2e, S2q, Nsipm, Xrms, Yrms) + resolution, S1w, S1h, S1e, S2w, S2h, S2e, S2q, Nsipm, 'Xrms, Yrms) for a given krypton distribution. Parameters @@ -493,7 +313,7 @@ def computing_kr_parameters(data : DataFrame, nbins_dv: int Number of bins in Z-coordinate to use in the histogram for the drift velocity calculation. - dt_range_dv: List + zrange_dv: List Range in Z-coordinate to use in the histogram for the drift velocity calculation. detector: string @@ -508,32 +328,27 @@ def computing_kr_parameters(data : DataFrame, # Computing E0, LT - fit_func, seed = get_fit_function_lt(fittype = fittype) - print('fit_func and seed ok') + fit_func, seed = get_fit_function_lt(fittype) + geo_correction_factor = e0_xy_correction(map = emaps, norm_strat = NormStrategy.max) - print('geo corr function ok') + x, y = select_fit_variables(fittype, data) # Importar de previo PR corr = geo_correction_factor(data.X.to_numpy(), data.Y.to_numpy()) - print('calculating corrections based on map ok') - y_corr = y - np.log(corr) if fittype == KrFitFunction.log_lin else y*corr - print('applying corrections ok') + + y_corr = y -np.log(corr) if fittype == KrFitFunction.log_lin else y*corr + # If we didn't care about using the specific fit function of the map computation, # this would be reduced to just the following: # y_corr = y*geo_correction_factor(data.X.to_numpy(), data.Y.to_numpy()) - print(x) - print(y) - print(data.X) - print(data.Y) - print(corr) - print(y_corr) + fit_output, _, _, _ = fit(func = fit_func, # Misma funcion que en el ajuste del mapa - x = x.to_numpy(), - y = y_corr.to_numpy(), + x = x, + y = y_corr, seed = seed(x, y_corr), - full_output = True) + full_output = False) if fittype == KrFitFunction.log_lin: @@ -551,7 +366,7 @@ def computing_kr_parameters(data : DataFrame, dv, dvu = compute_drift_v(zdata = data.DT.to_numpy(), nbins = nbins_dv, - zrange = dt_range_dv, + zrange = zrange_dv, seed = None, detector = detector) @@ -572,17 +387,6 @@ def computing_kr_parameters(data : DataFrame, res, err_res = resolution(values = par, errors = err) - # Computing Efficiencies - - # n0 = data.event.nunique() - # nS1 = data.query('nS1 == 1').event.nunique() - # nS2 = data.query('nS2 == 1').event.nunique() - # nband = nS2 # Best way to implement band? - - # S1eff = nS1 / n0, - # S2eff = nS2 / nS1 - # B_eff = nband / nS2 - # Averages of parameters parameters = ['S1w', 'S1h', 'S1e', @@ -599,8 +403,7 @@ def computing_kr_parameters(data : DataFrame, # Creating parameter evolution table - evol = DataFrame({'ts' : [ts] , # 'S1eff' : [S1eff] , - # 'S2eff': [S2eff] , 'B_eff' : [B_eff] , + evol = DataFrame({'ts' : [ts] , 'e0' : [e0] , 'e0u' : [e0u] , 'lt' : [lt] , 'ltu' : [ltu] , 'dv' : [dv] , 'dvu' : [dvu] , @@ -627,13 +430,13 @@ def computing_kr_parameters(data : DataFrame, return evol -def kr_time_evolution(ts : np.array, +def kr_time_evolution(ts : np.array[float], masks_time : List[np.array], dst : pd.DataFrame, emaps : pd.DataFrame, fittype : KrFitFunction, nbins_dv : int, - dt_range_dv : Tuple[float, float], + zrange_dv : Tuple[float, float], detector : str) -> pd.DataFrame: ''' @@ -657,7 +460,7 @@ def kr_time_evolution(ts : np.array, nbins_dv: int Number of bins in Z-coordinate for doing the histogram to compute the drift velocity. - dt_range_dv: Tuple + zrange_dv: Tuple Range in Z-coordinate for doing the histogram to compute the drift velocity. detector: string @@ -675,21 +478,23 @@ def kr_time_evolution(ts : np.array, # Create a DataFrame with the time slices (ts) mask_time_df = pd.DataFrame({'ts': ts}) - def compute_parameters(row, masks_time, dst, emaps, fittype, nbins_dv, dt_range_dv, detector): + def compute_parameters(row, masks_time, dst, emaps, fittype, nbins_dv, zrange_dv, detector): idx = row.name mask = masks_time[idx] time = row['ts'] - print(idx) + + # Select the data for the given mask sel_dst = dst[mask] - # Compute the krypton parameters using the time-binned data - pars = computing_kr_parameters(data = sel_dst, - ts = time, - emaps = emaps, - fittype = fittype, - nbins_dv = nbins_dv, - dt_range_dv = dt_range_dv, - detector = detector) + # Compute the krypton parameters using the selected data + pars = computing_kr_parameters(data = sel_dst, + ts = time, + emaps = emaps, + fittype = fittype, + nbins_dv = nbins_dv, + zrange_dv = zrange_dv, + detector = detector) + if pars is None or pars.empty: @@ -703,7 +508,7 @@ def compute_parameters(row, masks_time, dst, emaps, fittype, nbins_dv, dt_range_ # Apply the function to the previous df and compute krypton parameters for each time slice evol_pars = mask_time_df.apply(lambda row: compute_parameters(row, masks_time, dst, emaps, - fittype, nbins_dv, dt_range_dv, detector), + fittype, nbins_dv, zrange_dv, detector), axis=1) return evol_pars @@ -791,7 +596,7 @@ def add_krevol(r_fid : float, # Esto sería para meter en la ciudad de def add_krevol(map, kdst, mask_s1, # Más de lo mismo con la pregunta anterior mask_s2, mask_band, fittype, - nbins_dv, dt_range_dv, detector): + nbins_dv, zrange_dv, detector): fid_sel = (kdst.R < r_fid) & mask_s1 & mask_s2 & mask_band dstf = kdst[fid_sel] @@ -814,10 +619,16 @@ def add_krevol(map, kdst, mask_s1, # Más de lo mismo con la pregunta emaps = map, fittype = fittype, nbins_dv = nbins_dv, - dt_range_dv = dt_range_dv, + zrange_dv = zrange_dv, detector = detector) + evol_table_eff = cut_effs_evolution(masks_time = masks_time, + data = kdst, + mask_s1 = mask_s1, + mask_s2 = mask_s2, + mask_band = mask_band, + evol_table = evol_table) - return evol_table + return evol_table_eff return add_krevol From f6b1dedb3527148ad853b203635034c3861e6c8d Mon Sep 17 00:00:00 2001 From: carhc Date: Fri, 11 Oct 2024 11:58:09 +0200 Subject: [PATCH 22/48] Delete get_number_of_time_bins function --- invisible_cities/reco/krmap_evolution.py | 29 ------------------------ 1 file changed, 29 deletions(-) diff --git a/invisible_cities/reco/krmap_evolution.py b/invisible_cities/reco/krmap_evolution.py index ea51afa4be..c295bbf2ca 100644 --- a/invisible_cities/reco/krmap_evolution.py +++ b/invisible_cities/reco/krmap_evolution.py @@ -144,35 +144,6 @@ def quick_gauss_fit(data : np.array, return fit_result -def get_number_of_time_bins(nStimeprofile : int, - tstart : int, - tfinal : int)->int: - - ''' - Computes the number of time bins to use for a given time step - in seconds. - - Parameters - ---------- - nStimeprofile: int - Time step in seconds. - tstart: int - Initial timestamp for the dataset. - tfinal: int - Final timestamp for the dataset. - - Returns - ------- - ntimebins: int - Number of time bins. - ''' - - ntimebins = int(np.floor((tfinal - tstart) / nStimeprofile)) - ntimebins = np.max([ntimebins, 1]) - - return ntimebins - - def get_time_series_df(ntimebins : int, time_range : Tuple[float, float], dst : DataFrame)->Tuple[np.array, List[np.array]]: From 79c78e04509cba79747c4afc2bc905df77711754 Mon Sep 17 00:00:00 2001 From: carhc Date: Fri, 11 Oct 2024 14:12:47 +0200 Subject: [PATCH 23/48] Update from z to DT in some functions --- invisible_cities/reco/krmap_evolution.py | 83 +++++++++++++++++------- 1 file changed, 61 insertions(+), 22 deletions(-) diff --git a/invisible_cities/reco/krmap_evolution.py b/invisible_cities/reco/krmap_evolution.py index c295bbf2ca..6665123422 100644 --- a/invisible_cities/reco/krmap_evolution.py +++ b/invisible_cities/reco/krmap_evolution.py @@ -175,9 +175,9 @@ def get_time_series_df(ntimebins : int, return shift_to_bin_centers(time_bins), masks -def compute_drift_v(zdata : np.array, +def compute_drift_v(dtdata : np.array, nbins : int, - zrange : Tuple[float, float], + dtrange : Tuple[float, float], seed : Tuple[float, float, float, float], detector : str)->Tuple[float, float]: @@ -187,12 +187,12 @@ def compute_drift_v(zdata : np.array, Parameters ---------- - zdata: array_like - Values of Z coordinate. + dtdata: array_like + Values of DT coordinate. nbins: int (optional) The number of bins in the z coordinate for the binned fit. - zrange: length-2 tuple (optional) - Fix the range in z. + dtrange: length-2 tuple (optional) + Fix the range in DT. seed: length-4 tuple (optional) Seed for the fit. detector: string (optional) @@ -206,16 +206,16 @@ def compute_drift_v(zdata : np.array, Drift velocity uncertainty. ''' - y, x = np.histogram(zdata, nbins, zrange) + y, x = np.histogram(dtdata, nbins, dtrange) x = shift_to_bin_centers(x) - if seed is None: seed = np.max(y), np.mean(zrange), 0.5, np.min(y) + if seed is None: seed = np.max(y), np.mean(dtrange), 0.5, np.min(y) # At the moment there is not NEXT-100 DB so this won't work for that geometry z_cathode = DB.DetectorGeo(detector).ZMAX[0] try: - f = fit(sigmoid, x, y, seed, sigma = poisson_sigma(y), fit_range = zrange) + f = fit(sigmoid, x, y, seed, sigma = poisson_sigma(y), fit_range = dtrange) par = f.values err = f.errors @@ -263,7 +263,7 @@ def computing_kr_parameters(data : DataFrame, emaps : pd.DataFrame, fittype : KrFitFunction, nbins_dv : int, - zrange_dv : List[float, float], + dtrange_dv : List[float], detector : str)->DataFrame: ''' @@ -284,7 +284,7 @@ def computing_kr_parameters(data : DataFrame, nbins_dv: int Number of bins in Z-coordinate to use in the histogram for the drift velocity calculation. - zrange_dv: List + dtrange_dv: List Range in Z-coordinate to use in the histogram for the drift velocity calculation. detector: string @@ -335,9 +335,9 @@ def computing_kr_parameters(data : DataFrame, # Computing Drift Velocity - dv, dvu = compute_drift_v(zdata = data.DT.to_numpy(), + dv, dvu = compute_drift_v(dtdata = data.DT.to_numpy(), nbins = nbins_dv, - zrange = zrange_dv, + dtrange = dtrange_dv, seed = None, detector = detector) @@ -407,7 +407,7 @@ def kr_time_evolution(ts : np.array[float], emaps : pd.DataFrame, fittype : KrFitFunction, nbins_dv : int, - zrange_dv : Tuple[float, float], + dtrange_dv : Tuple[float, float], detector : str) -> pd.DataFrame: ''' @@ -431,7 +431,7 @@ def kr_time_evolution(ts : np.array[float], nbins_dv: int Number of bins in Z-coordinate for doing the histogram to compute the drift velocity. - zrange_dv: Tuple + dtrange_dv: Tuple Range in Z-coordinate for doing the histogram to compute the drift velocity. detector: string @@ -449,7 +449,7 @@ def kr_time_evolution(ts : np.array[float], # Create a DataFrame with the time slices (ts) mask_time_df = pd.DataFrame({'ts': ts}) - def compute_parameters(row, masks_time, dst, emaps, fittype, nbins_dv, zrange_dv, detector): + def compute_parameters(row, masks_time, dst, emaps, fittype, nbins_dv, dtrange_dv, detector): idx = row.name mask = masks_time[idx] time = row['ts'] @@ -463,7 +463,7 @@ def compute_parameters(row, masks_time, dst, emaps, fittype, nbins_dv, zrange_dv emaps = emaps, fittype = fittype, nbins_dv = nbins_dv, - zrange_dv = zrange_dv, + dtrange_dv = dtrange_dv, detector = detector) @@ -479,7 +479,7 @@ def compute_parameters(row, masks_time, dst, emaps, fittype, nbins_dv, zrange_dv # Apply the function to the previous df and compute krypton parameters for each time slice evol_pars = mask_time_df.apply(lambda row: compute_parameters(row, masks_time, dst, emaps, - fittype, nbins_dv, zrange_dv, detector), + fittype, nbins_dv, dtrange_dv, detector), axis=1) return evol_pars @@ -543,9 +543,48 @@ def cut_effs_evolution(masks_time : List[np.array], return evol_table_updated +def all_krevol(map, kdst, r_fid, nStimeprofile, + mask_s1, mask_s2, mask_band, fittype, + nbins_dv, dtrange_dv, detector): + + ''' COMPLETAR + + resumen: this function performs the whole proccess of evolution for all bins''' + + fid_sel = (kdst.R < r_fid) & mask_s1 & mask_s2 & mask_band + dstf = kdst[fid_sel] + t_start = dstf.time.min() + t_final = dstf.time.max() + ntimebins = np.max([int(np.floor((t_final - t_start) / nStimeprofile)), 1]) + + ts, masks_time = get_time_series_df(ntimebins = ntimebins, + time_range = (t_start, t_final), + dst = kdst) + + masks_timef = [mask[fid_sel] for mask in masks_time] + + evol_table = kr_time_evolution(ts = ts, + masks_time = masks_timef, + dst = dstf, + emaps = map, + fittype = fittype, + nbins_dv = nbins_dv, + dtrange_dv = dtrange_dv, + detector = detector) + + evol_table_eff = cut_effs_evolution(masks_time = masks_time, + dst = kdst, + mask_s1 = mask_s1, + mask_s2 = mask_s2, + mask_band = mask_band, + evol_table = evol_table) + + return evol_table_eff + + def add_krevol(r_fid : float, # Esto sería para meter en la ciudad de ICARO, flow.map(funcion, args, out) etc nStimeprofile : int, - **map_params): # PREGUNTA: Pongo aquí explícitamente todos los argumentos? O se pueden meter con un diccionario? + **map_params): ''' Adds the time evolution dataframe to the map. @@ -565,9 +604,9 @@ def add_krevol(r_fid : float, # Esto sería para meter en la ciudad de and returns the time evolution. ''' - def add_krevol(map, kdst, mask_s1, # Más de lo mismo con la pregunta anterior + def add_krevol(map, kdst, mask_s1, mask_s2, mask_band, fittype, - nbins_dv, zrange_dv, detector): + nbins_dv, dtrange_dv, detector): fid_sel = (kdst.R < r_fid) & mask_s1 & mask_s2 & mask_band dstf = kdst[fid_sel] @@ -590,7 +629,7 @@ def add_krevol(map, kdst, mask_s1, # Más de lo mismo con la pregunta emaps = map, fittype = fittype, nbins_dv = nbins_dv, - zrange_dv = zrange_dv, + dtrange_dv = dtrange_dv, detector = detector) evol_table_eff = cut_effs_evolution(masks_time = masks_time, From d01715f986cae1dfdd34ac6bbec210226117130c Mon Sep 17 00:00:00 2001 From: carhc Date: Fri, 11 Oct 2024 15:17:03 +0200 Subject: [PATCH 24/48] Cosmetics --- invisible_cities/reco/krmap_evolution.py | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/invisible_cities/reco/krmap_evolution.py b/invisible_cities/reco/krmap_evolution.py index 6665123422..9d886f01c7 100644 --- a/invisible_cities/reco/krmap_evolution.py +++ b/invisible_cities/reco/krmap_evolution.py @@ -175,9 +175,9 @@ def get_time_series_df(ntimebins : int, return shift_to_bin_centers(time_bins), masks -def compute_drift_v(dtdata : np.array, +def compute_drift_v(dtdata : np.array, nbins : int, - dtrange : Tuple[float, float], + dtrange : Tuple[float, float], seed : Tuple[float, float, float, float], detector : str)->Tuple[float, float]: @@ -263,7 +263,7 @@ def computing_kr_parameters(data : DataFrame, emaps : pd.DataFrame, fittype : KrFitFunction, nbins_dv : int, - dtrange_dv : List[float], + dtrange_dv : List[float], detector : str)->DataFrame: ''' @@ -325,7 +325,7 @@ def computing_kr_parameters(data : DataFrame, par, err, cov = transform_parameters(fit_output) - e0, lt = par + e0, lt = par e0u, ltu = err else: @@ -407,7 +407,7 @@ def kr_time_evolution(ts : np.array[float], emaps : pd.DataFrame, fittype : KrFitFunction, nbins_dv : int, - dtrange_dv : Tuple[float, float], + dtrange_dv : Tuple[float, float], detector : str) -> pd.DataFrame: ''' @@ -458,13 +458,13 @@ def compute_parameters(row, masks_time, dst, emaps, fittype, nbins_dv, dtrange_d sel_dst = dst[mask] # Compute the krypton parameters using the selected data - pars = computing_kr_parameters(data = sel_dst, - ts = time, - emaps = emaps, - fittype = fittype, - nbins_dv = nbins_dv, + pars = computing_kr_parameters(data = sel_dst, + ts = time, + emaps = emaps, + fittype = fittype, + nbins_dv = nbins_dv, dtrange_dv = dtrange_dv, - detector = detector) + detector = detector) @@ -629,7 +629,7 @@ def add_krevol(map, kdst, mask_s1, emaps = map, fittype = fittype, nbins_dv = nbins_dv, - dtrange_dv = dtrange_dv, + dtrange_dv = dtrange_dv, detector = detector) evol_table_eff = cut_effs_evolution(masks_time = masks_time, From fe215d7bb319f309a62216e20a8b9d62e1f3ce84 Mon Sep 17 00:00:00 2001 From: carhc Date: Mon, 14 Oct 2024 10:20:12 +0200 Subject: [PATCH 25/48] Upgrade explanations and delete useless line comments --- invisible_cities/reco/krmap_evolution.py | 143 +++++++++-------------- 1 file changed, 55 insertions(+), 88 deletions(-) diff --git a/invisible_cities/reco/krmap_evolution.py b/invisible_cities/reco/krmap_evolution.py index 9d886f01c7..174d2e9add 100644 --- a/invisible_cities/reco/krmap_evolution.py +++ b/invisible_cities/reco/krmap_evolution.py @@ -55,7 +55,9 @@ def gauss_seed(x : np.array, sigma_rel : Optional[int] = 0.05): ''' - This function estimates the seed for a gaussian fit. + This function estimates the seed for a gaussian fit. It looks for the + higher Y value and takes its corresponding X value as the center of the + Gaussian. Parameters ---------- @@ -72,11 +74,6 @@ def gauss_seed(x : np.array, Tuple with the seed estimation. ''' - # Looks for the higher Y value and takes its corresponding X value - # as the center of the Gaussian. Based on the "sigma_rel" parameter, - # applies some error to that X value and provides the estimation for the - # amplitude, the center and the sigma. - y_max = np.argmax(y) x_max = x[y_max] sigma = sigma_rel * x_max @@ -305,7 +302,7 @@ def computing_kr_parameters(data : DataFrame, norm_strat = NormStrategy.max) - x, y = select_fit_variables(fittype, data) # Importar de previo PR + x, y = select_fit_variables(fittype, data) corr = geo_correction_factor(data.X.to_numpy(), data.Y.to_numpy()) @@ -446,7 +443,6 @@ def kr_time_evolution(ts : np.array[float], parameters for a given time slice. ''' - # Create a DataFrame with the time slices (ts) mask_time_df = pd.DataFrame({'ts': ts}) def compute_parameters(row, masks_time, dst, emaps, fittype, nbins_dv, dtrange_dv, detector): @@ -454,10 +450,8 @@ def compute_parameters(row, masks_time, dst, emaps, fittype, nbins_dv, dtrange_d mask = masks_time[idx] time = row['ts'] - # Select the data for the given mask sel_dst = dst[mask] - # Compute the krypton parameters using the selected data pars = computing_kr_parameters(data = sel_dst, ts = time, emaps = emaps, @@ -466,17 +460,14 @@ def compute_parameters(row, masks_time, dst, emaps, fittype, nbins_dv, dtrange_d dtrange_dv = dtrange_dv, detector = detector) - - if pars is None or pars.empty: - return pd.Series(np.nan) # Return nan Series if no parameters found + return pd.Series(np.nan) if isinstance(pars, pd.DataFrame) and len(pars) == 1: - pars = pars.iloc[0] # Convert single-row DataFrame to Series + pars = pars.iloc[0] return pd.Series(pars) - # Apply the function to the previous df and compute krypton parameters for each time slice evol_pars = mask_time_df.apply(lambda row: compute_parameters(row, masks_time, dst, emaps, fittype, nbins_dv, dtrange_dv, detector), @@ -490,7 +481,7 @@ def cut_effs_evolution(masks_time : List[np.array], mask_s1 : np.array, mask_s2 : np.array, mask_band : np.array, - evol_table : pd.DataFrame): + evol_table : pd.DataFrame) -> pd.DataFrame: ''' Computes the efficiencies in time evolution for different time slices. @@ -543,102 +534,78 @@ def cut_effs_evolution(masks_time : List[np.array], return evol_table_updated -def all_krevol(map, kdst, r_fid, nStimeprofile, - mask_s1, mask_s2, mask_band, fittype, - nbins_dv, dtrange_dv, detector): +def all_krevol(emap : pd.DataFrame, + dst : pd.DataFrame, + r_fid : float, + nStimeprofile : int, + mask_s1 : np.array, + mask_s2 : np.array, + mask_band : np.array, + fittype : KrFitFunction, + nbins_dv : int, + dtrange_dv : Tuple[float, float], + detector : str) -> pd.DataFrame: + + ''' + Computes the whole krypton evolution parameters table. It applies some masks - ''' COMPLETAR + Parameters + ---------- + emap: pd.DataFrame + Map of e0, lt, etc. corrections for which the evolution is being computed. + dst: pd.DataFrame + Krypton dst containing the data to be analyzed. + r_fid: float + Fiducial radius for the cuts to exclude data outside the detector's inner region. + nStimeprofile: np.array + Duration (in seconds) of the time bins in which the evolution is computed. + mask_s1: np.array + Mask of S1 cut used to filter data. + mask_s2: np.array + Mask of S2 cut used to filter data. + mask_band: np.array + Mask of band cut used to filter data. + fittype: KrFitFunction + Type of fit to be used. + nbins_dv: int + Number of DT bins considered for the Drift Velocity computation + dtrange_dv: tuple + Range of DT considered for the Drift Velocity computation. + detector: str + Detector name. - resumen: this function performs the whole proccess of evolution for all bins''' + Returns + ------- + evol_table_updated: pd.DataFrame + Data frame containing Krypton evolution parameters including the efficiencies. + ''' - fid_sel = (kdst.R < r_fid) & mask_s1 & mask_s2 & mask_band - dstf = kdst[fid_sel] + fid_sel = (dst.R < r_fid) & mask_s1 & mask_s2 & mask_band + dstf = dst[fid_sel] t_start = dstf.time.min() t_final = dstf.time.max() ntimebins = np.max([int(np.floor((t_final - t_start) / nStimeprofile)), 1]) ts, masks_time = get_time_series_df(ntimebins = ntimebins, time_range = (t_start, t_final), - dst = kdst) + dst = dst) masks_timef = [mask[fid_sel] for mask in masks_time] evol_table = kr_time_evolution(ts = ts, masks_time = masks_timef, dst = dstf, - emaps = map, + emaps = emap, fittype = fittype, nbins_dv = nbins_dv, dtrange_dv = dtrange_dv, detector = detector) evol_table_eff = cut_effs_evolution(masks_time = masks_time, - dst = kdst, + dst = dst, mask_s1 = mask_s1, mask_s2 = mask_s2, mask_band = mask_band, evol_table = evol_table) return evol_table_eff - - -def add_krevol(r_fid : float, # Esto sería para meter en la ciudad de ICARO, flow.map(funcion, args, out) etc - nStimeprofile : int, - **map_params): - - ''' - Adds the time evolution dataframe to the map. - - Parameters - --------- - r_fid: float - Maximum radius for fiducial sample. - nStimeprofile: int - Number of seconds for each time bin. - map_params: dict - Dictionary containing the config file variables. - - Returns - --------- - Function which takes as input map, kr_data, and kr_mask - and returns the time evolution. - ''' - - def add_krevol(map, kdst, mask_s1, - mask_s2, mask_band, fittype, - nbins_dv, dtrange_dv, detector): - - fid_sel = (kdst.R < r_fid) & mask_s1 & mask_s2 & mask_band - dstf = kdst[fid_sel] - min_time = dstf.time.min() - max_time = dstf.time.max() - - ntimebins = get_number_of_time_bins(nStimeprofile = nStimeprofile, - tstart = min_time, - tfinal = max_time) - - ts, masks_time = get_time_series_df(ntimebins = ntimebins, - time_range = (min_time, max_time), - dst = kdst) - - masks_timef = [mask[fid_sel] for mask in masks_time] - - evol_table = kr_time_evolution(ts = ts, - masks_time = masks_timef, - dst = dstf, - emaps = map, - fittype = fittype, - nbins_dv = nbins_dv, - dtrange_dv = dtrange_dv, - detector = detector) - - evol_table_eff = cut_effs_evolution(masks_time = masks_time, - data = kdst, - mask_s1 = mask_s1, - mask_s2 = mask_s2, - mask_band = mask_band, - evol_table = evol_table) - - return evol_table_eff - - return add_krevol From ee2db0d73f005a7c03d71ed89f820d1f2d44e96b Mon Sep 17 00:00:00 2001 From: carhc Date: Mon, 14 Oct 2024 10:45:30 +0200 Subject: [PATCH 26/48] Rename get_fit_function_lt to get_function_and_seed_lt --- invisible_cities/reco/krmap_evolution.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/invisible_cities/reco/krmap_evolution.py b/invisible_cities/reco/krmap_evolution.py index 174d2e9add..59950336ad 100644 --- a/invisible_cities/reco/krmap_evolution.py +++ b/invisible_cities/reco/krmap_evolution.py @@ -14,7 +14,7 @@ from .. reco.corrections import apply_all_correction from .. core.stat_functions import poisson_sigma from .. database import load_db as DB -from .. reco.icaro_components import get_fit_function_lt, select_fit_variables, transform_parameters # Won't work until previous PR are approved +from .. reco.icaro_components import get_function_and_seed_lt, select_fit_variables, transform_parameters # Won't work until previous PR are approved def sigmoid(x : np.array, @@ -296,7 +296,7 @@ def computing_kr_parameters(data : DataFrame, # Computing E0, LT - fit_func, seed = get_fit_function_lt(fittype) + fit_func, seed = get_function_and_seed_lt(fittype) geo_correction_factor = e0_xy_correction(map = emaps, norm_strat = NormStrategy.max) From 1520822e00277ccaadf7641325ddcd03cd16de06 Mon Sep 17 00:00:00 2001 From: carhc Date: Mon, 14 Oct 2024 12:04:07 +0200 Subject: [PATCH 27/48] Remove blank spaces --- invisible_cities/reco/krmap_evolution.py | 23 ++++------------------- 1 file changed, 4 insertions(+), 19 deletions(-) diff --git a/invisible_cities/reco/krmap_evolution.py b/invisible_cities/reco/krmap_evolution.py index 59950336ad..8026de23f9 100644 --- a/invisible_cities/reco/krmap_evolution.py +++ b/invisible_cities/reco/krmap_evolution.py @@ -50,10 +50,7 @@ def sigmoid(x : np.array, return sigmoid -def gauss_seed(x : np.array, - y : np.array, - sigma_rel : Optional[int] = 0.05): - +def gauss_seed(x : np.array, y : np.array, sigma_rel : Optional[int] = 0.05): ''' This function estimates the seed for a gaussian fit. It looks for the higher Y value and takes its corresponding X value as the center of the @@ -83,9 +80,7 @@ def gauss_seed(x : np.array, return seed -def resolution(values : np.array, - errors : np.array): - +def resolution(values : np.array, errors : np.array): ''' Computes the resolution (% FWHM) from the Gaussian parameters. @@ -113,9 +108,7 @@ def resolution(values : np.array, return res, ures -def quick_gauss_fit(data : np.array, - bins : int): - +def quick_gauss_fit(data : np.array, bins : int): ''' This function histograms input data and then fits it to a Gaussian. @@ -144,7 +137,6 @@ def quick_gauss_fit(data : np.array, def get_time_series_df(ntimebins : int, time_range : Tuple[float, float], dst : DataFrame)->Tuple[np.array, List[np.array]]: - ''' Given a dst this function returns a time series (ts) and a list of masks which are used to divide the dst in time intervals. @@ -177,7 +169,6 @@ def compute_drift_v(dtdata : np.array, dtrange : Tuple[float, float], seed : Tuple[float, float, float, float], detector : str)->Tuple[float, float]: - ''' Computes the drift velocity for a given distribution using the sigmoid function to get the cathode edge. @@ -227,9 +218,7 @@ def compute_drift_v(dtdata : np.array, return dv, dvu -def e0_xy_correction(map : pd.DataFrame, - norm_strat : NormStrategy)->Callable: - +def e0_xy_correction(map : pd.DataFrame, norm_strat : NormStrategy)->Callable: ''' Provides the function to compute only the geometrical corrections. @@ -262,7 +251,6 @@ def computing_kr_parameters(data : DataFrame, nbins_dv : int, dtrange_dv : List[float], detector : str)->DataFrame: - ''' Computes some average parameters (e0, lt, drift v, energy resolution, S1w, S1h, S1e, S2w, S2h, S2e, S2q, Nsipm, 'Xrms, Yrms) @@ -406,7 +394,6 @@ def kr_time_evolution(ts : np.array[float], nbins_dv : int, dtrange_dv : Tuple[float, float], detector : str) -> pd.DataFrame: - ''' Computes some average parameters (e0, lt, drift v, S1w, S1h, S1e, S2w, S2h, S2e, S2q, Nsipm, Xrms, Yrms) @@ -482,7 +469,6 @@ def cut_effs_evolution(masks_time : List[np.array], mask_s2 : np.array, mask_band : np.array, evol_table : pd.DataFrame) -> pd.DataFrame: - ''' Computes the efficiencies in time evolution for different time slices. Returns the input DataFrame updated with S1eff, S2eff, Bandeff. @@ -545,7 +531,6 @@ def all_krevol(emap : pd.DataFrame, nbins_dv : int, dtrange_dv : Tuple[float, float], detector : str) -> pd.DataFrame: - ''' Computes the whole krypton evolution parameters table. It applies some masks From 36166d35a7ea68d23ffb2d1f4cac9dae7ad62aec Mon Sep 17 00:00:00 2001 From: carhc Date: Fri, 20 Dec 2024 09:46:25 +0100 Subject: [PATCH 28/48] Add sigma<0 protection in resolution function --- invisible_cities/reco/krmap_evolution.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/invisible_cities/reco/krmap_evolution.py b/invisible_cities/reco/krmap_evolution.py index 8026de23f9..ebfd43014e 100644 --- a/invisible_cities/reco/krmap_evolution.py +++ b/invisible_cities/reco/krmap_evolution.py @@ -102,6 +102,9 @@ def resolution(values : np.array, errors : np.array): amp , mu, sigma = values u_amp, u_mu, u_sigma = errors + if sigma < 0: + raise ValueError("Sigma cannot be negative") + res = 235.48 * sigma/mu ures = res * (u_mu**2/mu**2 + u_sigma**2/sigma**2)**0.5 @@ -197,7 +200,7 @@ def compute_drift_v(dtdata : np.array, 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) + 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 = DB.DetectorGeo(detector).ZMAX[0] From ac26eec8145127d1f53601deddb9445f56cedbeb Mon Sep 17 00:00:00 2001 From: carhc Date: Sun, 22 Dec 2024 19:27:44 +0100 Subject: [PATCH 29/48] Create krmap_evolution_test.py file with some tests --- invisible_cities/reco/krmap_evolution_test.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 invisible_cities/reco/krmap_evolution_test.py diff --git a/invisible_cities/reco/krmap_evolution_test.py b/invisible_cities/reco/krmap_evolution_test.py new file mode 100644 index 0000000000..e69de29bb2 From 63ea20ee1ffc6de056b716debe7c4b334c8b3962 Mon Sep 17 00:00:00 2001 From: carhc Date: Sun, 22 Dec 2024 19:32:57 +0100 Subject: [PATCH 30/48] Fix imports --- invisible_cities/reco/krmap_evolution.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/invisible_cities/reco/krmap_evolution.py b/invisible_cities/reco/krmap_evolution.py index ebfd43014e..d092e1c806 100644 --- a/invisible_cities/reco/krmap_evolution.py +++ b/invisible_cities/reco/krmap_evolution.py @@ -14,7 +14,7 @@ from .. reco.corrections import apply_all_correction from .. core.stat_functions import poisson_sigma from .. database import load_db as DB -from .. reco.icaro_components import get_function_and_seed_lt, select_fit_variables, transform_parameters # Won't work until previous PR are approved +from .. reco.krmap_functions import get_function_and_seed_lt, select_fit_variables, transform_parameters # Won't work until previous PR are approved def sigmoid(x : np.array, From 0d7bfe3cc908c77905ba2b2836fb4535f87ed3fe Mon Sep 17 00:00:00 2001 From: carhc Date: Sun, 22 Dec 2024 22:32:06 +0100 Subject: [PATCH 31/48] Fix function definition --- invisible_cities/reco/krmap_evolution.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/invisible_cities/reco/krmap_evolution.py b/invisible_cities/reco/krmap_evolution.py index d092e1c806..36b7377220 100644 --- a/invisible_cities/reco/krmap_evolution.py +++ b/invisible_cities/reco/krmap_evolution.py @@ -389,7 +389,7 @@ def computing_kr_parameters(data : DataFrame, return evol -def kr_time_evolution(ts : np.array[float], +def kr_time_evolution(ts : np.array, masks_time : List[np.array], dst : pd.DataFrame, emaps : pd.DataFrame, From 10db202eae12ef8788d6ff50bf13a74bccfb3206 Mon Sep 17 00:00:00 2001 From: carhc Date: Sun, 22 Dec 2024 22:37:37 +0100 Subject: [PATCH 32/48] Add tests in krmap_evolution_test.py --- invisible_cities/reco/krmap_evolution_test.py | 96 +++++++++++++++++++ 1 file changed, 96 insertions(+) diff --git a/invisible_cities/reco/krmap_evolution_test.py b/invisible_cities/reco/krmap_evolution_test.py index e69de29bb2..f7572fbfcb 100644 --- a/invisible_cities/reco/krmap_evolution_test.py +++ b/invisible_cities/reco/krmap_evolution_test.py @@ -0,0 +1,96 @@ +import numpy as np +import numpy.testing as npt + +from pytest import approx +from flaky import flaky + +from hypothesis import given +from hypothesis.strategies import floats +from .. reco.krmap_evolution import sigmoid, compute_drift_v, resolution + +from .. database import load_db as DB + +sensible_floats = floats(-1e4, +1e4) + + +def test_sigmoid_definition_fixed_values(): + + x = np.array(range(10)) + scale = 1 + slope = -5 + offset = 0 + inflection = 5 + + expected_output = scale / (1 + np.exp(-slope * (x - inflection))) + offset + actual_output = sigmoid(x, scale, inflection, slope, offset) + + npt.assert_allclose(actual_output, expected_output, rtol=1e-10) + + +def test_sigmoid_zero_slope(): + + x = np.array(range(10)) + scale = 1 + slope = 0 + offset = 0 + inflection = 0 + + expected_output = np.full(len(x), 0.5 * scale + offset) + actual_output = sigmoid(x, scale, inflection, slope, offset) + + npt.assert_allclose(actual_output, expected_output, rtol=1e-10) + + +@given(floats(min_value = 1e-2, max_value = 1e4), + floats(min_value = -1e2, max_value = 1e2)) +def test_sigmoid_limits(scale, offset): + + slope = 10 + inflection = 0 + aux_sigmoid = lambda x: sigmoid(x, scale, inflection, slope, offset) + + assert aux_sigmoid(-10) == approx(offset, rel=1e-2) + assert aux_sigmoid( 10) == approx(scale+offset, rel=1e-2) + + +@given(floats(min_value = 1e-2, max_value = 1e4), + floats(min_value = -10, max_value = 10), + floats(min_value = 1e-2, max_value = 1e2), + floats(min_value = -1e2, max_value = 1e2)) +def test_sigmoid_values_at_0(scale, inflection, slope, offset): + expected_output = scale / (1 + np.exp(slope*inflection)) + offset + actual_output = sigmoid(0, scale, inflection, slope, offset) + assert expected_output == actual_output + + +@flaky(max_runs=10, min_passes=9) +def test_compute_drift_v_when_moving_edge(): + edge = np.random.uniform(530, 570) + Nevents = 100 * 1000 + data = np.random.uniform(450, edge, Nevents) + data = np.random.normal(data, 1) + dv, dvu = compute_drift_v(data, 60, [500,600], + [1500, 550,1,0], 'next100') + dv_th = DB.DetectorGeo('next100').ZMAX[0]/edge + + assert dv_th == approx(dv, abs=5*dvu) + + +def test_compute_drift_v_failing_fit_return_nan(): + dst = np.random.rand(1000)*1200 + dv_vect = compute_drift_v(dst, 60, [500,600], [1500, 550, 1, 0], 'new') + nans = np.array([np.nan, np.nan]) + assert dv_vect == nans + + +def test_resolution_typical_case(): + + values = np.array([10, 100, 2]) + errors = np.array([0.1, 1, 0.05]) + expected_res = 235.48 * (values[2] / values[1]) # 235.48 * sigma/mu + expected_ures = expected_res * ((errors[1] / values[1])**2 + (errors[2] / values[2])**2)**0.5 + + res, ures = resolution(values, errors) + + assert res == approx(expected_res, rel=1e-5) + assert ures == approx(expected_ures, rel=1e-5) From 487612281d0d30125436e748a5523bee8c7ad624 Mon Sep 17 00:00:00 2001 From: carhc Date: Mon, 23 Dec 2024 09:18:23 +0100 Subject: [PATCH 33/48] Correct gauss_seed function and add test --- invisible_cities/reco/krmap_evolution.py | 12 ++++++------ invisible_cities/reco/krmap_evolution_test.py | 18 +++++++++++++++++- 2 files changed, 23 insertions(+), 7 deletions(-) diff --git a/invisible_cities/reco/krmap_evolution.py b/invisible_cities/reco/krmap_evolution.py index 36b7377220..2c27f18368 100644 --- a/invisible_cities/reco/krmap_evolution.py +++ b/invisible_cities/reco/krmap_evolution.py @@ -59,9 +59,9 @@ def gauss_seed(x : np.array, y : np.array, sigma_rel : Optional[int] = 0.05): Parameters ---------- x: np.array - Data to fit. - y: int - Number of bins for the histogram. + x data to fit. + y: np.array + y data to fit. sigma_rel (Optional): int Relative error, default 5%. @@ -71,9 +71,9 @@ def gauss_seed(x : np.array, y : np.array, sigma_rel : Optional[int] = 0.05): Tuple with the seed estimation. ''' - y_max = np.argmax(y) - x_max = x[y_max] - sigma = sigma_rel * x_max + y_max = max(y) + x_max = x[np.argmax(y)] + sigma = sigma_rel * (max(x)-min(x)) * 0.5 amp = y_max * (2 * np.pi)**0.5 * sigma * np.diff(x)[0] seed = amp, x_max, sigma diff --git a/invisible_cities/reco/krmap_evolution_test.py b/invisible_cities/reco/krmap_evolution_test.py index f7572fbfcb..d355806801 100644 --- a/invisible_cities/reco/krmap_evolution_test.py +++ b/invisible_cities/reco/krmap_evolution_test.py @@ -6,7 +6,7 @@ from hypothesis import given from hypothesis.strategies import floats -from .. reco.krmap_evolution import sigmoid, compute_drift_v, resolution +from .. reco.krmap_evolution import sigmoid, gauss_seed, compute_drift_v, resolution from .. database import load_db as DB @@ -63,6 +63,22 @@ def test_sigmoid_values_at_0(scale, inflection, slope, offset): assert expected_output == actual_output +@given(floats(min_value = 1e-3, max_value = 1)) +def test_gauss_seed_all_sigmas(sigma_rel): + + x = np.array(range(10)) + y = np.array([1, 3, 7, 12, 18, 25, 18, 12, 7, 3]) + + max_y = np.max(y) + max_x = x[np.argmax(y)] + + expected_amp = max_y * (2 * np.pi)**0.5 * (sigma_rel * (max(x)-min(x))*0.5) + expected_seed = (expected_amp, max_x, sigma_rel * (max(x)-min(x))*0.5) + + actual_seed = gauss_seed(x, y, sigma_rel=sigma_rel) + npt.assert_allclose(actual_seed, expected_seed, rtol=1e-5) + + @flaky(max_runs=10, min_passes=9) def test_compute_drift_v_when_moving_edge(): edge = np.random.uniform(530, 570) From acd115e1b7c1869869d5c39f64576cf0f745e8ca Mon Sep 17 00:00:00 2001 From: carhc Date: Mon, 23 Dec 2024 10:42:26 +0100 Subject: [PATCH 34/48] Add tests in krmap_evolution_test.py --- invisible_cities/reco/krmap_evolution_test.py | 66 +++++++++++++++---- 1 file changed, 55 insertions(+), 11 deletions(-) diff --git a/invisible_cities/reco/krmap_evolution_test.py b/invisible_cities/reco/krmap_evolution_test.py index d355806801..ed43db1f19 100644 --- a/invisible_cities/reco/krmap_evolution_test.py +++ b/invisible_cities/reco/krmap_evolution_test.py @@ -1,16 +1,19 @@ import numpy as np +import pandas as pd import numpy.testing as npt -from pytest import approx +from pytest import approx, mark from flaky import flaky from hypothesis import given -from hypothesis.strategies import floats +from hypothesis.strategies import floats, integers from .. reco.krmap_evolution import sigmoid, gauss_seed, compute_drift_v, resolution +from .. reco.krmap_evolution import quick_gauss_fit, get_time_series_df +from .. evm.ic_containers import FitFunction from .. database import load_db as DB - -sensible_floats = floats(-1e4, +1e4) +from .. io.dst_io import load_dst +from .. core.testing_utils import assert_dataframes_equal def test_sigmoid_definition_fixed_values(): @@ -79,27 +82,68 @@ def test_gauss_seed_all_sigmas(sigma_rel): npt.assert_allclose(actual_seed, expected_seed, rtol=1e-5) +@given(floats (min_value=-1e2, max_value=1e4), + floats (min_value=0.1, max_value=5), + integers(min_value=1e2, max_value=1e5), + integers(min_value=10, max_value = 1e2)) +def test_quick_gauss_fit(mean, sigma, size, bins): + + def generate_gaussian_data(mean, sigma, size): + return np.random.normal(loc=mean, scale=sigma, size=size) + + data = generate_gaussian_data(mean, sigma, size) + fit_result = quick_gauss_fit(data, bins) + + assert type(fit_result) is FitFunction + assert np.isclose(fit_result.values[1], mean, rtol=0.2) + assert np.isclose(fit_result.values[2], sigma, rtol=0.2) + + @flaky(max_runs=10, min_passes=9) def test_compute_drift_v_when_moving_edge(): edge = np.random.uniform(530, 570) Nevents = 100 * 1000 data = np.random.uniform(450, edge, Nevents) data = np.random.normal(data, 1) - dv, dvu = compute_drift_v(data, 60, [500,600], - [1500, 550,1,0], 'next100') - dv_th = DB.DetectorGeo('next100').ZMAX[0]/edge + dv, dvu = compute_drift_v(dtdata = data, nbins = 60, dtrange = [500, 600], + seed = [1500, 550, 1, 0], detector ='new') + dv_th = DB.DetectorGeo('new').ZMAX[0]/edge assert dv_th == approx(dv, abs=5*dvu) def test_compute_drift_v_failing_fit_return_nan(): dst = np.random.rand(1000)*1200 - dv_vect = compute_drift_v(dst, 60, [500,600], [1500, 550, 1, 0], 'new') - nans = np.array([np.nan, np.nan]) - assert dv_vect == nans + dv_vect = compute_drift_v(dtdata = dst, nbins = 50, dtrange = [500, 600], + seed = [1500, 550, 1, 0], detector ='new') + + assert np.all(np.isnan(dv_vect)) + + +def test_get_time_series_df_empty_data(): + + dst = pd.DataFrame({'time': []}) + ntimebins = 3 + time_range = (0, 10) + ts, masks = get_time_series_df(ntimebins, time_range, dst) + expected_ts = np.linspace(0, 10, ntimebins + 1)[:-1] + (10 / (2 * ntimebins)) + + npt.assert_allclose(ts, expected_ts, rtol=1e-5) + + for mask in masks: + assert np.sum(mask) == 0 + + +@mark.parametrize('n', [2, 5, 10]) +def test_get_time_series_df_pandas(KrMC_kdst, n): + kdst = load_dst(KrMC_kdst, 'DST', 'Events') + time_indx = pd.cut(kdst.time, n, labels=np.arange(n)) + ts, masks = get_time_series_df(n, (kdst.time.min(), kdst.time.max()), kdst) + for i, mask in enumerate(masks): + assert_dataframes_equal(kdst[mask], kdst[time_indx==i]) -def test_resolution_typical_case(): +def test_resolution_definition(): values = np.array([10, 100, 2]) errors = np.array([0.1, 1, 0.05]) From bfb4956462772f729f8ed6797ae9db5f0a9afd79 Mon Sep 17 00:00:00 2001 From: carhc Date: Mon, 23 Dec 2024 11:53:53 +0100 Subject: [PATCH 35/48] Uodate test_quick_gauss_fit --- invisible_cities/reco/krmap_evolution_test.py | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/invisible_cities/reco/krmap_evolution_test.py b/invisible_cities/reco/krmap_evolution_test.py index ed43db1f19..87daf53f5a 100644 --- a/invisible_cities/reco/krmap_evolution_test.py +++ b/invisible_cities/reco/krmap_evolution_test.py @@ -82,21 +82,26 @@ def test_gauss_seed_all_sigmas(sigma_rel): npt.assert_allclose(actual_seed, expected_seed, rtol=1e-5) -@given(floats (min_value=-1e2, max_value=1e4), - floats (min_value=0.1, max_value=5), - integers(min_value=1e2, max_value=1e5), - integers(min_value=10, max_value = 1e2)) +@mark.parametrize('mean sigma size bins'.split(), +( (0, 1, 1000, 50), + (10, 2, 10000, 100), +(-100, 1, 5000, 75),) +) def test_quick_gauss_fit(mean, sigma, size, bins): def generate_gaussian_data(mean, sigma, size): return np.random.normal(loc=mean, scale=sigma, size=size) - data = generate_gaussian_data(mean, sigma, size) + data = generate_gaussian_data(mean, sigma, size) + fit_result = quick_gauss_fit(data, bins) + print(fit_result.values[1]) + print(mean) + assert type(fit_result) is FitFunction - assert np.isclose(fit_result.values[1], mean, rtol=0.2) - assert np.isclose(fit_result.values[2], sigma, rtol=0.2) + assert np.isclose(fit_result.values[2], sigma, rtol=0.1) + assert np.isclose(fit_result.values[1], mean, atol=0.1) @flaky(max_runs=10, min_passes=9) From ea1fb6701a814b1b2adeaf5bb3a99d7b165081fb Mon Sep 17 00:00:00 2001 From: carhc Date: Mon, 23 Dec 2024 12:29:14 +0100 Subject: [PATCH 36/48] Update test_compute_drift_v_failing_fit_return_nan --- invisible_cities/reco/krmap_evolution_test.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/invisible_cities/reco/krmap_evolution_test.py b/invisible_cities/reco/krmap_evolution_test.py index 87daf53f5a..50e36596f6 100644 --- a/invisible_cities/reco/krmap_evolution_test.py +++ b/invisible_cities/reco/krmap_evolution_test.py @@ -118,10 +118,10 @@ def test_compute_drift_v_when_moving_edge(): def test_compute_drift_v_failing_fit_return_nan(): - dst = np.random.rand(1000)*1200 + '''Generating nonsense data to make sigmoid fit fail''' + dst = [500, 599] dv_vect = compute_drift_v(dtdata = dst, nbins = 50, dtrange = [500, 600], seed = [1500, 550, 1, 0], detector ='new') - assert np.all(np.isnan(dv_vect)) From c4fb5a499fc2bf4161bd68f7fe9f2821b7743e94 Mon Sep 17 00:00:00 2001 From: carhc Date: Mon, 23 Dec 2024 13:01:29 +0100 Subject: [PATCH 37/48] Update test_get_time_series_df_pandas --- invisible_cities/reco/krmap_evolution_test.py | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/invisible_cities/reco/krmap_evolution_test.py b/invisible_cities/reco/krmap_evolution_test.py index 50e36596f6..ca74499eff 100644 --- a/invisible_cities/reco/krmap_evolution_test.py +++ b/invisible_cities/reco/krmap_evolution_test.py @@ -138,14 +138,19 @@ def test_get_time_series_df_empty_data(): for mask in masks: assert np.sum(mask) == 0 - +@mark.parametrize('reverse', [False, True]) @mark.parametrize('n', [2, 5, 10]) -def test_get_time_series_df_pandas(KrMC_kdst, n): - kdst = load_dst(KrMC_kdst, 'DST', 'Events') - time_indx = pd.cut(kdst.time, n, labels=np.arange(n)) +def test_get_time_series_df_pandas(n, reverse): + + time = np.arange(n-1, -1, -1) if reverse else np.arange(n) + kdst = pd.DataFrame({'time': time}) ts, masks = get_time_series_df(n, (kdst.time.min(), kdst.time.max()), kdst) + assert len(ts) == len(masks) == n, f'{len(ts)}, {len(masks)}, {n}' for i, mask in enumerate(masks): - assert_dataframes_equal(kdst[mask], kdst[time_indx==i]) + assert np.count_nonzero(mask) == 1 + index = n-1-i if reverse else i + assert mask[index] == True + assert np.all(ts[:-1] < ts[1:]) def test_resolution_definition(): From e9e18f9801d8d5bb7482913ba3cd83b08ead9e4d Mon Sep 17 00:00:00 2001 From: carhc Date: Mon, 23 Dec 2024 13:30:06 +0100 Subject: [PATCH 38/48] Update test_resolution_definition --- invisible_cities/reco/krmap_evolution_test.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/invisible_cities/reco/krmap_evolution_test.py b/invisible_cities/reco/krmap_evolution_test.py index ca74499eff..03e46df9a3 100644 --- a/invisible_cities/reco/krmap_evolution_test.py +++ b/invisible_cities/reco/krmap_evolution_test.py @@ -154,11 +154,12 @@ def test_get_time_series_df_pandas(n, reverse): def test_resolution_definition(): - - values = np.array([10, 100, 2]) - errors = np.array([0.1, 1, 0.05]) - expected_res = 235.48 * (values[2] / values[1]) # 235.48 * sigma/mu - expected_ures = expected_res * ((errors[1] / values[1])**2 + (errors[2] / values[2])**2)**0.5 + expected_res = 3 # % FWHM + mean = 123 + sigma = mean*expected_res/235.48 + values = np.array([np.nan, mean, sigma]) + errors = np.array([np.nan, 0.33*mean, 0.33*sigma]) + expected_ures = expected_res * 0.33 * np.sqrt(2) res, ures = resolution(values, errors) From 548db286bdc40af22f3584ca53e59a36d1e3f2a8 Mon Sep 17 00:00:00 2001 From: carhc Date: Thu, 26 Dec 2024 11:05:36 +0100 Subject: [PATCH 39/48] Remove unused imports --- invisible_cities/reco/krmap_evolution_test.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/invisible_cities/reco/krmap_evolution_test.py b/invisible_cities/reco/krmap_evolution_test.py index 03e46df9a3..abeed49d1e 100644 --- a/invisible_cities/reco/krmap_evolution_test.py +++ b/invisible_cities/reco/krmap_evolution_test.py @@ -12,8 +12,6 @@ from .. evm.ic_containers import FitFunction from .. database import load_db as DB -from .. io.dst_io import load_dst -from .. core.testing_utils import assert_dataframes_equal def test_sigmoid_definition_fixed_values(): From 642533445b7b5d2726c096cac17ef83b2aeeca22 Mon Sep 17 00:00:00 2001 From: carhc Date: Thu, 26 Dec 2024 11:06:03 +0100 Subject: [PATCH 40/48] Cosmetics --- invisible_cities/reco/krmap_evolution.py | 12 +++++------- invisible_cities/reco/krmap_evolution_test.py | 6 +++--- 2 files changed, 8 insertions(+), 10 deletions(-) diff --git a/invisible_cities/reco/krmap_evolution.py b/invisible_cities/reco/krmap_evolution.py index 2c27f18368..c7f83608b0 100644 --- a/invisible_cities/reco/krmap_evolution.py +++ b/invisible_cities/reco/krmap_evolution.py @@ -303,21 +303,19 @@ def computing_kr_parameters(data : DataFrame, # this would be reduced to just the following: # y_corr = y*geo_correction_factor(data.X.to_numpy(), data.Y.to_numpy()) - fit_output, _, _, _ = fit(func = fit_func, # Misma funcion que en el ajuste del mapa - x = x, - y = y_corr, - seed = seed(x, y_corr), - full_output = False) + fit_output = fit(func = fit_func, # Misma funcion que en el ajuste del mapa + x = x, + y = y_corr, + seed = seed(x, y_corr), + full_output = False) if fittype == KrFitFunction.log_lin: - par, err, cov = transform_parameters(fit_output) e0, lt = par e0u, ltu = err else: - e0, lt = fit_output.values e0u, ltu = fit_output.errors diff --git a/invisible_cities/reco/krmap_evolution_test.py b/invisible_cities/reco/krmap_evolution_test.py index abeed49d1e..e5d0c16fd7 100644 --- a/invisible_cities/reco/krmap_evolution_test.py +++ b/invisible_cities/reco/krmap_evolution_test.py @@ -81,9 +81,9 @@ def test_gauss_seed_all_sigmas(sigma_rel): @mark.parametrize('mean sigma size bins'.split(), -( (0, 1, 1000, 50), - (10, 2, 10000, 100), -(-100, 1, 5000, 75),) +( (0, 1, 1000, 50), + (10, 2, 10000, 100), +(-100, 1, 5000, 75),) ) def test_quick_gauss_fit(mean, sigma, size, bins): From d627c8a621a4186fe35517896b87af4c612f5cee Mon Sep 17 00:00:00 2001 From: carhc Date: Fri, 27 Dec 2024 12:07:59 +0100 Subject: [PATCH 41/48] Create fixtures and add test for computing_kr_parameters --- invisible_cities/reco/krmap_evolution_test.py | 100 +++++++++++++++++- 1 file changed, 98 insertions(+), 2 deletions(-) diff --git a/invisible_cities/reco/krmap_evolution_test.py b/invisible_cities/reco/krmap_evolution_test.py index e5d0c16fd7..de583c4360 100644 --- a/invisible_cities/reco/krmap_evolution_test.py +++ b/invisible_cities/reco/krmap_evolution_test.py @@ -2,15 +2,16 @@ import pandas as pd import numpy.testing as npt -from pytest import approx, mark +from pytest import approx, mark, fixture from flaky import flaky from hypothesis import given from hypothesis.strategies import floats, integers from .. reco.krmap_evolution import sigmoid, gauss_seed, compute_drift_v, resolution -from .. reco.krmap_evolution import quick_gauss_fit, get_time_series_df +from .. reco.krmap_evolution import quick_gauss_fit, get_time_series_df, computing_kr_parameters from .. evm.ic_containers import FitFunction +from .. types.symbols import KrFitFunction from .. database import load_db as DB @@ -163,3 +164,98 @@ def test_resolution_definition(): assert res == approx(expected_res, rel=1e-5) assert ures == approx(expected_ures, rel=1e-5) + +@fixture +def dummy_kr_map(): + return pd.DataFrame(dict(bin = np.arange(10), + counts = 1, + e0 = 1, + ue0 = 0, + lt = 1, + ult = 0, + cov = 0, + res_std = 0, + pval = 1, + in_fid = True, + has_min_counts = True, + fit_success = True, + valid = True, + R = 1, + x_index = 1, + y_index = 1, + Y = 1, + X = 1)) + +@fixture +def dummy_kr_dst(): + n = 10 + res = 5.5 # % FWHM + lt = 1500 + dv = 1 + X = np.linspace(-200, 200, n) + Y = X + Z = np.linspace(400, 500, n) + + X, Y, Z = map(np.ravel, np.meshgrid(X, Y, Z, indexing='ij')) + + DT = Z/dv + time = np.linspace(0, 1, n**3) + event = np.arange(n**3) + mu = 15000 + sigma = res*mu/235.48 + energy = np.random.normal(mu, sigma, n**3) + energy = energy*np.exp(-DT/lt) + # FIJAR SEED + pars = dict(res = res, lt = lt, mu = mu, sigma = sigma, dv = dv) + + dst = pd.DataFrame(dict(S1t = 0, S1w = 1, S1h = 2, S1e = 3, S2t = 0, S2w = 4, S2h=5, S2e = energy, + S2q = 6, Nsipm = 7, Xrms = 0, Yrms = 0, X = X, Y = Y, DT = DT, time = time, + event = event, Z = Z, Zrms = 0, R = 0, Phi = 0, nS1 = 1, nS2 = 1, s1_peak = 0, + s2_peak = 0, qmax = 8)) + return dst, pars + + +def test_dummy_fixtures(dummy_kr_dst, dummy_kr_map): + + dst, pars_true = dummy_kr_dst + pars = computing_kr_parameters(data = dst, ts = 0.5, + emaps = dummy_kr_map, + fittype = KrFitFunction.log_lin, + nbins_dv = 10, + dtrange_dv = [470, 470+10*10], + detector = 'new') + + expected = dst.mean() + data_size = len(dst) + + assert pars['ts'] [0] == 0.5 + assert pars['s1w'] [0] == expected.S1w + assert pars['s1wu'] [0] == 0 + assert pars['s1h'] [0] == expected.S1h + assert pars['s1hu'] [0] == 0 + assert pars['s1e'] [0] == expected.S1e + assert pars['s1eu'] [0] == 0 + assert pars['s2w'] [0] == expected.S2w + assert pars['s2wu'] [0] == 0 + assert pars['s2h'] [0] == expected.S2h + assert pars['s2hu'] [0] == 0 + assert pars['s2e'] [0] == expected.S2e + assert pars['s2q'] [0] == expected.S2q + assert pars['s2qu'] [0] == 0 + assert pars['Nsipm'] [0] == expected.Nsipm + assert pars['Nsipmu'][0] == 0 + assert pars['Xrms'] [0] == expected.Xrms + assert pars['Xrmsu'] [0] == 0 + assert pars['Yrms'] [0] == expected.Yrms + assert pars['Yrmsu'] [0] == 0 + + assert np.isclose(pars['s2eu'][0], dst.S2e.std()/data_size**0.5, rtol=1e-1) + assert np.isclose(pars['e0'][0], pars_true['mu'], rtol=5e-2) # E0 + assert np.isclose(pars['lt'][0], pars_true['lt'], rtol=1e-1) # LT + assert np.isclose(pars['dv'][0], pars_true['dv'], rtol=1e-1) + assert np.isclose(pars['resol'][0], pars_true['res'], rtol=5e-2) + assert pars['e0u'][0]/pars['e0'][0] <= 0.1 + assert pars['ltu'][0]/pars['lt'][0] <= 0.1 + assert pars['resolu'][0]/pars['resol'][0] <= 0.1 + # assert pars['dvu'][0]/pars['dv'][0] <= 0.1 Uncertainty in dv for this dataset is too high + # but the validity of the drift velocity computation is tested elsewhere From 0c071848aa9a564b077dd38ff06e90c1beb64338 Mon Sep 17 00:00:00 2001 From: carhc Date: Fri, 27 Dec 2024 12:09:10 +0100 Subject: [PATCH 42/48] Update and correct docstrings --- invisible_cities/reco/krmap_evolution.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/invisible_cities/reco/krmap_evolution.py b/invisible_cities/reco/krmap_evolution.py index c7f83608b0..fc4c6e4e9a 100644 --- a/invisible_cities/reco/krmap_evolution.py +++ b/invisible_cities/reco/krmap_evolution.py @@ -256,7 +256,7 @@ def computing_kr_parameters(data : DataFrame, detector : str)->DataFrame: ''' Computes some average parameters (e0, lt, drift v, energy - resolution, S1w, S1h, S1e, S2w, S2h, S2e, S2q, Nsipm, 'Xrms, Yrms) + resolution, S1w, S1h, S1e, S2w, S2h, S2e, S2q, Nsipm, Xrms, Yrms) for a given krypton distribution. Parameters @@ -264,7 +264,7 @@ def computing_kr_parameters(data : DataFrame, data: DataFrame Kdst distribution to analyze. ts: float - Central time of the distribution. + Central time of the time bin. emaps: correction map Allows geometrical correction of the energy. fittype: KrFitFunction From 8b3ec759843334f0b3bb0b85f31e46e9494413eb Mon Sep 17 00:00:00 2001 From: carhc Date: Fri, 27 Dec 2024 12:10:04 +0100 Subject: [PATCH 43/48] Add temporary utils for debugging while testing --- invisible_cities/reco/krmap_evolution.py | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/invisible_cities/reco/krmap_evolution.py b/invisible_cities/reco/krmap_evolution.py index fc4c6e4e9a..0b3367ed23 100644 --- a/invisible_cities/reco/krmap_evolution.py +++ b/invisible_cities/reco/krmap_evolution.py @@ -204,7 +204,7 @@ def compute_drift_v(dtdata : np.array, # At the moment there is not NEXT-100 DB so this won't work for that geometry z_cathode = DB.DetectorGeo(detector).ZMAX[0] - + print(seed) try: f = fit(sigmoid, x, y, seed, sigma = poisson_sigma(y), fit_range = dtrange) @@ -289,6 +289,12 @@ def computing_kr_parameters(data : DataFrame, fit_func, seed = get_function_and_seed_lt(fittype) + # TEMPORARY: Remove after + def e0_xy_correction(map, norm_strat): + def function(x, y): + return np.ones_like(x) + return function + geo_correction_factor = e0_xy_correction(map = emaps, norm_strat = NormStrategy.max) @@ -306,8 +312,7 @@ def computing_kr_parameters(data : DataFrame, fit_output = fit(func = fit_func, # Misma funcion que en el ajuste del mapa x = x, y = y_corr, - seed = seed(x, y_corr), - full_output = False) + seed = seed(x, y_corr)) if fittype == KrFitFunction.log_lin: par, err, cov = transform_parameters(fit_output) @@ -329,6 +334,10 @@ def computing_kr_parameters(data : DataFrame, # Computing Resolution + def apply_all_correction(maps, apply_temp): + def function(x, y, z, t): + return np.exp(z/1500) + return function tot_corr_factor = apply_all_correction(maps = emaps, apply_temp = False) From e6e5bbb54e43c15cfdcb214de76c840263d5b965 Mon Sep 17 00:00:00 2001 From: carhc Date: Fri, 27 Dec 2024 13:22:35 +0100 Subject: [PATCH 44/48] Modify all_krevol function --- invisible_cities/reco/krmap_evolution.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/invisible_cities/reco/krmap_evolution.py b/invisible_cities/reco/krmap_evolution.py index 0b3367ed23..8239b999c1 100644 --- a/invisible_cities/reco/krmap_evolution.py +++ b/invisible_cities/reco/krmap_evolution.py @@ -579,8 +579,10 @@ def all_krevol(emap : pd.DataFrame, dstf = dst[fid_sel] t_start = dstf.time.min() t_final = dstf.time.max() - ntimebins = np.max([int(np.floor((t_final - t_start) / nStimeprofile)), 1]) + # ntimebins = np.max([int(np.floor((t_final - t_start) / nStimeprofile)), 1]) + # np.floor yet to be decided + ntimebins = np.max([int((t_final - t_start) / nStimeprofile), 1]) ts, masks_time = get_time_series_df(ntimebins = ntimebins, time_range = (t_start, t_final), dst = dst) From 55e78500dedbc9087ac23eb08150c88d26a9bad9 Mon Sep 17 00:00:00 2001 From: carhc Date: Fri, 27 Dec 2024 13:23:26 +0100 Subject: [PATCH 45/48] Update and rename test --- invisible_cities/reco/krmap_evolution_test.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/invisible_cities/reco/krmap_evolution_test.py b/invisible_cities/reco/krmap_evolution_test.py index de583c4360..4a3b105f27 100644 --- a/invisible_cities/reco/krmap_evolution_test.py +++ b/invisible_cities/reco/krmap_evolution_test.py @@ -215,7 +215,7 @@ def dummy_kr_dst(): return dst, pars -def test_dummy_fixtures(dummy_kr_dst, dummy_kr_map): +def test_computing_kr_parameters(dummy_kr_dst, dummy_kr_map): dst, pars_true = dummy_kr_dst pars = computing_kr_parameters(data = dst, ts = 0.5, @@ -253,7 +253,7 @@ def test_dummy_fixtures(dummy_kr_dst, dummy_kr_map): assert np.isclose(pars['e0'][0], pars_true['mu'], rtol=5e-2) # E0 assert np.isclose(pars['lt'][0], pars_true['lt'], rtol=1e-1) # LT assert np.isclose(pars['dv'][0], pars_true['dv'], rtol=1e-1) - assert np.isclose(pars['resol'][0], pars_true['res'], rtol=5e-2) + assert np.isclose(pars['resol'][0], pars_true['res'], rtol=1e-1) assert pars['e0u'][0]/pars['e0'][0] <= 0.1 assert pars['ltu'][0]/pars['lt'][0] <= 0.1 assert pars['resolu'][0]/pars['resol'][0] <= 0.1 From af2b995fcaaa0e6a3eedc4241c1e4559d5ee8a72 Mon Sep 17 00:00:00 2001 From: carhc Date: Fri, 27 Dec 2024 13:24:16 +0100 Subject: [PATCH 46/48] Add test_kr_time_evolution function --- invisible_cities/reco/krmap_evolution_test.py | 25 +++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/invisible_cities/reco/krmap_evolution_test.py b/invisible_cities/reco/krmap_evolution_test.py index 4a3b105f27..a15697a70d 100644 --- a/invisible_cities/reco/krmap_evolution_test.py +++ b/invisible_cities/reco/krmap_evolution_test.py @@ -9,6 +9,7 @@ from hypothesis.strategies import floats, integers from .. reco.krmap_evolution import sigmoid, gauss_seed, compute_drift_v, resolution from .. reco.krmap_evolution import quick_gauss_fit, get_time_series_df, computing_kr_parameters +from .. reco.krmap_evolution import kr_time_evolution from .. evm.ic_containers import FitFunction from .. types.symbols import KrFitFunction @@ -259,3 +260,27 @@ def test_computing_kr_parameters(dummy_kr_dst, dummy_kr_map): assert pars['resolu'][0]/pars['resol'][0] <= 0.1 # assert pars['dvu'][0]/pars['dv'][0] <= 0.1 Uncertainty in dv for this dataset is too high # but the validity of the drift velocity computation is tested elsewhere + +def test_kr_time_evolution(dummy_kr_dst, dummy_kr_map): + dst_1 = dummy_kr_dst[0] + dst_2 = dst_1.copy() + deltat = np.diff(dst_1.time)[0] + dst_2.time += deltat + max(dst_1.time) + dst = pd.concat([dst_1, dst_2], ignore_index=True) + t0 = min(dst.time) + t1 = max(dst.time) + + ts, masks_time = get_time_series_df(ntimebins = 2, + time_range = (t0, t1), + dst = dst) + + evol_pars = kr_time_evolution(ts = ts, masks_time = masks_time, dst = dst, + emaps = dummy_kr_map, fittype = KrFitFunction.log_lin, + nbins_dv = 10, dtrange_dv = [470, 570], detector = 'new') + + assert len(evol_pars) == 2 + for col in evol_pars.columns: + if col == 'ts': assert (evol_pars.ts == ts).all() + else: assert np.all(np.diff(evol_pars[col].values, axis=0) == 0), col + + From 2aa7deaba58ff4b9e8a1fabcc35e00f691eae233 Mon Sep 17 00:00:00 2001 From: carhc Date: Fri, 27 Dec 2024 13:43:45 +0100 Subject: [PATCH 47/48] Remove unused imports + Cosmetics --- invisible_cities/reco/krmap_evolution_test.py | 40 +++++++++---------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/invisible_cities/reco/krmap_evolution_test.py b/invisible_cities/reco/krmap_evolution_test.py index a15697a70d..b042aa6aa5 100644 --- a/invisible_cities/reco/krmap_evolution_test.py +++ b/invisible_cities/reco/krmap_evolution_test.py @@ -6,7 +6,7 @@ from flaky import flaky from hypothesis import given -from hypothesis.strategies import floats, integers +from hypothesis.strategies import floats from .. reco.krmap_evolution import sigmoid, gauss_seed, compute_drift_v, resolution from .. reco.krmap_evolution import quick_gauss_fit, get_time_series_df, computing_kr_parameters from .. reco.krmap_evolution import kr_time_evolution @@ -168,24 +168,24 @@ def test_resolution_definition(): @fixture def dummy_kr_map(): - return pd.DataFrame(dict(bin = np.arange(10), - counts = 1, - e0 = 1, - ue0 = 0, - lt = 1, - ult = 0, - cov = 0, - res_std = 0, - pval = 1, + return pd.DataFrame(dict(bin = np.arange(10), + counts = 1, + e0 = 1, + ue0 = 0, + lt = 1, + ult = 0, + cov = 0, + res_std = 0, + pval = 1, in_fid = True, has_min_counts = True, - fit_success = True, - valid = True, - R = 1, - x_index = 1, - y_index = 1, - Y = 1, - X = 1)) + fit_success = True, + valid = True, + R = 1, + x_index = 1, + y_index = 1, + Y = 1, + X = 1)) @fixture def dummy_kr_dst(): @@ -251,9 +251,9 @@ def test_computing_kr_parameters(dummy_kr_dst, dummy_kr_map): assert pars['Yrmsu'] [0] == 0 assert np.isclose(pars['s2eu'][0], dst.S2e.std()/data_size**0.5, rtol=1e-1) - assert np.isclose(pars['e0'][0], pars_true['mu'], rtol=5e-2) # E0 - assert np.isclose(pars['lt'][0], pars_true['lt'], rtol=1e-1) # LT - assert np.isclose(pars['dv'][0], pars_true['dv'], rtol=1e-1) + assert np.isclose(pars['e0'][0], pars_true['mu'], rtol=5e-2) # E0 + assert np.isclose(pars['lt'][0], pars_true['lt'], rtol=1e-1) # LT + assert np.isclose(pars['dv'][0], pars_true['dv'], rtol=1e-1) assert np.isclose(pars['resol'][0], pars_true['res'], rtol=1e-1) assert pars['e0u'][0]/pars['e0'][0] <= 0.1 assert pars['ltu'][0]/pars['lt'][0] <= 0.1 From 84d19ee7533fdf01e6d25f40d682f0e15f071f73 Mon Sep 17 00:00:00 2001 From: carhc Date: Fri, 3 Jan 2025 12:02:22 +0100 Subject: [PATCH 48/48] Add test_cut_effs_evolution function --- invisible_cities/reco/krmap_evolution_test.py | 35 ++++++++++++++++++- 1 file changed, 34 insertions(+), 1 deletion(-) diff --git a/invisible_cities/reco/krmap_evolution_test.py b/invisible_cities/reco/krmap_evolution_test.py index b042aa6aa5..d9ddc45a05 100644 --- a/invisible_cities/reco/krmap_evolution_test.py +++ b/invisible_cities/reco/krmap_evolution_test.py @@ -7,9 +7,10 @@ from hypothesis import given from hypothesis.strategies import floats + from .. reco.krmap_evolution import sigmoid, gauss_seed, compute_drift_v, resolution from .. reco.krmap_evolution import quick_gauss_fit, get_time_series_df, computing_kr_parameters -from .. reco.krmap_evolution import kr_time_evolution +from .. reco.krmap_evolution import kr_time_evolution, cut_effs_evolution from .. evm.ic_containers import FitFunction from .. types.symbols import KrFitFunction @@ -284,3 +285,35 @@ def test_kr_time_evolution(dummy_kr_dst, dummy_kr_map): else: assert np.all(np.diff(evol_pars[col].values, axis=0) == 0), col +def test_cut_effs_evolution(dummy_kr_dst): + dst_1 = dummy_kr_dst[0] + dst_2 = dst_1.copy() + deltat = np.diff(dst_1.time)[0] + dst_2.time += deltat + max(dst_1.time) + dst = pd.concat([dst_1, dst_2], ignore_index=True) + t0 = min(dst.time) + t1 = max(dst.time) + ts, masks_time = get_time_series_df(ntimebins=2, time_range=(t0, t1), dst=dst) + + n = len(dst) // 2 + mask = np.ones(2*n, dtype = bool) + mask_s1 = mask .copy(); mask_s1 [:1*n//10] = False; mask_s1 [-1*n//10:] = False + mask_s2 = mask_s1.copy(); mask_s2 [:2*n//10] = False; mask_s2 [-2*n//10:] = False + mask_band = mask_s2.copy(); mask_band[:3*n//10] = False; mask_band[-3*n//10:] = False + + evol_table = pd.DataFrame({'ts': ts}) + + evol_table_updated = cut_effs_evolution(masks_time = masks_time, + dst = dst, + mask_s1 = mask_s1, + mask_s2 = mask_s2, + mask_band = mask_band, + evol_table = evol_table) + + assert all(col in evol_table_updated.columns for col in ['S1eff', 'S2eff', 'Bandeff']) + + for _, row in evol_table_updated.iterrows(): + + assert np.isclose(row.S1eff, 1-0.1) + assert np.isclose(row.S2eff, (1-2*0.1)/(1-1*0.1)) + assert np.isclose(row.Bandeff, (1-3*0.1)/(1-2*0.1))