WIP: TDC measurement features

.github/workflows/codestyle.yml

@@ -28,4 +28,4 @@ jobs:
       - name: Install flake8
         run: pip install flake8
       - name: Check code
-        run: while read in; do flake8 --extend-ignore E501,W503 "$in"; done < <(git diff --name-only --diff-filter=ACMTUXB origin/master -- '*.py')
+        run: while read in; do flake8 --extend-ignore E501,W503 "$in"; done < <(git diff --name-only --diff-filter=ACMTUXB origin/development -- '*.py')

.github/workflows/tests.yml

@@ -12,15 +12,15 @@ on:
 jobs:
   tests:
     name: Python ${{matrix.python-version}} | ${{matrix.sim}} 
-    runs-on: ubuntu-20.04
+    runs-on: ubuntu-22.04
     env:
       SIM: ${{matrix.sim}}
     strategy:
       fail-fast: false
       matrix:
         include:
           - sim: verilator
-            sim-version: v4.106
+            sim-version: v5.020
             python-version: '3.10'
     steps:
     - uses: actions/checkout@v3
@@ -45,7 +45,7 @@ jobs:
     - name: Install Verilator
       if: matrix.sim == 'verilator'
       run: |
-        sudo apt install -y --no-install-recommends make g++ perl python3 autoconf flex bison libfl2 libfl-dev zlibc zlib1g zlib1g-dev
+        sudo apt install -y --no-install-recommends make g++ help2man perl python3 autoconf flex bison libfl2 libfl-dev zlib1g zlib1g-dev
         git clone https://github.com/verilator/verilator.git -b ${{matrix.sim-version}}
         cd verilator
         autoconf

firmware/src/tjmonopix2_core.v

@@ -533,7 +533,8 @@ rrp_arbiter
     .DATA_IN({
         RX_FIFO_DATA,
         TLU_FIFO_DATA,
-        TDC_FIFO_DATA}),
+        TDC_FIFO_DATA
+    }),
     .READ_GRANT({
         RX_FIFO_READ,
         TLU_FIFO_READ,
@@ -612,7 +613,7 @@ pulse_gen #(
 
 // ----- TDC ----- //
 localparam CLKDV = 4;  // division factor from 160 MHz clock to DV_CLK (here 40 MHz)
-wire [CLKDV * 4 - 1:0] FAST_TRIGGER_OUT;
+// wire [CLKDV * 4 - 1:0] FAST_TRIGGER_OUT;
 // wire LEMO_RX0_FROM_TDC;
 // wire HITOR_FROM_TDC;
 
@@ -624,12 +625,12 @@ tdc_s3 #(
     .DATA_IDENTIFIER(4'b0010),
     .FAST_TDC(1),
     .FAST_TRIGGER(1),
-    .BROADCAST(0)         // generate for first TDC module the 640MHz sampled trigger signal and share it with other modules using TRIGGER input
+    .BROADCAST(0)         // generate for LVDS TDC module the 640MHz sampled trigger signal and share it with other modules using TRIGGER input
 ) i_tdc (
     .CLK320(CLK320),      // 320 MHz
     .CLK160(CLK160),      // 160 MHz
     .DV_CLK(CLK40),       // 40 MHz
-    .TDC_IN(LVDS_HITOR),  // HITOR
+    .TDC_IN(LVDS_HITOR),  // LVDS HITOR (DP)
     .TDC_OUT(),
     .TRIG_IN(LEMO_RX[0]),
     .TRIG_OUT(),

tjmonopix2/analysis/analysis.py

@@ -19,7 +19,7 @@
 
 
 class Analysis(object):
-    def __init__(self, raw_data_file=None, analyzed_data_file=None,
+    def __init__(self, raw_data_file=None, analyzed_data_file=None, tot_calib_file=None,
                  store_hits=True, cluster_hits=False, analyze_tdc=False, use_tdc_trigger_dist=False,
                  build_events=False, chunk_size=1000000, **_):
         self.log = logger.setup_derived_logger('Analysis')
@@ -32,6 +32,7 @@ def __init__(self, raw_data_file=None, analyzed_data_file=None,
         self.chunk_size = chunk_size
         self.analyze_tdc = analyze_tdc
         self.use_tdc_trigger_dist = use_tdc_trigger_dist
+        self.tot_calib_file = tot_calib_file
 
         if self.build_events:
             self.cluster_hits = True
@@ -62,6 +63,7 @@ def _get_configs(self):
             self.scan_config = au.ConfigDict(in_file.root.configuration_in.scan.scan_config[:])
             self.chip_settings = au.ConfigDict(in_file.root.configuration_in.chip.settings[:])
             self.tlu_config = au.ConfigDict(in_file.root.configuration_in.bench.TLU[:])
+            self.tdc_config = au.ConfigDict(in_file.root.configuration_in.bench.TDC[:])
 
     def __enter__(self):
         return self
@@ -156,7 +158,7 @@ def _setup_clusterizer(self):
                            ('frame', 'u1'),
                            ('column', 'u2'),
                            ('row', 'u2'),
-                           ('charge', 'u1'),
+                           ('charge', 'u2'),
                            ('timestamp', 'i8')]
         cluster_fields = {'event_number': 'event_number',
                           'column': 'column',
@@ -172,7 +174,7 @@ def _setup_clusterizer(self):
         cluster_description = [('event_number', 'u4'),
                                ('id', '<u2'),
                                ('size', '<u2'),
-                               ('tot', '<u2'),
+                               ('tot', '<u4'),
                                ('seed_col', '<u2'),
                                ('seed_row', '<u2'),
                                ('mean_col', '<f4'),
@@ -254,14 +256,18 @@ def end_of_cluster_function(hits, clusters, cluster_size,
                                          noisy_pixels, disabled_pixels,
                                          seed_hit_index)
 
+            if self.tot_calib_file:
+                max_hit_charge = 2048
+            else:
+                max_hit_charge = 128
             # Initialize clusterizer with custom hit/cluster fields
             self.clz = HitClusterizer(
                 hit_fields=hit_fields,
                 hit_dtype=hit_dtype,
                 cluster_fields=cluster_fields,
                 cluster_dtype=self.cluster_dtype,
                 min_hit_charge=1,
-                max_hit_charge=128,
+                max_hit_charge=max_hit_charge,
                 column_cluster_distance=5,
                 row_cluster_distance=5,
                 frame_cluster_distance=1,
@@ -294,6 +300,9 @@ def analyze_data(self):
                 if self.build_events:
                     trigger_n, trigger_ts, event_n = 0, 0, 0
                     event_table = self._create_table(out_file, name='Hits', title='event_data', dtype=au.event_dtype)
+                if self.tot_calib_file is not None:
+                    with tb.open_file(self.tot_calib_file, 'r') as calib_file:
+                        self.tot_calib = calib_file.root.InjTotCalibration[:]
                 if self.cluster_hits:
                     cluster_table = out_file.create_table(
                         out_file.root, name='Cluster',
@@ -302,11 +311,21 @@ def analyze_data(self):
                         filters=tb.Filters(complib='blosc',
                                            complevel=5,
                                            fletcher32=False))
+                    if self.tot_calib_file:
+                        cs_tot_size = 2048
+                    else:
+                        cs_tot_size = 256
                     hist_cs_size = np.zeros(shape=(30, ), dtype=np.uint32)
-                    hist_cs_tot = np.zeros(shape=(256, ), dtype=np.uint32)
+                    hist_cs_tot = np.zeros(shape=(cs_tot_size, ), dtype=np.uint32)
                     hist_cs_shape = np.zeros(shape=(300, ), dtype=np.int32)
 
-                interpreter = RawDataInterpreter(n_scan_params=n_scan_params, trigger_data_format=self.tlu_config['DATA_FORMAT'])
+                # Setup interpreter
+                interpreter = RawDataInterpreter(
+                    n_scan_params=n_scan_params,
+                    trigger_data_format=self.tlu_config['DATA_FORMAT'],
+                    en_trigger_dist=self.tdc_config['EN_TRIGGER_DIST']
+                )
+
                 self.last_chunk = False
                 pbar = tqdm(total=n_words, unit=' Words', unit_scale=True)
                 upd = 0
@@ -336,6 +355,7 @@ def analyze_data(self):
                         if self.build_events:
                             data_to_clusterizer = event_dat
                         else:
+                            hit_dat = hit_dat[hit_dat['col'] < 1000]  # Can only call tot_calib for hit data
                             hit_data_cs_fmt = np.zeros(len(hit_dat), dtype=au.event_dtype)
                             hit_data_cs_fmt['event_number'][:] = hit_dat['timestamp'][:]
                             hit_data_cs_fmt['trigger_number'][:] = -1
@@ -346,11 +366,19 @@ def analyze_data(self):
                             hit_data_cs_fmt['timestamp'][:] = hit_dat['timestamp'][:]
                             data_to_clusterizer = hit_data_cs_fmt
 
+                        if self.tot_calib_file:
+                            data_to_clusterizer['charge'][:] = au._inv_tot_response_func(
+                                data_to_clusterizer['charge'][:],
+                                self.tot_calib[data_to_clusterizer[:]['column'], data_to_clusterizer[:]['row']][:, 0],
+                                self.tot_calib[data_to_clusterizer[:]['column'], data_to_clusterizer[:]['row']][:, 1],
+                                self.tot_calib[data_to_clusterizer[:]['column'], data_to_clusterizer[:]['row']][:, 2]
+                            )
+
                         _, cluster = self.clz.cluster_hits(data_to_clusterizer)
                         cluster_table.append(cluster)
                         # Create actual cluster hists
                         cs_size = np.bincount(cluster['size'], minlength=30)[:30]
-                        cs_tot = np.bincount(cluster['tot'], minlength=256)[:256]
+                        cs_tot = np.bincount(cluster['tot'], minlength=512)[:512]
                         sel = np.logical_and(cluster['cluster_shape'] > 0, cluster['cluster_shape'] < 300)
                         cs_shape = np.bincount(cluster['cluster_shape'][sel], minlength=300)[:300]
                         # Add to total hists
@@ -360,13 +388,13 @@ def analyze_data(self):
                     pbar.update(upd)
                 pbar.close()
 
-                hist_occ, hist_tot, hist_tdc = interpreter.get_histograms()
+                hist_occ, hist_tot, hist_tdc, hist_trigger_dist = interpreter.get_histograms()
 
-        self._create_additional_hit_data(hist_occ, hist_tot)
+        self._create_additional_hit_data(hist_occ, hist_tot, hist_tdc, hist_trigger_dist)
         if self.cluster_hits:
             self._create_additional_cluster_data(hist_cs_size, hist_cs_tot, hist_cs_shape)
 
-    def _create_additional_hit_data(self, hist_occ, hist_tot):
+    def _create_additional_hit_data(self, hist_occ, hist_tot, hist_tdc, hist_trigger_dist):
         with tb.open_file(self.analyzed_data_file, 'r+') as out_file:
             scan_id = self.run_config['scan_id']
 
@@ -385,7 +413,23 @@ def _create_additional_hit_data(self, hist_occ, hist_tot):
                                                       complevel=5,
                                                       fletcher32=False))
 
-            # if self.analyze_tdc:  # Only store if TDC analysis is used.
+            if self.analyze_tdc:  # Only store if TDC analysis is used.
+                out_file.create_carray(out_file.root,
+                                       name='HistTdc',
+                                       title='TDC Histogram',
+                                       obj=hist_tdc,
+                                       filters=tb.Filters(complib='blosc',
+                                                          complevel=5,
+                                                          fletcher32=False))
+
+                out_file.create_carray(out_file.root,
+                                       name='HistTriggerDist',
+                                       title='Trigger Dist Histogram',
+                                       obj=hist_trigger_dist,
+                                       filters=tb.Filters(complib='blosc',
+                                                          complevel=5,
+                                                          fletcher32=False))
+
             #     out_file.create_carray(out_file.root,
             #                            name='HistTdcStatus',
             #                            title='Tdc status Histogram',
@@ -394,11 +438,11 @@ def _create_additional_hit_data(self, hist_occ, hist_tot):
             #                                               complevel=5,
             #                                               fletcher32=False))
 
-            if scan_id in ['threshold_scan']:
+            if scan_id in ['threshold_scan', 'calibrate_tot']:
                 n_injections = self.scan_config['n_injections']
                 hist_scurve = hist_occ.reshape((self.rows * self.columns, -1))
 
-                if scan_id in ['threshold_scan']:
+                if scan_id in ['threshold_scan', 'calibrate_tot']:
                     scan_params = [self.scan_config['VCAL_HIGH'] - v for v in range(self.scan_config['VCAL_LOW_start'],
                                                                                     self.scan_config['VCAL_LOW_stop'], self.scan_config['VCAL_LOW_step'])]
                     self.threshold_map, self.noise_map, self.chi2_map = au.fit_scurves_multithread(hist_scurve, scan_params, n_injections, optimize_fit_range=False)

tjmonopix2/analysis/analysis_utils.py

@@ -35,7 +35,7 @@
     ("frame", "<u1"),
     ("column", "<u2"),
     ("row", "<u2"),
-    ("charge", "<u1"),
+    ("charge", "<u2"),
     ("timestamp", "<i8"),
 ])
 
@@ -78,6 +78,15 @@ def _type_cast(self, key, val):
             return key, val
 
 
+def _tot_response_func(x, a, b, d):
+    return (a / x + 1 / b) * (x - d)
+
+
+@numba.njit
+def _inv_tot_response_func(tot, a, b, d):
+    return (np.sqrt(b**2 * (a - tot)**2 + 2 * b * d * (a + tot) + d**2) - b * a + b * tot + d) * 0.5
+
+
 def scurve(x, A, mu, sigma):
     return 0.5 * A * erf((x - mu) / (np.sqrt(2) * sigma)) + 0.5 * A
 
@@ -424,3 +433,54 @@ def fit_scurves_multithread(scurves, scan_params, n_injections=None, invert_x=Fa
     sig2D = np.reshape(sig, (512, 512))
     chi2ndf2D = np.reshape(chi2ndf, (512, 512))
     return thr2D, sig2D, chi2ndf2D
+
+
+def fit_tot_response_multithread(tot_avg, scan_params):
+
+    scurves_masked = np.ma.masked_array(tot_avg)
+
+    logger.info("Start injection ToT calibration fit on %d CPU core(s)", mp.cpu_count())
+    partialfit_tot_inj_func = partial(_fit_tot_response, scan_params=scan_params)
+
+    result_list = imap_bar(partialfit_tot_inj_func, scurves_masked.tolist(), unit=' Fits', unit_scale=True)  # Masked array entries to list leads to NaNs
+    result_array = np.array(result_list)
+    logger.info("Fit finished")
+    return np.reshape(result_array, (512, 512, 4))
+
+
+def _fit_tot_response(data, scan_params):
+    '''
+        Fit one pixel data with injection Tot calibration function.
+        Has to be global function for the multiprocessing module.
+
+        Returns:
+            (m, b, c, d, chi2/ndf)
+    '''
+
+    # Typecast to working types
+    data = np.array(data, dtype=float)
+    # Scipy bug: fit does not work on float32 values, without any error message
+    scan_params = np.array(scan_params, dtype=float)
+
+    # Deselect masked values (== nan)
+    x = scan_params[~np.isnan(data)]
+    y = data[~np.isnan(data)]
+    yerr = np.ones(len(y)) / 2  # Assume +/- 0.5 because of integer values of ToT
+
+    # Only fit data that is fittable
+    if np.all(np.isnan(y)) or x.shape[0] < 3 or len(x[y > 0]) == 0:
+        return (0., 0., 0., 0.)
+
+    p0 = [40, 0.005, 0.1]
+
+    try:
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore", OptimizeWarning)
+            popt = curve_fit(f=lambda x, a, b, d: _tot_response_func(x, a, b, d),
+                             xdata=x[y > 0], ydata=y[y > 0], p0=p0, sigma=yerr[y > 0],
+                             absolute_sigma=True)[0]
+            chi2 = np.sum((y - _tot_response_func(x, *popt)) ** 2)
+    except RuntimeError:  # fit failed
+        return (0., 0., 0., 0.)
+
+    return (*popt, chi2 / (y.shape[0] - 3 - 1))