#  draw_medsig2.py -- reproduces Fig. 7 in
#  Cowan, Cranmer, Gross and Vitells, EPJC 71 (2011) 1554; arXiv:1007.1727.
#  G. Cowan, March 2025

import numpy as np
import matplotlib.pyplot as plt

# Load the data into a 2D NumPy array.
data = np.loadtxt('medsig.txt')

# Extract values
x = data[:, 0]
y_pois = data[:, [2, 3, 4]]
y_asimov = data[:, [6, 7, 8]]
y_exact = data[:, [10, 11, 12]]

# Create a square figure (e.g., 6x6 inches).
fig = plt.figure(figsize=(6, 6))
ax = fig.gca()
ax.tick_params(axis="x", direction="in")
ax.tick_params(axis="y", direction="in")
ax.tick_params(which="minor", axis="x", direction="in")
ax.tick_params(which="minor", axis="y", direction="in")
plt.tick_params(axis="x", which="major", pad=10)
plt.tick_params(axis="y", which="major", pad=10)

# Asimov formula (blue curves)
for i in range(y_asimov.shape[1]):
    plt.plot(x, y_asimov[:, i], 'b-', linewidth=1.5,
             label='$\\sqrt{q_{0,\\rm A}}$' if i == 0 else None)

# s/sqrt(b) (red dashed curves)
for i in range(y_pois.shape[1]):
    plt.plot(x, y_pois[:, i], 'r--', linewidth=1.5, 
             label='$s / \\sqrt{b}$' if i == 0 else None)

# "Exact" values (from Monte Carlo), black points
for i in range(y_exact.shape[1]):
    plt.plot(x, y_exact[:, i], 'k.', markersize=9, 
             label='exact' if i == 0 else None)

plt.xscale('log')
plt.xlim(0.1, 100)
plt.xlabel('$b$', labelpad=10)
plt.ylim(0, 8)
plt.yticks([0, 2, 4, 6, 8])
plt.ylabel(r'med$[Z_0|1]$', labelpad=10)
plt.legend(loc='upper right', frameon=False)

plt.text(0.2, 0.8, '$s = 2$')
plt.text(0.2, 5.2, '$s = 5$')
plt.text(20, 2.8, '$s = 10$')

plt.subplots_adjust(left=0.15)
plt.subplots_adjust(bottom=0.15)

plt.show()