Add Valdrada, the trigger simulation city

invisible_cities/cities/components.py

@@ -125,7 +125,7 @@ def create_timestamp(rate: float) -> float:
 
     Returns
     -------
-    Function to calculate timestamp for the given rate with 
+    Function to calculate timestamp for the given rate with
     event_number as parameter.
     """
 
@@ -142,12 +142,12 @@ def create_timestamp(rate: float) -> float:
     def create_timestamp_(event_number: Union[int, float]) -> float:
         """
         Calculates timestamp for a given Event Number and Rate.
-    
+
         Parameters
         ----------
         event_number : Union[int, float]
                        ID value of the current event.
-    
+
         Returns
         -------
         Calculated timestamp : float
@@ -326,6 +326,33 @@ def deconv_pmt(RWF):
     return deconv_pmt
 
 
+def deconv_pmt_fpga(dbfile     : str
+                   ,run_number : int
+                   ,selection  : List[int] = None
+                   ) -> Callable:
+    '''
+    Apply deconvolution as done in the FPGA.
+
+    Parameters
+    ----------
+    dbfile     : Database to be used
+    run_number : Run number of the database
+    selection  : List of PMT IDs to apply deconvolution to.
+
+    Returns
+    ----------
+    deconv_pmt : Function that will apply the deconvolution.
+    '''
+    DataPMT    = load_db.DataPMT(dbfile, run_number = run_number)
+    pmt_active = np.nonzero(DataPMT.Active.values)[0].tolist() if selection is None else selection
+    coeff_c    = DataPMT.coeff_c  .values.astype(np.double)[pmt_active]
+    coeff_blr  = DataPMT.coeff_blr.values.astype(np.double)[pmt_active]
+    def deconv_pmt(RWF):
+        return zip(*map(blr.deconvolve_signal_fpga, RWF[pmt_active],
+                        coeff_c                   , coeff_blr      ))
+    return deconv_pmt
+
+
 def get_run_number(h5in):
     if   "runInfo" in h5in.root.Run: return h5in.root.Run.runInfo[0]['run_number']
     elif "RunInfo" in h5in.root.Run: return h5in.root.Run.RunInfo[0]['run_number']

invisible_cities/cities/valdrada.py

@@ -0,0 +1,172 @@
+"""
+-----------------------------------------------------------------------
+                                 valdrada
+-----------------------------------------------------------------------
+
+This city finds simulates the trigger procedure at the FPGA over PMT waveforms.
+This includes a number of tasks:
+    - Use a moving average to obtain the baseline
+    - Truncate to integers the deconvolution output
+    - Find trigger candidates and compute their characteristics, taking into
+      account the unique circumstances derived from working on a continuous scheme
+      (for example, delay in signals)
+    - Evaluate coincidences between trigger candidates.
+
+A number of general configuration parameters are needed (input as a dict, trigger_config):
+- coincidence_window : Time window (in time bins) in which valid trigger coincidences are counted
+- discard_width      : Any trigger with width less than this parameter is discarded.
+- multipeak          : Dictionary with multipeak protection parameters, otherwise None.
+     - q_min      : Integrated charge threshold of a post-trigger peak to discard
+                    a trigger due to multipeak protection. In ADC counts.
+     - time_min   : Minimum width of a post-trigger peak to discard a trigger due to
+                    multipeak protection. In time bins.
+     - time_after : For how long is the multipeak protection evaluated after a
+                    valid trigger. In mus.
+
+A individual trigger configuration can be given per channel, through a dict
+with keys equal to PMT IDs, which marks the validity range of the peak
+characteristics:
+- q_min, q_max          : Range for the integrated charge of the peak (q_min < q < q_max).
+                          In ADC counts.
+- time_min, time_max    : Range for the peak width (time_min <= width < time_max).
+                          In time mus.
+- baseline_dev, amp_max : Range for peak height (baseline_dev < height < amp_max).
+                          In ADC counts.
+- pulse_valid_ext       : Time allowed for a pulse to go below baseline_dev. In ns.
+
+The result of the city is a dataframe containing the event ID and PMT ID of each
+trigger candidate. For each trigger candidate a number of parameters is computed:
+ - trigger_time : time bin at which the trigger candidate starts.
+ - q            : integrated ADC counts within the peak.
+ - width        : Length of the peak in time bins.
+ - height       : Maximum ADC values in a given time bien within the peak.
+ - baseline     : Value of the baseline at trigger_time.
+ - max_height   : Maximum (minimum due to PMT negative signals) height in the wvf.
+
+Additionally, a set of flags are assigned depending on wether the parameters are
+within range of the trigger configuration:
+ - valid_q    : If peak is within q range.
+ - valid_w    : If peak is within width range.
+ - valid_h    : If peak is within height range.
+ - valid_peak : Only if 'multipeak' is active, True if there isn't a post-trigger
+               candidate with the configuration parameters.
+ - valid_all  : boolean and of previous valid flags.
+
+ Finally, a series of coincidence-related values are also given:
+ - n_coinc       : Number of valid triggers within the coincidence_window,
+                   starting at the trigger trigger_time and including the trigger itself.
+                   -1 indicates no valid trigger (including the trigger itself).
+ - closest_ttime : Time difference to closest valid trigger.
+                   -1 if there are none aside from the trigger itself.
+ - closest_pmt   : PMT ID of the closest valid trigger.
+                   -1 if there are none aside from the trigger itself.
+"""
+
+import tables as tb
+import numpy  as np
+
+from .. reco                  import tbl_functions        as tbl
+from .. io  .run_and_event_io import run_and_event_writer
+from .. io  .trigger_io       import       trigger_writer, trigger_dst_writer
+
+from .. dataflow            import dataflow as fl
+from .. dataflow.dataflow   import push
+from .. dataflow.dataflow   import pipe
+from .. dataflow.dataflow   import sink
+
+from .  components import city
+from .  components import print_every
+from .  components import collect
+from .  components import copy_mc_info
+from .  components import deconv_pmt_fpga
+from .  components import WfType
+from .  components import wf_from_files
+from .  components import get_number_of_active_pmts
+
+from .. reco.trigger_functions import retrieve_trigger_information
+from .. reco.trigger_functions import get_trigger_candidates
+from .. reco.trigger_functions import check_trigger_coincidence
+
+from .. core       import system_of_units as units
+
+
+def check_empty_trigger(triggers) -> bool:
+    """
+    Filter for valdrada flow.
+    The flow stops if there are no candidate triggers/
+    """
+    return len(triggers) > 0
+
+@city
+def valdrada(files_in, file_out, compression, event_range, print_mod, detector_db, run_number,
+             trigger_config = dict(), channel_config = dict()):
+
+    # Change time units into number of waveform bins.
+    if trigger_config['multipeak'] is not None:
+        trigger_config['multipeak']['time_min'  ] /= 25*units.ns
+        trigger_config['multipeak']['time_after'] /= 25*units.ns
+    for k in channel_config.keys():
+        channel_config[k]['time_min']       /= 25*units.ns
+        channel_config[k]['time_max']       /= 25*units.ns
+        channel_config[k]['pulse_valid_ext'] = int(channel_config[k]['pulse_valid_ext']/25*units.ns)
+
+    #### Define data transformations
+    # Raw WaveForm to deconvolved waveforms
+    rwf_to_cwf           = fl.map(deconv_pmt_fpga(detector_db, run_number, list(channel_config.keys())),
+                                  args = "pmt",
+                                  out  = ("cwf", "baseline"))
+
+    # Extract all possible trigger candidates of each PMT
+    trigger_on_channels  = fl.map(retrieve_trigger_information(channel_config, trigger_config),
+                                  args = ("pmt", "cwf", "baseline", "event_number"),
+                                  out  = "triggers")
+
+    # Add coincidence between channels
+    check_coincidences   = fl.map(check_trigger_coincidence(trigger_config['coincidence_window']),
+                                  item = "triggers")
+
+    # Filter events with zero triggers
+    filter_empty_trigger = fl.map(check_empty_trigger,
+                                  args = "triggers",
+                                  out  = "empty_trigger")
+
+    event_count_in     = fl.spy_count()
+    event_count_out    = fl.spy_count()
+    events_no_triggers = fl.count_filter(bool, args = "empty_trigger")
+    evtnum_collect     = collect()
+
+    with tb.open_file(file_out, "w", filters = tbl.filters(compression)) as h5out:
+        # Define writers...
+        write_event_info_   = run_and_event_writer(h5out)
+        write_trigger_info_ = trigger_writer      (h5out, get_number_of_active_pmts(detector_db, run_number))
+        write_trigger_dst_  = trigger_dst_writer  (h5out)
+        # ... and make them sinks
+
+        write_event_info   = sink(write_event_info_  , args=( "run_number"  , "event_number"    , "timestamp"   ))
+        write_trigger_info = sink(write_trigger_info_, args=( "trigger_type", "trigger_channels"                ))
+        write_trigger_dst  = sink(write_trigger_dst_ , args=  "triggers"                                         )
+
+        result = push(source = wf_from_files(files_in, WfType.rwf),
+                      pipe   = pipe(fl.slice(*event_range, close_all=True),
+                                    print_every(print_mod),
+                                    event_count_in.spy,
+                                    rwf_to_cwf,
+                                    trigger_on_channels,
+                                    filter_empty_trigger,
+                                    events_no_triggers.filter,
+                                    check_coincidences,
+                                    event_count_out.spy,
+                                    fl.branch("event_number", evtnum_collect.sink),
+                                    fl.fork(write_trigger_dst ,
+                                            write_event_info  ,
+                                            write_trigger_info)),
+                      result = dict(events_in   = event_count_in    .future,
+                                    events_out  = event_count_out   .future,
+                                    evtnum_list = evtnum_collect    .future,
+                                    events_pass = events_no_triggers.future))
+
+        if run_number <= 0:
+            copy_mc_info(files_in, h5out, result.evtnum_list,
+                         detector_db, run_number)
+
+        return result

invisible_cities/cities/valdrada_test.py

@@ -0,0 +1,53 @@
+import os
+import tables as tb
+
+from .  valdrada           import valdrada
+from .. core.testing_utils import assert_tables_equality
+
+
+def test_valdrada_contains_all_tables(trigger_config):
+    conf, PATH_OUT = trigger_config
+    valdrada(**conf)
+    with tb.open_file(PATH_OUT) as h5out:
+        assert "Run"             in h5out.root
+        assert "Run/events"      in h5out.root
+        assert "Run/runInfo"     in h5out.root
+        assert "Trigger"         in h5out.root
+        assert "Trigger/events"  in h5out.root
+        assert "Trigger/trigger" in h5out.root
+        assert "Trigger/DST"     in h5out.root
+
+
+def test_valdrada_exact_result_multipeak(ICDATADIR, trigger_config):
+    true_out          = os.path.join(ICDATADIR, "exact_result_multipeak_valdrada.h5")
+    conf, PATH_OUT   = trigger_config
+    valdrada(**conf)
+
+    tables = ("Trigger/events", "Trigger/trigger", "Trigger/DST",
+              "Run/events"      , "Run/runInfo")
+
+    with tb.open_file(true_out)  as true_output_file:
+        with tb.open_file(PATH_OUT) as      output_file:
+            for table in tables:
+                assert hasattr(output_file.root, table)
+                got      = getattr(     output_file.root, table)
+                expected = getattr(true_output_file.root, table)
+                assert_tables_equality(got, expected)
+
+
+def test_valdrada_exact_result(ICDATADIR, trigger_config):
+    true_out          = os.path.join(ICDATADIR, "exact_result_valdrada.h5")
+    conf, PATH_OUT    = trigger_config
+    conf['trigger_config']['multipeak'] = None
+    valdrada(**conf)
+
+    tables = ("Trigger/events", "Trigger/trigger", "Trigger/DST",
+              "Run/events"      , "Run/runInfo")
+
+    with tb.open_file(true_out)  as true_output_file:
+        with tb.open_file(PATH_OUT) as      output_file:
+            for table in tables:
+                assert hasattr(output_file.root, table)
+                got      = getattr(     output_file.root, table)
+                expected = getattr(true_output_file.root, table)
+                assert_tables_equality(got, expected)

invisible_cities/config/valdrada.conf

@@ -0,0 +1,19 @@
+files_in    = '$ICDIR/database/test_data/blr_fpga_examples.h5'
+file_out   = '$ICDIR/database/test_data/valdrada.h5'
+
+compression = 'ZLIB4'
+run_number  = 8093
+detector_db = 'new'
+event_range = all
+
+trigger_config = {'coincidence_window':64, 'discard_width':40,
+                  'multipeak':{'q_min':100000, 'time_min':2*mus, 'time_after':800*mus}}
+
+channel_config = {0: {'q_min'       :5000, 'q_max'   :50000,
+                      'time_min'    :2*mus, 'time_max': 40*mus,
+                      'baseline_dev':  10, 'amp_max' : 1000,
+                      'pulse_valid_ext':50*ns},
+                  2: {'q_min'       :5000, 'q_max'   :50000,
+                      'time_min'    :2*mus, 'time_max': 40*mus,
+                      'baseline_dev':  10, 'amp_max' : 1000,
+                      'pulse_valid_ext':50*ns}}

invisible_cities/conftest.py

@@ -50,6 +50,11 @@ def example_blr_wfs_filename(ICDATADIR):
     return os.path.join(ICDATADIR, "blr_examples.h5")
 
 
+@pytest.fixture(scope='session')
+def example_blr_fpga_wfs_filename(ICDATADIR):
+    return os.path.join(ICDATADIR, "blr_fpga_examples.h5")
+
+
 @pytest.fixture(scope  = 'session',
                 params = ['electrons_40keV_z250_MCRD.h5'])
 def electron_MCRD_file(request, ICDATADIR):
@@ -699,6 +704,41 @@ def deconvolution_config(ICDIR, ICDATADIR, PSFDIR, config_tmpdir):
 
     return conf, PATH_OUT
 
+
+@pytest.fixture(scope='function') # Needs to be function as the config dict is modified when running
+def trigger_config(ICDIR, ICDATADIR, config_tmpdir):
+    PATH_IN     = os.path.join(ICDATADIR    ,    "blr_fpga_examples.h5")
+    PATH_OUT    = os.path.join(config_tmpdir,     "trigger_city_out.h5")
+    nevt_req    = 10
+    conf        = dict(files_in       = PATH_IN ,
+                       file_out       = PATH_OUT,
+                       event_range    = nevt_req,
+                       compression    = 'ZLIB4',
+                       print_mod      = 1000,
+                       run_number     = 8093,
+                       trigger_config = dict(coincidence_window = 64
+                                            ,discard_width      = 40
+                                            ,multipeak          = dict(q_min      = 100000
+                                                                      ,time_min   = 2  *units.mus
+                                                                      ,time_after = 800*units.mus)),
+                       channel_config = {0 : dict(q_min           = 5000
+                                                 ,q_max           = 50000
+                                                 ,time_min        = 2 *units.mus
+                                                 ,time_max        = 40*units.mus
+                                                 ,baseline_dev    = 10
+                                                 ,amp_max         = 1000
+                                                 ,pulse_valid_ext = 50*units.ns)
+                                        ,2 : dict(q_min           = 5000
+                                                 ,q_max           = 50000
+                                                 ,time_min        = 2 *units.mus
+                                                 ,time_max        = 40*units.mus
+                                                 ,baseline_dev    = 10
+                                                 ,amp_max         = 1000
+                                                 ,pulse_valid_ext = 50*units.ns)})
+
+    return conf, PATH_OUT
+
+
 ## To make very slow tests only run with specific option
 def pytest_addoption(parser):
     parser.addoption(

invisible_cities/evm/nh5.py

@@ -186,3 +186,22 @@ class VoxelsTable(tb.IsDescription):
 class EventPassedFilter(tb.IsDescription):
     event  = tb.Int32Col(pos=0)
     passed = tb. BoolCol(pos=1)
+
+
+class TriggerTable(tb.IsDescription):
+    event         = tb.UInt32Col(pos= 0)
+    pmt           = tb.UInt32Col(pos= 1)
+    trigger_time  = tb.UInt32Col(pos= 2)
+    q             = tb.UInt32Col(pos= 3)
+    width         = tb.UInt32Col(pos= 4)
+    height        = tb.UInt32Col(pos= 5)
+    valid_q       = tb.BoolCol  (pos= 6)
+    valid_w       = tb.BoolCol  (pos= 7)
+    valid_h       = tb.BoolCol  (pos= 8)
+    valid_peak    = tb.BoolCol  (pos= 9)
+    valid_all     = tb.BoolCol  (pos=10)
+    baseline      = tb.Int32Col (pos=11)
+    max_height    = tb.Int32Col (pos=12)
+    n_coinc       = tb.Int32Col (pos=13)
+    closest_ttime = tb.Int32Col (pos=14)
+    closest_pmt   = tb.Int32Col (pos=15)