Coverage for src / lstautorta / rta_var.py: 0%
346 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-08-10 11:56 +0000
1#!/usr/bin/env python
3import argparse
4import logging
5import os
6import smtplib
7import ssl
8from argparse import Namespace
9from datetime import datetime
10from pathlib import Path
12import astropy.units as u
13import matplotlib.pyplot as plt
14import numpy as np
15import scipy
16from astropy.coordinates import SkyCoord
17from astropy.time import Time
18from gammapy.data import DataStore, Observations
19from gammapy.datasets import Datasets, FluxPointsDataset, SpectrumDataset
20from gammapy.estimators import FluxPointsEstimator, LightCurveEstimator
21from gammapy.makers import ReflectedRegionsBackgroundMaker, SpectrumDatasetMaker
22from gammapy.maps import MapAxis, RegionGeom, WcsGeom
23from gammapy.modeling import Fit
24from gammapy.modeling.models import LogParabolaSpectralModel, PowerLawSpectralModel, SkyModel
25from regions import CircleSkyRegion
26from scipy.stats import chi2
29def build_argparser():
30 parser = argparse.ArgumentParser(
31 description="For DL3 files of one observationnight for a given source, search if the flux is variable."
32 )
34 parser.add_argument(
35 "--dl3_folder", default=None, dest="dl3_folder", type=str, required=True, help="Path of the DL3 files to use."
36 )
37 parser.add_argument(
38 "--output_dir", default=None, dest="output_dir", type=str, required=True, help="Path where results are saved."
39 )
40 parser.add_argument(
41 "--source",
42 "-s",
43 default="???",
44 dest="source_name",
45 type=str,
46 required=False,
47 help="Source name of the target.",
48 )
49 parser.add_argument(
50 "-ra", default=None, dest="ra", type=float, required=True, help="RA coordinate of the target."
51 )
52 parser.add_argument(
53 "-dec", default=None, dest="dec", type=float, required=True, help="DEC coordinate of the target."
54 )
55 parser.add_argument(
56 "--date", "-d", default=None, dest="date", type=str, required=True, help="Date of the observation night."
57 )
58 parser.add_argument(
59 "--runlist",
60 "-rl",
61 default=None,
62 dest="runlist",
63 type=str,
64 required=False,
65 help="Run list to be analysed from txt file (--rl file.txt).",
66 )
67 parser.add_argument(
68 "--distance",
69 "-dis",
70 default=3,
71 dest="distance",
72 type=float,
73 required=False,
74 help="Max distance in degrees between the target position and the run pointing position "
75 "for the run selection, negative value means no selection using this parameter.",
76 )
78 return parser
81class SourceData:
82 def __init__(self, args: Namespace) -> None:
83 self.path = args.dl3_folder
84 self.date = args.date
85 self.runlist = args.runlist
86 self.source_position = SkyCoord(args.ra, args.dec, unit="deg")
87 self.distance = args.distance
88 self.obs_ids = None
89 self.datastore = None
90 self.obs_collection = None
91 self.required_irf = ["aeff", "edisp"]
93 self.compute()
95 def get_runs_database(self) -> list:
96 if self.runlist: # and self.distance < 0.:
97 obs_ids = np.loadtxt(self.runlist, unpack=True, dtype=int, ndmin=1)
99 elif self.distance > 0.0: # and not self.runlist :
100 obs_table = self.datastore.obs_table
101 obs_table = obs_table[obs_table["LIVETIME"].data > 5 * 60] # remove runs<5min
102 obs_ids = obs_table.select_sky_circle(self.source_position, self.distance * u.deg, inverted=False)[
103 "OBS_ID"
104 ].tolist()
105 # logging.info('Selected obs_ids : %s'%obs_ids)
107 else:
108 raise RuntimeError("Cannot make a run selection. Either (and only) runlist or distance>0 is needed.")
110 if len(obs_ids) == 0:
111 raise RuntimeError("No run selected with current input.")
112 # logging.info(obs_ids)
114 return obs_ids
116 def compute(self) -> None:
117 self.datastore = DataStore.from_dir(self.path)
118 self.obs_ids = self.get_runs_database()
119 self.obs_collection = self.datastore.get_observations(self.obs_ids, required_irf=self.required_irf)
121 def get_datastore(self) -> DataStore:
122 return self.datastore
124 def get_obs_collection(self) -> Observations:
125 if not self.obs_collection:
126 self.compute()
127 return self.obs_collection
129 def get_obs_ids(self) -> list:
130 if not self.obs_ids:
131 self.compute()
132 return self.obs_ids
135#########################
138class Significance:
139 def __init__(
140 self, source_position: SkyCoord, e_min_reco: float = 0.01, e_max_reco: float = 10, on_radius: float = 0.2
141 ) -> None:
142 self.source_position = source_position
143 self.e_min_reco = e_min_reco * u.TeV
144 self.e_max_reco = e_max_reco * u.TeV
145 self.on_radius = on_radius
146 self.obs_collection = None
147 self.initialize()
149 def initialize(self) -> None:
150 on_radius = self.on_radius * u.deg
151 on_region = CircleSkyRegion(center=self.source_position, radius=on_radius)
152 exclusion_radius = 0.35 * u.deg
153 exclude_region = CircleSkyRegion(center=self.source_position, radius=exclusion_radius)
154 energy_axis = MapAxis.from_edges(
155 np.logspace(np.log10(self.e_min_reco.to_value(u.TeV)), np.log10(self.e_max_reco.to_value(u.TeV)), 2),
156 unit="TeV",
157 name="energy",
158 interp="log",
159 )
161 geom_on_region = RegionGeom.create(region=on_region, axes=[energy_axis])
162 geom = WcsGeom.create(
163 skydir=self.source_position, npix=(200, 200), binsz=0.02, frame="icrs", axes=[energy_axis]
164 )
165 geom_image = geom.to_image()
166 exclusion_mask = ~geom_image.region_mask([exclude_region])
168 self.dataset_maker_spectrum_significance = SpectrumDatasetMaker(
169 selection=["counts", "exposure"], use_region_center=True
170 )
171 self.spectrum_dataset_empty_significance = SpectrumDataset.create(geom=geom_on_region)
172 self.bkg_maker_spectrum_significance = ReflectedRegionsBackgroundMaker(exclusion_mask=exclusion_mask)
174 def compute(self, obs_collection: Observations) -> None:
175 datasets_significance = Datasets()
177 for obs in obs_collection:
178 dataset_spectrum_significance = self.dataset_maker_spectrum_significance.run(
179 self.spectrum_dataset_empty_significance.copy(name=f"obs-{obs.obs_id}"),
180 obs,
181 )
182 dataset_on_off_spectrum_significance = self.bkg_maker_spectrum_significance.run(
183 observation=obs,
184 dataset=dataset_spectrum_significance,
185 )
186 datasets_significance.append(dataset_on_off_spectrum_significance)
188 self.info_table = datasets_significance.info_table(cumulative=True)
190 def get_significance(self, obs_collection: Observations) -> float:
191 self.compute(obs_collection)
192 if self.info_table is not None:
193 significance = self.info_table["sqrt_ts"][-1]
194 if significance == np.nan:
195 return 0.0
196 report.add("Significance (Li&Ma) = %.2f" % significance)
197 return significance
198 return 0.0
200 def is_significant(self, obs_collection: Observations) -> bool:
201 bool_sign = self.get_significance(obs_collection) > 5.0
202 report.add("Is significant (>5σ) : %s" % bool_sign)
203 return bool_sign
206#########################
209class Flux:
210 def __init__(self, source_position: SkyCoord, e_min_reco: float = 0.2, e_max_reco: float = 10) -> None:
211 self.source_position = source_position
212 self.e_min_reco = e_min_reco * u.TeV
213 self.e_max_reco = e_max_reco * u.TeV
214 self.datasets_spectrum_joint_with_model = None
215 self.model = None
216 self.selected_model = None
217 self.energy_edges_flux_point = None
218 self.required_irf = ["aeff", "edisp"]
220 def initialize(self):
221 n_bin_per_decade_e_reco = 10
222 e_min_true = 0.005 * u.TeV
223 e_max_true = 40.0 * u.TeV
224 n_bin_per_decade_e_true = 20
226 e_min_flux_point = self.e_min_reco
227 e_max_flux_point = self.e_max_reco
228 n_bin_per_decade_flux_point = 10
230 n_bin_energy_reco = int(
231 (np.log10(self.e_max_reco.to_value(u.TeV)) - np.log10(self.e_min_reco.to_value(u.TeV)))
232 * n_bin_per_decade_e_reco
233 )
234 energy_axis_reco = MapAxis.from_edges(
235 np.logspace(
236 np.log10(self.e_min_reco.to_value(u.TeV)),
237 np.log10(self.e_max_reco.to_value(u.TeV)),
238 n_bin_energy_reco,
239 ),
240 unit="TeV",
241 name="energy",
242 interp="log",
243 )
245 n_bin_energy_true = int(
246 (np.log10(e_max_true.to_value(u.TeV)) - np.log10(e_min_true.to_value(u.TeV))) * n_bin_per_decade_e_true
247 )
248 energy_axis_true = MapAxis.from_edges(
249 np.logspace(
250 np.log10(e_min_true.to_value(u.TeV)), np.log10(e_max_true.to_value(u.TeV)), n_bin_energy_true
251 ),
252 unit="TeV",
253 name="energy_true",
254 interp="log",
255 )
257 on_radius = 0.2 * u.deg
258 on_region = CircleSkyRegion(center=self.source_position, radius=on_radius)
259 geom_on_region = RegionGeom.create(region=on_region, axes=[energy_axis_reco])
261 geom = WcsGeom.create(
262 skydir=self.source_position, npix=(200, 200), binsz=0.02, frame="icrs", axes=[energy_axis_reco]
263 )
264 geom_image = geom.to_image()
266 exclusion_radius = 0.35 * u.deg
267 exclude_region = CircleSkyRegion(center=self.source_position, radius=exclusion_radius)
268 exclusion_mask = ~geom_image.region_mask([exclude_region])
270 n_bin_energy_flux_point = int(
271 (np.log10(e_max_flux_point.to_value(u.TeV)) - np.log10(e_min_flux_point.to_value(u.TeV)))
272 * n_bin_per_decade_flux_point
273 )
274 energy_edges_flux_point = (
275 np.logspace(
276 np.log10(e_min_flux_point.to_value(u.TeV)),
277 np.log10(e_max_flux_point.to_value(u.TeV)),
278 n_bin_energy_flux_point,
279 )
280 * u.TeV
281 )
282 return geom_on_region, energy_axis_true, exclusion_mask, energy_axis_reco, energy_edges_flux_point
284 def get_datasets(self, obs_ids: list, obs_collection: Observations):
285 geom_on_region, energy_axis_true, exclusion_mask, energy_axis_reco, energy_edges_flux_point = (
286 self.initialize()
287 )
288 datasets_spectrum_joint = Datasets()
289 dataset_maker_spectrum = SpectrumDatasetMaker(
290 selection=["counts", "exposure", "edisp"], use_region_center=True
291 )
292 spectrum_dataset_empty = SpectrumDataset.create(geom=geom_on_region, energy_axis_true=energy_axis_true)
293 bkg_maker_spectrum = ReflectedRegionsBackgroundMaker(exclusion_mask=exclusion_mask)
295 for i, obs in enumerate(obs_collection):
296 dataset_spectrum = dataset_maker_spectrum.run(spectrum_dataset_empty.copy(name=f"obs-{i}".format(i)), obs)
297 dataset_on_off_spectrum = bkg_maker_spectrum.run(observation=obs, dataset=dataset_spectrum)
298 datasets_spectrum_joint.append(dataset_on_off_spectrum)
300 return datasets_spectrum_joint, energy_axis_reco, energy_edges_flux_point
302 def compute_models(self, energy_axis_reco):
303 reference_energy = energy_axis_reco.center[int(energy_axis_reco.nbin / 2)]
304 spectral_model = {
305 "power_law": PowerLawSpectralModel(
306 amplitude=1e-12 * u.Unit("cm-2 s-1 TeV-1"), index=2, reference=reference_energy
307 ),
308 "log_parabola": LogParabolaSpectralModel(
309 amplitude=1e-12 * u.Unit("cm-2 s-1 TeV-1"), alpha=2, beta=0, reference=reference_energy
310 ),
311 }
313 model = {}
314 for k in spectral_model:
315 model[k] = SkyModel(spectral_model=spectral_model[k], name=k)
317 return model, spectral_model
319 def select_model(self, khi2_flux) -> None:
320 """Select a model fit for the flux if pvalue >5 sigma with a power law as H0 hypothesis."""
321 challenger_model = min(khi2_flux, key=khi2_flux.get)
322 LR_statistic = khi2_flux["power_law"] - khi2_flux[challenger_model]
323 dof = 1
324 p_val = chi2.sf(LR_statistic, dof)
325 significance_challenger = np.sqrt(2) * scipy.special.erfinv(1 - p_val)
326 if significance_challenger > 5.0:
327 self.selected_model = challenger_model
328 else:
329 self.selected_model = "power_law"
331 def get_time_unix_from_lst_epoch(self, time):
332 if self.version == "lstchain":
333 return time
334 LST_EPOCH = Time("2018-10-01T00:00:00", scale="utc")
335 time_lst_unix = Time(LST_EPOCH, format="unix", scale="utc").value
336 if type(time) is list:
337 return [t + time_lst_unix for t in time]
338 return time + time_lst_unix
340 def get_dataset_joint(self, obs_collection: Observations, obs_ids: list) -> None:
341 datasets_spectrum_joint, energy_axis_reco, energy_edges_flux_point = self.get_datasets(
342 obs_ids, obs_collection
343 )
344 model, spectral_model = self.compute_models(energy_axis_reco)
345 datasets_spectrum_joint_with_model = {}
346 fit_joint = {}
347 result_fit_joint = {}
348 khi2_flux = {}
350 for k in spectral_model.keys():
351 datasets_spectrum_joint_with_model[k] = datasets_spectrum_joint.copy()
352 datasets_spectrum_joint_with_model[k].models = [model[k]]
353 fit_joint[k] = Fit()
354 result_fit_joint[k] = fit_joint[k].run(datasets=datasets_spectrum_joint_with_model[k])
355 khi2_flux[k] = result_fit_joint[k].total_stat
357 for k in spectral_model.keys():
358 if not result_fit_joint[k].success:
359 logging.warning("Fit of " + k + " did not converge.")
360 else:
361 self.select_model(khi2_flux)
362 self.datasets_spectrum_joint_with_model = datasets_spectrum_joint_with_model
363 self.model = model
364 self.energy_edges_flux_point = energy_edges_flux_point
365 return datasets_spectrum_joint_with_model
367 def get_plot(self, data: SourceData, output_file: Path):
368 self.required_irf = data.required_irf
369 obs_collection = data.get_obs_collection()
370 obs_ids = data.get_obs_ids()
371 self.get_dataset_joint(obs_collection, obs_ids)
372 fpe = FluxPointsEstimator(
373 energy_edges=self.energy_edges_flux_point,
374 source=self.selected_model,
375 selection_optional="all",
376 n_sigma_ul=3,
377 reoptimize=False,
378 )
379 flux_points = fpe.run(datasets=self.datasets_spectrum_joint_with_model[self.selected_model])
380 flux_point_dataset_loop = FluxPointsDataset(data=flux_points, models=self.model[self.selected_model])
381 flux_point_dataset_loop.plot_spectrum()
382 plt.legend([f"{self.selected_model} fit", "Stat. err.", "Flux"], loc=0, frameon=True, prop={"size": 12})
383 plt.tight_layout()
384 if output_file:
385 plt.savefig(output_file)
386 # plt.show()
387 plt.close()
390#######################
393class LightCurve:
394 def __init__(self, flux: Flux, obs_collection: Observations, obs_ids: list, output_path: Path = None) -> None:
395 self.output_path = output_path
396 self.obs_collection = obs_collection
397 self.obs_ids = obs_ids
398 self.flux = flux
399 self.energy_range = [0.01 * u.TeV, 10.0 * u.TeV]
400 self.name = "light_curve"
401 self.flux_points = None
402 self.flux_points_err = None
403 self.time = None
404 self.time_plot = None
405 self.tbin_l_plot = None
406 self.tbin_r_plot = None
407 self.up_lim = None
409 self.lc = self.compute()
411 def convert_mjd_unix(self, time):
412 return Time(time, format="mjd").to_value("unix")
414 def set_time_lc(self, t_init, time):
415 return self.convert_mjd_unix(time) - self.convert_mjd_unix(t_init)
417 def compute(self) -> None: # FluxPoints:
418 if self.flux.datasets_spectrum_joint_with_model is None:
419 datasets_spectrum_joint_with_model = self.flux.get_dataset_joint(self.obs_collection, self.obs_ids)
421 else:
422 datasets_spectrum_joint_with_model = self.flux.datasets_spectrum_joint_with_model
424 selected_model = self.flux.selected_model
426 light_curve_maker = LightCurveEstimator(
427 energy_edges=self.energy_range, reoptimize=False, n_sigma_ul=3, selection_optional="all"
428 )
429 lc = light_curve_maker.run(datasets_spectrum_joint_with_model[selected_model])
431 return lc
433 def get_plot(self) -> None:
434 plt.figure(figsize=(6, 4))
435 ax = self.lc.plot(sed_type="flux")
436 ax.tick_params(axis="x", labelsize=12)
437 ax.get_legend().remove()
439 plt.tight_layout()
440 if self.output_path:
441 plt.savefig(Path(self.output_path))
442 # plt.show()
443 plt.close()
445 def lc_extraction(self, lc) -> None:
446 lc_table = lc.to_table(sed_type="flux", format="lightcurve")
447 # len_table=len(lc_table)
448 lc_table = lc_table[~np.isnan(np.concatenate(lc_table["flux"]))] # sort nan values
449 # logging.info('Removed %i nan values from flux.'%(len_table-len(lc_table)))
451 self.flux_points = np.concatenate(lc_table["flux"].data)
452 flux_err = np.concatenate(lc_table["flux_err"].data)
453 self.up_lim = np.concatenate(lc_table["is_ul"].data)
455 # error of upper limit points is positive error
456 flux_points_errp = np.concatenate(lc_table["flux_errp"].data)
457 self.flux_points_err = [
458 flux_points_errp[i] if self.up_lim[i] == True else flux_err[i] for i in range(len(flux_err))
459 ]
461 tmin = lc_table["time_min"].data
462 tmax = lc_table["time_max"].data
463 self.time = (tmin + tmax) / 2
465 self.time_plot = self.set_time_lc(tmin[0], self.time)
466 self.tbin_l_plot = self.set_time_lc(tmin, self.time)
467 self.tbin_r_plot = self.set_time_lc(self.time, tmax)
469 def get_plot_var_fit(self, par: list, p_value: float, output_file: str) -> None:
470 def constant(x, a):
471 return a
473 fig, ax = plt.subplots(figsize=(6, 4))
474 ax.errorbar(
475 self.time_plot,
476 self.flux_points,
477 uplims=self.up_lim,
478 yerr=self.flux_points_err,
479 fmt="o",
480 label="flux error",
481 )
482 ax.errorbar(
483 self.time_plot,
484 self.flux_points,
485 uplims=self.up_lim,
486 xerr=np.array(list(zip(self.tbin_l_plot, self.tbin_r_plot))).T,
487 fmt="none",
488 capsize=4,
489 label="run",
490 color="darkblue",
491 )
492 ax.axhline(
493 constant(self.time_plot, par),
494 label=f"Fit with constant\np-value={np.format_float_scientific(p_value, precision=1)}",
495 color="orange",
496 )
497 ax.set_ylabel("Integrated flux (cm-2 s-1)")
498 ax.set_xlabel("Time (s)")
499 ax.legend(prop={"size": 12})
500 fig.tight_layout()
501 if output_file:
502 fig.savefig(Path(output_file))
503 report.add("saving variability lc in %s" % output_file)
504 # plt.show()
505 plt.close()
507 def variability(self) -> float:
508 self.lc_extraction(self.lc)
510 if len(self.flux_points) <= 1:
511 logging.error("Not enought lightcurve points for fitting.")
512 return np.nan, np.nan
514 def constant(x, a):
515 return a
517 par, cov = scipy.optimize.curve_fit(
518 constant, self.time, self.flux_points, sigma=self.flux_points_err, absolute_sigma=True
519 )
521 syst_lst_source_indep = 0.06 # crab paper 2023
522 syst_err = syst_lst_source_indep * np.array(self.flux_points)
523 khi2 = np.sum(
524 np.square(constant(self.time, par[0]) - self.flux_points)
525 / (np.square(self.flux_points_err) + np.square(syst_err))
526 )
528 ndof = len(self.time) - len(par)
529 p_value = chi2.sf(khi2, ndof)
531 if self.output_path:
532 self.get_plot_var_fit(par, p_value, self.output_path)
533 report.add("p-value = %f" % p_value)
535 return p_value
537 def is_variable(self) -> bool:
538 bool_pvalue = self.variability() < 5 * 10 ** (-7) # 5sigma
539 report.add("Is variable (<5σ) : %s" % bool_pvalue)
541 return bool_pvalue
544#######################
547class report_manager:
548 def __init__(self, data: SourceData, args: Namespace):
549 self.date = data.date
550 self.txt = (
551 "creation date : " + datetime.today().strftime("20%y-%m-%d %H:%M") + "\n\n"
552 "observation night : " + "%s-%s-%s" % (self.date[:4], self.date[4:6], self.date[6:]) + "\n\n"
553 "source : " + args.source_name + "\n"
554 "ra : " + str(args.ra) + "\n"
555 "dec : " + str(args.dec) + "\n"
556 "Max distance : " + str(args.distance) + "°\n"
557 "obs_ids : " + str(data.obs_ids) + "\n"
558 )
560 def add(self, txtline) -> None:
561 self.txt += txtline
562 self.txt += "\n"
564 def print(self) -> None:
565 print(self.txt)
567 def save(self, filename) -> None:
568 file = open(filename, "w")
569 file.write(self.txt)
570 file.close
572 def send(self) -> None:
573 smtp_server = "smtp.cnrs.fr"
574 port = 587 # For starttls
576 # sender_login = "login"
577 # sender_email = "mail"
578 # password = 'password'
579 # receiver_email = 'mail'
581 sender_login = "cyann.buisson@ods.services"
582 sender_email = "cyann.plard@cnrs.fr"
583 password = input("password")
584 receiver_email = "cyann.plard@lapp.in2p3.fr"
586 context = ssl.create_default_context()
588 theme = "RTA report of %s night" % self.date
589 msg = f'From: {sender_email}\r\nTo: {receiver_email}\r\nContent-Type: text/plain; charset="utf-8"\r\nSubject: {theme}\r\n\r\n'
590 msg += self.txt
592 try:
593 server = smtplib.SMTP(smtp_server, port)
594 server.ehlo() # Can be omitted
595 server.starttls(context=context) # Secure the connection
596 server.ehlo() # Can be omitted
597 server.login(sender_login, password)
598 server.sendmail(sender_email, receiver_email, msg.encode("utf8"))
600 except Exception as e:
601 print(e)
602 finally:
603 server.quit()
606#######################
609def main():
610 parser = build_argparser()
611 args = parser.parse_args()
612 data = SourceData(args)
613 global report
614 report = report_manager(data, args)
616 if Significance(data.source_position).is_significant(data.get_obs_collection()):
617 # if (0 > -1):
618 # Flux(data.source_position).get_plot(data, os.path.join(args.output_dir,'%s_flux.pdf'%args.date))
619 # LightCurve(Flux(data.source_position), data.obs_collection, data.obs_ids, os.path.join(args.output_dir,'%s_lc.pdf'%args.date)).get_plot()
620 LightCurve(
621 Flux(data.source_position),
622 data.obs_collection,
623 data.obs_ids,
624 os.path.join(args.output_dir, "%s_lc_variability.pdf" % args.date),
625 ).is_variable()
627 report.print()
630if __name__ == "__main__":
631 main()
633# python rta_var.py --dl3_folder /home/plard/Documents/code/gamma/data/DL3_RTA/20240111/16342/reco/dl3/ --source crab -ra 83.63 -dec 22.01 --distance 3 --date 20240111 --output_dir /home/plard/Documents/code/gamma/figures-results/tests/RTA/