Coverage for src / lstautorta / SourceAnalyse.py: 0%

1## Script 

2 

3# %matplotlib inline 

4### Licensed under a 3-clause BSD style license - see LICENSE.rst 

5import datetime as dt 

6import logging 

7 

8import astropy.units as u 

9import matplotlib 

10import matplotlib.pyplot as plt 

11import numpy as np 

12from acceptanceCalculation import create_radial_acceptance_map 

13from astropy.coordinates import SkyCoord 

14from gammapy.data import DataStore 

15from gammapy.datasets import ( 

16 Datasets, 

17 FluxPointsDataset, 

18 MapDataset, 

19 SpectrumDataset, 

20) 

21from gammapy.estimators import ExcessMapEstimator, FluxPointsEstimator 

22from gammapy.estimators.utils import find_peaks 

23from gammapy.makers import ( 

24 MapDatasetMaker, 

25 ReflectedRegionsBackgroundMaker, 

26 RingBackgroundMaker, 

27 SafeMaskMaker, 

28 SpectrumDatasetMaker, 

29) 

30from gammapy.maps import Map, MapAxis, WcsGeom, WcsNDMap 

31from gammapy.modeling import Fit 

32from gammapy.modeling.models import ( 

33 PowerLawSpectralModel, 

34 SkyModel, 

35 create_crab_spectral_model, 

36) 

37from gammapy.stats import WStatCountsStatistic 

38from gammapy.visualization import ( 

39 plot_spectrum_datasets_off_regions, 

40) 

41from regions import CircleAnnulusSkyRegion, CircleSkyRegion 

42from scipy.stats import norm 

43 

44# import astropy 

45 

46 

47__all__ = ["event"] 

48 

49log = logging.getLogger(__name__) 

50matplotlib.rcParams.update({"font.size": 17}) 

51 

52 

53class event: 

54 def __init__(self, datastore_path, runid): 

55 # self.table = table 

56 self.table = 0 

57 self.count_scatter = 0 

58 self.info_table = 0 

59 self.stats = 0 

60 self.sourcePos = 0 

61 self.ringbg_exclude_region = 0 

62 self.reflectedbg_exclude_region = 0 

63 reflected_exclusion_mask = 0 

64 self.skydir = 0 

65 self.exclusion_mask = 0 

66 self.on_region = 0 

67 self.datasets = 0 

68 self.datastore = 0 

69 self.obs_id = 0 

70 self.observations = 0 

71 self.create_datastore(datastore_path, runid) 

72 self.dataset_maker = 0 

73 self.dataset_empty = 0 

74 self.bkg_maker = 0 

75 self.safe_mask_masker = 0 

76 self.signal_table = 0 

77 model_best_joint = 0 

78 self.flux_points = 0 

79 self.model = 0 

80 self.model_best_stacked = 0 

81 self.obsCollection = 0 

82 self.stacked = 0 

83 self.npix_x = 0 

84 self.npix_y = 0 

85 self.ring_exclusion_mask = 0 

86 self.significance_map = 0 

87 self.excess_map = 0 

88 self.axis = 0 

89 self.exclusion_map = 0 

90 

91 def create_datastore(self, datastore_path, runid): 

92 datastore = DataStore.from_dir(datastore_path) 

93 obs_id = runid 

94 observations = datastore.get_observations(obs_id) 

95 

96 self.obs_id = obs_id 

97 self.datastore = datastore 

98 self.observations = observations 

99 self.table = [] 

100 for i in obs_id: 

101 self.table.append(datastore.obs(obs_id=i).events.table) 

102 

103 def plot_timee(self, id_=0, ax=None): 

104 """Plots an event rate time curve. 

105 

106 Parameters 

107 ---------- 

108 ax : `~matplotlib.axes.Axes` or None 

109 Axes 

110 

111 Returns 

112 ------- 

113 ax : `~matplotlib.axes.Axes` 

114 Axes 

115 i 

116 """ 

117 plt.figure(figsize=(16, 8)) 

118 for i in range(len(self.table)): 

119 ax = plt.gca() if ax is None else ax 

120 # Note the events are not necessarily in time order 

121 # print(self.table[i]) 

122 time = self.table[i]["TIME"] 

123 time = time - np.min(time) 

124 

125 ax.set_xlabel("Time (sec)") 

126 ax.set_ylabel("Counts") 

127 y, x_edges = np.histogram(time, bins=np.linspace(0, 2400, 100)) 

128 y = y[:-1] 

129 x_edges = x_edges[:-1] 

130 

131 xerr = np.diff(x_edges) / 2 

132 x = x_edges[:-1] + xerr 

133 yerr = np.sqrt(y) 

134 ax.errorbar(x=x, y=y, xerr=xerr, yerr=yerr, label="run" + str(self.obs_id[i])) 

135 

136 plt.legend( 

137 shadow=True, bbox_to_anchor=(1.05, 1), loc="upper left", borderaxespad=0, handlelength=1.5, fontsize=10 

138 ) 

139 prefix = (5 - len(str(id_))) * "0" 

140 plt.savefig("../plots/event_rate" + prefix + str(id_) + ".png") 

141 plt.show() 

142 

143 return x, y 

144 

145 def create_counts_scatter(self): 

146 time = self.table["TIME"] 

147 time = time - np.min(time) 

148 y, x_edges = np.histogram(time, bins=np.linspace(0, 1200, 20)) 

149 y = y[:-1] 

150 

151 self.count_scatter = y 

152 

153 def statisticals_parameters(self, runid): 

154 info_table = self.datastore.obs_table 

155 self.create_counts_scatter() 

156 

157 variation_max = ( 

158 np.abs((self.count_scatter.max() - self.count_scatter.mean()) / self.count_scatter.mean()) * 100 

159 ) 

160 variation_min = ( 

161 np.abs((self.count_scatter.min() - self.count_scatter.mean()) / self.count_scatter.mean()) * 100 

162 ) 

163 general_variation = self.count_scatter.std() / self.count_scatter.mean() * 100 

164 variability = self.count_scatter.std() / (np.sqrt(self.count_scatter.mean())) 

165 mean_trigger_rate = self.count_scatter.mean() / info_table["LIVETIME"][runid] 

166 zenith_angle = info_table["ZEN_PNT"][runid] 

167 livetime = info_table["LIVETIME"][runid] 

168 

169 print(f"Mean = {self.count_scatter.mean()}") 

170 print(f"Standard deviation = {self.count_scatter.std()}") 

171 print(f"Localized variation on the maximum = {variation_max} %") 

172 print(f"Localized variation on the maximum = {variation_min} %") 

173 print(f"General variation = {general_variation}") 

174 print(f"Variability = {variability}") 

175 print(f"Mean_trigger_rate = {mean_trigger_rate}") 

176 print(f"Zenith_angle = {zenith_angle}") 

177 print(f"Livetime = {livetime}") 

178 

179 self.info_table = info_table 

180 

181 return ( 

182 self.count_scatter.mean(), 

183 self.count_scatter.std(), 

184 variation_max, 

185 variation_min, 

186 general_variation, 

187 variability, 

188 mean_trigger_rate, 

189 zenith_angle, 

190 livetime, 

191 ) 

192 

193 def create_on_region(self, ra=83.633, dec=22.014, on_radius_angle=0.2, radius_excluded=0.35): 

194 # def create_on_exclusion_regions(self, ra=238.929, dec=11.190, on_radius_angle=0.14, radius_excluded=0.2): 

195 # def create_on_exclusion_regions(self, ra=187.706, dec=12.391, on_radius_angle=0.14, radius_excluded=0.2): 

196 

197 sourcePos = SkyCoord(ra=ra * u.deg, dec=dec * u.deg, frame="icrs") 

198 

199 on_radius = on_radius_angle * u.deg 

200 on_region = CircleSkyRegion(center=sourcePos, radius=on_radius) 

201 

202 # Create exclusion region 

203 reflectedbg_exclude_region = CircleSkyRegion(center=sourcePos, radius=radius_excluded * u.deg) 

204 ringbg_exclude_region = CircleSkyRegion(center=sourcePos, radius=on_radius_angle * 3 * u.deg) 

205 

206 self.sourcePos = sourcePos 

207 self.on_region = on_region 

208 self.reflectedbg_exclude_region = reflectedbg_exclude_region 

209 self.ringbg_exclude_region = ringbg_exclude_region 

210 

211 # Reflected Background Method 

212 

213 def create_exclusion_mask(self, plot=False, binsz=0.02, npix_x=400, npix_y=400): 

214 skydir = self.sourcePos.galactic 

215 

216 reflected_exclusion_mask = Map.create( 

217 npix=(npix_x, npix_y), binsz=binsz, skydir=skydir, proj="TAN", frame="icrs" 

218 ) 

219 mask = reflected_exclusion_mask.geom.region_mask([self.reflectedbg_exclude_region], inside=False) 

220 reflected_exclusion_mask.data = mask 

221 

222 self.skydir = skydir 

223 self.reflected_exclusion_mask = reflected_exclusion_mask 

224 

225 if plot: 

226 mask.plot() 

227 plt.show() 

228 

229 def reduction_chain( 

230 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 

231 ): 

232 # Run data reduction chain, e_=energy 

233 

234 e_reco = MapAxis.from_energy_bounds(e_reco_min, e_reco_max, e_reco_bin, unit="TeV", name="energy") 

235 e_true = MapAxis.from_energy_bounds(e_true_min, e_true_max, e_true_bin, unit="TeV", name="energy_true") 

236 dataset_empty = SpectrumDataset.create(e_reco=e_reco, e_true=e_true, region=self.on_region) 

237 self.dataset_empty = dataset_empty 

238 

239 def data_maker(self): 

240 dataset_maker = SpectrumDatasetMaker(containment_correction=False, selection=["counts", "exposure", "edisp"]) 

241 bkg_maker = ReflectedRegionsBackgroundMaker(exclusion_mask=self.reflected_exclusion_mask) 

242 safe_mask_masker = SafeMaskMaker(methods=["aeff-max"], aeff_percent=10) 

243 

244 self.dataset_maker = dataset_maker 

245 self.bkg_maker = bkg_maker 

246 self.safe_mask_masker = safe_mask_masker 

247 

248 def data_run(self): 

249 datasets = Datasets() 

250 

251 observations = self.datastore.get_observations(self.obs_id) 

252 

253 for obs_id, observation in zip(self.obs_id, observations): 

254 dataset = self.dataset_maker.run(self.dataset_empty.copy(name=str(obs_id)), observation) 

255 dataset_on_off = self.bkg_maker.run(dataset, observation) 

256 dataset_on_off = self.safe_mask_masker.run(dataset_on_off, observation) 

257 datasets.append(dataset_on_off) 

258 

259 self.datasets = datasets 

260 

261 def plot_off_region(self): 

262 plt.figure(figsize=(8, 8)) 

263 _, ax, _ = self.reflected_exclusion_mask.plot() 

264 self.on_region.to_pixel(ax.wcs).plot(ax=ax, edgecolor="k") 

265 plot_spectrum_datasets_off_regions(ax=ax, datasets=self.datasets) 

266 

267 plt.show() 

268 

269 def signal_info(self): 

270 info_table = self.datastore.obs_table 

271 print(info_table[0]) 

272 signal_table = self.datasets.info_table(cumulative=True) 

273 alpha = signal_table["alpha"] 

274 

275 excess = signal_table["excess"] 

276 uncertainty_excess = np.sqrt( 

277 signal_table["counts"] + ((1 / alpha) * signal_table["background"]) / ((1 / alpha) ** 2) 

278 ) 

279 

280 bg = signal_table["background"] 

281 uncertainty_bg = np.sqrt((1 / alpha) * signal_table["background"]) / (1 / alpha) 

282 

283 self.signal_table = signal_table 

284 

285 print(f"reflected_excess = {excess[-1]}") 

286 print(f"reflected_uncertainty_excess = {uncertainty_excess[-1]}") 

287 print(f"reflected_uncertainty_excess = {self.signal_table['livetime'].to('h')[-1]} hours") 

288 print(f"sqrt_sts = {self.signal_table['sqrt_ts'][-1]}") 

289 # print(f"reflected_bg = {bg}") 

290 # print(f"reflected_uncertainty_bg = {uncertainty_bg}") 

291 

292 return ( 

293 excess[-1], 

294 uncertainty_excess[-1], 

295 bg[-1], 

296 uncertainty_bg[-1], 

297 self.signal_table["sqrt_ts"][-1], 

298 self.signal_table["livetime"][-1], 

299 ) 

300 

301 def statistic_source_excess(self, plot=False): 

302 if plot: 

303 plt.plot( 

304 # self.signal_table["livetime"].to("h"), 

305 self.signal_table["name"], 

306 self.signal_table["excess"], 

307 marker="o", 

308 ls="none", 

309 ) 

310 plt.xlabel("Livetime [h]") 

311 plt.ylabel("Excess") 

312 

313 plt.show() 

314 

315 def statistic_source_ts(self, plot=False): 

316 if plot: 

317 plt.plot(self.signal_table["livetime"].to("h"), self.signal_table["sqrt_ts"], marker="o", ls="none") 

318 plt.xlabel("Livetime [h]") 

319 plt.ylabel("Sqrt(TS)") 

320 

321 plt.show() 

322 

323 def fit_spectrum(self): 

324 # Fit spectrum 

325 spectral_model = PowerLawSpectralModel( 

326 index=2, amplitude=2e-11 * u.Unit("cm-2 s-1 TeV-1"), reference=1 * u.TeV 

327 ) 

328 model = SkyModel(spectral_model=spectral_model, name="crab") 

329 

330 for self.dataset in self.datasets: 

331 self.dataset.models = model 

332 

333 fit_joint = Fit(self.datasets) 

334 result_joint = fit_joint.run() 

335 

336 # Make a copy here to compare it later 

337 model_best_joint = model.copy() 

338 

339 self.model = model 

340 self.model_best_joint = model_best_joint 

341 

342 def fit_quality(self, plot=False): 

343 if plot: 

344 ax_spectrum, ax_residuals = self.datasets[0].plot_fit() 

345 ax_spectrum.set_ylim(0.1, 40) 

346 plt.show() 

347 

348 def flux_points_plot(self, e_min=0.7, e_max=30, plot=False): 

349 # Compute Flux Points 

350 energy_edges = np.logspace(np.log10(e_min), np.log10(e_max), 11) * u.TeV 

351 

352 # Create an instance of the FluxPointsEstimator 

353 fpe = FluxPointsEstimator(energy_edges=energy_edges, source="crab") 

354 flux_points = fpe.run(datasets=self.datasets) 

355 

356 flux_points.table_formatted 

357 

358 # Plot the flux points and their likelihood profile with Treshold < 4 

359 if plot: 

360 plt.figure(figsize=(8, 5)) 

361 flux_points.table["is_ul"] = flux_points.table["ts"] < 4 

362 ax = flux_points.plot(energy_power=2, flux_unit="erg-1 cm-2 s-1", color="darkorange") 

363 flux_points.to_sed_type("e2dnde").plot_ts_profiles(ax=ax) 

364 plt.show() 

365 

366 self.flux_points = flux_points 

367 

368 def model_fit_plot(self, plot=False): 

369 # Final plot with the best fit model 

370 

371 flux_points_dataset = FluxPointsDataset(data=self.flux_points, models=self.model_best_joint) 

372 

373 if plot: 

374 flux_points_dataset.plot_fit() 

375 plt.show() 

376 

377 def flux_plot(self, plot=False, id_=0): 

378 dataset_stacked = Datasets(self.datasets).stack_reduce() 

379 

380 dataset_stacked.models = self.model 

381 stacked_fit = Fit([dataset_stacked]) 

382 result_stacked = stacked_fit.run() 

383 

384 # Make a copy to compare later 

385 model_best_stacked = self.model.copy() 

386 

387 if plot: 

388 plt.figure(figsize=(16, 8)) 

389 plot_kwargs = { 

390 "energy_range": [0.1, 30] * u.TeV, 

391 "energy_power": 2, 

392 "flux_unit": "erg-1 cm-2 s-1", 

393 } 

394 

395 # courb stacked model 

396 model_best_stacked.spectral_model.plot(**plot_kwargs, label="run" + str(self.obs_id)) 

397 model_best_stacked.spectral_model.plot_error(**plot_kwargs) 

398 

399 # courb reference 

400 create_crab_spectral_model("hess_pl").plot(**plot_kwargs, label="Crab reference") 

401 

402 plt.legend() 

403 prefix = (5 - len(str(id_))) * "0" 

404 plt.savefig("./Plots/flux_plot" + prefix + str(id_) + ".png") 

405 plt.show() 

406 

407 self.model_best_stacked = model_best_stacked 

408 

409 def flux_parameters(self): 

410 parameters_table = self.model_best_stacked.parameters.to_table() 

411 

412 spectral_index = parameters_table[0][1] 

413 uncertainty_spectral_index = parameters_table[0][6] 

414 amplitude = parameters_table[1][1] 

415 uncertainty_amplitude = parameters_table[1][6] 

416 

417 return spectral_index, uncertainty_spectral_index, amplitude, uncertainty_amplitude 

418 

419 """Ring Background Method""" 

420 

421 def calculate_acceptance_model( 

422 self, 

423 energyaxis_min=np.log10(0.1), 

424 energyaxis_max=np.log10(2.0), 

425 energy_nbin=5, 

426 offset_min=0.0, 

427 offset_max=5.0, 

428 offset_nbin=8, 

429 plot=False, 

430 ): 

431 obsCollection = self.datastore.get_observations(self.obs_id) 

432 energyAxisAcceptance = MapAxis.from_edges( 

433 np.logspace(energyaxis_min, energyaxis_max, energy_nbin), unit="TeV", name="energy", interp="log" 

434 ) 

435 offsetAxisAcceptance = MapAxis.from_edges( 

436 np.linspace(offset_min, offset_max, offset_nbin), unit="deg", name="offset", interp="lin" 

437 ) 

438 background = create_radial_acceptance_map( 

439 obsCollection, 

440 energyAxisAcceptance, 

441 offsetAxisAcceptance, 

442 exclude_regions=[self.ringbg_exclude_region], 

443 oversample_map=10, 

444 ) 

445 

446 if plot: 

447 background.peek() 

448 plt.show 

449 

450 hduBackground = background.to_table_hdu() 

451 hduBackground.writeto("background.fits", overwrite=True) 

452 

453 self.obsCollection = obsCollection 

454 

455 def add_acceptance_map(self): 

456 # Add acceptance map to observations 

457 

458 listIDObs = [] 

459 for obs in self.obsCollection: 

460 listIDObs.append(obs.obs_id) 

461 self.datastore.hdu_table.add_row( 

462 { 

463 "OBS_ID": listIDObs[-1], 

464 "HDU_TYPE": "bkg", 

465 "HDU_CLASS": "bkg_2d", 

466 "FILE_DIR": "", 

467 "FILE_NAME": "background.fits", 

468 "HDU_NAME": "BACKGROUND", 

469 } 

470 ) 

471 self.datastore.hdu_table = self.datastore.hdu_table.copy() 

472 self.datastore.hdu_table 

473 

474 self.obsCollection = self.datastore.get_observations(listIDObs) 

475 

476 def map_geometry( 

477 self, 

478 RA, 

479 DEC, 

480 plot=False, 

481 npix_x=200, 

482 npix_y=200, 

483 binsize=0.02, 

484 axis_min=np.log10(0.1), 

485 axis_max=np.log10(2.0), 

486 axis_nbin=5, 

487 ): 

488 sourcePos = SkyCoord(ra=RA * u.deg, dec=DEC * u.deg, frame="icrs") 

489 self.axis = MapAxis.from_edges( 

490 np.logspace(axis_min, axis_max, axis_nbin), unit="TeV", name="energy", interp="log" 

491 ) 

492 geom = WcsGeom.create(skydir=sourcePos, npix=(npix_x, npix_y), binsz=binsize, frame="icrs", axes=[self.axis]) 

493 geom 

494 

495 ring_exclusion_mask = geom.to_image().region_mask([self.ringbg_exclude_region], inside=False) 

496 self.exclusion_map = WcsNDMap(geom.to_image(), ring_exclusion_mask) 

497 

498 self.ring_exclusion_mask = ring_exclusion_mask 

499 

500 if plot: 

501 self.exclusion_map.plot() 

502 plt.show() 

503 

504 stacked = MapDataset.create(geom=geom) 

505 unstacked = Datasets() 

506 maker = MapDatasetMaker(selection=["counts", "background"]) 

507 # maker = MapDatasetMaker(selection=["counts"]) 

508 maker_safe_mask = SafeMaskMaker(methods=["offset-max"], offset_max=3.0 * u.deg) 

509 

510 for obs in self.obsCollection: 

511 cutout = stacked.cutout(obs.pointing_radec, width="10 deg") 

512 dataset = maker.run(cutout, obs) 

513 dataset = maker_safe_mask.run(dataset, obs) 

514 stacked.stack(dataset) 

515 unstacked.append(dataset) 

516 

517 self.stacked = stacked 

518 self.npix_x = npix_x 

519 self.npix_y = npix_y 

520 

521 def ring_data_estimation(self): 

522 ring_maker = RingBackgroundMaker(r_in="0.5 deg", width="0.3 deg", exclusion_mask=self.ring_exclusion_mask) 

523 

524 estimator = ExcessMapEstimator(0.2 * u.deg) 

525 lima_maps = estimator.run(self.stacked) 

526 

527 significance_map = lima_maps["sqrt_ts"] 

528 excess_map = lima_maps["excess"] 

529 

530 npix_x = int(self.npix_x / 2) 

531 npix_y = int(self.npix_y / 2) 

532 

533 ring_sqrt_ts = lima_maps["sqrt_ts"].data[0][npix_x][npix_y] 

534 

535 ring_excess = lima_maps["excess"].data[0][npix_x][npix_y] 

536 ring_bg = lima_maps["background"].data[0][npix_x][npix_y] 

537 

538 uncertainty_excess_ring = lima_maps["err"].data[0][npix_x][npix_y] 

539 # uncertainty_ring_bg = np.sqrt((1/self.alpha)*ring_bg)/(1/self.alpha) 

540 

541 self.significance_map = significance_map 

542 self.excess_map = excess_map 

543 

544 print(f"ring_excess = {ring_excess}") 

545 print(f"ring_sqrt_ts = {ring_sqrt_ts}") 

546 print(f"ring_bg = {ring_bg}") 

547 print(f"uncertainty_excess_ring = {uncertainty_excess_ring}") 

548 # print(f"uncertainty_ring_bg{uncertainty_ring_bg}") 

549 

550 return ring_excess, ring_sqrt_ts, ring_bg, uncertainty_excess_ring 

551 # uncertainty_ring_bg 

552 

553 def excess_significance_plot(self, directory, plot=False, id_=0): 

554 if plot: 

555 plt.figure(figsize=(16, 8)) 

556 ax1 = plt.subplot(121, projection=self.significance_map.geom.wcs) 

557 ax2 = plt.subplot(122, projection=self.excess_map.geom.wcs) 

558 

559 ax2.set_title("Significance map") 

560 self.significance_map.plot(ax=ax2, add_cbar=True) 

561 sources = find_peaks( 

562 self.significance_map.get_image_by_idx((0,)), 

563 threshold=5, 

564 min_distance="0.2 deg", 

565 ) 

566 print("Found sources") 

567 print(sources) 

568 f = open(directory + "../plots/results.txt", "w") 

569 f.write("Found sources") 

570 f.write(str(sources)) 

571 f.close() 

572 sources_2 = find_peaks( 

573 self.significance_map.get_image_by_idx((0,)), 

574 threshold=7, 

575 min_distance="0.2 deg", 

576 ) 

577 

578 now = dt.datetime.now() 

579 timestamp_str = now.strftime("%Y-%m-%d %H:%M:%S") 

580 ax1.text( 

581 0.02, 

582 0.98, 

583 timestamp_str, 

584 transform=ax1.transAxes, 

585 fontsize=11, 

586 fontweight="bold", 

587 va="top", 

588 ha="left", 

589 ) 

590 ax2.text( 

591 0.02, 

592 0.98, 

593 timestamp_str, 

594 transform=ax2.transAxes, 

595 fontsize=11, 

596 fontweight="bold", 

597 va="top", 

598 ha="left", 

599 ) 

600 if len(sources) > 0: 

601 ax2.scatter( 

602 sources["ra"], 

603 sources["dec"], 

604 transform=plt.gca().get_transform("icrs"), 

605 color="none", 

606 edgecolor="white", 

607 marker="o", 

608 s=300, 

609 lw=1.5, 

610 ) 

611 if len(sources_2) > 0: 

612 ax2.scatter( 

613 sources_2["ra"], 

614 sources_2["dec"], 

615 transform=plt.gca().get_transform("icrs"), 

616 color="none", 

617 edgecolor="blue", 

618 marker="o", 

619 s=300, 

620 lw=1.5, 

621 ) 

622 

623 ax1.set_title("Excess map") 

624 self.excess_map.plot(ax=ax1, add_cbar=True) 

625 

626 prefix = (5 - len(str(id_))) * "0" 

627 plt.savefig(directory + "../plots/sig_excess_plot" + prefix + str(id_) + ".png") 

628 plt.show() 

629 return sources 

630 

631 def off_distribution_plot(self, plot=False): 

632 # create a 2D mask for the images 

633 exclusion_map_ring = self.significance_map.geom.region_mask([self.ringbg_exclude_region], inside=False) 

634 significance_map_off = self.significance_map * exclusion_map_ring 

635 significance_all = self.significance_map.data[np.isfinite(self.significance_map.data)] 

636 significance_off = significance_map_off.data[np.isfinite(significance_map_off.data)] 

637 

638 bins = np.linspace( 

639 np.min(significance_all), 

640 np.max(significance_all), 

641 num=int(np.max(significance_all - np.min(significance_all)) * 3), 

642 ) 

643 mu, std = norm.fit(significance_off) 

644 

645 if plot: 

646 plt.hist( 

647 significance_all, 

648 density=True, 

649 alpha=0.5, 

650 color="red", 

651 label="all bins", 

652 bins=bins, 

653 ) 

654 

655 plt.hist( 

656 significance_off, 

657 density=True, 

658 alpha=0.5, 

659 color="blue", 

660 label="off bins", 

661 bins=bins, 

662 ) 

663 

664 # Now, fit the off distribution with a Gaussian 

665 

666 x = np.linspace(-8, 8, 50) 

667 p = norm.pdf(x, mu, std) 

668 plt.plot(x, p, lw=2, color="black") 

669 plt.legend() 

670 plt.xlabel("Significance") 

671 plt.yscale("log") 

672 plt.ylim(1e-5, 1) 

673 xmin, xmax = np.min(significance_all), np.max(significance_all) 

674 plt.xlim(xmin, xmax) 

675 

676 plt.show() 

677 print(mu, std) 

678 return mu, std 

679 

680 def thetaSquarePlot(self): 

681 # theta2Edge = np.linspace(0.0, 0.35**2, num=15) 

682 theta2Edge = np.linspace(0.0, 0.35**2, num=15) 

683 thetaEdge = np.sqrt(theta2Edge) 

684 theta2 = (theta2Edge[1:] + theta2Edge[:-1]) / 2.0 

685 count = np.zeros(theta2.shape) 

686 countBack = np.zeros(theta2.shape) 

687 alpha = np.zeros(theta2.shape) 

688 significance = np.zeros(theta2.shape) 

689 sb = np.zeros(theta2.shape) 

690 

691 for i in range(len(theta2)): 

692 on_region_theta2_ring = CircleAnnulusSkyRegion( 

693 center=self.sourcePos, inner_radius=thetaEdge[i] * u.deg, outer_radius=thetaEdge[i + 1] * u.deg 

694 ) 

695 on_region_theta2_circle = CircleSkyRegion(center=self.sourcePos, radius=thetaEdge[i + 1] * u.deg) 

696 dataset_maker_spectrum_significance_theta2 = SpectrumDatasetMaker( 

697 selection=["counts", "exposure", "edisp"] 

698 ) 

699 spectrum_dataset_empty_significance_theta2_ring = SpectrumDataset.create( 

700 e_reco=self.axis, region=on_region_theta2_ring 

701 ) 

702 spectrum_dataset_empty_significance_theta2_circle = SpectrumDataset.create( 

703 e_reco=self.axis, region=on_region_theta2_circle 

704 ) 

705 bkg_maker_spectrum_significance_theta2 = ReflectedRegionsBackgroundMaker( 

706 exclusion_mask=self.exclusion_map 

707 ) 

708 

709 countsRing = np.zeros(len(self.obsCollection)) 

710 countsOffRing = np.zeros(len(self.obsCollection)) 

711 alphaRing = np.zeros(len(self.obsCollection)) 

712 exposureRing = np.zeros(len(self.obsCollection)) 

713 countsCircle = np.zeros(len(self.obsCollection)) 

714 countsOffCircle = np.zeros(len(self.obsCollection)) 

715 alphaCircle = np.zeros(len(self.obsCollection)) 

716 exposureCircle = np.zeros(len(self.obsCollection)) 

717 for j, obs in enumerate(self.obsCollection): 

718 print(i) 

719 dataset_spectrum_significance_theta2_ring = dataset_maker_spectrum_significance_theta2.run( 

720 spectrum_dataset_empty_significance_theta2_ring.copy(name=f"obs-{obs.obs_id}"), obs 

721 )