test_python.py
3.39 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
#! /usr/bin/env python
import gammalib
import math
import os
# ================================== #
# Extract data points from FITS file #
# ================================== #
def extract_data(filename):
"""
Extract data from FITS file.
"""
# Set plot styles
ecolor = ['r', 'b', 'g']
styles = ['ro', 'bo', 'go']
# Open FITS file
fits = gammalib.GFits(filename)
# Loop over all table HDUs
k = 0
for extno in range(fits.size()):
if fits.at(extno).exttype() == gammalib.GFitsHDU.HT_ASCII_TABLE or \
fits.at(extno).exttype() == gammalib.GFitsHDU.HT_BIN_TABLE:
# Load spectrum
spectrum = gammalib.GMWLSpectrum()
spectrum.load(filename+"["+str(extno)+"]")
# Add spectral points
x = []
y = []
dy = []
for i in range(spectrum.size()):
xval = spectrum[i].energy().MeV()
yval = spectrum[i].flux()
yerr = spectrum[i].flux_err()
conv = xval * xval * 1.6021765e-6
x.append(xval)
y.append(yval * conv)
dy.append(yerr * conv)
plt.loglog(x, y, styles[k])
plt.errorbar(x, y, dy, ecolor=ecolor[k])
# Increment style counter
k = k + 1
# Add axes
plt.xlabel("Energy (MeV)")
plt.ylabel("Flux (erg/cm2/s)")
# Close FITS file
fits.close()
# ============ #
# Fit spectrum #
# ============ #
def fit_spectrum(filename, xmlname):
"""
Fit spectrum.
"""
# Load observations
obs = gammalib.GObservations()
comptel = gammalib.GMWLObservation(filename+"[COMPTEL]")
lat = gammalib.GMWLObservation(filename+"[LAT]")
hess = gammalib.GMWLObservation(filename+"[HESS]")
obs.append(comptel)
obs.append(lat)
obs.append(hess)
# Setup model
models = gammalib.GModels(xmlname)
obs.models(models)
# Perform optimization
opt = gammalib.GOptimizerLM()
opt.max_iter(1000)
obs.optimize(opt)
print(obs)
# Get model values
x = np.power(10., np.arange(-1., 8., 0.1))
y = []
d = gammalib.GSkyDir()
d.radec_deg(83.6331, 22.0145)
e = gammalib.GEnergy()
t = gammalib.GTime()
for energy in x:
e.MeV(energy)
yval = 0.0
for model in obs.models():
if model.type() == "PointSource":
yval += model.value(gammalib.GPhoton(d, e, t)) * \
energy * energy * 1.6021765e-6
y.append(yval)
plt.plot(x, y)
#==========================#
# Main routine entry point #
#==========================#
if __name__ == '__main__':
"""
Perform testing.
"""
# Dump result
print("")
print("****************************")
print("* Python interface testing *")
print("****************************")
# Only proceed if matplotlib is available
try:
# Import matplotlib
import numpy as np
import matplotlib.pyplot as plt
# Open plot
plt.figure(1)
plt.title('GammaLib fit of the Crab PWN')
# Extract data from FITS file
extract_data("data/crab_mwl.fits")
# Fit spectrum
fit_spectrum("data/crab_mwl.fits", "data/crab_mwl2.xml")
# Show spectrum
plt.show()
except ImportError:
print("Matplotlib is not (correctly) installed on your system.")