diff --git a/ana/README.md b/ana/README.md
index 1051dd9ec..dbba48beb 100644
--- a/ana/README.md
+++ b/ana/README.md
@@ -1,12 +1,12 @@
 # pflib Analysis
 These analysis scripts are meant to be _simple_ and _fast_.
 
-They use `pandas` and `matplotlib`.
+They use `pandas`, `matplotlib` and `astropy`.
 It is **not recommended** to do these analyses on the ZCU itself
 due to disk space constraints, but you've been warned.
 
 ```
 python3 -m venv venv
 . venv/bin/activate
-pip install pandas matplotlib
+pip install pandas matplotlib astropy
 ```
diff --git a/ana/charge/parameter-time-scan.py b/ana/charge/parameter-time-scan.py
index b8fbfdfe7..2112cdea7 100644
--- a/ana/charge/parameter-time-scan.py
+++ b/ana/charge/parameter-time-scan.py
@@ -10,6 +10,7 @@
 import update_path
 from plt_gen import plt_gen
 from get_params import get_params
+from astropy.stats import sigma_clipped_stats
 
 import os
 
@@ -17,7 +18,11 @@
 import argparse
 
 import matplotlib.pyplot as plt
+import matplotlib.ticker as ticker
+
+from scipy.optimize import curve_fit
 import numpy as np
+import pandas as pd
 
 from read import read_pflib_csv
 
@@ -30,7 +35,7 @@
 parser.add_argument('-od', '--output_directory', type=Path, help='directory to which to print.')
 parser.add_argument('-r', '--repository', type=bool, help='True if you would like to indicate that the input is a repository with csv files rather than a single csv file.')
 plot_types = ['SCATTER', 'HEATMAP']
-plot_funcs = ['ADC-TIME', 'TOT-TIME', 'TOT', 'TOT-EFF', 'PARAMS', 'MULTI-CHANNEL']
+plot_funcs = ['ADC-TIME', 'ADC-ALL-CHANNELS', 'TOT-TIME', 'TOT', 'TOT-EFF', 'PARAMS', 'MULTI-CHANNEL']
 parser.add_argument('-pt','--plot_type', choices=plot_types, default=plot_types[0], type=str, help=f'Plotting type. Options: {", ".join(plot_types)}')
 parser.add_argument('-pf','--plot_function', choices=plot_funcs, default=plot_types[0], type=str, help=f'Plotting function. Options: {", ".join(plot_types)}')
 parser.add_argument('-xl','--xlabel', default='time [ns]', type=str, help=f'What to label the x-axis with.')
@@ -71,10 +76,197 @@ def set_xticks (
     xmin = 25*np.floor(xmin/25)
     xmax = 25*np.ceil(xmax/25)
     ax.set_xticks(ticks = np.arange(xmin, xmax+1, 25), minor=False)
-    ax.set_xticks(ticks = np.arange(xmin, xmax, 25/16), minor=True) 
+    ax.set_xticks(ticks = np.arange(xmin, xmax, 25/16), minor=True)
+
+#################### HELPER FUNCTIONS ########################
+
+def linear(x, k, m):
+    return k*x + m
 
 #################### PLOTTING FUNCTIONS ########################
 
+def adc_all_channels(
+        samples, 
+        run_params, 
+        ax_holder,
+        plot_params="all",   # "all" or "one"
+        param_index=6        # used only when plot_params="one"
+):
+    """
+    Plot max ADC vs channels, compute RMS trends, and fit dependence on parameters.
+    """
+
+    def compute_rms(ch_df, pedestal):
+        """RMS for a single channel dataframe."""
+        return np.sqrt(((ch_df['adc'] - pedestal)**2).mean())
+
+    def compute_all_rms(df, pedestal):
+        """Return {channel: rms}"""
+        rms = df.groupby('channel').apply(lambda c: compute_rms(c, pedestal))
+        return rms.reset_index(name='rms')
+
+    def conditional_legend(ax, max_labels=37, **kwargs):
+        """Only draw legend if below threshold."""
+        handles, labels = ax.get_legend_handles_labels()
+        if len(labels) < max_labels:
+            ax.legend(**kwargs)
+
+    def save_and_close(fig, fname):
+        fig.savefig(fname, dpi=400)
+        plt.close(fig)
+
+    charges_group = samples.groupby("nr channels")
+    fig_rms, ax_rms = plt.subplots(1, 1)
+    ax_rms.set_xlabel('Channels activated')
+    ax_rms.set_ylabel('Average RMS of non-activated channels')
+
+    link = 0
+    central_val = 54 if link == 1 else 18
+
+    activated_channels_list = []
+    avg_rms_list = []
+
+    runs = len(charges_group)
+    activated_channels_list = [[] for _ in range(runs - 1)]
+    avg_rms_list = [[] for _ in range(runs - 1)]
+
+    parameter_values = set()
+
+    # Loop over charge groups
+    for i, (charge_id, charge_df) in enumerate(charges_group):
+        print(f'Running {i+1} out of {runs}')
+        # Pedestal case
+        if i == 0:
+            fig, ax = plt.subplots()
+            df = charge_df[charge_df['nr channels'] == 0]
+
+            df = df[df['channel'] < 36] if link == 0 else df[df['channel'] >= 36]
+
+            mean, med, std = sigma_clipped_stats(df['adc'], sigma=3)
+            pedestal = mean
+
+            ax.scatter(df['channel'], df['adc'], s=5, color='r', lw=1)
+            ax.set_ylabel('ADC')
+            ax.set_xlabel('Channel')
+
+            save_and_close(fig, 'Pedestal.png')
+            continue
+
+        fig, ax1 = plt.subplots()
+        fig_time, ax_time = plt.subplots()
+
+        ax1.set_ylabel(r'$\Delta$ADC')
+        ax1.set_xlabel('Channel')
+        ax1.xaxis.set_minor_locator(ticker.AutoMinorLocator())
+
+        ax_time.set_ylabel('ADC')
+        ax_time.set_xlabel('Time')
+
+        # Parameter grouping
+        try:
+            param_group, param_name = get_params(charge_df, 0)
+        except Exception:
+            param_group = [(None, charge_df)]
+
+        cmap = plt.get_cmap('viridis')
+
+        for j, (param_id, param_df) in enumerate(param_group):
+
+            # --- Parameter selection ---
+            if plot_params == "one":
+                if j != param_index:
+                    continue
+
+            try:
+                key = param_name.split('.')[1]
+                val = param_df[param_name].iloc[0]
+            except Exception:  
+                key = "Channels flashed, voltage [mVpp]"
+                val = (4, 2000)
+
+            parameter_values.add(val)
+
+            max_diff = (
+                param_df
+                .groupby('channel')['adc']
+                .apply(lambda s: (s - pedestal).abs().max())
+                .reset_index(name='max_diff')
+            )
+            
+            # Get RMS
+            rms_df = compute_all_rms(param_df, pedestal)
+            rms_vals = rms_df['rms'].tolist()
+            avg_rms = 0
+            count = 0
+            for k in range(central_val - 18, central_val + 18):
+                if central_val - (i - 1) <= k <= central_val + (i - 1):
+                    continue
+                avg_rms += rms_vals[k]
+                count += 1
+
+            if count > 0:
+                avg_rms /= count
+                avg_rms_list[i - 1].append(avg_rms)
+                activated_channels_list[i - 1].append(i)
+
+            merged = max_diff.merge(rms_df, on='channel')
+            ax1.set_ylim(-10, merged['max_diff'].max() + 10)
+            ax1.scatter(
+                merged['channel'], merged['max_diff'],
+                label=f'{key} = {val}',
+                s=5, color=cmap(j / len(param_group)), lw=1
+            )
+            if link == 0:
+                link_df = param_df[param_df['channel'] < 36]
+            elif link == 1:
+                link_df = param_df[param_df['channel'] >= 36]
+            for k, (ch_id, ch_df) in enumerate(link_df.groupby('channel')):
+                ax_time.scatter(
+                    ch_df['time'], ch_df['adc'],
+                    label=f'ch{ch_id}',
+                    s=5, color=plt.get_cmap('tab20')(k / 20)
+                    )
+
+        # Pedestal and link lines
+        ax1.axhline(y=0, color='k', linestyle='--', linewidth=0.8)
+        ax1.axvline(x=35.5, color='r', linestyle='--', alpha=0.5)
+
+        conditional_legend(ax1, fontsize=8)
+        conditional_legend(ax_time, fontsize=4, ncols=3)
+
+        save_and_close(fig, f'adc_channels_{i}.png')
+        save_and_close(fig_time, f'adc_time_{i}.png')
+
+    # slope/intercept plots
+    activated_T = list(map(list, zip(*activated_channels_list)))
+    rms_T = list(map(list, zip(*avg_rms_list)))
+    parameter_values = sorted(parameter_values)
+
+    if len(activated_T) > 1:
+        slopes = []
+        intercepts = []
+
+        for xvals, yvals in zip(activated_T, rms_T):
+            popt, _ = curve_fit(linear, xvals, yvals)
+            slopes.append(popt[0])
+            intercepts.append(popt[1])
+
+            ax_rms.scatter(xvals, yvals, s=8)
+            xfit = np.linspace(min(xvals), max(xvals), 100)
+            ax_rms.plot(xfit, linear(xfit, *popt),
+                        label=f'Fit: y={popt[0]:.3f}x + {popt[1]:.3f}, CALIB={parameter_values[len(slopes)-1]}')
+
+        conditional_legend(ax_rms, fontsize=8)
+        save_and_close(fig_rms, 'avg_rms.png')
+
+        fig_par, ax_par = plt.subplots()
+        ax_par.scatter(parameter_values, slopes, label='slopes')
+        ax_par.scatter(parameter_values, intercepts, label='intercepts')
+        ax_par.legend()
+        ax_par.set_xlabel('CALIB')
+        ax_par.set_title('Slope and intercepts from fits')
+        save_and_close(fig_par, 'slope_intercept.png')
+
 def time(
     samples,
     run_params,
@@ -374,6 +566,10 @@ def multiparams(
 
 ############################## MAIN ###################################
 
+if args.plot_function == 'ADC-ALL-CHANNELS':
+    plt_gen(adc_all_channels, samples, run_params, args.output,
+            xlabel = args.xlabel, ylabel = args.ylabel)
+
 if args.plot_function == 'ADC-TIME' and args.plot_type == 'SCATTER':
     plt_gen(partial(time, yval = 'adc'), samples, run_params, args.output,
             xlabel = args.xlabel, ylabel = args.ylabel)
diff --git a/app/tool/tasks/multi_channel_scan.cxx b/app/tool/tasks/multi_channel_scan.cxx
new file mode 100644
index 000000000..01188536c
--- /dev/null
+++ b/app/tool/tasks/multi_channel_scan.cxx
@@ -0,0 +1,261 @@
+#include "multi_channel_scan.h"
+
+#include <filesystem>
+#include <nlohmann/json.hpp>
+
+#include "../daq_run.h"
+#include "../econ_links.h"
+#include "load_parameter_points.h"
+#include "pflib/utility/string_format.h"
+
+ENABLE_LOGGING();
+
+template <class EventPacket>
+static void multi_channel_scan_writer(Target* tgt, pflib::ROC& roc,
+                                      size_t n_events, int calib, bool isLED,
+                                      int highrange, int link,
+                                      std::string fname, int start_bx,
+                                      int n_bx) {
+  int ch0{0};
+  link == 0 ? ch0 = 18 : ch0 = 54;
+  int n_links = determine_n_links(tgt);
+
+  if (isLED) {
+    auto refvol_page =
+        pflib::utility::string_format("REFERENCEVOLTAGE_%d", link);
+    auto test_param_handle = roc.testParameters();
+    test_param_handle.add(refvol_page, "INTCTEST", 1);
+    auto test_param = test_param_handle.apply();
+
+    int phase_strobe{0};
+    int charge_to_l1a{0};
+    double time{0};
+    double clock_cycle{25.0};
+    int n_phase_strobe{16};
+    int offset{1};
+    int nr_channels{-1};
+    DecodeAndWriteToCSV<EventPacket> writer{
+        fname,
+        [&](std::ofstream& f) {
+          nlohmann::json header;
+          header["highrange"] = highrange;
+          header["preCC"] = !highrange;
+          header["ledflash"] = isLED;
+          f << std::boolalpha << "# " << header << '\n'
+            << "nr channels," << "time,";
+          f << "channel,";
+          f << pflib::packing::Sample::to_csv_header << '\n';
+        },
+        [&](std::ofstream& f, const EventPacket& ep) {
+          for (int ch{0}; ch < 72; ch++) {
+            f << nr_channels + 1 << ',';
+            f << time << ',';
+            f << ch << ',';
+            if constexpr (std::is_same_v<
+                              EventPacket,
+                              pflib::packing::MultiSampleECONDEventPacket>) {
+              ep.samples[ep.i_soi].channel(link, ch).to_csv(f);
+            } else if constexpr (std::is_same_v<
+                                     EventPacket,
+                                     pflib::packing::SingleROCEventPacket>) {
+              ep.channel(ch).to_csv(f);
+            } else {
+              PFEXCEPTION_RAISE("BadConf",
+                                "Unable to do all_channels_to_csv for the "
+                                "currently configured format.");
+            }
+            f << '\n';
+          }
+        },
+        n_links  // number of links
+    };
+    tgt->setup_run(1 /* dummy - not stored */, pftool::state.daq_format_mode,
+                   1 /* dummy */);
+    // Do the scan for increasing amount of channels
+    for (nr_channels; nr_channels < 1; nr_channels++) {
+      auto test_param_handle = roc.testParameters();
+      if (nr_channels == -1) {  // In case of -1, just take pedestal data
+        daq_run(tgt, "PEDESTAL", writer, 1, pftool::state.daq_rate);
+        continue;
+      }
+      int central_charge_to_l1a = tgt->fc().fc_get_setup_led();
+      for (charge_to_l1a = central_charge_to_l1a + start_bx;
+           charge_to_l1a < central_charge_to_l1a + start_bx + n_bx;
+           charge_to_l1a++) {
+        tgt->fc().fc_setup_led(charge_to_l1a);
+        pflib_log(info) << "led_to_l1a = " << tgt->fc().fc_get_setup_led();
+        for (phase_strobe = 0; phase_strobe < n_phase_strobe; phase_strobe++) {
+          auto phase_strobe_test_handle =
+              roc.testParameters()
+                  .add("TOP", "PHASE_STROBE", phase_strobe)
+                  .apply();
+          pflib_log(info) << "TOP.PHASE_STROBE = " << phase_strobe;
+          usleep(10);  // make sure parameters are applied
+          time =
+              (charge_to_l1a - central_charge_to_l1a + offset) * clock_cycle -
+              phase_strobe * clock_cycle / n_phase_strobe;
+          daq_run(tgt, "LED", writer, 1, pftool::state.daq_rate);
+        }
+      }
+      // reset charge_to_l1a to central value
+      tgt->fc().fc_setup_led(central_charge_to_l1a);
+    }
+  } else {
+    auto refvol_page =
+        pflib::utility::string_format("REFERENCEVOLTAGE_%d", link);
+    auto test_param_handle = roc.testParameters();
+    test_param_handle.add(refvol_page, "CALIB", highrange ? calib : 0)
+        .add(refvol_page, "CALIB_2V5", highrange ? 0 : calib)
+        .add(refvol_page, "INTCTEST", 1)
+        .add(refvol_page, "CHOICE_CINJ", highrange ? 1 : 0);
+    auto test_param = test_param_handle.apply();
+
+    std::filesystem::path parameter_points_file =
+        pftool::readline("File of parameter points: ");
+
+    pflib_log(info) << "loading parameter points file...";
+    auto [param_names, param_values] =
+        load_parameter_points(parameter_points_file);
+    pflib_log(info) << "successfully loaded parameter points";
+
+    int phase_strobe{0};
+    int charge_to_l1a{0};
+    double time{0};
+    double clock_cycle{25.0};
+    int n_phase_strobe{16};
+    int offset{1};
+    std::size_t i_param_point{0};
+    int nr_channels{-1};
+    DecodeAndWriteToCSV<EventPacket> writer{
+        fname,
+        [&](std::ofstream& f) {
+          nlohmann::json header;
+          header["highrange"] = highrange;
+          header["precc"] = !highrange;
+          header["ledflash"] = isLED;
+          f << std::boolalpha << "# " << header << '\n'
+            << "nr channels," << "time,";
+          for (const auto& [page, parameter] : param_names) {
+            f << page << '.' << parameter << ',';
+          }
+          f << "channel,";
+          f << pflib::packing::Sample::to_csv_header << '\n';
+        },
+        [&](std::ofstream& f, const EventPacket& ep) {
+          for (int ch{0}; ch < 72; ch++) {
+            f << nr_channels + 1 << ',';
+            f << time << ',';
+
+            for (const auto& val : param_values[i_param_point]) {
+              f << val << ',';
+            }
+            f << ch << ',';
+
+            if constexpr (std::is_same_v<
+                              EventPacket,
+                              pflib::packing::MultiSampleECONDEventPacket>) {
+              ep.samples[ep.i_soi].channel(link, ch).to_csv(f);
+            } else if constexpr (std::is_same_v<
+                                     EventPacket,
+                                     pflib::packing::SingleROCEventPacket>) {
+              ep.channel(ch).to_csv(f);
+            } else {
+              PFEXCEPTION_RAISE("BadConf",
+                                "Unable to do all_channels_to_csv for the "
+                                "currently configured format.");
+            }
+            f << '\n';
+          }
+        },
+        n_links  // number of links
+    };
+    tgt->setup_run(1 /* dummy - not stored */, pftool::state.daq_format_mode,
+                   1 /* dummy */);
+    // Do the scan for increasing amount of channels
+    for (nr_channels; nr_channels < 18; nr_channels++) {
+      pflib_log(info) << "running scan " << nr_channels + 1 << " out of 19";
+      auto test_param_handle = roc.testParameters();
+      if (nr_channels == -1) {  // In case of -1, just take pedestal data
+        daq_run(tgt, "PEDESTAL", writer, n_events, pftool::state.daq_rate);
+        continue;
+      }
+      for (int ch = ch0 - nr_channels; ch <= ch0 + nr_channels; ch++) {
+        auto channel_page = pflib::utility::string_format("CH_%d", ch);
+        test_param_handle.add(channel_page, "HIGHRANGE", 1);
+      }
+      auto test_param_two = test_param_handle.apply();
+      i_param_point = 0;
+      for (; i_param_point < param_values.size(); i_param_point++) {
+        // set parameters
+        auto test_param_builder = roc.testParameters();
+        for (std::size_t i_param{0}; i_param < param_names.size(); i_param++) {
+          const auto& full_page = param_names[i_param].first;
+          std::string str_page = full_page.substr(0, full_page.find("_"));
+          const auto& param = param_names[i_param].second;
+          int value = param_values[i_param_point][i_param];
+
+          if (str_page == "CH" || str_page == "TOP") {
+            test_param_builder.add(full_page, param, value);
+          } else {
+            test_param_builder.add(str_page + "_" + std::to_string(link), param,
+                                   value);
+          }
+          pflib_log(info) << param << " = " << value;
+        }
+        auto test_param_three = test_param_builder.apply();
+        // timescan
+        int central_charge_to_l1a;
+        central_charge_to_l1a = tgt->fc().fc_get_setup_calib();
+        for (charge_to_l1a = central_charge_to_l1a + start_bx;
+             charge_to_l1a < central_charge_to_l1a + start_bx + n_bx;
+             charge_to_l1a++) {
+          tgt->fc().fc_setup_calib(charge_to_l1a);
+          pflib_log(info) << "charge_to_l1a = "
+                          << tgt->fc().fc_get_setup_calib();
+          for (phase_strobe = 0; phase_strobe < n_phase_strobe;
+               phase_strobe++) {
+            auto phase_strobe_test_handle =
+                roc.testParameters()
+                    .add("TOP", "PHASE_STROBE", phase_strobe)
+                    .apply();
+            pflib_log(info) << "TOP.PHASE_STROBE = " << phase_strobe;
+            usleep(10);  // make sure parameters are applied
+            time =
+                (charge_to_l1a - central_charge_to_l1a + offset) * clock_cycle -
+                phase_strobe * clock_cycle / n_phase_strobe;
+            daq_run(tgt, "CHARGE", writer, n_events, pftool::state.daq_rate);
+          }
+        }
+        // reset charge_to_l1a to central value
+        tgt->fc().fc_setup_calib(central_charge_to_l1a);
+      }
+    }
+  }
+}
+
+void multi_channel_scan(Target* tgt) {
+  size_t n_events = pftool::readline_int("How many events per time point? ", 1);
+  int calib = pftool::readline_int("Setting for calib pulse amplitude? ", 256);
+  bool isLED = pftool::readline_bool("Use external LED flashes?", false);
+  bool highrange =
+      pftool::readline_bool("Use highrange (Y) or preCC (N)? ", false);
+  int start_bx = pftool::readline_int("Starting BX? ", -1);
+  int n_bx = pftool::readline_int("Number of BX? ", 3);
+  int link{0};
+  pftool::readline_bool("Link 0 [Y] or link 1 [N]", "true") ? link = 0
+                                                            : link = 1;
+  std::string fname = pftool::readline_path("multi-channel-scan", ".csv");
+
+  pflib::ROC roc{tgt->roc(pftool::state.iroc)};
+
+  if (pftool::state.daq_format_mode == Target::DaqFormat::SIMPLEROC) {
+    multi_channel_scan_writer<pflib::packing::SingleROCEventPacket>(
+        tgt, roc, n_events, calib, isLED, highrange, link, fname, start_bx,
+        n_bx);
+  } else if (pftool::state.daq_format_mode ==
+             Target::DaqFormat::ECOND_SW_HEADERS) {
+    multi_channel_scan_writer<pflib::packing::MultiSampleECONDEventPacket>(
+        tgt, roc, n_events, calib, isLED, highrange, link, fname, start_bx,
+        n_bx);
+  }
+}
diff --git a/app/tool/tasks/multi_channel_scan.h b/app/tool/tasks/multi_channel_scan.h
new file mode 100644
index 000000000..dbcd48f7d
--- /dev/null
+++ b/app/tool/tasks/multi_channel_scan.h
@@ -0,0 +1,14 @@
+#pragma once
+
+#include "../pftool.h"
+
+/*
+ * TASKS.MULTI_CHANNEL_SCAN
+ *
+ * Scan a internal calibration pulse in time by varying the charge_to_l1a
+ * and top.phase_strobe parameters, and then change given parameters pulse-wise.
+ * In essence, an implementation of gen_scan into charge_timescan, to see how
+ * pulse shapes transform with varying parameters.
+ *
+ */
+void multi_channel_scan(Target* tgt);
diff --git a/app/tool/tasks/tasks.cxx b/app/tool/tasks/tasks.cxx
index a677517dc..596ec0f7f 100644
--- a/app/tool/tasks/tasks.cxx
+++ b/app/tool/tasks/tasks.cxx
@@ -11,6 +11,7 @@
 #include "inv_vref_scan.h"
 #include "level_pedestals.h"
 #include "load_parameter_points.h"
+#include "multi_channel_scan.h"
 #include "noinv_vref_scan.h"
 #include "parameter_timescan.h"
 #include "sampling_phase_scan.h"
@@ -40,6 +41,9 @@ auto menu_tasks =
         ->line("SAMPLING_PHASE_SCAN",
                "scan phase_ck, pedestal for clock phase alignment",
                sampling_phase_scan)
+        ->line("MULTI_CHANNEL_SCAN",
+               "scans multiple channels to look for cross-talk",
+               multi_channel_scan)
         ->line("VT50_SCAN",
                "Hones in on the vt50 with a binary or bisectional scan",
                vt50_scan)