 codex/pulse-baseline-covariance-cancellation.py

"""Illustrative correlated-baseline uncertainty; writes PNG in working directory."""
from pathlib import Path
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import numpy as np


def main():
    rho = np.linspace(-1, 1, 401)
    residual_variance = .3**2 + .4**2
    offset_energy = 2.0  # nJ per reading
    variance = residual_variance + 2*offset_energy**2*(1-rho)
    uncertainty = np.sqrt(variance)
    independent = np.sqrt(residual_variance + 2*offset_energy**2)
    assert np.isclose(uncertainty[-1], .5)
    assert np.isclose(independent, np.sqrt(8.25))
    assert np.isclose(np.sqrt(residual_variance+8*(1-.9)), np.sqrt(1.05))
    fig, ax = plt.subplots(figsize=(9, 4), dpi=100, facecolor='white')
    fig.subplots_adjust(left=.10, right=.97, bottom=.20, top=.79)
    ax.plot(rho, uncertainty, color='#126e91', lw=2.5, label='Correct covariance included')
    ax.axhline(independent, color='#b65c26', ls='--', lw=1.8, label='Offsets incorrectly treated as independent')
    ax.scatter([.9, 1], [np.sqrt(1.05), .5], color='#126e91', zorder=3)
    ax.annotate('ρ = 0.90: 1.025 nJ', (.9, np.sqrt(1.05)), xytext=(.14, 1.30),
                arrowprops={'arrowstyle':'->','color':'#126e91'}, fontsize=10)
    ax.annotate('ρ = 1: 0.50 nJ', (1, .5), xytext=(.44, .30),
                arrowprops={'arrowstyle':'->','color':'#126e91'}, fontsize=10)
    ax.set(xlim=(-1.04, 1.04), ylim=(0,4.5), xlabel='Correlation ρ between offsets in the two acquisitions',
           ylabel='Corrected energy uncertainty (nJ)')
    ax.grid(alpha=.22)
    ax.legend(loc='upper right', framealpha=1, fontsize=10)
    fig.suptitle('Baseline subtraction can cancel a shared error', fontsize=14, y=.97)
    fig.text(.5,.85,'Illustrative offset uncertainty: 2 nJ per reading; independent residuals: 0.3 nJ and 0.4 nJ',ha='center',fontsize=10)
    fig.savefig(Path.cwd()/'pulse-baseline-covariance-cancellation.png',dpi=100,facecolor='white')
    plt.close(fig)


if __name__ == '__main__':
    main()