diff --git a/.gitignore b/.gitignore
index 5fb506c..16031b7 100644
--- a/.gitignore
+++ b/.gitignore
@@ -6,5 +6,6 @@ notebooks/unions_shapepipe_extended_2022_W3_v1.0.3.fits
 __pycache__
 __pycache__/*
 shear_psf_leakage/__pycache__/*
+shear_psf_leakage/__pycache__
 */__pycache__/*
 build
diff --git a/citation.cff b/citation.cff
new file mode 100644
index 0000000..8143e34
--- /dev/null
+++ b/citation.cff
@@ -0,0 +1,16 @@
+cff-version: 1.2.0
+message: "If you use this software, please cite it as below."
+authors:
+- family-names: "Kilbinger"
+  given-names: "Martin"
+  orcid: "https://orcid.org/0000-0001-9513-7138"
+- family-names: "Guerrini"
+  given-names: "Sacha"
+  orcid: "https://orcid.org/0009-0004-3655-4870
+- family-names: "Bonini"
+  given-names: "Clara"
+title: "Shear PSF leakage"
+version: 0.1.6
+doi: 00.0000/zenodo.0000
+date-released: 2024-12-12
+url: "https://github.com/CosmoStat/shear_psf_leakage"
diff --git a/notebooks/leakage_autocorr_psf.py b/notebooks/leakage_autocorr_psf.py
deleted file mode 100644
index 6938690..0000000
--- a/notebooks/leakage_autocorr_psf.py
+++ /dev/null
@@ -1,206 +0,0 @@
-# ---
-# jupyter:
-#   jupytext:
-#     formats: ipynb,py:light
-#     text_representation:
-#       extension: .py
-#       format_name: light
-#       format_version: '1.5'
-#       jupytext_version: 1.16.1
-#   kernelspec:
-#     display_name: shear_psf_leakage
-#     language: python
-#     name: shear_psf_leakage
-# ---
-
-# # PSF leakage: PSF auto-correlation contributions
-#
-# Martin Kilbinger <martin.kilbinger@cea.fr>
-
-# %reload_ext autoreload
-# %autoreload 2
-
-# +
-import os
-import matplotlib.pylab as plt
-import numpy as np
-
-from astropy import units
-
-from uncertainties import ufloat
-from uncertainties import unumpy                                                 
-
-from cs_util import canfar
-from cs_util import plots as cs_plots
-
-import shear_psf_leakage.run_object as run_object
-import shear_psf_leakage.run_scale as run_scale
-from shear_psf_leakage import leakage
-# -
-
-# ## Set input parameters
-
-params_in = {}
-
-# ### Paths
-
-# +
-# Input galaxy shear catalogue
-params_in["input_path_shear"] = "unions_shapepipe_extended_2022_v1.0.fits"
-
-# Output directory
-params_in["output_dir"] = "leakage_PSF_autocorr"
-# -
-
-# ### Job control
-
-# Compute (spin-preserving) PSF ellipticity leakage
-params_in["PSF_leakage"] = True
-
-# ### Other parameters
-
-# +
-# PSF ellipticity column names
-params_in["e1_PSF_col"] = "e1_PSF"
-params_in["e2_PSF_col"] = "e2_PSF"
-
-# Set verbose output
-params_in["verbose"] = True
-# -
-
-# ### Retrieve test catalogue from VOspace if not yet downloade
-
-vos_dir = "vos:cfis/XXXX/"
-canfar.download(
-    f"{vos_dir}/{params_in['input_path_shear']}",
-    params_in["input_path_shear"],
-    verbose=params_in["verbose"],
-)
-
-# ## Compute leakage
-
-# Create leakage instance
-obj_object = run_object.LeakageObject()
-
-# Set instance parameters, copy from above
-for key in params_in:
-    obj_object._params[key] = params_in[key]
-
-# Run commands as given in LeakageObject.run()
-
-# +
-# Check parameter validity
-obj_object.check_params()
-
-# Update parameters (here: strings to list)
-obj_object.update_params()
-
-# Prepare output directory
-obj_object.prepare_output()
-# -
-
-# Read input catalogue
-obj_object.read_data()
-
-# Object-by-object spin-consistent PSF leakage
-mix = True
-order = "lin"
-obj_object.PSF_leakage(mix=mix, order=order)
-
-p_dp = {}                                                                    
-for p in obj_object.par_best_fit:                                                         
-    p_dp[p] = ufloat(obj_object.par_best_fit[p].value, obj_object.par_best_fit[p].stderr)
-print(p_dp)
-
-s_ds = leakage.param_order2spin(p_dp, order, mix)
-print(s_ds)
-
-# Get spin elements
-x0 = s_ds["x0"]                                                                     
-x4 = s_ds["x4"]
-y4 = s_ds["y4"]
-y0 = s_ds["y0"]
-
-
-# Create scale-dependent leakage instance
-obj_scale = run_scale.LeakageScale()
-
-# Set instance parameters, copy from above                                   
-for key in params_in:                                                        
-    obj_scale._params[key] = params_in[key]
-
-
-obj_scale._params["input_path_PSF"] = "unions_shapepipe_psf_2022_v1.0.2.fits"
-obj_scale._params["dndz_path"] = "dndz_SP_A.txt"
-
-# ### Run                                                                        
-
-# Check parameter validity                                                       
-obj_scale.check_params()                                                               
-
-# Prepare output directory and stats file                                        
-obj_scale.prepare_output()                                                             
-
-# Prepare input                                                                  
-obj_scale.read_data()                                                                  
-
-# #### Compute correlation function and alpha matrices                           
-# The following command calls `treecorr` to compute auto- and cross-correlation functions.
-# This can take a few minutes.                                                   
-
-obj_scale.compute_corr_gp_pp_alpha_matrix()                                            
-
-obj_scale.alpha_matrix()                                                               
-
-#### For comparison: scalar alpha leakage                                        
-obj_scale.compute_corr_gp_pp_alpha()                                                   
-obj_scale.do_alpha()
-
-
-# Plot terms
-theta_arcmin = obj_scale.get_theta()
-
-term_1 = (
-    (x0 ** 2 + x4 **2 + y0 **2 + y4 ** 2)
-    * (obj_scale.xi_pp_m[0][0] + obj_scale.xi_pp_m[1][1])
-)
-term_2 = 4 * (x0 * y4  - x4 * y0) * obj_scale.xi_pp_m[0][1]
-
-theta = theta_arcmin * units.arcmin
-xi_theo_p, xi_theo_m = run_scale.get_theo_xi(theta, obj_scale._params["dndz_path"])
-
-# +
-y = [
-    unumpy.nominal_values(term_1), 
-    -unumpy.nominal_values(term_2),
-    xi_theo_p,
-]
-dy = [
-    unumpy.std_devs(term_1), 
-    unumpy.std_devs(term_2),
-    [np.nan] * len(xi_theo_p), 
-]
-x = [theta_arcmin] * len(y)
-
-title = "PSF auto-correlation additive terms"
-xlabel = r"$\theta$ [arcmin]" 
-ylabel = "terms"
-out_path = f"{obj_scale._params['output_dir']}/auto_corr_terms.png"                                  
-labels = ["t1" , "-t2", r"$\xi_p$"]
-
-cs_plots.plot_data_1d(
-    x,
-    y,
-    dy,
-    title,
-    xlabel,
-    ylabel,
-    out_path,
-    labels=labels,
-    shift_x=True,
-    xlog=True,
-    ylog=True,
-)
-# -
-
-
diff --git a/notebooks/leakage_scale_spin.py b/notebooks/leakage_scale_spin.py
index ed20567..456354c 100644
--- a/notebooks/leakage_scale_spin.py
+++ b/notebooks/leakage_scale_spin.py
@@ -7,23 +7,22 @@
 #       extension: .py
 #       format_name: light
 #       format_version: '1.5'
-#       jupytext_version: 1.16.1
+#       jupytext_version: 1.15.1
 #   kernelspec:
-#     display_name: shear-psf-leakage-cHl61fDi-py3.9
+#     display_name: sp_validation
 #     language: python
 #     name: python3
 # ---
 
-# # Scale-dependent PSF leakage
+# # PSF contamination: Scale-dependent PSF leakage
 #
-# ## Spin-consistent calculations
-#
-# Computes and plots scale-dependent spin-consistent PSF leakage functions.
+# ## Spin-consistent formalism
 #
 # Martin Kilbinger <martin.kilbinger@cea.fr>
 
 # %reload_ext autoreload
 # %autoreload 2
+# %matplotlib inline
 
 # +
 import os
@@ -35,45 +34,53 @@
 from cs_util import plots as cs_plots
 from cs_util import args
 
-import shear_psf_leakage.run_scale as run
+from shear_psf_leakage import run_scale
 from shear_psf_leakage.leakage import *
 # -
 
-# ## Compute leakage
+# ## Set up
 
-# Create leakage instance
-obj = run.LeakageScale()
+# Create instance of scale-dependent leakage object
+obj_scale = run_scale.LeakageScale()
 
 
 # ### Set input parameters
 
 # Read python parameter file or get user input
-params_upd = args.read_param_script("params_leakage_scale.py", obj._params, verbose=True)
+params_upd = args.read_param_script(
+    "params_leakage_scale.py",
+    obj_scale._params,
+    verbose=True
+)
 for key in params_upd:
-    obj._params[key] = params_upd[key]
+    obj_scale._params[key] = params_upd[key]
 
 # ### Run
 
+# +
 # Check parameter validity
-obj.check_params()
+obj_scale.check_params()
 
 # Prepare output directory and stats file
-obj.prepare_output()
+obj_scale.prepare_output()
+# -
+
+# ### Read input catalogues
 
-# Prepare input
-obj.read_data()
+# Read input galaxy and star catalogue
+obj_scale.read_data()
 
 # #### Compute correlation function and alpha matrices
 # The following command calls `treecorr` to compute auto- and cross-correlation functions.
 # This can take a few minutes.
 
-obj.compute_corr_gp_pp_alpha_matrix()
+obj_scale.compute_corr_gp_pp_alpha_matrix()
 
-obj.alpha_matrix()
+obj_scale.alpha_matrix()
 
 #### For comparison: scalar alpha leakage
-obj.compute_corr_gp_pp_alpha()
-obj.do_alpha()
+obj_scale.compute_corr_gp_pp_alpha()
+obj_scale.do_alpha()
 
 # #### PSF auto-correlation function correlation matrix
 #
@@ -85,25 +92,26 @@
 # $$
 
 # Check symmetry of PSF auto-correlation matrix
-diff = obj.Xi_pp_m[0][1] - obj.Xi_pp_m[1][0]
+diff = obj_scale.Xi_pp_m[0][1] - obj_scale.Xi_pp_m[1][0]
 print(
     "Is r symmetrical? max abs (rel) diff ="
-    +f" {max(np.abs(diff)):.3e} ({max(np.abs(diff / obj.Xi_pp_m[0][1])):.3e})",
+    + f" {max(np.abs(diff)):.3e}"
+    + f" ({max(np.abs(diff / obj_scale.Xi_pp_m[0][1])):.3e})",
 )
 
 # +
 # Plot r and ratios of r
 
-theta = obj.get_theta()
+theta = obj_scale.get_theta()
 
-# Exact calculation
+# Exact: Using centered correlation functions
 r = []
 r_ratio_1 = []
 r_ratio_2 = []
 for ndx in range(len(theta)):
     my_r = (
-        obj.Xi_pp_ufloat[ndx][0, 1] ** 2
-        / (obj.Xi_pp_ufloat[ndx][0, 0] * obj.Xi_pp_ufloat[ndx][1, 1])
+        obj_scale.Xi_pp_ufloat[ndx][0, 1] ** 2
+        / (obj_scale.Xi_pp_ufloat[ndx][0, 0] * obj_scale.Xi_pp_ufloat[ndx][1, 1])
     )
     r.append(my_r)
     r_ratio_1.append(1 / (1 - my_r))
@@ -111,8 +119,8 @@
 
 print("min max mean r = ", np.min(r), np.max(r), np.mean(r))
 
-# Approximate
-r_fast = obj.xi_pp_m[0][1] ** 2 / (obj.xi_pp_m[0][0] * obj.xi_pp_m[1][1])
+# Approximate: Using uncentered correlation functions
+r_fast = obj_scale.xi_pp_m[0][1] ** 2 / (obj_scale.xi_pp_m[0][0] * obj_scale.xi_pp_m[1][1])
 
 n = 6
 theta_arr = [theta] * n
@@ -140,14 +148,14 @@
 xlabel = r"$\theta$ [arcmin]"
 ylabel = r"functions of $r(\theta)$"
 
-fac = 1.1
+fac = 0.9
 xlim = [
-    obj._params["theta_min_amin"] / fac,
-    obj._params["theta_max_amin"] * fac
+    obj_scale._params["theta_min_amin"] * fac,
+    obj_scale._params["theta_max_amin"],
 ]
 ylim = (-0.5, 2)
 
-out_path = f"{obj._params['output_dir']}/r.png"
+out_path = f"{obj_scale._params['output_dir']}/r.png"
 
 title = ""
 
@@ -173,19 +181,17 @@
 # +
 # Plot alpha matrix elements
 
-theta = obj.get_theta()
+ylim = obj_scale._params["leakage_alpha_ylim"]
 
-xlim = [obj._params["theta_min_amin"] / fac, obj._params["theta_max_amin"] * fac]
-ylim = obj._params["leakage_alpha_ylim"]
-
-theta_arr = [theta] * 4
+n = 4
+theta_arr = [theta] * n
     
 alpha = []
 yerr = []
 labels = []
 for idx in (0, 1):
     for jdx in (0, 1):   
-        alpha_ufloat = obj.get_alpha_ufloat(idx, jdx)
+        alpha_ufloat = obj_scale.get_alpha_ufloat(idx, jdx)
         alpha.append(unumpy.nominal_values(alpha_ufloat))
         yerr.append(unumpy.std_devs(alpha_ufloat))
         labels.append(rf"$\alpha_{{{idx+1}{jdx+1}}}$")
@@ -197,7 +203,7 @@
 xlabel = r"$\theta$ [arcmin]"
 ylabel = r"$\alpha_{ij}(\theta)$"
 title = ""
-out_path = f"{obj._params['output_dir']}/alpha_ij.png"
+out_path = f"{obj_scale._params['output_dir']}/alpha_ij.png"
 
 cs_plots.plot_data_1d(
     theta_arr,
@@ -219,49 +225,39 @@
 )
 
 # +
-# Plot spin elements: Spin-0 and spin-4
-
-x0 = 0.5 * (
-     obj.get_alpha_ufloat(0, 0)
-    + obj.get_alpha_ufloat(1, 1)
-)
-x4 = 0.5 * (
-    obj.get_alpha_ufloat(0, 0)
-    - obj.get_alpha_ufloat(1, 1)
-)
-y4 = 0.5 * (obj.get_alpha_ufloat(0, 1) + obj.get_alpha_ufloat(1, 0))
-y0 = 0.5 * (obj.get_alpha_ufloat(1, 0) - obj.get_alpha_ufloat(0, 1))
+# Plot spin coefficients
 
+obj_scale.compute_alpha_spin_coeffs()
 
-ylim = obj._params["leakage_alpha_ylim"]
+n = 4
+theta_arr = [theta] * n
 
 y_arr = [
-    unumpy.nominal_values(x0),
-    unumpy.nominal_values(x4),
-    unumpy.nominal_values(y4),
-    unumpy.nominal_values(y0),
+    unumpy.nominal_values(obj_scale._alpha_0_r),
+    unumpy.nominal_values(obj_scale._alpha_0_i),
+    unumpy.nominal_values(obj_scale._alpha_4_r),
+    unumpy.nominal_values(obj_scale._alpha_4_i),
 ]
 dy_arr = [
-    unumpy.std_devs(x0),
-    unumpy.std_devs(x4),
-    unumpy.std_devs(y4),
-    unumpy.std_devs(y0),
+    unumpy.std_devs(obj_scale._alpha_0_r),
+    unumpy.std_devs(obj_scale._alpha_0_i),
+    unumpy.std_devs(obj_scale._alpha_4_r),
+    unumpy.std_devs(obj_scale._alpha_4_i),
 ]
-theta_arr = [obj.get_theta()] * len(y_arr)
 labels = [
-    r"$\alpha^\Re_0 = (\alpha_{11} + \alpha_{22})/2$",
-    r"$\alpha^\Re_4 = (\alpha_{11} - \alpha_{22})/2$",
-    r"$\alpha^\Im_4 = (\alpha_{12} + \alpha_{21})/2$",
-    r"$\alpha^\Im_0 = (-\alpha_{12} + \alpha_{21})/2$"
+    r"$\alpha^\Re_0$",
+    r"$\alpha^\Im_0$",
+    r"$\alpha^\Re_r$",
+    r"$\alpha^\Im_r$"
 ]
 colors = ["blue", "orange", "green", "magenta"]
-markers = ["o", "s", "h", "^"]
+markers = ["o", "s", "^", "v"]
 linestyles = ["-"] * 4 
 
 xlabel = r"$\theta$ [arcmin]"
 ylabel = r"Components of leakage matrix"
 title = ""
-out_path = f"{obj._params['output_dir']}/alpha_leakage_m_s0_s4.png"
+out_path = f"{obj_scale._params['output_dir']}/alpha_leakage_m_s0_s4.png"
 
 cs_plots.plot_data_1d(
     theta_arr,
@@ -283,13 +279,13 @@
 
 # Including scalar leakage for comparison
 theta_arr.append(theta)
-y_arr.append(obj.alpha_leak)
-dy_arr.append(obj.sig_alpha_leak)
-labels.append(r"$\alpha$ (scalar approx.)")
-colors.append("blue")
-markers.append("p")
+y_arr.append(obj_scale.alpha_leak)
+dy_arr.append(obj_scale.sig_alpha_leak)
+labels.append(r"$\alpha^{\rm s}_{+}$")
+colors.append("cyan")
+markers.append("x")
 linestyles.append("--")
-out_path = f"{obj._params['output_dir']}/alpha_leakage_m_s0_s4_as.png"
+out_path = f"{obj_scale._params['output_dir']}/alpha_leakage_m_s0_s4_as.png"
 
 cs_plots.plot_data_1d(
     theta_arr,
@@ -309,47 +305,49 @@
     close_fig=False,
     shift_x=True,
 )
-
+# + \xi_ {"incorrectly_encoded_metadata": "{22}^\\textrm{g,p} \\, \\xi_{11}^\\textrm{p,p}"}
 # #### Consistency relations for scalar leakage
 
 # If the leakage is a scalar function, it can be expressed in three different ways.
 
 # +
-xlim = [
-    obj._params["theta_min_amin"] / fac,
-    obj._params["theta_max_amin"] * fac
-]
-ylim = obj._params["leakage_alpha_ylim"]
-
 alpha_1 = []
 alpha_2 = []
 for ndx in range(len(theta)):
-    my_a1 = obj.Xi_gp_ufloat[ndx][0, 0] / obj.Xi_pp_ufloat[ndx][0, 0]
+    my_a1 = obj_scale.Xi_gp_ufloat[ndx][0, 0] / obj_scale.Xi_pp_ufloat[ndx][0, 0]
     alpha_1.append(my_a1)
-    my_a2 = obj.Xi_gp_ufloat[ndx][1, 1] / obj.Xi_pp_ufloat[ndx][1, 1]
+    my_a2 = obj_scale.Xi_gp_ufloat[ndx][1, 1] / obj_scale.Xi_pp_ufloat[ndx][1, 1]
     alpha_2.append(my_a2)
 
+# TODO: Use centered functions for all cases
+
 y = [
-    unumpy.nominal_values(x0),
-    obj.alpha_leak,
+    obj_scale.alpha_leak,
     unumpy.nominal_values(alpha_1),
     unumpy.nominal_values(alpha_2),
+    unumpy.nominal_values(obj_scale._alpha_0_r),
 ]
 dy = [
-    unumpy.std_devs(x0),
-    obj.sig_alpha_leak,
+    obj_scale.sig_alpha_leak,
     unumpy.std_devs(alpha_1),
     unumpy.std_devs(alpha_2),
+    unumpy.std_devs(obj_scale._alpha_0_r),
 ]
-theta_arr = [obj.get_theta()] * len(y)
+theta_arr = [theta] * len(y)
 
-
-labels = [r"$\alpha_0^\Re$", r"$\alpha_+$", r"$\alpha_1$", r"$\alpha_2$"]
+labels = [
+    r"$\alpha^{\rm s}_+$",
+    r"$\alpha^{\rm s}_1$",
+    r"$\alpha^{\rm s}_2$",
+    r"$\alpha^\Re_0$",
+]
 xlabel = r"$\theta$ [arcmin]"
 ylabel = r"$\alpha(\theta)$"
-linestyles = ["-", "--", "--", "--"] 
 title = ""
-out_path = f"{obj._params['output_dir']}/alpha_leakage_scalar_consistency.png"
+out_path = f"{obj_scale._params['output_dir']}/alpha_leakage_scalar_consistency.png"
+colors = ["cyan", "black", "grey", "blue"]
+markers = ["x", "h", "o", "o"]
+linestyles = ["--", "--", "--", "-"]
 
 cs_plots.plot_data_1d(
     theta_arr,
@@ -364,25 +362,23 @@
     xlim=xlim,
     ylim=ylim,
     close_fig=False,
-    linestyles=linestyles,
     shift_x=True,
+    colors=colors,
+    markers=markers,
+    linestyles=linestyles,
 )
 # -
 
 # If alpha is a scalar, the mixed-component centered cross-correlation functions should be identical.
 
 # +
-xlim = [obj._params["theta_min_amin"], obj._params["theta_max_amin"]]
-
-theta_arr = [theta] * 3
-
 Xi_12 = []
 Xi_21 = []
 Xi_tr = []
 for ndx in range(len(theta)):
-    Xi_12.append(obj.Xi_gp_ufloat[ndx][0, 1])
-    Xi_21.append(obj.Xi_gp_ufloat[ndx][1, 0])
-    Xi_tr.append(obj.Xi_gp_ufloat[ndx][0, 0] + obj.Xi_gp_ufloat[ndx][0, 0])
+    Xi_12.append(obj_scale.Xi_gp_ufloat[ndx][0, 1])
+    Xi_21.append(obj_scale.Xi_gp_ufloat[ndx][1, 0])
+    Xi_tr.append(obj_scale.Xi_gp_ufloat[ndx][0, 0] + obj_scale.Xi_gp_ufloat[ndx][0, 0])
 
 y = [
     unumpy.nominal_values(Xi_12),
@@ -394,12 +390,15 @@
     unumpy.std_devs(Xi_21),
     unumpy.std_devs(Xi_tr),
 ]
+theta_arr = [theta] * len(y)
 
-labels = [r"$\tau_{12,0}$", r"$\tau_{21,0}$", r"tr$\tau_0$"]
+labels = [r"$\tau_{0, 12}$", r"$\tau_{0,21}$", r"tr$\tau_{0}$"]
 xlabel = r"$\theta$ [arcmin]"
-ylabel = r"correlation functions $\tau_{ij,0}$"
+ylabel = r"correlation functions $\tau_{0,ij}$"
 title = ""
-out_path = f"{obj._params['output_dir']}/Xi_mixed_consistency.png"
+out_path = f"{obj_scale._params['output_dir']}/Xi_mixed_consistency.png"
+markers = ["o", "s", "d"]
+linestyles = ["-", "--", ":"]
 
 cs_plots.plot_data_1d(
     theta_arr,
@@ -414,22 +413,20 @@
     xlim=xlim,
     close_fig=False,
     shift_x=True,
+    markers=markers,
+    linestyles=linestyles,
 )
 
 # +
 # For comparison, plot the same for the PSF - PSF correlations
 
-xlim = [obj._params["theta_min_amin"], obj._params["theta_max_amin"]]
-
-theta_arr = [theta] * 3
-
 Xi_12 = []
 Xi_21 = []
 Xi_tr = []
 for ndx in range(len(theta)):
-    Xi_12.append(obj.Xi_pp_ufloat[ndx][0, 1])
-    Xi_21.append(obj.Xi_pp_ufloat[ndx][1, 0])
-    Xi_tr.append(obj.Xi_pp_ufloat[ndx][0, 0] + obj.Xi_pp_ufloat[ndx][0, 0])
+    Xi_12.append(obj_scale.Xi_pp_ufloat[ndx][0, 1])
+    Xi_21.append(obj_scale.Xi_pp_ufloat[ndx][1, 0])
+    Xi_tr.append(obj_scale.Xi_pp_ufloat[ndx][0, 0] + obj_scale.Xi_pp_ufloat[ndx][0, 0])
 
 y = [
     unumpy.nominal_values(Xi_12),
@@ -441,13 +438,15 @@
     unumpy.std_devs(Xi_21),
     unumpy.std_devs(Xi_tr),
 ]
+theta_arr = [theta] * len(y)
 
-
-labels = [r"$\rho_{12,0}$", r"$\rho_{21,0}$", r"tr$\rho_0$"]
+labels = [r"$\rho_{0,12}$", r"$\rho_{0,21}$", r"tr$\rho_0$"]
 xlabel = r"$\theta$ [arcmin]"
 ylabel = r"Centered correlation functions"
 title = ""
-out_path = f"{obj._params['output_dir']}/Xi_pp_mixed_consistency.png"
+out_path = f"{obj_scale._params['output_dir']}/Xi_pp_mixed_consistency.png"
+markers = ["o", "s", "d"]
+linestyles = ["-", "--", ":"]
 
 cs_plots.plot_data_1d(
     theta_arr,
@@ -462,4 +461,9 @@
     xlim=xlim,
     close_fig=False,
     shift_x=True,
+    markers=markers,
+    linestyles=linestyles,
 )
+# -
+
+
diff --git a/notebooks/leakage_xi_sys.py b/notebooks/leakage_xi_sys.py
new file mode 100644
index 0000000..1e25819
--- /dev/null
+++ b/notebooks/leakage_xi_sys.py
@@ -0,0 +1,455 @@
+# ---
+# jupyter:
+#   jupytext:
+#     formats: ipynb,py:light
+#     text_representation:
+#       extension: .py
+#       format_name: light
+#       format_version: '1.5'
+#       jupytext_version: 1.15.1
+#   kernelspec:
+#     display_name: sp_validation
+#     language: python
+#     name: python3
+# ---
+
+# # PSF contamination: Star-galaxy correlation function $\xi_{\rm{sys}}$
+#
+# Martin Kilbinger <martin.kilbinger@cea.fr>
+
+# The PSF contamination diagnostic function $\xi_{\rm{sys}} was introduced in Bacon et al. (2003) in the scalar PSF approximation.
+# It is defined as
+# \begin{equation}
+#     \xi_{\textrm{sys}, +}(\theta)
+#         = \frac{\left( \xi_+^\textrm{e, p} \right)^2(\theta)}
+#                 {\xi_+^\textrm{p, p}(\theta)}
+#         = \frac{\tau_{+, 0}^2(\theta)}{\rho_{0, +}(\theta)}
+#         = \alpha^2(\theta) \rho_{0, +}(\theta).
+# \end{equation}
+#
+# The generalisation of this function to spin-consistent notation is the matrix
+# \begin{equation}
+#     \xi_{\textrm{sys}, ij}(\theta)
+#     = \sum_{k,l=1}^2 \alpha_{ik} \alpha_{jl} \, \rho_{0, kl}(\theta).
+# \end{equation}
+# The ‘+’ components can be computed as
+# \begin{align}
+#     \xi_{\textrm{sys}, +}(\theta)
+#     = & \sum_{k,l=1}^2 \sum_{i=1}^2 \alpha_{ik} \alpha_{il} \rho_{0, kl}(\theta)
+#     \nonumber \\
+#     = & \,
+#         \left(
+#             \alpha_{11}^2 + \alpha_{21}^2
+#         \right) \rho_{0, 11}(\theta)
+#         + \left(
+#             \alpha_{12}^2 + \alpha_{22}
+#         \right)^2 \rho_{0, 22}(\theta)
+#         \nonumber \\
+#         & \,
+#         + 2 \left(
+#             \alpha_{11} \alpha_{12} + \alpha_{21} \alpha_{22}
+#         \right) \rho_{0, 12}(\theta)
+#         \nonumber \\
+#     = & \,
+#         \left| \alpha_0 + \alpha_4 \right|^2 \rho_{0, 11}(\theta)
+#         + \left| \alpha_0 - \alpha_4 \right|^2 \rho_{0, 22}(\theta)
+#         \nonumber \\
+#         & \,
+#         + 4 \left( \alpha_0^\Re \alpha_4^\Im + \alpha_4^\Re \alpha_0^\Im \right)
+#         \rho_{0, 12}(\theta)
+#         \nonumber \\
+#     = & \,
+#         2 \left[
+#             \left(\alpha_0^\Re \right)^2 + \left(\alpha_4^\Im\right)^2
+#             + \left(\alpha_0^\Re \right)^2 + \left(\alpha_4^\Im\right)^2
+#         \right]
+#         \rho_{0, +}(\theta)
+#         \nonumber \\
+#       &
+#       + 2 \left( \alpha_0^\Re \alpha_4^\Re + \alpha_0^\Im \alpha_4^\Im \right)
+#         \left[ \rho_{0, 11}(\theta) - \rho_{0, 22}(\theta) \right]
+#        \nonumber \\
+#       &
+#       + 4 \left( \alpha_0^\Re \alpha_4^\Im - \alpha_4^\Re \alpha_0^\Im \right)
+#         \rho_{0, 12}(\theta)
+#         ,
+# \end{align}
+
+# %reload_ext autoreload
+# %autoreload 2
+# %matplotlib inline
+
+# +
+import os
+import matplotlib.pylab as plt
+import numpy as np
+
+from astropy import units
+
+from uncertainties import ufloat
+from uncertainties import unumpy                                                 
+
+from cs_util import canfar
+from cs_util import plots as cs_plots
+from cs_util import args
+
+from shear_psf_leakage import run_object
+from shear_psf_leakage import run_scale
+from shear_psf_leakage import leakage
+# -
+
+# ## Set up
+
+# Create instance of scale-dependent leakage object
+obj_scale = run_scale.LeakageScale()
+
+# Read python parameter file or get user input
+params_upd = args.read_param_script("params_xi_sys.py", obj_scale._params, verbose=True)
+for key in params_upd:
+    obj_scale._params[key] = params_upd[key]
+
+
+# +
+# Check parameter validity
+obj_scale.check_params()
+
+# Prepare output directory
+obj_scale.prepare_output()
+# -
+
+# ### Read input catalogues
+
+# Read input galaxy and star catalogue
+obj_scale.read_data()
+
+# #### Compute correlation function and alpha matrices                           
+# The following command calls `treecorr` to compute auto- and cross-correlation functions.
+# This can take a few minutes.                                                   
+
+obj_scale.compute_corr_gp_pp_alpha_matrix()                                            
+
+obj_scale.alpha_matrix()                                                               
+
+# +
+# Compute spin coefficients from matrix elements
+
+alpha_0_r = (
+    0.5 * (
+        obj_scale.get_alpha_ufloat(0, 0)
+        + obj_scale.get_alpha_ufloat(1, 1)
+    )
+)
+alpha_0_i = (
+    0.5 * (
+        -obj_scale.get_alpha_ufloat(0, 1)
+        + obj_scale.get_alpha_ufloat(1, 0)
+    )
+)
+alpha_4_r = (
+    0.5 * (
+        obj_scale.get_alpha_ufloat(0, 0)
+        - obj_scale.get_alpha_ufloat(1, 1)
+    )
+)
+alpha_4_i = (
+    0.5 * (
+        obj_scale.get_alpha_ufloat(0, 1)
+        + obj_scale.get_alpha_ufloat(1, 0)
+    )
+)
+
+
+# -
+
+def get_rho_0(self, idx, jdx):                                        
+	"""Get Rho 0.                                                     
+
+	Return (i, j) matrix element of rho_0.                
+
+	Parameters                                                               
+	----------                                                               
+	idx : int                                                                
+		line index, allowed are 0 or 1                                       
+	jdx : int                                                                
+		column index, allowed are 0 or 1                                     
+
+	Returns                                                                  
+	-------                                                                  
+	numpy.ndarray                                                            
+		matrix element as array over scales, each entry is                   
+		of type ufloat                                                       
+
+	"""  
+	mat = []                                                                 
+	n_theta = self._params["n_theta"]                                        
+	for ndx in range(n_theta):                                               
+		mat.append(self.Xi_pp_ufloat[ndx][idx, jdx])
+	return np.array(mat)
+
+
+# +
+# xi_sys terms
+
+
+# Split into +, '-', and mixed parts
+
+xi_sys_term_p = (
+    (alpha_0_r ** 2 + alpha_0_i ** 2 + alpha_4_r ** 2 + alpha_4_i ** 2)
+    * (get_rho_0(obj_scale, 0, 0) + get_rho_0(obj_scale, 1, 1))
+)
+xi_sys_term_m = (
+    2 * (alpha_0_r * alpha_4_r + alpha_0_i * alpha_4_i)
+     * (get_rho_0(obj_scale, 0, 0) - get_rho_0(obj_scale, 1, 1))
+)
+xi_sys_term_mixed = (
+    4 * (
+        alpha_0_r * alpha_4_i - alpha_4_r * alpha_0_i
+    ) * get_rho_0(obj_scale, 0, 1)
+)
+
+xi_sys_tot = xi_sys_term_p + xi_sys_term_m + xi_sys_term_mixed
+# -
+
+
+#### For comparison: scalar xi sys.
+obj_scale.compute_corr_gp_pp_alpha()
+obj_scale.do_alpha()
+obj_scale.compute_xi_sys()
+
+
+# ## Plotting
+
+# angular scales in arcmin
+theta_arcmin = obj_scale.get_theta()
+
+# +
+# Mean ellipticities for centered correlation functions
+
+e1_gal = obj_scale.dat_shear["e1"]
+e2_gal = obj_scale.dat_shear["e2"]
+weights_gal = obj_scale.dat_shear["w"]
+
+complex_gal = (                                                         
+    np.average(e1_gal, weights=weights_gal)                             
+    + np.average(e2_gal, weights=weights_gal) * 1j
+)
+
+e1_star = obj_scale.dat_PSF["E1_STAR_HSM"]
+e2_star = obj_scale.dat_PSF["E2_STAR_HSM"]
+complex_psf = np.mean(e1_star) + np.mean(e2_star) * 1j
+
+# +
+# At the moment we compute the scalar xi_sys here using centered correlation functions.
+# TODO: Implement consistent centering in shear_psf_leakage classes. 
+
+# xi_sys = tau_0^2/rho_0
+xi_sys_scalar = (obj_scale.r_corr_gp.xip - np.real(np.conj(complex_gal) * complex_psf)) ** 2 / (obj_scale.r_corr_pp.xip - np.abs((complex_psf) ** 2))
+std_xi_sys_scalar = obj_scale.C_sys_std_p
+
+y = [
+    unumpy.nominal_values(xi_sys_term_p),
+    unumpy.nominal_values(xi_sys_term_m),
+    unumpy.nominal_values(xi_sys_term_mixed),
+    unumpy.nominal_values(xi_sys_tot),
+    xi_sys_scalar,
+    
+]
+dy = [
+    unumpy.std_devs(xi_sys_term_p),
+    unumpy.std_devs(xi_sys_term_m),
+    unumpy.std_devs(xi_sys_term_mixed),
+    unumpy.std_devs(xi_sys_tot),
+    std_xi_sys_scalar,
+]
+x = [theta_arcmin] * len(y)
+
+title = r"Bacon et al. (2003) $\xi_{sys}$"
+xlabel = r"$\theta$ [arcmin]" 
+ylabel = "terms"
+out_path = f"{obj_scale._params['output_dir']}/xi_sys_terms.png"                                  
+labels = ["$t_+$", "$t_-$", r"$t_{\rm mixed}$", r"$\sum t_i$", "scalar"]
+markers = ["o", "d", "^", "x", "s"]
+
+factor = 0.9
+xlim = [obj_scale._params["theta_min_amin"] ** factor, obj_scale._params["theta_max_amin"]]
+
+#ylim = [-1e-6, 2e-6]
+
+cs_plots.plot_data_1d(
+    x,
+    y,
+    dy,
+    title,
+    xlabel,
+    ylabel,
+    out_path,
+    labels=labels,
+    shift_x=True,
+    xlog=True,
+    ylog=False,
+    xlim=xlim,
+    markers=markers,
+)
+# +
+# Consistency test. Different computations of xi_sys terms
+
+# Using a_ij matrix
+xi_sys_term_m_11 = (
+    (obj_scale.get_alpha_ufloat(0, 0) ** 2 + obj_scale.get_alpha_ufloat(1, 0) ** 2)
+    * get_rho_0(obj_scale, 0, 0)
+)
+xi_sys_term_m_22 = (
+    (obj_scale.get_alpha_ufloat(0, 1) ** 2 + obj_scale.get_alpha_ufloat(1, 1) ** 2)
+    * get_rho_0(obj_scale, 1, 1)
+)
+xi_sys_term_m_12 = (
+    2 * (obj_scale.get_alpha_ufloat(0, 0) * obj_scale.get_alpha_ufloat(0, 1)
+     + obj_scale.get_alpha_ufloat(1, 0) * obj_scale.get_alpha_ufloat(1, 1))
+     * get_rho_0(obj_scale, 0, 1)
+)
+
+# Sum of the three terms
+xi_sys_m_tot =  xi_sys_term_m_11 + xi_sys_term_m_22 + xi_sys_term_m_12
+
+# Matrix sum
+xi_sys_m2_tot = 0
+for idx in (0, 1):
+    for kdx in (0, 1):
+        for ldx in (0, 1):
+            xi_sys_m2_tot +=  (
+                obj_scale.get_alpha_ufloat(idx, kdx)
+                * obj_scale.get_alpha_ufloat(idx, ldx)
+                * get_rho_0(obj_scale, kdx, ldx)
+            )
+
+y = [
+    unumpy.nominal_values(xi_sys_term_m_11),
+    unumpy.nominal_values(xi_sys_term_m_22),
+    unumpy.nominal_values(xi_sys_term_m_12),
+    unumpy.nominal_values(xi_sys_m_tot),
+    unumpy.nominal_values(xi_sys_m2_tot),
+    unumpy.nominal_values(xi_sys_tot),
+    
+]
+dy = [
+    unumpy.std_devs(xi_sys_term_m_11),
+    unumpy.std_devs(xi_sys_term_m_22),
+    unumpy.std_devs(xi_sys_term_m_12),
+    unumpy.std_devs(xi_sys_m_tot),
+    unumpy.std_devs(xi_sys_m2_tot),
+    unumpy.std_devs(xi_sys_tot),
+]
+x = [theta_arcmin] * len(y)
+
+title = r"Bacon et al. (2003) $\xi_{sys}$"
+xlabel = r"$\theta$ [arcmin]" 
+ylabel = "terms with matrix coeffs"
+out_path = f"{obj_scale._params['output_dir']}/xi_sys_terms_m.png"                                  
+labels = ["$t_{11}$", "$t_{22}$", r"$t_{12}$", r"$\sum t_i$ (matrix terms)", r"$\sum t_i$ (matrix terms 2)", "$\sum t_i$ (spin terms)"]
+markers = ["o", "d", "^", "x", "h", "s"]
+linestyles = [":"] * 3
+linestyles.extend(["-", "-.", "--"])
+
+ylim = [-1e-6, 2e-6]
+
+cs_plots.plot_data_1d(
+    x,
+    y,
+    dy,
+    title,
+    xlabel,
+    ylabel,
+    out_path,
+    labels=labels,
+    shift_x=False,
+    xlog=True,
+    ylog=False,
+    xlim=xlim,
+    ylim=ylim,
+    markers=markers,
+    linestyles=linestyles,
+)
+
+# +
+# Plot difference, should be consistent with zero
+
+d1 = xi_sys_m_tot - xi_sys_tot
+d2 = xi_sys_m2_tot - xi_sys_tot
+d3 = xi_sys_m_tot - xi_sys_m2_tot
+
+d4 = xi_sys_term_m_11 + xi_sys_term_m_22 - (xi_sys_term_p + xi_sys_term_m)
+d5 = xi_sys_term_m_12 - xi_sys_term_mixed
+
+y = [
+    unumpy.nominal_values(d1),
+    unumpy.nominal_values(d2),
+    unumpy.nominal_values(d3),
+    unumpy.nominal_values(d4),
+    unumpy.nominal_values(d5),
+]
+dy = [
+    unumpy.std_devs(d1),
+    unumpy.std_devs(d2),
+    unumpy.std_devs(d3),
+    #unumpy.std_devs(d4),
+    #unumpy.std_devs(d5),
+]
+x = [theta_arcmin] * len(y)
+
+title = r"Bacon et al. (2003) $\xi_{sys}$"
+xlabel = r"$\theta$ [arcmin]" 
+ylabel = "difference"
+out_path = f"{obj_scale._params['output_dir']}/xi_sys_diff.png"                                  
+labels = [r"(matrix - spin) tot", r"(matrix 2 - spin) tot", "(matrix - matrix 2) tot", "11 22", "mixed"]
+markers = ["o", "d", "s", "x", "v"]
+linestyles = ["-", "--", ":", ":", ":"]
+
+
+cs_plots.plot_data_1d(
+    x[:3],
+    y[:3],
+    dy[:3],
+    title,
+    xlabel,
+    ylabel,
+    out_path,
+    labels=labels,
+    shift_x=True,
+    xlog=True,
+    ylog=False,
+    xlim=xlim,
+    markers=markers,
+    linestyles=linestyles,
+)
+
+# +
+# Test: compare different estimates of xi_sys. Some use uncentered, some centered correlation functions.
+
+plt.clf()
+plt.loglog(theta_arcmin, obj_scale.C_sys_p, "p:", label="scalar uncentered")
+
+# spin t_+
+plt.loglog(
+    theta_arcmin,
+    unumpy.nominal_values(xi_sys_term_p), 
+    "d-",
+    label="spin $t_+$ centered",
+)
+
+# tau^2 / rho
+plt.loglog(theta_arcmin, obj_scale.r_corr_gp.xip ** 2 / obj_scale.r_corr_pp.xip, "v:", label="scalar centered")
+plt.loglog(
+    theta_arcmin,
+    (obj_scale.r_corr_gp.xip - np.real(np.conj(complex_gal) * complex_psf)) ** 2 / (obj_scale.r_corr_pp.xip - np.abs((complex_psf) ** 2)),
+    "v-",
+    label="3c"
+)
+
+plt.legend()
+plt.show()
+plt.savefig(f"{obj_scale._params['output_dir']}/xi_sys_test.png")
+
+# -
+
+
+
diff --git a/shear_psf_leakage/run_scale.py b/shear_psf_leakage/run_scale.py
index d5783e6..ca0f5c0 100644
--- a/shear_psf_leakage/run_scale.py
+++ b/shear_psf_leakage/run_scale.py
@@ -771,15 +771,17 @@ def plot_alpha_leakage(self, close_fig=True, xlinlog="log"):
 
         """
         # Set some x-values and limits for plot
+        factor = 0.9
+        x0 = self._params["theta_min_amin"]
         if xlinlog == "log":
-            x0 = self._params["theta_min_amin"]
+            x0 = x0 * factor
             x_affine = np.geomspace(x0, self._params["theta_max_amin"])
             xlim = [x0, self._params["theta_max_amin"]]
             xlog = True
         else:
-            x0 = -10
-            x_affine = np.linspace(0, self._params["theta_max_amin"])
-            xlim = [x0 * 2, self._params["theta_max_amin"]]
+            x0 = x0 ** factor
+            x_affine = np.linspace(x0, self._params["theta_max_amin"])
+            xlim = [x0, self._params["theta_max_amin"]]
             xlog = False
 
         # alpha(theta) data points
@@ -1155,6 +1157,38 @@ def get_alpha_ufloat(self, idx, jdx):
             mat.append(self.alpha_leak_ufloat[ndx][idx, jdx])
 
         return np.array(mat)
+    
+    def compute_alpha_spin_coeffs(self):
+        """Compute Alpha Spin Coefficients.
+        
+        Compute the spin coefficients of the PSF leakage alpha(theta) from
+        the matrix elements.
+        
+        """
+        self._alpha_0_r = (
+            0.5 * (
+                self.get_alpha_ufloat(0, 0)
+                + self.get_alpha_ufloat(1, 1)
+            )
+        )
+        self._alpha_0_i = (
+            0.5 * (
+                -self.get_alpha_ufloat(0, 1)
+                + self.get_alpha_ufloat(1, 0)
+            )
+        )
+        self._alpha_4_r = (
+            0.5 * (
+                self.get_alpha_ufloat(0, 0)
+                - self.get_alpha_ufloat(1, 1)
+            )
+        )
+        self._alpha_4_i = (
+            0.5 * (
+                self.get_alpha_ufloat(0, 1)
+                + self.get_alpha_ufloat(1, 0)
+            )
+        )
 
     def do_alpha_matrix(self):
         """Do Alpha Matrix.