Selection city

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

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)
 

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):

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")