diff --git a/invisible_cities/cities/components.py b/invisible_cities/cities/components.py index 2a339e1598..4bcad93cc3 100644 --- a/invisible_cities/cities/components.py +++ b/invisible_cities/cities/components.py @@ -17,38 +17,41 @@ import pandas as pd import inspect - -from .. dataflow import dataflow as fl -from .. evm .ic_containers import SensorData -from .. evm .event_model import KrEvent -from .. evm .event_model import Hit -from .. evm .event_model import Cluster -from .. evm .event_model import HitCollection -from .. evm .event_model import MCInfo -from .. core.system_of_units_c import units -from .. core.exceptions import XYRecoFail -from .. core.exceptions import NoInputFiles -from .. core.exceptions import NoOutputFile -from .. core.exceptions import InvalidInputFileStructure -from .. core.configure import EventRange -from .. core.configure import event_range_help -from .. reco import calib_functions as cf -from .. reco import calib_sensors_functions as csf -from .. reco import peak_functions as pkf -from .. reco import pmaps_functions as pmf -from .. reco import hits_functions as hif -from .. reco.tbl_functions import get_mc_info -from .. reco.xy_algorithms import corona -from .. filters.s1s2_filter import S12Selector -from .. filters.s1s2_filter import pmap_filter -from .. database import load_db -from .. sierpe import blr -from .. io.pmaps_io import load_pmaps -from .. io. hits_io import load_hits -from .. io. dst_io import load_dst -from .. types.ic_types import xy -from .. types.ic_types import NN -from .. types.ic_types import NNN +from .. dataflow import dataflow as fl +from .. evm .ic_containers import SensorData +from .. evm .event_model import KrEvent +from .. evm .event_model import Hit +from .. evm .event_model import Cluster +from .. evm .event_model import HitCollection +from .. evm .event_model import MCInfo +from .. evm .pmaps import SiPMCharge +from .. core .system_of_units_c import units +from .. core .exceptions import XYRecoFail +from .. core .exceptions import NoInputFiles +from .. core .exceptions import NoOutputFile +from .. core .exceptions import InvalidInputFileStructure +from .. core .configure import EventRange +from .. core .configure import event_range_help +from .. core .random_sampling import NoiseSampler +from .. reco import calib_functions as cf +from .. reco import calib_sensors_functions as csf +from .. reco import peak_functions as pkf +from .. reco import pmaps_functions as pmf +from .. reco import hits_functions as hif +from .. reco .tbl_functions import get_mc_info +from .. reco .xy_algorithms import corona +from .. filters.s1s2_filter import S12Selector +from .. filters.s1s2_filter import pmap_filter +from .. database import load_db +from .. sierpe import blr +from .. io .pmaps_io import load_pmaps +from .. io .hits_io import hits_from_df +from .. io .dst_io import load_dst +from .. io .hits_io import load_hits +from .. io .dst_io import load_dst +from .. types .ic_types import xy +from .. types .ic_types import NN +from .. types .ic_types import NNN NoneType = type(None) @@ -284,8 +287,8 @@ def hits_and_kdst_from_files(paths: List[str]) -> Iterator[Dict[str,Union[HitCol """Reader of the files, yields HitsCollection, pandas DataFrame with kdst info, mc_info, run_number, event_number and timestamp""" for path in paths: try: - hits = load_hits(path) - kdst_df = load_dst (path, 'DST', 'Events') + hits_df = load_dst (path, 'RECO', 'Events') + kdst_df = load_dst (path, 'DST' , 'Events') except tb.exceptions.NoSuchNodeError: continue @@ -297,10 +300,11 @@ def hits_and_kdst_from_files(paths: List[str]) -> Iterator[Dict[str,Union[HitCol except (tb.exceptions.NoSuchNodeError, IndexError): continue - check_lengths(event_info, hits) + check_lengths(event_info, hits_df.event.unique()) for evtinfo in event_info: event_number, timestamp = evtinfo.fetch_all_fields() + hits = hits_from_df(hits_df.loc[hits_df.event == event_number]) yield dict(hits = hits[event_number], kdst = kdst_df.loc[kdst_df.event==event_number], mc=mc_info, run_number=run_number, event_number=event_number, timestamp=timestamp) # NB, the monte_carlo writer is different from the others: @@ -389,11 +393,14 @@ def compute_z_and_dt(t_s2, t_s1, drift_v): return z, dt -def build_pointlike_event(dbfile, run_number, drift_v, reco): - datasipm = load_db.DataSiPM(dbfile, run_number) - sipm_xs = datasipm.X.values - sipm_ys = datasipm.Y.values - sipm_xys = np.stack((sipm_xs, sipm_ys), axis=1) +def build_pointlike_event(dbfile, run_number, drift_v, + reco, charge_type = SiPMCharge.raw): + datasipm = load_db.DataSiPM(dbfile, run_number) + sipm_xs = datasipm.X.values + sipm_ys = datasipm.Y.values + sipm_xys = np.stack((sipm_xs, sipm_ys), axis=1) + + sipm_noise = NoiseSampler(dbfile, run_number).signal_to_noise def build_pointlike_event(pmap, selector_output, event_number, timestamp): evt = KrEvent(event_number, timestamp * 1e-3) @@ -420,7 +427,8 @@ def build_pointlike_event(pmap, selector_output, event_number, timestamp): evt.S2t.append(peak.time_at_max_energy) xys = sipm_xys[peak.sipms.ids ] - qs = peak.sipms.sum_over_times + qs = peak.sipm_charge_array(sipm_noise, charge_type, + single_point = True) try: clusters = reco(xys, qs) except XYRecoFail: @@ -450,12 +458,16 @@ def build_pointlike_event(pmap, selector_output, event_number, timestamp): return build_pointlike_event -def hit_builder(dbfile, run_number, drift_v, reco, rebin_slices, rebin_method): +def hit_builder(dbfile, run_number, drift_v, reco, + rebin_slices, rebin_method, + charge_type = SiPMCharge.raw): datasipm = load_db.DataSiPM(dbfile, run_number) sipm_xs = datasipm.X.values sipm_ys = datasipm.Y.values sipm_xys = np.stack((sipm_xs, sipm_ys), axis=1) + sipm_noise = NoiseSampler(dbfile, run_number).signal_to_noise + barycenter = partial(corona, all_sipms = datasipm, Qthr = 0 * units.pes, @@ -492,25 +504,30 @@ def build_hits(pmap, selector_output, event_number, timestamp): peak = pmf.rebin_peak(peak, rebin_slices, rebin_method) - xys = sipm_xys[peak.sipms.ids ] - qs = peak.sipms.sum_over_times + xys = sipm_xys[peak.sipms.ids] + qs = peak.sipm_charge_array(sipm_noise, charge_type, + single_point = True) try : cluster = barycenter(xys, qs)[0] except XYRecoFail: xy_peak = xy(NN, NN) else : xy_peak = xy(cluster.X, cluster.Y) - for slice_no, t_slice in enumerate(peak.times): + sipm_charge = peak.sipm_charge_array(sipm_noise , + charge_type , + single_point=False) + for slice_no, (t_slice, qs) in enumerate(zip(peak.times , + sipm_charge)): z_slice = (t_slice - s1_t) * units.ns * drift_v e_slice = peak.pmts.sum_over_sensors[slice_no] try: - xys = sipm_xys[peak.sipms.ids ] - qs = peak.sipms.time_slice(slice_no) + xys = sipm_xys[peak.sipms.ids] clusters = reco(xys, qs) es = hif.split_energy(e_slice, clusters) for c, e in zip(clusters, es): hit = Hit(peak_no, c, z_slice, e, xy_peak) hitc.hits.append(hit) except XYRecoFail: - hit = Hit(peak_no, empty_cluster(), z_slice, e_slice, xy_peak) + hit = Hit(peak_no, empty_cluster(), z_slice, + e_slice, xy_peak) hitc.hits.append(hit) return hitc diff --git a/invisible_cities/cities/diomira.py b/invisible_cities/cities/diomira.py index 07fc161cb2..8e29ac247a 100644 --- a/invisible_cities/cities/diomira.py +++ b/invisible_cities/cities/diomira.py @@ -25,6 +25,7 @@ """ from functools import partial +from itertools import compress import numpy as np import tables as tb @@ -197,6 +198,9 @@ def do_nothing(*args, **kwargs): max = s2_params["s2_lmax" ]), rebin_stride = s2_params["s2_rebin_stride"]) + for channel in channels: + if channel not in compress(datapmt.ChannelID, datapmt.Active.values): + raise ValueError("Cannot trigger on a masked PMT") deconvolver = deconv_pmt(detector_db, run_number, n_baseline, IC_ids) diff --git a/invisible_cities/cities/diomira_test.py b/invisible_cities/cities/diomira_test.py index 6e561105f7..3a219f77c7 100644 --- a/invisible_cities/cities/diomira_test.py +++ b/invisible_cities/cities/diomira_test.py @@ -16,38 +16,6 @@ from . diomira import diomira -@mark.skip(reason="Table not written in liquid cities anymore") -def test_diomira_fee_table(ICDATADIR): - "Test that FEE table reads back correctly with expected values." - RWF_file = os.path.join(ICDATADIR, 'electrons_40keV_z250_RWF.h5') - - with tb.open_file(RWF_file, 'r') as e40rwf: - fee = tbl.read_FEE_table(e40rwf.root.FEE.FEE) - feep = fee.fee_param - eps = 1e-04 - # Ignoring PEP8 to improve readability by making symmetry explicit. - assert len(fee.adc_to_pes) == e40rwf.root.RD.pmtrwf.shape[1] - assert len(fee.coeff_blr) == e40rwf.root.RD.pmtrwf.shape[1] - assert len(fee.coeff_c) == e40rwf.root.RD.pmtrwf.shape[1] - assert len(fee.pmt_noise_rms) == e40rwf.root.RD.pmtrwf.shape[1] - assert feep.NBITS == FEE.NBITS - assert abs(feep.FEE_GAIN - FEE.FEE_GAIN) < eps - assert abs(feep.LSB - FEE.LSB) < eps - assert abs(feep.NOISE_I - FEE.NOISE_I) < eps - assert abs(feep.NOISE_DAQ - FEE.NOISE_DAQ) < eps - assert abs(feep.C2/units.nF - FEE.C2/units.nF) < eps - assert abs(feep.C1/units.nF - FEE.C1/units.nF) < eps - assert abs(feep.R1/units.ohm - FEE.R1/units.ohm) < eps - assert abs(feep.ZIN/units.ohm - FEE.Zin/units.ohm) < eps - assert abs(feep.t_sample - FEE.t_sample) < eps - assert abs(feep.f_sample - FEE.f_sample) < eps - assert abs(feep.f_mc - FEE.f_mc) < eps - assert abs(feep.f_LPF1 - FEE.f_LPF1) < eps - assert abs(feep.f_LPF2 - FEE.f_LPF2) < eps - assert abs(feep.OFFSET - FEE.OFFSET) < eps - assert abs(feep.CEILING - FEE.CEILING) < eps - - def test_diomira_identify_bug(ICDATADIR): """Read a one-event file in which the energy of PMTs is equal to zero and asset it must be son. This test would fail for a normal file where there @@ -112,16 +80,6 @@ def test_diomira_copy_mc_and_offset(ICDATADIR, config_tmpdir): last_evt_number = h5out.root.MC.extents[-1][0] assert first_evt_number != last_evt_number - -@mark.skip(reason="Feature removed") -@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), - ('electrons_40keV_z250_MCRD.h5' , 0))) -def test_event_number_from_input_file_name(filename, first_evt): - assert Diomira.event_number_from_input_file_name(filename) == first_evt - - @mark.slow def test_diomira_mismatch_between_input_and_database(ICDATADIR, output_tmpdir): file_in = os.path.join(ICDATADIR , 'electrons_40keV_z250_MCRD.h5') @@ -140,18 +98,26 @@ def test_diomira_mismatch_between_input_and_database(ICDATADIR, output_tmpdir): assert cnt.events_in == 1 -@mark.skip(reason="Trigger not implemented in liquid cities") @mark.slow def test_diomira_trigger_on_masked_pmt_raises_ValueError(ICDATADIR, output_tmpdir): file_in = os.path.join(ICDATADIR , 'electrons_40keV_z250_MCRD.h5') file_out = os.path.join(output_tmpdir, 'electrons_40keV_z250_RWF_test_trigger.h5') conf = configure('diomira invisible_cities/config/diomira.conf'.split()) - conf.update(dict(run_number = -4500, # Must be a run number with dead pmts - files_in = file_in, - file_out = file_out, - trigger_type = "S2", - tr_channels = (0,), # This is a masked PMT for this run + conf.update(dict(run_number = -4500, # Must be a run number with dead pmts + files_in = file_in, + file_out = file_out, + trigger_type = "S2", + trigger_params = dict( + tr_channels = (0,), # This is a masked PMT for this run + min_number_channels = 1 , + data_mc_ratio = 1 , + min_height = 0 , + max_height = 1e6 , + min_charge = 0 , + max_charge = 1e6 , + min_width = 0 , + max_width = 1e6 ), event_range = (0, 1))) with raises(ValueError): diff --git a/invisible_cities/cities/dorothea_test.py b/invisible_cities/cities/dorothea_test.py index 7b071f5e36..c95b6d0d9e 100644 --- a/invisible_cities/cities/dorothea_test.py +++ b/invisible_cities/cities/dorothea_test.py @@ -103,41 +103,6 @@ def test_dorothea_filter_events(config_tmpdir, Kr_pmaps_run4628_filename): assert np.all(dst.s2_peak.values == s2_peak_pass) - -@mark.skip(reason="The configuration of corona in dorothea can only return 1 cluster." - "The assertion of more than one cluster has been removed") -def test_dorothea_issue_347(Kr_pmaps_run4628_filename, config_tmpdir): - PATH_IN = Kr_pmaps_run4628_filename - PATH_OUT = os.path.join(config_tmpdir, 'KrDST.h5') - conf = configure('dummy invisible_cities/config/dorothea_with_corona.conf'.split()) - - # with this parameters Corona will find several clusters - conf.update(dict(run_number = 4628, - files_in = PATH_IN, - file_out = PATH_OUT, - lm_radius = 10.0, - new_lm_radius = 13.0, - msipm = 1)) - cnt = dorothea(**conf) - assert cnt.n_events_more_than_1_cluster == 3 - - -@mark.skip("This scenario is not possible in liquid cities") -def test_dorothea_event_not_found(ICDATADIR, output_tmpdir): - file_in = os.path.join(ICDATADIR , "kr_rwf_0_0_7bar_NEXT_v1_00_05_v0.9.2_20171011_krmc_irene_3evt.h5") - file_out = os.path.join(output_tmpdir, "test_dorothea_event_not_found.h5") - - conf = configure('dummy invisible_cities/config/dorothea.conf'.split()) - nevt = 3 - - conf.update(dict(files_in = file_in, - file_out = file_out, - event_range = (0, nevt))) - - cnt = dorothea(**conf) - assert cnt.n_empty_pmaps == 1 - - def test_dorothea_exact_result(ICDATADIR, output_tmpdir): file_in = os.path.join(ICDATADIR , "Kr83_nexus_v5_03_00_ACTIVE_7bar_3evts.PMP.h5") file_out = os.path.join(output_tmpdir, "exact_result_dorothea.h5") diff --git a/invisible_cities/cities/esmeralda.py b/invisible_cities/cities/esmeralda.py index 90a64750ff..66478da235 100644 --- a/invisible_cities/cities/esmeralda.py +++ b/invisible_cities/cities/esmeralda.py @@ -6,24 +6,25 @@ This city corrects hits energy and extracts topology information. The input is penthesilea output containing hits, kdst global information and mc info. The city outputs : - - RECO corrected hits table, according to the given map and hits threshold and merging parameters + - CHITS corrected hits table, + - lowTh - contains corrected hits table that passed h.Q >= charge_threshold_low constrain + - highTh - contains corrected hits table that passed h.Q >= charge_threshold_high constrain. + it also contains: + - Ep field that is the energy of a hit after applying drop_end_point_voxel algorithm. + - track_id denoting to which track from Tracking/Tracks dataframe the hit belong to - MC info (if run number <=0) - - PAOLINA Summary of topological anlaysis containing: - - corrected (paolina) hits, hits table after performing paolina drop_voxels - - summary of per track information - - summary of per event information + - Tracking/Tracks - summary of per track information + - Summary/events - summary of per event information + - DST/Events - copy of kdst information from penthesilea """ import os import tables as tb import numpy as np import pandas as pd -import warnings -from functools import partial -from typing import Tuple +from collections import OrderedDict from typing import Callable -from typing import Optional from .. reco import tbl_functions as tbl from .. reco import paolina_functions as plf @@ -38,7 +39,6 @@ from . components import print_every from . components import hits_and_kdst_from_files -from .. types. ic_types import NN from .. types. ic_types import xy from .. io. hits_io import hits_writer from .. io. mcinfo_io import mc_info_writer @@ -46,11 +46,31 @@ from .. io. event_filter_io import event_filter_writer from .. io. dst_io import _store_pandas_as_tables +types_dict_summary = OrderedDict({'event' : np.int32 , 'evt_energy' : np.float64, 'evt_charge' : np.float64, + 'evt_ntrks' : np.int , 'evt_nhits' : np.int , 'evt_x_avg' : np.float64, + 'evt_y_avg' : np.float64, 'evt_z_avg' : np.float64, 'evt_r_avg' : np.float64, + 'evt_x_min' : np.float64, 'evt_y_min' : np.float64, 'evt_z_min' : np.float64, + 'evt_r_min' : np.float64, 'evt_x_max' : np.float64, 'evt_y_max' : np.float64, + 'evt_z_max' : np.float64, 'evt_r_max' : np.float64, 'evt_out_of_map': bool }) + +types_dict_tracks = OrderedDict({'event' : np.int32 , 'trackID' : np.int , 'energy' : np.float64, + 'length' : np.float64, 'numb_of_voxels': np.int , 'numb_of_hits': np.int , + 'numb_of_tracks' : np.int , 'x_min' : np.float64, 'y_min' : np.float64, + 'z_min' : np.float64, 'r_min' : np.float64, 'x_max' : np.float64, + 'y_max' : np.float64, 'z_max' : np.float64, 'r_max' : np.float64, + 'x_ave' : np.float64, 'y_ave' : np.float64, 'z_ave' : np.float64, + 'r_ave' : np.float64, 'extreme1_x' : np.float64, 'extreme1_y' : np.float64, + 'extreme1_z' : np.float64, 'extreme2_x' : np.float64, 'extreme2_y' : np.float64, + 'extreme2_z' : np.float64, 'blob1_x' : np.float64, 'blob1_y' : np.float64, + 'blob1_z' : np.float64, 'blob2_x' : np.float64, 'blob2_y' : np.float64, + 'blob2_z' : np.float64, 'eblob1' : np.float64, 'eblob2' : np.float64, + 'ovlp_blob_energy': np.float64, + 'vox_size_x' : np.float64, 'vox_size_y' : np.float64, 'vox_size_z' : np.float64}) + def hits_threshold_and_corrector(map_fname : str , threshold_charge : float, same_peak : bool , - apply_temp : bool , - norm_strat : cof.norm_strategy + apply_temp : bool ) -> Callable: """ For a given correction map and hit threshold/ merging parameters returns a function that applies thresholding, merging and @@ -83,7 +103,7 @@ class norm_strategy(AutoNameEnumBase): """ map_fname = os.path.expandvars(map_fname) maps = cof.read_maps(map_fname) - get_coef = cof.apply_all_correction(maps, apply_temp = apply_temp, norm_strat = norm_strat) + get_coef = cof.apply_all_correction(maps, apply_temp = apply_temp, norm_strat = cof.norm_strategy.kr) if maps.t_evol is not None: time_to_Z = cof.get_df_to_z_converter(maps) else: @@ -92,8 +112,6 @@ def threshold_and_correct_hits(hitc : evm.HitCollection) -> evm.HitCollection: t = hitc.time thr_hits = hif.threshold_hits(hitc.hits, threshold_charge ) mrg_hits = hif.merge_NN_hits ( thr_hits, same_peak = same_peak) - if len(mrg_hits) == 0: - return None X = np.fromiter((h.X for h in mrg_hits), float) Y = np.fromiter((h.Y for h in mrg_hits), float) Z = np.fromiter((h.Z for h in mrg_hits), float) @@ -110,30 +128,25 @@ def threshold_and_correct_hits(hitc : evm.HitCollection) -> evm.HitCollection: return new_hitc return threshold_and_correct_hits - -def events_filter(allow_nans : bool) -> Callable: +def copy_Ec_to_Ep_hit_attribute_(hitc : evm.HitCollection) -> evm.HitCollection: """ - Filter for Esmeralda flow. The flow stops if: - 1. there are no hits - 2. there are hits with NaN energy and allow_nans is set to False + The functions copies values of Ec attributes into Ep attributes. Takes as input HitCollection and returns a copy. """ - def filter_events(hitc : Optional[evm.HitCollection]) -> bool: - if hitc == None: - return False - elif allow_nans == False: - nans = np.isnan([h.Ec for h in hitc.hits]) - return not(any(nans)) - else: - return True - return filter_events + mod_hits = [] + for hit in hitc.hits: + hit = evm.Hit(hit.npeak, evm.Cluster(hit.Q, xy(hit.X, hit.Y), hit.var, hit.nsipm), + hit.Z, hit.E, xy(hit.Xpeak, hit.Ypeak), s2_energy_c=hit.Ec, Ep=hit.Ec) + mod_hits.append(hit) + mod_hitc = evm.HitCollection(hitc.event, hitc.time, hits=mod_hits) + return mod_hitc def track_blob_info_creator_extractor(vox_size : [float, float, float], - energy_type : evm.HitEnergy , strict_vox_size : bool , energy_threshold : float , min_voxels : int , - blob_radius : float + blob_radius : float , + max_num_hits : int ) -> Callable: """ For a given paolina parameters returns a function that extract tracks / blob information from a HitCollection. @@ -142,12 +155,6 @@ def track_blob_info_creator_extractor(vox_size : [float, float, float], ---------- vox_size : [float, float, float] (maximum) size of voxels for track reconstruction - energy_type : HitEnergy - class HitEnergy(AutoNameEnumBase): - E = auto() - Ec = auto() - energy attribute to use for voxelization/ tracks - strict_vox_size : bool if False allows per event adaptive voxel size, smaller of equal thatn vox_size @@ -164,99 +171,89 @@ class HitEnergy(AutoNameEnumBase): A function that from a given HitCollection returns a pandas DataFrame with per track information. """ def create_extract_track_blob_info(hitc): - voxels = plf.voxelize_hits(hitc.hits, vox_size, strict_vox_size, energy_type) - mod_voxels = plf.drop_end_point_voxels(voxels, energy_threshold, min_voxels) - tracks = plf.make_track_graphs(mod_voxels) - - vox_size_x = voxels[0].size[0] - vox_size_y = voxels[0].size[1] - vox_size_z = voxels[0].size[2] - - def get_track_energy(track): - return sum([vox.Ehits for vox in track.nodes()]) - #sort tracks in energy - tracks = sorted(tracks, key = get_track_energy, reverse = True) - - track_hits = [] - df = pd.DataFrame(columns=['event', 'trackID', 'energy', 'length', 'numb_of_voxels', - 'numb_of_hits', 'numb_of_tracks', 'x_min', 'y_min', 'z_min', - 'x_max', 'y_max', 'z_max', 'r_max', 'x_ave', 'y_ave', 'z_ave', - 'extreme1_x', 'extreme1_y', 'extreme1_z', - 'extreme2_x', 'extreme2_y', 'extreme2_z', - 'blob1_x', 'blob1_y', 'blob1_z', - 'blob2_x', 'blob2_y', 'blob2_z', - 'eblob1', 'eblob2', 'ovlp_blob_energy', - 'vox_size_x', 'vox_size_y', 'vox_size_z']) - for c, t in enumerate(tracks, 0): - tID = c - energy = get_track_energy(t) - length = plf.length(t) - numb_of_hits = len([h for vox in t.nodes() for h in vox.hits]) - numb_of_voxels = len(t.nodes()) - numb_of_tracks = len(tracks ) - - min_x = min([h.X for v in t.nodes() for h in v.hits]) - max_x = max([h.X for v in t.nodes() for h in v.hits]) - min_y = min([h.Y for v in t.nodes() for h in v.hits]) - max_y = max([h.Y for v in t.nodes() for h in v.hits]) - min_z = min([h.Z for v in t.nodes() for h in v.hits]) - max_z = max([h.Z for v in t.nodes() for h in v.hits]) - max_r = max([np.sqrt(h.X*h.X + h.Y*h.Y) for v in t.nodes() for h in v.hits]) - - pos = [h.pos for v in t.nodes() for h in v.hits] - e = [getattr(h, energy_type.value) for v in t.nodes() for h in v.hits] - ave_pos = np.average(pos, weights=e, axis=0) - - extr1, extr2 = plf.find_extrema(t) - extr1_pos = extr1.XYZ - extr2_pos = extr2.XYZ - - blob_pos1, blob_pos2 = plf.blob_centres(t, blob_radius) - - e_blob1, e_blob2, hits_blob1, hits_blob2 = plf.blob_energies_and_hits(t, blob_radius) - overlap = False - overlap = sum([h.Ec for h in set(hits_blob1).intersection(hits_blob2)]) - list_of_vars = [hitc.event, tID, energy, length, numb_of_voxels, - numb_of_hits, numb_of_tracks, - min_x, min_y, min_z, max_x, max_y, max_z, max_r, - ave_pos[0], ave_pos[1], ave_pos[2], - extr1_pos[0], extr1_pos[1], extr1_pos[2], - extr2_pos[0], extr2_pos[1], extr2_pos[2], - blob_pos1[0], blob_pos1[1], blob_pos1[2], - blob_pos2[0], blob_pos2[1], blob_pos2[2], - e_blob1, e_blob2, overlap, - vox_size_x, vox_size_y, vox_size_z] - - df.loc[c] = list_of_vars - try: - types_dict - except NameError: - types_dict = dict(zip(df.columns, [type(x) for x in list_of_vars])) - - for vox in t.nodes(): - for hit in vox.hits: - hit.track_id = tID - track_hits.append(hit) - - + df = pd.DataFrame(columns=list(types_dict_tracks.keys())) + 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 = evm.HitCollection(hitc.event, hitc.time) - track_hitc.hits = track_hits - #change dtype of columns to match type of variables - df = df.apply(lambda x : x.astype(types_dict[x.name])) - - return df, track_hitc + 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 = evm.HitCollection(hitc.event, hitc.time) + hitc.hits = hits_without_nan + + if len(hitc.hits) > 0: + voxels = plf.voxelize_hits(hitc.hits, vox_size, strict_vox_size, evm.HitEnergy.Ep) + ( mod_voxels, + dropped_voxels) = plf.drop_end_point_voxels(voxels, energy_threshold, min_voxels) + tracks = plf.make_track_graphs(mod_voxels) + + for v in dropped_voxels: + track_hitc.hits.extend(v.hits) + + vox_size_x = voxels[0].size[0] + vox_size_y = voxels[0].size[1] + vox_size_z = voxels[0].size[2] + del(voxels) + #sort tracks in energy + tracks = sorted(tracks, key=plf.get_track_energy, reverse=True) + + track_hits = [] + for c, t in enumerate(tracks, 0): + tID = c + energy = plf.get_track_energy(t) + length = plf.length(t) + numb_of_hits = len([h for vox in t.nodes() for h in vox.hits]) + numb_of_voxels = len(t.nodes()) + numb_of_tracks = len(tracks ) + pos = [h.pos for v in t.nodes() for h in v.hits] + x, y, z = map(np.array, zip(*pos)) + r = np.sqrt(x**2 + y**2) + + e = [h.Ep for v in t.nodes() for h in v.hits] + ave_pos = np.average(pos, weights=e, axis=0) + ave_r = np.average(r , weights=e, axis=0) + extr1, extr2 = plf.find_extrema(t) + extr1_pos = extr1.XYZ + extr2_pos = extr2.XYZ + + blob_pos1, blob_pos2 = plf.blob_centres(t, blob_radius) + + e_blob1, e_blob2, hits_blob1, hits_blob2 = plf.blob_energies_and_hits(t, blob_radius) + overlap = float(sum(h.Ep for h in set(hits_blob1).intersection(set(hits_blob2)))) + list_of_vars = [hitc.event, tID, energy, length, numb_of_voxels, + numb_of_hits, numb_of_tracks, + min(x), min(y), min(z), min(r), max(x), max(y), max(z), max(r), + *ave_pos, ave_r, *extr1_pos, + *extr2_pos, *blob_pos1, *blob_pos2, + e_blob1, e_blob2, overlap, + vox_size_x, vox_size_y, vox_size_z] + + df.loc[c] = list_of_vars + + for vox in t.nodes(): + for hit in vox.hits: + hit.track_id = tID + track_hits.append(hit) + + #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 create_extract_track_blob_info def make_event_summary(event_number : int , - timestamp : int , topology_info : pd.DataFrame , paolina_hits : evm.HitCollection, - kdst : pd.DataFrame + out_of_map : bool ) -> pd.DataFrame: """ - For a given event number, timestamp, topology info dataframe, paolina hits and kdst information returns a + For a given event number, timestamp, topology info dataframe, paolina hits and kdst information returns a dataframe with the whole event summary. Parameters @@ -277,78 +274,59 @@ def make_event_summary(event_number : int , ---------- DataFrame containing relevant per event information. """ - es = pd.DataFrame(columns=['event', 'time', 'S1e', 'S1t', - 'nS2', 'ntrks', 'nhits', 'S2e0', 'S2ec', - 'S2q0', 'S2qc', 'x_avg', 'y_avg', 'z_avg', - 'r_avg', 'x_min', 'y_min', 'z_min', 'r_min', - 'x_max', 'y_max', 'z_max', 'r_max']) - - if(len(kdst.s1_peak.unique()) != 1): - warnings.warn("Number of recorded S1 energies differs from 1 in event {}.Choosing first S1".format(event_number), UserWarning) - - S1e = kdst.S1e.values[0] - S1t = kdst.S1t.values[0] - - nS2 = kdst.nS2.values[0] + es = pd.DataFrame(columns=list(types_dict_summary.keys())) ntrks = len(topology_info.index) nhits = len(paolina_hits.hits) - S2e0 = np.sum(kdst.S2e.values) - S2ec = max(-1, sum([h.Ec for h in paolina_hits.hits])) - - S2q0 = np.sum(kdst.S2q.values) - S2qc = max(-1, sum([h.Qc for h in paolina_hits.hits])) + S2ec = sum(h.Ec for h in paolina_hits.hits) + S2qc = -1 #not implemented yet - x_avg = sum([h.X*h.Ec for h in paolina_hits.hits]) - y_avg = sum([h.Y*h.Ec for h in paolina_hits.hits]) - z_avg = sum([h.Z*h.Ec for h in paolina_hits.hits]) - r_avg = sum([(h.X**2 + h.Y**2)**0.5*h.Ec for h in paolina_hits.hits]) - if(S2ec > 0): - x_avg /= S2ec - y_avg /= S2ec - z_avg /= S2ec - r_avg /= S2ec + pos = [h.pos for h in paolina_hits.hits] + x, y, z = map(np.array, zip(*pos)) + r = np.sqrt(x**2 + y**2) - x_min = min([h.X for h in paolina_hits.hits]) - y_min = min([h.Y for h in paolina_hits.hits]) - z_min = min([h.Z for h in paolina_hits.hits]) - r_min = min([(h.X**2 + h.Y**2)**0.5 for h in paolina_hits.hits]) + e = [h.Ec for h in paolina_hits.hits] + ave_pos = np.average(pos, weights=e, axis=0) + ave_r = np.average(r , weights=e, axis=0) - x_max = max([h.X for h in paolina_hits.hits]) - y_max = max([h.Y for h in paolina_hits.hits]) - z_max = max([h.Z for h in paolina_hits.hits]) - r_max = max([(h.X**2 + h.Y**2)**0.5 for h in paolina_hits.hits]) - list_of_vars = [event_number, int(timestamp), S1e, S1t, nS2, ntrks, nhits, S2e0, - S2ec, S2q0, S2qc, x_avg, y_avg, z_avg, r_avg, x_min, y_min, - z_min, r_min, x_max, y_max, z_max, r_max] + 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] es.loc[0] = list_of_vars #change dtype of columns to match type of variables - types_dict = dict(zip(es.columns, [type(x) for x in list_of_vars])) - es = es.apply(lambda x : x.astype(types_dict[x.name])) - + es = es.apply(lambda x : x.astype(types_dict_summary[x.name])) return es -def track_writer(h5out, compression='ZLIB4', group_name='PAOLINA', table_name='Tracks', descriptive_string='Track information', str_col_length=32): +def track_writer(h5out, compression='ZLIB4', group_name='Tracking', table_name='Tracks', descriptive_string='Track information', str_col_length=32): """ For a given open table returns a writer for topology info dataframe """ def write_tracks(df): - return _store_pandas_as_tables(h5out = h5out, df = df, compression = compression, group_name = group_name, table_name = table_name, descriptive_string = descriptive_string, str_col_length = str_col_length) + return _store_pandas_as_tables(h5out=h5out, df=df, compression=compression, group_name=group_name, table_name=table_name, descriptive_string=descriptive_string, str_col_length=str_col_length) return write_tracks -def summary_writer(h5out, compression='ZLIB4', group_name='PAOLINA', table_name='Summary', descriptive_string='Event summary information', str_col_length=32): +def summary_writer(h5out, compression='ZLIB4', group_name='Summary', table_name='Events', descriptive_string='Event summary information', str_col_length=32): """ For a given open table returns a writer for summary info dataframe """ def write_summary(df): - return _store_pandas_as_tables(h5out = h5out, df = df, compression = compression, group_name = group_name, table_name = table_name, descriptive_string = descriptive_string, str_col_length = str_col_length) + return _store_pandas_as_tables(h5out=h5out, df=df, compression=compression, group_name=group_name, table_name=table_name, descriptive_string=descriptive_string, str_col_length=str_col_length) return write_summary +def kdst_from_df_writer(h5out, compression='ZLIB4', group_name='DST', table_name='Events', descriptive_string='KDST Events', str_col_length=32): + """ + For a given open table returns a writer for KDST dataframe info + """ + def write_kdst(df): + return _store_pandas_as_tables(h5out=h5out, df=df, compression=compression, group_name=group_name, table_name=table_name, descriptive_string=descriptive_string, str_col_length=str_col_length) + return write_kdst + @city def esmeralda(files_in, file_out, compression, event_range, print_mod, run_number, @@ -375,23 +353,19 @@ def esmeralda(files_in, file_out, compression, event_range, print_mod, run_numbe cor_hits_params : dict map_fname : string (filepath) filename of the map - threshold_charge_NN : float - minimum pes for a RECO hit - threshold_charge_paolina : float - minimum pes for a PAOLINA hit + threshold_charge_low : float + minimum pes for a lowTh hit + threshold_charge_high : float + minimum pes for a highTh hit same_peak : bool if True energy of NN hits is assigned only to the hits from the same peak apply_temp : bool whether to apply temporal corrections must be set to False if no temporal correction dataframe exists in map file - norm_strat : 'max', 'mean' or 'kr' - strategy to normalize the energy paolina_params :dict vox_size : [float, float, float] (maximum) size of voxels for track reconstruction - energy_type : 'corrected' ot 'uncorrected' - energy attribute to use for voxelization/ tracks strict_vox_size : bool if False allows per event adaptive voxel size, smaller of equal thatn vox_size @@ -402,6 +376,8 @@ def esmeralda(files_in, file_out, compression, event_range, print_mod, run_numbe after min_voxel number of voxels is reached no dropping will happen. blob_radius : float radius of blob + max_num_hits : int + maximum number of hits allowed per event to run paolina functions. ---------- Input ---------- @@ -409,91 +385,103 @@ def esmeralda(files_in, file_out, compression, event_range, print_mod, run_numbe ---------- Output ---------- - RECO : corrected hits table, according to the given map and hits threshold and merging parameters - MC info : (if run number <=0) - PAOLINA : Summary of topological anlaysis containing: - Events : corrected (paolina) hits, hits table after performing paolina drop_voxels - Tracks : summary of per track topological information - Summary : summary of per event information + - CHITS corrected hits table, + - lowTh - contains corrected hits table that passed h.Q >= charge_threshold_low constrain + - highTh - contains corrected hits table that passed h.Q >= charge_threshold_high constrain. + it also contains: + - Ep field that is the energy of a hit after applying drop_end_point_voxel algorithm. + - track_id denoting to which track from Tracking/Tracks dataframe the hit belong to + - MC info (if run number <=0) + - Tracking/Tracks - summary of per track information + - Summary/events - summary of per event information + - DST/Events - copy of kdst information from penthesilea """ - if paolina_params ['energy_type'] == 'corrected' : energy_type = evm.HitEnergy.Ec - elif paolina_params ['energy_type'] == 'uncorrected' : energy_type = evm.HitEnergy.E - else : raise ValueError(f"Unrecognized processing mode: {paolina_params['energy_type']}") - if cor_hits_params['norm_strat'] == 'max' : norm_strat = cof.norm_strategy.max - elif cor_hits_params['norm_strat'] == 'mean' : norm_strat = cof.norm_strategy.mean - elif cor_hits_params['norm_strat'] == 'kr' : norm_strat = cof.norm_strategy.kr - else : raise ValueError(f"Unrecognized processing mode: {cor_hits_params['norm_strat']}") + cor_hits_params_ = {value : cor_hits_params.get(value) for value in ['map_fname', 'same_peak', 'apply_temp']} - cor_hits_params_ = {value : cor_hits_params.get(value) for value in ['map_fname', 'same_peak' , 'apply_temp' ]} - paolina_params_ = {value : paolina_params .get(value) for value in ['vox_size' , 'strict_vox_size', 'energy_threshold', 'min_voxels', 'blob_radius']} + threshold_and_correct_hits_low = fl.map(hits_threshold_and_corrector(threshold_charge=cor_hits_params['threshold_charge_low' ], **cor_hits_params_), + args = 'hits', + out = 'cor_low_th_hits') - threshold_and_correct_hits_NN = fl.map(hits_threshold_and_corrector(**cor_hits_params_, threshold_charge = cor_hits_params['threshold_charge_NN' ], norm_strat = norm_strat), - args = 'hits', - out = 'NN_hits') + threshold_and_correct_hits_high = fl.map(hits_threshold_and_corrector(threshold_charge=cor_hits_params['threshold_charge_high'], **cor_hits_params_), + args = 'hits', + out = 'cor_high_th_hits') - threshold_and_correct_hits_paolina = fl.map(hits_threshold_and_corrector(**cor_hits_params_, threshold_charge = cor_hits_params['threshold_charge_paolina'], norm_strat = norm_strat), - args = 'hits', - out = 'corrected_hits') + filter_events_low_th = fl.map(lambda x : len(x.hits) > 0, + args = 'cor_low_th_hits', + out = 'low_th_hits_passed') - filter_events_NN = fl.map(events_filter(allow_nans = True), - args = 'NN_hits', - out = 'NN_hits_passed') + filter_events_high_th = fl.map(lambda x : len(x.hits) > 0, + args = 'cor_high_th_hits', + out = 'high_th_hits_passed') - filter_events_paolina = fl.map(events_filter(allow_nans = False), - args = 'corrected_hits', - out = 'paolina_hits_passed') + hits_passed_low_th = fl.count_filter(bool, args="low_th_hits_passed") + hits_passed_high_th = fl.count_filter(bool, args="high_th_hits_passed") - hits_passed_NN = fl.count_filter(bool, args = "NN_hits_passed") - hits_passed_paolina = fl.count_filter(bool, args = "paolina_hits_passed") + copy_Ec_to_Ep_hit_attribute = fl.map(copy_Ec_to_Ep_hit_attribute_, + args = 'cor_high_th_hits', + out = 'cor_Ep_high_th_hits') - create_extract_track_blob_info = fl.map(track_blob_info_creator_extractor(**paolina_params_, energy_type = energy_type), - args = 'corrected_hits', - out = ('topology_info', 'paolina_hits')) + create_extract_track_blob_info = fl.map(track_blob_info_creator_extractor(**paolina_params), + args = 'cor_Ep_high_th_hits', + out = ('topology_info', 'paolina_hits', 'out_of_map')) + filter_events_topology = fl.map(lambda x : len(x) > 0, + args = 'topology_info', + out = 'topology_passed') + events_passed_topology = fl.count_filter(bool, args="topology_passed") - make_final_summary = fl.map(make_event_summary, - args = ('event_number', 'timestamp', 'topology_info', 'paolina_hits', 'kdst'), - out = 'event_info') + make_final_summary = fl.map(make_event_summary, + args = ('event_number', 'topology_info', 'paolina_hits', 'out_of_map'), + out = 'event_info') event_count_in = fl.spy_count() event_count_out = fl.spy_count() - with tb.open_file(file_out, "w", filters = tbl.filters(compression)) as h5out: + with tb.open_file(file_out, "w", filters=tbl.filters(compression)) as h5out: # Define writers... write_event_info = fl.sink(run_and_event_writer(h5out), args=("run_number", "event_number", "timestamp")) write_mc_ = mc_info_writer(h5out) if run_number <= 0 else (lambda *_: None) - write_mc = fl.sink( write_mc_ , args = ("mc", "event_number")) - write_hits_NN = fl.sink( hits_writer (h5out) , args = "NN_hits" ) - write_hits_paolina = fl.sink( hits_writer (h5out, group_name = 'PAOLINA'), args = "paolina_hits" ) - write_tracks = fl.sink( track_writer (h5out=h5out) , args = "topology_info" ) - write_summary = fl.sink( summary_writer (h5out=h5out) , args = "event_info" ) - write_paolina_filter = fl.sink( event_filter_writer(h5out, "paolina_select") , args = ("event_number", "paolina_hits_passed")) - write_NN_filter = fl.sink( event_filter_writer(h5out, "NN_select") , args = ("event_number", "NN_hits_passed")) - - return push(source = hits_and_kdst_from_files(files_in), - pipe = pipe( - fl.slice(*event_range, close_all=True) , - print_every(print_mod) , - event_count_in .spy , - fl.branch(threshold_and_correct_hits_NN , - filter_events_NN , - fl.branch(write_NN_filter) , - hits_passed_NN .filter , - fl.fork(write_mc , - write_hits_NN) ), - threshold_and_correct_hits_paolina , - filter_events_paolina , - fl.branch(write_paolina_filter) , - hits_passed_paolina .filter , - event_count_out .spy , - create_extract_track_blob_info , - fl.fork(write_hits_paolina , - write_tracks , - (make_final_summary, write_summary), - write_event_info)) , - result = dict(events_in = event_count_in .future, - events_out = event_count_out.future)) + write_hits_low_th = fl.sink( hits_writer (h5out, group_name='CHITS', table_name='lowTh'), + args="cor_low_th_hits") + write_hits_paolina = fl.sink( hits_writer (h5out, group_name='CHITS', table_name='highTh' ), + args="paolina_hits" ) + write_mc = fl.sink( write_mc_ , args=("mc", "event_number")) + write_tracks = fl.sink( track_writer (h5out=h5out) , args="topology_info" ) + write_summary = fl.sink( summary_writer (h5out=h5out) , args="event_info" ) + write_high_th_filter = fl.sink( event_filter_writer(h5out, "high_th_select" ) , args=("event_number", "high_th_hits_passed")) + write_low_th_filter = fl.sink( event_filter_writer(h5out, "low_th_select" ) , args=("event_number", "low_th_hits_passed" )) + write_topology_filter = fl.sink( event_filter_writer(h5out, "topology_select") , args=("event_number", "topology_passed" )) + write_kdst_table = fl.sink( kdst_from_df_writer(h5out) , args="kdst" ) + + return push(source = hits_and_kdst_from_files(files_in), + pipe = pipe( + fl.slice(*event_range, close_all=True) , + print_every(print_mod) , + event_count_in .spy , + fl.branch(fl.fork(write_kdst_table , + write_event_info , + write_mc )), + fl.branch(threshold_and_correct_hits_low , + filter_events_low_th , + fl.branch(write_low_th_filter) , + hits_passed_low_th.filter , + write_hits_low_th ), + threshold_and_correct_hits_high , + filter_events_high_th , + fl.branch(write_high_th_filter) , + hits_passed_high_th .filter , + copy_Ec_to_Ep_hit_attribute , + create_extract_track_blob_info , + filter_events_topology , + fl.branch(make_final_summary, write_summary) , + fl.branch(write_topology_filter) , + fl.branch(write_hits_paolina) , + events_passed_topology.filter , + event_count_out .spy , + write_tracks ), + result = dict(events_in =event_count_in .future , + events_out=event_count_out.future )) diff --git a/invisible_cities/cities/esmeralda_test.py b/invisible_cities/cities/esmeralda_test.py index e9d1dbb560..c861e5e7d9 100644 --- a/invisible_cities/cities/esmeralda_test.py +++ b/invisible_cities/cities/esmeralda_test.py @@ -3,8 +3,8 @@ import tables as tb import pandas as pd -from pytest import mark - +from . components import get_event_info +from . components import length_of from .. core.system_of_units_c import units from .. core.configure import configure from .. core.configure import all as all_events @@ -13,22 +13,21 @@ from .. core.testing_utils import assert_dataframes_close from .. core.testing_utils import assert_tables_equality -@mark.serial + def test_esmeralda_contains_all_tables(KrMC_hdst_filename, correction_map_MC_filename, config_tmpdir): PATH_IN = KrMC_hdst_filename PATH_OUT = os.path.join(config_tmpdir, "Kr_tracks_with_MC.h5") conf = configure('dummy invisible_cities/config/esmeralda.conf'.split()) nevt_req = all_events - conf.update(dict(files_in = PATH_IN , - file_out = PATH_OUT , - event_range = nevt_req , - cor_hits_params = dict( - map_fname = correction_map_MC_filename, - threshold_charge_NN = 6 * units.pes , - threshold_charge_paolina = 30 * units.pes , - same_peak = True , - norm_strat = 'kr' , - apply_temp = False ))) + conf.update(dict(files_in = PATH_IN , + file_out = PATH_OUT , + event_range = nevt_req , + cor_hits_params = dict( + map_fname = correction_map_MC_filename, + threshold_charge_low = 6 * units.pes , + threshold_charge_high = 30 * units.pes , + same_peak = True , + apply_temp = False ))) esmeralda(**conf) with tb.open_file(PATH_OUT) as h5out: @@ -36,171 +35,343 @@ def test_esmeralda_contains_all_tables(KrMC_hdst_filename, correction_map_MC_fil assert "MC/extents" in h5out.root assert "MC/hits" in h5out.root assert "MC/particles" in h5out.root - assert "PAOLINA" in h5out.root - assert "PAOLINA/Events" in h5out.root - assert "PAOLINA/Summary" in h5out.root - assert "PAOLINA/Tracks" in h5out.root - assert "RECO/Events" in h5out.root + assert "Tracking/Tracks" in h5out.root + assert "Summary/Events" in h5out.root + assert "CHITS" in h5out.root + assert "CHITS/highTh" in h5out.root + assert "CHITS/lowTh" in h5out.root assert "Run" in h5out.root assert "Run/events" in h5out.root assert "Run/runInfo" in h5out.root - assert "Filters/NN_select" in h5out.root - assert "Filters/paolina_select" in h5out.root + assert "Filters/low_th_select" in h5out.root + assert "Filters/high_th_select" in h5out.root + assert "Filters/topology_select" in h5out.root + assert "DST/Events" in h5out.root + -@mark.serial def test_esmeralda_filters_events(KrMC_hdst_filename_toy, correction_map_MC_filename, config_tmpdir): PATH_IN = KrMC_hdst_filename_toy PATH_OUT = os.path.join(config_tmpdir, "Kr_tracks_with_MC_filtered.h5") conf = configure('dummy invisible_cities/config/esmeralda.conf'.split()) nevt_req = 8 - conf.update(dict(files_in = PATH_IN , - file_out = PATH_OUT , - event_range = nevt_req , - cor_hits_params = dict( - map_fname = correction_map_MC_filename, - threshold_charge_NN = 150 * units.pes , - threshold_charge_paolina = 200 * units.pes , - same_peak = True , - norm_strat = 'kr' , - apply_temp = False ))) + conf.update(dict(files_in = PATH_IN , + file_out = PATH_OUT , + event_range = nevt_req , + cor_hits_params = dict( + map_fname = correction_map_MC_filename, + threshold_charge_low = 150 * units.pes , + threshold_charge_high = 200 * units.pes , + same_peak = True , + apply_temp = False ))) cnt = esmeralda(**conf) - events_pass_NN = [0, 1, 2, 3, 4, 5, 6] + events_pass_low_th = [0, 1, 2, 3, 4, 5, 6] events_pass_paolina = [3, 4, 5, 6] nevt_in = cnt.events_in nevt_out = cnt.events_out assert nevt_req == nevt_in assert nevt_out == len(set(events_pass_paolina)) - df_hits_NN = dio.load_dst(PATH_OUT, 'RECO' , 'Events' ) - df_hits_paolina = dio.load_dst(PATH_OUT, 'PAOLINA', 'Events' ) - df_tracks_paolina = dio.load_dst(PATH_OUT, 'PAOLINA', 'Tracks' ) - df_summary_paolina = dio.load_dst(PATH_OUT, 'PAOLINA', 'Summary') + df_hits_low_th = dio.load_dst(PATH_OUT, 'CHITS' , 'lowTh' ) + df_hits_paolina = dio.load_dst(PATH_OUT, 'CHITS' , 'highTh') + df_tracks_paolina = dio.load_dst(PATH_OUT, 'Tracking', 'Tracks') + df_summary_paolina = dio.load_dst(PATH_OUT, 'Summary' , 'Events') - assert set(df_hits_NN .event.unique()) == set(events_pass_NN ) + assert set(df_hits_low_th .event.unique()) == set(events_pass_low_th ) assert set(df_hits_paolina.event.unique()) == set(events_pass_paolina) assert set(df_hits_paolina.event.unique()) == set(df_tracks_paolina .event.unique()) assert set(df_hits_paolina.event.unique()) == set(df_summary_paolina.event.unique()) + #assert event number in EventInfo and MC/Extents iqual to nevt_req + with tb.open_file(PATH_OUT) as h5out: + event_info = get_event_info(h5out) + assert length_of(event_info) == nevt_req + MC_num_evs = h5out.root.MC.extents[:]['evt_number'] + assert len(MC_num_evs) == nevt_req + + -@mark.serial def test_esmeralda_with_out_of_map_hits(KrMC_hdst_filename_toy, correction_map_MC_filename, config_tmpdir): PATH_IN = KrMC_hdst_filename_toy PATH_OUT = os.path.join(config_tmpdir, "Kr_tracks_with_MC_out_of_map.h5") conf = configure('dummy invisible_cities/config/esmeralda.conf'.split()) nevt_req = 8 - conf.update(dict(files_in = PATH_IN , - file_out = PATH_OUT , - event_range = nevt_req , - cor_hits_params = dict( - map_fname = correction_map_MC_filename, - threshold_charge_NN = 20 * units.pes , - threshold_charge_paolina = 20 * units.pes , - same_peak = True , - norm_strat = 'kr' , - apply_temp = False ))) + conf.update(dict(files_in = PATH_IN , + file_out = PATH_OUT , + event_range = nevt_req , + cor_hits_params = dict( + map_fname = correction_map_MC_filename, + threshold_charge_low = 20 * units.pes , + threshold_charge_high = 20 * units.pes , + same_peak = True , + apply_temp = False ))) cnt = esmeralda(**conf) - events_pass_NN = [0, 1, 2, 3, 4, 5, 6, 7] - events_pass_paolina = [0, 1, 2, 4, 5, 6, 7] + events_pass_low_th = [0, 1, 2, 3, 4, 5, 6, 7] + events_pass_paolina = [0, 1, 2, 3, 4, 5, 6, 7] nevt_in = cnt.events_in nevt_out = cnt.events_out assert nevt_req == nevt_in assert nevt_out == len(set(events_pass_paolina)) - df_hits_NN = dio.load_dst(PATH_OUT, 'RECO' , 'Events') - df_hits_paolina = dio.load_dst(PATH_OUT, 'PAOLINA', 'Events') + df_hits_low_th = dio.load_dst(PATH_OUT, 'CHITS', 'lowTh' ) + df_hits_paolina = dio.load_dst(PATH_OUT, 'CHITS', 'highTh') - assert set(df_hits_NN .event.unique()) == set(events_pass_NN ) + assert set(df_hits_low_th .event.unique()) == set(events_pass_low_th ) assert set(df_hits_paolina.event.unique()) == set(events_pass_paolina) + summary_table = dio.load_dst(PATH_OUT, 'Summary', 'Events') + #assert event with nan energy labeled in summary_table + events_energy = df_hits_paolina.groupby('event').Ec.apply(pd.Series.sum, skipna=False) + np.testing.assert_array_equal(summary_table.evt_out_of_map, np.isnan(events_energy.values)) + -@mark.serial def test_esmeralda_tracks_exact(data_hdst, esmeralda_tracks, correction_map_filename, config_tmpdir): PATH_IN = data_hdst PATH_OUT = os.path.join(config_tmpdir, "exact_tracks_esmeralda.h5") conf = configure('dummy invisible_cities/config/esmeralda.conf'.split()) nevt_req = all_events - conf.update(dict(files_in = PATH_IN , - file_out = PATH_OUT , - event_range = nevt_req , - run_number = 6822 , - cor_hits_params = dict( - map_fname = correction_map_filename, - threshold_charge_NN = 10 * units.pes , - threshold_charge_paolina = 30 * units.pes , - same_peak = True , - norm_strat = 'kr' , - apply_temp = False ), - paolina_params = dict( - vox_size = [15 * units.mm] * 3 , - energy_type = 'corrected' , - strict_vox_size = False , - energy_threshold = 0 * units.keV , - min_voxels = 2 , - blob_radius = 21 * units.mm) )) + conf.update(dict(files_in = PATH_IN , + file_out = PATH_OUT , + event_range = nevt_req , + run_number = 6822 , + cor_hits_params = dict( + map_fname = correction_map_filename, + threshold_charge_low = 10 * units.pes , + threshold_charge_high = 30 * units.pes , + same_peak = True , + apply_temp = False ), + paolina_params = dict( + vox_size = [15 * units.mm] * 3 , + strict_vox_size = False , + energy_threshold = 0 * units.keV , + min_voxels = 2 , + blob_radius = 21 * units.mm , + max_num_hits = 10000 ))) cnt = esmeralda(**conf) - df_tracks = dio.load_dst(PATH_OUT, 'PAOLINA', 'Tracks') + df_tracks = dio.load_dst(PATH_OUT, 'Tracking', 'Tracks' ) df_tracks_exact = pd.read_hdf(esmeralda_tracks, key = 'Tracks') - columns1 = df_tracks .columns columns2 = df_tracks_exact.columns - assert(sorted(columns1) == sorted(columns2)) - assert_dataframes_close (df_tracks[sorted(columns1)], df_tracks_exact[sorted(columns2)]) + #some events are not in df_tracks_exact + events = df_tracks_exact.event.unique() + df_tracks_cut = df_tracks[df_tracks.event.isin(events)] + assert_dataframes_close (df_tracks_cut[columns2], df_tracks_exact[columns2]) + #make sure out_of_map is true for events not in df_tracks_exact + diff_events = list(set(df_tracks.event.unique()).difference(events)) + df_summary = dio.load_dst(PATH_OUT, 'Summary', 'Events') + assert np.all(df_summary.loc[df_summary.event.isin(diff_events),'evt_out_of_map']) + + +def test_esmeralda_empty_input_file(config_tmpdir, ICDATADIR): + # Esmeralda must run on an empty file without raising any exception + # The input file has the complete structure of a PMAP but no events. + + PATH_IN = os.path.join(ICDATADIR , 'empty_hdst.h5') + PATH_OUT = os.path.join(config_tmpdir, 'empty_voxels.h5') + conf = configure('dummy invisible_cities/config/esmeralda.conf'.split()) + conf.update(dict(files_in = PATH_IN, + file_out = PATH_OUT)) + + esmeralda(**conf) + +#if the first analyzed events has no overlap in blob buggy esmeralda will cast all overlap energy to integers +def test_esmeralda_blob_overlap_bug(data_hdst, correction_map_filename, config_tmpdir): + PATH_IN = data_hdst + PATH_OUT = os.path.join(config_tmpdir, "exact_tracks_esmeralda_bug.h5") + conf = configure('dummy invisible_cities/config/esmeralda.conf'.split()) + nevt_req = 4, 8 + + conf.update(dict(files_in = PATH_IN , + file_out = PATH_OUT , + event_range = nevt_req , + run_number = 6822 , + cor_hits_params = dict( + map_fname = correction_map_filename, + threshold_charge_low = 10 * units.pes , + threshold_charge_high = 30 * units.pes , + same_peak = True , + apply_temp = False ), + paolina_params = dict( + vox_size = [15 * units.mm] * 3 , + strict_vox_size = False , + energy_threshold = 0 * units.keV , + min_voxels = 2 , + blob_radius = 21 * units.mm , + max_num_hits = 10000 ))) -def test_esmeralda_exact_result(ICDATADIR, KrMC_hdst_filename, correction_map_MC_filename, config_tmpdir): + cnt = esmeralda(**conf) + + df_tracks = dio.load_dst(PATH_OUT, 'Tracking', 'Tracks') + assert df_tracks['ovlp_blob_energy'].dtype == float + +def test_esmeralda_exact_result_all_events(ICDATADIR, KrMC_hdst_filename, correction_map_MC_filename, config_tmpdir): file_in = KrMC_hdst_filename file_out = os.path.join(config_tmpdir, "exact_Kr_tracks_with_MC.h5") conf = configure('dummy invisible_cities/config/esmeralda.conf'.split()) - true_out = os.path.join(ICDATADIR, "exact_Kr_tracks_with_MC.h5") + true_out = os.path.join(ICDATADIR, "exact_Kr_tracks_with_MC_KDST_no_filter.h5") nevt_req = all_events - conf.update(dict(files_in = file_in , - file_out = file_out , - event_range = nevt_req , - cor_hits_params = dict( - map_fname = correction_map_MC_filename, - threshold_charge_NN = 10 * units.pes , - threshold_charge_paolina = 20 * units.pes , - same_peak = True , - norm_strat = 'kr' , - apply_temp = False ), - paolina_params = dict( - vox_size = [15 * units.mm] * 3 , - energy_type = 'corrected' , - strict_vox_size = False , - energy_threshold = 0 * units.keV , - min_voxels = 2 , - blob_radius = 21 * units.mm ))) + conf.update(dict(files_in = file_in , + file_out = file_out , + event_range = nevt_req , + cor_hits_params = dict( + map_fname = correction_map_MC_filename, + threshold_charge_low = 10 * units.pes , + threshold_charge_high = 20 * units.pes , + same_peak = True , + apply_temp = False ), + paolina_params = dict( + vox_size = [15 * units.mm] * 3 , + strict_vox_size = False , + energy_threshold = 0 * units.keV , + min_voxels = 2 , + blob_radius = 21 * units.mm , + max_num_hits = 10000 ))) cnt = esmeralda(**conf) - tables = ( "MC/extents" , "MC/hits" , "MC/particles" , "MC/generators" , - "RECO/Events" , "PAOLINA/Events", "PAOLINA/Tracks" , "PAOLINA/Summary" , - "Run/events" , "Run/runInfo" , "Filters/NN_select", "Filters/paolina_select") + + + tables = ["MC/extents", "MC/hits", "MC/particles", + "Tracking/Tracks", "CHITS/lowTh", "CHITS/highTh", + "Run/events", "Run/runInfo", "DST/Events", "Summary/Events", + "Filters/low_th_select", "Filters/high_th_select", "Filters/topology_select"] with tb.open_file(true_out) as true_output_file: - with tb.open_file(file_out) as output_file: + with tb.open_file(file_out) as output_file: for table in tables: got = getattr( output_file.root, table) expected = getattr(true_output_file.root, table) assert_tables_equality(got, expected) -def test_esmeralda_empty_input_file(config_tmpdir, ICDATADIR): - # Esmeralda must run on an empty file without raising any exception - # The input file has the complete structure of a PMAP but no events. - PATH_IN = os.path.join(ICDATADIR , 'empty_hdst.h5') - PATH_OUT = os.path.join(config_tmpdir, 'empty_voxels.h5') +#test showing that all events that pass charge threshold are contained in hits output +def test_esmeralda_bug_duplicate_hits(data_hdst, correction_map_filename, config_tmpdir): + PATH_IN = data_hdst + PATH_OUT = os.path.join(config_tmpdir, "exact_tracks_esmeralda_drop_voxels_bug.h5") + conf = configure('dummy invisible_cities/config/esmeralda.conf'.split()) + nevt_req = 1 + conf.update(dict(files_in = PATH_IN , + file_out = PATH_OUT , + event_range = nevt_req , + run_number = 6822 , + cor_hits_params = dict( + map_fname = correction_map_filename, + threshold_charge_low = 10 * units.pes , + threshold_charge_high = 30 * units.pes , + same_peak = True , + apply_temp = False ), + paolina_params = dict( + vox_size = [15 * units.mm] * 3 , + strict_vox_size = False , + energy_threshold = 0 * units.keV , + min_voxels = 2 , + blob_radius = 21 * units.mm , + max_num_hits = 10000 ))) - conf = configure('dummy invisible_cities/config/esmeralda.conf'.split()) - conf.update(dict(files_in = PATH_IN, - file_out = PATH_OUT)) + cnt = esmeralda(**conf) - esmeralda(**conf) + df_tracks = dio.load_dst(PATH_OUT, 'Tracking', 'Tracks') + df_phits = dio.load_dst(PATH_OUT, 'CHITS' , 'highTh') + + for (event_num, ev_phits) in df_phits.groupby('event'): + assert sum(df_tracks[df_tracks.event==event_num].numb_of_hits) == len(ev_phits) + + +#test showing that all events that pass charge threshold are contained in hits output +def test_esmeralda_all_hits_after_drop_voxels(data_hdst, esmeralda_tracks, correction_map_filename, config_tmpdir): + PATH_IN = data_hdst + PATH_OUT = os.path.join(config_tmpdir, "exact_tracks_esmeralda_drop_voxels.h5") + conf = configure('dummy invisible_cities/config/esmeralda.conf'.split()) + nevt_req = all_events + th_p = 30 * units.pes + conf.update(dict(files_in = PATH_IN , + file_out = PATH_OUT , + event_range = nevt_req , + run_number = 6822 , + cor_hits_params = dict( + map_fname = correction_map_filename, + threshold_charge_low = 10 * units.pes , + threshold_charge_high = th_p , + same_peak = True , + apply_temp = False ), + paolina_params = dict( + vox_size = [15 * units.mm] * 3 , + strict_vox_size = False , + energy_threshold = 20 * units.keV , + min_voxels = 2 , + blob_radius = 21 * units.mm , + max_num_hits = 10000 ))) + cnt = esmeralda(**conf) + + df_phits = dio.load_dst(PATH_OUT, 'CHITS' , 'highTh') + df_in_hits = dio.load_dst(PATH_IN , 'RECO' , 'Events') + + num_pass_th_in_hits = sum(df_in_hits.Q >= th_p) + num_pass_th_p_hits = len(df_phits) + assert num_pass_th_in_hits == num_pass_th_p_hits + + #the number of finite Ep should be equal to Ec if no voxels were dropped. + num_paolina_hits = sum(np.isfinite(df_phits.Ep)) + assert num_paolina_hits <= num_pass_th_p_hits + + #check that the sum of Ep and Ec energies is the same + assert np.nansum(df_phits.Ec) == np.nansum(df_phits.Ep) + +def test_esmeralda_filters_events_with_too_many_hits(data_hdst, correction_map_filename, config_tmpdir): + PATH_IN = data_hdst + PATH_OUT = os.path.join(config_tmpdir, "esmeralda_filters_long_events.h5") + conf = configure('dummy invisible_cities/config/esmeralda.conf'.split()) + nevt_req = 9 + all_hits = 601 + nhits_max = 50 + paolina_events = {3021898, 3021914, 3021930, 3020951, 3020961} + conf.update(dict(files_in = PATH_IN , + file_out = PATH_OUT , + event_range = nevt_req , + run_number = 6822 , + cor_hits_params = dict( + map_fname = correction_map_filename, + threshold_charge_low = 10 * units.pes , + threshold_charge_high = 30 * units.pes , + same_peak = True , + apply_temp = False ), + paolina_params = dict( + vox_size = [15 * units.mm] * 3 , + strict_vox_size = False , + energy_threshold = 20 * units.keV , + min_voxels = 2 , + blob_radius = 21 * units.mm , + max_num_hits = nhits_max ))) + cnt = esmeralda(**conf) + + summary = dio.load_dst(PATH_OUT, 'Summary' , 'Events') + tracks = dio.load_dst(PATH_OUT, 'Tracking', 'Tracks') + hits = dio.load_dst(PATH_OUT, 'CHITS' , 'highTh') + + #assert only paolina_events inside tracks + assert set(tracks .event.unique()) == paolina_events + + #assert all events in summary table + assert summary.event.nunique() == nevt_req + #assert ntrk is 0 for non_paolina events + assert np.all(summary[ summary.event.isin(paolina_events)].evt_ntrks > 0 ) + assert np.all(summary[~summary.event.isin(paolina_events)].evt_ntrks == 0 ) + assert np.all(summary[~summary.event.isin(paolina_events)].evt_nhits > nhits_max) + + #assert all hits and events are in hits table + assert len(hits) == 601 + assert hits.event.nunique() == nevt_req + + #assert all events in topology filter with corresponding bool + topology_filter = dio.load_dst(PATH_OUT, 'Filters', 'topology_select') + assert len(topology_filter) == nevt_req + assert np.all(topology_filter[ topology_filter.event.isin(paolina_events)].passed == 1) + assert np.all(topology_filter[~topology_filter.event.isin(paolina_events)].passed == 0) diff --git a/invisible_cities/cities/irene.py b/invisible_cities/cities/irene.py index c69f978ec3..1d66fca786 100644 --- a/invisible_cities/cities/irene.py +++ b/invisible_cities/cities/irene.py @@ -22,14 +22,15 @@ from .. types.ic_types import minmax from .. database import load_db -from .. reco import tbl_functions as tbl -from .. reco import peak_functions as pkf -from .. core.random_sampling import NoiseSampler as SiPMsNoiseSampler -from .. io . pmaps_io import pmap_writer -from .. io. mcinfo_io import mc_info_writer -from .. io .run_and_event_io import run_and_event_writer -from .. io .trigger_io import trigger_writer -from .. io .event_filter_io import event_filter_writer +from .. reco import tbl_functions as tbl +from .. reco import peak_functions as pkf +from .. core.random_sampling import NoiseSampler as SiPMsNoiseSampler +from .. core.system_of_units_c import units +from .. io .pmaps_io import pmap_writer +from .. io .mcinfo_io import mc_info_writer +from .. io .run_and_event_io import run_and_event_writer +from .. io .trigger_io import trigger_writer +from .. io .event_filter_io import event_filter_writer from .. dataflow import dataflow as fl from .. dataflow.dataflow import push @@ -51,7 +52,8 @@ def irene(files_in, file_out, compression, event_range, print_mod, detector_db, run_number, n_baseline, n_mau, thr_mau, thr_sipm, thr_sipm_type, s1_lmin, s1_lmax, s1_tmin, s1_tmax, s1_rebin_stride, s1_stride, thr_csum_s1, - s2_lmin, s2_lmax, s2_tmin, s2_tmax, s2_rebin_stride, s2_stride, thr_csum_s2, thr_sipm_s2): + s2_lmin, s2_lmax, s2_tmin, s2_tmax, s2_rebin_stride, s2_stride, thr_csum_s2, thr_sipm_s2, + pmt_samp_wid=25*units.ns, sipm_samp_wid=1*units.mus): if thr_sipm_type.lower() == "common": # In this case, the threshold is a value in pes sipm_thr = thr_sipm @@ -87,7 +89,7 @@ def irene(files_in, file_out, compression, event_range, print_mod, detector_db, item = "sipm") # Build the PMap - compute_pmap = fl.map(build_pmap(detector_db, run_number, + compute_pmap = fl.map(build_pmap(detector_db, run_number, pmt_samp_wid, sipm_samp_wid, s1_lmax, s1_lmin, s1_rebin_stride, s1_stride, s1_tmax, s1_tmin, s2_lmax, s2_lmin, s2_rebin_stride, s2_stride, s2_tmax, s2_tmin, thr_sipm_s2), args = ("ccwfs", "s1_indices", "s2_indices", "sipm"), @@ -152,7 +154,7 @@ def irene(files_in, file_out, compression, event_range, print_mod, detector_db, -def build_pmap(detector_db, run_number, +def build_pmap(detector_db, run_number, pmt_samp_wid, sipm_samp_wid, s1_lmax, s1_lmin, s1_rebin_stride, s1_stride, s1_tmax, s1_tmin, s2_lmax, s2_lmin, s2_rebin_stride, s2_stride, s2_tmax, s2_tmin, thr_sipm_s2): s1_params = dict(time = minmax(min = s1_tmin, @@ -174,7 +176,8 @@ def build_pmap(detector_db, run_number, def build_pmap(ccwf, s1_indx, s2_indx, sipmzs): # -> PMap return pkf.get_pmap(ccwf, s1_indx, s2_indx, sipmzs, - s1_params, s2_params, thr_sipm_s2, pmt_ids) + s1_params, s2_params, thr_sipm_s2, pmt_ids, + pmt_samp_wid, sipm_samp_wid) return build_pmap diff --git a/invisible_cities/cities/irene_test.py b/invisible_cities/cities/irene_test.py index e3da969562..9ba965427a 100644 --- a/invisible_cities/cities/irene_test.py +++ b/invisible_cities/cities/irene_test.py @@ -4,6 +4,7 @@ import tables as tb import numpy as np +from pytest import raises from pytest import mark from pytest import fixture @@ -16,8 +17,9 @@ from .. io.run_and_event_io import read_run_and_event from .. evm.ic_containers import S12Params as S12P -from .. database import load_db -from .. database.load_db import DetDB +from .. database import load_db +from .. database.load_db import DetDB +from .. io .pmaps_io import load_pmaps from . irene import irene @@ -82,8 +84,8 @@ def test_irene_electrons_40keV(config_tmpdir, ICDATADIR, s12params, # since they are in general test-specific # NB: avoid taking defaults for run number (test-specific) - PATH_IN = os.path.join(ICDATADIR , 'electrons_40keV_z250_RWF.h5') - PATH_OUT = os.path.join(config_tmpdir, 'electrons_40keV_z250_CWF.h5') + PATH_IN = os.path.join(ICDATADIR , 'electrons_40keV_ACTIVE_10evts_RWF.h5') + PATH_OUT = os.path.join(config_tmpdir, 'electrons_40keV_ACTIVE_10evts_CWF.h5') nrequired = 2 @@ -210,23 +212,6 @@ def test_empty_events_issue_81(config_tmpdir, ICDATADIR, s12params): assert cnt. over_thr.n_failed == 1 -@mark.skip -def test_irene_electrons_40keV_pmt_active_is_correctly_set(job_info_missing_pmts, s12params): - "Check that PMT active correctly describes the PMT configuration of the detector" - nrequired = 1 - conf = configure('dummy invisible_cities/config/irene.conf'.split()) - conf.update(dict(run_number = job_info_missing_pmts.run_number, - files_in = job_info_missing_pmts. input_filename, - file_out = job_info_missing_pmts.output_filename, - event_range = (0, nrequired), - **unpack_s12params(s12params))) # s12params are just dummy values in this test - - irene = Irene(**conf) - - assert irene.pmt_active == job_info_missing_pmts.pmt_active - - -@mark.skip def test_irene_empty_pmap_output(ICDATADIR, output_tmpdir, s12params): file_in = os.path.join(ICDATADIR , "kr_rwf_0_0_7bar_NEXT_v1_00_05_v0.9.2_20171011_krmc_diomira_3evt.h5") file_out = os.path.join(output_tmpdir, "kr_rwf_0_0_7bar_NEXT_v1_00_05_v0.9.2_20171011_pmaps_3evt.h5") @@ -333,64 +318,6 @@ def test_irene_trigger_channels(config_tmpdir, ICDATADIR, s12params): np.testing.assert_array_equal(trigger_in, trigger_out) -@mark.skip(reason="Trigger split not implemented in liquid cities") -def test_irene_split_trigger(config_tmpdir, ICDATADIR, s12params): - PATH_IN = os.path.join(ICDATADIR , '6229_000_split_trigger.h5') - PATH_OUT1 = os.path.join(config_tmpdir, 'pmaps_6229_000_trg1.h5') - PATH_OUT2 = os.path.join(config_tmpdir, 'pmaps_6229_000_trg2.h5') - - nrequired = 3 - run_number = 6229 - - conf = configure('dummy invisible_cities/config/irene.conf'.split()) - conf.update(dict(run_number = run_number, - files_in = PATH_IN, - file_out = PATH_OUT1, - file_out2 = PATH_OUT2, - split_triggers = True, - trg1_code = 1, - trg2_code = 9, - event_range = (0, nrequired), - **unpack_s12params(s12params))) - - with warns(UserWarning) as record: - cnt = irene(**conf) - assert cnt.events_in == nrequired - - #check there is a warning for unknown trigger types - assert len(record) >= 1 - - evt_warn = 2 - trg_warn = 5 - message = "Event {} has an unknown trigger type ({})".format(evt_warn, trg_warn) - assert sum(1 for r in record if r.message.args[0] == message) == 1 - - # Check events has been properly redirected to their files - with tb.open_file(PATH_IN , mode='r') as h5in , \ - tb.open_file(PATH_OUT1, mode='r') as h5out1, \ - tb.open_file(PATH_OUT2, mode='r') as h5out2: - # There is only one event per file - assert h5out1.root.Trigger.events.shape[0] == 1 - assert h5out2.root.Trigger.events.shape[0] == 1 - - # Event number and trigger type are correct - evt1 = 0 - evt2 = 2 - trigger_in1 = h5in .root.Trigger.events[evt1] - trigger_in2 = h5in .root.Trigger.events[evt2] - trigger_out1 = h5out1.root.Trigger.events[0] - trigger_out2 = h5out2.root.Trigger.events[0] - np.testing.assert_array_equal(trigger_in1, trigger_out1) - np.testing.assert_array_equal(trigger_in2, trigger_out2) - - evt_in1 = h5in .root.Run.events[evt1][0] - evt_in2 = h5in .root.Run.events[evt2][0] - evt_out1 = h5out1.root.Run.events[0][0] - evt_out2 = h5out2.root.Run.events[0][0] - np.testing.assert_array_equal(evt_in1, evt_out1) - np.testing.assert_array_equal(evt_in2, evt_out2) - - def test_irene_exact_result(ICDATADIR, output_tmpdir): file_in = os.path.join(ICDATADIR , "Kr83_nexus_v5_03_00_ACTIVE_7bar_3evts.RWF.h5") file_out = os.path.join(output_tmpdir, "exact_result_irene.h5") @@ -461,3 +388,88 @@ def test_irene_empty_input_file(config_tmpdir, ICDATADIR): file_out = PATH_OUT)) irene(**conf) + + + +def test_irene_sequential_times(config_tmpdir, ICDATADIR): + + PATH_IN = os.path.join(ICDATADIR , 'single_evt_nonseqtime_rwf.h5') + PATH_OUT = os.path.join(config_tmpdir, 'test_pmaps.h5') + + conf = configure('dummy invisible_cities/config/irene.conf'.split()) + conf.update(dict(files_in = PATH_IN , + file_out = PATH_OUT , + run_number = 6351 , + n_baseline = 48000 , + thr_sipm = 1 * units.pes, + s1_tmin = 0 * units.mus, + s1_tmax = 640 * units.mus, + s1_lmin = 5 , + s1_lmax = 30 , + s2_tmin = 645 * units.mus, + s2_tmax = 1300 * units.mus, + s2_lmin = 80 , + s2_lmax = 200000 , + thr_sipm_s2 = 5 * units.pes)) + + irene(**conf) + + pmaps_out = load_pmaps(PATH_OUT) + + assert np.all(np.diff(pmaps_out[3348].s2s[1].times) > 0) + +def test_error_when_different_sample_widths(ICDATADIR, config_tmpdir): + ## Tests that ValueError message is correct + PATH_IN = os.path.join(ICDATADIR, 'run_2983.h5') + PATH_OUT = os.path.join(config_tmpdir, 'r2983_pmaps_diffparams.h5') + conf = configure('dummy invisible_cities/config/irene.conf'.split()) + + pmt_samp_wid = 25 * units.ns + sipm_samp_wid = 2 * units.mus + conf.update(dict(run_number = 2983, + files_in = PATH_IN, + file_out = PATH_OUT, + pmt_samp_wid = pmt_samp_wid, + sipm_samp_wid = sipm_samp_wid)) + + msg = "Shapes don't match!\n" + msg += "times has length 6\n" + msg += "pmts has length 6 \n" + msg += "sipms has length 3\n" + + with raises(ValueError) as error: + irene(**conf) + assert str(error.value) == msg + + +def test_irene_other_sample_widths(ICDATADIR, config_tmpdir): + ## Tests irene works when running with other sample frequencies + file_in = os.path.join(ICDATADIR, 'run_2983.h5') + file_out = os.path.join(config_tmpdir, 'r2983_pmaps_output.h5') + true_output = os.path.join(ICDATADIR, 'run_2983_pmaps_2evts_sfreq_5ns_200ns.h5') + conf = configure('dummy invisible_cities/config/irene.conf'.split()) + + pmt_samp_wid = 5 * units.ns + sipm_samp_wid = 0.2 * units.mus + n_events = 2 + + conf.update(dict(run_number = 2983, + files_in = file_in, + file_out = file_out, + event_range = (0, n_events), + pmt_samp_wid = pmt_samp_wid, + sipm_samp_wid = sipm_samp_wid)) + irene(**conf) + + tables = ( "PMAPS/S1" , "PMAPS/S2" , "PMAPS/S2Si" , + "PMAPS/S1Pmt" , "PMAPS/S2Pmt" , + "Run/events" , "Run/runInfo" , + "Trigger/events" , "Trigger/trigger" , + "Filters/s12_indices", "Filters/empty_pmap") + with tb.open_file(true_output) as true_output_file: + with tb.open_file(file_out) as output_file: + for table in tables: + got = getattr( output_file.root, table) + expected = getattr(true_output_file.root, table) + assert_tables_equality(got, expected) + diff --git a/invisible_cities/cities/paola.py b/invisible_cities/cities/paola.py new file mode 100644 index 0000000000..77bf85b478 --- /dev/null +++ b/invisible_cities/cities/paola.py @@ -0,0 +1,53 @@ +import os +import numpy as np +import tables as tb + +from .. dataflow import dataflow as fl + +from . components import city +from . components import hits_and_kdst_from_files + +from .. io.hits_io import hits_writer +from .. io.dst_io import _store_pandas_as_tables +from .. io.run_and_event_io import run_and_event_writer + + +def kdst_from_df_writer(h5out, compression='ZLIB4', group_name='DST', table_name='Events', descriptive_string='KDST Events', str_col_length=32): + """ + For a given open table returns a writer for KDST dataframe info + """ + def write_kdst(df): + return _store_pandas_as_tables(h5out=h5out, df=df, compression=compression, group_name=group_name, table_name=table_name, descriptive_string=descriptive_string, str_col_length=str_col_length) + return write_kdst + +@city +def paola(files_in, file_out, event_range, sel_event_file): + + selected_evts = np.loadtxt(sel_event_file, dtype=int) + + def event_sel(evt): + return evt in selected_evts + + evt_selector = fl.filter(event_sel, args='event_number') + + event_count_in = fl.spy_count() + event_count_out = fl.spy_count() + + with tb.open_file(file_out, "w") as h5out: + write_kdst = fl.sink(kdst_from_df_writer(h5out), args='kdst') + write_hits = fl.sink(hits_writer(h5out), args='hits') + write_event_info = fl.sink(run_and_event_writer(h5out), args=("run_number", "event_number", "timestamp")) + + result = fl.push(source = hits_and_kdst_from_files(files_in), + pipe = fl.pipe(#fl.slice(*event_range), + event_count_in.spy, + evt_selector, + event_count_out.spy, + fl.fork(write_kdst, + write_hits, + write_event_info) + ), + result = dict (events_in = event_count_in.future, + events_out = event_count_out.future )) + + return result diff --git a/invisible_cities/cities/paola_test.py b/invisible_cities/cities/paola_test.py new file mode 100644 index 0000000000..ae16211f35 --- /dev/null +++ b/invisible_cities/cities/paola_test.py @@ -0,0 +1,78 @@ +import pytest + +import os +import numpy as np +import tables as tb +import pandas as pd + +from .. core.configure import configure +from .. core.testing_utils import assert_dataframes_equal +from . paola import paola + +@pytest.fixture(scope = 'module') +def paola_file_input_data(ICDATADIR): + return os.path.join(ICDATADIR, 'hdst_paola_test.h5') + + +@pytest.fixture(scope = 'module') +def paola_sel_event_file(ICDATADIR): + return os.path.join(ICDATADIR, 'events_paola.txt') + + +@pytest.fixture(scope = 'module') +def evts_to_sel(ICDATADIR, paola_sel_event_file): + return np.loadtxt(paola_sel_event_file, dtype=int) + + +def test_selection_is_correct(paola_file_input_data, paola_sel_event_file, evts_to_sel, output_tmpdir): + + file_out = os.path.join(output_tmpdir,'hdst_paola_out.h5') + conf = configure('dummy invisible_cities/config/paola.conf'.split()) + conf.update(dict(files_in = paola_file_input_data, + file_out = file_out, + sel_event_file = paola_sel_event_file)) + res = paola(**conf) + + with tb.open_file(file_out) as h5out: + selected_evts = h5out.root.Run.events[:]['evt_number'] + + assert np.all(evts_to_sel == selected_evts) + + +def test_content_is_the_same(paola_file_input_data, paola_sel_event_file, evts_to_sel, output_tmpdir): + + file_out = os.path.join(output_tmpdir,'hdst_paola_out.h5') + conf = configure('dummy invisible_cities/config/paola.conf'.split()) + conf.update(dict(files_in = paola_file_input_data, + file_out = file_out, + sel_event_file = paola_sel_event_file)) + res = paola(**conf) + + with tb.open_file(paola_file_input_data) as h5in: + + in_table = getattr(getattr(h5in.root, 'RECO'), 'Events').read() + in_hits = pd.DataFrame.from_records(in_table) + sel_in_hits = in_hits[in_hits.event.isin(evts_to_sel)].reset_index(drop=True) + + in_table = getattr(getattr(h5in.root, 'DST'), 'Events').read() + in_kdsts = pd.DataFrame.from_records(in_table) + sel_in_kdsts = in_kdsts[in_kdsts.event.isin(evts_to_sel)].reset_index(drop=True) + + in_table = getattr(getattr(h5in.root, 'Run'), 'events').read() + in_times = pd.DataFrame.from_records(in_table) + sel_in_times = in_times[in_times.evt_number.isin(evts_to_sel)].reset_index(drop=True) + + with tb.open_file(file_out) as h5out: + + out_table = getattr(getattr(h5out.root, 'RECO'), 'Events').read() + out_hits = pd.DataFrame.from_records(out_table).reset_index(drop=True) + + out_table = getattr(getattr(h5out.root, 'DST'), 'Events').read() + out_kdsts = pd.DataFrame.from_records(out_table).reset_index(drop=True) + + out_table = getattr(getattr(h5out.root, 'Run'), 'events').read() + out_times = pd.DataFrame.from_records(out_table).reset_index(drop=True) + + assert_dataframes_equal(sel_in_hits, out_hits) + assert_dataframes_equal(sel_in_kdsts, out_kdsts) + assert_dataframes_equal(sel_in_times, out_times) diff --git a/invisible_cities/cities/penthesilea.py b/invisible_cities/cities/penthesilea.py index 35171b7c8f..a2504b8787 100644 --- a/invisible_cities/cities/penthesilea.py +++ b/invisible_cities/cities/penthesilea.py @@ -30,24 +30,25 @@ from .. database import load_db from .. core.system_of_units_c import units from .. reco import tbl_functions as tbl -from .. reco. pmaps_functions import RebinMethod -from .. io . hits_io import hits_writer -from .. io . mcinfo_io import mc_info_writer -from .. io .run_and_event_io import run_and_event_writer -from .. io. kdst_io import kr_writer -from .. io. event_filter_io import event_filter_writer - -from .. dataflow import dataflow as df -from .. dataflow.dataflow import push -from .. dataflow.dataflow import pipe - -from . components import city -from . components import print_every -from . components import peak_classifier -from . components import compute_xy_position -from . components import pmap_from_files -from . components import hit_builder -from . components import build_pointlike_event as build_pointlike_event_ +from .. reco. pmaps_functions import RebinMethod +from .. io . hits_io import hits_writer +from .. io . mcinfo_io import mc_info_writer +from .. io . run_and_event_io import run_and_event_writer +from .. io . kdst_io import kr_writer +from .. io . event_filter_io import event_filter_writer +from .. evm . pmaps import SiPMCharge + +from .. dataflow import dataflow as df +from .. dataflow.dataflow import push +from .. dataflow.dataflow import pipe + +from . components import city +from . components import print_every +from . components import peak_classifier +from . components import compute_xy_position +from . components import pmap_from_files +from . components import hit_builder +from . components import build_pointlike_event as build_pointlike_event_ @city def penthesilea(files_in, file_out, compression, event_range, print_mod, detector_db, run_number, @@ -56,7 +57,8 @@ def penthesilea(files_in, file_out, compression, event_range, print_mod, detecto s2_nmin, s2_nmax, s2_emin, s2_emax, s2_wmin, s2_wmax, s2_hmin, s2_hmax, s2_ethr, s2_nsipmmin, s2_nsipmmax, slice_reco_params = dict(), global_reco_params = dict(), - rebin_method = 'stride'): + rebin_method = 'stride', + sipm_charge_type = 'raw'): # slice_reco_params are qth, qlm, lm_radius, new_lm_radius, msipm used for hits reconstruction # global_reco_params are qth, qlm, lm_radius, new_lm_radius, msipm used for overall global (pointlike event) reconstruction @@ -69,7 +71,7 @@ def penthesilea(files_in, file_out, compression, event_range, print_mod, detecto pmap_select = df.count_filter(bool, args="pmap_passed") reco_algo_slice = compute_xy_position(detector_db, run_number, **slice_reco_params) - build_hits = df.map(hit_builder(detector_db, run_number, drift_v, reco_algo_slice, rebin, RebinMethod[rebin_method]), + build_hits = df.map(hit_builder(detector_db, run_number, drift_v, reco_algo_slice, rebin, RebinMethod[rebin_method], SiPMCharge[sipm_charge_type]), args = ("pmap", "selector_output", "event_number", "timestamp"), out = "hits" ) reco_algo_global = compute_xy_position(detector_db, run_number, **global_reco_params) diff --git a/invisible_cities/cities/penthesilea_test.py b/invisible_cities/cities/penthesilea_test.py index c896966b57..5a376c53d1 100644 --- a/invisible_cities/cities/penthesilea_test.py +++ b/invisible_cities/cities/penthesilea_test.py @@ -125,12 +125,42 @@ def test_penthesilea_threshold_rebin(ICDATADIR, output_tmpdir): conf = configure('dummy invisible_cities/config/penthesilea.conf'.split()) rebin_thresh = 4000 - conf.update(dict(run_number = -6340, - files_in = file_in, - file_out = file_out, - event_range = all_events, - rebin = rebin_thresh, - rebin_method = 'threshold')) + conf.update(dict(run_number = -6340, + files_in = file_in, + file_out = file_out, + event_range = all_events, + rebin = rebin_thresh, + rebin_method = 'threshold')) + + cnt = penthesilea(**conf) + + output_dst = dio.load_dst(file_out , 'RECO', 'Events') + expected_dst = dio.load_dst(true_output, 'RECO', 'Events') + + assert len(set(output_dst.event)) == len(set(expected_dst.event)) + assert_dataframes_close(output_dst[columns], expected_dst[columns], check_types=False) + + +def test_penthesilea_signal_to_noise(ICDATADIR, output_tmpdir): + file_in = os.path.join(ICDATADIR , "Kr83_nexus_v5_03_00_ACTIVE_7bar_3evts.PMP.h5") + file_out = os.path.join(output_tmpdir, "exact_result_penthesilea_SN.h5") + true_output = os.path.join(ICDATADIR , "Kr83_nexus_v5_03_00_ACTIVE_7bar_3evts_SN.HDST.h5") + + conf = configure('dummy invisible_cities/config/penthesilea.conf'.split()) + + reco_params = dict(Qthr = 2 , + Qlm = 6 , + lm_radius = 0 * units.mm, + new_lm_radius = 15 * units.mm, + msipm = 1 ) + + conf.update(dict(run_number = -6340, + files_in = file_in, + file_out = file_out, + event_range = all_events, + rebin = 2, + sipm_charge_type = 'signal_to_noise', + slice_reco_params = reco_params)) cnt = penthesilea(**conf) @@ -176,23 +206,6 @@ def test_penthesilea_true_hits_are_correct(KrMC_true_hits, config_tmpdir): (assert_MChit_equality(p_hit, t_hit) for p_hit, t_hit in zip(penthesilea_hits, true_hits)) -@mark.skip("This scenario is not possible in liquid cities") -def test_penthesilea_event_not_found(ICDATADIR, output_tmpdir): - file_in = os.path.join(ICDATADIR , "kr_rwf_0_0_7bar_NEXT_v1_00_05_v0.9.2_20171011_krmc_irene_3evt.h5") - file_out = os.path.join(output_tmpdir, "test_penthesilea_event_not_found.h5") - - conf = configure('dummy invisible_cities/config/penthesilea.conf'.split()) - nevt = 3 - - conf.update(dict(run_number = 4714, - files_in = file_in, - file_out = file_out, - event_range = (0, nevt))) - - cnt = penthesilea(**conf) - assert cnt.n_empty_pmaps == 1 - - def test_penthesilea_read_multiple_files(ICDATADIR, output_tmpdir): file_in = os.path.join(ICDATADIR , "Tl_v1_00_05_nexus_v5_02_08_7bar_pmaps_5evts_*.h5") file_out = os.path.join(output_tmpdir, "Tl_v1_00_05_nexus_v5_02_08_7bar_hits.h5") diff --git a/invisible_cities/config/esmeralda.conf b/invisible_cities/config/esmeralda.conf index 2840bd6fac..973174dd7d 100644 --- a/invisible_cities/config/esmeralda.conf +++ b/invisible_cities/config/esmeralda.conf @@ -8,19 +8,18 @@ run_number = 0 # How frequently to print events print_mod = 1 -cor_hits_params = dict( - map_fname = '$ICDIR/database/test_data/kr_emap_xy_100_100_r_6573_time.h5', - threshold_charge_NN = 6 * pes, - threshold_charge_paolina = 30 * pes, - same_peak = True, - norm_strat = 'kr', - apply_temp = True) +cor_hits_params = dict( + map_fname = '$ICDIR/database/test_data/kr_emap_xy_100_100_r_6573_time.h5', + threshold_charge_low = 6 * pes, + threshold_charge_high = 30 * pes, + same_peak = True, + apply_temp = True) paolina_params = dict( vox_size = [10 * mm, 10 * mm, 10 * mm], - energy_type = 'corrected', strict_vox_size = True, energy_threshold = 10 * keV, min_voxels = 2, - blob_radius = 21 * mm) + blob_radius = 21 * mm, + max_num_hits = 10000) diff --git a/invisible_cities/config/irene.conf b/invisible_cities/config/irene.conf index 45f6e793de..5c34cce170 100644 --- a/invisible_cities/config/irene.conf +++ b/invisible_cities/config/irene.conf @@ -30,7 +30,6 @@ thr_csum_s2 = 1.0 * pes thr_sipm = 3.5 * pes thr_sipm_type = "common" - # Set parameters to search for S1 # Notice that in MC file S1 is in t=100 mus s1_tmin = 99 * mus # position of S1 in MC files at 100 mus diff --git a/invisible_cities/config/paola.conf b/invisible_cities/config/paola.conf new file mode 100644 index 0000000000..5699adf5e5 --- /dev/null +++ b/invisible_cities/config/paola.conf @@ -0,0 +1,11 @@ +files_in = '/Users/paola/Software/IC-crash-course/data/hdst.h5' +file_out = '/Users/paola/Software/IC-crash-course/data/paola.h5' +#compression = 'ZLIB4' +event_range = 10 + +sel_event_file = '/Users/paola/Software/IC-crash-course/data/events.txt' +# run number 0 is for MC +#run_number = 0 + +# How frequently to print events +#print_mod = 1 diff --git a/invisible_cities/core/fit_functions_test.py b/invisible_cities/core/fit_functions_test.py index 86289650cb..6a8392cad8 100644 --- a/invisible_cities/core/fit_functions_test.py +++ b/invisible_cities/core/fit_functions_test.py @@ -54,12 +54,12 @@ def test_get_chi2_and_pvalue_when_data_equals_fit(): assert pvalue == approx(1., rel=1e-3) -@mark.skip(reason="Delaying elimination of solid cities") +@flaky(max_runs=5, min_passes=4) @given(floats(min_value = -2500, max_value = +2500), floats(min_value = + 100, max_value = + 300)) -@settings(max_examples=500) +@settings(max_examples=100) def test_get_chi2_and_pvalue_gauss_errors(mean, sigma): Nevt = int(1e6) ydata = np.random.normal(mean, sigma, Nevt) diff --git a/invisible_cities/core/random_sampling.py b/invisible_cities/core/random_sampling.py index 8d7c6d133d..5d964a6f45 100644 --- a/invisible_cities/core/random_sampling.py +++ b/invisible_cities/core/random_sampling.py @@ -2,12 +2,14 @@ import numpy as np -from typing import Tuple +from scipy.signal import fftconvolve -from functools import partial -from functools import lru_cache +from typing import Tuple -from .. database import load_db as DB +from functools import partial +from functools import lru_cache + +from .. database import load_db as DB class DarkModel(Enum): @@ -216,6 +218,9 @@ def signal_to_noise(self, ids : np.array, charges : np.array, An array of S/N values for the sipms in the slice. """ + if sample_width < 1: + error_msg = 'Unphysical sample width = {}'.format(sample_width) + raise ValueError(error_msg) dark_levels = self.dark_expectation(sample_width, dark_model) return charges / np.sqrt(charges + dark_levels[ids]) @@ -273,7 +278,7 @@ def multi_sample_distributions(self, sample_width : int) -> np.array: if sample_width == 1: return pad_pdfs(self.xbins, self.probs)[1] - mapping = map(np.convolve , + mapping = map(fftconvolve , self.multi_sample_distributions( 1), self.multi_sample_distributions(sample_width - 1), np.full(self.probs.shape[0], "same") ) diff --git a/invisible_cities/core/random_sampling_test.py b/invisible_cities/core/random_sampling_test.py index 96c3a2f911..fe2306663e 100644 --- a/invisible_cities/core/random_sampling_test.py +++ b/invisible_cities/core/random_sampling_test.py @@ -1,9 +1,10 @@ import numpy as np -from flaky import flaky -from pytest import mark +from flaky import flaky +from pytest import mark from pytest import fixture -from pytest import approx +from pytest import approx +from pytest import raises from hypothesis import given from hypothesis.strategies import composite @@ -142,7 +143,7 @@ def test_pad_pdfs(bins, expected_length): spectra = np.full((2, len(bins)), 0.1) _, padded_spectra = pad_pdfs(bins, spectra) assert padded_spectra.shape[1] == expected_length - + ## @mark.parametrize("bins", ## (np.arange(-1, 1, 0.01), @@ -153,7 +154,7 @@ def test_pad_pdfs(bins, expected_length): ## bin_diffs = np.diff(padded_bins) ## assert np.allclose(bin_diffs, bin_diffs[0]) - + @mark.parametrize("domain frequencies percentile true_value".split(), ((np.arange( 10), np.linspace(0, 1, 10), 0.6, 6), @@ -320,3 +321,13 @@ def test_noise_sampler_signal_to_noise(noise_sampler, assert len(signal_to_noise) == len(ids) assert np.allclose(np.round(signal_to_noise, 2), expected) + + +def test_signal_to_noise_zero_sample_raises_error(noise_sampler): + + noise_sampler, *_ = noise_sampler + ids = [1000, 1010, 1023] + qs = np.array([0.2, 1, 1.03]) + with raises(ValueError): + + noise_sampler.signal_to_noise(ids, qs, 0) diff --git a/invisible_cities/database/localdb.DEMOPPDB.sqlite3 b/invisible_cities/database/localdb.DEMOPPDB.sqlite3 index 4219f0faf2..5ec4354ca5 100644 --- a/invisible_cities/database/localdb.DEMOPPDB.sqlite3 +++ b/invisible_cities/database/localdb.DEMOPPDB.sqlite3 @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:27bca7c34a8b1c8d06f5806a0e24e5f9920bcfa86f3cf777788cb4aed8babcbb -size 2482176 +oid sha256:543c744f37ad6c61fefea492dfb7f99117bd74b7cfab2f706cf34283de714ef1 +size 5324800 diff --git a/invisible_cities/database/localdb.NEWDB.sqlite3 b/invisible_cities/database/localdb.NEWDB.sqlite3 index eb915fae09..35dd5af9d9 100644 --- a/invisible_cities/database/localdb.NEWDB.sqlite3 +++ b/invisible_cities/database/localdb.NEWDB.sqlite3 @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:b0ff72e35811fff697fe46cf1aa32e070b9e7ce5e8d82b79481b3712f81ac800 -size 218296320 +oid sha256:9bf7b0eef6cfe49d787116383538a1fa0b9e2105e41ca0e8299e8a94450b8ddb +size 300511232 diff --git a/invisible_cities/database/test_data/Kr83_nexus_v5_03_00_ACTIVE_7bar_3evts.HDST.h5 b/invisible_cities/database/test_data/Kr83_nexus_v5_03_00_ACTIVE_7bar_3evts.HDST.h5 index 638a4416eb..e24a02a776 100644 --- a/invisible_cities/database/test_data/Kr83_nexus_v5_03_00_ACTIVE_7bar_3evts.HDST.h5 +++ b/invisible_cities/database/test_data/Kr83_nexus_v5_03_00_ACTIVE_7bar_3evts.HDST.h5 @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:2dbcbe23e749112af9c4519ecfd9487e39b06132b30bfeff9e4f85089efccf61 -size 80343 +oid sha256:66eea9f1090974baa45e4b07f0f074a77eb8afa0c32efcca189843760c58cd9a +size 80218 diff --git a/invisible_cities/database/test_data/Kr83_nexus_v5_03_00_ACTIVE_7bar_3evts.PMP.h5 b/invisible_cities/database/test_data/Kr83_nexus_v5_03_00_ACTIVE_7bar_3evts.PMP.h5 index c8c5e5f20f..e79d89fe50 100644 --- a/invisible_cities/database/test_data/Kr83_nexus_v5_03_00_ACTIVE_7bar_3evts.PMP.h5 +++ b/invisible_cities/database/test_data/Kr83_nexus_v5_03_00_ACTIVE_7bar_3evts.PMP.h5 @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:e9c5b40c16295bcbd25ac4b1be1021b10a3b63e2d9b37a6dbd7437ac88766aa1 +oid sha256:7284c833efc80b2879bffae58fdeb85804f0c2665ff54595ecbe659b383b8625 size 130492 diff --git a/invisible_cities/database/test_data/Kr83_nexus_v5_03_00_ACTIVE_7bar_3evts_SN.HDST.h5 b/invisible_cities/database/test_data/Kr83_nexus_v5_03_00_ACTIVE_7bar_3evts_SN.HDST.h5 new file mode 100644 index 0000000000..278956d7ce --- /dev/null +++ b/invisible_cities/database/test_data/Kr83_nexus_v5_03_00_ACTIVE_7bar_3evts_SN.HDST.h5 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:159a3940a58bf3b4e66c2f71f37ff37b36ccb74a3fbc0b7b4bf228dc46ab48f6 +size 79636 diff --git a/invisible_cities/database/test_data/Kr83_nexus_v5_03_00_ACTIVE_7bar_noS1.HDST.h5 b/invisible_cities/database/test_data/Kr83_nexus_v5_03_00_ACTIVE_7bar_noS1.HDST.h5 index 4c0c04857a..260e897d5c 100644 --- a/invisible_cities/database/test_data/Kr83_nexus_v5_03_00_ACTIVE_7bar_noS1.HDST.h5 +++ b/invisible_cities/database/test_data/Kr83_nexus_v5_03_00_ACTIVE_7bar_noS1.HDST.h5 @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:bacc03d7bca51fbe59932c04b9851e536296d75e7b44b16051537e5462e92a16 -size 89333 +oid sha256:0cb15ade893c41f23851e5e3feecc42c3f3f4f183980f8f869a124e3580ebb84 +size 89541 diff --git a/invisible_cities/database/test_data/electrons_40keV_ACTIVE_10evts_RWF.h5 b/invisible_cities/database/test_data/electrons_40keV_ACTIVE_10evts_RWF.h5 new file mode 100644 index 0000000000..c7d68c9527 --- /dev/null +++ b/invisible_cities/database/test_data/electrons_40keV_ACTIVE_10evts_RWF.h5 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:ca43fc96b2393385371b54cdf8af38456b516c85f0097a8f257f757b9921372c +size 1497328 diff --git a/invisible_cities/database/test_data/events_paola.txt b/invisible_cities/database/test_data/events_paola.txt new file mode 100644 index 0000000000..401d2ef92e --- /dev/null +++ b/invisible_cities/database/test_data/events_paola.txt @@ -0,0 +1,2 @@ +1046622 +1046748 diff --git a/invisible_cities/database/test_data/exact_Kr_tracks_with_MC.h5 b/invisible_cities/database/test_data/exact_Kr_tracks_with_MC.h5 index ec07c7677e..9eaeedf358 100644 --- a/invisible_cities/database/test_data/exact_Kr_tracks_with_MC.h5 +++ b/invisible_cities/database/test_data/exact_Kr_tracks_with_MC.h5 @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:bdb76e4c655a97d83088f097c6b71cf301a82b5b13d390869297a0867ea6be64 -size 110403 +oid sha256:b30a5f6c70335436a32e512d755cd4ff0112deec6b85e0864b7e06747c8e6a1a +size 122010 diff --git a/invisible_cities/database/test_data/exact_Kr_tracks_with_MC_KDST.h5 b/invisible_cities/database/test_data/exact_Kr_tracks_with_MC_KDST.h5 new file mode 100644 index 0000000000..b27da6ef7e --- /dev/null +++ b/invisible_cities/database/test_data/exact_Kr_tracks_with_MC_KDST.h5 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:fc3b1553d8a69f8afb8b4c6beebd4f1f0b8a9eb55b23e97ccdf0b925738afb6d +size 118040 diff --git a/invisible_cities/database/test_data/exact_Kr_tracks_with_MC_KDST_no_filter.h5 b/invisible_cities/database/test_data/exact_Kr_tracks_with_MC_KDST_no_filter.h5 new file mode 100644 index 0000000000..cacdd914c1 --- /dev/null +++ b/invisible_cities/database/test_data/exact_Kr_tracks_with_MC_KDST_no_filter.h5 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:25880bc28e98b1d8c8098a5d37f872470408f090deded3aba73e2bf25b229e18 +size 123577 diff --git a/invisible_cities/database/test_data/hdst_paola_test.h5 b/invisible_cities/database/test_data/hdst_paola_test.h5 new file mode 100644 index 0000000000..c10dea8085 --- /dev/null +++ b/invisible_cities/database/test_data/hdst_paola_test.h5 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:124e1663f7462ba834b0669800ca38f2046d660d080f03f0f8e6963ce26e083e +size 90645 diff --git a/invisible_cities/database/test_data/pmaps_negative_bins.h5 b/invisible_cities/database/test_data/pmaps_negative_bins.h5 new file mode 100644 index 0000000000..427b0d1fba --- /dev/null +++ b/invisible_cities/database/test_data/pmaps_negative_bins.h5 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:72b66655b67d3c542c931686745ae7dca70326b87bd96648aaf17a247c90733b +size 141417 diff --git a/invisible_cities/database/test_data/run_2983_pmaps_2evts_sfreq_5ns_200ns.h5 b/invisible_cities/database/test_data/run_2983_pmaps_2evts_sfreq_5ns_200ns.h5 new file mode 100644 index 0000000000..8aed94f464 --- /dev/null +++ b/invisible_cities/database/test_data/run_2983_pmaps_2evts_sfreq_5ns_200ns.h5 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:c4bb3cbcb78833133174e30a93a60215dfe8143451f98ebc1364ba946686b796 +size 87877 diff --git a/invisible_cities/database/test_data/single_evt_nonseqtime_rwf.h5 b/invisible_cities/database/test_data/single_evt_nonseqtime_rwf.h5 new file mode 100644 index 0000000000..fe93cfa961 --- /dev/null +++ b/invisible_cities/database/test_data/single_evt_nonseqtime_rwf.h5 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:8fcd19eab10e9a298eb60110b0f96b616934f6183293fc87c8cb74e5ab6a61f0 +size 591126 diff --git a/invisible_cities/evm/event_model.py b/invisible_cities/evm/event_model.py index dbf877f1f7..614a07e28e 100644 --- a/invisible_cities/evm/event_model.py +++ b/invisible_cities/evm/event_model.py @@ -116,6 +116,7 @@ def __str__(self): class HitEnergy(AutoNameEnumBase): E = auto() Ec = auto() + Ep = auto() class BHit: @@ -225,7 +226,7 @@ def __str__(self): class Hit(Cluster): """Represents a reconstructed hit (cluster + z + energy)""" def __init__(self, peak_number, cluster, z, s2_energy, peak_xy, - s2_energy_c=-1, track_id=-1): + s2_energy_c=-1, track_id=-1, Ep=-1): super().__init__(cluster.Q, @@ -237,13 +238,14 @@ def __init__(self, peak_number, cluster, z, s2_energy, peak_xy, self.Ypeak = peak_xy.y self.Ec = s2_energy_c self.track_id = track_id + self.Ep = Ep @property def npeak(self): return self.peak_number def __str__(self): - return """<{} : npeak = {} z = {} XYpeak = {}, {} E = {} cluster ={} >""".format(self.__class__.__name__, - self.npeak, self.Z, self.Xpeak, self.Ypeak, self.E, super().__str__()) + return """<{} : npeak = {} z = {} XYpeak = {}, {} E = {} Ec = {} Ep = {} trackid = {} cluster ={} >""".format(self.__class__.__name__, + self.npeak, self.Z, self.Xpeak, self.Ypeak, self.E, self.Ec, self.Ep, self.track_id, super().__str__()) __repr__ = __str__ @@ -268,7 +270,7 @@ def __repr__(self): return self.__str__() -class Blob(): +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], @@ -352,9 +354,9 @@ def __str__(self): class HitCollection(Event): """A Collection of hits""" - def __init__(self, event_number, event_time): + def __init__(self, event_number, event_time, hits=None): Event.__init__(self, event_number, event_time) - self.hits = [] + self.hits = [] if hits is None else hits def store(self, table): row = table.row @@ -375,6 +377,7 @@ def store(self, table): row["Qc" ] = hit .Qc row["Ec" ] = hit .Ec row["track_id"] = hit .track_id + row["Ep" ] = hit .Ep row.append() def __str__(self): diff --git a/invisible_cities/evm/event_model_test.py b/invisible_cities/evm/event_model_test.py index 3e7554d1f8..189f80fb2e 100644 --- a/invisible_cities/evm/event_model_test.py +++ b/invisible_cities/evm/event_model_test.py @@ -4,6 +4,9 @@ from pytest import mark from hypothesis import given +from hypothesis.strategies import just +from hypothesis.strategies import lists +from hypothesis.strategies import one_of from hypothesis.strategies import floats from hypothesis.strategies import integers from hypothesis.strategies import composite @@ -18,7 +21,7 @@ from . event_model import Hit from . event_model import Voxel from . event_model import HitCollection -from . event_model import HitCollection +from . event_model import HitEnergy from . event_model import KrEvent @@ -37,8 +40,9 @@ def event_input(draw): time = draw(floats (allow_nan=False)) return evt_no, time + @composite -def voxel_input(draw, min_value=0, max_value=100): +def voxel_input(draw): x = draw(floats( 1, 5)) y = draw(floats(-10, 10)) z = draw(floats(.01, .5)) @@ -47,7 +51,7 @@ def voxel_input(draw, min_value=0, max_value=100): return x, y, z, E, size @composite -def cluster_input(draw, min_value=0, max_value=100): +def cluster_input(draw): x = draw(floats ( 1, 5)) y = draw(floats (-10, 10)) xvar = draw(floats (.01, .5)) @@ -58,13 +62,26 @@ def cluster_input(draw, min_value=0, max_value=100): @composite -def hit_input(draw, min_value=0, max_value=100): +def hit_input(draw): z = draw(floats (.1, .9)) s2_energy = draw(floats (50, 100)) peak_number = draw(integers( 1, 20)) x_peak = draw(floats (-10., 2.)) y_peak = draw(floats (-20., 5.)) - return peak_number, s2_energy, z, x_peak, y_peak + s2_energy_c = draw(one_of(just(-1), floats (50, 100))) + track_id = draw(one_of(just(-1), integers( 0, 10))) + Ep = draw(one_of(just(-1), floats (50, 100))) + return peak_number, s2_energy, z, x_peak, y_peak, s2_energy_c, track_id, Ep + + +@composite +def hits(draw): + Q, x, y, xvar, yvar, nsipm = draw(cluster_input()) + peak_number, E, z, x_peak, y_peak, s2ec, track_id, ep = draw( hit_input()) + + c = Cluster(Q, xy(x,y), xy(xvar,yvar), nsipm) + h = Hit(peak_number, c, z, E, xy(x_peak, y_peak), s2ec, track_id, ep) + return h @given(sensor_params_input()) @@ -80,7 +97,6 @@ def test_sensor_params(sensor_pars): @mark.parametrize("test_class", (Event, - HitCollection, HitCollection, KrEvent)) @given(event_input()) @@ -92,8 +108,7 @@ def test_event(test_class, event_pars): assert evt.time == time - -@given(cluster_input(1)) +@given(cluster_input()) def test_cluster(ci): Q, x, y, xvar, yvar, nsipm = ci xrms = np.sqrt(xvar) @@ -105,61 +120,73 @@ def test_cluster(ci): c = Cluster(Q, xy(x,y), xy(xvar,yvar), nsipm) assert c.nsipm == nsipm - np.isclose (c.Q , Q, rtol=1e-4) - np.isclose (c.X , x, rtol=1e-4) - np.isclose (c.Y , y, rtol=1e-4) - np.isclose (c.Xrms, xrms, rtol=1e-3) - np.isclose (c.Yrms, yrms, rtol=1e-3) - np.isclose (c.var.XY, varar, rtol=1e-4) - np.allclose(c.XY, xyar, rtol=1e-4) - np.isclose (c.R , r, rtol=1e-4) - np.isclose (c.Phi , phi, rtol=1e-4) - np.allclose(c.posxy , pos, rtol=1e-4) - - -@given(cluster_input(1), hit_input(1)) + np.isclose (c.Q , Q , rtol=1e-4) + np.isclose (c.X , x , rtol=1e-4) + np.isclose (c.Y , y , rtol=1e-4) + np.isclose (c.Xrms , xrms , rtol=1e-4) + np.isclose (c.Yrms , yrms , rtol=1e-4) + np.isclose (c.var.XY, varar, rtol=1e-4) + np.allclose(c.XY , xyar , rtol=1e-4) + np.isclose (c.R , r , rtol=1e-4) + np.isclose (c.Phi , phi , rtol=1e-4) + np.allclose(c.posxy , pos , rtol=1e-4) + + +@mark.parametrize("value", "E Ec Ep".split()) +def test_hitenergy_value(value): + assert getattr(HitEnergy, value).value == value + + +@given(cluster_input(), hit_input()) def test_hit(ci, hi): - Q, x, y, xvar, yvar, nsipm = ci - peak_number, E, z, x_peak, y_peak = hi + Q, x, y, xvar, yvar, nsipm = ci + peak_number, E, z, x_peak, y_peak, s2ec, track_id, ep = hi xyz = x, y, z c = Cluster(Q, xy(x,y), xy(xvar,yvar), nsipm) - h = Hit(peak_number, c, z, E, xy(x_peak, y_peak)) + h = Hit(peak_number, c, z, E, xy(x_peak, y_peak), s2ec, track_id, ep) assert h.peak_number == peak_number assert h.npeak == peak_number - np.isclose (h.X , x , rtol=1e-4) - np.isclose (h.Y , y , rtol=1e-4) - np.isclose (h.Z , z, rtol=1e-4) - np.isclose (h.E , E, rtol=1e-4) - np.isclose (h.Xpeak , x_peak, rtol=1e-4) - np.isclose (h.Ypeak , y_peak, rtol=1e-4) - np.allclose(h.XYZ , xyz, rtol=1e-4) - np.allclose(h.pos , np.array(xyz), rtol=1e-4) + np.isclose (h.X , x , rtol=1e-4) + np.isclose (h.Y , y , rtol=1e-4) + np.isclose (h.Z , z , rtol=1e-4) + np.isclose (h.E , E , rtol=1e-4) + np.isclose (h.Xpeak , x_peak , rtol=1e-4) + np.isclose (h.Ypeak , y_peak , rtol=1e-4) + np.allclose(h.XYZ , xyz , rtol=1e-4) + np.allclose(h.pos , xyz , rtol=1e-4) + np.allclose(h.Ec , s2ec , rtol=1e-4) + np.allclose(h.track_id, track_id, rtol=1e-4) + np.allclose(h.Ep , ep , rtol=1e-4) -@given(voxel_input(1)) + +@given(voxel_input()) def test_voxel(vi): x, y, z, E, size = vi xyz = x, y, z v = Voxel(x, y, z, E, size) - np.allclose(v.XYZ , xyz, rtol=1e-4) - np.allclose(v.pos , np.array(xyz), rtol=1e-4) - np.isclose (v.E , E, rtol=1e-4) - np.isclose (v.X , x , rtol=1e-4) - np.isclose (v.Y , y , rtol=1e-4) - np.isclose (v.Z , z , rtol=1e-4) + np.allclose(v.XYZ, xyz, rtol=1e-4) + np.allclose(v.pos, xyz, rtol=1e-4) + np.isclose (v.E , E , rtol=1e-4) + np.isclose (v.X , x , rtol=1e-4) + np.isclose (v.Y , y , rtol=1e-4) + np.isclose (v.Z , z , rtol=1e-4) -@mark.parametrize("test_class", - (HitCollection, - HitCollection)) -def test_hit_collection_empty(test_class): - hc = test_class(-1, -1) +def test_hit_collection_empty(): + hc = HitCollection(-1, -1) assert hc.hits == [] +@given(lists(hits())) +def test_hit_collection_nonempty(hits): + hc = HitCollection(-1, -1, hits=hits) + assert hc.hits == hits + + def test_kr_event_attributes(): evt = KrEvent(-1, -1) diff --git a/invisible_cities/evm/nh5.py b/invisible_cities/evm/nh5.py index f0ac04ab09..9998962927 100644 --- a/invisible_cities/evm/nh5.py +++ b/invisible_cities/evm/nh5.py @@ -76,17 +76,17 @@ class MCHitInfo(tb.IsDescription): class MCParticleInfo(tb.IsDescription): """Stores the simulated particles as metadata using Pytables. """ - particle_indx = tb. Int16Col( pos=0) - particle_name = tb. StringCol(20, pos=1) - primary = tb. Int16Col( pos=2) - mother_indx = tb. Int16Col( pos=3) - initial_vertex = tb.Float32Col( pos=4, shape=4) - final_vertex = tb.Float32Col( pos=5, shape=4) - initial_volume = tb. StringCol(20, pos=6) - final_volume = tb. StringCol(20, pos=7) - momentum = tb.Float32Col( pos=8, shape=3) - kin_energy = tb.Float32Col( pos=9) - creator_proc = tb. StringCol(20, pos=10) + particle_indx = tb. Int16Col( pos= 0) + particle_name = tb. StringCol( 20, pos= 1) + primary = tb. Int16Col( pos= 2) + mother_indx = tb. Int16Col( pos= 3) + initial_vertex = tb.Float32Col( pos= 4, shape=4) + final_vertex = tb.Float32Col( pos= 5, shape=4) + initial_volume = tb. StringCol( 20, pos= 6) + final_volume = tb. StringCol( 20, pos= 7) + momentum = tb.Float32Col( pos= 8, shape=3) + kin_energy = tb.Float32Col( pos= 9) + creator_proc = tb. StringCol(100, pos=10) class SENSOR_WF(tb.IsDescription): @@ -233,6 +233,7 @@ class HitsTable(tb.IsDescription): Qc = tb.Float64Col(pos=13) Ec = tb.Float64Col(pos=14) track_id = tb. Int32Col(pos=15) + Ep = tb.Float64Col(pos=16) class VoxelsTable(tb.IsDescription): diff --git a/invisible_cities/evm/pmaps.py b/invisible_cities/evm/pmaps.py index 96475d6867..0a2abfcc30 100644 --- a/invisible_cities/evm/pmaps.py +++ b/invisible_cities/evm/pmaps.py @@ -2,8 +2,16 @@ import numpy as np +from enum import auto + from .. core.system_of_units_c import units from .. core.core_functions import weighted_mean_and_std +from .. types.ic_types import AutoNameEnumBase + + +class SiPMCharge(AutoNameEnumBase): + raw = auto() + signal_to_noise = auto() class PMap: @@ -98,7 +106,32 @@ def _check_valid_input(self, times, pmts, sipms): raise ValueError(msg) class S1(_Peak): pass -class S2(_Peak): pass +class S2(_Peak): + + def sipm_charge_array(self, noise_func, charge_type, + single_point=False): + + if charge_type is SiPMCharge.raw: + if single_point: + return self.sipms.sum_over_times + return self.sipms.all_waveforms.T + + ## Should generalise to use base sample width rather + ## than mus in case of electronics changes in the + ## future. + if single_point: + charges = self.sipms.sum_over_times + sample_wid = int(np.ceil(np.round(self.width) / units.mus)) + return noise_func(self.sipms.ids, charges, sample_wid) + + sample_widths = np.ceil(np.round(self.bin_widths) / units.mus) + id_sample_repeat = np.repeat(self.sipms.ids[np.newaxis], + len(self.bin_widths), 0) + mapping = map(noise_func , + id_sample_repeat , + self.sipms.all_waveforms.T, + sample_widths.astype(int) ) + return tuple(mapping) class _SensorResponses: diff --git a/invisible_cities/evm/pmaps_test.py b/invisible_cities/evm/pmaps_test.py index 625e215cd1..241ff73aa6 100644 --- a/invisible_cities/evm/pmaps_test.py +++ b/invisible_cities/evm/pmaps_test.py @@ -1,28 +1,36 @@ import numpy as np -from pytest import approx -from pytest import raises -from pytest import mark - -from hypothesis import assume -from hypothesis import given -from hypothesis.strategies import integers -from hypothesis.strategies import floats -from hypothesis.strategies import sampled_from -from hypothesis.strategies import composite -from hypothesis.extra.numpy import arrays - -from .. core.core_functions import weighted_mean_and_std -from .. core.testing_utils import exactly -from .. core.testing_utils import assert_SensorResponses_equality -from .. core.testing_utils import assert_Peak_equality -from .. core.testing_utils import previous_float +from numpy.testing import assert_allclose +from pytest import approx +from pytest import raises +from pytest import mark +from pytest import fixture + +from hypothesis import assume +from hypothesis import given +from hypothesis import settings +from hypothesis.strategies import integers +from hypothesis.strategies import floats +from hypothesis.strategies import sampled_from +from hypothesis.strategies import composite +from hypothesis.extra.numpy import arrays + +from .. core.core_functions import weighted_mean_and_std +from .. core.random_sampling import NoiseSampler +from .. core.system_of_units_c import units +from .. core.testing_utils import exactly +from .. core.testing_utils import assert_SensorResponses_equality +from .. core.testing_utils import assert_Peak_equality +from .. core.testing_utils import previous_float + +from invisible_cities.database import load_db as DB from . pmaps import PMTResponses from . pmaps import SiPMResponses -from . pmaps import S1 -from . pmaps import S2 -from . pmaps import PMap +from . pmaps import S1 +from . pmaps import S2 +from . pmaps import PMap +from . pmaps import SiPMCharge wf_min = 0 @@ -45,7 +53,7 @@ def sensor_responses(draw, n_samples=None, subtype=None, ids=None): @composite def peaks(draw, subtype=None, pmt_ids=None, with_sipms=True): - nsamples = draw(integers(1, 50)) + nsamples = draw(integers(1, 20)) _, pmt_r = draw(sensor_responses(nsamples, PMTResponses, pmt_ids)) sipm_r = SiPMResponses.build_empty_instance() assume(pmt_r.sum_over_sensors[ 0] != 0) @@ -56,13 +64,16 @@ def peaks(draw, subtype=None, pmt_ids=None, with_sipms=True): if with_sipms: _, sipm_r = draw(sensor_responses(nsamples, SiPMResponses)) - times = draw(arrays(float, nsamples, - floats(min_value=0, max_value=1e3), - unique = True)) - bin_widths = [1] + times = draw(arrays(float, nsamples, + floats(min_value=0, max_value=1e3), + unique = True).map(sorted)) + + bin_widths = np.array([1]) if len(times) > 1: - bin_widths = np.append(np.diff(times), max(np.diff(times))) - args = np.sort(times), bin_widths, pmt_r, sipm_r + time_differences = np.diff(times) + bin_widths = np.append(time_differences, max(time_differences)) + + args = times, bin_widths, pmt_r, sipm_r return args, subtype(*args) @@ -305,3 +316,161 @@ def test_PMap_s2s(pmps): assert len(pmp.s2s) == len(s2s) for kept_s2, true_s2 in zip(pmp.s2s, s2s): assert_Peak_equality(kept_s2, true_s2) + + + +@fixture(scope='module') +def signal_to_noise_6400(): + return NoiseSampler('new', 6400).signal_to_noise + + +@fixture(scope='module') +def s2_peak(): + times = np.arange(20) * units.mus + bin_widths = np.full_like(times, units.mus) + + pmt_ids = DB.DataPMT ('new', 6400).SensorID.values + sipm_ids = DB.DataSiPM('new', 6400).index .values + + np.random.seed(22) + pmts = PMTResponses(pmt_ids, + np.random.uniform(0, 100, + (len(pmt_ids), + len(times)))) + + sipms = SiPMResponses(sipm_ids, + np.random.uniform(0, 10, + (len(sipm_ids), + len(times)))) + + ## Values for comparison in tests + chan_slice = slice(3, 9) + raw_arr = [[7.51759755, 5.60509619, 6.42950506, + 9.20429423, 9.70312128, 2.02362046], + [7.05972334, 2.22846472, 5.46749176, + 7.31111882, 3.40960906, 4.56294121], + [6.70472842, 3.59272265, 5.16309049, + 1.56536978, 9.81538459, 6.82850322], + [4.15689334, 3.41990664, 5.79051761, + 0.81210519, 1.0304272 , 0.40297561], + [4.73786661, 0.6207543 , 1.6485871 , + 2.52693314, 8.91532864, 1.52840296], + [5.92994762, 3.13080909, 2.56228467, + 3.18354286, 1.55118873, 0.60703854], + [3.51052921, 6.40861592, 4.52017217, + 1.08859183, 5.93003437, 7.38501235], + [8.57725314, 5.32329684, 4.3160815 , + 6.77629621, 7.49288644, 9.25124645], + [8.0475847 , 6.4163836 , 3.62396899, + 5.04696816, 8.86944591, 0.08356193], + [3.0936211 , 3.95499815, 5.96360799, + 0.6605126 , 3.86466816, 9.6719478 ], + [1.96541611, 5.36489548, 3.70399066, + 1.12084339, 7.96734339, 3.24741959], + [6.14385189, 1.83113282, 2.0741162 , + 4.83746891, 3.2362641 , 0.17227037], + [1.67604863, 5.59126381, 2.09007632, + 5.46945561, 9.71342408, 1.6268169 ], + [9.61221953, 5.27643349, 2.7304223 , + 7.30143289, 2.37011787, 7.83264795], + [9.90983548, 5.79429284, 3.23532663, + 9.17294685, 5.97553658, 4.4452637 ], + [0.35683435, 4.30303963, 1.90266857, + 4.17023551, 4.06186979, 2.7965879 ], + [0.83376644, 9.7300586 , 6.05654181, + 5.09992868, 8.97882656, 4.6280693 ], + [8.48222031, 0.12398542, 7.06665329, + 7.90802805, 8.97576213, 1.98215824], + [0.45055814, 9.90187098, 0.7664873 , + 4.74533329, 8.9772245 , 6.48808184], + [0.88080893, 5.19986295, 6.27365681, + 4.46247949, 4.42179112, 0.3344096 ]] + sn_arr = [[2.55346379, 2.1727307 , 2.35093766, + 2.86042287, 2.97285019, 1.15426079], + [2.46390195, 1.23075187, 2.14164701, + 2.51344824, 1.6329798 , 1.92597048], + [2.3922998 , 1.6673905 , 2.07136149, + 0.94995374, 2.99152438, 2.43278461], + [1.80446599, 1.61757723, 2.21397598, + 0.57976386, 0.73221374, 0.33431018], + [1.95255757, 0.48024565, 1.00356077, + 1.31780419, 2.83845215, 0.95183676], + [2.22861903, 1.5311806 , 1.34814435, + 1.52928755, 0.98082474, 0.47157412], + [1.62612892, 2.34657869, 1.91522882, + 0.72757474, 2.26079876, 2.54277753], + [2.75012934, 2.10862141, 1.86320656, + 2.40689549, 2.57871696, 2.88240738], + [2.65355845, 2.34820047, 1.67626069, + 2.02740696, 2.83043378, 0.07848485], + [1.50175408, 1.76784557, 2.25184422, + 0.49088564, 1.76128254, 2.95377289], + [1.11351605, 2.11819761, 1.69879401, + 0.74379173, 2.66806122, 1.56651489], + [2.27489977, 1.07879322, 1.17346325, + 1.97695394, 1.58172378, 0.15582134], + [0.99700377, 2.16962471, 1.17948397, + 2.1258517 , 2.97456874, 0.99432725], + [2.93003466, 2.09778509, 1.40429925, + 2.51155634, 1.30074808, 2.62807274], + [2.97983427, 2.21480833, 1.56288192, + 2.85500604, 2.27064999, 1.89628347], + [0.29069255, 1.859838 , 1.10727576, + 1.80798637, 1.81437105, 1.42590517], + [0.59196924, 2.96349975, 2.27193322, + 2.03998311, 2.84951285, 1.94222486], + [2.73304224, 0.11442983, 2.48042035, + 2.62755231, 2.84898002, 1.13831482], + [0.35613537, 2.99207933, 0.56870727, + 1.95439588, 2.8492343 , 2.36312065], + [0.61808354, 2.07996798, 2.3182313 , + 1.8836441 , 1.90782557, 0.28421464]] + expected_array = {SiPMCharge.raw : raw_arr, + SiPMCharge.signal_to_noise: sn_arr} + expected_single = {SiPMCharge.raw : [ 99.6473, 93.8179, + 81.3852, 92.4639, + 125.2603, 75.8990], + SiPMCharge.signal_to_noise: [ 9.6651, 9.3927, + 8.7087, 9.3396, + 10.9941, 8.4203]} + + return (S2(times, bin_widths, pmts, sipms), chan_slice, + expected_single, expected_array) + + +@mark.parametrize("charge_type", SiPMCharge) +def test_sipm_charge_array(charge_type , + s2_peak , + signal_to_noise_6400): + s2_peak, chan_slice, _, expected_array = s2_peak + charge_tpl = s2_peak.sipm_charge_array(signal_to_noise_6400, + charge_type , + single_point = False) + + all_wf = s2_peak.sipms.all_waveforms + orig_zeros = np.count_nonzero(all_wf) + charge_arr = np.array(charge_tpl) + calc_zeros = np.count_nonzero(charge_arr) + assert charge_arr.shape == all_wf.T.shape + assert calc_zeros == orig_zeros + + calc_slice = np.array(charge_tpl)[:, chan_slice] + assert_allclose(calc_slice, expected_array[charge_type]) + + +@mark.parametrize("charge_type", SiPMCharge) +def test_sipm_charge_array_single(charge_type , + s2_peak , + signal_to_noise_6400): + s2_peak, chan_slice, expected_single, _ = s2_peak + charge_tpl = s2_peak.sipm_charge_array(signal_to_noise_6400, + charge_type , + single_point = True) + + assert charge_tpl.shape == s2_peak.sipms.ids.shape + orig_zeros = np.count_nonzero(s2_peak.sipms.sum_over_times) + calc_zeros = np.count_nonzero(charge_tpl) + assert calc_zeros == orig_zeros + + assert_allclose(charge_tpl[chan_slice], expected_single[charge_type], + atol=5e-5) diff --git a/invisible_cities/filters/trigger_filters.py b/invisible_cities/filters/trigger_filters.py index a9532ae266..961a1343c3 100644 --- a/invisible_cities/filters/trigger_filters.py +++ b/invisible_cities/filters/trigger_filters.py @@ -20,9 +20,6 @@ def trigger_filter(peak_data : '{channel_no: s2}'): for channel_no, s2s in peak_data.items(): for peak_number, peak in enumerate(s2s): - # print(peak.width) - # print(peak.total_energy) - # print(peak.height) if min_width < peak.width <= max_width: if min_charge < peak.total_energy <= max_charge: if min_height < peak.height <= max_height: diff --git a/invisible_cities/io/dst_io.py b/invisible_cities/io/dst_io.py index 71167760fa..23ce1d5583 100644 --- a/invisible_cities/io/dst_io.py +++ b/invisible_cities/io/dst_io.py @@ -27,7 +27,7 @@ def load_dst(filename, group, node): def load_dsts(dst_list, group, node): dsts = [load_dst(filename, group, node) for filename in dst_list] - return pd.concat(dsts) + return pd.concat(dsts, ignore_index=True) def _make_tabledef(column_types : pd.Series, str_col_length : int=32) -> dict: tabledef = {} diff --git a/invisible_cities/io/dst_io_test.py b/invisible_cities/io/dst_io_test.py index 9b02df3b70..db9cb96d98 100644 --- a/invisible_cities/io/dst_io_test.py +++ b/invisible_cities/io/dst_io_test.py @@ -2,6 +2,7 @@ import string import pandas as pd import tables as tb +import numpy as np from ..core.testing_utils import assert_dataframes_close from ..core.testing_utils import assert_dataframes_equal from ..core.exceptions import TableMismatch @@ -61,10 +62,12 @@ def test_load_dsts_double_file(KrMC_kdst): tbl = KrMC_kdst[0].file_info df_true = KrMC_kdst[0].true df_read = load_dsts([tbl.filename]*2, tbl.group, tbl.node) - df_true = pd.concat([df_true, df_true]) + df_true = pd.concat([df_true, df_true], ignore_index=True) assert_dataframes_close(df_read, df_true , False , rtol=1e-5) + #assert index number unique (important for saving it to pytable) + assert all(~df_read.index.duplicated()) def test_load_dsts_reads_good_kdst(ICDATADIR): good_file = "Kr83_nexus_v5_03_00_ACTIVE_7bar_10evts_KDST.h5" diff --git a/invisible_cities/io/hits_io.py b/invisible_cities/io/hits_io.py index 86c1939736..83ba1f3daf 100644 --- a/invisible_cities/io/hits_io.py +++ b/invisible_cities/io/hits_io.py @@ -39,15 +39,19 @@ def hits_from_df (dst : pd.DataFrame, skip_NN : bool = False) -> Dict[int, HitCo continue hit = Hit(row.npeak, Cluster(row.Q, xy(row.X, row.Y), xy(row.Xrms**2, row.Yrms**2), - row.nsipm, row.Z, row.E, Qc = getattr(row, 'Qc', -1)), - row.Z, row.E, xy(getattr(row, 'Xpeak', -1000) , getattr(row, 'Ypeak', -1000)), - s2_energy_c = getattr(row, 'Ec', -1), track_id = getattr(row, 'track_id', -1)) + row.nsipm, row.Z, row.E, + Qc = getattr(row, 'Qc', -1)), # for backwards compatibility + row.Z, row.E, + xy(getattr(row, 'Xpeak', -1000), # for backwards compatibility + getattr(row, 'Ypeak', -1000)), # for backwards compatibility + s2_energy_c = getattr(row, 'Ec' , -1), # for backwards compatibility + track_id = getattr(row, 'track_id', -1), # for backwards compatibility + Ep = getattr(row, "Ep" , -1)) # for backwards compatibility hits.append(hit) - if len(hits)>0: - all_events.update({event : HitCollection(event, time)}) - all_events[event].hits.extend(hits) + if len(hits): + all_events[event] = HitCollection(event, time, hits=hits) return all_events diff --git a/invisible_cities/io/mcinfo_io.py b/invisible_cities/io/mcinfo_io.py index 760b034370..28155125bd 100644 --- a/invisible_cities/io/mcinfo_io.py +++ b/invisible_cities/io/mcinfo_io.py @@ -1,6 +1,7 @@ import tables as tb import numpy as np +import pandas as pd from .. reco import tbl_functions as tbl from .. core import system_of_units as units @@ -15,8 +16,11 @@ from .. evm.nh5 import MCHitInfo from .. evm.nh5 import MCParticleInfo +from .. database import load_db as DB + from typing import Mapping from typing import Sequence +from typing import Tuple # use Mapping (duck type) rather than dict @@ -168,12 +172,64 @@ def load_mchits(file_name: str, event_range=(0, int(1e9))) -> Mapping[int, Sequence[MCHit]]: with tb.open_file(file_name, mode='r') as h5in: - mcevents = read_mcinfo(h5in, event_range) - mchits_dict = compute_mchits_dict(mcevents) + mchits_dict = read_mchit_info(h5in, event_range) return mchits_dict +def load_mchits_df(file_name : str) -> pd.DataFrame: + """ + Opens file and calls read_mchits_df + file_name : str + The name of the file to be read + """ + extents = pd.read_hdf(file_name, 'MC/extents') + with tb.open_file(file_name) as h5in: + hits = read_mchits_df(h5in, extents) + + return hits + + +def read_mchits_df(h5in : tb.file.File, + extents : pd.DataFrame) -> pd.DataFrame: + """ + Loads the MC hit information into a pandas DataFrame. + h5in : pytables file + A file already read into a pytables object + extents : pd.DataFrame + The extents table from the file which gives + information about the event tracture. + """ + + hits_tb = h5in.root.MC.hits + + # Generating hits DataFrame + hits = pd.DataFrame({'hit_id' : hits_tb.col('hit_indx'), + 'particle_id' : hits_tb.col('particle_indx'), + 'label' : hits_tb.col('label').astype('U13'), + 'time' : hits_tb.col('hit_time'), + 'x' : hits_tb.col('hit_position')[:, 0], + 'y' : hits_tb.col('hit_position')[:, 1], + 'z' : hits_tb.col('hit_position')[:, 2], + 'energy' : hits_tb.col('hit_energy')}) + + evt_hit_df = extents[['last_hit', 'evt_number']] + evt_hit_df.set_index('last_hit', inplace = True) + + hits = hits.merge(evt_hit_df , + left_index = True, + right_index = True, + how = 'left') + hits.rename(columns={"evt_number": "event_id"}, inplace = True) + hits.event_id.fillna(method='bfill', inplace = True) + hits.event_id = hits.event_id.astype(int) + + # Setting the indexes + hits.set_index(['event_id', 'particle_id', 'hit_id'], inplace=True) + + return hits + + def load_mcparticles(file_name: str, event_range=(0, int(1e9))) -> Mapping[int, Mapping[int, MCParticle]]: @@ -181,6 +237,68 @@ def load_mcparticles(file_name: str, return read_mcinfo(h5in, event_range) +def load_mcparticles_df(file_name: str) -> pd.DataFrame: + """ + Opens file and calls read_mcparticles_df + file_name : str + The name of the file to be read + """ + extents = pd.read_hdf(file_name, 'MC/extents') + with tb.open_file(file_name, mode='r') as h5in: + particles = read_mcparticles_df(h5in, extents) + + return particles + + +def read_mcparticles_df(h5in : tb.file.File, + extents : pd.DataFrame) -> pd.DataFrame: + """ + A reader for the MC particle output based + on pandas DataFrames. + + file_name: string + Name of the file to be read + """ + p_tb = h5in.root.MC.particles + + # Generating parts DataFrame + parts = pd.DataFrame({'particle_id' : p_tb.col('particle_indx'), + 'particle_name' : p_tb.col('particle_name').astype('U20'), + 'primary' : p_tb.col('primary').astype('bool'), + 'mother_id' : p_tb.col('mother_indx'), + 'initial_x' : p_tb.col('initial_vertex')[:, 0], + 'initial_y' : p_tb.col('initial_vertex')[:, 1], + 'initial_z' : p_tb.col('initial_vertex')[:, 2], + 'initial_t' : p_tb.col('initial_vertex')[:, 3], + 'final_x' : p_tb.col('final_vertex')[:, 0], + 'final_y' : p_tb.col('final_vertex')[:, 1], + 'final_z' : p_tb.col('final_vertex')[:, 2], + 'final_t' : p_tb.col('final_vertex')[:, 3], + 'initial_volume' : p_tb.col('initial_volume').astype('U20'), + 'final_volume' : p_tb.col('final_volume').astype('U20'), + 'initial_momentum_x': p_tb.col('momentum')[:, 0], + 'initial_momentum_y': p_tb.col('momentum')[:, 1], + 'initial_momentum_z': p_tb.col('momentum')[:, 2], + 'kin_energy' : p_tb.col('kin_energy'), + 'creator_proc' : p_tb.col('creator_proc').astype('U20')}) + + # Adding event info + evt_part_df = extents[['last_particle', 'evt_number']] + evt_part_df.set_index('last_particle', inplace = True) + parts = parts.merge(evt_part_df , + left_index = True, + right_index = True, + how = 'left') + parts.rename(columns={"evt_number": "event_id"}, inplace = True) + parts.event_id.fillna(method='bfill', inplace = True) + parts.event_id = parts.event_id.astype(int) + + # Setting the indexes + parts.set_index(['event_id', 'particle_id'], inplace=True) + + return parts + + def load_mcsensor_response(file_name: str, event_range=(0, int(1e9))) -> Mapping[int, Mapping[int, Waveform]]: @@ -188,8 +306,64 @@ def load_mcsensor_response(file_name: str, return read_mcsns_response(h5in, event_range) -def read_mcinfo_evt (mctables: (tb.Table, tb.Table, tb.Table, tb.Table), - event_number: int, last_row=0) -> ([tb.Table], [tb.Table], [tb.Table]): +def get_sensor_binning(file_name : str) -> Tuple: + """ + Looks in the configuration table of the + input file and extracts the binning used + for both types of sensitive detector. + """ + config = pd.read_hdf(file_name, 'MC/configuration') + bins = config[config.param_key.str.contains('time_binning')] + pmt_bin = bins.param_value[bins.param_key.str.contains('Pmt')].iloc[0] + pmt_bin = float(pmt_bin.split()[0]) * units_dict[pmt_bin.split()[1]] + sipm_bin = bins.param_value[bins.param_key.str.contains('SiPM')].iloc[0] + sipm_bin = float(sipm_bin.split()[0]) * units_dict[sipm_bin.split()[1]] + + return pmt_bin, sipm_bin + + +def load_mcsensor_response_df(file_name : str, + db_file : str, + run_no : int) -> Tuple: + """ + A reader for the MC sensor output based + on pandas DataFrames. + + file_name: string + Name of the file to be read + db_file : string + Name of the detector database to be accessed + run_no : int + Run number for database access + """ + pmt_ids = DB.DataPMT(db_file, run_no).SensorID + + pmt_bin, sipm_bin = get_sensor_binning(file_name) + + extents = pd.read_hdf(file_name, 'MC/extents') + + sns = pd.read_hdf(file_name, 'MC/waveforms') + evt_sns = extents[['last_sns_data', 'evt_number']] + evt_sns.set_index('last_sns_data', inplace = True) + + sns = sns.merge(evt_sns , + left_index = True, + right_index = True, + how = 'left') + sns.evt_number.fillna(method='bfill', inplace = True) + + sns['time'] = sns[sns.sensor_id.isin(pmt_ids)].time_bin * pmt_bin + sns.time.fillna(sns.time_bin * sipm_bin, inplace = True) + + sns.evt_number = sns.evt_number.astype(int) + sns.rename(columns = {'evt_number': 'event_id'}, inplace = True) + sns.set_index(['event_id', 'sensor_id', 'time_bin'], inplace = True) + + return extents.evt_number.unique(), pmt_bin, sipm_bin, sns + + +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]): h5extents = mctables[0] h5hits = mctables[1] h5particles = mctables[2] @@ -210,16 +384,18 @@ def read_mcinfo_evt (mctables: (tb.Table, tb.Table, tb.Table, tb.Table), previous_row = h5extents[iext-1] ihit = int(previous_row['last_hit']) + 1 - ipart = int(previous_row['last_particle']) + 1 + if not return_only_hits: + ipart = int(previous_row['last_particle']) + 1 ihit_end = this_row['last_hit'] - ipart_end = this_row['last_particle'] - if len(h5hits) != 0: while ihit <= ihit_end: hit_rows.append(h5hits[ihit]) ihit += 1 + if return_only_hits: break + + ipart_end = this_row['last_particle'] while ipart <= ipart_end: particle_rows.append(h5particles[ipart]) ipart += 1 @@ -235,7 +411,7 @@ def read_mcinfo_evt (mctables: (tb.Table, tb.Table, tb.Table, tb.Table), def read_mcinfo(h5f, event_range=(0, int(1e9))) -> Mapping[int, Mapping[int, Sequence[MCParticle]]]: mc_info = tbl.get_mc_info(h5f) - + h5extents = mc_info.extents events_in_file = len(h5extents) @@ -292,6 +468,35 @@ def compute_mchits_dict(mcevents:Mapping[int, Mapping[int, MCParticle]]) -> Mapp return mchits_dict +def read_mchit_info(h5f, event_range=(0, int(1e9))) -> Mapping[int, Sequence[MCHit]]: + """Returns all hits in the event""" + mc_info = tbl.get_mc_info(h5f) + h5extents = mc_info.extents + events_in_file = len(h5extents) + + all_events = {} + + for iext in range(*event_range): + if iext >= events_in_file: + break + + current_event = {} + evt_number = h5extents[iext]['evt_number'] + hit_rows, _, _ = read_mcinfo_evt(mc_info, evt_number, iext, True) + + hits = [] + for h5hit in hit_rows: + hit = MCHit(h5hit['hit_position'], + h5hit['hit_time'], + h5hit['hit_energy'], + h5hit['label'].decode('utf-8','ignore')) + hits.append(hit) + + all_events[evt_number] = hits + + return all_events + + def read_mcsns_response(h5f, event_range=(0, 1e9)) -> Mapping[int, Mapping[int, Waveform]]: h5config = h5f.root.MC.configuration diff --git a/invisible_cities/io/mcinfo_io_test.py b/invisible_cities/io/mcinfo_io_test.py index 533de8f38c..65c6278155 100644 --- a/invisible_cities/io/mcinfo_io_test.py +++ b/invisible_cities/io/mcinfo_io_test.py @@ -1,18 +1,26 @@ import os import numpy as np +import pandas as pd import tables as tb -from glob import glob -from os.path import expandvars +from glob import glob +from os.path import expandvars +from numpy.testing import assert_allclose from . mcinfo_io import load_mchits +from . mcinfo_io import load_mchits_df +from . mcinfo_io import read_mchits_df from . mcinfo_io import load_mcparticles +from . mcinfo_io import load_mcparticles_df +from . mcinfo_io import read_mcparticles_df +from . mcinfo_io import get_sensor_binning from . mcinfo_io import load_mcsensor_response +from . mcinfo_io import load_mcsensor_response_df from . mcinfo_io import mc_info_writer from . mcinfo_io import read_mcinfo_evt from .. core import system_of_units as units -from ..core.exceptions import NoParticleInfoInFile +from .. core.exceptions import NoParticleInfoInFile from .. reco.tbl_functions import get_mc_info @@ -193,6 +201,34 @@ def test_load_mchits(mc_particle_and_hits_nexus_data): assert np.allclose(t, ht) +def test_read_mchits_df(mc_particle_and_hits_nexus_data): + efile, _, _, _, _, _, _, X, Y, Z, E, t = mc_particle_and_hits_nexus_data + + with tb.open_file(efile) as h5in: + extents = pd.read_hdf(efile, 'MC/extents') + hit_df = read_mchits_df(h5in, extents) + + evt = 0 + assert np.allclose(X, hit_df.loc[evt].x .values) + assert np.allclose(Y, hit_df.loc[evt].y .values) + assert np.allclose(Z, hit_df.loc[evt].z .values) + assert np.allclose(E, hit_df.loc[evt].energy.values) + assert np.allclose(t, hit_df.loc[evt].time .values) + + +def test_load_mchits_df(mc_particle_and_hits_nexus_data): + efile, _, _, _, _, _, _, X, Y, Z, E, t = mc_particle_and_hits_nexus_data + + hit_df = load_mchits_df(efile) + + evt = 0 + assert np.allclose(X, hit_df.loc[evt].x .values) + assert np.allclose(Y, hit_df.loc[evt].y .values) + assert np.allclose(Z, hit_df.loc[evt].z .values) + assert np.allclose(E, hit_df.loc[evt].energy.values) + assert np.allclose(t, hit_df.loc[evt].time .values) + + def test_load_mcparticles(mc_particle_and_hits_nexus_data): efile, name, vi, vf, p, Ep, nhits, X, Y, Z, E, t = mc_particle_and_hits_nexus_data @@ -218,6 +254,58 @@ def test_load_mcparticles(mc_particle_and_hits_nexus_data): assert np.allclose(t, ht) +def test_load_mcparticles_df(mc_particle_and_hits_nexus_data): + efile, name, vi, vf, p, k_eng, *_ = mc_particle_and_hits_nexus_data + + mcparticle_df = load_mcparticles_df(efile) + + evt = 0 + p_id = 1 + particle = mcparticle_df.loc[evt].loc[p_id] + assert particle.particle_name == name + assert np.isclose(particle.kin_energy, k_eng) + + init_vtx = particle[['initial_x', 'initial_y', + 'initial_z', 'initial_t']] + assert_allclose(init_vtx.tolist(), vi) + + fin_vtx = particle[['final_x', 'final_y', + 'final_z', 'final_t']] + assert_allclose(fin_vtx.tolist(), vf) + + init_mom = particle[['initial_momentum_x', + 'initial_momentum_y', + 'initial_momentum_z']] + assert_allclose(init_mom.tolist(), p) + + +def test_read_mcparticles_df(mc_particle_and_hits_nexus_data): + efile, name, vi, vf, p, k_eng, *_ = mc_particle_and_hits_nexus_data + + with tb.open_file(efile) as h5in: + extents = pd.read_hdf(efile, 'MC/extents') + mcparticle_df = read_mcparticles_df(h5in, extents) + + evt = 0 + p_id = 1 + particle = mcparticle_df.loc[evt].loc[p_id] + assert particle.particle_name == name + assert np.isclose(particle.kin_energy, k_eng) + + init_vtx = particle[['initial_x', 'initial_y', + 'initial_z', 'initial_t']] + assert_allclose(init_vtx.tolist(), vi) + + fin_vtx = particle[['final_x', 'final_y', + 'final_z', 'final_t']] + assert_allclose(fin_vtx.tolist(), vf) + + init_mom = particle[['initial_momentum_x', + 'initial_momentum_y', + 'initial_momentum_z']] + assert_allclose(init_mom.tolist(), p) + + def test_load_sensors_data(mc_sensors_nexus_data): efile, pmt0_first, pmt0_last, pmt0_tot_samples, sipm_id, sipm = mc_sensors_nexus_data @@ -241,6 +329,46 @@ def test_load_sensors_data(mc_sensors_nexus_data): assert np.allclose(samples, sipm) +def test_get_sensor_binning(mc_sensors_nexus_data): + fullsim_data, *_ = mc_sensors_nexus_data + + sensor_binning = 100 * units.nanosecond, 1 * units.microsecond + + binning = get_sensor_binning(fullsim_data) + + assert len(binning) == 2 + assert binning[0] in sensor_binning + assert binning[1] in sensor_binning + + +def test_load_mcsensor_response_df(mc_sensors_nexus_data): + efile, pmt0_first, pmt0_last, pmt0_tot_samples, sipm_id, sipm = mc_sensors_nexus_data + + ## Should be generalised for other detectors but data + ## available at the moment is for new. + ## Should consider adding to test data! + list_evt, _, _, wfs = load_mcsensor_response_df(efile, 'new', -6400) + + ## Check first pmt + pmt0_id = 0 + wf = wfs.loc[list_evt[0]].loc[pmt0_id] + n_samp = len(wf) + indx0 = (wf.index[0], wf.iloc[0].charge) + indx_last = (wf.index[n_samp - 1], wf.iloc[n_samp - 1].charge) + + assert n_samp == pmt0_tot_samples + assert indx0 == pmt0_first + assert indx_last == pmt0_last + + ## Check chosen SiPM + sipm_wf = wfs.loc[list_evt[0]].loc[sipm_id] + bin_q = [(sipm_wf.index[i], sipm_wf.iloc[i].charge) + for i in range(len(sipm_wf))] + + assert np.all(bin_q == sipm) + + + def test_read_last_sensor_response(mc_sensors_nexus_data): efile, _, _, _, _, _ = mc_sensors_nexus_data diff --git a/invisible_cities/reco/corrections_new.py b/invisible_cities/reco/corrections_new.py index 4273dddfe5..a444fbb535 100644 --- a/invisible_cities/reco/corrections_new.py +++ b/invisible_cities/reco/corrections_new.py @@ -9,10 +9,11 @@ from typing import Optional from enum import auto -from .. core import system_of_units as units -from .. core.exceptions import TimeEvolutionTableMissing -from .. types.ic_types import AutoNameEnumBase -from .. evm.event_model import Hit +from .. core.core_functions import in_range +from .. core import system_of_units as units +from .. core.exceptions import TimeEvolutionTableMissing +from .. types.ic_types import AutoNameEnumBase +from .. evm.event_model import Hit @dataclass @@ -93,12 +94,20 @@ def maps_coefficient_getter(mapinfo : Series, for a given (X,Y) position """ - binsx = np.linspace(mapinfo.xmin,mapinfo.xmax,mapinfo.nx+1) - binsy = np.linspace(mapinfo.ymin,mapinfo.ymax,mapinfo.ny+1) + binsx = np.linspace(mapinfo.xmin, mapinfo.xmax, mapinfo.nx + 1) + binsy = np.linspace(mapinfo.ymin, mapinfo.ymax, mapinfo.ny + 1) + def get_maps_coefficient(x : np.array, y : np.array) -> np.array: - ix = np.digitize(x, binsx)-1 - iy = np.digitize(y, binsy)-1 - return np.array([map_df.get(j, {}).get(i, np.nan) for i, j in zip(iy,ix)]) + ix = np.digitize(x, binsx) - 1 + iy = np.digitize(y, binsy) - 1 + + valid = in_range(x, mapinfo.xmin, mapinfo.xmax) + valid &= in_range(y, mapinfo.ymin, mapinfo.ymax) + output = np.full(len(valid), np.nan, dtype=np.float) + + output[valid] = map_df.values[iy[valid], ix[valid]] + return output + return get_maps_coefficient @@ -168,6 +177,7 @@ def time_coefs_corr(time_evt : np.array, par_factor = par_i/par_mean return par_factor + def get_df_to_z_converter(map_te: ASectorMap) -> Callable: """ For given map, it returns a function that provides the conversion @@ -202,6 +212,7 @@ class norm_strategy(AutoNameEnumBase): kr = auto() custom = auto() + def get_normalization_factor(map_e0 : ASectorMap, norm_strat: norm_strategy = norm_strategy.max, norm_value: Optional[float] = None @@ -242,6 +253,7 @@ def get_normalization_factor(map_e0 : ASectorMap, return norm_value + def apply_all_correction_single_maps(map_e0 : ASectorMap, map_lt : ASectorMap, map_te : Optional[ASectorMap] = None, diff --git a/invisible_cities/reco/hits_functions.py b/invisible_cities/reco/hits_functions.py index b6588cec31..e208d2ade7 100644 --- a/invisible_cities/reco/hits_functions.py +++ b/invisible_cities/reco/hits_functions.py @@ -31,7 +31,7 @@ def merge_NN_hits(hits : List[evm.Hit], same_peak : bool = True) -> List[evm.Hit z_closest = min(abs(h.Z-nn_h.Z) for h in hits_to_merge) except ValueError: continue - h_closest = [h for h in hits_to_merge if abs(h.Z-nn_h.Z)==z_closest] + h_closest = [h for h in hits_to_merge if np.isclose(abs(h.Z-nn_h.Z), z_closest)] en_tot = sum([h.E for h in h_closest]) for h in h_closest: hits_to_correct.append([h,nn_h.E*(h.E/en_tot)]) diff --git a/invisible_cities/reco/paolina_functions.py b/invisible_cities/reco/paolina_functions.py index b30369efd7..e6d80076b0 100644 --- a/invisible_cities/reco/paolina_functions.py +++ b/invisible_cities/reco/paolina_functions.py @@ -307,74 +307,69 @@ def drop_voxel(voxels: Sequence[Voxel], the_vox: Voxel) -> int: """Eliminate an individual voxel from a set of voxels and give its energy to the hit that is closest to the barycenter of the eliminated voxel hits, provided that it belongs to a neighbour voxel.""" - - ### be sure that the voxel to be eliminated has at least one neighbour - ### beyond itself - the_neighbours = np.array([neighbours(the_vox, v) for v in voxels]) - if len(the_neighbours[the_neighbours>0]) <= 1: - return 0 - - ### remove voxel from list of voxels - voxels.remove(the_vox) + the_neighbour_voxels = [v for v in voxels if neighbours(the_vox, v)] pos = [h.pos for h in the_vox.hits] qs = [getattr(h, e_type) for h in the_vox.hits] + #if there are no hits associated to voxels the pos will be an empty list if len(pos) == 0: - min_dist = 1e+06 - min_v = voxels[0] - for v in voxels: - if neighbours(the_vox, v): - dist = np.linalg.norm(the_vox.pos - v.pos) - if dist < min_dist: - min_dist = dist - min_v = v + min_dist = min(np.linalg.norm(the_vox.pos-v.pos) for v in the_neighbour_voxels) + min_v = [v for v in the_neighbour_voxels if np.isclose(np.linalg.norm(the_vox.pos-v.pos), min_dist)] - ### add voxel energy to voxel - min_v.E += the_vox.E - - return 1 + ### add dropped voxel energy to closest voxels, proportional to the voxels energy + sum_en_v = sum(v.E for v in min_v) + for v in min_v: + v.E += the_vox.E/sum_en_v * v.E + return bary_pos = np.average(pos, weights=qs, axis=0) ### find hit with minimum distance, only among neighbours - min_dist = 1e+06 - min_hit = voxels[0].hits[0] - min_v = voxels[0] - for v in voxels: - if neighbours(the_vox, v): - for hh in v.hits: - dist = np.linalg.norm(bary_pos - hh.pos) - if dist < min_dist: - min_dist = dist - min_hit = hh - min_v = v - - ### add voxel energy to hit and to voxel, separately - setattr(min_hit, e_type, getattr(min_hit, e_type) + the_vox.E) - min_v.E += the_vox.E - - return 1 - - mod_voxels = copy.deepcopy(voxels) - - while True: - n_modified_voxels = 0 + min_dist = min(np.linalg.norm(bary_pos-hh.pos) for v in the_neighbour_voxels for hh in v.hits) + min_h_v = [(h, v) for v in the_neighbour_voxels for h in v.hits if np.isclose(np.linalg.norm(bary_pos-h.pos), min_dist)] + min_hs = set(h for (h,v) in min_h_v) + min_vs = set(v for (h,v) in min_h_v) + + ### add dropped voxel energy to closest hits/voxels, proportional to the hits/voxels energy + sum_en_h = sum(getattr(h, e_type) for h in min_hs) + sum_en_v = sum(v.E for v in min_vs) + for h in min_hs: + setattr(h, e_type, getattr(h, e_type) + getattr(h, e_type) * the_vox.E/sum_en_h) + for v in min_vs: + v.E = sum(getattr(h, e_type) for h in v.hits) + + def nan_energy(voxel): + voxel.E = np.nan + for hit in voxel.hits: + setattr(hit, e_type, np.nan) + + mod_voxels = copy.deepcopy(voxels) + dropped_voxels = [] + + modified = True + while modified: + modified = False trks = make_track_graphs(mod_voxels) for t in trks: if len(t.nodes()) < min_vxls: continue - extr1, extr2 = find_extrema(t) - if extr1.E < energy_threshold: - n_mod = drop_voxel(mod_voxels, extr1) - n_modified_voxels += n_mod - if extr2.E < energy_threshold: - n_mod = drop_voxel(mod_voxels, extr2) - n_modified_voxels += n_mod + for extreme in find_extrema(t): + if extreme.E < energy_threshold: + ### be sure that the voxel to be eliminated has at least one neighbour + ### beyond itself + n_neighbours = sum(neighbours(extreme, v) for v in mod_voxels) + if n_neighbours > 1: + mod_voxels .remove(extreme) + dropped_voxels.append(extreme) + drop_voxel(mod_voxels, extreme) + nan_energy(extreme) + modified = True + + return mod_voxels, dropped_voxels - if n_modified_voxels == 0: - break - return mod_voxels +def get_track_energy(track): + return sum([vox.Ehits for vox in track.nodes()]) diff --git a/invisible_cities/reco/paolina_functions_test.py b/invisible_cities/reco/paolina_functions_test.py index 26fa731051..b7a7477464 100644 --- a/invisible_cities/reco/paolina_functions_test.py +++ b/invisible_cities/reco/paolina_functions_test.py @@ -8,6 +8,7 @@ import networkx as nx from itertools import combinations +from operator import attrgetter from numpy.testing import assert_almost_equal @@ -51,6 +52,7 @@ from . paolina_functions import Contiguity from . paolina_functions import drop_end_point_voxels from . paolina_functions import make_tracks +from . paolina_functions import get_track_energy from .. core.exceptions import NoHits from .. core.exceptions import NoVoxels @@ -86,8 +88,17 @@ def hit(draw, min_value=1, max_value=100): E_c = draw(floats ( 50, 100)) x_peak = draw(floats (-10, 10)) y_peak = draw(floats (-10, 10)) + track = draw(integers( 0, 10)) + E_p = draw(floats ( 50, 100)) - return Hit(npeak,Cluster(Q,xy(x,y),xy(xvar,yvar),nsipm),z,E,xy(x_peak,y_peak),s2_energy_c=E_c) + return Hit(npeak, + Cluster(Q, xy(x,y), xy(xvar,yvar), nsipm), + z, + E, + xy(x_peak, y_peak), + s2_energy_c = E_c, + track_id = track, + Ep = E_p) @composite @@ -374,16 +385,16 @@ def test_voxelize_single_hit(): @fixture(scope='module') def voxels_without_hits(): voxel_spec = ((10,10,10,1), - (10,10,11,1), - (10,10,12,1), - (10,10,13,1), - (10,10,14,1), - (10,10,15,1), - (10,11,15,1), - (10,12,15,1), - (10,13,15,1), - (10,14,15,1), - (10,15,15,1) + (10,10,11,2), + (10,10,12,2), + (10,10,13,2), + (10,10,14,2), + (10,10,15,2), + (10,11,15,2), + (10,12,15,2), + (10,13,15,2), + (10,14,15,2), + (10,15,15,2) ) voxels = [Voxel(x,y,z, E, np.array([1,1,1])) for (x,y,z,E) in voxel_spec] @@ -491,7 +502,7 @@ def test_energy_is_conserved_with_dropped_voxels(hits, requested_voxel_dimension energies = [v.E for v in voxels] e_thr = min(energies) + fraction_zero_one * (max(energies) - min(energies)) - mod_voxels = drop_end_point_voxels(voxels, e_thr, min_voxels) + mod_voxels, _ = drop_end_point_voxels(voxels, e_thr, min_voxels) tot_final_energy = sum(v.E for v in mod_voxels) final_trks = make_track_graphs(mod_voxels) final_trk_energies = [sum(vox.E for vox in t.nodes()) for t in final_trks] @@ -501,6 +512,78 @@ def test_energy_is_conserved_with_dropped_voxels(hits, requested_voxel_dimension assert np.allclose(ini_trk_energies, final_trk_energies) +@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) + for voxel in dropped_voxels: + assert np.isnan(voxel.E) + for hit in voxel.hits: + assert np.isnan(getattr(hit, energy_type.value)) + + +@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_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] + e_thr = min(energies) + fraction_zero_one * (max(energies) - min(energies)) + mod, drop = drop_end_point_voxels(voxels, e_thr, min_voxels) + new_voxels = sorted(mod + drop, key=attrgetter("xyz")) + + for e_type in HitEnergy: + if e_type is energy_type: continue + + for v_before, v_after in zip(voxels, new_voxels): + for h_before, h_after in zip(v_before.hits, v_after.hits): + #assert sum(energy_from_hits(v_before, e_type.value)) == sum(energy_from_hits(v_after, e_type.value)) + assert np.isclose(getattr(h_before, e_type.value), getattr(h_after, e_type.value)) + + +@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_drop_voxels_voxel_energy_is_sum_of_hits_general(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)) + mod_voxels, _ = drop_end_point_voxels(voxels, e_thr, min_voxels) + for v in mod_voxels: + assert np.isclose(v.E, sum(getattr(h, energy_type.value) for h in v.hits)) + + + +@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_drop_end_point_voxels_constant_number_of_voxels_and_hits(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)) + new_voxels = drop_end_point_voxels(voxels, e_thr, min_voxels) + (mod_voxels, + dropped_voxels) = new_voxels + assert len(mod_voxels) + len(dropped_voxels) == len(voxels) + assert sum(1 for vs in new_voxels for v in vs for h in v.hits) == len(hits) + + def test_initial_voxels_are_the_same_after_dropping_voxels(ICDATADIR): # Get some test data: nothing interesting to see here @@ -517,7 +600,7 @@ def test_initial_voxels_are_the_same_after_dropping_voxels(ICDATADIR): # This is the core of the test: collect data before/after ... ante_energies = [v.E for v in voxels] ante_positions = [v.XYZ for v in voxels] - mod_voxels = drop_end_point_voxels(voxels, e_thr, min_voxels) + mod_voxels, _ = drop_end_point_voxels(voxels, e_thr, min_voxels) post_energies = [v.E for v in voxels] post_positions = [v.XYZ for v in voxels] @@ -549,7 +632,7 @@ def test_tracks_with_dropped_voxels(ICDATADIR): ini_energies = [sum(vox.E for vox in t.nodes()) for t in ini_trks] ini_n_voxels = np.array([len(t.nodes()) for t in ini_trks]) - mod_voxels = drop_end_point_voxels(voxels, e_thr, min_voxels) + mod_voxels, _ = drop_end_point_voxels(voxels, e_thr, min_voxels) trks = make_track_graphs(mod_voxels) n_of_tracks = len(trks) @@ -566,17 +649,58 @@ def test_tracks_with_dropped_voxels(ICDATADIR): assert np.all(ini_n_voxels - n_voxels == expected_diff_n_voxels) +def test_drop_voxels_deterministic(ICDATADIR): + hit_file = os.path.join(ICDATADIR, 'tracks_0000_6803_trigger2_v0.9.9_20190111_krth1600.h5') + evt_number = 19 + e_thr = 5867.92 + min_voxels = 3 + vox_size = [15.] * 3 + + all_hits = load_hits(hit_file) + hits = all_hits[evt_number].hits + voxels = voxelize_hits(hits, vox_size, strict_voxel_size=False) + mod_voxels , _ = drop_end_point_voxels(sorted(voxels, key = lambda v:v.E, reverse = False), e_thr, min_voxels) + mod_voxels_r, _ = drop_end_point_voxels(sorted(voxels, key = lambda v:v.E, reverse = True ), e_thr, min_voxels) + + for v1, v2 in zip(sorted(mod_voxels, key = lambda v:v.E), sorted(mod_voxels_r, key = lambda v:v.E)): + assert np.isclose(v1.E, v2.E) + + def test_voxel_drop_in_short_tracks(): hits = [BHit(10, 10, 10, 1), BHit(26, 10, 10, 1)] voxels = voxelize_hits(hits, [15,15,15], strict_voxel_size=True) e_thr = sum(v.E for v in voxels) + 1. min_voxels = 0 - mod_voxels = drop_end_point_voxels(voxels, e_thr, min_voxels) + mod_voxels, _ = drop_end_point_voxels(voxels, e_thr, min_voxels) assert len(mod_voxels) >= 1 +def test_drop_voxels_voxel_energy_is_sum_of_hits(): + def make_hit(x, y, z, e): + return Hit(peak_number = 0, + cluster = Cluster(0, xy(x, y), xy(0, 0), 1), + z = z, + s2_energy = e, + peak_xy = xy(0, 0), + Ep = e) + + # Create a track with an extreme to be dropped and two hits at the same + # distance from the barycenter of the the voxel to be dropped with + # different energies in the hits *and* in the voxels + voxels = [Voxel( 0, 0, 0, 0.1, size=5, e_type=HitEnergy.Ep, hits = [make_hit( 0, 0, 0, 0.1) ]), + Voxel( 5, -5, 0, 1.0, size=5, e_type=HitEnergy.Ep, hits = [make_hit( 5, -5, 0, 0.7), make_hit( 5, -8, 0, 0.3)]), + Voxel( 5, 5, 0, 1.5, size=5, e_type=HitEnergy.Ep, hits = [make_hit( 5, 5, 0, 0.9), make_hit( 5, 8, 0, 0.6)]), + Voxel(10, 0, 0, 2.0, size=5, e_type=HitEnergy.Ep, hits = [make_hit(10, 5, 0, 1.2), make_hit(11, 5, 0, 0.8)]), + Voxel(15, 0, 0, 2.5, size=5, e_type=HitEnergy.Ep, hits = [make_hit(15, 0, 0, 1.8), make_hit(11, 0, 0, 0.7)]), + Voxel(20, 0, 0, 3.0, size=5, e_type=HitEnergy.Ep, hits = [make_hit(20, 0, 0, 1.5), make_hit(11, 0, 0, 1.5)])] + + modified_voxels, _ = drop_end_point_voxels(voxels, energy_threshold = 0.5, min_vxls = 1) + for v in modified_voxels: + assert np.isclose(v.E, sum(h.Ep for h in v.hits)) + + @parametrize('radius, expected', ((10., ( 60, 20)), (12., ( 60, 60)), @@ -654,7 +778,7 @@ def test_paolina_functions_with_voxels_without_associated_hits(blob_radius, min_ energies = [v.E for v in voxels] e_thr = min(energies) + fraction_zero_one * (max(energies) - min(energies)) - mod_voxels = drop_end_point_voxels(voxels, e_thr, min_voxels) + mod_voxels, _ = drop_end_point_voxels(voxels, e_thr, min_voxels) trks = make_track_graphs(mod_voxels) for t in tracks: a, b = find_extrema(t) @@ -663,11 +787,12 @@ def test_paolina_functions_with_voxels_without_associated_hits(blob_radius, min_ assert np.allclose(blob_centre(b), b.pos) -@given(bunch_of_corrected_hits(), box_sizes, radius, fraction_zero_one) -def test_paolina_functions_with_hit_energy_different_from_default_value(hits, requested_voxel_dimensions, blob_radius, fraction_zero_one): - - energy_type = HitEnergy.Ec - +@mark.parametrize("energy_type", (HitEnergy.Ec, HitEnergy.Ep)) +@given(hits = bunch_of_corrected_hits(), + requested_voxel_dimensions = box_sizes, + blob_radius = radius, + fraction_zero_one = fraction_zero_one) +def test_paolina_functions_with_hit_energy_different_from_default_value(hits, requested_voxel_dimensions, blob_radius, fraction_zero_one, energy_type): voxels = voxelize_hits(hits, requested_voxel_dimensions, strict_voxel_size=False) voxels_c = voxelize_hits(hits, requested_voxel_dimensions, strict_voxel_size=False, energy_type=energy_type) @@ -682,10 +807,10 @@ def test_paolina_functions_with_hit_energy_different_from_default_value(hits, re energies_c = [v.E for v in voxels_c] e_thr = min(energies_c) + fraction_zero_one * (max(energies_c) - min(energies_c)) # Test that this function doesn't fail - mod_voxels_c = drop_end_point_voxels(voxels_c, e_thr, min_vxls=0) + mod_voxels_c, _ = drop_end_point_voxels(voxels_c, e_thr, min_vxls=0) - tot_energy = sum(getattr(h, energy_type.value) for v in voxels_c for h in v.hits) - tot_mod_energy = sum(getattr(h, energy_type.value) for v in mod_voxels_c for h in v.hits) + tot_energy = sum(getattr(h, energy_type.value) for v in voxels_c for h in v.hits) + tot_mod_energy = sum(getattr(h, energy_type.value) for v in mod_voxels_c for h in v.hits) assert np.isclose(tot_energy, tot_mod_energy) @@ -739,3 +864,11 @@ def test_make_tracks_function(ICDATADIR): tc_blob_energies = (tc.blobs[0].E, tc.blobs[1].E) assert np.allclose(blob_energies(t, blob_radius), tc_blob_energies) + + +@given(bunch_of_hits, box_sizes) +def test_make_voxel_graph_keeps_all_voxels(hits, voxel_dimensions): + voxels = voxelize_hits (hits , voxel_dimensions) + tracks = make_track_graphs(voxels) + # assert sum of track energy equal to sum of hits energies + assert_almost_equal(sum(get_track_energy(track) for track in tracks), sum(h.E for h in hits)) diff --git a/invisible_cities/reco/peak_functions.py b/invisible_cities/reco/peak_functions.py index dc0a47f1d1..e5d13c5f16 100644 --- a/invisible_cities/reco/peak_functions.py +++ b/invisible_cities/reco/peak_functions.py @@ -36,22 +36,23 @@ def split_in_peaks(indices, stride): return np.split(indices, where + 1) -def select_peaks(peaks, time, length): +def select_peaks(peaks, time, length, pmt_samp_wid=25*units.ns): def is_valid(indices): - return (time .contains(indices[ 0] * 25 * units.ns) and - time .contains(indices[-1] * 25 * units.ns) and + return (time .contains(indices[ 0] * pmt_samp_wid) and + time .contains(indices[-1] * pmt_samp_wid) and length.contains(indices[-1] + 1 - indices[0])) return tuple(filter(is_valid, peaks)) def pick_slice_and_rebin(indices, times, widths, - wfs, rebin_stride, pad_zeros=False): + wfs, rebin_stride, pad_zeros=False, + sipm_pmt_bin_ratio=40): slice_ = slice(indices[0], indices[-1] + 1) times_ = times [ slice_] widths_ = widths[ slice_] wfs_ = wfs [:, slice_] if pad_zeros: - n_miss = indices[0] % 40 + n_miss = indices[0] % sipm_pmt_bin_ratio n_wfs = wfs.shape[0] times_ = np.concatenate([np.zeros( n_miss) , times_]) widths_ = np.concatenate([np.zeros( n_miss) , widths_]) @@ -63,12 +64,14 @@ def pick_slice_and_rebin(indices, times, widths, def build_pmt_responses(indices, times, widths, ccwf, - pmt_ids, rebin_stride, pad_zeros): + pmt_ids, rebin_stride, pad_zeros, + sipm_pmt_bin_ratio): (pk_times , pk_widths, pmt_wfs ) = pick_slice_and_rebin(indices, times, widths, ccwf , rebin_stride, - pad_zeros = pad_zeros) + pad_zeros = pad_zeros, + sipm_pmt_bin_ratio = sipm_pmt_bin_ratio) return pk_times, pk_widths, PMTResponses(pmt_ids, pmt_wfs) @@ -87,17 +90,23 @@ def build_peak(indices, times, widths, ccwf, pmt_ids, rebin_stride, with_sipms, Pk, - sipm_wfs = None, - thr_sipm_s2 = 0): + pmt_samp_wid = 25 * units.ns, + sipm_samp_wid = 1 * units.mus, + sipm_wfs = None, + thr_sipm_s2 = 0): + sipm_pmt_bin_ratio = int(sipm_samp_wid/pmt_samp_wid) (pk_times , pk_widths, pmt_r ) = build_pmt_responses(indices, times, widths, ccwf, pmt_ids, - rebin_stride, pad_zeros = with_sipms) + rebin_stride, pad_zeros = with_sipms, + sipm_pmt_bin_ratio = sipm_pmt_bin_ratio) if with_sipms: - sipm_r = build_sipm_responses(indices // 40, times // 40, - widths * 40, sipm_wfs, - rebin_stride // 40, + sipm_r = build_sipm_responses(indices // sipm_pmt_bin_ratio, + times // sipm_pmt_bin_ratio, + widths * sipm_pmt_bin_ratio, + sipm_wfs, + rebin_stride // sipm_pmt_bin_ratio, thr_sipm_s2) else: sipm_r = SiPMResponses.build_empty_instance() @@ -109,14 +118,16 @@ def find_peaks(ccwfs, index, time, length, stride, rebin_stride, Pk, pmt_ids, + pmt_samp_wid = 25*units.ns, + sipm_samp_wid = 1*units.mus, sipm_wfs=None, thr_sipm_s2=0): ccwfs = np.array(ccwfs, ndmin=2) peaks = [] - times = np.arange (ccwfs.shape[1]) * 25 * units.ns - widths = np.full (ccwfs.shape[1], 25 * units.ns) + times = np.arange (ccwfs.shape[1]) * pmt_samp_wid + widths = np.full (ccwfs.shape[1], pmt_samp_wid) indices_split = split_in_peaks(index, stride) - selected_splits = select_peaks (indices_split, time, length) + selected_splits = select_peaks (indices_split, time, length, pmt_samp_wid) with_sipms = Pk is S2 and sipm_wfs is not None for indices in selected_splits: @@ -124,17 +135,23 @@ def find_peaks(ccwfs, index, widths, ccwfs, pmt_ids, rebin_stride, with_sipms, Pk, + pmt_samp_wid, sipm_samp_wid, sipm_wfs, thr_sipm_s2) peaks.append(pk) return peaks def get_pmap(ccwf, s1_indx, s2_indx, sipm_zs_wf, - s1_params, s2_params, thr_sipm_s2, pmt_ids): - return PMap(find_peaks(ccwf, s1_indx, Pk=S1, pmt_ids=pmt_ids, **s1_params), + s1_params, s2_params, thr_sipm_s2, pmt_ids, + pmt_samp_wid, sipm_samp_wid): + return PMap(find_peaks(ccwf, s1_indx, Pk=S1, pmt_ids=pmt_ids, + pmt_samp_wid=pmt_samp_wid, + **s1_params), find_peaks(ccwf, s2_indx, Pk=S2, pmt_ids=pmt_ids, - sipm_wfs = sipm_zs_wf, - thr_sipm_s2 = thr_sipm_s2, + sipm_wfs = sipm_zs_wf, + thr_sipm_s2 = thr_sipm_s2, + pmt_samp_wid = pmt_samp_wid, + sipm_samp_wid = sipm_samp_wid, **s2_params)) @@ -142,7 +159,7 @@ def rebin_times_and_waveforms(times, widths, waveforms, rebin_stride=2, slices=None): if rebin_stride < 2: return times, widths, waveforms - if not slices: + if slices is None: n_bins = int(np.ceil(len(times) / rebin_stride)) reb = rebin_stride slices = [slice(reb * i, reb * (i + 1)) for i in range(n_bins)] @@ -152,11 +169,15 @@ def rebin_times_and_waveforms(times, widths, waveforms, rebinned_widths = np.zeros( len(slices) ) rebinned_wfs = np.zeros((n_sensors, len(slices))) - for i, s in enumerate(slices): - t = times [ s] - e = waveforms[:, s] - w = np.sum(e, axis=0) if np.any(e) else None + for i, sl in enumerate(slices): + t = times [ sl] + e = waveforms[:, sl] + ## Weight with the charge sum per slice + ## if positive and unweighted if all + ## negative. + s = np.sum(e, axis=0).clip(0) + w = s if np.any(s) else None rebinned_times [ i] = np.average(t, weights=w) - rebinned_widths[ i] = np.sum ( widths[s]) + rebinned_widths[ i] = np.sum ( widths[sl]) rebinned_wfs [:, i] = np.sum (e, axis=1) return rebinned_times, rebinned_widths, rebinned_wfs diff --git a/invisible_cities/reco/peak_functions_test.py b/invisible_cities/reco/peak_functions_test.py index 292567031d..bea7a6cae1 100644 --- a/invisible_cities/reco/peak_functions_test.py +++ b/invisible_cities/reco/peak_functions_test.py @@ -1,3 +1,5 @@ +import os + import numpy as np from pytest import approx @@ -22,6 +24,7 @@ from ..evm .pmaps import S1 from ..evm .pmaps import S2 from ..evm .pmaps import PMap +from ..io .pmaps_io import load_pmaps from ..types.ic_types_c import minmax from . import peak_functions as pf @@ -334,7 +337,7 @@ def test_build_pmt_responses(wf_with_indices): times, widths, wfs, indices = wf_with_indices ids = np.arange(wfs.shape[0]) ts, wid, pmt_r = pf.build_pmt_responses(indices, times, widths, - wfs, ids, 1, False) + wfs, ids, 1, False, 40) assert ts == approx (times[indices]) assert wid == approx (widths[indices]) assert pmt_r.ids == exactly(ids) @@ -495,11 +498,14 @@ def test_get_pmap(s1_and_s2_with_indices): pmt_ids = np.arange( pmt_wfs.shape[0]) sipm_ids = np.arange(sipm_wfs.shape[0]) - + pmt_samp_wid = 25 * units.ns + sipm_samp_wid = 1 * units.mus pmap = pf.get_pmap(pmt_wfs, s1_indx, s2_indx, sipm_wfs, s1_params, s2_params, thr_sipm_s2 = -1, - pmt_ids = pmt_ids) + pmt_ids = pmt_ids, + pmt_samp_wid = pmt_samp_wid , + sipm_samp_wid = sipm_samp_wid) (rebinned_times , rebinned_widths, @@ -597,3 +603,14 @@ def test_rebin_times_and_waveforms_times_are_consistent(t_and_wf, stride): np.ones_like(wfs), stride) assert np.sum(rb_times) * stride == approx(np.sum(times)) + + +def test_rebin_times_and_waveforms_negative_bins(ICDATADIR): + # Tests that rebin doesn't raise an error + # even with multiple adjacent negative bins + pmap_file = os.path.join(ICDATADIR, 'pmaps_negative_bins.h5') + pmaps = load_pmaps(pmap_file) + s2 = pmaps[48490].s2s[0] + rebinned_info = pf.rebin_times_and_waveforms(s2.times , + s2.bin_widths , + s2.pmts.all_waveforms) diff --git a/manage.sh b/manage.sh index b246372744..6b375e09aa 100644 --- a/manage.sh +++ b/manage.sh @@ -103,6 +103,7 @@ dependencies: - flaky = 3.4.0 - hypothesis = 3.68.0 - pytest-xdist = 1.23.2 +- coverage = 4.5.4 - pip: - pytest-instafail==0.4.0 EOF