diff --git a/.coveragerc b/.coveragerc new file mode 100644 index 0000000000..e33d0d9918 --- /dev/null +++ b/.coveragerc @@ -0,0 +1,11 @@ +# replace by exclude_also after bumping to >=7.2 +# with that, "pragma: no cover" can be removed +[report] +exclude_lines = + def __str__ + def __repr__ + pragma: no cover + +[run] +omit = + invisible_cities/database/* diff --git a/.github/workflows/coverage.yml b/.github/workflows/coverage.yml new file mode 100644 index 0000000000..eacfa9622d --- /dev/null +++ b/.github/workflows/coverage.yml @@ -0,0 +1,43 @@ +# This workflow will install IC and run all the tests with pytest + +name: Test suite + +on: + pull_request: + +jobs: + check-coverage: + runs-on: ubuntu-latest + + steps: + - uses: actions/checkout@v4 + - name: Get LFS files + run: git lfs pull + + - name: Install IC + run: | + source $CONDA/etc/profile.d/conda.sh + source manage.sh work_in_python_version_no_tests 3.8 + + - name: Produce coverage report + run: | + source $CONDA/etc/profile.d/conda.sh + source manage.sh work_in_python_version_no_tests 3.8 + coverage run -m pytest + coverage report -m > coverage_report + + - name: Check coverage report + run: | + echo "Showing only files with coverage < 100%" + head -2 coverage_report + incomplete_coverage=$(head -n -2 coverage_report | tail -n +3 | grep -v "100%") + + while read line; do + printf '\033[0;31m%s\033[0m\n' "$line" # red color ... reset color + done <<<${incomplete_coverage} + + tail -2 coverage_report + + if [ ! -z "$incomplete_coverage" ]; then + echo "::warning ::Test suite coverage incomplete" + fi diff --git a/.github/workflows/test_suite.yml b/.github/workflows/test_suite.yml index 48da250784..330b51f138 100644 --- a/.github/workflows/test_suite.yml +++ b/.github/workflows/test_suite.yml @@ -7,7 +7,7 @@ on: pull_request: jobs: - build: + build-and-test: strategy: fail-fast: false @@ -30,7 +30,6 @@ jobs: source $CONDA/etc/profile.d/conda.sh source manage.sh work_in_python_version_no_tests ${{ matrix.python-version }} - - name: Run tests run: | set -o pipefail # necessary for test failure to propagate to `tee` diff --git a/invisible_cities/calib/calib_functions_test.py b/invisible_cities/calib/calib_functions_test.py index a8e377f289..176a0eea39 100644 --- a/invisible_cities/calib/calib_functions_test.py +++ b/invisible_cities/calib/calib_functions_test.py @@ -15,6 +15,7 @@ from .. core import tbl_functions as tbl from .. core import fit_functions as fitf from .. core import system_of_units as units +from .. core.core_functions import in_range from .. core.stat_functions import poisson_sigma from .. evm.nh5 import SensorTable from .. types.symbols import SensorType @@ -148,9 +149,9 @@ def test_copy_sensor_table(config_tmpdir, sensor_type, sensors): @mark.parametrize('sensor_type , n_channel, gain_seed, gain_sigma_seed', ((SensorType.SIPM, 1, 16.5622, 2.5), (SensorType.PMT , 5, 24.9557, 9.55162))) -def test_seeds_db(sensor_type, n_channel, gain_seed, gain_sigma_seed, dbnew): +def test_seeds_db(sensor_type, n_channel, gain_seed, gain_sigma_seed): run_number = 6217 - result = cf.seeds_db(sensor_type, dbnew, run_number, n_channel) + result = cf.seeds_db(sensor_type, "new", run_number, n_channel) assert result == (gain_seed, gain_sigma_seed) @@ -190,14 +191,16 @@ def test_sensor_values(sensor_type, scaler, expected_range, min_b, max_b, half_w @mark.parametrize('sensor_type, run_number, n_chann, scaler', (( None, 6217, 1023, _dark_scaler_sipm), ( None, 6217, 0, _dark_scaler_pmt))) -def test_incorrect_sensor_type_raises_ValueError(sensor_type, dbnew, run_number, n_chann, scaler): +def test_incorrect_sensor_type_raises_ValueError(sensor_type, run_number, n_chann, scaler): bins = np.array([ -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7]) spec = np.array([ 28, 539, 1072, 1845, 2805, 3251, 3626, 3532, 3097, 2172, 1299, 665, 371, 174]) ped_vals = np.array([2.65181178e+04, 1.23743445e-01, 2.63794236e+00]) with raises(ValueError): - cf. seeds_db(sensor_type, dbnew, run_number, n_chann) + cf. seeds_db(sensor_type, "new", run_number, n_chann) + with raises(ValueError): cf.poisson_mu_seed(sensor_type, scaler, bins, spec, ped_vals) + with raises(ValueError): cf. sensor_values(sensor_type, scaler, bins, spec, ped_vals) @@ -214,6 +217,7 @@ def test_pedestal_values(): assert_approx_equal(ped_values.sigma_min, 0.001) +@mark.slow def test_compute_seeds_from_spectrum(ICDATADIR): PATH_IN = os.path.join(ICDATADIR, 'sipmcalspectra_R6358.h5') # Suppress warnings from division by zero in some bins. @@ -227,12 +231,10 @@ def test_compute_seeds_from_spectrum(ICDATADIR): for ich, (led, dar) in enumerate(zip(specsL, specsD)): b1 = 0 b2 = len(dar) - try: - valid_bins = np.argwhere(led>=min_stat) - b1 = valid_bins[ 0][0] - b2 = valid_bins[-1][0] - except IndexError: - continue + + valid_bins = np.argwhere(led>=min_stat) + b1 = valid_bins[ 0][0] + b2 = valid_bins[-1][0] peaks_dark = find_peaks_cwt(dar, np.arange(2, 20), min_snr=2) if len(peaks_dark) == 0: @@ -260,7 +262,53 @@ def test_compute_seeds_from_spectrum(ICDATADIR): assert gain_sigma_seed != 0 -def test_seeds_without_using_db(ICDATADIR, dbnew): +@mark.filterwarnings("ignore:Covariance of the parameters could not be estimated") +def test_seeds_and_bounds_with_db(): + values = np.array([1, 0.5, .1], dtype=float) + errors = (values + 1) / 20 + x = np.arange(-50, 300) + y = fitf.gauss(x, *values) + sel = in_range(x, -5, 5, right_closed=True) # magic + scaler_func = cf.dark_scaler(y[sel]) + seeds, bounds = cf.seeds_and_bounds(SensorType.SIPM, 8000, 0, + scaler_func, x, y, + values, "new", errors, + use_db_gain_seeds=True) + + + # this test doesn't do much, to be improved + assert all(seeds) + assert len(seeds) == 4 + assert len(bounds) == 2 + assert len(bounds[0]) == 4 + assert len(bounds[1]) == 4 + + +@mark.filterwarnings("ignore:Covariance of the parameters could not be estimated") +def test_seeds_and_bounds_with_gau(): + values = np.array([1, 0.5, .1], dtype=float) + errors = (values + 1) / 20 + x = np.arange(-50, 300) + y = fitf.gauss(x, *values) + sel = in_range(x, -5, 5, right_closed=True) # magic + scaler_func = cf.dark_scaler(y[sel]) + seeds, bounds = cf.seeds_and_bounds(SensorType.SIPM, 8000, 0, + scaler_func, x, y, + values, "new", errors, + "gau", + use_db_gain_seeds=True) + + + # this test doesn't do much, to be improved + assert all(seeds) + assert len(seeds) == 6 + assert len(bounds) == 2 + assert len(bounds[0]) == 6 + assert len(bounds[1]) == 6 + + +@mark.slow +def test_seeds_without_using_db(ICDATADIR): PATH_IN = os.path.join(ICDATADIR, 'sipmcalspectra_R6358.h5') # Suppress warnings from division by zero in some bins. with warnings.catch_warnings(), tb.open_file(PATH_IN) as h5in: @@ -275,12 +323,10 @@ def test_seeds_without_using_db(ICDATADIR, dbnew): for ich, (led, dar) in enumerate(zip(specsL, specsD)): b1 = 0 b2 = len(dar) - try: - valid_bins = np.argwhere(led>=min_stat) - b1 = valid_bins[ 0][0] - b2 = valid_bins[-1][0] - except IndexError: - continue + + valid_bins = np.argwhere(led>=min_stat) + b1 = valid_bins[ 0][0] + b2 = valid_bins[-1][0] peaks_dark = find_peaks_cwt(dar, np.arange(2, 20), min_snr=2) if len(peaks_dark) == 0: @@ -300,8 +346,23 @@ def test_seeds_without_using_db(ICDATADIR, dbnew): scaler_func = cf.dark_scaler(dar[p_range][p_bins]) seeds, bounds = cf.seeds_and_bounds(SensorType.SIPM, run_no, ich, scaler_func, bins[p_range], - led[p_range], ped_vals, dbnew, + led[p_range], ped_vals, "new", gfitRes.errors, use_db_gain_seeds=False) assert all(seeds) assert bounds == ((0, 0, 0, 0.001), (np.inf, 10000, 10000, 10000)) + + +@mark.parametrize( "f nargs".split() + , ( (cf. seeds_db, 3) + , (cf.poisson_mu_seed, 4) + , (cf. sensor_values, 4))) +def test_sensor_type_input(f, nargs): + """ + Test that if the first argument (sensor_type) is not a symbol, the + functions complain. The rest of the arguments are irrelevant, so + they are just dummy values. + """ + dummy = [None]*nargs + with raises(ValueError): f( "PMT", *dummy) + with raises(ValueError): f("SiPM", *dummy) diff --git a/invisible_cities/calib/calib_sensors_functions.py b/invisible_cities/calib/calib_sensors_functions.py index f8db736c3b..8ae9ea7dbc 100644 --- a/invisible_cities/calib/calib_sensors_functions.py +++ b/invisible_cities/calib/calib_sensors_functions.py @@ -10,10 +10,6 @@ from ..types.symbols import SiPMCalibMode -def mask_sensors(wfs, active): - return wfs * active.astype(wfs.dtype).reshape(active.size, 1) - - def scipy_mode(x, axis=0): """ Scipy implementation of the mode (runs very slow). diff --git a/invisible_cities/cities/beersheba.py b/invisible_cities/cities/beersheba.py index c30cd50048..412862d91d 100644 --- a/invisible_cities/cities/beersheba.py +++ b/invisible_cities/cities/beersheba.py @@ -187,15 +187,6 @@ def deconvolve_signal(det_db : pd.DataFrame, **satellite_params, inter_method = inter_method) - if not isinstance(energy_type , HitEnergy ): - raise ValueError(f'energy_type {energy_type} is not a valid energy type.') - if not isinstance(inter_method, InterpolationMethod): - raise ValueError(f'inter_method {inter_method} is not a valid interpolation method.') - if not isinstance(cut_type , CutType ): - raise ValueError(f'cut_type {cut_type} is not a valid cut type.') - if not isinstance(deconv_mode , DeconvolutionMode ): - raise ValueError(f'deconv_mode {deconv_mode} is not a valid deconvolution mode.') - def deconvolve_hits(df, z): ''' Given an slice, applies deconvolution using the PSF @@ -259,7 +250,14 @@ def apply_deconvolution(df): return apply_deconvolution -def create_deconvolution_df(hits, deconv_e, pos, cut_type, e_cut, n_dim): +@check_annotations +def create_deconvolution_df( hits : pd.DataFrame + , deconv_e : np.ndarray + , pos : Tuple[np.ndarray, np.ndarray, np.ndarray] + , cut_type : CutType + , e_cut : float + , n_dim : int + ) -> pd.DataFrame: ''' Given the output of the deconvolution, it cuts the low energy voxels and creates a dataframe object with the resulting output. @@ -280,12 +278,8 @@ def create_deconvolution_df(hits, deconv_e, pos, cut_type, e_cut, n_dim): df = pd.DataFrame(columns=['event', 'npeak', 'X', 'Y', 'Z', 'E']) - if cut_type is CutType.abs: - sel_deconv = deconv_e > e_cut - elif cut_type is CutType.rel: - sel_deconv = deconv_e / deconv_e.max() > e_cut - else: - raise ValueError(f'cut_type {cut_type} is not a valid cut type.') + if cut_type is CutType.abs: sel_deconv = deconv_e > e_cut + else : sel_deconv = deconv_e / deconv_e.max() > e_cut df['E'] = deconv_e[sel_deconv] df['event'] = hits.event.unique()[0] @@ -497,9 +491,10 @@ def beersheba( files_in : OneOrManyFiles for p in ['sample_width', 'bin_size', 'diffusion']: if len(deconv_params[p]) != deconv_params['n_dim']: - raise ValueError (f"Parameter {p} dimensions do not match n_dim parameter") + raise ValueError(f"Parameter {p} dimensions do not match n_dim parameter") + if deconv_params['n_dim'] > 2: - raise NotImplementedError(f"{deconv_params['n_dim']}-dimensional PSF not yet implemented") + raise NotImplementedError(f"{deconv_params['n_dim']}-dimensional PSF not yet implemented") cut_sensors = fl.map(cut_over_Q (deconv_params.pop("q_cut") , ['E', 'Ec']), item = 'hits') diff --git a/invisible_cities/cities/beersheba_test.py b/invisible_cities/cities/beersheba_test.py index a843790b52..f71d13b20f 100644 --- a/invisible_cities/cities/beersheba_test.py +++ b/invisible_cities/cities/beersheba_test.py @@ -20,28 +20,22 @@ from .. types.symbols import CutType -def test_create_deconvolution_df(ICDATADIR): +@mark.parametrize("cut_type", CutType) +def test_create_deconvolution_df(ICDATADIR, cut_type): true_in = os.path.join(ICDATADIR, "exact_Kr_deconvolution_with_MC.h5") true_dst = dio.load_dst(true_in, 'DECO', 'Events') ecut = 1e-2 new_dst = pd.concat([create_deconvolution_df(t, t.E.values, (t.X.values, t.Y.values, t.Z.values), - CutType.abs, ecut, 3) for _, t in true_dst.groupby('event')]) - true_dst = true_dst.loc[true_dst.E > ecut, :].reset_index(drop=True) - # compare only existing columns - true_dst = true_dst.loc[:, new_dst.columns.values.tolist()] - assert_dataframes_close(new_dst .reset_index(drop=True), true_dst.reset_index(drop=True)) + cut_type, ecut, 3) for _, t in true_dst.groupby('event')]) + for event, df in true_dst.groupby("event"): + sel = (df.E > ecut if cut_type is CutType.abs else + df.E > ecut * df.E.max()) + df = df.loc[sel].reset_index(drop=True) + # compare only existing columns + df = df.loc[:, new_dst.columns.values.tolist()] - -@mark.parametrize("cut_type", CutType.__members__) -def test_create_deconvolution_df_cuttype(ICDATADIR, cut_type): - true_in = os.path.join(ICDATADIR, "exact_Kr_deconvolution_with_MC.h5") - true_dst = dio.load_dst(true_in, 'DECO', 'Events') - ecut = 1e-2 - - with raises(ValueError): - create_deconvolution_df(true_dst, true_dst.E.values, - (true_dst.X.values, true_dst.Y.values, true_dst.Z.values), - cut_type, ecut, 3) + new_event = new_dst.loc[new_dst.event==event] + assert_dataframes_close(new_event.reset_index(drop=True), df.reset_index(drop=True)) def test_distribute_energy(ICDATADIR): @@ -56,6 +50,18 @@ def test_distribute_energy(ICDATADIR): assert np.isclose (true_dst1.E.sum(), true_dst2.E.sum()) +def test_beersheba_ndim(beersheba_config, config_tmpdir): + path_out = os.path.join(config_tmpdir, "beersheba_ndim.h5") + beersheba_config.update(dict(file_out = path_out)) + beersheba_config["deconv_params"].update(dict( n_dim = 3 + , sample_width = [1]*3 + , bin_size = [1]*3 + , diffusion = [1]*3 )) + + with raises(NotImplementedError): + beersheba(**beersheba_config) + + @ignore_warning.no_config_group @ignore_warning.str_length @ignore_warning.not_kdst @@ -107,6 +113,23 @@ def test_beersheba_exact_result( deco expected = getattr(true_output_file.root, table) assert_tables_equality(got, expected, rtol=1e-6) +@mark.parametrize( "updates" + , ( dict(sample_width = [1]), dict(sample_width = [1,1,1]) + , dict(bin_size = [1]), dict(bin_size = [1,1,1]) + , dict(diffusion = [1]), dict(diffusion = [1,1,1]) + )) +def test_beersheba_checks_dimensions(config_tmpdir, beersheba_config, updates): + # check that if these quantities do not have the correct number of + # dimensions (2), the city raises an error + par = next(iter(updates)) + n = len(updates[par]) + path_out = os.path.join(config_tmpdir, f"beersheba_only_check_dims_{par}_{n}.h5") + beersheba_config.update(dict(file_out = path_out)) + beersheba_config['deconv_params'].update(updates) + + with raises(ValueError, match="Parameter .* dimensions do not match n_dim parameter"): + beersheba(**beersheba_config) + @ignore_warning.no_config_group @ignore_warning.str_length diff --git a/invisible_cities/cities/cities_test.py b/invisible_cities/cities/cities_test.py index 0e66874cb3..b86aaccc11 100644 --- a/invisible_cities/cities/cities_test.py +++ b/invisible_cities/cities/cities_test.py @@ -22,6 +22,7 @@ @ignore_warning.no_config_group @ignore_warning.no_mc_tables @ignore_warning.not_kdst +@ignore_warning.no_pmaps @ignore_warning.str_length @mark.parametrize("city", online_cities) def test_city_empty_input_file(config_tmpdir, ICDATADIR, city): @@ -118,6 +119,7 @@ def test_city_output_contains_configuration(config_tmpdir, city): @ignore_warning.no_config_group +@ignore_warning.no_pmaps def test_cities_carry_on_configuration_from_previous_ones(ICDATADIR, config_tmpdir): rwf_file = os.path.join( ICDATADIR, "electrons_40keV_ACTIVE_10evts_RWF.h5") pmaps_file = os.path.join(config_tmpdir, f"test_cities_carry_on_configuration_pmaps.h5") diff --git a/invisible_cities/cities/command_line_execution_test.py b/invisible_cities/cities/command_line_execution_test.py index 9476e8f623..02419bafa9 100644 --- a/invisible_cities/cities/command_line_execution_test.py +++ b/invisible_cities/cities/command_line_execution_test.py @@ -26,7 +26,7 @@ def test_command_line_run(city, tmpdir_factory): .format(**locals())) try: check_output(command, shell = True, stderr=STDOUT) - except CalledProcessError as e: + except CalledProcessError as e: # pragma: no cover # Ensure that stdout and stderr are visible when test fails print(e.stdout.decode()) raise diff --git a/invisible_cities/cities/components.py b/invisible_cities/cities/components.py index 1bb886c93d..9557c1ee20 100644 --- a/invisible_cities/cities/components.py +++ b/invisible_cities/cities/components.py @@ -261,7 +261,7 @@ def print_every(N): fl.sink (lambda data: print(f"events processed: {data['index']}, event number: {data['event_number']}"))) -def print_every_alternative_implementation(N): +def print_every_alternative_implementation(N): # pragma: no cover @fl.coroutine def print_every_loop(target): with fl.closing(target): @@ -273,7 +273,11 @@ def print_every_loop(target): return print_every_loop -def get_actual_sipm_thr(thr_sipm_type, thr_sipm, detector_db, run_number): +@check_annotations +def get_actual_sipm_thr( thr_sipm_type : SiPMThreshold + , thr_sipm : float + , detector_db : str + , run_number : int): if thr_sipm_type is SiPMThreshold.common: # In this case, the threshold is a value in pes sipm_thr = thr_sipm @@ -283,10 +287,6 @@ def get_actual_sipm_thr(thr_sipm_type, thr_sipm, detector_db, run_number): noise_sampler = NoiseSampler(detector_db, run_number) sipm_thr = noise_sampler.compute_thresholds(thr_sipm) - else: - raise ValueError(f"Unrecognized thr type: {thr_sipm_type}. " - "Only valid options are `common` and `individual`") - return sipm_thr @@ -407,21 +407,27 @@ def deconv_pmt(RWF): def get_run_number(h5in): - if "runInfo" in h5in.root.Run: return h5in.root.Run.runInfo[0]['run_number'] - elif "RunInfo" in h5in.root.Run: return h5in.root.Run.RunInfo[0]['run_number'] + group = "Run" in h5in.root + table = None + if group and "runInfo" in h5in.root.Run: table = h5in.root.Run.runInfo + elif group and "RunInfo" in h5in.root.Run: table = h5in.root.Run.RunInfo + + if table is not None and len(table): + return table[0]['run_number'] raise tb.exceptions.NoSuchNodeError(f"No node runInfo or RunInfo in file {h5in}") -def get_pmt_wfs(h5in, wf_type): - if wf_type is WfType.rwf : return h5in.root.RD.pmtrwf - elif wf_type is WfType.mcrd: return h5in.root. pmtrd - else : raise TypeError(f"Invalid WfType: {type(wf_type)}") +@check_annotations +def get_pmt_wfs(h5in : tb.File, wf_type : WfType): + if wf_type is WfType.rwf: return h5in.root.RD.pmtrwf + else : return h5in.root. pmtrd -def get_sipm_wfs(h5in, wf_type): - if wf_type is WfType.rwf : return h5in.root.RD.sipmrwf - elif wf_type is WfType.mcrd: return h5in.root. sipmrd - else : raise TypeError(f"Invalid WfType: {type(wf_type)}") + +@check_annotations +def get_sipm_wfs(h5in : tb.File, wf_type : WfType): + if wf_type is WfType.rwf: return h5in.root.RD.sipmrwf + else : return h5in.root. sipmrd def get_trigger_info(h5in): @@ -553,18 +559,13 @@ def wf_from_files(paths, wf_type): def pmap_from_files(paths): for path in paths: - try: - pmaps = load_pmaps(path, lazy=True) - except tb.exceptions.NoSuchNodeError: - continue - with tb.open_file(path, "r") as h5in: try: - run_number = get_run_number(h5in) - event_info = get_event_info(h5in) + pmaps = load_pmaps(path, lazy=True) + run_number = get_run_number(h5in) + event_info = get_event_info(h5in) except tb.exceptions.NoSuchNodeError: - continue - except IndexError: + warnings.warn("No PMAPs in input file", UserWarning) continue for evtinfo, (evt, pmap) in zip(event_info, pmaps): @@ -589,11 +590,8 @@ def hits_and_kdst_from_files( paths : List[str] continue with tb.open_file(path, "r") as h5in: - try: - run_number = get_run_number(h5in) - event_info = get_event_info(h5in) - except (tb.exceptions.NoSuchNodeError, IndexError): - continue + run_number = get_run_number(h5in) + event_info = get_event_info(h5in) check_lengths(event_info, kdst_df.event.unique()) @@ -620,13 +618,18 @@ def dst_from_files(paths: List[str], group: str, node:str) -> Iterator[Dict[str, for path in paths: try: df = load_dst(path, group, node) - with tb.open_file(path, "r") as h5in: - run_number = get_run_number(h5in) - evt_numbers = get_event_info(h5in).col("evt_number") - timestamps = get_event_info(h5in).col("timestamp") except (tb.exceptions.NoSuchNodeError, IndexError): continue + if not len(df): + warnings.warn(f"No data in node /{group}/{node} in input file", UserWarning) + continue + + with tb.open_file(path, "r") as h5in: + run_number = get_run_number(h5in) + evt_numbers = get_event_info(h5in).col("evt_number") + timestamps = get_event_info(h5in).col("timestamp") + yield dict( dst = df , run_number = run_number , event_numbers = evt_numbers @@ -638,10 +641,7 @@ def dst_from_files(paths: List[str], group: str, node:str) -> Iterator[Dict[str, def MC_hits_from_files(files_in : List[str], rate: float) -> Generator: timestamp = create_timestamp(rate) for filename in files_in: - try: - hits_df = load_mchits_df(filename) - except tb.exceptions.NoSuchNodeError: - continue + hits_df = load_mchits_df(filename) l_type = hits_df.dtypes['label'] map_df = load_mcstringmap(filename) if l_type == np.int32 else None @@ -669,15 +669,15 @@ def dhits_from_files(paths: List[str]) -> Iterator[Dict[str,Union[HitCollection, except tb.exceptions.NoSuchNodeError: continue + if not len(dhits_df): + warnings.warn(f"No data in node /DECO/Events in input file", UserWarning) + continue + kdst_df = load_dst(path, 'DST', 'Events', ignore_errors=True) with tb.open_file(path, "r") as h5in: - try: - run_number = get_run_number(h5in) - except (tb.exceptions.NoSuchNodeError, IndexError): - continue - event_info = load_dst(path, 'Run', 'events') + run_number = get_run_number(h5in) for evt in dhits_df.event.unique(): this_event = lambda df: df.event==evt timestamp = event_info[event_info.evt_number==evt].timestamp.values[0] @@ -696,11 +696,11 @@ def dhits_from_files(paths: List[str]) -> Iterator[Dict[str,Union[HitCollection, timestamp = timestamp ) -def sensor_data(path, wf_type): +@check_annotations +def sensor_data(path : str, wf_type : WfType): with tb.open_file(path, "r") as h5in: - if wf_type is WfType.rwf : (pmt_wfs, sipm_wfs) = (h5in.root.RD .pmtrwf, h5in.root.RD .sipmrwf) - elif wf_type is WfType.mcrd: (pmt_wfs, sipm_wfs) = (h5in.root. pmtrd , h5in.root. sipmrd ) - else : raise TypeError(f"Invalid WfType: {type(wf_type)}") + if wf_type is WfType.rwf: (pmt_wfs, sipm_wfs) = (h5in.root.RD .pmtrwf, h5in.root.RD .sipmrwf) + else : (pmt_wfs, sipm_wfs) = (h5in.root. pmtrd , h5in.root. sipmrd ) _, NPMT , PMTWL = pmt_wfs.shape _, NSIPM, SIPMWL = sipm_wfs.shape return SensorData(NPMT=NPMT, PMTWL=PMTWL, NSIPM=NSIPM, SIPMWL=SIPMWL) @@ -761,21 +761,6 @@ def calibrate_sipms(rwf): return calibrate_sipms -def calibrate_with_mean(dbfile, run_number): - DataSiPM = load_db.DataSiPM(dbfile, run_number) - adc_to_pes = np.abs(DataSiPM.adc_to_pes.values) - def calibrate_with_mean(wfs): - return csf.subtract_baseline_and_calibrate(wfs, adc_to_pes) - return calibrate_with_mean - -def calibrate_with_maw(dbfile, run_number, n_maw_sipm): - DataSiPM = load_db.DataSiPM(dbfile, run_number) - adc_to_pes = np.abs(DataSiPM.adc_to_pes.values) - def calibrate_with_maw(wfs): - return csf.subtract_baseline_maw_and_calibrate(wfs, adc_to_pes, n_maw_sipm) - return calibrate_with_maw - - def zero_suppress_wfs(thr_csum_s1, thr_csum_s2): def ccwfs_to_zs(ccwf_sum, ccwf_sum_maw): return (pkf.indices_and_wf_above_threshold(ccwf_sum_maw, thr_csum_s1).indices, @@ -1212,7 +1197,6 @@ def copy_Efield(hitc : HitCollection) -> HitCollection: def make_event_summary(event_number : int , topology_info : pd.DataFrame , paolina_hits : HitCollection, - out_of_map : bool ) -> pd.DataFrame: """ For a given event number, timestamp, topology info dataframe, paolina hits and kdst information returns a @@ -1256,7 +1240,7 @@ def make_event_summary(event_number : int , list_of_vars = [event_number, S2ec, S2qc, ntrks, nhits, *ave_pos, ave_r, min(x), min(y), min(z), min(r), max(x), max(y), max(z), max(r), - out_of_map] + False] # out_of_map es.loc[0] = list_of_vars #change dtype of columns to match type of variables @@ -1328,17 +1312,9 @@ def create_extract_track_blob_info(hitc): if len(hitc.hits) > max_num_hits: return df, hitc, True #track_hits is a new Hitcollection object that contains hits belonging to tracks, and hits that couldnt be corrected - track_hitc = HitCollection(hitc.event, hitc.time) - out_of_map = np.any(np.isnan([h.Ep for h in hitc.hits])) - if out_of_map: - #add nan hits to track_hits, the track_id will be -1 - track_hitc.hits.extend ([h for h in hitc.hits if np.isnan (h.Ep)]) - hits_without_nan = [h for h in hitc.hits if np.isfinite(h.Ep)] - #create new Hitcollection object but keep the name hitc - hitc = HitCollection(hitc.event, hitc.time) - hitc.hits = hits_without_nan - + track_hitc = HitCollection(hitc.event, hitc.time) hit_energies = np.array([getattr(h, HitEnergy.Ep.value) for h in hitc.hits]) + assert not np.isnan(hit_energies).any() if len(hitc.hits) > 0 and (hit_energies>0).any(): voxels = plf.voxelize_hits(hitc.hits, vox_size, strict_vox_size, HitEnergy.Ep) @@ -1396,7 +1372,7 @@ def create_extract_track_blob_info(hitc): #change dtype of columns to match type of variables df = df.apply(lambda x : x.astype(types_dict_tracks[x.name])) track_hitc.hits.extend(track_hits) - return df, track_hitc, out_of_map + return df, track_hitc return create_extract_track_blob_info @@ -1424,7 +1400,7 @@ def compute_and_write_tracks_info(paolina_params, h5out, # Create tracks and compute topology-related information create_extract_track_blob_info = fl.map(track_blob_info_creator_extractor(**paolina_params), args = 'Ep_hits', - out = ('topology_info', 'paolina_hits', 'out_of_map')) + out = ('topology_info', 'paolina_hits')) sort_hits_ = fl.map(sort_hits, item="paolina_hits") @@ -1436,7 +1412,7 @@ def compute_and_write_tracks_info(paolina_params, h5out, # Create table with summary information make_final_summary = fl.map(make_event_summary, - args = ('event_number', 'topology_info', 'paolina_hits', 'out_of_map'), + args = ('event_number', 'topology_info', 'paolina_hits'), out = 'event_info') diff --git a/invisible_cities/cities/components_test.py b/invisible_cities/cities/components_test.py index 7fb441fb8e..7e7ed95c00 100644 --- a/invisible_cities/cities/components_test.py +++ b/invisible_cities/cities/components_test.py @@ -17,32 +17,41 @@ from .. core.exceptions import InvalidInputFileStructure from .. core.exceptions import SensorIDMismatch from .. core.exceptions import NoInputFiles +from .. core.exceptions import MCEventNotFound from .. core.testing_utils import assert_tables_equality from .. core.testing_utils import ignore_warning from .. core import system_of_units as units from .. evm.event_model import Cluster from .. evm.event_model import Hit from .. evm.event_model import HitCollection +from .. reco.xy_algorithms import barycenter from .. types.ic_types import xy from .. types.symbols import WfType from .. types.symbols import EventRange as ER from .. types.symbols import NormStrategy from .. types.symbols import XYReco +from .. types.symbols import SiPMThreshold from . components import event_range +from . components import get_actual_sipm_thr from . components import collect from . components import copy_mc_info from . components import wf_from_files from . components import pmap_from_files from . components import compute_xy_position from . components import city +from . components import dst_from_files from . components import hits_and_kdst_from_files +from . components import dhits_from_files from . components import mcsensors_from_file from . components import create_timestamp from . components import check_max_time from . components import hits_corrector from . components import write_city_configuration from . components import copy_cities_configuration +from . components import try_global_reco +from . components import length_of +from . components import get_run_number from .. dataflow import dataflow as fl @@ -206,8 +215,7 @@ def test_city_fails_if_bad_input_file(config_tmpdir, ICDATADIR): @city def dummy_city( files_in : Union[str, list] , file_out : str - , event_range : tuple): - pass + , event_range : tuple): pass with raises(FileNotFoundError): dummy_city(files_in=files_in, file_out=file_out, event_range=(0, 1)) @@ -249,7 +257,7 @@ def test_compute_xy_position_depends_on_actual_run_number(): @mark.skip(reason="there should not be a default detector db") -def test_city_adds_default_detector_db(config_tmpdir): +def test_city_adds_default_detector_db(config_tmpdir): # pragma: no cover default_detector_db = 'new' args = {'files_in' : 'dummy_in', 'file_out' : os.path.join(config_tmpdir, 'dummy_out')} @@ -264,7 +272,7 @@ def dummy_city(files_in, file_out, event_range, detector_db): @mark.skip(reason="there should not be a default detector db") -def test_city_does_not_overwrite_detector_db(config_tmpdir): +def test_city_does_not_overwrite_detector_db(config_tmpdir): # pragma: no cover args = {'detector_db' : 'some_detector', 'files_in' : 'dummy_in', 'file_out' : os.path.join(config_tmpdir, 'dummy_out')} @@ -279,7 +287,7 @@ def dummy_city(files_in, file_out, event_range, detector_db): @mark.skip(reason="there should not be a default detector db") -def test_city_only_pass_default_detector_db_when_expected(config_tmpdir): +def test_city_only_pass_default_detector_db_when_expected(config_tmpdir): # pragma: no cover args = {'files_in' : 'dummy_in', 'file_out' : os.path.join(config_tmpdir, 'dummy_out')} @city @@ -289,6 +297,7 @@ def dummy_city(files_in, file_out, event_range): dummy_city(**args) + def test_pmap_from_files_event_number_mismatch_raises(KrMC_pmaps_filename, output_tmpdir): filename = os.path.join(output_tmpdir, "test_pmap_from_files_event_number_mismatch_raises.h5") @@ -303,6 +312,12 @@ def test_pmap_from_files_event_number_mismatch_raises(KrMC_pmaps_filename, outpu next(generator) +def test_dst_from_files_empty_input(ICDATADIR): + empty_file = os.path.join(ICDATADIR, "empty_kdst.h5") + with warns(UserWarning, match="No data in node /DST/Events in input file"): + tuple(dst_from_files([empty_file], "DST", "Events")) # consume iterator + + def test_hits_and_kdst_from_files(ICDATADIR): event_number = 1 timestamp = 0. @@ -327,6 +342,27 @@ def test_hits_and_kdst_from_files_missing_hits(Th228_hits_missing, config_tmpdir assert n_events == n_events_true +def test_dhits_from_files_empty_input(ICDATADIR): + empty_file = os.path.join(ICDATADIR, "empty_ddst.h5") + with warns(UserWarning, match="No data in node /DECO/Events in input file"): + tuple(dhits_from_files([empty_file])) # consume iterator + + +@mark.parametrize("thr", (0, 1.5, 3)) +def test_get_actual_sipm_threshold_common(thr): + got = get_actual_sipm_thr(SiPMThreshold.common, thr, "new", 0) + assert np.isclose(got, thr) + + +def test_get_actual_sipm_threshold_individual(): + """ + 99.9% of the noise pds should be above 1 pe + """ + thr = 0.999 + got = get_actual_sipm_thr(SiPMThreshold.individual, thr, "new", 8000) + assert np.all(got > 1) + + def test_collect(): the_source = list(range(0,10)) the_collector = collect() @@ -345,7 +381,7 @@ def test_copy_mc_info_noMC(ICDATADIR, config_tmpdir): @mark.xfail -def test_copy_mc_info_repeated_event_numbers(ICDATADIR, config_tmpdir): +def test_copy_mc_info_repeated_event_numbers(ICDATADIR, config_tmpdir): # pragma: no cover file_in = os.path.join(ICDATADIR, "Kr83_nexus_v5_03_00_ACTIVE_7bar_10evts.sim.h5") file_out = os.path.join(config_tmpdir, "dummy_out.h5") @@ -373,6 +409,15 @@ def test_copy_mc_info_split_nexus_events(ICDATADIR, config_tmpdir): assert_tables_equality(got, expected) +def test_copy_mc_info_raises_missing_event(ICDATADIR, config_tmpdir): + file_in = os.path.join(ICDATADIR, "Kr83_nexus_v5_03_00_ACTIVE_7bar_10evts.sim.h5") + file_out = os.path.join(config_tmpdir, "copy_mc_info_raises_missing_event.h5") + + event_that_does_not_exist = 99999 + with tb.open_file(file_out, 'w') as h5out: + with raises(MCEventNotFound, match="Some events not found in MC tables"): + copy_mc_info([file_in], h5out, [event_that_does_not_exist], 'new', 0) + def test_mcsensors_from_file_fast_returns_empty(ICDATADIR): rate = 0.5 @@ -475,6 +520,7 @@ def test_check_max_time_units(): with raises(ValueError): check_max_time(max_time, buffer_length) + def test_read_wrong_pmt_ids(ICDATADIR): """ The input file of this test contains sensor IDs that are not present in the database. @@ -494,12 +540,12 @@ def test_hits_Z_uncorrected( correction_map_filename ''' Test to ensure that z is uncorrected when `apply_z` is False ''' - + hc = random_hits_toy_data hz = [h.Z for h in hc.hits] - + correct = hits_corrector(correction_map_filename, - apply_temp = False, + apply_temp = False, norm_strat = NormStrategy.kr, norm_value = None, apply_z = False) @@ -514,12 +560,12 @@ def test_hits_Z_corrected_when_flagged( correction_map_filename ''' Test to ensure that the correction is applied when `apply_z` is True ''' - + hc = random_hits_toy_data hz = [h.Z for h in hc.hits] - + correct = hits_corrector(correction_map_filename, - apply_temp = False, + apply_temp = False, norm_strat = NormStrategy.kr, norm_value = None, apply_z = True) @@ -527,7 +573,7 @@ def test_hits_Z_corrected_when_flagged( correction_map_filename # raise assertion error as expected assert_raises(AssertionError, assert_equal, corrected_z, hz) - + @mark.parametrize( "norm_strat norm_value".split(), ( (NormStrategy.kr , None) # None marks the default value @@ -574,6 +620,35 @@ def test_hits_corrector_invalid_normalization_options_raises( correction_map_fil hits_corrector(correction_map_filename, False, norm_strat, norm_value) +@mark.parametrize( "pos qs thr".split() + , ( ([ ], [ ], 0) # no hits + , ([[1,1]], [1], 2) # no hits above threshold + )) +def test_try_global_reco_creates_empty_cluster(pos, qs, thr): + pos, qs = np.array(pos), np.array(qs) + reco = lambda *args: barycenter(*args, Qthr=thr) + out = try_global_reco(reco, pos, qs) + assert out.x == xy.empty().x + assert out.y == xy.empty().y + + +@mark.parametrize("t", (int, float, object)) +def test_length_of_incorrect_type(t): + with raises(TypeError, match="Cannot determine size of .*"): + length_of(t()) + + +def test_get_run_number_raises(ICDATADIR, config_tmpdir): + filename = os.path.join(config_tmpdir, "test_get_run_number_raises.h5") + with tb.open_file(filename, "a") as file: + with raises(tb.exceptions.NoSuchNodeError, match="No node runInfo or RunInfo .*"): + get_run_number(file) + + file.create_group(file.root, "Run") + with raises(tb.exceptions.NoSuchNodeError, match="No node runInfo or RunInfo .*"): + get_run_number(file) + + def test_write_city_configuration(config_tmpdir): filename = os.path.join(config_tmpdir, "test_write_configuration.h5") city_name = "acity" @@ -597,31 +672,58 @@ def test_write_city_configuration(config_tmpdir): def test_copy_cities_configuration(config_tmpdir): - filename1 = os.path.join(config_tmpdir, "test_copy_cities_configuration_1.h5") - filename2 = os.path.join(config_tmpdir, "test_copy_cities_configuration_2.h5") + filename12 = os.path.join(config_tmpdir, "test_copy_cities_configuration_12.h5") + filename3 = os.path.join(config_tmpdir, "test_copy_cities_configuration_3.h5") city_name1 = "acity" city_name2 = "bcity" + city_name3 = "ccity" args = dict( a = 1, b = 2.3, c = "a_string", ) - write_city_configuration(filename1, city_name1, args) - write_city_configuration(filename2, city_name2, args) + write_city_configuration(filename12, city_name1, args) + write_city_configuration(filename12, city_name2, args) + write_city_configuration(filename3 , city_name3, args) - copy_cities_configuration(filename1, filename2) - with tb.open_file(filename2, "r") as file: + copy_cities_configuration(filename12, filename3) + with tb.open_file(filename3, "r") as file: assert "config" in file.root assert city_name1 in file.root.config assert city_name2 in file.root.config + assert city_name3 in file.root.config - df1 = pd.read_hdf(filename1, "/config/" + city_name1).set_index("variable") - df2 = pd.read_hdf(filename2, "/config/" + city_name2).set_index("variable") + df1 = pd.read_hdf(filename3, f"/config/{city_name1}").set_index("variable") + df2 = pd.read_hdf(filename3, f"/config/{city_name2}").set_index("variable") + df3 = pd.read_hdf(filename3, f"/config/{city_name3}").set_index("variable") for var, value in args.items(): assert var in df1.index assert var in df2.index + assert var in df3.index assert str(value) == df1.value.loc[var] assert str(value) == df2.value.loc[var] + assert str(value) == df3.value.loc[var] + + +def test_copy_cities_configuration_without_existing_node(config_tmpdir): + filename1 = os.path.join(config_tmpdir, "test_copy_cities_configuration_without_existing_node_in.h5") + filename2 = os.path.join(config_tmpdir, "test_copy_cities_configuration_without_existing_node_out.h5") + city_name = "acity" + args = dict( + a = 1, + b = 2.3, + c = "a_string", + ) + write_city_configuration(filename1, city_name, args) + copy_cities_configuration(filename1, filename2) + with tb.open_file(filename2, "r") as file: + assert "config" in file.root + assert city_name in file.root.config + + df = pd.read_hdf(filename1, f"/config/{city_name}").set_index("variable") + for var, value in args.items(): + assert var in df.index + assert str(value) == df.value.loc[var] def test_copy_cities_configuration_warns_when_nothing_to_copy(ICDATADIR, config_tmpdir): diff --git a/invisible_cities/cities/diomira_test.py b/invisible_cities/cities/diomira_test.py index eee76cebeb..221c062043 100644 --- a/invisible_cities/cities/diomira_test.py +++ b/invisible_cities/cities/diomira_test.py @@ -23,6 +23,8 @@ from .. types.symbols import all_events from . diomira import diomira +from . diomira import select_trigger_filter +from . diomira import emulate_trigger def test_diomira_identify_bug(ICDATADIR): @@ -215,7 +217,7 @@ def test_diomira_can_fix_random_seed(output_tmpdir): ## to run the following test, use the --runslow option with pytest @mark.veryslow -def test_diomira_reproduces_singlepe(ICDATADIR, output_tmpdir): +def test_diomira_reproduces_singlepe(ICDATADIR, output_tmpdir): # pragma: no cover file_in = os.path.join(ICDATADIR , 'single_pe_pmts.h5') file_out = os.path.join(output_tmpdir, 'single_pe_elec_sim.h5') @@ -246,3 +248,15 @@ def test_diomira_reproduces_singlepe(ICDATADIR, output_tmpdir): ## show the mean is within error of reproducing 1 pe ADC for fit_res in fits: assert np.abs(fit_res.values[1]) < fit_res.errors[1] + + +def test_select_trigger_filter_invalid_trigger_type(): + ftype = "DOES_NOT_EXIST" + with raises(ValueError, match=f"Invalid trigger type: '{ftype}'"): + select_trigger_filter(ftype, dict(), dict()) + + +def test_emulate_trigger_invalid_trigger_type(): + ftype = "DOES_NOT_EXIST" + with raises(ValueError, match=f"Invalid trigger type: '{ftype}'"): + emulate_trigger("new", 0, ftype, dict(), dict()) diff --git a/invisible_cities/cities/dorothea_test.py b/invisible_cities/cities/dorothea_test.py index 8b86469d5d..b92cb4427a 100644 --- a/invisible_cities/cities/dorothea_test.py +++ b/invisible_cities/cities/dorothea_test.py @@ -3,7 +3,9 @@ import numpy as np import tables as tb -from pytest import mark +from numpy.testing import assert_allclose +from pytest import mark +from pytest import warns from .. io.dst_io import load_dst from .. core.testing_utils import assert_dataframes_close @@ -12,6 +14,7 @@ from .. core.configure import configure from .. core.system_of_units import pes, mm, mus, ns from .. types.symbols import all_events +from .. types.ic_types import NN from . dorothea import dorothea @@ -176,4 +179,41 @@ def test_dorothea_empty_input_file(config_tmpdir, ICDATADIR): conf.update(dict(files_in = PATH_IN, file_out = PATH_OUT)) - dorothea(**conf) + with warns(UserWarning, match="No PMAPs in input file"): + dorothea(**conf) + + +@ignore_warning.no_config_group +def test_dorothea_high_threshold(config_tmpdir, KrMC_pmaps_filename, KrMC_kdst): + # set a ver high Qthr in reconstruction to verify proper behaviour + config = KrMC_kdst[0].config + reco_config = config["global_reco_params"] + reco_config["Qthr"] = 1e6 + + PATH_IN = KrMC_pmaps_filename + PATH_OUT = os.path.join(config_tmpdir, 'test_dorothea_high_threshold.h5') + nrequired = 3 + + conf = configure('dummy invisible_cities/config/dorothea.conf'.split()) + conf.update(dict(files_in = PATH_IN, + file_out = PATH_OUT, + event_range = (0, nrequired), + **config)) + + cnt = dorothea(**conf) + nevt_in = cnt.events_in + nevt_out = cnt.events_out + + if nrequired > 0: + assert nrequired == nevt_in + + dst = load_dst(PATH_OUT, + group = "DST", + node = "Events") + assert len(set(dst.event)) == nevt_out + print(nevt_out, dst.columns) + assert_allclose(dst.S2q, NN) + assert_allclose(dst.X , NN) + assert_allclose(dst.Y , NN) + assert_allclose(dst.Z , NN) + assert_allclose(dst.DT , NN) diff --git a/invisible_cities/cities/eutropia.py b/invisible_cities/cities/eutropia.py index b25b623a1c..92ef1c8a1d 100644 --- a/invisible_cities/cities/eutropia.py +++ b/invisible_cities/cities/eutropia.py @@ -32,6 +32,8 @@ - nevt : number of entries at this bin """ +import warnings + import numpy as np import tables as tb import pandas as pd @@ -135,10 +137,13 @@ def eutropia( files_in : OneOrManyFiles , result = accumulate_psf.future ) - df_writer(h5out, result - , "PSF", "PSFs" - , descriptive_string = f"PSF with {bin_size_xy} mm bin size" - ) + if result is None: + warnings.warn("No output file produced because no events were processed", UserWarning) + else: + df_writer(h5out, result + , "PSF", "PSFs" + , descriptive_string = f"PSF with {bin_size_xy} mm bin size" + ) return result diff --git a/invisible_cities/cities/eutropia_test.py b/invisible_cities/cities/eutropia_test.py index 5ea82dea37..51ab98e72e 100644 --- a/invisible_cities/cities/eutropia_test.py +++ b/invisible_cities/cities/eutropia_test.py @@ -5,6 +5,7 @@ import pandas as pd from pytest import fixture +from pytest import warns from .. cities.eutropia import eutropia from .. core .configure import configure @@ -13,11 +14,6 @@ from .. core .testing_utils import ignore_warning -@fixture(scope="module") -def kr_hits(ICDATADIR): - return os.path.join(ICDATADIR, "Kr_hits_for_psf_*.h5") - - @fixture(scope="module") def fast_psf_config(ICDATADIR): return dict( files_in = os.path.join(ICDATADIR, "Kr_hits_for_psf_0.h5") @@ -119,3 +115,17 @@ def test_eutropia_exact_result(ICDATADIR, output_tmpdir): got = getattr( output_file.root, table) expected = getattr(true_output_file.root, table) assert_tables_equality(got, expected) + + +def test_eutropia_skips_bad_files(ICDATADIR, config_tmpdir): + PATH_IN = os.path.join(ICDATADIR , 'empty_file.h5') + PATH_OUT = os.path.join(config_tmpdir, 'test_eutropia_skips_bad_files.h5') + + config_file = 'dummy invisible_cities/config/eutropia.conf' + conf = configure(config_file.split()) + conf.update(dict(files_in = PATH_IN, + file_out = PATH_OUT)) + + # just check that it doesn't complain + with warns(UserWarning, match="No output file produced because no events were processed"): + eutropia(**conf) diff --git a/invisible_cities/cities/phyllis.py b/invisible_cities/cities/phyllis.py index ae63a435ce..68cf9449e6 100644 --- a/invisible_cities/cities/phyllis.py +++ b/invisible_cities/cities/phyllis.py @@ -87,7 +87,6 @@ def phyllis( files_in : OneOrManyFiles if proc_mode is PMTCalibMode.gain : proc = pmt_deconvolver (detector_db, run_number, n_baseline ) elif proc_mode is PMTCalibMode.gain_maw : proc = pmt_deconvolver_maw(detector_db, run_number, n_baseline, n_maw) elif proc_mode is PMTCalibMode.gain_nodeconv: proc = mode_subtractor (detector_db, run_number) - else : raise ValueError(f"Unrecognized processing mode: {proc_mode}") bin_edges = np.arange(min_bin, max_bin, bin_width) bin_centres = shift_to_bin_centers(bin_edges) diff --git a/invisible_cities/cities/sophronia_test.py b/invisible_cities/cities/sophronia_test.py index 3b8be279e7..73f6231fc0 100644 --- a/invisible_cities/cities/sophronia_test.py +++ b/invisible_cities/cities/sophronia_test.py @@ -10,15 +10,24 @@ @ignore_warning.no_config_group -def test_sophronia_runs(sophronia_config, config_tmpdir): - path_out = os.path.join(config_tmpdir, 'sophronia_runs.h5') +@mark.parametrize("corrections", (True, False)) +def test_sophronia_runs(sophronia_config, config_tmpdir, corrections): + path_out = os.path.join(config_tmpdir, f'sophronia_runs_{corrections}.h5') nevt_req = 1 config = dict(**sophronia_config) config.update(dict(file_out = path_out, event_range = nevt_req)) + # ensure that `sophronia_config` uses corrections and that we + # remove them in the other case + if corrections: + assert "corrections" in config + assert config["corrections"] is not None + else: + del config["corrections"] + cnt = sophronia(**config) - assert cnt.events_in == nevt_req + assert cnt.events_in == nevt_req @ignore_warning.no_config_group diff --git a/invisible_cities/cities/trude.py b/invisible_cities/cities/trude.py index 8fa8343c8a..9ea60a8b01 100644 --- a/invisible_cities/cities/trude.py +++ b/invisible_cities/cities/trude.py @@ -73,9 +73,6 @@ def trude( files_in : OneOrManyFiles , integral_width : float , integrals_period : float ): - if proc_mode not in SiPMCalibMode: - raise ValueError(f"Unrecognized processing mode: {proc_mode}") - bin_edges = np.arange(min_bin, max_bin, bin_width) bin_centres = shift_to_bin_centers(bin_edges) sd = sensor_data(files_in[0], WfType.rwf) diff --git a/invisible_cities/conftest.py b/invisible_cities/conftest.py index 9080d759d6..d21930ad0d 100644 --- a/invisible_cities/conftest.py +++ b/invisible_cities/conftest.py @@ -77,15 +77,6 @@ def electron_MCRD_file(request, ICDATADIR): def krypton_MCRD_file(request, ICDATADIR): return os.path.join(ICDATADIR, request.param) -@pytest.fixture(scope = 'session', - params = ['mcfile_nohits.sim.h5']) -def nohits_sim_file(request, ICDATADIR): - return os.path.join(ICDATADIR, request.param) - -@pytest.fixture(scope = 'session', - params = ['mcfile_sns_only.sim.h5']) -def sns_only_sim_file(request, ICDATADIR): - return os.path.join(ICDATADIR, request.param) @pytest.fixture(scope = 'session', params = ['Kr83_full_nexus_v5_03_01_ACTIVE_7bar_1evt.sim.h5']) @@ -93,14 +84,6 @@ def full_sim_file(request, ICDATADIR): return os.path.join(ICDATADIR, request.param) -@pytest.fixture(scope='session') -def mc_all_hits_data(krypton_MCRD_file): - number_of_hits = 8 - evt_number = 2 - efile = krypton_MCRD_file - return efile, number_of_hits, evt_number - - def _get_pmaps_dict_and_event_numbers(filename): dict_pmaps = load_pmaps(filename) @@ -154,33 +137,6 @@ def correction_map_filename(ICDATADIR): test_file = os.path.join(ICDATADIR, test_file) return test_file -@pytest.fixture(scope='session') -def correction_map_MC_filename(ICDATADIR): - test_file = "kr_emap_xy_100_100_mc.h5" - test_file = os.path.join(ICDATADIR, test_file) - return test_file - - -@pytest.fixture(scope='session') -def KrMC_hdst_filename(ICDATADIR): - test_file = "Kr83_nexus_v5_03_00_ACTIVE_7bar_10evts_HDST.h5" - test_file = os.path.join(ICDATADIR, test_file) - return test_file - - -@pytest.fixture(scope='session') -def KrMC_hdst_filename_toy(ICDATADIR): - test_file = "toy_MC_HDST.h5" - test_file = os.path.join(ICDATADIR, test_file) - return test_file - - -@pytest.fixture(scope='session') -def esmeralda_tracks(ICDATADIR): - test_file = "esmeralda_tracks.hdf5" - test_file = os.path.join(ICDATADIR, test_file) - return test_file - @pytest.fixture(scope='session') def data_hdst(ICDATADIR): @@ -603,13 +559,6 @@ def TlMC_hits_skipping_NN(ICDATADIR): hits = load_hits_skipping_NN(hits_file_name) return hits -@pytest.fixture(scope='session') -def TlMC_hits_merged(ICDATADIR): - hits_file_name = "dst_NEXT_v1_00_05_Tl_ACTIVE_100_0_7bar_DST_10_merged.h5" - hits_file_name = os.path.join(ICDATADIR, hits_file_name) - hits = load_hits(hits_file_name) - return hits - @pytest.fixture(scope="session") def Th228_pmaps(ICDATADIR): @@ -803,22 +752,6 @@ def voxels_toy_data(ICDATADIR): return voxels_filename, (event, X, Y, Z, E, size) -@pytest.fixture(scope='session') -def dbdemopp(): - return 'demopp' - -@pytest.fixture(scope='session') -def dbnew(): - return 'new' - -@pytest.fixture(scope='session') -def dbnext100(): - return 'next100' - -@pytest.fixture(scope='session') -def dbflex100(): - return 'flex100' - @pytest.fixture(scope='session', params=[db_data('demopp' , 3, 256, 3, 79), db_data('new' , 12, 1792, 3, 79), @@ -881,7 +814,7 @@ def pytest_configure(config): def pytest_collection_modifyitems(config, items): - if config.getoption("--runslow"): + if config.getoption("--runslow"): # pragma: no cover # --runslow given in cli: do not skip slow tests return skip_slow = pytest.mark.skip(reason="need --runslow option to run") diff --git a/invisible_cities/core/configure.py b/invisible_cities/core/configure.py index fad09f1215..3f62b7af2a 100644 --- a/invisible_cities/core/configure.py +++ b/invisible_cities/core/configure.py @@ -183,7 +183,7 @@ def compare_signature_to_values( function : Callable """ def get_name(type_): try : return type_.__name__ - except: return str(type_) + except: return str(type_) # pragma: no cover signature = inspect.signature(function) function_name = function.__name__ @@ -219,7 +219,7 @@ def get_name(type_): msg = (f"The function `{function_name}` is missing an argument " f"`{name}` of type `{get_name(type_expected)}`") raise ValueError(msg) - else: + else: # pragma: no cover continue value = kw_values[name] @@ -305,7 +305,7 @@ def push_file(self, file_name): def pop_file(self): self._file_stack.pop() - def display(self): + def display(self): # pragma: no cover conf = self.as_namespace longest = max(self._data, key=len) width = len(longest) + 3 diff --git a/invisible_cities/core/configure_test.py b/invisible_cities/core/configure_test.py index 40d21617f1..a4a942c3b4 100644 --- a/invisible_cities/core/configure_test.py +++ b/invisible_cities/core/configure_test.py @@ -8,6 +8,7 @@ from pytest import raises from pytest import warns +from . configure import event_range from . configure import configure from . configure import Configuration from . configure import make_config_file_reader @@ -84,8 +85,8 @@ """ # The values that will be fed into the above. -config_file_spec = dict(files_in = 'electrons_40keV_z25_RWF.h5', - file_out = 'electrons_40keV_z25_PMP.h5', +config_file_spec = dict(files_in = '$ICDATADIR/electrons_40keV_z25_RWF.h5', + file_out = '/tmp/electrons_40keV_z25_PMP.h5', compression = 'ZLIB4', run_number = 23, nprint = 24, @@ -135,6 +136,20 @@ def default_conf(config_tmpdir): return conf_file_name +@mark.parametrize( "input_s result".split() + , ( ("0", 0) + , ("1", 1) + , ("all", EventRange.all) + , ("last", EventRange.last))) +def test_event_range(input_s, result): + assert event_range(input_s) is result + + +def test_event_range_raises(): + with raises(ValueError, match="`--event-range` must be an int, all or last"): + event_range("a_string_that_is_not_an_integer_or_an_EventRange") + + # This test is run repeatedly with different specs. Each spec is a # sequence of command-line variable specifications. The spec should # contain one entry for each variable that is being set on the command @@ -347,7 +362,7 @@ def simple_conf_file_name(tmpdir_factory): @city def dummy(**kwds): - pass + pass # pragma: no cover def test_config_city_fails_without_config_file(): argv = 'dummy'.split() @@ -522,7 +537,7 @@ class D(Enum): member = 0 def f(a: int, b:float, c:str, d: D): - return + return # pragma: no cover pos_values = (1, 2.0, "a_str", D.member) compare_signature_to_values(f, pos_values, {}) @@ -533,7 +548,7 @@ class D(Enum): member = 0 def f(a: int, b:float, c:str, d: D): - return + return # pragma: no cover kwd_values = dict(a = 1, b = 2.0, c = "a_str", d = D.member) compare_signature_to_values(f, (), kwd_values) @@ -544,7 +559,7 @@ class D(Enum): member = 0 def f(a: int, b:float, c:str, d: D): - return + return # pragma: no cover pos_values = (1, 2.0) kwd_values = dict(c = "a_str", d = D.member) @@ -553,7 +568,7 @@ def f(a: int, b:float, c:str, d: D): def test_compare_signature_to_values_duck_match(): def f(a: int, b:float): - return + return # pragma: no cover values = dict(a = int(1), b = int(2)) compare_signature_to_values(f, (), values) @@ -562,7 +577,7 @@ def f(a: int, b:float): @mark.parametrize("seq", (list, tuple, np.array)) def test_compare_signature_to_values_sequences(seq): def f(a: Sequence): - return + return # pragma: no cover values = (seq([0, 1, 2]),) compare_signature_to_values(f, values, {}) @@ -570,7 +585,7 @@ def f(a: Sequence): def test_compare_signature_to_values_missing_without_default(): def f(a: int, b:int, c:int): - return + return # pragma: no cover pos_values = (1,) for arg_name in "bc": @@ -583,7 +598,7 @@ def f(a: int, b:int, c:int): @mark.parametrize("mode", "positional keyword".split()) def test_compare_signature_to_values_missing_with_default(mode): def f(a: int = 1, b:int = 2): - return + return # pragma: no cover pos_values = (1,) if mode == "positional" else () kwd_values = {} if mode == "positional" else dict(b = 1) @@ -592,7 +607,7 @@ def f(a: int = 1, b:int = 2): def test_compare_signature_to_values_unused_arguments(): def f(a: int): - return + return # pragma: no cover values = dict(a=1, b=2) match_str = "Argument .* is not being used by .*" @@ -608,11 +623,11 @@ def f(a: int): , (tuple, (3,)) , (set , {4,}))) @mark.parametrize("type2", (int, str, list, dict, tuple, set)) -def test_compare_signature_to_values_raises(mode, type1, value, type2): +def test_compare_signature_to_values_raises_wrong_type(mode, type1, value, type2): if type1 is type2: return def f(a : type1): - pass + pass # pragma: no cover pos_values = (type2(),) if mode == "positional" else () kwd_values = {} if mode == "positional" else dict(a = type2()) @@ -622,6 +637,15 @@ def f(a : type1): compare_signature_to_values(f, pos_values, kwd_values) +def test_compare_signature_to_values_raises_wrong_number_of_args(): + def f(a: int, b: int): + return # pragma: no cover + + values = 1, 2, 3 + with raises(ValueError, match=f"The function `f` received {len(values)} positional arguments, but it only accepts 2"): + compare_signature_to_values(f, values, {}) + + @mark.parametrize("do_check", (True, False)) def test_check_annotations(do_check): def f(a: int): diff --git a/invisible_cities/core/core_functions.py b/invisible_cities/core/core_functions.py index 7265748147..05ea7cd012 100644 --- a/invisible_cities/core/core_functions.py +++ b/invisible_cities/core/core_functions.py @@ -265,24 +265,6 @@ def farray_from_string(sfl): return np.array(list(map(float, sfl.split()))) -def _rebin_array(arr, stride, met=np.sum, remainder=False): - """ - rebin arr by a factor stride, using method (ex: np.sum or np.mean), keep the remainder in the - last bin or not - """ - lenb = int(len(arr) / int(stride)) - if remainder and len(arr) % stride != 0: - rebinned = np.empty(lenb + 1) - rebinned[-1] = met(arr[lenb*stride:]) - else: - rebinned = np.empty(lenb) - for i in range(lenb): - s = i * stride - f = s + stride - rebinned[i] = met(arr[s:f]) - return rebinned - - def define_window(wf, window_size): """Define a window based on a peak. Takes max plus/minus *window_size*.""" peak = np.argmax(abs(wf - np.mean(wf))) diff --git a/invisible_cities/core/fit_functions_test.py b/invisible_cities/core/fit_functions_test.py index 8a20398c92..f49a0df531 100644 --- a/invisible_cities/core/fit_functions_test.py +++ b/invisible_cities/core/fit_functions_test.py @@ -68,10 +68,19 @@ def test_get_chi2_and_pvalue_gauss_errors(mean, sigma): assert chi2 == approx(1, rel=1e-2) +def test_get_chi2_and_pvalue_auto_errors(): + ndf = 10 + ydata = np.linspace(1, 2, ndf+1) + yfit = ydata + 1 + chi2, pval = fitf.get_chi2_and_pvalue(ydata, yfit, ndf, sigma=None) + + assert np.isclose(chi2, 0.7687714032) + assert np.isclose(pval, 0.6593105139) + + @mark.parametrize("ey", (0, -1, -100)) def test_fit_raises_ValueError_when_negative_or_zero_value_in_sigma(ey): - def dummy(x, m): - return m + def dummy(x, m): pass # no need to do anything x = np.array([0 , 1 ]) y = np.array([4.1, 4.2]) @@ -313,14 +322,28 @@ def test_fixed_parameters(): @mark.parametrize("pars", - ({'amp' : 2, 'mu' : 0, 'sigma' : 2}, - {'fake1' : 2}, + ({'fake1' : 2}, {'fake1' : 1, 'mu' : 0})) -def test_fix_wrong_parameters_raises_error(pars): +def test_fixed_parameters_raises_wrong_parameters(pars): with raises(ValueError): fitf.fixed_parameters(fitf.gauss, **pars) +def test_fixed_parameters_raises_fixed_all(): + with raises(ValueError): + fitf.fixed_parameters(fitf.gauss, area=1, mu=2, sigma=3) + + +def test_fixed_parameters_n_gaussians(): + def f(x, arg1, arg2): + return x + arg1 + arg2 + f.n_gaussians = 0 + fixed_f = fitf.fixed_parameters(f, arg1=10) + + assert fixed_f(-25, 5) == -10 + assert hasattr(fixed_f, "n_gaussians") + + @mark.parametrize(["func"], ((fitf.profileX,), (fitf.profileY,))) @@ -394,7 +417,7 @@ def test_profile1D_full_range_x(func, xdata, ydata, xrange): # (fitf.profileY, # FLOAT_ARRAY(10000, -100, 100), # FLOAT_ARRAY(10000, -500, 0)))) -def test_profile1D_one_bin_missing_x(func, xdata, ydata): +def test_profile1D_one_bin_missing_x(func, xdata, ydata): # pragma: no cover xdata[core.in_range(xdata, -2, 0)] += 5 xp, yp, ye = func(xdata, ydata) assert xp.size == 99 @@ -464,7 +487,8 @@ def test_profile_data_in_edges(func, xdata, ydata, drop_nan): @mark.slow -def test_profileXY_full_range(): +@mark.parametrize("zrange", (None, (.25, .75))) +def test_profileXY(zrange): N = 10000 Nbin = 100 rms = 12**-0.5 @@ -477,7 +501,8 @@ def test_profileXY_full_range(): ydata, zdata, Nbin, - Nbin) + Nbin, + zrange=zrange) assert np.all(abs(zp - 0.5) < 3.00*rms) assert np.all(abs(ze - rms*(Nbin**2/N)**0.5) < eps) diff --git a/invisible_cities/core/random_sampling_test.py b/invisible_cities/core/random_sampling_test.py index 3a024dcb7d..3a6f206cc4 100644 --- a/invisible_cities/core/random_sampling_test.py +++ b/invisible_cities/core/random_sampling_test.py @@ -182,12 +182,12 @@ def run_number(): return 4714 @fixture(scope="module") -def datasipm(run_number, dbnew): - return DataSiPM(dbnew, run_number) +def datasipm(run_number ): + return DataSiPM("new", run_number) @fixture(scope="module", params=[False, True]) -def noise_sampler(request, dbnew, run_number): +def noise_sampler(request, run_number): nsamples = 1000 smear = request.param thr = 0.99 @@ -201,7 +201,7 @@ def noise_sampler(request, dbnew, run_number): 1.45 : 6, 1.55 : 3, np.inf: 9} - return (NoiseSampler(dbnew, run_number, nsamples, smear), + return (NoiseSampler("new", run_number, nsamples, smear), nsamples, smear, thr, true_threshold_counts) diff --git a/invisible_cities/core/system_of_units.py b/invisible_cities/core/system_of_units.py index 247e208630..06b84d16cf 100644 --- a/invisible_cities/core/system_of_units.py +++ b/invisible_cities/core/system_of_units.py @@ -366,5 +366,5 @@ adc = 1 -def celsius(tKelvin): +def celsius(tKelvin): # pragma: no cover return tKelvin - 273.15 diff --git a/invisible_cities/core/tbl_functions_test.py b/invisible_cities/core/tbl_functions_test.py index 7d6ef55c6d..b429e7ffd0 100644 --- a/invisible_cities/core/tbl_functions_test.py +++ b/invisible_cities/core/tbl_functions_test.py @@ -1,8 +1,40 @@ from pytest import mark +from pytest import raises from . import tbl_functions as tbl +@mark.parametrize("name", "NOCOMPR ZLIB1 ZLIB4 ZLIB5 ZLIB9 BLOSC5 BLZ4HC5".split()) +def test_filters(name): + assert tbl.filters(name) is not None + + +def test_filters_None(): + assert tbl.filters(None) is None + + +def test_filters_raises(): + with raises(ValueError): + tbl.filters("this_filter_does_not_exist") + + +""" + if name is None: return None + + try: + level, lib = {"NOCOMPR": (0, None) , + "ZLIB1" : (1, 'zlib') , + "ZLIB4" : (4, 'zlib') , + "ZLIB5" : (5, 'zlib') , + "ZLIB9" : (9, 'zlib') , + "BLOSC5" : (5, 'blosc') , + "BLZ4HC5": (5, 'blosc:lz4hc'), + }[name] + return tb.Filters(complevel=level, complib=lib) + except KeyError: + raise ValueError("Compression option {} not found.".format(name)) +""" + @mark.parametrize('filename first_evt'.split(), (('dst_NEXT_v0_08_09_Co56_INTERNALPORTANODE_74_0_7bar_MCRD_10000.root.h5', 740000), ('NEXT_v0_08_09_Co56_2_0_7bar_MCRD_1000.root.h5' , 2000), diff --git a/invisible_cities/core/testing_utils.py b/invisible_cities/core/testing_utils.py index 7ce18ac0ae..1a0126ad49 100644 --- a/invisible_cities/core/testing_utils.py +++ b/invisible_cities/core/testing_utils.py @@ -58,7 +58,7 @@ def previous_float(x): return np.nextafter(x, -np.inf) -def next_float(x): +def next_float(x): # pragma: no cover """ Return the next float towards +inf. """ @@ -82,10 +82,6 @@ def float_arrays(size = 100, unique = unique) -def FLOAT_ARRAY(*args, **kwargs): - return float_arrays(*args, **kwargs).example() - - def random_length_float_arrays(min_length = 0, max_length = 100, **kwargs ): @@ -155,7 +151,8 @@ def assert_tables_equality(got_table, expected_table, rtol=1e-7, atol=0): is_float = got.dtype.kind == 'f' if is_float: assert_allclose(got, expected, rtol=rtol, atol=atol) else : assert_equal (got, expected) - except: + except: # pragma: no cover + # I think this branch is unreachable print(f"Mismatch in column {col_name} of tables {names}") raise else: @@ -219,3 +216,4 @@ class ignore_warning: unphysical_rate = mark.filterwarnings("ignore:(Zero|Negative) rate") max_time_short = mark.filterwarnings("ignore:`max_time` shorter than `buffer_length`") no_mc_tables = mark.filterwarnings("ignore:File does not contain MC tables.( *)Use positve run numbers for data") + no_pmaps = mark.filterwarnings("ignore:No PMAPs in input file") diff --git a/invisible_cities/core/testing_utils_test.py b/invisible_cities/core/testing_utils_test.py index bce70bdfa1..2d1747ad7d 100644 --- a/invisible_cities/core/testing_utils_test.py +++ b/invisible_cities/core/testing_utils_test.py @@ -75,7 +75,7 @@ def test_assert_tables_equality_fails_different_names(): assert_tables_equality(table1, table2) -def test_assert_tables_equality_fails_different_values(): +def test_assert_tables_equality_fails_different_types(): # modify the type of a column and check that the function picks up # the difference dtypes1 = [('text', 'U10'), ('integer', 'i4'), ('float', 'f4')] diff --git a/invisible_cities/database/db_connection.py b/invisible_cities/database/db_connection.py index b401cdbe39..e05ea080d3 100644 --- a/invisible_cities/database/db_connection.py +++ b/invisible_cities/database/db_connection.py @@ -1,5 +1,3 @@ -from pytest import mark - import sqlite3 import pymysql pymysql.install_as_MySQLdb() @@ -11,7 +9,6 @@ def connect_sqlite(dbfile): return conn_sqlite, cursor_sqlite -@mark.skip(reason='server timeouts cause too many spurious test failures') def connect_mysql(dbname): conn_mysql = pymysql.connect(host="next.ific.uv.es", user='nextreader',passwd='readonly', db=dbname) diff --git a/invisible_cities/database/download_test.py b/invisible_cities/database/download_test.py index 63ef56a254..f9d0c1b1f1 100644 --- a/invisible_cities/database/download_test.py +++ b/invisible_cities/database/download_test.py @@ -8,7 +8,8 @@ from . import download as db -@mark.skip(reason='server timeouts cause too many spurious test failures') + +@mark.xfail(raises=pymysql.err.OperationalError) @mark.parametrize('dbname', 'DEMOPPDB NEWDB NEXT100DB Flex100DB'.split()) def test_create_table_sqlite(dbname, output_tmpdir): dbfile = os.path.join(output_tmpdir, 'db.sqlite3') @@ -16,9 +17,6 @@ def test_create_table_sqlite(dbname, output_tmpdir): if os.path.isfile(dbfile): os.remove(dbfile) - dbname = 'NEXT100DB' - table = 'PmtBlr' - connSqlite = sqlite3.connect(dbfile) connMySql = pymysql.connect(host="next.ific.uv.es", user='nextreader',passwd='readonly', db=dbname) @@ -26,7 +24,7 @@ def test_create_table_sqlite(dbname, output_tmpdir): cursorMySql = connMySql .cursor() cursorSqlite = connSqlite.cursor() - for table in db.tables: + for table in db.common_tables: db.create_table_sqlite(cursorSqlite, cursorMySql, table) @@ -39,6 +37,7 @@ def test_table_assignment(dbname): assert name in db.table_dict[dbname] +@mark.xfail(raises=pymysql.err.OperationalError) @mark.parametrize('dbname', db.dbnames) def test_tables_exist(dbname): connMySql = pymysql.connect(host="next.ific.uv.es", @@ -50,5 +49,3 @@ def test_tables_exist(dbname): for name in db.table_dict[dbname]: assert (name,) in available - - diff --git a/invisible_cities/database/load_db.py b/invisible_cities/database/load_db.py index a7dec46f45..d40f41b4b2 100644 --- a/invisible_cities/database/load_db.py +++ b/invisible_cities/database/load_db.py @@ -146,6 +146,7 @@ def PMTLowFrequencyNoise(db_file, run_number=1e5): @lru_cache(maxsize=10) def RadioactivityData(db_file, version=None): + if db_file != "next100": return conn = sqlite3.connect(get_db(db_file)) diff --git a/invisible_cities/database/load_db_test.py b/invisible_cities/database/load_db_test.py index 1524e96155..ac9a5de928 100644 --- a/invisible_cities/database/load_db_test.py +++ b/invisible_cities/database/load_db_test.py @@ -49,9 +49,9 @@ def test_SiPMNoise(db): assert noise.shape[1] == energy.shape[0] -def test_DetectorGeometry(dbnew): +def test_DetectorGeometry(): """Check Detector Geometry.""" - geo = DB.DetectorGeo(dbnew) + geo = DB.DetectorGeo("new") assert geo['XMIN'][0] == -198 assert geo['XMAX'][0] == 198 assert geo['YMIN'][0] == -198 @@ -204,15 +204,29 @@ def test_RadioactivityData_max_version(): np.testing.assert_allclose(req[column], max_version[column]) -@mark.parametrize("version", (0, 1, 2)) -def test_RadioactivityData_get_version(version): +@mark.parametrize( "version exists".split() + , ( ( 0, False) # these versions don't exist + , ( 1, False) + , ( 2, False) + , ( 3, True) # these one exist + , ( 4, True) + , ( 5, True) # this one doesn't, but it picks up the highest one + )) +def test_RadioactivityData_get_version(version, exists): '''Check that RadioactivityData returns the desired version of data''' db_file = "next100" activity, efficiency = DB.RadioactivityData(db_file, version=version) - conn = sqlite3.connect(DB.get_db(db_file)) + if not exists: + assert len(activity ) == 0 + assert len(efficiency) == 0 + return + + assert len(activity ) > 0 + assert len(efficiency) > 0 + conn = sqlite3.connect(DB.get_db(db_file)) for table, table_name, column in zip(( activity , efficiency ), ( "Activity", "Efficiency"), ("TotalActivity", "MCEfficiency")): diff --git a/invisible_cities/database/test_data/corrupted_file.h5 b/invisible_cities/database/test_data/corrupted_file.h5 new file mode 100644 index 0000000000..3d19141534 --- /dev/null +++ b/invisible_cities/database/test_data/corrupted_file.h5 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:2fb8f205bb95631f926f2c0c408d281c00938ff7b6867da9cc27bc5ecf3497b7 +size 19409 diff --git a/invisible_cities/database/test_data/dst_NEXT_v1_00_05_Tl_ACTIVE_100_0_7bar_DST_10_merged.h5 b/invisible_cities/database/test_data/dst_NEXT_v1_00_05_Tl_ACTIVE_100_0_7bar_DST_10_merged.h5 deleted file mode 100644 index f3339c9a8c..0000000000 --- a/invisible_cities/database/test_data/dst_NEXT_v1_00_05_Tl_ACTIVE_100_0_7bar_DST_10_merged.h5 +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:3fe2354fee3aa164df4ff4ed5424ce2cd79b84c5df030ed048e340e7bc5b83c6 -size 12682 diff --git a/invisible_cities/database/test_data/empty_ddst.h5 b/invisible_cities/database/test_data/empty_ddst.h5 new file mode 100644 index 0000000000..d4571fe8ad --- /dev/null +++ b/invisible_cities/database/test_data/empty_ddst.h5 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:25e3bf62f88422e323514a7e69cb0742ee32391a2704d73a2f0a465c20efb11a +size 839714 diff --git a/invisible_cities/database/test_data/empty_kdst.h5 b/invisible_cities/database/test_data/empty_kdst.h5 new file mode 100644 index 0000000000..4b29eea6e5 --- /dev/null +++ b/invisible_cities/database/test_data/empty_kdst.h5 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:81d51054e988db49c50b02bc82e4120b40ad38f1ce60a243a9c7b8172f675de9 +size 30658 diff --git a/invisible_cities/database/test_data/esmeralda_tracks.hdf5 b/invisible_cities/database/test_data/esmeralda_tracks.hdf5 deleted file mode 100644 index e9bfd09882..0000000000 Binary files a/invisible_cities/database/test_data/esmeralda_tracks.hdf5 and /dev/null differ diff --git a/invisible_cities/database/test_data/kr_emap_xy_100_100_mc.h5 b/invisible_cities/database/test_data/kr_emap_xy_100_100_mc.h5 deleted file mode 100644 index c55faf25ed..0000000000 --- a/invisible_cities/database/test_data/kr_emap_xy_100_100_mc.h5 +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:88637e141f5e21d595c7f3e72c120f7fdb968c5bc1f31f03f7413969e40d87e9 -size 440296 diff --git a/invisible_cities/database/test_data/mcfile_nohits.sim.h5 b/invisible_cities/database/test_data/mcfile_nohits.sim.h5 deleted file mode 100644 index 9483fc8ed3..0000000000 --- a/invisible_cities/database/test_data/mcfile_nohits.sim.h5 +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:10eb6ae70a2431c3253f99e16e4d6b2011ba0dd7acd98f161cca1e457ecbb40d -size 12465048 diff --git a/invisible_cities/database/test_data/mcfile_sns_only.sim.h5 b/invisible_cities/database/test_data/mcfile_sns_only.sim.h5 deleted file mode 100644 index f2436e3d39..0000000000 --- a/invisible_cities/database/test_data/mcfile_sns_only.sim.h5 +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:5b04deb6be000113f06054359f424b7f7cceaa0e8ca5eea7ff9618274e3a425a -size 7875432 diff --git a/invisible_cities/database/test_data/toy_MC_HDST.h5 b/invisible_cities/database/test_data/toy_MC_HDST.h5 deleted file mode 100644 index 260c8aab6e..0000000000 --- a/invisible_cities/database/test_data/toy_MC_HDST.h5 +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:2c844c04c180bd12978a9e600f0ceffcbe313cd9f3e65808703a4a220e250316 -size 86395 diff --git a/invisible_cities/dataflow/dataflow.py b/invisible_cities/dataflow/dataflow.py index 1b41496aeb..d1cc93706d 100644 --- a/invisible_cities/dataflow/dataflow.py +++ b/invisible_cities/dataflow/dataflow.py @@ -371,7 +371,7 @@ def if_tuple_make_pipe(thing): ###################################################################### -if __name__ == '__main__': +if __name__ == '__main__': # pragma: no cover show = sink(print) count_2_fut, count2 = count(); every2 = filter(lambda n:not n%2)(count2) diff --git a/invisible_cities/dataflow/dataflow_test.py b/invisible_cities/dataflow/dataflow_test.py index 96d2da1d86..0c5ef0ff1c 100644 --- a/invisible_cities/dataflow/dataflow_test.py +++ b/invisible_cities/dataflow/dataflow_test.py @@ -68,6 +68,21 @@ def the_operation(n): return n*n assert result == list(map(the_operation, the_source)) +def test_flatmap_item(): + the_source = [{"the_item": i} for i in range(10)] + + def the_operation(n): return range(n) + square = df.flatmap(the_operation, item="the_item") + + result = [] + the_sink = df.sink(result.append, args="the_item") + + df.push(source = the_source, + pipe = square(the_sink)) + + assert result == [j for i in range(10) for j in range(i)] + + def test_pipe(): # The basic syntax requires any element of a pipeline to be passed @@ -194,6 +209,15 @@ def collect(a, b): assert result_b == expected_b +@mark.parametrize("comp", (df.map, df.flatmap)) +def test_map_mutually_exclusive_arguments(comp): + with raises(ValueError): + comp(lambda x:x, item="something", args="something else") + + with raises(ValueError): + comp(lambda x:x, item="something", out="something else") + + def test_count(): # 'count' is an example of a sink which only produces a result @@ -212,6 +236,20 @@ def test_count(): assert count.future.result() == len(the_source) +def test_count_filter(): + count = df.count_filter(lambda x: x>5) + the_source = list(range(10)) + out = [] + the_sink = df.sink(out.append) + df.push(source = the_source, + pipe = df.pipe(count.filter, the_sink)) + + filter_out = count.future.result() + assert out == list(range(6,10)) + assert filter_out.n_passed == 4 + assert filter_out.n_failed == 6 + + def test_push_futures(): # 'push' provides a higher-level interface to using such futures: diff --git a/invisible_cities/detsim/buffer_functions_test.py b/invisible_cities/detsim/buffer_functions_test.py index c6643ab420..5e8f9d7072 100644 --- a/invisible_cities/detsim/buffer_functions_test.py +++ b/invisible_cities/detsim/buffer_functions_test.py @@ -8,6 +8,7 @@ from . buffer_functions import bin_sensors from . buffer_functions import buffer_calculator from . buffer_functions import find_signal_start +from . buffer_functions import pad_safe def test_bin_sensors(mc_waveforms, pmt_ids, sipm_ids): @@ -90,6 +91,26 @@ def test_find_signal_start_correct_index(): assert pulses[0] == thr_bin +def test_pad_safe(): + sensors = np.arange(12).reshape(6,2) + pad_before = 3 + pad_after = 5 + padded = pad_safe(sensors, (pad_before, pad_after)) + + np.testing.assert_array_equal(padded[:, :pad_before ], 0) + np.testing.assert_array_equal(padded[:, -pad_after: ], 0) + np.testing.assert_array_equal(padded[:, pad_before:-pad_after], sensors) + + +def test_pad_safe_empty(): + shape = (0, 10) + sensors = np.empty(shape) + padding = 3, 5 + padded = pad_safe(sensors, padding) + + assert padded.shape == (0, sum(padding) + 1) + + @mark.parametrize("pre_trigger signal_thresh".split(), ((100 * units.mus, 2), (400 * units.mus, 10))) diff --git a/invisible_cities/detsim/light_tables_c_test.py b/invisible_cities/detsim/light_tables_c_test.py index c1d7a245bf..81d1215af0 100644 --- a/invisible_cities/detsim/light_tables_c_test.py +++ b/invisible_cities/detsim/light_tables_c_test.py @@ -3,6 +3,8 @@ import numpy as np +from pytest import mark +from pytest import fixture from pytest import warns from .. database.load_db import DataSiPM @@ -11,6 +13,7 @@ from . light_tables_c import LT_SiPM from . light_tables_c import LT_PMT +from . light_tables import create_lighttable_function from hypothesis import given from hypothesis import settings @@ -18,116 +21,157 @@ from hypothesis.strategies import integers -few_examples = settings(deadline=None, max_examples=100) +@fixture(scope="session") +def default_lt_pmt(get_dfs): + fname, *_ = get_dfs['lt'] + lt = LT_PMT(fname=fname) + npmt = lt.num_sensors + return lt, npmt + + +@fixture(scope="session") +def default_lt_sipm(get_dfs): + fname, *_ = get_dfs['psf'] + datasipm = DataSiPM("new") + lt = LT_SiPM(fname=fname, sipm_database=datasipm) + nsipm = lt.num_sensors + xy = datasipm[list("XY")].values + return lt, nsipm, xy + def test_LT_SiPM_optional_arguments(get_dfs): datasipm = DataSiPM('new') fname, psf_df, psf_conf = get_dfs['psf'] lt = LT_SiPM(fname=fname, sipm_database=datasipm) + #check the values are read from the table assert lt.el_gap_width == psf_conf.loc['EL_GAP' ].astype(float).value assert lt.active_radius == psf_conf.loc['ACTIVE_rad'].astype(float).value + #check optional arguments are set with User Warning with warns(UserWarning): lt = LT_SiPM(fname=fname, sipm_database=datasipm, el_gap_width=2, active_radius=150) assert lt.el_gap_width == 2 assert lt.active_radius == 150 -@few_examples -@given(xs=floats(min_value=-500, max_value=500), - ys=floats(min_value=-500, max_value=500), - sipm_indx=integers(min_value=0, max_value=1500)) -def test_LT_SiPM_values(get_dfs, xs, ys, sipm_indx): - datasipm = DataSiPM('new') + +@given(xf=floats(min_value=-.25, max_value=.25), + yf=floats(min_value=-.25, max_value=.25)) +def test_LT_SiPM_values_within_active(get_dfs, xf, yf, default_lt_sipm): fname, psf_df, psf_conf = get_dfs['psf'] + r_active = psf_conf.loc['ACTIVE_rad'].astype(float).value + lt, nsipm, xy_sipm = default_lt_sipm + + xs = r_active * xf + ys = r_active * yf + xys = np.array([xs, ys]) + r = np.sum(xys**2)**0.5 - r_active = psf_conf.loc['ACTIVE_rad'].astype(float).value - r = np.sqrt(xs**2 + ys**2) psfbins = psf_df.index.values - lenz = psf_df.shape[1] - psf_df = psf_df /lenz - lt = LT_SiPM(fname=fname, sipm_database=datasipm) - x_sipm, y_sipm = datasipm.iloc[sipm_indx][['X', 'Y']] - dist = np.sqrt((xs-x_sipm)**2+(ys-y_sipm)**2) - psf_bin = np.digitize(dist, psfbins)-1 - max_psf = psf_df.index.max() - if (dist>=max_psf) or (r>=r_active): - values = np.zeros(psf_df.shape[1]) - else: - values = (psf_df.loc[psf_bin].values) + psf_df = psf_df / psf_df.shape[1] # some sort of normalization - ltvals = lt.get_values(xs, ys, sipm_indx) - np.testing.assert_allclose(values, ltvals) + indices = np.arange(nsipm) + dist = np.sum((xy_sipm - xys)**2, axis=1)**0.5 + for sipm_indx, d in zip(indices, dist): + if d < psfbins.max(): + psf_bin = np.digitize(d, psfbins) - 1 + expected = psf_df.loc[psf_bin].values + else: + expected = np.zeros(psf_df.shape[1]) + ltvals = lt.get_values(xs, ys, sipm_indx) + assert np.allclose(ltvals, expected) -@few_examples -@given(xs=floats(min_value=-500, max_value=500), - ys=floats(min_value=-500, max_value=500), - pmt_indx=integers(min_value=0, max_value=11)) -def test_LT_PMTs_values(get_dfs, xs, ys, pmt_indx): + +def test_LT_SiPM_values_outside_active(get_dfs, default_lt_sipm): + fname, psf_df, psf_conf = get_dfs['psf'] + r_active = psf_conf.loc['ACTIVE_rad'].astype(float).value + + lt, nsipm, _ = default_lt_sipm + for sipm_indx in range(nsipm): + ltvals = lt.get_values(r_active, r_active, sipm_indx) + assert np.allclose(ltvals, 0) + + +@given(xf=floats(min_value=-.25, max_value=.25), + yf=floats(min_value=-.25, max_value=.25)) +def test_LT_PMTs_values_within_active(get_dfs, xf, yf, default_lt_pmt): fname, lt_df, lt_conf = get_dfs['lt'] r_active = lt_conf.loc['ACTIVE_rad'].astype(float).value - r = np.sqrt(xs**2 + ys**2) - lt = LT_PMT(fname=fname) + xs = r_active * xf + ys = r_active * yf - xs_lt = find_nearest(np.sort(np.unique(lt_df.index.get_level_values('x'))), xs) - ys_lt = find_nearest(np.sort(np.unique(lt_df.index.get_level_values('y'))), ys) - if (r>=r_active): - values = np.array([0]) #the values are one dimension only - else: - values = lt_df.loc[xs_lt, ys_lt].values[pmt_indx] - ltvals = lt.get_values(xs, ys, pmt_indx) - np.testing.assert_allclose(values, ltvals) + lt, npmt = default_lt_pmt + lt_xs = np.sort(np.unique(lt_df.index.get_level_values('x'))) + lt_ys = np.sort(np.unique(lt_df.index.get_level_values('y'))) + xs_lt = find_nearest(lt_xs, xs) + ys_lt = find_nearest(lt_ys, ys) -@few_examples -@given(xs=floats(min_value=-500, max_value=500), - ys=floats(min_value=-500, max_value=500), - pmt_indx=integers(min_value=0, max_value=11)) -def test_LT_PMTs_values_extended_r(get_dfs, xs, ys, pmt_indx): + expected = lt_df.loc[xs_lt, ys_lt].values[:npmt] # discard sum + ltvals = [lt.get_values(xs, ys, pmt)[0] for pmt in range(npmt)] + assert np.allclose(ltvals, expected) + + +def test_LT_PMTs_values_outside_active(get_dfs, default_lt_pmt): fname, lt_df, lt_conf = get_dfs['lt'] r_active = lt_conf.loc['ACTIVE_rad'].astype(float).value - r_new = 2*r_active - r = np.sqrt(xs**2 + ys**2) + + lt, npmt = default_lt_pmt + + expected = np.zeros(npmt) + ltvals = [lt.get_values(r_active, r_active, pmt)[0] for pmt in range(npmt)] + assert np.allclose(ltvals, expected) + + +@mark.parametrize( "xf yf zero".split(), + ( (0.5, 0.5, False), # inside `r_active_lt` + (1.1, 1.1, False), # outside `r_active_lt` but inside `r_active_forced` + (2.0, 2.0, True), # outside `r_active_forced` + )) +def test_LT_PMTs_values_extended_r(get_dfs, xf, yf, zero): + fname, lt_df, lt_conf = get_dfs['lt'] + r_active_lt = lt_conf.loc['ACTIVE_rad'].astype(float).value + r_active_forced = 2 * r_active_lt + with warns(UserWarning): - lt = LT_PMT(fname=fname,active_radius=r_new) - xs_lt = find_nearest(np.sort(np.unique(lt_df.index.get_level_values('x'))), xs) - ys_lt = find_nearest(np.sort(np.unique(lt_df.index.get_level_values('y'))), ys) - if (r>=r_new): - values = np.array([0]) #the values are one dimension only + lt = LT_PMT(fname=fname,active_radius=r_active_forced) + npmt = lt.num_sensors + + xs = r_active_lt * xf + ys = r_active_lt * yf + if zero: + expected = np.zeros(npmt) else: - values = lt_df.loc[xs_lt, ys_lt].values[pmt_indx] - ltvals = lt.get_values(xs, ys, pmt_indx) - np.testing.assert_allclose(values, ltvals) + lt_xs = np.sort(np.unique(lt_df.index.get_level_values('x'))) + lt_ys = np.sort(np.unique(lt_df.index.get_level_values('y'))) + xs_lt = find_nearest(lt_xs, xs) + ys_lt = find_nearest(lt_ys, ys) + expected = lt_df.loc[xs_lt, ys_lt].values[:npmt] # discard sum + ltvals = [lt.get_values(xs, ys, pmt)[0] for pmt in range(npmt)] + assert np.allclose(ltvals, expected) -def test_LTs_non_physical_sensor(get_dfs): - datasipm = DataSiPM('new') - fname, psf_df, psf_conf = get_dfs['psf'] - lt = LT_SiPM(fname=fname, sipm_database=datasipm) - xs = 0 - ys = 0 - sipm_id = len(datasipm) - values = lt.get_values(xs, ys, sipm_id) - assert all(values==0) - fname, lt_df, lt_conf = get_dfs['lt'] - lt = LT_PMT(fname=fname) - pmt_id = 12 - values = lt.get_values(xs, ys, pmt_id) - assert all(values==0) +def test_LTs_non_physical_sensor(get_dfs, default_lt_pmt, default_lt_sipm): + lt, nsipm, xy_sipm = default_lt_sipm + values = lt.get_values(0, 0, nsipm) # nsipm sensor ID does not exist + assert np.all(values==0) + lt, npmt = default_lt_pmt + values = lt.get_values(0, 0, npmt) # npmt sensor ID does not exist + assert np.all(values==0) -from .light_tables import create_lighttable_function -@few_examples + +@settings(deadline=None) @given(xs=floats(min_value=-500, max_value=500), ys=floats(min_value=-500, max_value=500)) -def test_light_tables_pmt_equal(get_dfs, xs, ys): - fname, lt_df, lt_conf = get_dfs['lt'] - lt_c = LT_PMT(fname=fname) - s2_lt = create_lighttable_function(fname) - pmts_ids = np.arange(lt_c.num_sensors) - vals_lt_c = np.concatenate([lt_c.get_values(xs,ys, pmtid)for pmtid in pmts_ids]).flatten() - vals_lt = s2_lt(np.array([xs]), np.array([ys])).flatten() - np.testing.assert_allclose(vals_lt_c, vals_lt) +def test_light_tables_pmt_equal_to_function(get_dfs, default_lt_pmt, xs, ys): + fname, *_ = get_dfs['lt'] + lt, npmt = default_lt_pmt + lt_fun = create_lighttable_function(fname) + + vals_lt = [lt.get_values(xs,ys, pmtid)[0] for pmtid in range(npmt)] + vals_fun = lt_fun(np.array([xs]), np.array([ys])).flatten() + assert np.allclose(vals_fun, vals_lt) diff --git a/invisible_cities/detsim/light_tables_test.py b/invisible_cities/detsim/light_tables_test.py index 2c91adcb40..b9b3c3ea52 100644 --- a/invisible_cities/detsim/light_tables_test.py +++ b/invisible_cities/detsim/light_tables_test.py @@ -8,6 +8,7 @@ from .light_tables import create_lighttable_function from pytest import fixture +from pytest import raises from hypothesis import given, settings from hypothesis.strategies import floats @@ -80,3 +81,13 @@ def test_get_lt_values_s2(lighttable_filenames, xs, ys): expected = np.zeros((1, 12)) np.testing.assert_allclose(S2_LT(np.array([xs]), np.array([ys])), expected) + + +def test_get_lt_values_complains_if_input_mismatch(lighttable_filenames): + s2_lighttable = lighttable_filenames["s2"] + S2_LT = create_lighttable_function(s2_lighttable) + + with raises(Exception, match="input arrays must be of same shape"): + xs = np.arange(5) + ys = np.arange(6) + S2_LT(xs, ys) diff --git a/invisible_cities/detsim/sensor_functions_test.py b/invisible_cities/detsim/sensor_functions_test.py index dcb9d24a3f..dc9c7378dc 100644 --- a/invisible_cities/detsim/sensor_functions_test.py +++ b/invisible_cities/detsim/sensor_functions_test.py @@ -14,11 +14,11 @@ from .. reco import wfm_functions as wfm -def test_cwf_blr(dbnew, electron_MCRD_file): +def test_cwf_blr(electron_MCRD_file): """Check that CWF -> (deconv) (RWF) == BLR within 1 %. """ run_number = 0 - DataPMT = load_db.DataPMT(dbnew, run_number) + DataPMT = load_db.DataPMT("new", run_number) pmt_active = np.nonzero(DataPMT.Active.values)[0].tolist() coeff_blr = abs(DataPMT.coeff_blr.values) coeff_c = abs(DataPMT.coeff_c .values) @@ -42,14 +42,14 @@ def test_cwf_blr(dbnew, electron_MCRD_file): assert diff[0] < 1 @mark.slow -def test_sipm_noise_sampler(dbnew, electron_MCRD_file): +def test_sipm_noise_sampler(electron_MCRD_file): """This test checks that the number of SiPMs surviving a hard energy cut (50 pes) is small (<10). The test exercises the full construction of the SiPM vectors as well as the noise suppression. """ run_number = 0 - DataSiPM = load_db.DataSiPM(dbnew, run_number) + DataSiPM = load_db.DataSiPM("new", run_number) sipm_adc_to_pes = DataSiPM.adc_to_pes.values.astype(np.double) cal_min = 13 @@ -71,7 +71,7 @@ def test_sipm_noise_sampler(dbnew, electron_MCRD_file): assert np.mean(sipm_adc_to_pes[sipm_adc_to_pes>0]) < cal_max sipms_thresholds = sipm_noise_cut * sipm_adc_to_pes - noise_sampler = SiPMsNoiseSampler(dbnew, run_number, SIPMWL, True) + noise_sampler = SiPMsNoiseSampler("new", run_number, SIPMWL, True) # signal in sipm with noise sipmrwf = e40rd.root.sipmrd[event] + noise_sampler.sample() @@ -80,7 +80,7 @@ def test_sipm_noise_sampler(dbnew, electron_MCRD_file): n_sipm = 0 for j in range(sipmzs.shape[0]): if np.sum(sipmzs[j] > 0): - n_sipm+=1 + n_sipm+=1 # pragma: no cover assert n_sipm < max_sipm_with_signal def test_channel_id_to_IC_id(): @@ -92,8 +92,8 @@ def test_channel_id_to_IC_id(): np.array ([0,4,2])) -def test_channel_id_to_IC_id_with_real_data(dbnew): - pmt_df = load_db.DataPMT(dbnew, 0) +def test_channel_id_to_IC_id_with_real_data(): + pmt_df = load_db.DataPMT("new", 0) assert np.array_equal( convert_channel_id_to_IC_id(pmt_df, pmt_df.ChannelID.values), pmt_df.SensorID .values) diff --git a/invisible_cities/evm/event_model.py b/invisible_cities/evm/event_model.py index 40f1719ba4..173583925c 100644 --- a/invisible_cities/evm/event_model.py +++ b/invisible_cities/evm/event_model.py @@ -185,9 +185,6 @@ def __init__(self, voxels : List[Voxel]): self.voxels = voxels self.E = sum(v.E for v in voxels) - @property - def number_of_voxels(self): - return len(self.voxels) def __str__(self): s = "VoxelCollection: (number of voxels = {})\n".format(self.number_of_voxels) @@ -199,88 +196,6 @@ def __repr__(self): return self.__str__() -class Blob: - """A Blob is a collection of Hits with a seed and a radius. """ - def __init__(self, seed: Tuple[float, float, float], - hits : List[BHit], - radius : float, - e_type : HitEnergy = HitEnergy.E) ->None: - self.seed = seed - self.hits = hits - self.E = sum(getattr(h, e_type.value) for h in hits) - self.radius = radius - self.e_type = e_type.value - - @property - def Etype(self): return self.e_type - - def __str__(self): - s = """Blob: (hits = {} \n - seed = {} \n - blob energy = {} \n - blob radius = {} - """.format(self.hits, self.seed, self.energy, self.radius) - - return s - - def __repr__(self): - return self.__str__() - - -class Track(VoxelCollection): - """A track is a collection of linked voxels. """ - def __init__(self, voxels : List[Voxel], blobs: Tuple[Blob, Blob]) ->None: - super().__init__(voxels) - self.blobs = blobs - - def __str__(self): - s = """Track: (number of voxels = {})\n, - voxels = {} \n - blob_a = {} \n - blob_b = {} - """.format(self.number_of_voxels, self.voxels, self.blobs[0], self.blobs[1]) - - return s - - def __repr__(self): - return self.__str__() - - -class TrackCollection(Event): - """A Collection of tracks""" - def __init__(self, event_number, event_time): - Event.__init__(self, event_number, event_time) - self.tracks = [] - - @property - def number_of_tracks(self): - return len(self.tracks) - - def store(self, table): - row = table.row - for i, t in enumerate(self.tracks): - row["event"] = self.event - row["time" ] = self.time - row["track_no"] = i - - for j, voxel in enumerate(t.voxels): - row["voxel_no"] = j - row["X" ] = voxel.X - row["Y" ] = voxel.Y - row["Z" ] = voxel.Z - row["E" ] = voxel.E - - row.append() - - def __str__(self): - s = "{}".format(self.__class__.__name__) - s+= "Track list:" - s2 = [str(trk) for trk in self.tracks] - return s + ''.join(s2) - - __repr__ = __str__ - - class HitCollection(Event): """A Collection of hits""" def __init__(self, event_number, event_time, hits=None): @@ -356,8 +271,8 @@ def fill_defaults(self): for attribute in ["w", "h", "e", "t", "q"]: setattr(self, "S2" + attribute, [np.nan]) - self.Nsipm = 0 - self.qmax = 0 + self.Nsipm = [0] + self.qmax = [0] for attribute in ["X", "Y", "R", "Phi", "Xrms", "Yrms", "Zrms"]: setattr(self, attribute, [np.nan]) diff --git a/invisible_cities/evm/event_model_test.py b/invisible_cities/evm/event_model_test.py index 3b0dfeb57e..c57f4c447c 100644 --- a/invisible_cities/evm/event_model_test.py +++ b/invisible_cities/evm/event_model_test.py @@ -175,3 +175,19 @@ def test_kr_event_attributes(): "X", "Y", "R", "Phi", "Xrms", "Yrms"]: assert getattr(evt, attr) == [] + +def test_kr_event_fill_defaults(): + evt = KrEvent(-1, -1) + evt.nS1 = 0 + evt.nS2 = 0 + evt.fill_defaults() + + for attr in ["S1w", "S1h", "S1e", "S1t", + "S2w", "S2h", "S2e", "S2t", "S2q", + "X", "Y", "R", "Phi", + "Xrms", "Yrms"]: + assert getattr(evt, attr) == [np.nan], attr + + assert evt.Nsipm == [0] + assert evt.DT == [[np.nan]] + assert evt.Z == [[np.nan]] diff --git a/invisible_cities/evm/pmaps.py b/invisible_cities/evm/pmaps.py index 518dc6a9f6..f5c7f39eb8 100644 --- a/invisible_cities/evm/pmaps.py +++ b/invisible_cities/evm/pmaps.py @@ -91,7 +91,7 @@ def _check_valid_input(self, times, pmts, sipms): if length_times == 0: raise ValueError("Attempt to initialize an empty" - f"{self.__class__.__name} instance.") + f"{self.__class__.__name__} instance.") if ((length_times != length_pmts ) or (length_times != length_sipms != 0)): msg = "Shapes don't match!\n" diff --git a/invisible_cities/evm/pmaps_test.py b/invisible_cities/evm/pmaps_test.py index 085de04de9..1d4a9b6f60 100644 --- a/invisible_cities/evm/pmaps_test.py +++ b/invisible_cities/evm/pmaps_test.py @@ -277,16 +277,15 @@ def test_Peak_width_above_threshold_max_zero(pks): assert peak.width_above_threshold(peak.height) == 0 -@given(peaks(), floats(wf_min, wf_max)) -def test_Peak_width_above_threshold(pks, thr): +@given(peaks(), floats(0.1, 0.9)) +def test_Peak_width_above_threshold_regular(pks, thr_fraction): _, peak = pks - i_above_thr = _get_indices_above_thr(peak.pmts, thr) - if np.size(i_above_thr) > 0: - imin = i_above_thr[0] - imax = i_above_thr[-1] - expected = np.sum(peak.bin_widths[imin:imax+1]) - else: - expected = 0 + thr = thr_fraction * peak.pmts.sum_over_sensors.max() + i_above_thr = _get_indices_above_thr(peak.pmts, thr) + + imin = i_above_thr[0] + imax = i_above_thr[-1] + expected = np.sum(peak.bin_widths[imin:imax+1]) assert peak.width_above_threshold(thr) == approx(expected) @@ -325,6 +324,9 @@ def test_Peak_raises_exception_when_shapes_dont_match(PK, sr1, sr2): PK(np.empty(n_samples + 1), np.empty(n_samples + 1), sr1, sr2) + with raises(ValueError): + PK([], [], sr1, sr2) + @given(pmaps()) def test_PMap_s1s(pmps): diff --git a/invisible_cities/filters/s1s2_filter_test.py b/invisible_cities/filters/s1s2_filter_test.py index 0d491a5369..74c1df8729 100644 --- a/invisible_cities/filters/s1s2_filter_test.py +++ b/invisible_cities/filters/s1s2_filter_test.py @@ -1,5 +1,7 @@ import numpy as np +from pytest import mark +from pytest import raises from pytest import fixture from hypothesis import given from hypothesis.strategies import lists @@ -203,6 +205,55 @@ def test_s12selector_select_s2(selector_conf, truth = (True, False, False, False, False) assert tuple(selector_output) == truth + +@mark.parametrize( "left right".split() + , ( (False, False) + , (False, True) + , ( True, False) + , ( True, True))) +def test_s12selectoroutput_andor_overall(left, right): + output1 = S12SelectorOutput( left, [], []) + output2 = S12SelectorOutput(right, [], []) + + assert (output1 & output2).passed == (left & right) + assert (output1 | output2).passed == (left | right) + + +def test_s12selectoroutput_andor_peaks(): + output1 = S12SelectorOutput( True + , ( True, False, False, True) + , ( True, True, False, True)) + output2 = S12SelectorOutput( True + , (False, True, False, True) + , (False, True, True, False)) + out_and = S12SelectorOutput( True + , (False, False, False, True) + , (False, True, False, False)) + out_or = S12SelectorOutput( True + , ( True, True, False, True) + , ( True, True, True, True)) + + anded = output1 & output2 + ored = output1 | output2 + assert anded.s1_peaks == out_and.s1_peaks + assert anded.s2_peaks == out_and.s2_peaks + assert ored.s1_peaks == out_or .s1_peaks + assert ored.s2_peaks == out_or .s2_peaks + + +@mark.parametrize( "s1 s2".split(), + ( (2, 1) # same s2, different s1 + , (0, 1) + , (1, 0) # same s1 different s2 + , (1, 2))) +def test_s12selectoroutput_invalid_operations(s1, s2): + output1 = S12SelectorOutput(True, (True,)*1 , (False,)*1 ) + output2 = S12SelectorOutput(True, (True,)*s1, (False,)*s2) + + with raises(ValueError): output1 & output2 + with raises(ValueError): output1 | output2 + + def test_pmap_filter(selector_conf, true_s1_peak, small_s1_peak, diff --git a/invisible_cities/filters/trigger_filters.py b/invisible_cities/filters/trigger_filters.py index 961a1343c3..91943b6e10 100644 --- a/invisible_cities/filters/trigger_filters.py +++ b/invisible_cities/filters/trigger_filters.py @@ -7,7 +7,7 @@ last revised: JJGC, July-2017 """ -def TriggerFilter(trigger_params): +def TriggerFilter(trigger_params): # pragma: no cover """Trigger Filter module""" def trigger_filter(peak_data : '{channel_no: s2}'): min_charge, max_charge = trigger_params.charge diff --git a/invisible_cities/io/channel_param_io.py b/invisible_cities/io/channel_param_io.py index 790d79a7e5..5c3d860f1b 100644 --- a/invisible_cities/io/channel_param_io.py +++ b/invisible_cities/io/channel_param_io.py @@ -97,7 +97,7 @@ def basic_param_reader(h5in): return table_names, param_names, param_tables except tb.NoSuchNodeError: print('File does not contain FITPARAMS node') - exit() + raise def generator_param_reader(h5in, table_name): @@ -106,29 +106,26 @@ def generator_param_reader(h5in, table_name): loops over the table requested yielding sensor number, (parameter value dict, parameter error dict) """ - table_names, param_names, param_tables = basic_param_reader(h5in) try: indx = table_names.index(table_name) - for row in param_tables[indx].iterrows(): - yield row['SensorID'], parameters_and_errors(row, param_names[indx][1:]) except ValueError: print('Requested table not present') - exit() + raise + return all_channel_value_reader(param_tables[indx], param_names[indx][1:]) def subset_param_reader(h5in, table_name, param_names): - table_names, _, param_tables = basic_param_reader(h5in) try: indx = table_names.index(table_name) - for row in param_tables[indx].iterrows(): - yield row['SensorID'], parameters_and_errors(row, param_names) except ValueError: print('Requested table not present') - exit() + raise + + return all_channel_value_reader(param_tables[indx], param_names) def all_channel_value_reader(param_table, param_names): @@ -147,11 +144,11 @@ def single_channel_value_reader(channel, param_table, param_names): """ channel_params = param_table.read_where('SensorID=='+str(channel)) if channel_params.size == 0: - print('Sensor info not found') - exit() - elif channel_params.size > 1: - print('Ambiguous call, more than one sensor entry found') - exit() + raise ValueError('Sensor info not found') + + elif channel_params.size > 1: # pragma: no cover + raise ValueError('Ambiguous call, more than one sensor entry found') + return parameters_and_errors(channel_params[0], param_names) @@ -163,4 +160,3 @@ def parameters_and_errors(table_row, parameters): error_dict[pname] = table_row[pname][1] return param_dict, error_dict - diff --git a/invisible_cities/io/channel_param_io_test.py b/invisible_cities/io/channel_param_io_test.py index 42cbde8fd4..c5631a8ee0 100644 --- a/invisible_cities/io/channel_param_io_test.py +++ b/invisible_cities/io/channel_param_io_test.py @@ -3,20 +3,26 @@ import numpy as np import tables as tb +from pytest import fixture +from pytest import raises + from numpy.testing import assert_allclose from .. core import tbl_functions as tbl -from . channel_param_io import generic_params -from . channel_param_io import store_fit_values -from . channel_param_io import channel_param_writer -from . channel_param_io import make_table_dictionary -from . channel_param_io import basic_param_reader -from . channel_param_io import generator_param_reader +from . channel_param_io import generic_params +from . channel_param_io import store_fit_values +from . channel_param_io import channel_param_writer +from . channel_param_io import make_table_dictionary +from . channel_param_io import basic_param_reader +from . channel_param_io import generator_param_reader +from . channel_param_io import subset_param_reader +from . channel_param_io import single_channel_value_reader -def test_generic_parameters(config_tmpdir): - filename = os.path.join(config_tmpdir, 'test_param.h5') +@fixture(scope="session") +def fake_channel_data(config_tmpdir): + filename = os.path.join(config_tmpdir, 'fake_channel_data.h5') out_dict = {} val_list = [] @@ -33,38 +39,34 @@ def test_generic_parameters(config_tmpdir): pWrite(sens, out_dict) val_list.append(list(out_dict.values())) + val_list = np.array(val_list) + return filename, n_rows, out_dict, val_list + + +def test_basic_param_reader_raises(config_tmpdir): + filename = os.path.join(config_tmpdir, "test_basic_param_reader.h5") + with tb.open_file(filename, "w"): pass + + with tb.open_file(filename) as file: + with raises(tb.NoSuchNodeError): + basic_param_reader(file) + + +def test_basic_param_reader(fake_channel_data): + filename, n_rows, out_dict, val_list = fake_channel_data col_list = ['SensorID'] + list(out_dict.keys()) - val_list = np.concatenate(val_list) with tb.open_file(filename) as data_in: - tbl_names, param_names, tbls = basic_param_reader(data_in) assert len(tbls) == 1 assert tbls[0].nrows == n_rows assert len(param_names[0]) == len(col_list) assert param_names[0] == col_list all_values = [ list(x)[1:] for x in tbls[0][:] ] - assert_allclose(np.concatenate(all_values), val_list) + assert_allclose(all_values, val_list) -def test_generator_param_reader(config_tmpdir): - filename = os.path.join(config_tmpdir, 'test_param.h5') - - out_dict = {} - val_list = [] - n_rows = 5 - with tb.open_file(filename, 'w') as data_out: - pWrite = channel_param_writer(data_out, sensor_type="test", - func_name="generic", - param_names=generic_params) - - for sens in range(n_rows): - for i, par in enumerate(generic_params): - out_dict[par] = [i, (i + sens) / 10] - - pWrite(sens, out_dict) - val_list.append(list(out_dict.values())) - - val_list = np.array(val_list) +def test_generator_param_reader(fake_channel_data): + filename, n_rows, out_dict, val_list = fake_channel_data with tb.open_file(filename) as data_in: counter = 0 for sens, (vals, errs) in generator_param_reader(data_in, 'FIT_test_generic'): @@ -76,9 +78,49 @@ def test_generator_param_reader(config_tmpdir): assert counter == n_rows +def test_generator_param_reader_raises(fake_channel_data): + filename, *_ = fake_channel_data + with tb.open_file(filename) as file: + with raises(ValueError): + tuple(generator_param_reader(file, "this_table_does_not_exist")) # consume iterator + + +def test_subset_param_reader(fake_channel_data): + filename, n_rows, _, val_list = fake_channel_data + with tb.open_file(filename) as file: + values_and_errors = tuple(subset_param_reader(file, "FIT_test_generic", ("gain",))) + assert len(values_and_errors) == n_rows + for got, expected in zip(values_and_errors, val_list[:, 4]): + assert np.isclose(got[1][0]["gain"], expected[0]) + assert np.isclose(got[1][1]["gain"], expected[1]) + + +def test_subset_param_reader_raises(fake_channel_data): + filename, *_ = fake_channel_data + with tb.open_file(filename) as file: + with raises(ValueError): + tuple(subset_param_reader(file, "this_table_does_not_exist", ())) # consume iterator + + +def test_single_channel_value_reader(fake_channel_data): + filename, *_ = fake_channel_data + with tb.open_file(filename) as file: + params_table = file.root.FITPARAMS.FIT_test_generic + params, errors = single_channel_value_reader(0, params_table, generic_params) + assert len(params) == 8 + assert len(errors) == 8 + + +def test_single_channel_value_reader_wrong_id(fake_channel_data): + filename, *_ = fake_channel_data + with tb.open_file(filename) as file: + params_table = file.root.FITPARAMS.FIT_test_generic + with raises(ValueError): + single_channel_value_reader(-1, params_table, generic_params) + def test_simple_parameters_with_covariance(config_tmpdir): - filename = os.path.join(config_tmpdir, 'test_param.h5') + filename = os.path.join(config_tmpdir, 'test_param_cov.h5') simple = ["par0", "par1", "par2"] cov = np.array([[0, 1, 2], [3, 4, 5]]) @@ -95,13 +137,11 @@ def test_simple_parameters_with_covariance(config_tmpdir): pWrite(0, out_dict) with tb.open_file(filename) as data_in: - file_cov = data_in.root.FITPARAMS.FIT_test_simple[0]["covariance"] assert_allclose(file_cov, cov) def test_make_table_dictionary(): - param_names = ["par0", "par1", "par2"] par_dict = make_table_dictionary(param_names) @@ -113,7 +153,7 @@ def test_make_table_dictionary(): def test_store_fit_values(config_tmpdir): - filename = os.path.join(config_tmpdir, 'test_param.h5') + filename = os.path.join(config_tmpdir, 'test_param_fit_values.h5') dummy_dict = make_table_dictionary(['par0']) @@ -134,4 +174,3 @@ def test_store_fit_values(config_tmpdir): assert tblRead.nrows == 1 assert len(tblRead.colnames) == len(dummy_dict) assert tblRead.colnames == list(dummy_dict.keys()) - diff --git a/invisible_cities/io/dst_io_test.py b/invisible_cities/io/dst_io_test.py index 3be551aa75..e33a034c15 100644 --- a/invisible_cities/io/dst_io_test.py +++ b/invisible_cities/io/dst_io_test.py @@ -134,6 +134,12 @@ def test_load_dst_converts_from_bytes(ICDATADIR, fixed_dataframe): assert all(config[col].apply(type) != bytes) +def test_load_dst_warns_corrupted_file(ICDATADIR): + wrong_file = os.path.join(ICDATADIR, "corrupted_file.h5") + with pytest.warns(UserWarning, match='.* corrupted: file .*'): + load_dst(wrong_file, "DST", "Events", ignore_errors=True) + + @given(df=dataframe) def test_df_writer_exact(config_tmpdir, df): filename = config_tmpdir + 'dataframe_to_table_exact.h5' diff --git a/invisible_cities/io/mcinfo_io.py b/invisible_cities/io/mcinfo_io.py index 3a913029a6..499d3a9d93 100644 --- a/invisible_cities/io/mcinfo_io.py +++ b/invisible_cities/io/mcinfo_io.py @@ -203,8 +203,6 @@ def read_mc_tables(file_in : str , elif tbl is MCTableType.string_map : str_map = load_mcstringmap(file_in) tbl_dict[tbl] = str_map - else : - raise TypeError("MC table has no reader") return tbl_dict @@ -736,7 +734,7 @@ def get_sensor_types(file_name : str) -> pd.DataFrame: Sensor position info for the MC sensors which saw light in this simulation. """ - if is_oldformat_file(file_name): + if is_oldformat_file(file_name): # pragma: no cover raise TypeError('Old format files not valid for get_sensor_types') sns_pos = load_dst(file_name, 'MC', 'sns_positions') sns_pos.drop(['x', 'y', 'z'], axis=1, inplace=True) @@ -822,7 +820,7 @@ def load_mcsensors_dfold(file_name : str) -> pd.DataFrame: return sns -def get_mc_info(h5in): +def get_mc_info(h5in): # pragma: no cover """Return MC info bank""" extents = h5in.root.MC.extents @@ -879,7 +877,7 @@ def convert_timebin_to_time(sns_resp : pd.DataFrame, def read_mcinfo_evt (mctables: (tb.Table, tb.Table, tb.Table, tb.Table), event_number: int, last_row=0, - return_only_hits: bool=False) -> ([tb.Table], [tb.Table], [tb.Table]): + return_only_hits: bool=False) -> ([tb.Table], [tb.Table], [tb.Table]): # pragma: no cover h5extents = mctables[0] h5hits = mctables[1] h5particles = mctables[2] @@ -925,7 +923,7 @@ def read_mcinfo_evt (mctables: (tb.Table, tb.Table, tb.Table, tb.Table), event_n return hit_rows, particle_rows, generator_rows -def _read_mchit_info(h5f, event_range=(0, int(1e9))) -> Mapping[int, Sequence[MCHit]]: +def _read_mchit_info(h5f, event_range=(0, int(1e9))) -> Mapping[int, Sequence[MCHit]]: # pragma: no cover """Returns all hits in the event""" mc_info = get_mc_info(h5f) h5extents = mc_info.extents diff --git a/invisible_cities/reco/corrections_test.py b/invisible_cities/reco/corrections_test.py index 3311caef06..d40f54a889 100644 --- a/invisible_cities/reco/corrections_test.py +++ b/invisible_cities/reco/corrections_test.py @@ -1,5 +1,6 @@ import os import numpy as np +import tables as tb from . corrections import maps_coefficient_getter from . corrections import read_maps @@ -155,6 +156,26 @@ def test_read_maps_maps_are_correct(map_filename): assert np.all(maps.e0u == 2 ) assert np.all(maps.lt == 5000 ) assert np.all(maps.ltu == 3 ) + assert maps.t_evol is not None + + +def test_read_maps_data_without_time_evolution(map_filename, config_tmpdir): + tmp_filename = os.path.join(config_tmpdir, "map_without_tevol.h5") + with tb.open_file(map_filename) as inputf: + with tb.open_file(tmp_filename, "w") as outputf: + for group in "chi2 e0 e0u lt ltu mapinfo".split(): + outputf.copy_node(getattr(inputf.root, group) + , outputf.root + , recursive = True) + + maps = read_maps(tmp_filename) + assert np.all(maps.chi2 == 1 ) + assert np.all(maps.e0 == 13000) + assert np.all(maps.e0u == 2 ) + assert np.all(maps.lt == 5000 ) + assert np.all(maps.ltu == 3 ) + assert maps.t_evol is None + @given(random_length_float_arrays(min_value = 1e-4, max_value = 3e4 )) diff --git a/invisible_cities/reco/deconv_functions.py b/invisible_cities/reco/deconv_functions.py index 799d8e08d2..85286fc7bc 100644 --- a/invisible_cities/reco/deconv_functions.py +++ b/invisible_cities/reco/deconv_functions.py @@ -408,7 +408,7 @@ def richardson_lucy(image, psf, satellite_start_iter, satellite_max_size, e_cut, if time_ratio <= 1 or len(image.shape) > 2: convolve_method = fftconvolve else: - convolve_method = convolve + convolve_method = convolve # pragma: no cover image = image.astype(float) psf = psf.astype(float) diff --git a/invisible_cities/reco/deconv_functions_test.py b/invisible_cities/reco/deconv_functions_test.py index b9c9b584f1..407dbc5a7d 100644 --- a/invisible_cities/reco/deconv_functions_test.py +++ b/invisible_cities/reco/deconv_functions_test.py @@ -36,6 +36,7 @@ from .. database.load_db import DataSiPM from scipy.stats import multivariate_normal +from scipy.signal import fftconvolve @pytest.fixture(scope='function') @@ -84,23 +85,41 @@ def test_cut_and_redistribute_df(df): def test_drop_isolated_sensors(): - size = 20 - dist = [10.1, 10.1] - x = random.choices(np.linspace(-200, 200, 41), k=size) - y = random.choices(np.linspace(-200, 200, 41), k=size) - q = np.random.uniform(0, 20, size) - e = np.random.uniform(0, 200, size) - df = pd.DataFrame({'X':x, 'Y':y, 'Q':q, 'E':e}) + """ + Generate a cluster of 9 adjacent sensors + 1 isolated sensor and + check that the function drops it. + """ + p = 10.1 + x = [-p, -p, -p, 0, 0, 0, p, p, p, 3*p] + y = [-p, 0, p, -p, 0, p, -p, 0, p, 3*p] + q = [100] * 10 + e = [100] * 10 + dist = [p, p] + df = pd.DataFrame({'X':x, 'Y':y, 'Q':q, 'E':e}) + drop_function = drop_isolated_sensors(dist, ['E']) df_cut = drop_function(df) - if len(df_cut) > 0: - assert np.isclose(df_cut.E.sum(), df.E.sum()) + assert len(df_cut) == 9 + assert np.isclose(df_cut.E.sum(), df.E.sum()) + assert np.all(df_cut.X >= -p) + assert np.all(df_cut.X <= p) + assert np.all(df_cut.Y >= -p) + assert np.all(df_cut.Y <= p) + + +def test_drop_isolated_sensors_drop_all(): + x = [0] * 10 + y = [0] * 10 + q = [1] * 10 + e = [1] * 10 + dist = [1, 1] + df = pd.DataFrame({'X':x, 'Y':y, 'Q':q, 'E':e}) + + drop_function = drop_isolated_sensors(dist, ['E']) + df_cut = drop_function(df) - for row in df_cut.itertuples(index=False): - n_neighbours = len(df_cut[in_range(df_cut.X, row.X - dist[0], row.X + dist[0]) & - in_range(df_cut.Y, row.Y - dist[1], row.Y + dist[1])]) - assert n_neighbours > 1 + assert df_cut.shape == (0, 4) def test_interpolate_signal(): @@ -335,6 +354,20 @@ def test_nonexact_binning(data_hdst, data_hdst_deconvolved): assert(np.all(check_y % 9 == 0)) +#TODO: implement timeout +def test_richardson_lucy_convergence(): + """ + The test passes if it ends quickly, an indication that it was + terminated by threshold, not by number of iterations. + """ + shape = 50, 50 + obj = np.ones (shape, dtype=float) * 0.5 + psf = np.ones (shape, dtype=float) / shape[0] / shape[1] + image = fftconvolve(obj, psf, "same") + out = richardson_lucy(image, psf, iterations=1000000000, iter_thr=1e-8, **no_satellite_killer) + assert np.all(out>0) + + @mark.parametrize("cut_type", CutType) def test_removing_satellites(sat_arr, cut_type): ''' diff --git a/invisible_cities/reco/dst_functions.py b/invisible_cities/reco/dst_functions.py deleted file mode 100644 index a9ac9cbed0..0000000000 --- a/invisible_cities/reco/dst_functions.py +++ /dev/null @@ -1,21 +0,0 @@ -import numpy as np - -def dst_event_id_selection(data, event_ids): - """Filter a DST by a list of event IDs. - Parameters - ---------- - data : pandas DataFrame - DST to be filtered. - event_ids : list - List of event ids that will be kept. - Returns - ------- - filterdata : pandas DataFrame - Filtered DST - """ - if 'event' in data: - sel = np.isin(data.event.values, event_ids) - return data[sel] - else: - print(r'DST does not have an "event" field. Data returned is unfiltered.') - return data diff --git a/invisible_cities/reco/dst_functions_test.py b/invisible_cities/reco/dst_functions_test.py deleted file mode 100644 index 52494d35c5..0000000000 --- a/invisible_cities/reco/dst_functions_test.py +++ /dev/null @@ -1,27 +0,0 @@ -import numpy as np -import pandas as pd - -from hypothesis import given -from hypothesis.strategies import integers -from hypothesis.strategies import lists -from hypothesis.extra.pandas import columns, data_frames - -from . dst_functions import dst_event_id_selection - - -@given(data_frames(columns=columns(['event'], elements=integers(min_value=-1e5, max_value=1e5))), - lists(integers(min_value=-1e5, max_value=1e5))) -def test_dst_event_id_selection(dst, events): - filtered_dst = dst_event_id_selection(dst, events) - assert set(filtered_dst.event.values) == set(dst.event.values) & set(events) - - -def test_dst_event_id_selection_2(): - data = {'event': [1, 1, 3, 6, 7], 'values': [3, 4, 2, 5, 6]} - filt_data = {'event': [1, 1, 6], 'values': [3, 4, 5]} - - df_data = pd.DataFrame(data=data) - df_filt_data = pd.DataFrame(data=filt_data) - df_real_filt = dst_event_id_selection(df_data, [1, 2, 6, 10]) - - assert np.all(df_real_filt.reset_index(drop=True) == df_filt_data.reset_index(drop=True)) diff --git a/invisible_cities/reco/hits_functions_test.py b/invisible_cities/reco/hits_functions_test.py index 415220a763..c053ba9e04 100644 --- a/invisible_cities/reco/hits_functions_test.py +++ b/invisible_cities/reco/hits_functions_test.py @@ -47,11 +47,6 @@ def list_of_hits(draw): assume(sum((h.Q > 0 for h in list_of_hits if h.Q != NN)) >= 1) return list_of_hits -@composite -def thresholds(draw, min_value=1, max_value=1): - th1 = draw (integers( 10 , 20)) - th2 = draw (integers( th1+1, 30)) - return th1, th2 @given(list_of_hits()) def test_merge_NN_does_not_modify_input(hits): @@ -79,6 +74,20 @@ def test_merge_nn_hits_does_not_leave_nn_hits(hits): assert hit.Q != NN +def test_merge_nn_hits_different_peaks(): + zero = xy(0,0) + # peakno Q X Y z E + hits = [ Hit(0, Cluster(NN, xy( 0, 0), zero, 1), 0, 1, zero) + , Hit(0, Cluster( 1, xy(100, 100), zero, 1), -2, 2, zero) + , Hit(1, Cluster( 1, xy( 0, 0), zero, 1), 1, 3, zero)] + merged = merge_NN_hits(hits, same_peak=False) + + assert len(merged) == 2 + assert_hit_equality(merged[0], hits[1]) # second hit unchanged + assert merged[1].E == hits[0].E + hits[2].E + assert merged[1].Q == hits[2].Q + + @given(list_of_hits(), floats()) def test_threshold_hits_does_not_modify_input(hits, th): hits_org = deepcopy(hits) diff --git a/invisible_cities/reco/paolina_functions.py b/invisible_cities/reco/paolina_functions.py index 99b9552f20..e2db7b3dec 100644 --- a/invisible_cities/reco/paolina_functions.py +++ b/invisible_cities/reco/paolina_functions.py @@ -11,9 +11,6 @@ from .. core.exceptions import NoHits from .. core.exceptions import NoVoxels from .. evm.event_model import BHit -from .. evm.event_model import Track -from .. evm.event_model import Blob -from .. evm.event_model import TrackCollection from .. core import system_of_units as units from .. types.symbols import Contiguity from .. types.symbols import HitEnergy @@ -277,27 +274,6 @@ def blob_centres(track_graph : Graph, radius : float) -> Tuple[Tuple[float, floa return (c1, c2) -def make_tracks(evt_number : float, - evt_time : float, - voxels : List[Voxel], - voxel_dimensions : np.ndarray, - contiguity : Contiguity = Contiguity.CORNER, - blob_radius : float = 30 * units.mm, - energy_type : HitEnergy = HitEnergy.E) -> TrackCollection: - """Make a track collection.""" - tc = TrackCollection(evt_number, evt_time) # type: TrackCollection - track_graphs = make_track_graphs(voxels, contiguity) # type: Sequence[Graph] - for trk in track_graphs: - energy_a, energy_b, hits_a, hits_b = blob_energies_and_hits(trk, blob_radius) - a, b = blob_centres(trk, blob_radius) - blob_a = Blob(a, hits_a, blob_radius, energy_type) # type: Blob - blob_b = Blob(b, hits_b, blob_radius, energy_type) - blobs = (blob_a, blob_b) - track = Track(voxels_from_track_graph(trk), blobs) - tc.tracks.append(track) - return tc - - def drop_end_point_voxels(voxels : Sequence[Voxel], energy_threshold : float, min_vxls : int = 3, diff --git a/invisible_cities/reco/paolina_functions_test.py b/invisible_cities/reco/paolina_functions_test.py index 2149fe6b18..f1d7536e6e 100644 --- a/invisible_cities/reco/paolina_functions_test.py +++ b/invisible_cities/reco/paolina_functions_test.py @@ -2,6 +2,7 @@ from math import sqrt from functools import partial +from copy import deepcopy import numpy as np import networkx as nx @@ -33,6 +34,7 @@ from .. evm.event_model import Voxel from . paolina_functions import bounding_box +from . paolina_functions import voxels_within_radius from . paolina_functions import energy_of_voxels_within_radius from . paolina_functions import find_extrema from . paolina_functions import find_extrema_and_length @@ -47,7 +49,6 @@ from . paolina_functions import voxels_from_track_graph from . paolina_functions import length from . paolina_functions import drop_end_point_voxels -from . paolina_functions import make_tracks from . paolina_functions import get_track_energy from .. core import system_of_units as units @@ -128,11 +129,9 @@ def test_voxelize_hits_should_detect_no_hits(): voxelize_hits([], None) +# TODO: deal with no hits in the input @given(bunch_of_hits) def test_bounding_box(hits): - if not len(hits): # TODO: deal with empty sequences - return - lo, hi = bounding_box(hits) mins = [float(' inf')] * 3 @@ -152,13 +151,11 @@ def test_bounding_box(hits): assert_almost_equal(maxs[d], hi[d]) +# TODO: deal with empty hits in the input @given(bunch_of_hits, box_sizes) def test_voxelize_hits_does_not_lose_energy(hits, voxel_dimensions): voxels = voxelize_hits(hits, voxel_dimensions, strict_voxel_size=False) - if not hits: - assert voxels == [] - def sum_energy(seq): return sum(e.E for e in seq) @@ -314,41 +311,64 @@ def test_find_extrema_no_voxels(): find_extrema_and_length({}) -@fixture(scope='module') -def track_extrema(): - voxel_spec = ((10,10,10, 1000), - (10,10,11, 1), - (10,10,12, 2), - (10,10,13, 4), - (10,10,14, 8), - (10,10,15, 16), - (10,11,15, 32), - (10,12,15, 64), - (10,13,15, 128), - (10,14,15, 265), - (10,15,15, 512), - (11,15,15, 256), - (12,15,15, 128), - (13,15,15, 64), - (14,15,15, 32), - (15,15,15, 16), - (16,16,16, 8), - (17,17,17, 4), - (18,18,18, 2), - (19,19,19, 1), - (20,20,20, 2000), - ) - vox_size = np.array([1,1,1]) - voxels = [Voxel(x,y,z, E, vox_size) for (x,y,z,E) in voxel_spec] - tracks = make_track_graphs(voxels) - - assert len(tracks) == 1 - extrema = find_extrema(tracks[0]) - - assert voxels[ 0] in extrema - assert voxels[-1] in extrema - - return tracks[0], extrema +def test_voxels_within_radius(): + dummy_voxel = lambda e: Voxel(0, 0, 0, e, 1) + r = 3.5 + es = np.arange(10, 20) + ds = np.arange( 0, 10) + + vdist = {dummy_voxel(e): d for e, d in zip(es, ds)} + selected = voxels_within_radius(vdist, r) + assert len(selected) == 4 + assert all(v.E in es[:4] for v in selected) + + +def test_energy_of_voxels_within_radius(): + dummy_voxel = lambda e: Voxel(0, 0, 0, e, 1) + r = 3.5 + es = np.arange(10, 20) + ds = np.arange( 0, 10) + + vdist = {dummy_voxel(e): d for e, d in zip(es, ds)} + tot_energy = energy_of_voxels_within_radius(vdist, r) + assert np.isclose(tot_energy, sum(es[:4])) + + +# @fixture(scope='module') +# def track_extrema(): +# voxel_spec = ((10,10,10, 1000), +# (10,10,11, 1), +# (10,10,12, 2), +# (10,10,13, 4), +# (10,10,14, 8), +# (10,10,15, 16), +# (10,11,15, 32), +# (10,12,15, 64), +# (10,13,15, 128), +# (10,14,15, 265), +# (10,15,15, 512), +# (11,15,15, 256), +# (12,15,15, 128), +# (13,15,15, 64), +# (14,15,15, 32), +# (15,15,15, 16), +# (16,16,16, 8), +# (17,17,17, 4), +# (18,18,18, 2), +# (19,19,19, 1), +# (20,20,20, 2000), +# ) +# vox_size = np.array([1,1,1]) +# voxels = [Voxel(x,y,z, E, vox_size) for (x,y,z,E) in voxel_spec] +# tracks = make_track_graphs(voxels) + +# assert len(tracks) == 1 +# extrema = find_extrema(tracks[0]) + +# assert voxels[ 0] in extrema +# assert voxels[-1] in extrema + +# return tracks[0], extrema def test_voxelize_single_hit(): @@ -488,16 +508,28 @@ def test_energy_is_conserved_with_dropped_voxels(hits, requested_voxel_dimension @mark.parametrize("energy_type", HitEnergy) -@given(hits = bunch_of_corrected_hits(), - requested_voxel_dimensions = box_sizes, - min_voxels = min_n_of_voxels, - fraction_zero_one = fraction_zero_one) -def test_dropped_voxels_have_nan_energy(hits, requested_voxel_dimensions, min_voxels, fraction_zero_one, energy_type): - voxels = voxelize_hits(hits, requested_voxel_dimensions, strict_voxel_size=False, energy_type=energy_type) - energies = [v.E for v in voxels] - e_thr = min(energies) + fraction_zero_one * (max(energies) - min(energies)) - _, dropped_voxels = drop_end_point_voxels(voxels, e_thr, min_voxels) +@given(hits=bunch_of_hits) +def test_dropped_voxels_have_nan_energy(hits, energy_type): + """ + Generate 5 voxels in a straight line. The ones in each end and the + central one have lower energy. Check that only those in the end are + dropped. The hits in these voxels are meaningless, we just enforce + their energy to match the voxel energy. + """ + def make_hits(e): + hs = deepcopy(hits) + for h in hs: + h.E = h.Ec = h.Ep = e/len(hs) + return hs + + vox_size = np.ones(3) + es = [1, 10, 1, 10, 1] + voxels = [Voxel(i, 0, 0, e, vox_size, make_hits(e), energy_type) for i, e in enumerate(es)] + _, dropped_voxels = drop_end_point_voxels(voxels, energy_threshold=5, min_vxls=3) + + assert len(dropped_voxels) == 2 for voxel in dropped_voxels: + assert voxel.X in [voxels[0].X, voxels[-1].X] assert np.isnan(voxel.E) for hit in voxel.hits: assert np.isnan(getattr(hit, energy_type.value)) @@ -509,9 +541,6 @@ def test_dropped_voxels_have_nan_energy(hits, requested_voxel_dimensions, min_vo min_voxels = min_n_of_voxels, fraction_zero_one = fraction_zero_one) def test_drop_end_point_voxels_doesnt_modify_other_energy_types(hits, requested_voxel_dimensions, min_voxels, fraction_zero_one, energy_type): - def energy_from_hits(voxel, e_type): - return [getattr(hit, e_type) for hit in voxel.hits] - voxels = voxelize_hits(hits, requested_voxel_dimensions, strict_voxel_size=False, energy_type=energy_type) voxels = sorted(voxels, key=attrgetter("xyz")) energies = [v.E for v in voxels] @@ -790,56 +819,13 @@ def test_paolina_functions_with_hit_energy_different_from_default_value(hits, re assert np.isclose(tot_energy, tot_mod_energy) - tot_default_energy = sum(h.E for v in voxels_c for h in v.hits) - tot_mod_default_energy = sum(h.E for v in mod_voxels_c for h in v.hits) - # We don't want to modify the default energy of hits, if the voxels are made with energy_c - if len(mod_voxels_c) < len(voxels_c): - assert tot_default_energy > tot_mod_default_energy - - -def test_make_tracks_function(ICDATADIR): - - # Get some test data - hit_file = os.path.join(ICDATADIR, 'tracks_0000_6803_trigger2_v0.9.9_20190111_krth1600.h5') - evt_number = 19 - size = 15. - voxel_size = np.array([size,size,size], dtype=np.float16) - blob_radius = 21*units.mm - - # Read the hits and voxelize - all_hits = load_hits(hit_file) - - for evt_number, hit_coll in all_hits.items(): - evt_hits = hit_coll.hits - evt_time = hit_coll.time - voxels = voxelize_hits(evt_hits, voxel_size, strict_voxel_size=False, energy_type=HitEnergy.E) - - tracks = list(make_track_graphs(voxels)) - - track_coll = make_tracks(evt_number, evt_time, voxels, voxel_size, - contiguity=Contiguity.CORNER, - blob_radius=blob_radius, - energy_type=HitEnergy.E) - tracks_from_coll = track_coll.tracks - - tracks.sort (key=lambda x : len(x.nodes())) - tracks_from_coll.sort(key=lambda x : x.number_of_voxels) - - # Compare the two sets of tracks - assert len(tracks) == len(tracks_from_coll) - for i in range(len(tracks)): - t = tracks[i] - tc = tracks_from_coll[i] - - assert len(t.nodes()) == tc.number_of_voxels - assert sum(v.E for v in t.nodes()) == tc.E - - tc_blobs = list(tc.blobs) - tc_blobs.sort(key=lambda x : x.E) - tc_blob_energies = (tc.blobs[0].E, tc.blobs[1].E) + tot_default_energy = sum(h.E for v in voxels_c for h in v.hits) + tot_mod_default_energy = sum(h.E for v in mod_voxels_c for h in v.hits) - assert np.allclose(blob_energies(t, blob_radius), tc_blob_energies) + no_drop = len(mod_voxels_c) == len(voxels_c) + lost_energy = tot_default_energy > tot_mod_default_energy + assert no_drop ^ lost_energy @given(bunch_of_hits, box_sizes) diff --git a/invisible_cities/reco/peak_functions_test.py b/invisible_cities/reco/peak_functions_test.py index 16efa9bcf7..0da55acf58 100644 --- a/invisible_cities/reco/peak_functions_test.py +++ b/invisible_cities/reco/peak_functions_test.py @@ -51,14 +51,14 @@ def times_and_waveforms(draw): waveforms = draw(multiple_waveforms()) n_samples = waveforms.shape[1] times = draw(arrays(float, n_samples, elements=floats(0, 10*n_samples), unique=True)) - return times, waveforms + widths = np.diff(np.append(times, times[-1] + 1)) + return times, widths, waveforms @composite def rebinned_sliced_waveforms(draw): - times, wfs = draw(times_and_waveforms()) + times, widths, wfs = draw(times_and_waveforms()) assume(times.size >= 5) - widths = np.append(np.diff(times), max(np.diff(times))) indices = np.arange(times.size) first = draw(integers( 0, times.size - 2)) @@ -77,6 +77,7 @@ def rebinned_sliced_waveforms(draw): return times, widths, wfs, indices, times_slice, widths_slice, wfs_slice, rebin +#TODO: rethink this composite @composite def peak_indices(draw): size = draw(integers(10, 50)) @@ -267,17 +268,31 @@ def test_select_peaks_without_bounds(peak_data): assert got == exactly(expected) -@given(peak_indices(), - floats(0, 5), floats( 6, 10), - floats(0, 10), floats(11, 20)) -def test_select_peaks_filtered_out(peak_data, i0, i1, l0, l1): - _, peaks , _ = peak_data - i_limits = minmax(i0, i1) - l_limits = minmax(l0, l1) +def _peak_length(i): + return i[-1] - i[0] + 1 + +@given(peak_indices()) +def test_select_peaks_filters_window(peak_data): + indices, peaks , _ = peak_data + i_limits = minmax(indices.min() + 1, indices.max() - 1) + l_limits = minmax(0, np.inf) # don't filter by length selected = pf.select_peaks(peaks, i_limits * 25 * units.ns, l_limits) for peak in selected: - assert i0 <= peak[0] <= i1 - assert l0 <= peak[-1] + 1 - peak[0] <= l1 + assert i_limits.min <= peak[0] <= i_limits.max + + +@given(peak_indices()) +def test_select_peaks_filters_length(peak_data): + indices, peaks , _ = peak_data + lengths = sorted(map(_peak_length, peaks)) + if max(lengths) - 2 <= min(lengths): + lengths.append(max(lengths)+2) + + i_limits = minmax(0, indices.max() + 1) # don't filter by window + l_limits = minmax(lengths[0] + 1, lengths[-1] - 1) + selected = pf.select_peaks(peaks, i_limits * 25 * units.ns, l_limits) + for peak in selected: + assert l_limits.min <= _peak_length(peak) <= l_limits.max def test_select_peaks_right_length_with_holes(): @@ -298,24 +313,45 @@ def test_select_peaks_right_length_with_holes(): assert selected_wrong == () -@given(peak_indices(), - floats(0, 5), floats( 6, 10), - floats(0, 10), floats(11, 20)) -def test_select_peaks(peak_data, t0, t1, l0, l1): +@given(peak_indices().filter(lambda data: len(data[1]) > 2)) +def test_select_peaks_tlimits(peak_data): + # peaks are time-ordered. Reject the first and the last one _, peaks , _ = peak_data + t0 = peaks[ 0][ 0] + 1 + t1 = peaks[-1][-1] - 1 + i_limits = minmax(t0, t1) t_limits = i_limits * 25 * units.ns + l_limits = minmax(0, 1000000) + selected = pf.select_peaks(peaks, t_limits, l_limits) + + assert len(selected) == len(peaks) - 2 + for got, expected in zip(selected, peaks[1:-1]): + assert got == exactly(expected) + + +@given(peak_indices() + .filter(lambda data: np.unique(list(map(_peak_length, data[1]))).size > 2)) +def test_select_peaks_llimits(peak_data): + _, peaks , _ = peak_data + ls = np.array(list(map(_peak_length, peaks))) + l0 = ls.min() + 1 + l1 = ls.max() - 1 + + i_limits = minmax(0, peaks[-1][-1] + 1) + t_limits = i_limits * 25 * units.ns l_limits = minmax(l0, l1) selected = pf.select_peaks(peaks, t_limits, l_limits) - select = lambda ids: (i_limits.contains(ids[ 0]) and - i_limits.contains(ids[-1]) and - l_limits.contains(ids[-1] + 1 - ids[0])) - expected_peaks = tuple(filter(select, peaks)) + selected_ls = np.array(list(map(_peak_length, selected))) - assert len(selected) == len(expected_peaks) - for got, expected in zip(selected, expected_peaks): - assert got == exactly(expected) + # at least one left, condition set in the input it should have + # filtered at least 2 peaks: the shortest and the longest peaks + assert len(selected) >= 1 + assert len(peaks) - len(selected) >= 2 + + assert np.all(selected_ls >= l0) + assert np.all(selected_ls <= l1) @given(rebinned_sliced_waveforms()) @@ -527,44 +563,32 @@ def test_get_pmap(s1_and_s2_with_indices): assert_PMap_equality(pmap, expected_pmap) -@given(times_and_waveforms(), integers(2, 10)) +@given( times_and_waveforms(), integers(2, 10)) def test_rebin_times_and_waveforms_sum_axis_1_does_not_change(t_and_wf, stride): - times, wfs = t_and_wf - widths = [1] - if len(times) > 1: - widths = np.append(np.diff(times), max(np.diff(times))) - _, _, rb_wfs = pf.rebin_times_and_waveforms(times, widths, wfs, stride) + times, widths, wfs = t_and_wf + _, _, rb_wfs = pf.rebin_times_and_waveforms(times, widths, wfs, stride) assert np.sum(wfs, axis=1) == approx(np.sum(rb_wfs, axis=1)) @given(times_and_waveforms(), integers(2, 10)) def test_rebin_times_and_waveforms_sum_axis_0_does_not_change(t_and_wf, stride): - times, wfs = t_and_wf - widths = [1] - if len(times) > 1: - widths = np.append(np.diff(times), max(np.diff(times))) - sum_wf = np.stack([np.sum(wfs, axis=0)]) - _, _, rb_wfs = pf.rebin_times_and_waveforms(times, widths, wfs, stride) - _, _, rb_sum = pf.rebin_times_and_waveforms(times, widths, sum_wf , stride) + times, widths, wfs = t_and_wf + sum_wf = np.stack([np.sum(wfs, axis=0)]) + _, _, rb_wfs = pf.rebin_times_and_waveforms(times, widths, wfs, stride) + _, _, rb_sum = pf.rebin_times_and_waveforms(times, widths, sum_wf , stride) assert rb_sum[0] == approx(np.sum(rb_wfs, axis=0)) @given(times_and_waveforms(), integers(2, 10)) def test_rebin_times_and_waveforms_number_of_wfs_does_not_change(t_and_wf, stride): - times, wfs = t_and_wf - widths = [1] - if len(times) > 1: - widths = np.append(np.diff(times), max(np.diff(times))) - _, _, rb_wfs = pf.rebin_times_and_waveforms(times, widths, wfs, stride) + times, widths, wfs = t_and_wf + _, _, rb_wfs = pf.rebin_times_and_waveforms(times, widths, wfs, stride) assert len(wfs) == len(rb_wfs) @given(times_and_waveforms(), integers(2, 10)) def test_rebin_times_and_waveforms_number_of_bins_is_correct(t_and_wf, stride): - times, wfs = t_and_wf - widths = [1] - if len(times) > 1: - widths = np.append(np.diff(times), max(np.diff(times))) + times, widths, wfs = t_and_wf rb_times, _, rb_wfs = pf.rebin_times_and_waveforms(times, widths, wfs, stride) expected_n_bins = times.size // stride @@ -577,10 +601,7 @@ def test_rebin_times_and_waveforms_number_of_bins_is_correct(t_and_wf, stride): @given(times_and_waveforms()) def test_rebin_times_and_waveforms_stride_1_does_not_rebin(t_and_wf): - times, wfs = t_and_wf - widths = [1] - if len(times) > 1: - widths = np.append(np.diff(times), max(np.diff(times))) + times, widths, wfs = t_and_wf rb_times, _, rb_wfs = pf.rebin_times_and_waveforms(times, widths, wfs, 1) assert np.all(times == rb_times) @@ -589,10 +610,7 @@ def test_rebin_times_and_waveforms_stride_1_does_not_rebin(t_and_wf): @given(times_and_waveforms(), integers(2, 10)) def test_rebin_times_and_waveforms_times_are_consistent(t_and_wf, stride): - times, wfs = t_and_wf - widths = [1] - if len(times) > 1: - widths = np.append(np.diff(times), max(np.diff(times))) + times, widths, wfs = t_and_wf # The samples falling in the last bin cannot be so easily # compared as the other ones so I remove them. diff --git a/invisible_cities/reco/pmaps_functions.py b/invisible_cities/reco/pmaps_functions.py index 804ac5a269..721af29a65 100644 --- a/invisible_cities/reco/pmaps_functions.py +++ b/invisible_cities/reco/pmaps_functions.py @@ -1,7 +1,5 @@ import numpy as np -from typing import List, Union - from ..evm.pmaps import _Peak from ..evm.pmaps import PMTResponses from ..evm.pmaps import SiPMResponses @@ -9,6 +7,8 @@ from ..core.core_functions import dict_filter_by_key from ..types.symbols import RebinMethod +from typing import List, Union + def get_even_slices(bins : int, stride : int) -> List[slice]: new_bins = int(np.ceil(bins / stride)) diff --git a/invisible_cities/reco/pmaps_functions_test.py b/invisible_cities/reco/pmaps_functions_test.py index 6444d04601..33db099085 100644 --- a/invisible_cities/reco/pmaps_functions_test.py +++ b/invisible_cities/reco/pmaps_functions_test.py @@ -1,3 +1,5 @@ +from itertools import combinations + import numpy as np from pytest import approx @@ -91,31 +93,54 @@ def test_rebin_peak(pk, fraction): assert_Peak_equality(rebinned_pk, expected_pk) -@mark.parametrize("threshold", (4000, 10000)) @given(peaks(subtype=S2)) -def test_rebin_peak_threshold(threshold, pk): +def test_rebin_peak_threshold_below(pk): + """With a threshold below the minimum slice energy, there should be no rebinning""" + _, pk = pk + + threshold = pk.pmts.sum_over_sensors.min() - 1 + rebinned = pmf.rebin_peak(pk, threshold, pmf.RebinMethod.threshold) + assert_Peak_equality(pk, rebinned) + + +@given(peaks(subtype=S2)) +def test_rebin_peak_threshold_above_sum(pk): + """With a threshold above the integrated signal, there should be only one slice""" + _, pk = pk + + threshold = pk.pmts.sum_over_sensors.sum() + rebinned = pmf.rebin_peak(pk, threshold, pmf.RebinMethod.threshold) + + assert len(rebinned.times) == 1 + assert rebinned.total_energy == approx(pk.total_energy) + assert rebinned.total_charge == approx(pk.total_charge) + + +@given(peaks(subtype=S2).filter(lambda p: len(p[1].times)>2)) +def test_rebin_peak_threshold(pk): _, pk = pk pk_eng = pk.total_energy pk_char = pk.total_charge + threshold = pk.pmts.sum_over_sensors.mean() + threshold = np.round(threshold, 6) # avoid precision problems with floats rebinned_pk = pmf.rebin_peak(pk, threshold, pmf.RebinMethod.threshold) assert rebinned_pk.total_energy == approx(pk_eng) assert rebinned_pk.total_charge == approx(pk_char) - if pk_eng < threshold: - assert rebinned_pk.times.shape[0] == 1 - else: - assert np.all(rebinned_pk.pmts.sum_over_sensors >= threshold) + assert np.all(rebinned_pk.pmts.sum_over_sensors >= threshold) -@given(dictionaries(keys = integers(min_value=-1e5, max_value=1e5), - values = pmaps(), - max_size = 5), - lists(integers(min_value=-1e5, max_value=1e5))) +@given(pmaps(), pmaps(), pmaps(), pmaps()) @settings(suppress_health_check=(HealthCheck.too_slow,)) -def test_pmap_event_id_selection(pmaps, events): - filtered_pmaps = pmf.pmap_event_id_selection(pmaps, events) - assert set(filtered_pmaps) == set(pmaps) & set(events) - for evt, pmap in filtered_pmaps.items(): - assert pmap == pmaps[evt] +def test_pmap_event_id_selection(pmap0, pmap1, pmap2, pmap3): + pmaps = {0:pmap0, 12:pmap1, 345:pmap2, 6789:pmap3} + events = list(pmaps) + for chosen in combinations(events, 2): + pm0, pm1 = chosen + filtered_pmaps = pmf.pmap_event_id_selection(pmaps, chosen) + assert len(filtered_pmaps) == 2 + assert pm0 in filtered_pmaps and pm1 in filtered_pmaps + assert filtered_pmaps[pm0] == pmaps[pm0] + assert filtered_pmaps[pm1] == pmaps[pm1] diff --git a/invisible_cities/reco/wfm_functions.py b/invisible_cities/reco/wfm_functions.py index 4ecced9e21..cf6dfa63c1 100644 --- a/invisible_cities/reco/wfm_functions.py +++ b/invisible_cities/reco/wfm_functions.py @@ -27,7 +27,7 @@ def to_adc(wfs, adc_to_pes): return wfs * adc_to_pes.reshape(wfs.shape[0], 1) -def to_pes(wfs, adc_to_pes): +def to_pes(wfs, adc_to_pes): # pragma: no cover """ Convert waveform in adc to pes. diff --git a/invisible_cities/reco/wfm_functions_test.py b/invisible_cities/reco/wfm_functions_test.py index 35ab09a885..41ceccbe91 100644 --- a/invisible_cities/reco/wfm_functions_test.py +++ b/invisible_cities/reco/wfm_functions_test.py @@ -15,7 +15,7 @@ @mark.slow -def test_compare_cwf_blr(dbnew, ICDATADIR): +def test_compare_cwf_blr(ICDATADIR): """Test functions cwf_from_rwf() and compare_cwf_blr(). The test: 1) Computes CWF from RWF (function cwf_from_rwf()) @@ -29,8 +29,8 @@ def test_compare_cwf_blr(dbnew, ICDATADIR): thr_trigger = 5) run_number = 0 - DataPMT = load_db.DataPMT (dbnew, run_number) - DataSiPM = load_db.DataSiPM(dbnew, run_number) + DataPMT = load_db.DataPMT ("new", run_number) + DataSiPM = load_db.DataSiPM("new", run_number) calib = CalibVectors(channel_id = DataPMT .ChannelID .values , coeff_blr = abs(DataPMT .coeff_blr .values), @@ -49,3 +49,26 @@ def test_compare_cwf_blr(dbnew, ICDATADIR): event_list=range(NEVT), window_size=300) assert max(diff) < 0.15 + + +def test_noise_suppression_constant_threshold(): + wfms = np.array([ np.arange(0, 5) + , np.arange(1, 6) + , np.arange(2, 7) ]) + threshold = 3.5 + suppressed = wfm.noise_suppression(wfms, threshold) + expected = np.where(wfms np.ndarray: +def sawteeth_waveform(wfp: WfmPar) -> np.ndarray: # pragma: no cover t_tot = wfp.t_tot t_pre = wfp.t_pre t_rise_s2 = wfp.t_rise_s2 @@ -119,7 +119,7 @@ def simulate_pmt_response(fee : FEE, wf : np.ndarray) -> WaveformPmt: def waveform_generator(fee : FEE, wfp : WfmPar, nsensors=5, pedestal=1024, - random_t0=True) -> WaveformPmt: + random_t0=True) -> WaveformPmt: # pragma: no cover """ Generate PMT-like (e.g, negative swing) and SiPM-like (positive) waveforms. diff --git a/invisible_cities/types/ic_types_test.py b/invisible_cities/types/ic_types_test.py index fe5580b2fc..b7cbefeb8f 100644 --- a/invisible_cities/types/ic_types_test.py +++ b/invisible_cities/types/ic_types_test.py @@ -1,4 +1,5 @@ from string import ascii_letters +from copy import deepcopy import numpy as np @@ -75,12 +76,26 @@ def test_minmax_sub(mm, f): np.isclose (lowered.max , hi - f, rtol=1e-4) +@given(minmaxes) +def test_minmax_eq(mm): + mm2 = deepcopy(mm) + mm3 = deepcopy(mm) + mm2.min = mm.min + 1 + mm3.max = mm.max + 1 + + assert mm == mm + assert mm != mm2 + assert mm != mm3 + + @given(xys, sensible_floats, sensible_floats) def test_xy(xy, a, b): ab = a, b r = np.sqrt(a ** 2 + b ** 2) phi = np.arctan2(b, a) pos = np.stack(([a], [b]), axis=1) + np.isclose (xy[0] , a, rtol=1e-4) + np.isclose (xy[1] , b, rtol=1e-4) np.isclose (xy.x , a, rtol=1e-4) np.isclose (xy.y , b, rtol=1e-4) np.isclose (xy.X , a, rtol=1e-4) @@ -91,6 +106,11 @@ def test_xy(xy, a, b): np.allclose(xy.pos, pos, rtol=1e-3, atol=1e-03) +def test_xy_getitem_raises(): + with raises(IndexError): + xy(0,0)[2] + + @given(text(min_size=1, max_size=10, alphabet=ascii_letters)) def test_NNN_generates_NN_for_every_attribute_name(name): nnn = NNN()