726 Change rounding in PMT simulation

#726

[author: andLaing]

Changes the way the output of the PMT simulation used in diomira is
converted to integer to better simulate the ADC.

This change should remove or significantly improve the ~15% loss of
signal observed in passing MC through the diomira city.

[reviewer: mmkekic]

This PR fixes ~15% signal loss after passing through diomira. The main
changes are:

- Using `round` instead of `floor` for adc converter simulation (a test is
added to show that for each PMT the mean ADC is equivalent to the
gain).
- Change `daq_decimate` to use `scipy` default (a test is added to show
that the signal is maintained after the down-sampling).
- Change ADC central values.

All the algorithm changes are discussed and approved by
V. Herrero. Good job Andrew in recovering our MC signal!

invisible_cities/cities/diomira.py

@@ -159,7 +159,7 @@ def simulate_pmt_response(pmtrd):
         rwf, blr = sf.simulate_pmt_response(0, pmtrd[np.newaxis],
                                             adc_to_pes, pe_resolution,
                                             detector, run_number)
-        return rwf.astype(np.int16), blr.astype(np.int16)
+        return np.round(rwf).astype(np.int16), np.round(blr).astype(np.int16)
     return simulate_pmt_response
 
 

invisible_cities/cities/diomira_test.py

@@ -10,11 +10,15 @@
 from .. core.configure     import                 configure
 from .. core.testing_utils import   assert_dataframes_close
 from .. core.testing_utils import    assert_tables_equality
+from .. database           import                   load_db
 from .. reco               import             tbl_functions as tbl
 from .. io  .mcinfo_io     import get_event_numbers_in_file
 from .. io  .mcinfo_io     import            load_mchits_df
 from .. io  .mcinfo_io     import       load_mcparticles_df
 
+from .. core                import fit_functions as fitf
+from .. core.core_functions import shift_to_bin_centers
+
 from .  diomira import diomira
 
 
@@ -166,7 +170,7 @@ def test_diomira_exact_result(ICDATADIR, output_tmpdir):
     file_out    = os.path.join(output_tmpdir                                  ,
                                "exact_result_diomira.h5"                      )
     true_output = os.path.join(ICDATADIR                                      ,
-                               "Kr83_nexus_v5_03_00_ACTIVE_7bar_3evts.RWF.h5" )
+                               "Kr83_nexus_v5_03_00_ACTIVE_7bar_3evts.CRWF.h5")
 
     conf = configure("diomira invisible_cities/config/diomira.conf".split())
     conf.update(dict(run_number   = -6340,
@@ -182,12 +186,8 @@ def test_diomira_exact_result(ICDATADIR, output_tmpdir):
     diomira(**conf)
     np.random.set_state(original_random_state)
 
-    ## tables = (     "MC/extents",  "MC/hits"   , "MC/particles", "MC/generators",
-    ##                "RD/pmtrwf" ,  "RD/pmtblr" , "RD/sipmrwf"  ,
-    ##               "Run/events" , "Run/runInfo",
-    ##           "Filters/trigger")
-    tables = (     "RD/pmtrwf" ,  "RD/pmtblr" , "RD/sipmrwf"  ,
-                  "Run/events" , "Run/runInfo",
+    tables = ("RD/pmtrwf"      ,  "RD/pmtblr" , "RD/sipmrwf",
+              "Run/events"     , "Run/runInfo",
               "Filters/trigger")
     with tb.open_file(true_output)  as true_output_file:
         with tb.open_file(file_out) as      output_file:
@@ -204,3 +204,40 @@ def test_diomira_can_fix_random_seed(output_tmpdir):
     conf.update(dict(random_seed = 123),
                 file_out     = file_out)
     diomira(**conf)
+
+
+## to run the following test, use the --runslow option with pytest
+@mark.veryslow
+def test_diomira_reproduces_singlepe(ICDATADIR, output_tmpdir):
+    file_in  = os.path.join(ICDATADIR    ,     'single_pe_pmts.h5')
+    file_out = os.path.join(output_tmpdir, 'single_pe_elec_sim.h5')
+
+    run_no = 7951
+    nevt   =  400
+    conf   = configure("diomira invisible_cities/config/diomira.conf".split())
+    conf.update(dict(files_in     = file_in ,
+                     file_out     = file_out,
+                     run_number   =   run_no,
+                     event_range  =     nevt,
+                     print_mod    =      100,
+                     trigger_type =     None,
+                     random_seed  =     1847))
+    diomira(**conf)
+
+    pmt_gain = load_db.DataPMT('new', run_no).adc_to_pes.values
+    with tb.open_file(file_out) as h5saved:
+        pmt_signal  = h5saved.root.RD.pmtblr
+
+        pmt_sum_adc = np.sum(h5saved.root.RD.pmtblr, axis=2)
+
+        bins        = np.arange(-50, 50, 0.5)
+        bin_centres = np.repeat(shift_to_bin_centers(bins)[np.newaxis, :],
+                                len(pmt_gain), 0)
+        gain_diff   = pmt_gain[:, np.newaxis] - pmt_sum_adc.T
+        histos      = [np.histogram(diffs, bins=bins)[0] for diffs in gain_diff]
+        seeds       = [(x.sum(), x.mean(), x.std(ddof=1)) for x in histos]
+        fits        = tuple(map(fitf.fit, [fitf.gauss]*len(pmt_gain),
+                                bin_centres, histos, seeds))
+        ## show the mean is within error of reproducing 1 pe ADC
+        for fit_res in fits:
+            assert np.abs(fit_res.values[1]) < fit_res.errors[1]

invisible_cities/conftest.py

@@ -707,3 +707,23 @@ def deconvolution_config(ICDIR, ICDATADIR, PSFDIR, config_tmpdir):
                                             inter_method  =    'cubic'))
 
     return conf, PATH_OUT
+
+## To make very slow tests only run with specific option
+def pytest_addoption(parser):
+    parser.addoption(
+        "--runslow", action="store_true", default=False, help="run slow tests"
+        )
+
+
+def pytest_configure(config):
+    config.addinivalue_line("markers", "veryslow: mark test as very slow to run")
+
+
+def pytest_collection_modifyitems(config, items):
+    if config.getoption("--runslow"):
+        # --runslow given in cli: do not skip slow tests
+        return
+    skip_slow = pytest.mark.skip(reason="need --runslow option to run")
+    for item in items:
+        if "veryslow" in item.keywords:
+            item.add_marker(skip_slow)

invisible_cities/sierpe/fee.py

@@ -29,7 +29,7 @@
 f_LPF1 =  3 * units.MHZ
 f_LPF2 = 10 * units.MHZ
 ADC_TO_PES_LPF = 24.1  # After LPF, comes out from spe area
-ADC_TO_PES     = 18.71
+ADC_TO_PES     = 16.241
 OFFSET = 2500  # offset adc
 CEILING = 4096  # ceiling of adc
 
@@ -399,5 +399,13 @@ def daq_decimator(f_sample1, f_sample2, signal_in):
     Includes anti-aliasing filter
     """
 
+    ## The default 30 point FIR filter with Hamming window of order
+    ## 20 * scale is used as an acceptable compromise between charge
+    ## losses and speed as approved by V. Herrero. A slightly better
+    ## charge conservation (~0.001% loss vs ~0.01%) is achieved using:
+    ## signal.decimate(signal_i, scale,
+    ##                 ftype=signal.dlti(*signal.butter(1, 0.8/scale)),
+    ##                 zero_phase=True)
+    ## but for an order 35 increase in processing time.
     scale = int(f_sample1 / f_sample2)
-    return signal.decimate(signal_in, scale, n=30, ftype='fir', zero_phase=False)
+    return signal.decimate(signal_in, scale, ftype='fir', zero_phase=True)

invisible_cities/sierpe/fee_test.py

@@ -1,11 +1,13 @@
-import numpy as np
+import numpy  as np
+import tables as tb
 
 import pytest
 from scipy import signal
 from flaky import flaky
 
-from .. core import system_of_units as units
-from .       import fee             as FE
+from .. core     import system_of_units as units
+from .. database import load_db
+from .           import fee             as FE
 
 def signal_i_th():
     """Generates a "theoretical" current signal (signal_i)"""
@@ -74,7 +76,29 @@ def test_show_signal_decimate_signature():
     assert 23 < adc_to_pes     < 23.1
 
 
+def test_signal_maintained(electron_MCRD_file):
+    """Test that there is no appreciable charge loss in daq_decimate"""
+    detector_db = 'new'
+    run_number  = -7951
+    pmt_id      = 0
+    with tb.open_file(electron_MCRD_file) as h5in:
+        evt0_pmt   = h5in.root.pmtrd[0]
+        adc_to_pes = load_db.DataPMT(detector_db, run_number).adc_to_pes.values
+        spe        = FE.SPE()
+        fee        = FE.FEE(detector_db, run_number,
+                            noise_FEEPMB_rms = FE.NOISE_I  ,
+                            noise_DAQ_rms    = FE.NOISE_DAQ)
 
+        cc         = adc_to_pes[pmt_id] / FE.ADC_TO_PES
+        signal_i   = FE.spe_pulse_from_vector(spe, evt0_pmt[pmt_id], norm=cc)
+        signal_d   = FE.daq_decimator(FE.f_mc, FE.f_sample, signal_i)
+        signal_blr = FE.signal_v_lpf(fee, signal_d) * FE.v_to_adc()
+
+        input_adc = adc_to_pes[pmt_id] * np.sum(evt0_pmt[pmt_id])
+        assert np.isclose(np.sum(signal_blr), input_adc, rtol=2e-4)
+
+
+@pytest.mark.skip('Skipped as uses outdated functions not used in code')
 def test_spe_to_adc():
     """Convert SPE to adc values with the current FEE Parameters must be."""
     ipmt = 0