diff --git a/data_processing.ipynb b/data_processing.ipynb
index 8691005..a8a1939 100644
--- a/data_processing.ipynb
+++ b/data_processing.ipynb
@@ -172,7 +172,7 @@
     "print(f\"path_data: {path_data}\")\n",
     "path_qc = os.path.join(path_data, \"qc\")\n",
     "shim_modes = [\"CP\", \"patient\", \"volume\", \"phase\", \"CoV\", \"target\", \"SAReff\"]\n",
-    "shim_modes_MPRAGE = [\"CP\", \"CoV\"]\n",
+    "shim_modes_MPRAGE = [\"CP\", \"CoV\"]  # TODO: change variable name PEP8\n",
     "# shim_modes = [\"CP\", \"CoV\"]  # debugging\n",
     "print(f\"shim_modes: {shim_modes}\")\n",
     "subjects = sorted(glob.glob(\"sub-*\"))\n",
@@ -183,120 +183,88 @@
     "os.makedirs(path_results, exist_ok=True)"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "e95d87a4-3194-4369-b6cb-933ad1d5cfa3",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# Run segmentation on GRE scan\n",
-    "# ℹ️ The \"CoV reduction\" RF shimming scenario was chosen as the segmentation baseline due to the more \n",
-    "# homogeneous signal intensity in the I-->S direction, which results in a better segmentation peformance\n",
-    "# in the C7-T2 region\n",
-    "\n",
-    "for subject in subjects:\n",
-    "    os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
-    "    # Use another syntax for sub-04. See: https://github.com/shimming-toolbox/rf-shimming-7t/issues/31\n",
-    "    if subject == 'sub-04':\n",
-    "        !sct_deepseg_sc -i {subject}_acq-CoV_T2starw.nii.gz -c t2 -qc {path_qc}\n",
-    "    else:\n",
-    "        !sct_deepseg_sc -i {subject}_acq-CoV_T2starw.nii.gz -c t2s -qc {path_qc}"
-   ]
-  },
   {
    "cell_type": "markdown",
-   "id": "aac9d889",
+   "id": "cc07e3e6",
    "metadata": {},
    "source": [
-    "## Verify QC report (GRE segmentation)\n",
-    "\n",
-    "Open the quality control (QC) report located under `ds004906/qc/index.html`. Make sure the spinal cord segmentations are correct before resuming the analysis.\n",
-    "\n",
-    ">If you run this notebook on Google Colab, you can skip as the QC report cannot easily be viewed from Google Colab. If you *really* want to see the QC report, you can create a cell where you zip the 'qc/' folder, then you can download it on your local station and open the index.html file. \n"
+    "## Process anat/T1w (MPRAGE)\n"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "862dc562",
+   "id": "49a99b31",
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Crop GRE scan for faster processing and better registration results\n",
+    "# Run segmentation on CoV_T1w scan\n",
     "\n",
     "for subject in subjects:\n",
     "    os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
-    "    !sct_crop_image -i {subject}_acq-CoV_T2starw.nii.gz -m {subject}_acq-CoV_T2starw_seg.nii.gz -dilate 20x20x0 -o {subject}_acq-CoV_T2starw_crop.nii.gz\n",
-    "    !sct_crop_image -i {subject}_acq-CoV_T2starw_seg.nii.gz -m {subject}_acq-CoV_T2starw_seg.nii.gz -dilate 20x20x0 -o {subject}_acq-CoV_T2starw_crop_seg.nii.gz"
+    "    !sct_deepseg_sc -i {subject}_acq-CoV_T1w.nii.gz -c t1 -qc {path_qc}"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "6137bb2b",
+   "id": "663bb98c",
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Label vertebrae on GRE scan\n",
+    "# Crop the images for faster processing\n",
     "\n",
-    "# Given the low resolution of the GRE scan, the automatic detection of C2-C3 disc is unreliable. Therefore we need to use the manual disc labels that are part of the dataset.\n",
     "for subject in subjects:\n",
     "    os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
-    "    fname_label_discs = os.path.join(path_data, \"derivatives\", \"labels\", subject, \"anat\", f\"{subject}_acq-CoV_T2starw_label-discs_dseg.nii.gz\")\n",
-    "    !sct_label_utils -i {subject}_acq-CoV_T2starw_crop_seg.nii.gz -disc {fname_label_discs} -o {subject}_acq-CoV_T2starw_crop_seg_labeled.nii.gz\n",
-    "    # Generate QC report to assess labeled segmentation\n",
-    "    !sct_qc -i {subject}_acq-CoV_T2starw_crop.nii.gz -s {subject}_acq-CoV_T2starw_crop_seg_labeled.nii.gz -p sct_label_vertebrae -qc {path_qc} -qc-subject {subject}"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "c46e8d9a",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# Register the other shim methods to the GRE CoV scan\n",
-    "\n",
-    "# TODO"
+    "    !sct_crop_image -i {subject}_acq-CoV_T1w.nii.gz -m {subject}_acq-CoV_T1w_seg.nii.gz -dilate 34x34x0 -o {subject}_acq-CoV_T1w_crop.nii.gz\n",
+    "    !sct_crop_image -i {subject}_acq-CoV_T1w_seg.nii.gz -m {subject}_acq-CoV_T1w_seg.nii.gz -dilate 34x34x0 -o {subject}_acq-CoV_T1w_crop_seg.nii.gz"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "90deb5e5",
+   "id": "75071e07",
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Extract the signal intensity on the GRE scan within the spinal cord between levels C3 and T2 (included), which correspond to the region where RF shimming was prescribed\n",
+    "# Vertebral labeling on the CoV_T1w scan\n",
     "\n",
     "for subject in subjects:\n",
-    "    # TODO: loop across other shim methods    \n",
     "    os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
-    "    !sct_extract_metric -i {subject}_acq-CoV_T2starw_crop.nii.gz -f {subject}_acq-CoV_T2starw_crop_seg.nii.gz -method wa -vert 3:9 -vertfile {subject}_acq-CoV_T2starw_crop_seg_labeled.nii.gz -append 1 -perslice 1 -o gre_CoV.csv\n"
+    "    !sct_label_vertebrae -i {subject}_acq-CoV_T1w_crop.nii.gz -s {subject}_acq-CoV_T1w_crop_seg.nii.gz -c t1 -qc {path_qc}\n",
+    "    \n",
+    "# TODO: replace by code below\n",
+    "# https://github.com/shimming-toolbox/rf-shimming-7t/issues/59\n",
+    "# Given the low resolution of the GRE scan, the automatic detection of C2-C3 disc is unreliable. Therefore we need to use the manual disc labels that are part of the dataset.\n",
+    "# for subject in subjects:\n",
+    "#     os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
+    "#     fname_label_discs = os.path.join(path_data, \"derivatives\", \"labels\", subject, \"anat\", f\"{subject}_acq-CoV_T2starw_label-discs_dseg.nii.gz\")\n",
+    "#     !sct_label_utils -i {subject}_acq-CoV_T2starw_crop_seg.nii.gz -disc {fname_label_discs} -o {subject}_acq-CoV_T2starw_crop_seg_labeled.nii.gz\n",
+    "#     # Generate QC report to assess labeled segmentation\n",
+    "#     !sct_qc -i {subject}_acq-CoV_T2starw_crop.nii.gz -s {subject}_acq-CoV_T2starw_crop_seg_labeled.nii.gz -p sct_label_vertebrae -qc {path_qc} -qc-subject {subject}"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "e0cde907",
+   "id": "436a2728",
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Register TFL flip angle maps to the GRE scan ⏳\n",
+    "# Register the CP_T1w with the CoV_T1w\n",
     "\n",
     "for subject in subjects:\n",
-    "    os.chdir(os.path.join(path_data, subject, \"fmap\"))\n",
-    "    for shim_mode in shim_modes:\n",
-    "        !sct_register_multimodal -i {subject}_acq-anat{shim_mode}_TB1TFL.nii.gz -d ../anat/{subject}_acq-CoV_T2starw_crop.nii.gz -dseg ../anat/{subject}_acq-CoV_T2starw_crop_seg.nii.gz -param step=1,type=im,algo=slicereg,metric=CC -qc {path_qc}"
+    "    os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
+    "    !sct_register_multimodal -i {subject}_acq-CP_T1w.nii.gz -d {subject}_acq-CoV_T1w_crop.nii.gz -dseg {subject}_acq-CoV_T1w_crop_seg.nii.gz -qc {path_qc}"
    ]
   },
   {
    "cell_type": "markdown",
-   "id": "39de4623",
+   "id": "a2074892",
    "metadata": {},
    "source": [
-    "## Verify QC report (B1maps to GRE registration)\n",
+    "### Verify QC report (T1w segmentation, labeling and registration)\n",
+    "\n",
     "\n",
     "Open the QC report located under `ds004906/qc/index.html`. Make sure the registration are correct before resuming the analysis."
    ]
@@ -304,93 +272,91 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "6a6d7670",
+   "id": "75c43da0",
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Warping spinal cord segmentation and vertebral level to each flip angle map\n",
+    "# Extract the signal intensity on the MPRAGE scan within the spinal cord between levels C3 and T2 (included)\n",
     "\n",
-    "for subject in subjects:\n",
-    "    os.chdir(os.path.join(path_data, subject, \"fmap\"))\n",
-    "    for shim_mode in shim_modes:\n",
-    "        # Use linear interpolation to preserve partial volume information (needed when extracting values along the cord)\n",
-    "        !sct_apply_transfo -i ../anat/{subject}_acq-CoV_T2starw_crop_seg.nii.gz -d {subject}_acq-anat{shim_mode}_TB1TFL.nii.gz -w warp_{subject}_acq-CoV_T2starw_crop2{subject}_acq-anat{shim_mode}_TB1TFL.nii.gz -x linear -o {subject}_acq-anat{shim_mode}_TB1TFL_seg.nii.gz\n",
-    "        # Use nearest neighbour (nn) interpolation because we are dealing with non-binary discrete labels\n",
-    "        !sct_apply_transfo -i ../anat/{subject}_acq-CoV_T2starw_crop_seg_labeled.nii.gz -d {subject}_acq-anat{shim_mode}_TB1TFL.nii.gz -w warp_{subject}_acq-CoV_T2starw_crop2{subject}_acq-anat{shim_mode}_TB1TFL.nii.gz -x nn -o {subject}_acq-anat{shim_mode}_TB1TFL_seg_labeled.nii.gz"
+    "for subject in subjects:    \n",
+    "    os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
+    "    #CP\n",
+    "    !sct_extract_metric -i {subject}_acq-CP_T1w_reg.nii.gz -f {subject}_acq-CoV_T1w_crop_seg.nii.gz -method wa -vert 3:9 -vertfile {subject}_acq-CoV_T1w_crop_seg_labeled.nii.gz -append 1 -perslice 1 -o {subject}_acq-CP_T1w.csv\n",
+    "    #CoV\n",
+    "    !sct_extract_metric -i {subject}_acq-CoV_T1w_crop.nii.gz -f {subject}_acq-CoV_T1w_crop_seg.nii.gz -method wa -vert 3:9 -vertfile {subject}_acq-CoV_T1w_crop_seg_labeled.nii.gz -append 1 -perslice 1 -o {subject}_acq-CoV_T1w.csv"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "a835cdee",
+   "id": "63e6181a",
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Convert the flip angle maps to B1+ efficiency maps [nT/V] (inspired by code from Kyle Gilbert)\n",
-    "# The approach consists in calculating the B1+ efficiency using a 1ms, pi-pulse at the acquisition voltage,\n",
-    "# then scale the efficiency by the ratio of the measured flip angle to the requested flip angle in the pulse sequence.\n",
-    "\n",
-    "GAMMA = 2.675e8;  # [rad / (s T)]\n",
-    "requested_fa = 90  # saturation flip angle -- hard-coded in sequence\n",
-    "\n",
+    "# Flattening sagittaly\n",
     "for subject in subjects:\n",
-    "    os.chdir(os.path.join(path_data, subject, \"fmap\"))\n",
-    "    for shim_mode in shim_modes:\n",
-    "        # Fetch the reference voltage from the JSON sidecar to the TFL B1map sequence\n",
-    "        with open(f\"{subject}_acq-famp{shim_mode}_TB1TFL.json\", \"r\") as f:\n",
-    "            metadata = json.load(f)\n",
-    "            ref_voltage = metadata.get(\"TxRefAmp\", \"N/A\")\n",
-    "            print(f\"ref_voltage [V]: {ref_voltage} ({subject}_acq-famp{shim_mode}_TB1TFL)\")\n",
-    "\n",
-    "        # Open flip angle map with nibabel\n",
-    "        nii = nib.load(f\"{subject}_acq-famp{shim_mode}_TB1TFL.nii.gz\")\n",
-    "        acquired_fa = nii.get_fdata()\n",
-    "\n",
-    "        # Siemens maps are in units of flip angle * 10 (in degrees)\n",
-    "        acquired_fa = acquired_fa / 10\n",
-    "\n",
-    "        # Account for the power loss between the coil and the socket. That number was given by Siemens.\n",
-    "        voltage_at_socket = ref_voltage * 10 ** -0.095\n",
-    "\n",
-    "        # Compute B1 map in [T/V]\n",
-    "        b1_map = (acquired_fa / requested_fa) * (np.pi / (GAMMA * 1e-3 * voltage_at_socket))\n",
-    "\n",
-    "        # Convert to [nT/V]\n",
-    "        b1_map = b1_map * 1e9\n",
-    "\n",
-    "        # Save as NIfTI file\n",
-    "        nii_b1 = nib.Nifti1Image(b1_map, nii.affine, nii.header)\n",
-    "        nib.save(nii_b1, f\"{subject}_acq-{shim_mode}_TB1map.nii.gz\")"
+    "    os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
+    "    #flatten \n",
+    "    !sct_flatten_sagittal -i {subject}_acq-CoV_T1w_crop.nii.gz -s {subject}_acq-CoV_T1w_crop_seg.nii.gz\n",
+    "    !sct_flatten_sagittal -i {subject}_acq-CP_T1w_reg.nii.gz -s {subject}_acq-CoV_T1w_crop_seg.nii.gz"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "740fda91",
+   "id": "77d24f64",
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Extract B1+ value along the spinal cord between levels C3 and T2 (included)\n",
+    "# Display images\n",
     "\n",
     "for subject in subjects:\n",
-    "    os.chdir(os.path.join(path_data, subject, \"fmap\"))\n",
-    "    for shim_mode in shim_modes:\n",
-    "        !sct_extract_metric -i {subject}_acq-{shim_mode}_TB1map.nii.gz -f {subject}_acq-anat{shim_mode}_TB1TFL_seg.nii.gz -method wa -vert 3:9 -vertfile {subject}_acq-anat{shim_mode}_TB1TFL_seg_labeled.nii.gz -perslice 1 -o TB1map_{shim_mode}.csv"
+    "    os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
+    "    CoV_flatten = os.path.join(path_data, subject, \"anat\", f\"{subject}_acq-CoV_T1w_crop_flatten.nii.gz\")\n",
+    "    CoV_flatten_data = nib.load(CoV_flatten).get_fdata()\n",
+    "    CP_flatten = os.path.join(path_data, subject, \"anat\", f\"{subject}_acq-CP_T1w_reg_flatten.nii.gz\")\n",
+    "    CP_flatten_data=nib.load(CP_flatten).get_fdata()\n",
+    "\n",
+    "    # Data manipulation to get the central slice\n",
+    "    # Already cropped inf-sup and right-left, no need to crop more\n",
+    "    CoV_flatten_data_cslice=CoV_flatten_data[CoV_flatten_data.shape[0] // 2,:,:]\n",
+    "    CP_flatten_data_cslice=CP_flatten_data[CP_flatten_data.shape[0] // 2,:,:]\n",
+    "    #And now plot side by side\n",
+    "    plt.figure(figsize=(3, 4))\n",
+    "    vmin = -1  \n",
+    "    vmax = -0.70\n",
+    "\n",
+    "    plt.subplot(1, 2, 1)\n",
+    "    plt.imshow(CP_flatten_data_cslice.T, cmap='gray', origin='lower',vmin=vmin,vmax=vmax)  \n",
+    "    plt.title('CP mode',color='blue')\n",
+    "    plt.axis('off')\n",
+    "\n",
+    "    plt.subplot(1, 2, 2)\n",
+    "    plt.imshow(CoV_flatten_data_cslice.T, cmap='gray', origin='lower',vmin=vmin, vmax=vmax)\n",
+    "    plt.title('CoV reduction',color='orange')\n",
+    "    plt.axis('off')\n",
+    "\n",
+    "    # Add the subject title at the top\n",
+    "    plt.suptitle(subject, y=0.9)  # Adjust 'y' as needed for position\n",
+    "    plt.tight_layout()\n",
+    "    plt.subplots_adjust(wspace=0.1, hspace=0)\n",
+    "    plt.show()    "
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "b0185bb0",
+   "id": "cd2a3c9c",
    "metadata": {},
    "outputs": [],
    "source": [
+    "# Generate signal intensity curves\n",
     "# Make figure of B1+ values along the spinal cord across shim methods\n",
     "\n",
     "# Go back to root data folder\n",
     "os.chdir(os.path.join(path_data))\n",
     "\n",
-    "def smooth_data(data, window_size=20):\n",
+    "def smooth_data(data, window_size=11):\n",
     "    \"\"\" Apply a simple moving average to smooth the data. \"\"\"\n",
     "    return uniform_filter1d(data, size=window_size, mode='nearest')\n",
     "\n",
@@ -398,7 +364,7 @@
     "x_grid = np.linspace(0, 1, 100)\n",
     "\n",
     "# z-slices corresponding to levels C1 to T2 on the PAM50 template. These will be used to scale the x-label of each subject.\n",
-    "original_vector = np.array([984, 938, 907, 870, 833, 800, 769, 735, 692, 646])\n",
+    "original_vector =np.array([907, 870, 833, 800, 769, 735, 692, 646]) # np.array([984, 938, 907, 870, 833, 800, 769, 735, 692, 646])\n",
     "\n",
     "# Normalize the PAM50 z-slice numbers to the 1-0 range (to show inferior-superior instead of superior-inferior)\n",
     "min_val = original_vector.min()\n",
@@ -423,13 +389,13 @@
     "for i, subject in enumerate(subjects):\n",
     "    ax = axes[i]\n",
     "    \n",
-    "    os.chdir(os.path.join(path_data, subject, \"fmap\"))\n",
+    "    os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
     "\n",
-    "    for shim_method in shim_modes:\n",
+    "    for shim_method in shim_modes_MPRAGE:\n",
     "        # Initialize list to collect data for this shim method\n",
     "        method_data = []\n",
     "\n",
-    "        for file_path in glob.glob(f\"TB1map_*{shim_method}*.csv\"):\n",
+    "        for file_path in glob.glob(f\"{subject}_acq-{shim_method}_T1w.csv\"):\n",
     "            df = pd.read_csv(file_path)\n",
     "            wa_data = df['WA()']\n",
     "\n",
@@ -456,7 +422,7 @@
     "    ax.set_xticklabels([''] * len(original_vector))\n",
     "\n",
     "    ax.set_title(f'{subject}', fontsize=font_size)\n",
-    "    ax.set_ylabel('B1+ efficiency [nT/V]', fontsize=font_size)\n",
+    "    ax.set_ylabel('Signal intensity [a.u.]', fontsize=font_size)\n",
     "    ax.tick_params(axis='y', which='major', labelsize=font_size-4)\n",
     "\n",
     "    # Add legend only to the first subplot\n",
@@ -473,11 +439,11 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "f8b90a09",
+   "id": "1a690f03",
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Create tables and perform statistics\n",
+    "# Create mean and CoV \n",
     "\n",
     "# Go back to root data folder\n",
     "os.chdir(path_data)\n",
@@ -485,242 +451,270 @@
     "# Data storage\n",
     "data_summary = []\n",
     "\n",
-    "# Compute mean and SD of B1+ for each subject and each shim mode\n",
+    "# Compute mean and SD of MPRAGE sigint for each subject and each shim mode\n",
     "for subject in subjects:\n",
-    "    for shim_mode in shim_modes:\n",
+    "    for shim_mode in shim_modes_MPRAGE:\n",
     "        all_data = []\n",
-    "        for file_path in glob.glob(os.path.join(path_data, subject, \"fmap\", f\"TB1map_*{shim_mode}*.csv\")):\n",
+    "        for file_path in glob.glob(os.path.join(path_data, subject, \"anat\", f\"{subject}_acq-{shim_mode}_T1w.csv\")):\n",
+    "           # print(file_path)\n",
     "            df = pd.read_csv(file_path)\n",
-    "            # Drop rows where VertLevel is 1 or 2 (because RF shimming was not performed at these levels)\n",
-    "            df = df[df['VertLevel'].isin(['1', '2']) == False]\n",
     "            wa_data = df['WA()']\n",
     "            all_data.extend(wa_data)\n",
     "\n",
     "        if all_data:\n",
     "            mean_data = np.mean(all_data)\n",
     "            sd_data = np.std(all_data)\n",
-    "            data_summary.append([subject, shim_mode, mean_data, sd_data])\n",
+    "            CoV_data=np.round((sd_data/mean_data)*100,2)\n",
+    "            data_summary.append([subject, shim_mode, mean_data, sd_data, CoV_data])\n",
     "\n",
     "# Convert to DataFrame and save to CSV\n",
-    "df_summary = pd.DataFrame(data_summary, columns=['Subject', 'Shim_Mode', 'Average', 'Standard_Deviation'])\n",
+    "df_summary = pd.DataFrame(data_summary, columns=['Subject', 'Shim_Mode', 'Average', 'Standard_Deviation', 'CoV'])\n",
     "df_summary.to_csv(os.path.join(path_results, 'subject_shim_mode_summary.csv'), index=False)\n",
     "\n",
-    "# Compute statistics across subjects\n",
-    "grouped_means = df_summary.groupby('Shim_Mode').agg({'Average': ['mean', 'std'], 'Standard_Deviation': ['mean', 'std']})\n",
-    "\n",
-    "# Format mean ± standard deviation\n",
-    "grouped_means['Average_formatted'] = grouped_means['Average']['mean'].map(\"{:.2f}\".format) + \" ± \" + grouped_means['Average']['std'].map(\"{:.2f}\".format)\n",
-    "grouped_means['Standard_Deviation_formatted'] = grouped_means['Standard_Deviation']['mean'].map(\"{:.2f}\".format) + \" ± \" + grouped_means['Standard_Deviation']['std'].map(\"{:.2f}\".format)\n",
-    "\n",
-    "# Drop multi-level index and only keep formatted columns\n",
-    "grouped_means = grouped_means.drop(columns=['Average', 'Standard_Deviation'])\n",
-    "grouped_means.columns = ['Average', 'Standard_Deviation']  # Rename columns for clarity\n",
-    "grouped_means.reset_index().to_csv(os.path.join(path_results, 'average_across_subjects.csv'), index=False)\n",
-    "\n",
-    "# Perform ANOVA and Posthoc Tests\n",
-    "anova_result = f_oneway(*[group[\"Average\"].values for name, group in df_summary.groupby(\"Shim_Mode\")])\n",
-    "print(\"ANOVA Result:\", anova_result)\n",
-    "\n",
-    "if anova_result.pvalue < 0.05:\n",
-    "    posthoc_result = pairwise_tukeyhsd(df_summary['Average'], df_summary['Shim_Mode'])\n",
-    "    print(\"Posthoc Tukey HSD Result:\\n\", posthoc_result)\n"
+    "df_summary_onlyCoV=df_summary[['Subject', 'Shim_Mode', 'CoV']]\n",
+    "print(df_summary_onlyCoV)\n"
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "a916a6f7",
+   "cell_type": "markdown",
+   "id": "b0c4e225",
    "metadata": {},
-   "outputs": [],
    "source": [
-    "# Visualize RF maps to generate Figure 2\n",
-    "\n",
-    "# TODO"
+    "## Process anat/T2starw (GRE)"
    ]
   },
   {
-   "cell_type": "markdown",
-   "id": "cc07e3e6",
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e95d87a4-3194-4369-b6cb-933ad1d5cfa3",
    "metadata": {},
+   "outputs": [],
    "source": [
-    "# MPRAGE PROCESSING BLOCK\n"
+    "# Run segmentation on GRE scan\n",
+    "# ℹ️ The \"CoV reduction\" RF shimming scenario was chosen as the segmentation baseline due to the more \n",
+    "# homogeneous signal intensity in the I-->S direction, which results in a better segmentation peformance\n",
+    "# in the C7-T2 region\n",
+    "\n",
+    "for subject in subjects:\n",
+    "    os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
+    "    # Use another syntax for sub-04. See: https://github.com/shimming-toolbox/rf-shimming-7t/issues/31\n",
+    "    if subject == 'sub-04':\n",
+    "        !sct_deepseg_sc -i {subject}_acq-CoV_T2starw.nii.gz -c t2 -qc {path_qc}\n",
+    "    else:\n",
+    "        !sct_deepseg_sc -i {subject}_acq-CoV_T2starw.nii.gz -c t2s -qc {path_qc}"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "49a99b31",
+   "id": "862dc562",
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Run segmentation on MPRAGE CoV scan\n",
+    "# Crop GRE scan for faster processing and better registration results\n",
     "\n",
     "for subject in subjects:\n",
     "    os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
-    "    !sct_deepseg_sc -i {subject}_acq-CoV_T1w.nii.gz -c t1 -qc {path_qc}"
+    "    !sct_crop_image -i {subject}_acq-CoV_T2starw.nii.gz -m {subject}_acq-CoV_T2starw_seg.nii.gz -dilate 20x20x0 -o {subject}_acq-CoV_T2starw_crop.nii.gz\n",
+    "    !sct_crop_image -i {subject}_acq-CoV_T2starw_seg.nii.gz -m {subject}_acq-CoV_T2starw_seg.nii.gz -dilate 20x20x0 -o {subject}_acq-CoV_T2starw_crop_seg.nii.gz"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "663bb98c",
+   "id": "87030540",
    "metadata": {},
    "outputs": [],
    "source": [
-    "#Crop the scans for easier and faster coregistration\n",
-    "#Using a more permissive bounding box than GRE\n",
+    "# Register CoV_T2starw to CoV_T1w to bring vertebral labels\n",
     "\n",
     "for subject in subjects:\n",
-    "    os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
-    "    !sct_crop_image -i {subject}_acq-CoV_T1w.nii.gz -m {subject}_acq-CoV_T1w_seg.nii.gz -dilate 34x34x0 -o {subject}_acq-CoV_T1w_crop.nii.gz\n",
-    "    !sct_crop_image -i {subject}_acq-CoV_T1w_seg.nii.gz -m {subject}_acq-CoV_T1w_seg.nii.gz -dilate 34x34x0 -o {subject}_acq-CoV_T1w_crop_seg.nii.gz"
+    "    os.chdir(os.path.join(path_data, subject, \"anat\"))    \n",
+    "    !sct_register_multimodal -i {subject}_acq-CoV_T1w_crop.nii.gz -iseg {subject}_acq-CoV_T1w_crop_seg.nii.gz -d {subject}_acq-CoV_T2starw_crop.nii.gz -dseg {subject}_acq-CoV_T2starw_crop_seg.nii.gz -param step=1,type=seg,algo=slicereg -qc {path_qc}"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "75071e07",
+   "id": "c46e8d9a",
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Vertebral labeling on the CoV MPRAGE scan\n",
+    "# Register *_T2starw to CoV_T2starw\n",
     "\n",
-    "for subject in subjects:\n",
-    "    os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
-    "    !sct_label_vertebrae -i {subject}_acq-CoV_T1w_crop.nii.gz -s {subject}_acq-CoV_T1w_crop_seg.nii.gz -c t1 -qc {path_qc}"
+    "# Commenting for now, due to https://github.com/shimming-toolbox/rf-shimming-7t/issues/35\n",
+    "# for subject in subjects:\n",
+    "#     os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
+    "#     for shim_mode in shim_modes:\n",
+    "#         # Don't do it for CoV_T2starw\n",
+    "#         if shim_mode != 'CoV':\n",
+    "#             !sct_register_multimodal -i {subject}_acq-{shim_mode}_T2starw.nii.gz -d {subject}_acq-CoV_T2starw_crop.nii.gz -dseg {subject}_acq-CoV_T2starw_crop_seg.nii.gz -param step=1,type=im,algo=slicereg,metric=CC -qc {path_qc}"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "436a2728",
+   "id": "0ffcfb28",
    "metadata": {},
    "outputs": [],
    "source": [
-    "#Coregistering the CP mode to the CoV case\n",
-    "#Seems to be working good enough\n",
+    "# Warp vertebral level to *_T2starw\n",
+    "\n",
     "for subject in subjects:\n",
-    "    os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
-    "    !sct_register_multimodal -i {subject}_acq-CP_T1w.nii.gz -d {subject}_acq-CoV_T1w_crop.nii.gz -dseg {subject}_acq-CoV_T1w_crop_seg.nii.gz -qc {path_qc}"
+    "    os.chdir(os.path.join(path_data, subject, \"anat\"))    \n",
+    "    !sct_apply_transfo -i {subject}_acq-CoV_T1w_crop_seg_labeled.nii.gz -d {subject}_acq-CoV_T2starw_crop.nii.gz -w warp_{subject}_acq-CoV_T1w_crop2{subject}_acq-CoV_T2starw_crop.nii.gz -x nn -o {subject}_acq-CoV_T2starw_crop_seg_labeled.nii.gz"
    ]
   },
   {
    "cell_type": "markdown",
-   "id": "a2074892",
+   "id": "aac9d889",
    "metadata": {},
    "source": [
-    "## Verify QC report (MPRAGE segmentation and coregistration)\n",
+    "## Verify QC report (GRE segmentation)\n",
     "\n",
+    "Open the quality control (QC) report located under `ds004906/qc/index.html`. Make sure the spinal cord segmentations are correct before resuming the analysis.\n",
     "\n",
-    "Open the QC report located under `ds004906/qc/index.html`. Make sure the registration are correct before resuming the analysis."
+    ">If you run this notebook on Google Colab, you can skip as the QC report cannot easily be viewed from Google Colab. If you *really* want to see the QC report, you can create a cell where you zip the 'qc/' folder, then you can download it on your local station and open the index.html file. \n"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "0ffcfb28",
+   "id": "90deb5e5",
    "metadata": {},
    "outputs": [],
    "source": [
-    "# UNNECESSARY, but keeping code it in case we want to register the GRE scans to the T1w scans.\n",
-    "# https://github.com/shimming-toolbox/rf-shimming-7t/issues/49\n",
-    "\n",
-    "# Warping spinal cord segmentation and vertebral level to the CP mode scan for each subject\n",
+    "# Extract the signal intensity on the GRE scan within the spinal cord between levels C3 and T2 (included), which correspond to the region where RF shimming was prescribed\n",
     "\n",
-    "# for subject in subjects:\n",
-    "#     os.chdir(os.path.join(path_data, subject, \"anat\"))    \n",
-    "#     !sct_apply_transfo -i {subject}_acq-CoV_T1w_crop_seg.nii.gz -d {subject}_acq-CP_T1w.nii.gz  -w warp_{subject}_acq-CoV_T1w_crop2{subject}_acq-CP_T1w.nii.gz -x linear -o {subject}_acq-CP_T1w_seg_reg.nii.gz\n",
-    "#     !sct_apply_transfo -i {subject}_acq-CoV_T1w_crop_seg_labeled.nii.gz -d {subject}_acq-CP_T1w.nii.gz  -w warp_{subject}_acq-CoV_T1w_crop2{subject}_acq-CP_T1w.nii.gz -x linear -o {subject}_acq-CP_T1w_seg_labeled_reg.nii.gz"
+    "for subject in subjects:\n",
+    "    # TODO: loop across other shim methods    \n",
+    "    os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
+    "    !sct_extract_metric -i {subject}_acq-CoV_T2starw_crop.nii.gz -f {subject}_acq-CoV_T2starw_crop_seg.nii.gz -method wa -vert 3:9 -vertfile {subject}_acq-CoV_T2starw_crop_seg_labeled.nii.gz -append 1 -perslice 1 -o gre_CoV.csv\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "bda7c3ac",
+   "metadata": {},
+   "source": [
+    "## Process fmap/TFL (flip angle maps)"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "75c43da0",
+   "id": "e0cde907",
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Extract the signal intensity on the MPRAGE scan within the spinal cord between levels C3 and T2 (included)\n",
+    "# Register TFL flip angle maps to the GRE scan ⏳\n",
     "\n",
-    "for subject in subjects:    \n",
-    "    os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
-    "    #CP\n",
-    "    !sct_extract_metric -i {subject}_acq-CP_T1w_reg.nii.gz -f {subject}_acq-CoV_T1w_crop_seg.nii.gz -method wa -vert 3:9 -vertfile {subject}_acq-CoV_T1w_crop_seg_labeled.nii.gz -append 1 -perslice 1 -o {subject}_acq-CP_T1w.csv\n",
-    "    #CoV\n",
-    "    !sct_extract_metric -i {subject}_acq-CoV_T1w_crop.nii.gz -f {subject}_acq-CoV_T1w_crop_seg.nii.gz -method wa -vert 3:9 -vertfile {subject}_acq-CoV_T1w_crop_seg_labeled.nii.gz -append 1 -perslice 1 -o {subject}_acq-CoV_T1w.csv"
+    "for subject in subjects:\n",
+    "    os.chdir(os.path.join(path_data, subject, \"fmap\"))\n",
+    "    for shim_mode in shim_modes:\n",
+    "        !sct_register_multimodal -i {subject}_acq-anat{shim_mode}_TB1TFL.nii.gz -d ../anat/{subject}_acq-CoV_T2starw_crop.nii.gz -dseg ../anat/{subject}_acq-CoV_T2starw_crop_seg.nii.gz -param step=1,type=im,algo=slicereg,metric=CC -qc {path_qc}"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "39de4623",
+   "metadata": {},
+   "source": [
+    "## Verify QC report (B1maps to GRE registration)\n",
+    "\n",
+    "Open the QC report located under `ds004906/qc/index.html`. Make sure the registration are correct before resuming the analysis."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "63e6181a",
+   "id": "6a6d7670",
    "metadata": {},
    "outputs": [],
    "source": [
-    "#Flattening sagittaly\n",
+    "# Warping spinal cord segmentation and vertebral level to each flip angle map\n",
+    "\n",
     "for subject in subjects:\n",
-    "    os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
-    "    #flatten \n",
-    "    !sct_flatten_sagittal -i {subject}_acq-CoV_T1w_crop.nii.gz -s {subject}_acq-CoV_T1w_crop_seg.nii.gz\n",
-    "    !sct_flatten_sagittal -i {subject}_acq-CP_T1w_reg.nii.gz -s {subject}_acq-CoV_T1w_crop_seg.nii.gz"
+    "    os.chdir(os.path.join(path_data, subject, \"fmap\"))\n",
+    "    for shim_mode in shim_modes:\n",
+    "        # Use linear interpolation to preserve partial volume information (needed when extracting values along the cord)\n",
+    "        !sct_apply_transfo -i ../anat/{subject}_acq-CoV_T2starw_crop_seg.nii.gz -d {subject}_acq-anat{shim_mode}_TB1TFL.nii.gz -w warp_{subject}_acq-CoV_T2starw_crop2{subject}_acq-anat{shim_mode}_TB1TFL.nii.gz -x linear -o {subject}_acq-anat{shim_mode}_TB1TFL_seg.nii.gz\n",
+    "        # Use nearest neighbour (nn) interpolation because we are dealing with non-binary discrete labels\n",
+    "        !sct_apply_transfo -i ../anat/{subject}_acq-CoV_T2starw_crop_seg_labeled.nii.gz -d {subject}_acq-anat{shim_mode}_TB1TFL.nii.gz -w warp_{subject}_acq-CoV_T2starw_crop2{subject}_acq-anat{shim_mode}_TB1TFL.nii.gz -x nn -o {subject}_acq-anat{shim_mode}_TB1TFL_seg_labeled.nii.gz"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "77d24f64",
+   "id": "a835cdee",
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Display images\n",
+    "# Convert the flip angle maps to B1+ efficiency maps [nT/V] (inspired by code from Kyle Gilbert)\n",
+    "# The approach consists in calculating the B1+ efficiency using a 1ms, pi-pulse at the acquisition voltage,\n",
+    "# then scale the efficiency by the ratio of the measured flip angle to the requested flip angle in the pulse sequence.\n",
+    "\n",
+    "GAMMA = 2.675e8;  # [rad / (s T)]\n",
+    "requested_fa = 90  # saturation flip angle -- hard-coded in sequence\n",
     "\n",
     "for subject in subjects:\n",
-    "    os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
-    "    CoV_flatten = os.path.join(path_data, subject, \"anat\", f\"{subject}_acq-CoV_T1w_crop_flatten.nii.gz\")\n",
-    "    CoV_flatten_data = nib.load(CoV_flatten).get_fdata()\n",
-    "    CP_flatten = os.path.join(path_data, subject, \"anat\", f\"{subject}_acq-CP_T1w_reg_flatten.nii.gz\")\n",
-    "    CP_flatten_data=nib.load(CP_flatten).get_fdata()\n",
+    "    os.chdir(os.path.join(path_data, subject, \"fmap\"))\n",
+    "    for shim_mode in shim_modes:\n",
+    "        # Fetch the reference voltage from the JSON sidecar to the TFL B1map sequence\n",
+    "        with open(f\"{subject}_acq-famp{shim_mode}_TB1TFL.json\", \"r\") as f:\n",
+    "            metadata = json.load(f)\n",
+    "            ref_voltage = metadata.get(\"TxRefAmp\", \"N/A\")\n",
+    "            print(f\"ref_voltage [V]: {ref_voltage} ({subject}_acq-famp{shim_mode}_TB1TFL)\")\n",
     "\n",
-    "    # Data manipulation to get the central slice\n",
-    "    # Already cropped inf-sup and right-left, no need to crop more\n",
-    "    CoV_flatten_data_cslice=CoV_flatten_data[CoV_flatten_data.shape[0] // 2,:,:]\n",
-    "    CP_flatten_data_cslice=CP_flatten_data[CP_flatten_data.shape[0] // 2,:,:]\n",
-    "    #And now plot side by side\n",
-    "    plt.figure(figsize=(3, 4))\n",
-    "    vmin = -1  \n",
-    "    vmax = -0.70\n",
+    "        # Open flip angle map with nibabel\n",
+    "        nii = nib.load(f\"{subject}_acq-famp{shim_mode}_TB1TFL.nii.gz\")\n",
+    "        acquired_fa = nii.get_fdata()\n",
     "\n",
-    "    plt.subplot(1, 2, 1)\n",
-    "    plt.imshow(CP_flatten_data_cslice.T, cmap='gray', origin='lower',vmin=vmin,vmax=vmax)  \n",
-    "    plt.title('CP mode',color='blue')\n",
-    "    plt.axis('off')\n",
+    "        # Siemens maps are in units of flip angle * 10 (in degrees)\n",
+    "        acquired_fa = acquired_fa / 10\n",
     "\n",
-    "    plt.subplot(1, 2, 2)\n",
-    "    plt.imshow(CoV_flatten_data_cslice.T, cmap='gray', origin='lower',vmin=vmin, vmax=vmax)\n",
-    "    plt.title('CoV reduction',color='orange')\n",
-    "    plt.axis('off')\n",
+    "        # Account for the power loss between the coil and the socket. That number was given by Siemens.\n",
+    "        voltage_at_socket = ref_voltage * 10 ** -0.095\n",
     "\n",
-    "    # Add the subject title at the top\n",
-    "    plt.suptitle(subject, y=0.9)  # Adjust 'y' as needed for position\n",
-    "    plt.tight_layout()\n",
-    "    plt.subplots_adjust(wspace=0.1, hspace=0)\n",
-    "    plt.show()    "
+    "        # Compute B1 map in [T/V]\n",
+    "        b1_map = (acquired_fa / requested_fa) * (np.pi / (GAMMA * 1e-3 * voltage_at_socket))\n",
+    "\n",
+    "        # Convert to [nT/V]\n",
+    "        b1_map = b1_map * 1e9\n",
+    "\n",
+    "        # Save as NIfTI file\n",
+    "        nii_b1 = nib.Nifti1Image(b1_map, nii.affine, nii.header)\n",
+    "        nib.save(nii_b1, f\"{subject}_acq-{shim_mode}_TB1map.nii.gz\")"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "cd2a3c9c",
+   "id": "740fda91",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Extract B1+ value along the spinal cord between levels C3 and T2 (included)\n",
+    "\n",
+    "for subject in subjects:\n",
+    "    os.chdir(os.path.join(path_data, subject, \"fmap\"))\n",
+    "    for shim_mode in shim_modes:\n",
+    "        !sct_extract_metric -i {subject}_acq-{shim_mode}_TB1map.nii.gz -f {subject}_acq-anat{shim_mode}_TB1TFL_seg.nii.gz -method wa -vert 3:9 -vertfile {subject}_acq-anat{shim_mode}_TB1TFL_seg_labeled.nii.gz -perslice 1 -o TB1map_{shim_mode}.csv"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b0185bb0",
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Generate signal intensity curves\n",
     "# Make figure of B1+ values along the spinal cord across shim methods\n",
     "\n",
     "# Go back to root data folder\n",
     "os.chdir(os.path.join(path_data))\n",
     "\n",
-    "def smooth_data(data, window_size=11):\n",
+    "def smooth_data(data, window_size=20):\n",
     "    \"\"\" Apply a simple moving average to smooth the data. \"\"\"\n",
     "    return uniform_filter1d(data, size=window_size, mode='nearest')\n",
     "\n",
@@ -728,7 +722,7 @@
     "x_grid = np.linspace(0, 1, 100)\n",
     "\n",
     "# z-slices corresponding to levels C1 to T2 on the PAM50 template. These will be used to scale the x-label of each subject.\n",
-    "original_vector =np.array([907, 870, 833, 800, 769, 735, 692, 646]) # np.array([984, 938, 907, 870, 833, 800, 769, 735, 692, 646])\n",
+    "original_vector = np.array([984, 938, 907, 870, 833, 800, 769, 735, 692, 646])\n",
     "\n",
     "# Normalize the PAM50 z-slice numbers to the 1-0 range (to show inferior-superior instead of superior-inferior)\n",
     "min_val = original_vector.min()\n",
@@ -753,13 +747,13 @@
     "for i, subject in enumerate(subjects):\n",
     "    ax = axes[i]\n",
     "    \n",
-    "    os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
+    "    os.chdir(os.path.join(path_data, subject, \"fmap\"))\n",
     "\n",
-    "    for shim_method in shim_modes_MPRAGE:\n",
+    "    for shim_method in shim_modes:\n",
     "        # Initialize list to collect data for this shim method\n",
     "        method_data = []\n",
     "\n",
-    "        for file_path in glob.glob(f\"{subject}_acq-{shim_method}_T1w.csv\"):\n",
+    "        for file_path in glob.glob(f\"TB1map_*{shim_method}*.csv\"):\n",
     "            df = pd.read_csv(file_path)\n",
     "            wa_data = df['WA()']\n",
     "\n",
@@ -786,7 +780,7 @@
     "    ax.set_xticklabels([''] * len(original_vector))\n",
     "\n",
     "    ax.set_title(f'{subject}', fontsize=font_size)\n",
-    "    ax.set_ylabel('Signal intensity [a.u.]', fontsize=font_size)\n",
+    "    ax.set_ylabel('B1+ efficiency [nT/V]', fontsize=font_size)\n",
     "    ax.tick_params(axis='y', which='major', labelsize=font_size-4)\n",
     "\n",
     "    # Add legend only to the first subplot\n",
@@ -803,11 +797,11 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "1a690f03",
+   "id": "f8b90a09",
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Create mean and CoV \n",
+    "# Create tables and perform statistics\n",
     "\n",
     "# Go back to root data folder\n",
     "os.chdir(path_data)\n",
@@ -815,28 +809,57 @@
     "# Data storage\n",
     "data_summary = []\n",
     "\n",
-    "# Compute mean and SD of MPRAGE sigint for each subject and each shim mode\n",
+    "# Compute mean and SD of B1+ for each subject and each shim mode\n",
     "for subject in subjects:\n",
-    "    for shim_mode in shim_modes_MPRAGE:\n",
+    "    for shim_mode in shim_modes:\n",
     "        all_data = []\n",
-    "        for file_path in glob.glob(os.path.join(path_data, subject, \"anat\", f\"{subject}_acq-{shim_mode}_T1w.csv\")):\n",
-    "           # print(file_path)\n",
+    "        for file_path in glob.glob(os.path.join(path_data, subject, \"fmap\", f\"TB1map_*{shim_mode}*.csv\")):\n",
     "            df = pd.read_csv(file_path)\n",
+    "            # Drop rows where VertLevel is 1 or 2 (because RF shimming was not performed at these levels)\n",
+    "            df = df[df['VertLevel'].isin(['1', '2']) == False]\n",
     "            wa_data = df['WA()']\n",
     "            all_data.extend(wa_data)\n",
     "\n",
     "        if all_data:\n",
     "            mean_data = np.mean(all_data)\n",
     "            sd_data = np.std(all_data)\n",
-    "            CoV_data=np.round((sd_data/mean_data)*100,2)\n",
-    "            data_summary.append([subject, shim_mode, mean_data, sd_data, CoV_data])\n",
+    "            data_summary.append([subject, shim_mode, mean_data, sd_data])\n",
     "\n",
     "# Convert to DataFrame and save to CSV\n",
-    "df_summary = pd.DataFrame(data_summary, columns=['Subject', 'Shim_Mode', 'Average', 'Standard_Deviation', 'CoV'])\n",
+    "df_summary = pd.DataFrame(data_summary, columns=['Subject', 'Shim_Mode', 'Average', 'Standard_Deviation'])\n",
     "df_summary.to_csv(os.path.join(path_results, 'subject_shim_mode_summary.csv'), index=False)\n",
     "\n",
-    "df_summary_onlyCoV=df_summary[['Subject', 'Shim_Mode', 'CoV']]\n",
-    "print(df_summary_onlyCoV)\n"
+    "# Compute statistics across subjects\n",
+    "grouped_means = df_summary.groupby('Shim_Mode').agg({'Average': ['mean', 'std'], 'Standard_Deviation': ['mean', 'std']})\n",
+    "\n",
+    "# Format mean ± standard deviation\n",
+    "grouped_means['Average_formatted'] = grouped_means['Average']['mean'].map(\"{:.2f}\".format) + \" ± \" + grouped_means['Average']['std'].map(\"{:.2f}\".format)\n",
+    "grouped_means['Standard_Deviation_formatted'] = grouped_means['Standard_Deviation']['mean'].map(\"{:.2f}\".format) + \" ± \" + grouped_means['Standard_Deviation']['std'].map(\"{:.2f}\".format)\n",
+    "\n",
+    "# Drop multi-level index and only keep formatted columns\n",
+    "grouped_means = grouped_means.drop(columns=['Average', 'Standard_Deviation'])\n",
+    "grouped_means.columns = ['Average', 'Standard_Deviation']  # Rename columns for clarity\n",
+    "grouped_means.reset_index().to_csv(os.path.join(path_results, 'average_across_subjects.csv'), index=False)\n",
+    "\n",
+    "# Perform ANOVA and Posthoc Tests\n",
+    "anova_result = f_oneway(*[group[\"Average\"].values for name, group in df_summary.groupby(\"Shim_Mode\")])\n",
+    "print(\"ANOVA Result:\", anova_result)\n",
+    "\n",
+    "if anova_result.pvalue < 0.05:\n",
+    "    posthoc_result = pairwise_tukeyhsd(df_summary['Average'], df_summary['Shim_Mode'])\n",
+    "    print(\"Posthoc Tukey HSD Result:\\n\", posthoc_result)\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a916a6f7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Visualize RF maps to generate Figure 2\n",
+    "\n",
+    "# TODO"
    ]
   }
  ],