calibration.py

Go to the documentation of this file.

r"""!
@file  calibration.py
@brief `main()` shows how to
       [minimize()](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html)
       the [objective function()](#calibration.objective_function) for antenna phase center
       calibration.
"""
import ROOT
import re
import numpy as np
import polars as pl
from pathlib import Path
from initialise import load_pueoEvent_Dataset

 
 
def generate_antenna_phase_center_pairs_with_placeholders(qrh_dot_dat: Path) -> pl.DataFrame:
    r"""!Creates [antenna pairs](#antenna_pairs.generate_MI_antenna_pairs) with placeholder
         integers as antenna phase centers.
@ingroup CALPC
@param[in] qrh_dot_dat              Path to `qrh.dat`
@retval nominal_phase_center_pairs
 
*   Output schema:
```
$ A1_AntNum            <enum>
$ A1_Ring                <u8>
$ A1_Rho (placeholder)  <u16>
$ A1_Phi (placeholder)  <u16>
$ A1_Z (placeholder)    <u16>
$ A1_Rho                <f64>
$ A1_Phi                <f64>
$ A1_Z                  <f64>
$ A2_AntNum            <enum>
$ A2_Ring                <u8>
$ A2_Rho (placeholder)  <u16>
$ A2_Phi (placeholder)  <u16>
$ A2_Z (placeholder)    <u16>
$ A2_Rho                <f64>
$ A2_Phi                <f64>
$ A2_Z                  <f64>
```
*   Details:
    -#  Calls #antenna_attributes.read_MI_antenna_geometry() to read `qrh.dat` to obtain
        antenna face center positions
    -#  Calls #antenna_attributes.get_MI_nominal_phase_center() to convert face centers to
        nominal phase centers
    -#  Adds three columns to store placeholder values of the phase centers; these integer placeholders
        serve as dictionary keys.
```
shape: (96, 9)
┌────────┬────────┬──────┬────────────────┬────────────────┬───────────────┬──────────┬───────────┬───────┐
│ AntNum ┆ AntIdx ┆ Ring ┆ Rho            ┆ Phi            ┆ Z             ┆ Rho      ┆ Phi       ┆ Z     │
│ ---    ┆ ---    ┆ ---  ┆ (placeholder)  ┆ (placeholder)  ┆ (placeholder) ┆ ---      ┆ ---       ┆ ---   │
│ enum   ┆ u32    ┆ u8   ┆ ---            ┆ ---            ┆ ---           ┆ f64      ┆ f64       ┆ f64   │
│        ┆        ┆      ┆ u16            ┆ u16            ┆ u16           ┆          ┆           ┆       │
╞════════╪════════╪══════╪════════════════╪════════════════╪═══════════════╪══════════╪═══════════╪═══════╡
│ 101    ┆ 0      ┆ 1    ┆ 0              ┆ 96             ┆ 192           ┆ 1.188    ┆ 0.0       ┆ 5.114 │
│ 201    ┆ 1      ┆ 2    ┆ 1              ┆ 97             ┆ 193           ┆ 2.296    ┆ 0.0       ┆ 1.457 │
│ 301    ┆ 2      ┆ 3    ┆ 2              ┆ 98             ┆ 194           ┆ 2.296    ┆ 0.0       ┆ 0.728 │
│ 401    ┆ 3      ┆ 4    ┆ 3              ┆ 99             ┆ 195           ┆ 2.296    ┆ 0.0       ┆ 0.0   │
│ 102    ┆ 4      ┆ 1    ┆ 4              ┆ 100            ┆ 196           ┆ 1.076821 ┆ 0.261538  ┆ 4.476 │
│ …      ┆ …      ┆ …    ┆ …              ┆ …              ┆ …             ┆ …        ┆ …         ┆ …     │
│ 423    ┆ 91     ┆ 4    ┆ 91             ┆ 187            ┆ 283           ┆ 2.295825 ┆ -0.523638 ┆ 0.0   │
│ 124    ┆ 92     ┆ 1    ┆ 92             ┆ 188            ┆ 284           ┆ 1.07708  ┆ -0.26224  ┆ 4.476 │
│ 224    ┆ 93     ┆ 2    ┆ 93             ┆ 189            ┆ 285           ┆ 2.295502 ┆ -0.261893 ┆ 1.457 │
│ 324    ┆ 94     ┆ 3    ┆ 94             ┆ 190            ┆ 286           ┆ 2.295502 ┆ -0.261893 ┆ 0.728 │
│ 424    ┆ 95     ┆ 4    ┆ 95             ┆ 191            ┆ 287           ┆ 2.295502 ┆ -0.261893 ┆ 0.0   │
└────────┴────────┴──────┴────────────────┴────────────────┴───────────────┴──────────┴───────────┴───────┘
```
    -#  Calls #antenna_pairs.generate_MI_antenna_pairs with the above dataframe as the input.
    """
    from antenna_pairs import generate_MI_antenna_pairs
    from antenna_attributes import read_MI_antenna_geometry, get_MI_nominal_phase_center
 
    fc: pl.DataFrame = read_MI_antenna_geometry(qrh_dot_dat=qrh_dot_dat)
 
    nominal_and_placeholders: pl.DataFrame = (
        get_MI_nominal_phase_center(face_centers=fc, coordinates="cylindrical")
        .with_columns(
            pl.arange(pl.len()).alias("Rho (placeholder)").cast(pl.UInt16),
            pl.arange(pl.len(), pl.len() * 2).alias("Phi (placeholder)").cast(pl.UInt16),
            pl.arange(pl.len() * 2, pl.len() * 3).alias("Z (placeholder)").cast(pl.UInt16)
        )
        .select(
            "AntNum", "AntIdx", "Ring", "Rho (placeholder)", "Phi (placeholder)", "Z (placeholder)",
            pl.col("Rho[m]").alias("Rho"), pl.col("Phi[rad]").alias("Phi"), pl.col("Z[m]").alias("Z")
        )
    )
 
    nominal_phase_center_pairs: pl.DataFrame = (
        generate_MI_antenna_pairs(antennas=nominal_and_placeholders)
        .select(pl.all().exclude(r"^A[12]_AntIdx$"))
    )
 
    return nominal_phase_center_pairs
 
 
def __get_all_run_numbers(pueo_mc_data: Path) -> [int]:
    """! A temporary helper function of #load_dataset_and_measure_time_delays()
    @ingroup CALPC
    @param[in] pueo_mc_data  See the parameter of #initialise.load_pueoEvent_Dataset
 
    *   This function retrieves the run numbers (negative sign allowed) of the `run/` subdirectories
        in `pueo_mc_data`
 
    *   The runs are returned as a list of integers.
 
    *   The function is here only for testing purposes,
        as it makes little sense to use every single run for antenna calibration.
    """
 
    # regular expression, allowing positive and negative signs in run numbers
    run_number_pattern_ = re.compile(r'^run([+\-]?\d+)$')
 
    return [int(run_number_pattern_.match(run_directories_.name).group(1))
            for run_directories_ in pueo_mc_data.rglob("run*")]
 
 
def load_dataset_and_measure_time_delays(dataset: ROOT.pueo.Dataset,
                                         phase_center_pairs: pl.DataFrame) -> pl.DataFrame:
    r"""! Uses [SciPy's correlate()](https://docs.scipy.org/doc//scipy-1.16.2/reference/generated/scipy.signal.correlate.html)
    to "measure" the delays of [waveforms](#waveform_plots.load_waveforms) between antenna pairs
    @ingroup CALPC
    @param[in] dataset              The output of #initialise.load_pueoEvent_Dataset
    @param[in] phase_center_pairs   The output of #generate_antenna_phase_center_pairs_with_placeholders
    @retval    pulser_directions_and_measured_time_delays
 
    *   The input `dataset` should have been initialized with a (valid) run number,
        but that is not the only run that gets loaded.
        For now, [all runs](#__get_all_run_numbers) are used.
    *   Output schema:
```
$ measured time delays (ns)                      <f64>
$ max correlation                                <f64>
$ A1_AntNum                                     <enum>
$ A1_Ring                                         <u8>
$ A1_Rho (placeholder)                           <u16>
$ A1_Phi (placeholder)                           <u16>
$ A1_Z (placeholder)                             <u16>
$ A1_Rho                                         <f64>
$ A1_Phi                                         <f64>
$ A1_Z                                           <f64>
$ A2_AntNum                                     <enum>
$ A2_Ring                                         <u8>
$ A2_Rho (placeholder)                           <u16>
$ A2_Phi (placeholder)                           <u16>
$ A2_Z (placeholder)                             <u16>
$ A2_Rho                                         <f64>
$ A2_Phi                                         <f64>
$ A2_Z                                           <f64>
$ incoming pulse direction (Phi Theta) <array[f64, 2]>
```
    *   The `dataset` is also used to load the directions of the incoming pulses
        (\f$^\circ\phi\f$ and \f$^\circ\theta\f$ in payload coordinates).
    *   `measured time delays (ns)` is computed by the helper function
        #_scipy_cross_correlate_all_pairs_for_one_event().
    """
 
    # Note: eventually will have to get pulser direction somehow, instead of from the truth file
    from locate_signal import _get_true_direction
    from waveform_plots import load_waveforms, upsample_waveforms
 
    # initialize an empty dataframe that will be in-place updated
    pulser_directions_and_measured_time_delays = pl.DataFrame()
 
    # Note: probably doesn't make sense to use all runs, for now it's okay
    all_runs: [int] = __get_all_run_numbers(
        pueo_mc_data=Path(
            dataset.getDataDir(ROOT.pueo.Dataset.PUEO_MC_DATA)
        )
    )
    for run in all_runs:
        dataset.loadRun(run, ROOT.pueo.Dataset.PUEO_MC_DATA)
        for i in range(dataset.N()):  # loop through all events in the run
            dataset.getEntry(i)
 
            wf: pl.DataFrame = (
                load_waveforms(dataset=dataset)
                # .filter(pl.col("Pol") == "V")
                .select("AntNum", "Pol", "waveforms (volts)", "step size (nanoseconds)")
            )
            up: pl.DataFrame = upsample_waveforms(waveforms=wf, upsample_factor=5)
 
            pulser_directions_and_measured_time_delays.vstack(
                _scipy_cross_correlate_all_pairs_for_one_event(
                    phase_center_pairs=phase_center_pairs,
                    waveforms=up
                ).with_columns(
                    (pl.lit(_get_true_direction(dataset=dataset)) * np.pi / 180)
                    .cast(pl.Array(pl.Float64, 2))
                    .alias("incoming pulse direction (Phi Theta)")
                ),
                in_place=True
            )
    return pulser_directions_and_measured_time_delays
 
 
def _scipy_cross_correlate_all_pairs_for_one_event(phase_center_pairs: pl.DataFrame,
                                                   waveforms: pl.DataFrame) -> pl.DataFrame:
    r"""! A helper function of #load_dataset_and_measure_time_delays()
    @ingroup CALPC
    @param[in] phase_center_pairs   The output of #generate_antenna_phase_center_pairs_with_placeholders
    @param[in] waveforms            The output of #waveform_plots.load_waveforms
    @retval    measured_time_delays
 
    *   Output schema:
```
$ measured time delays (ns)  <f64>
$ max correlation            <f64>
$ A1_AntNum                 <enum>
$ A1_Ring                     <u8>
$ A1_Rho (placeholder)       <u16>
$ A1_Phi (placeholder)       <u16>
$ A1_Z (placeholder)         <u16>
$ A1_Rho                     <f64>
$ A1_Phi                     <f64>
$ A1_Z                       <f64>
$ A2_AntNum                 <enum>
$ A2_Ring                     <u8>
$ A2_Rho (placeholder)       <u16>
$ A2_Phi (placeholder)       <u16>
$ A2_Z (placeholder)         <u16>
$ A2_Rho                     <f64>
$ A2_Phi                     <f64>
$ A2_Z                       <f64>
```
    *   Details:
        -#  [Load the waveforms](#waveform_plots.load_waveforms) of both antennas in the antenna pair
        -#  Cross correlate the waveform pair
        -#  Figure out by how much one needs to shift the second waveform against the first
            in order to achieve the maximum correlation.
        -#  The above is called `measured time delays (ns)` in the output dataframe.
        -#  For more details, see [Explanation](#scipy_corr_expl).
    """
    from scipy.signal import correlate, correlation_lags
 
    # the fact that the waveforms are in the same column guarantees they have the same length,
    # which is probably an integer multiple of 1024 (if upsampled)
    signal_length = len(waveforms["waveforms (volts)"][0])
    CORRLAGS: np.ndarray = correlation_lags(signal_length, signal_length)
 
    STEP_SIZE = waveforms["step size (nanoseconds)"][0]
    # For peace of mind, check that all signals were sampled at the same step size :)
    assert np.mean(waveforms["step size (nanoseconds)"].to_numpy()) - STEP_SIZE < 1e-10
 
    measured_time_delays = (
        phase_center_pairs
        .join(waveforms, left_on="A1_AntNum", right_on="AntNum")
        .join(waveforms, left_on=["A2_AntNum", "Pol"], right_on=["AntNum", "Pol"])
        .with_columns(
            pl.struct("waveforms (volts)", "waveforms (volts)_right")
            .map_batches(
                lambda s: np.array([
                    correlate(  # zero-center normalized cross-correlation (ZNCC), bounded b/w [-1, 1]
                        (wf1 - np.mean(wf1)) / np.std(wf1),
                        (wf2 - np.mean(wf2)) / np.std(wf2),
                        method='fft'
                    ) / len(wf1)
                    for wf2, wf1 in s.to_numpy()
                ])
            ).alias("correlation"),
        )
        .with_columns(
            pl.col("correlation").arr.max().alias("max correlation"),
            pl.col("correlation").arr.arg_max().alias("max correlation idx"),
        )
        .with_columns(  # Find the time where the maximum correlation occurs and convert to [ns]
            pl.col("max correlation idx").map_batches(lambda idx: CORRLAGS[idx] * STEP_SIZE)
            .alias("measured time delays (ns)")
        )
        .select(
            "measured time delays (ns)", "max correlation",
            r"^A[12]_(?:AntNum|Ring|Rho|Phi|Z)(?: \‍(placeholder\‍))?$",
        )
 
    )
    return measured_time_delays
 
 
def retrieve_relevant_antennas_from_all_pairs(phase_center_pairs: pl.DataFrame):
    r"""! Retrieves the (unique) antennas, given a bunch of antenna pairs.
    @ingroup CALPC
    @param[in] phase_center_pairs   The output of #load_dataset_and_measure_time_delays
    @retval    relevant_antennas
 
    *   This step is not trivial, as explained below.
    *   It is conceivable that not all antennas receive signals clear enough for calibration purposes
        (resulting in weak correlations).
    *   As such, the weakly correlated antenna pairs should be filtered out.
        For example,
```python3
    pulser_directions_and_pairwise_time_delays = (
        load_dataset_and_measure_time_delays(dataset=ds, phase_center_pairs=nominal_pairs)
        .filter(pl.col("max correlation") > 0.8)
    )
```
    *   It is possible that some antennas end up getting filtered out.
    *   Sample output (note that only 76 out of 96 MI antennas survived the filter in this case):
```
shape: (76, 8)
┌────────┬──────┬───────────────────┬───────────────────┬─────────────────┬──────────┬───────────┬───────┐
│ AntNum ┆ Ring ┆ Rho (placeholder) ┆ Phi (placeholder) ┆ Z (placeholder) ┆ Rho      ┆ Phi       ┆ Z     │
│ ---    ┆ ---  ┆ ---               ┆ ---               ┆ ---             ┆ ---      ┆ ---       ┆ ---   │
│ enum   ┆ u8   ┆ u16               ┆ u16               ┆ u16             ┆ f64      ┆ f64       ┆ f64   │
╞════════╪══════╪═══════════════════╪═══════════════════╪═════════════════╪══════════╪═══════════╪═══════╡
│ 101    ┆ 1    ┆ 0                 ┆ 96                ┆ 192             ┆ 1.188    ┆ 0.0       ┆ 5.114 │
│ 201    ┆ 2    ┆ 1                 ┆ 97                ┆ 193             ┆ 2.296    ┆ 0.0       ┆ 1.457 │
│ 301    ┆ 3    ┆ 2                 ┆ 98                ┆ 194             ┆ 2.296    ┆ 0.0       ┆ 0.728 │
│ 401    ┆ 4    ┆ 3                 ┆ 99                ┆ 195             ┆ 2.296    ┆ 0.0       ┆ 0.0   │
│ 102    ┆ 1    ┆ 4                 ┆ 100               ┆ 196             ┆ 1.076821 ┆ 0.261538  ┆ 4.476 │
│ …      ┆ …    ┆ …                 ┆ …                 ┆ …               ┆ …        ┆ …         ┆ …     │
│ 423    ┆ 4    ┆ 91                ┆ 187               ┆ 283             ┆ 2.295825 ┆ -0.523638 ┆ 0.0   │
│ 124    ┆ 1    ┆ 92                ┆ 188               ┆ 284             ┆ 1.07708  ┆ -0.26224  ┆ 4.476 │
│ 224    ┆ 2    ┆ 93                ┆ 189               ┆ 285             ┆ 2.295502 ┆ -0.261893 ┆ 1.457 │
│ 324    ┆ 3    ┆ 94                ┆ 190               ┆ 286             ┆ 2.295502 ┆ -0.261893 ┆ 0.728 │
│ 424    ┆ 4    ┆ 95                ┆ 191               ┆ 287             ┆ 2.295502 ┆ -0.261893 ┆ 0.0   │
└────────┴──────┴───────────────────┴───────────────────┴─────────────────┴──────────┴───────────┴───────┘
```
    *   If we still pass these filtered antennas to the #objective_function,
        SciPy's minimizer would act as if these antennas are relevant and try to still modify their
        phase centers during the minization.
    """
    a1 = (
        phase_center_pairs
        .select(pl.col(r"^A1_(?:AntNum|Ring|Rho|Phi|Z)(?: \‍(placeholder\‍))?$"))
        .rename(lambda name: re.sub(r"^A1_(.*)$", r"\1", name))
    )
    a2 = (
        phase_center_pairs
        .select(pl.col(r"^A2_(?:AntNum|Ring|Rho|Phi|Z)(?: \‍(placeholder\‍))?$"))
        .rename(lambda name: re.sub(r"^A2_(.*)$", r"\1", name))
    )
 
    relevant_antennas = (
        a1.vstack(a2).unique().sort("AntNum")
    )
    return relevant_antennas
 
 
def objective_function(phase_center_guess: pl.Series,
                       phase_center_placeholders: pl.Series,
                       pulser_directions_and_measured_time_delays: pl.DataFrame,
                       axis: str) -> int:
    r"""!Computes the expected time delays and compares them with the measured time delays,
    designed to be run by SciPy's minimizer.
    @ingroup CALPC
    @param[in] phase_center_guess               The output of #retrieve_relevant_antennas_from_all_pairs
    @param[in] phase_center_placeholders        The output of #retrieve_relevant_antennas_from_all_pairs
    @param[in] pulser_directions_and_measured_time_delays  The output of #load_dataset_and_measure_time_delays
    @param[in] axis                             Axis to fit: \f$\rho\f$, \f$\phi\f$, or \f$z\f$
    @retval    total_difference_bw_expected_time_delay_and_measured_time_delay  A scalar
 
    * In #load_dataset_and_measure_time_delays(), the `measured_time_delay` has already been computed.
    * The job of this function is to
        -#  Given `phase_center_guess`, replace the placeholder antenna phase centers in
            `pulser_directions_and_pairwise_time_delays`
        -#  Compute the pairwise antenna displacement vectors (call these \f$\bf{d}\f$),
            using the updated antenna positions.
        -#  Compute the `expected time delay` based on the dot product between \f$\bf{d}\f$ and
            `incoming pulse direction (Phi Theta)`,
            see [this figure](@ref timedelaybydotproduct) for an illustration.
 ```
 ┌───────────────────────────┬───────────────────────────┐
 │ measured time delays (ns) ┆ expected time delays (ns) │
 │ ---                       ┆ ---                       │
 │ f64                       ┆ f64                       │
 ╞═══════════════════════════╪═══════════════════════════╡
 │ -0.866667                 ┆ -0.853384                 │
 │ -0.866667                 ┆ -0.853384                 │
 │ -7.466667                 ┆ -7.487111                 │
 │ 0.6                       ┆ 0.611216                  │
 │ -9.466667                 ┆ -9.48138                  │
 │ …                         ┆ …                         │
 │ 8.333333                  ┆ 8.33698                   │
 │ 0.2                       ┆ 0.15054                   │
 │ -0.866667                 ┆ -0.855732                 │
 │ 0.8                       ┆ 0.760589                  │
 │ 1.466667                  ┆ 1.463433                  │
 └───────────────────────────┴───────────────────────────┘
 ```
        -#  Quantify the difference between the expectation and the measurement; this is done by
            summing the difference over all rows above.
    * SciPy's minimizer tries to minimize the total difference by updating `phase_center_guess`.
    """
    from scipy.constants import c as speed_of_light
 
    phi_p = pl.col("incoming pulse direction (Phi Theta)").arr.get(0)
    theta_p = pl.col("incoming pulse direction (Phi Theta)").arr.get(1)
    x_p = phi_p.cos() * theta_p.sin()
    y_p = phi_p.sin() * theta_p.sin()
    z_p = theta_p.cos()
 
    r_1 = pl.col("A1_Rho")
    r_2 = pl.col("A2_Rho")
    z_1 = pl.col("A1_Z")
    z_2 = pl.col("A2_Z")
    phi_1 = pl.col("A1_Phi")
    phi_2 = pl.col("A2_Phi")
 
    A1minusA2_dotProduct_signalDir = (
        x_p * (r_1 * phi_1.cos() - r_2 * phi_2.cos()) +
        y_p * (r_1 * phi_1.sin() - r_2 * phi_2.sin()) + z_p * (z_1 - z_2)
    )
 
    total_difference_bw_expected_time_delay_and_measured_time_delay = (
        pulser_directions_and_measured_time_delays
        .with_columns(
            pl.col(f"A1_{axis} (placeholder)").replace_strict(phase_center_placeholders, phase_center_guess).alias(f"A1_{axis}"),
            pl.col(f"A2_{axis} (placeholder)").replace_strict(phase_center_placeholders, phase_center_guess).alias(f"A2_{axis}")
        )
        .with_columns(
            (A1minusA2_dotProduct_signalDir / speed_of_light * 1e9).alias("expected time delays (ns)")
        )
        .select(
            (pl.col("measured time delays (ns)") - pl.col("expected time delays (ns)"))
            .abs().sum().alias("time delay difference")
        )
    )
    return total_difference_bw_expected_time_delay_and_measured_time_delay.item()
 
 
def _plot_calibration_result(axis: str, before: pl.DataFrame, after: np.ndarray):
    """!
    @param[in] axis     Axis to fit: \f$\rho\f$, \f$\phi\f$, or \f$z\f$
    @param[in] before   The output of of #retrieve_relevant_antennas_from_all_pairs
    @param[in] after    The `restul.x` of running `scipy.optimize.minimize` on #objective_function
    @ingroup CALPC
    Nothing crazy here!
 
    ![](example_phase_center_fit_result.png)
    """
    import matplotlib.pyplot as plt
    plt.rcParams.update({'font.size': 13})
    fig1, ax1 = plt.subplots(nrows=4, figsize=(10, 15), tight_layout=True)
 
    before_and_after = (
        before.with_columns(
            pl.Series(after).alias(f"{axis} minimization result")
        )
    )
 
    for i in range(before['Ring'].unique().len()):
        ring = before_and_after.filter(pl.col("Ring") == (i + 1))
 
        ax1[i].scatter(ring['AntNum'], ring[f"{axis}"] * 1e2,
                       label=f"{axis} initial guess (face center from qrh.dat)")
        ax1[i].scatter(ring['AntNum'], ring[f"{axis} minimization result"] * 1e2,
                       label=f"{axis} minimization result")
        ax1[i].set_title(f"Ring {i + 1}")
        ax1[i].tick_params(axis='x', labelrotation=45, labelright=True)
        ax1[i].tick_params(axis='y', labelright=True)
        ax1[i].grid(alpha=.3)
    ax1[0].legend()
 
    fig1.supylabel(f"{axis} [cm]")
    fig1.supxlabel('Antenna Label')
    fig1.suptitle('Antenna Phase Centers (Payload Coordinates)\n')
    fig1.savefig("phase_center_fit_result.png")
 
 
if __name__ == "__main__":
    import os
    from scipy.optimize import minimize
 
    # use nominal phase center as the initial guesss of the minimzation
    j25: Path = os.environ.get("PUEO_UTIL_INSTALL_DIR") / Path("share/pueo/geometry/jun25/qrh.dat")
    nominal_pairs = generate_antenna_phase_center_pairs_with_placeholders(qrh_dot_dat=j25)
 
    # start with some run number to initialize the Dataset instance.
    ds: ROOT.pueo.Dataset = load_pueoEvent_Dataset(pueo_mc_data=Path("/tmp"), run_number=2)
    pulser_directions_and_pairwise_time_delays = (
        load_dataset_and_measure_time_delays(dataset=ds, phase_center_pairs=nominal_pairs)
        .filter(pl.col("max correlation") > 0.8)
    )
 
    relevant_antennas = (
        retrieve_relevant_antennas_from_all_pairs(
            phase_center_pairs=pulser_directions_and_pairwise_time_delays
        )
    )
 
    AXIS = "Rho"  # "Rho", "Phi", or "Z"
    result = minimize(fun=objective_function,
                      x0=relevant_antennas[AXIS],
                      args=(relevant_antennas[f"{AXIS} (placeholder)"],
                            pulser_directions_and_pairwise_time_delays, AXIS),
                      method='Nelder-Mead',
                      options={'maxiter': 500_000})
 
    print(result)
    _plot_calibration_result(axis=AXIS, before=relevant_antennas, after=result.x)