class event:
def __init__(self, datastore_path, runid):
# self.table = table
self.table = 0
self.count_scatter = 0
self.info_table = 0
self.stats = 0
self.sourcePos = 0
self.ringbg_exclude_region = 0
self.reflectedbg_exclude_region = 0
reflected_exclusion_mask = 0
self.skydir = 0
self.exclusion_mask = 0
self.on_region = 0
self.datasets = 0
self.datastore = 0
self.obs_id = 0
self.observations = 0
self.create_datastore(datastore_path, runid)
self.dataset_maker = 0
self.dataset_empty = 0
self.bkg_maker = 0
self.safe_mask_masker = 0
self.signal_table = 0
model_best_joint = 0
self.flux_points = 0
self.model = 0
self.model_best_stacked = 0
self.obsCollection = 0
self.stacked = 0
self.npix_x = 0
self.npix_y = 0
self.ring_exclusion_mask = 0
self.significance_map = 0
self.excess_map = 0
self.axis = 0
self.exclusion_map = 0
def create_datastore(self, datastore_path, runid):
datastore = DataStore.from_dir(datastore_path)
obs_id = runid
observations = datastore.get_observations(obs_id)
self.obs_id = obs_id
self.datastore = datastore
self.observations = observations
self.table = []
for i in obs_id:
self.table.append(datastore.obs(obs_id=i).events.table)
def plot_timee(self, id_=0, ax=None):
"""Plots an event rate time curve.
Parameters
----------
ax : `~matplotlib.axes.Axes` or None
Axes
Returns
-------
ax : `~matplotlib.axes.Axes`
Axes
i
"""
plt.figure(figsize=(16, 8))
for i in range(len(self.table)):
ax = plt.gca() if ax is None else ax
# Note the events are not necessarily in time order
# print(self.table[i])
time = self.table[i]["TIME"]
time = time - np.min(time)
ax.set_xlabel("Time (sec)")
ax.set_ylabel("Counts")
y, x_edges = np.histogram(time, bins=np.linspace(0, 2400, 100))
y = y[:-1]
x_edges = x_edges[:-1]
xerr = np.diff(x_edges) / 2
x = x_edges[:-1] + xerr
yerr = np.sqrt(y)
ax.errorbar(x=x, y=y, xerr=xerr, yerr=yerr, label="run" + str(self.obs_id[i]))
plt.legend(
shadow=True, bbox_to_anchor=(1.05, 1), loc="upper left", borderaxespad=0, handlelength=1.5, fontsize=10
)
prefix = (5 - len(str(id_))) * "0"
plt.savefig("../plots/event_rate" + prefix + str(id_) + ".png")
plt.show()
return x, y
def create_counts_scatter(self):
time = self.table["TIME"]
time = time - np.min(time)
y, x_edges = np.histogram(time, bins=np.linspace(0, 1200, 20))
y = y[:-1]
self.count_scatter = y
def statisticals_parameters(self, runid):
info_table = self.datastore.obs_table
self.create_counts_scatter()
variation_max = (
np.abs((self.count_scatter.max() - self.count_scatter.mean()) / self.count_scatter.mean()) * 100
)
variation_min = (
np.abs((self.count_scatter.min() - self.count_scatter.mean()) / self.count_scatter.mean()) * 100
)
general_variation = self.count_scatter.std() / self.count_scatter.mean() * 100
variability = self.count_scatter.std() / (np.sqrt(self.count_scatter.mean()))
mean_trigger_rate = self.count_scatter.mean() / info_table["LIVETIME"][runid]
zenith_angle = info_table["ZEN_PNT"][runid]
livetime = info_table["LIVETIME"][runid]
print(f"Mean = {self.count_scatter.mean()}")
print(f"Standard deviation = {self.count_scatter.std()}")
print(f"Localized variation on the maximum = {variation_max} %")
print(f"Localized variation on the maximum = {variation_min} %")
print(f"General variation = {general_variation}")
print(f"Variability = {variability}")
print(f"Mean_trigger_rate = {mean_trigger_rate}")
print(f"Zenith_angle = {zenith_angle}")
print(f"Livetime = {livetime}")
self.info_table = info_table
return (
self.count_scatter.mean(),
self.count_scatter.std(),
variation_max,
variation_min,
general_variation,
variability,
mean_trigger_rate,
zenith_angle,
livetime,
)
def create_on_region(self, ra=83.633, dec=22.014, on_radius_angle=0.2, radius_excluded=0.35):
# def create_on_exclusion_regions(self, ra=238.929, dec=11.190, on_radius_angle=0.14, radius_excluded=0.2):
# def create_on_exclusion_regions(self, ra=187.706, dec=12.391, on_radius_angle=0.14, radius_excluded=0.2):
sourcePos = SkyCoord(ra=ra * u.deg, dec=dec * u.deg, frame="icrs")
on_radius = on_radius_angle * u.deg
on_region = CircleSkyRegion(center=sourcePos, radius=on_radius)
# Create exclusion region
reflectedbg_exclude_region = CircleSkyRegion(center=sourcePos, radius=radius_excluded * u.deg)
ringbg_exclude_region = CircleSkyRegion(center=sourcePos, radius=on_radius_angle * 3 * u.deg)
self.sourcePos = sourcePos
self.on_region = on_region
self.reflectedbg_exclude_region = reflectedbg_exclude_region
self.ringbg_exclude_region = ringbg_exclude_region
# Reflected Background Method
def create_exclusion_mask(self, plot=False, binsz=0.02, npix_x=400, npix_y=400):
skydir = self.sourcePos.galactic
reflected_exclusion_mask = Map.create(
npix=(npix_x, npix_y), binsz=binsz, skydir=skydir, proj="TAN", frame="icrs"
)
mask = reflected_exclusion_mask.geom.region_mask([self.reflectedbg_exclude_region], inside=False)
reflected_exclusion_mask.data = mask
self.skydir = skydir
self.reflected_exclusion_mask = reflected_exclusion_mask
if plot:
mask.plot()
plt.show()
def reduction_chain(
self, e_reco_min=0.1, e_reco_max=40, e_reco_bin=40, e_true_min=0.05, e_true_max=100, e_true_bin=200
):
# Run data reduction chain, e_=energy
e_reco = MapAxis.from_energy_bounds(e_reco_min, e_reco_max, e_reco_bin, unit="TeV", name="energy")
e_true = MapAxis.from_energy_bounds(e_true_min, e_true_max, e_true_bin, unit="TeV", name="energy_true")
dataset_empty = SpectrumDataset.create(e_reco=e_reco, e_true=e_true, region=self.on_region)
self.dataset_empty = dataset_empty
def data_maker(self):
dataset_maker = SpectrumDatasetMaker(containment_correction=False, selection=["counts", "exposure", "edisp"])
bkg_maker = ReflectedRegionsBackgroundMaker(exclusion_mask=self.reflected_exclusion_mask)
safe_mask_masker = SafeMaskMaker(methods=["aeff-max"], aeff_percent=10)
self.dataset_maker = dataset_maker
self.bkg_maker = bkg_maker
self.safe_mask_masker = safe_mask_masker
def data_run(self):
datasets = Datasets()
observations = self.datastore.get_observations(self.obs_id)
for obs_id, observation in zip(self.obs_id, observations):
dataset = self.dataset_maker.run(self.dataset_empty.copy(name=str(obs_id)), observation)
dataset_on_off = self.bkg_maker.run(dataset, observation)
dataset_on_off = self.safe_mask_masker.run(dataset_on_off, observation)
datasets.append(dataset_on_off)
self.datasets = datasets
def plot_off_region(self):
plt.figure(figsize=(8, 8))
_, ax, _ = self.reflected_exclusion_mask.plot()
self.on_region.to_pixel(ax.wcs).plot(ax=ax, edgecolor="k")
plot_spectrum_datasets_off_regions(ax=ax, datasets=self.datasets)
plt.show()
def signal_info(self):
info_table = self.datastore.obs_table
print(info_table[0])
signal_table = self.datasets.info_table(cumulative=True)
alpha = signal_table["alpha"]
excess = signal_table["excess"]
uncertainty_excess = np.sqrt(
signal_table["counts"] + ((1 / alpha) * signal_table["background"]) / ((1 / alpha) ** 2)
)
bg = signal_table["background"]
uncertainty_bg = np.sqrt((1 / alpha) * signal_table["background"]) / (1 / alpha)
self.signal_table = signal_table
print(f"reflected_excess = {excess[-1]}")
print(f"reflected_uncertainty_excess = {uncertainty_excess[-1]}")
print(f"reflected_uncertainty_excess = {self.signal_table['livetime'].to('h')[-1]} hours")
print(f"sqrt_sts = {self.signal_table['sqrt_ts'][-1]}")
# print(f"reflected_bg = {bg}")
# print(f"reflected_uncertainty_bg = {uncertainty_bg}")
return (
excess[-1],
uncertainty_excess[-1],
bg[-1],
uncertainty_bg[-1],
self.signal_table["sqrt_ts"][-1],
self.signal_table["livetime"][-1],
)
def statistic_source_excess(self, plot=False):
if plot:
plt.plot(
# self.signal_table["livetime"].to("h"),
self.signal_table["name"],
self.signal_table["excess"],
marker="o",
ls="none",
)
plt.xlabel("Livetime [h]")
plt.ylabel("Excess")
plt.show()
def statistic_source_ts(self, plot=False):
if plot:
plt.plot(self.signal_table["livetime"].to("h"), self.signal_table["sqrt_ts"], marker="o", ls="none")
plt.xlabel("Livetime [h]")
plt.ylabel("Sqrt(TS)")
plt.show()
def fit_spectrum(self):
# Fit spectrum
spectral_model = PowerLawSpectralModel(
index=2, amplitude=2e-11 * u.Unit("cm-2 s-1 TeV-1"), reference=1 * u.TeV
)
model = SkyModel(spectral_model=spectral_model, name="crab")
for self.dataset in self.datasets:
self.dataset.models = model
fit_joint = Fit(self.datasets)
result_joint = fit_joint.run()
# Make a copy here to compare it later
model_best_joint = model.copy()
self.model = model
self.model_best_joint = model_best_joint
def fit_quality(self, plot=False):
if plot:
ax_spectrum, ax_residuals = self.datasets[0].plot_fit()
ax_spectrum.set_ylim(0.1, 40)
plt.show()
def flux_points_plot(self, e_min=0.7, e_max=30, plot=False):
# Compute Flux Points
energy_edges = np.logspace(np.log10(e_min), np.log10(e_max), 11) * u.TeV
# Create an instance of the FluxPointsEstimator
fpe = FluxPointsEstimator(energy_edges=energy_edges, source="crab")
flux_points = fpe.run(datasets=self.datasets)
flux_points.table_formatted
# Plot the flux points and their likelihood profile with Treshold < 4
if plot:
plt.figure(figsize=(8, 5))
flux_points.table["is_ul"] = flux_points.table["ts"] < 4
ax = flux_points.plot(energy_power=2, flux_unit="erg-1 cm-2 s-1", color="darkorange")
flux_points.to_sed_type("e2dnde").plot_ts_profiles(ax=ax)
plt.show()
self.flux_points = flux_points
def model_fit_plot(self, plot=False):
# Final plot with the best fit model
flux_points_dataset = FluxPointsDataset(data=self.flux_points, models=self.model_best_joint)
if plot:
flux_points_dataset.plot_fit()
plt.show()
def flux_plot(self, plot=False, id_=0):
dataset_stacked = Datasets(self.datasets).stack_reduce()
dataset_stacked.models = self.model
stacked_fit = Fit([dataset_stacked])
result_stacked = stacked_fit.run()
# Make a copy to compare later
model_best_stacked = self.model.copy()
if plot:
plt.figure(figsize=(16, 8))
plot_kwargs = {
"energy_range": [0.1, 30] * u.TeV,
"energy_power": 2,
"flux_unit": "erg-1 cm-2 s-1",
}
# courb stacked model
model_best_stacked.spectral_model.plot(**plot_kwargs, label="run" + str(self.obs_id))
model_best_stacked.spectral_model.plot_error(**plot_kwargs)
# courb reference
create_crab_spectral_model("hess_pl").plot(**plot_kwargs, label="Crab reference")
plt.legend()
prefix = (5 - len(str(id_))) * "0"
plt.savefig("./Plots/flux_plot" + prefix + str(id_) + ".png")
plt.show()
self.model_best_stacked = model_best_stacked
def flux_parameters(self):
parameters_table = self.model_best_stacked.parameters.to_table()
spectral_index = parameters_table[0][1]
uncertainty_spectral_index = parameters_table[0][6]
amplitude = parameters_table[1][1]
uncertainty_amplitude = parameters_table[1][6]
return spectral_index, uncertainty_spectral_index, amplitude, uncertainty_amplitude
"""Ring Background Method"""
def calculate_acceptance_model(
self,
energyaxis_min=np.log10(0.1),
energyaxis_max=np.log10(2.0),
energy_nbin=5,
offset_min=0.0,
offset_max=5.0,
offset_nbin=8,
plot=False,
):
obsCollection = self.datastore.get_observations(self.obs_id)
energyAxisAcceptance = MapAxis.from_edges(
np.logspace(energyaxis_min, energyaxis_max, energy_nbin), unit="TeV", name="energy", interp="log"
)
offsetAxisAcceptance = MapAxis.from_edges(
np.linspace(offset_min, offset_max, offset_nbin), unit="deg", name="offset", interp="lin"
)
background = create_radial_acceptance_map(
obsCollection,
energyAxisAcceptance,
offsetAxisAcceptance,
exclude_regions=[self.ringbg_exclude_region],
oversample_map=10,
)
if plot:
background.peek()
plt.show
hduBackground = background.to_table_hdu()
hduBackground.writeto("background.fits", overwrite=True)
self.obsCollection = obsCollection
def add_acceptance_map(self):
# Add acceptance map to observations
listIDObs = []
for obs in self.obsCollection:
listIDObs.append(obs.obs_id)
self.datastore.hdu_table.add_row(
{
"OBS_ID": listIDObs[-1],
"HDU_TYPE": "bkg",
"HDU_CLASS": "bkg_2d",
"FILE_DIR": "",
"FILE_NAME": "background.fits",
"HDU_NAME": "BACKGROUND",
}
)
self.datastore.hdu_table = self.datastore.hdu_table.copy()
self.datastore.hdu_table
self.obsCollection = self.datastore.get_observations(listIDObs)
def map_geometry(
self,
RA,
DEC,
plot=False,
npix_x=200,
npix_y=200,
binsize=0.02,
axis_min=np.log10(0.1),
axis_max=np.log10(2.0),
axis_nbin=5,
):
sourcePos = SkyCoord(ra=RA * u.deg, dec=DEC * u.deg, frame="icrs")
self.axis = MapAxis.from_edges(
np.logspace(axis_min, axis_max, axis_nbin), unit="TeV", name="energy", interp="log"
)
geom = WcsGeom.create(skydir=sourcePos, npix=(npix_x, npix_y), binsz=binsize, frame="icrs", axes=[self.axis])
geom
ring_exclusion_mask = geom.to_image().region_mask([self.ringbg_exclude_region], inside=False)
self.exclusion_map = WcsNDMap(geom.to_image(), ring_exclusion_mask)
self.ring_exclusion_mask = ring_exclusion_mask
if plot:
self.exclusion_map.plot()
plt.show()
stacked = MapDataset.create(geom=geom)
unstacked = Datasets()
maker = MapDatasetMaker(selection=["counts", "background"])
# maker = MapDatasetMaker(selection=["counts"])
maker_safe_mask = SafeMaskMaker(methods=["offset-max"], offset_max=3.0 * u.deg)
for obs in self.obsCollection:
cutout = stacked.cutout(obs.pointing_radec, width="10 deg")
dataset = maker.run(cutout, obs)
dataset = maker_safe_mask.run(dataset, obs)
stacked.stack(dataset)
unstacked.append(dataset)
self.stacked = stacked
self.npix_x = npix_x
self.npix_y = npix_y
def ring_data_estimation(self):
ring_maker = RingBackgroundMaker(r_in="0.5 deg", width="0.3 deg", exclusion_mask=self.ring_exclusion_mask)
estimator = ExcessMapEstimator(0.2 * u.deg)
lima_maps = estimator.run(self.stacked)
significance_map = lima_maps["sqrt_ts"]
excess_map = lima_maps["excess"]
npix_x = int(self.npix_x / 2)
npix_y = int(self.npix_y / 2)
ring_sqrt_ts = lima_maps["sqrt_ts"].data[0][npix_x][npix_y]
ring_excess = lima_maps["excess"].data[0][npix_x][npix_y]
ring_bg = lima_maps["background"].data[0][npix_x][npix_y]
uncertainty_excess_ring = lima_maps["err"].data[0][npix_x][npix_y]
# uncertainty_ring_bg = np.sqrt((1/self.alpha)*ring_bg)/(1/self.alpha)
self.significance_map = significance_map
self.excess_map = excess_map
print(f"ring_excess = {ring_excess}")
print(f"ring_sqrt_ts = {ring_sqrt_ts}")
print(f"ring_bg = {ring_bg}")
print(f"uncertainty_excess_ring = {uncertainty_excess_ring}")
# print(f"uncertainty_ring_bg{uncertainty_ring_bg}")
return ring_excess, ring_sqrt_ts, ring_bg, uncertainty_excess_ring
# uncertainty_ring_bg
def excess_significance_plot(self, directory, plot=False, id_=0):
if plot:
plt.figure(figsize=(16, 8))
ax1 = plt.subplot(121, projection=self.significance_map.geom.wcs)
ax2 = plt.subplot(122, projection=self.excess_map.geom.wcs)
ax2.set_title("Significance map")
self.significance_map.plot(ax=ax2, add_cbar=True)
sources = find_peaks(
self.significance_map.get_image_by_idx((0,)),
threshold=5,
min_distance="0.2 deg",
)
print("Found sources")
print(sources)
f = open(directory + "../plots/results.txt", "w")
f.write("Found sources")
f.write(str(sources))
f.close()
sources_2 = find_peaks(
self.significance_map.get_image_by_idx((0,)),
threshold=7,
min_distance="0.2 deg",
)
now = dt.datetime.now()
timestamp_str = now.strftime("%Y-%m-%d %H:%M:%S")
ax1.text(
0.02,
0.98,
timestamp_str,
transform=ax1.transAxes,
fontsize=11,
fontweight="bold",
va="top",
ha="left",
)
ax2.text(
0.02,
0.98,
timestamp_str,
transform=ax2.transAxes,
fontsize=11,
fontweight="bold",
va="top",
ha="left",
)
if len(sources) > 0:
ax2.scatter(
sources["ra"],
sources["dec"],
transform=plt.gca().get_transform("icrs"),
color="none",
edgecolor="white",
marker="o",
s=300,
lw=1.5,
)
if len(sources_2) > 0:
ax2.scatter(
sources_2["ra"],
sources_2["dec"],
transform=plt.gca().get_transform("icrs"),
color="none",
edgecolor="blue",
marker="o",
s=300,
lw=1.5,
)
ax1.set_title("Excess map")
self.excess_map.plot(ax=ax1, add_cbar=True)
prefix = (5 - len(str(id_))) * "0"
plt.savefig(directory + "../plots/sig_excess_plot" + prefix + str(id_) + ".png")
plt.show()
return sources
def off_distribution_plot(self, plot=False):
# create a 2D mask for the images
exclusion_map_ring = self.significance_map.geom.region_mask([self.ringbg_exclude_region], inside=False)
significance_map_off = self.significance_map * exclusion_map_ring
significance_all = self.significance_map.data[np.isfinite(self.significance_map.data)]
significance_off = significance_map_off.data[np.isfinite(significance_map_off.data)]
bins = np.linspace(
np.min(significance_all),
np.max(significance_all),
num=int(np.max(significance_all - np.min(significance_all)) * 3),
)
mu, std = norm.fit(significance_off)
if plot:
plt.hist(
significance_all,
density=True,
alpha=0.5,
color="red",
label="all bins",
bins=bins,
)
plt.hist(
significance_off,
density=True,
alpha=0.5,
color="blue",
label="off bins",
bins=bins,
)
# Now, fit the off distribution with a Gaussian
x = np.linspace(-8, 8, 50)
p = norm.pdf(x, mu, std)
plt.plot(x, p, lw=2, color="black")
plt.legend()
plt.xlabel("Significance")
plt.yscale("log")
plt.ylim(1e-5, 1)
xmin, xmax = np.min(significance_all), np.max(significance_all)
plt.xlim(xmin, xmax)
plt.show()
print(mu, std)
return mu, std
def thetaSquarePlot(self):
# theta2Edge = np.linspace(0.0, 0.35**2, num=15)
theta2Edge = np.linspace(0.0, 0.35**2, num=15)
thetaEdge = np.sqrt(theta2Edge)
theta2 = (theta2Edge[1:] + theta2Edge[:-1]) / 2.0
count = np.zeros(theta2.shape)
countBack = np.zeros(theta2.shape)
alpha = np.zeros(theta2.shape)
significance = np.zeros(theta2.shape)
sb = np.zeros(theta2.shape)
for i in range(len(theta2)):
on_region_theta2_ring = CircleAnnulusSkyRegion(
center=self.sourcePos, inner_radius=thetaEdge[i] * u.deg, outer_radius=thetaEdge[i + 1] * u.deg
)
on_region_theta2_circle = CircleSkyRegion(center=self.sourcePos, radius=thetaEdge[i + 1] * u.deg)
dataset_maker_spectrum_significance_theta2 = SpectrumDatasetMaker(
selection=["counts", "exposure", "edisp"]
)
spectrum_dataset_empty_significance_theta2_ring = SpectrumDataset.create(
e_reco=self.axis, region=on_region_theta2_ring
)
spectrum_dataset_empty_significance_theta2_circle = SpectrumDataset.create(
e_reco=self.axis, region=on_region_theta2_circle
)
bkg_maker_spectrum_significance_theta2 = ReflectedRegionsBackgroundMaker(
exclusion_mask=self.exclusion_map
)
countsRing = np.zeros(len(self.obsCollection))
countsOffRing = np.zeros(len(self.obsCollection))
alphaRing = np.zeros(len(self.obsCollection))
exposureRing = np.zeros(len(self.obsCollection))
countsCircle = np.zeros(len(self.obsCollection))
countsOffCircle = np.zeros(len(self.obsCollection))
alphaCircle = np.zeros(len(self.obsCollection))
exposureCircle = np.zeros(len(self.obsCollection))
for j, obs in enumerate(self.obsCollection):
print(i)
dataset_spectrum_significance_theta2_ring = dataset_maker_spectrum_significance_theta2.run(
spectrum_dataset_empty_significance_theta2_ring.copy(name=f"obs-{obs.obs_id}"), obs
)
dataset_on_off_spectrum_significance_theta2_ring = bkg_maker_spectrum_significance_theta2.run(
observation=obs, dataset=dataset_spectrum_significance_theta2_ring
)
countsRing[j] = dataset_on_off_spectrum_significance_theta2_ring.counts.get_by_idx([0])[0, 0, 0]
countsOffRing[j] = dataset_on_off_spectrum_significance_theta2_ring.counts_off.get_by_idx([0])[
0, 0, 0
]
alphaRing[j] = dataset_on_off_spectrum_significance_theta2_ring.alpha.get_by_idx([0])[0, 0, 0]
exposureRing[j] = dataset_on_off_spectrum_significance_theta2_ring.exposure.get_by_idx([0])[0, 0, 0]
dataset_spectrum_significance_theta2_circle = dataset_maker_spectrum_significance_theta2.run(
spectrum_dataset_empty_significance_theta2_circle.copy(name=f"obs-{obs.obs_id}"), obs
)
dataset_on_off_spectrum_significance_theta2_circle = bkg_maker_spectrum_significance_theta2.run(
observation=obs, dataset=dataset_spectrum_significance_theta2_circle
)
countsCircle[j] = dataset_on_off_spectrum_significance_theta2_circle.counts.get_by_idx([0])[0, 0, 0]
countsOffCircle[j] = dataset_on_off_spectrum_significance_theta2_circle.counts_off.get_by_idx([0])[
0, 0, 0
]
alphaCircle[j] = dataset_on_off_spectrum_significance_theta2_circle.alpha.get_by_idx([0])[0, 0, 0]
exposureCircle[j] = dataset_on_off_spectrum_significance_theta2_circle.exposure.get_by_idx([0])[
0, 0, 0
]
countsRing = np.sum(countsRing)
countsOffRing = np.sum(countsOffRing)
alphaRing = np.sum(alphaRing * exposureRing / np.sum(exposureRing))
countsCircle = np.sum(countsCircle)
countsOffCircle = np.sum(countsOffCircle)
alphaCircle = np.sum(alphaCircle * exposureCircle / np.sum(exposureCircle))
wstatCircle = WStatCountsStatistic(n_on=countsCircle, n_off=countsOffCircle, alpha=alphaCircle)
count[i] = countsRing
countBack[i] = countsOffRing
alpha[i] = alphaRing
significance[i] = wstatCircle.sqrt_ts
sb[i] = wstatCircle.n_sig / wstatCircle.n_bkg
plt.figure(figsize=(30, 10))
plt.subplot(1, 3, 1)
plt.errorbar(
x=theta2,
xerr=[theta2 - theta2Edge[:-1], theta2Edge[1:] - theta2],
y=count,
yerr=np.sqrt(count),
label="On",
fmt=".",
)
plt.errorbar(
x=theta2,
xerr=[theta2 - theta2Edge[:-1], theta2Edge[1:] - theta2],
y=countBack * alpha,
yerr=np.sqrt(countBack) * alpha,
label="Off",
fmt=".",
)
plt.legend()
plt.xlabel("Theta^2")
plt.ylabel("Count")
plt.axvline(x=0.14**2, ls="--", c="k")
plt.subplot(1, 3, 2)
plt.plot(theta2Edge[:-1], significance, "+", mew=3.0)
plt.xlabel("Theta^2")
plt.ylabel("Significance")
plt.axvline(x=0.14**2, ls="--", c="k")
plt.subplot(1, 3, 3)
plt.plot(theta2Edge[:-1], sb, "+", mew=3.0)
plt.xlabel("Theta^2")
plt.ylabel("S/B")
plt.axvline(x=0.14**2, ls="--", c="k")
return 0