Coverage for lst_auto_rta/Theta_square_ACADA.py: 0%

1#!/usr/bin/env python 

2 

3import astropy 

4import gammapy 

5import matplotlib 

6import numpy as np 

7import regions 

8 

9print("gammapy:", gammapy.__version__) 

10print("numpy:", np.__version__) 

11print("astropy", astropy.__version__) 

12print("regions", regions.__version__) 

13print("matplotlib", matplotlib.__version__) 

14 

15import os 

16 

17import astropy.units as u 

18import matplotlib.pyplot as plt 

19import matplotlib.style as style 

20import numpy as np 

21from astropy.coordinates import SkyCoord 

22 

23style.use("tableau-colorblind10") 

24import argparse 

25from pathlib import Path 

26 

27from gammapy.data import DataStore 

28from gammapy.datasets import ( 

29 Datasets, 

30 MapDataset, 

31 SpectrumDataset, 

32) 

33from gammapy.makers import ( 

34 MapDatasetMaker, 

35 ReflectedRegionsBackgroundMaker, 

36 SafeMaskMaker, 

37 SpectrumDatasetMaker, 

38) 

39from gammapy.makers.utils import make_theta_squared_table 

40from gammapy.maps import Map, MapAxis, RegionGeom, WcsGeom 

41from gammapy.visualization import plot_spectrum_datasets_off_regions, plot_theta_squared_table 

42from matplotlib.offsetbox import AnchoredText 

43from regions import CircleSkyRegion 

44from scipy.optimize import curve_fit 

45 

46parser = argparse.ArgumentParser( 

47 description="Automatic Script for the DL1 check", formatter_class=argparse.ArgumentDefaultsHelpFormatter 

48) 

49parser.add_argument("-d", "--directory", default="/fefs/onsite/pipeline/rta/data/", help="Directory for data") 

50parser.add_argument("-da", "--date", default="20230705", help="Date of the run to check") 

51parser.add_argument("-r", "--run-id", default="13600", help="run id to check") 

52# parser.add_argument("-add", "--add-string", default="" , help="add a string to the path") 

53parser.add_argument("-add", "--add-string", default="reco/", help="add a string to the path") 

54parser.add_argument("-RA", "--right-ascension", default="82.0", help="right-ascension in deg") 

55parser.add_argument("-DEC", "--declination", default="23.0", help="declination in deg") 

56 

57args = parser.parse_args() 

58config = vars(args) 

59 

60print(config) 

61 

62location_data = ( 

63 config["directory"] + config["date"] + "/" + config["run_id"] + "/" + config["add_string"] + "/plots" 

64) # path to DL3 folder 

65source_name = config["run_id"] # e.g., Crab, GRB210807A 

66cut_type = "standard" # e.g., loose, hard, ... 

67filename_output = "{name}_{cut}".format(name=source_name, cut=cut_type) 

68 

69source_position = SkyCoord(ra=config["right_ascension"], dec=config["declination"], unit="deg", frame="icrs") 

70max_offset_run = 5 * u.deg 

71work_directory = location_data 

72path_plot = Path(work_directory + "/../plots") 

73print(work_directory + "/../plots") 

74path_plot.mkdir(exist_ok=True) 

75 

76 

77e_min = 0.05 * u.TeV 

78e_max = 20.0 * u.TeV 

79on_radius = 0.2 * u.deg 

80exclusion_radius = 0.35 * u.deg 

81max_theta_squared = (1.5 * exclusion_radius.to(u.deg).value) ** 2 

82n_bin_per_square_deg = 120 

83 

84data_store = DataStore.from_dir(location_data) 

85 

86print(data_store.obs_table) 

87 

88obs_ids = data_store.obs_table[source_position.separation(data_store.obs_table.pointing_radec) < max_offset_run][ 

89 "OBS_ID" 

90] 

91print(len(obs_ids)) 

92obs_collection = data_store.get_observations(obs_ids, required_irf="point-like") 

93 

94print("OBS_ID ", obs_ids) 

95 

96for i in range(len(obs_collection)): 

97 obs_collection[i].events.peek() 

98 plt.savefig(os.path.join(path_plot, "{}__obs.png".format(filename_output)), dpi=300) 

99 plt.close() 

100 

101on_region = CircleSkyRegion(center=source_position, radius=on_radius) 

102exclude_region = CircleSkyRegion(center=source_position, radius=exclusion_radius) 

103 

104energy_axis = MapAxis.from_edges( 

105 np.logspace(np.log10(e_min.to_value(u.TeV)), np.log10(e_max.to_value(u.TeV)), 2), 

106 unit="TeV", 

107 name="energy", 

108 interp="log", 

109) 

110geom_on_region = RegionGeom.create(region=on_region, axes=[energy_axis]) 

111geom = WcsGeom.create(skydir=source_position, npix=(200, 200), binsz=0.05, frame="icrs", axes=[energy_axis]) 

112 

113geom_image = geom.to_image() 

114exclusion_mask = ~geom_image.region_mask([exclude_region]) 

115# exclusion_mask.sum_over_axes().plot(); 

116 

117stacked = MapDataset.create(geom=geom, name=source_name + "_stacked") 

118unstacked = Datasets() 

119maker = MapDatasetMaker(selection=["counts"]) 

120maker_safe_mask = SafeMaskMaker(methods=["offset-max"], offset_max="3 deg") 

121 

122for obs in obs_collection: 

123 cutout = stacked.cutout(obs.pointing_radec, width="5 deg") 

124 # cutout = stacked.cutout(obs.get_pointing_icrs(), width="5 deg") 

125 dataset = maker.run(cutout, obs) 

126 dataset = maker_safe_mask.run(dataset, obs) 

127 stacked.stack(dataset) 

128 unstacked.append(dataset) 

129 

130stacked.to_image().counts.smooth("0.12 deg", kernel="disk").plot(add_cbar=True, cmap="coolwarm") 

131cbar = plt.gca().images[-1].colorbar 

132cbar.ax.get_yaxis().labelpad = 5 

133cbar.ax.set_ylabel(r"events/0.0004 deg$^2$", rotation=90) 

134plt.gca().scatter( 

135 source_position.ra, 

136 source_position.dec, 

137 transform=plt.gca().get_transform("world"), 

138 marker="+", 

139 c="k", 

140 label=source_name, 

141 s=100, 

142) 

143plt.legend() 

144 

145plt.savefig(os.path.join(path_plot, "{}__countmap.png".format(filename_output)), dpi=300) 

146plt.close() 

147 

148dataset_maker_spectrum_significance = SpectrumDatasetMaker(selection=["counts", "exposure"], use_region_center=True) 

149spectrum_dataset_empty_significance = SpectrumDataset.create(geom=geom_on_region) 

150bkg_maker_spectrum_significance = ReflectedRegionsBackgroundMaker(exclusion_mask=exclusion_mask) 

151 

152datasets_significance = Datasets() 

153for obs in obs_collection: 

154 dataset_spectrum_significance = dataset_maker_spectrum_significance.run( 

155 spectrum_dataset_empty_significance.copy(name=f"obs-{obs.obs_id}"), obs 

156 ) 

157 dataset_on_off_spectrum_significance = bkg_maker_spectrum_significance.run( 

158 observation=obs, dataset=dataset_spectrum_significance 

159 ) 

160 datasets_significance.append(dataset_on_off_spectrum_significance) 

161 

162info_table = datasets_significance.info_table(cumulative=True) 

163 

164 

165print("Livetime :", info_table["livetime"].to("s")[-1]) 

166print("Significance :", info_table["sqrt_ts"][-1]) 

167print("Excess :", info_table["excess"][-1]) 

168print("S/B :", info_table["excess"][-1] / info_table["background"][-1]) 

169print("Rate :", info_table["excess_rate"].to(1.0 / u.min)[-1].value, "per min") 

170 

171theta2_axis = MapAxis.from_bounds( 

172 0, max_theta_squared, nbin=int(n_bin_per_square_deg * max_theta_squared), interp="lin", unit="deg2" 

173) 

174 

175theta2_table = make_theta_squared_table( 

176 observations=obs_collection, 

177 position=source_position, 

178 theta_squared_axis=theta2_axis, 

179) 

180 

181# Values to be used for the plots 

182counts_on = theta2_table["counts"] 

183counts_on_err = np.sqrt(counts_on) 

184 

185counts_off = theta2_table["counts_off"] 

186counts_off_err = np.sqrt(counts_off) 

187 

188excess = theta2_table["excess"] 

189excess_err = np.append( 

190 np.abs(theta2_table["excess_errn"].data[np.newaxis, :]), theta2_table["excess_errp"].data[np.newaxis, :], axis=0 

191) 

192 

193bin_center = theta2_axis.center 

194theta2_cut = on_radius.to_value(u.deg) ** 2 

195 

196N_off_err = np.sqrt(info_table["counts_off"][-1]) 

197significance_lima = info_table["sqrt_ts"][-1] 

198text = ( 

199 r"Livetime = {:.1f}" 

200 f"\n" 

201 r"N$_{{\rm on}}$ = {:.0f}" 

202 f"\n" 

203 r"1/$\alpha$ = {:.1f}" 

204 f"\n" 

205 r"N$_{{\rm off}}$ = {:.1f} $\pm$ {:.1f}" 

206 f"\n" 

207 r"N_excess = {:.1f}" 

208 f"\n" 

209 r"Significance (Li&Ma) = {:.1f} $\sigma$ ".format( 

210 info_table["livetime"].to("s")[-1], 

211 info_table["counts"][-1], 

212 1 / info_table["alpha"][-1], 

213 info_table["counts_off"][-1], 

214 N_off_err, 

215 info_table["excess"][-1], 

216 significance_lima, 

217 ) 

218) 

219 

220fig, ax = plt.subplots(figsize=(8, 6)) 

221 

222ax.errorbar( 

223 bin_center, 

224 counts_on, 

225 xerr=theta2_axis.bin_width / 2, 

226 yerr=counts_on_err, 

227 fmt="o", 

228 label="ON data", 

229 ms=10, 

230 color="tab:orange", 

231) 

232ax.errorbar( 

233 bin_center, 

234 counts_off, 

235 xerr=theta2_axis.bin_width / 2, 

236 yerr=counts_off_err, 

237 fmt="s", 

238 label="OFF data", 

239 ms=10, 

240 color="tab:blue", 

241) 

242ax.set_xlim(left=0) 

243ax.grid(ls="-") 

244ax.axvline(theta2_cut, color="black", ls="--", alpha=0.75) 

245ax.set_xlabel("$\\theta^2$ [deg$^2$]") 

246ax.set_ylabel("Counts") 

247ax.legend(loc="center right") 

248ax.set_title(source_name) 

249 

250box_prop = dict(boxstyle="Round", facecolor="white", alpha=0.5) 

251text_prop = dict(fontsize="x-large", bbox=box_prop) 

252txt = AnchoredText(text, loc=1, transform=ax.transAxes, prop=text_prop, frameon=False) 

253ax.add_artist(txt) 

254 

255plt.savefig(os.path.join(path_plot, "{}_theta2.png".format(filename_output)), dpi=300) 

256plt.close() 

257 

258# Excess plot 

259fig, ax = plt.subplots(figsize=(12, 10)) 

260 

261ax.errorbar(bin_center, excess, yerr=excess_err, fmt="o", color="tab:blue", label="Excess") 

262ax.bar(bin_center, excess, theta2_axis.bin_width.to_value(u.deg**2), color="lightblue", ec="lightblue", alpha=0.5) 

263# ax.axhline(0) 

264ax.set_xlim(left=0) 

265ax.axvline(theta2_cut, color="black", ls="--", alpha=0.75) 

266ax.grid(ls="-") 

267ax.set_xlabel("$\\theta^2$ [deg$^2$]") 

268ax.set_ylabel("Excess") 

269ax.legend(loc="upper right", title=f"Significance = {significance_lima:.1f} $\sigma$") 

270 

271plt.savefig(os.path.join(path_plot, "{}_theta2_excess.png".format(filename_output)), dpi=300) 

272plt.close() 

273 

274excess_err_livetime = np.sqrt( 

275 np.sqrt(info_table["counts"]) ** 2 + (np.sqrt(info_table["counts_off"]) * info_table["alpha"]) ** 2 

276) 

277 

278excess_livetime_function = lambda x, a: a * x 

279significance_livetime_function = lambda x, a: a * np.sqrt(x) 

280 

281popt_excess_livetime, pcov_excess_livetime = curve_fit( 

282 excess_livetime_function, info_table["livetime"].to(u.h), info_table["excess"], sigma=excess_err_livetime 

283) 

284popt_significance_livetime, pcov_significance_livetime = curve_fit( 

285 significance_livetime_function, info_table["livetime"].to(u.h), info_table["sqrt_ts"] 

286) 

287 

288popt_excess_livetime = popt_excess_livetime[0] * u.Unit("h^-1") 

289pcov_excess_livetime = pcov_excess_livetime[0, 0] * u.Unit("h^-2") 

290popt_significance_livetime = popt_significance_livetime[0] * u.Unit("h^(-1/2)") 

291pcov_significance_livetime = pcov_significance_livetime[0, 0] * u.Unit("h^(-1)") 

292 

293x_plot_livetime = np.linspace(0, np.max(info_table["livetime"].to(u.h)) * 1.2, num=len(info_table["livetime"]) * 10) 

294 

295y_plot_fit_excess_livetime = excess_livetime_function(x_plot_livetime, popt_excess_livetime) 

296y_err_plot_fit_excess_livetime = x_plot_livetime * np.sqrt(pcov_excess_livetime) 

297 

298y_plot_fit_significance_livetime = significance_livetime_function(x_plot_livetime, popt_significance_livetime) 

299y_err_plot_fit_significance_livetime = np.sqrt(x_plot_livetime) * np.sqrt(pcov_significance_livetime) 

300 

301fig, ax = plt.subplots(figsize=(15, 7), ncols=2) 

302ax1, ax2 = ax[0], ax[1] 

303ax1.fill_between( 

304 x_plot_livetime.to_value(u.h), 

305 y_plot_fit_excess_livetime - y_err_plot_fit_excess_livetime, 

306 y_plot_fit_excess_livetime + y_err_plot_fit_excess_livetime, 

307 color="silver", 

308 alpha=0.5, 

309) 

310ax1.plot(x_plot_livetime, y_plot_fit_excess_livetime, c="silver") 

311ax1.errorbar( 

312 info_table["livetime"].to(u.h), 

313 info_table["excess"], 

314 yerr=excess_err_livetime, 

315 fmt="o", 

316 color="tab:blue", 

317 label="Excess", 

318) 

319ax1.set_xlabel("Livetime [h]") 

320ax1.set_ylabel("Excess") 

321txt_ax1 = AnchoredText( 

322 "Excess rate : {:.1f} +/- {:.1f} ".format(popt_excess_livetime.value, np.sqrt(pcov_excess_livetime).value) 

323 + str(popt_excess_livetime.unit), 

324 loc=2, 

325 transform=ax1.transAxes, 

326 prop=text_prop, 

327 frameon=False, 

328) 

329ax1.add_artist(txt_ax1) 

330 

331ax2.fill_between( 

332 x_plot_livetime.to_value(u.h), 

333 y_plot_fit_significance_livetime - y_err_plot_fit_significance_livetime, 

334 y_plot_fit_significance_livetime + y_err_plot_fit_significance_livetime, 

335 color="silver", 

336 alpha=0.5, 

337) 

338ax2.plot(x_plot_livetime, y_plot_fit_significance_livetime, c="silver") 

339ax2.errorbar(info_table["livetime"].to(u.h), info_table["sqrt_ts"], fmt="o", color="tab:blue", label="Excess") 

340txt_ax2 = AnchoredText( 

341 "Significance rate : {:.1f} +/- {:.1f} ".format( 

342 popt_significance_livetime.value, np.sqrt(pcov_significance_livetime).value 

343 ) 

344 + str(popt_significance_livetime.unit), 

345 loc=2, 

346 transform=ax2.transAxes, 

347 prop=text_prop, 

348 frameon=False, 

349) 

350ax2.add_artist(txt_ax2) 

351ax2.set_xlabel("Livetime [h]") 

352ax2.set_ylabel("Significance") 

353 

354plt.savefig(os.path.join(path_plot, "{}__excess_significance_over_time.png".format(filename_output)), dpi=300) 

355plt.close() 

356 

357ax = stacked.to_image().counts.smooth("0.12 deg", kernel="disk").plot(add_cbar=True, cmap="coolwarm") 

358on_region.to_pixel(ax.wcs).plot(ax=ax, color="black") 

359plot_spectrum_datasets_off_regions(datasets=datasets_significance, ax=ax) 

360plt.gca().scatter( 

361 source_position.ra, 

362 source_position.dec, 

363 transform=plt.gca().get_transform("world"), 

364 marker="+", 

365 c="k", 

366 label=source_name, 

367 s=[100], 

368) 

369plt.savefig(os.path.join(path_plot, "{}__countmap_with_OFF.png".format(filename_output)), dpi=300) 

370plt.close()