|
1
|
|
|
2
|
import gammalib
|
|
3
|
import ctools as ct
|
|
4
|
import numpy as np
|
|
5
|
import matplotlib as mp
|
|
6
|
import matplotlib.pyplot as plt
|
|
7
|
from astropy.io import fits
|
|
8
|
|
|
9
|
|
|
10
|
def _scatter(array1, array2, xtitle="", ytitle="", maintitle=""):
|
|
11
|
plt.scatter(array1, array2, marker=".")
|
|
12
|
plt.ylabel(ytitle)
|
|
13
|
plt.xlabel(xtitle)
|
|
14
|
plt.show()
|
|
15
|
|
|
16
|
def _plot_energy(filename):
|
|
17
|
events = gammalib.GCTAEventList(filename)
|
|
18
|
event_num = [events[i].event_id() for i in range(events.size())]
|
|
19
|
energy_TeV = [events[i].energy().TeV() for i in range(events.size())]
|
|
20
|
_scatter(event_num, energy_TeV, "Event Number", "Energy [TeV]")
|
|
21
|
|
|
22
|
def _plot_phase(filename):
|
|
23
|
events = gammalib.GCTAEventList(filename)
|
|
24
|
hdu = fits.open("phased_events.fits")
|
|
25
|
event_num = [events[i].event_id() for i in range(events.size())]
|
|
26
|
times = [events[i].time().mjd() for i in range(events.size())]
|
|
27
|
phases = [ hdu]
|
|
28
|
phases = hdu[1].data.field("PHASE")
|
|
29
|
_scatter(times, phases, "MJD", "Phase")
|
|
30
|
|
|
31
|
def _plot_energy_time(filename):
|
|
32
|
events = gammalib.GCTAEventList(filename)
|
|
33
|
energy_TeV = [events[i].energy().TeV() for i in range(events.size())]
|
|
34
|
times = [events[i].time().mjd() for i in range(events.size())]
|
|
35
|
_scatter(times, energy_TeV, "MJD", "Energy [TeV]")
|
|
36
|
|
|
37
|
def _plot_phase_distribution(filename):
|
|
38
|
hdu = fits.open(filename)
|
|
39
|
phases = hdu[1].data.field("PHASE")
|
|
40
|
plt.hist(phases, 10)
|
|
41
|
plt.xlim(0.,1.)
|
|
42
|
plt.ylabel("Counts")
|
|
43
|
plt.xlabel("Phase")
|
|
44
|
plt.show()
|
|
45
|
|
|
46
|
|
|
47
|
obs = ct.ctobssim()
|
|
48
|
obs["inobs"] = ""
|
|
49
|
obs["inmodel"] = "test_model.xml"
|
|
50
|
obs["outevents"] = "events.fits"
|
|
51
|
obs["caldb"] = "prod2"
|
|
52
|
obs["irf"] = "South_0.5h"
|
|
53
|
obs["ra"] = 83.63
|
|
54
|
obs["dec"] = 22.01
|
|
55
|
obs["rad"] = 2.5
|
|
56
|
obs["tmin"] = 0
|
|
57
|
obs["tmax"] = 1800
|
|
58
|
obs["emin"] = 0.1
|
|
59
|
obs["emax"] = 100.0
|
|
60
|
obs.execute()
|
|
61
|
|
|
62
|
|
|
63
|
phase = ct.ctphase()
|
|
64
|
phase["inobs"] = "events.fits"
|
|
65
|
phase["outobs"] = "phased_events.fits"
|
|
66
|
phase["inmodel"] = "test_model.xml"
|
|
67
|
phase["source"] = "Crab"
|
|
68
|
phase["mjd"] = 51544.5
|
|
69
|
phase["p0"] = 0.0
|
|
70
|
phase["f0"] = 0.2345
|
|
71
|
phase["f1"] = -1.0e-10
|
|
72
|
phase["f2"] = -1.0e-13
|
|
73
|
phase.logFileOpen()
|
|
74
|
phase.execute()
|
|
75
|
|
|
76
|
|
|
77
|
|
|
78
|
|
|
79
|
|
|
80
|
|
|
81
|
_plot_phase("phased_events.fits")
|
|
82
|
|
|
83
|
_plot_phase_distribution("phased_events.fits")
|
|
84
|
|
|
85
|
|
|
86
|
|