|
1
|
from gammalib import *
|
|
2
|
from ctools import *
|
|
3
|
import os
|
|
4
|
import sys
|
|
5
|
import pyfits
|
|
6
|
|
|
7
|
def fits_header(fname,*args,**kwargs):
|
|
8
|
hdu = 1
|
|
9
|
for name,value in kwargs.items():
|
|
10
|
if name == "hdu":
|
|
11
|
hdu = value
|
|
12
|
fdata = pyfits.open(fname)
|
|
13
|
arr = []
|
|
14
|
for arg in args:
|
|
15
|
if isinstance(arg,str):
|
|
16
|
if arg in fdata[hdu].header:
|
|
17
|
arr.append(fdata[hdu].header[arg])
|
|
18
|
else:
|
|
19
|
arr.append(float('NaN'))
|
|
20
|
print fname+": Argument",arg,"not found in fits header (append NaN to array)"
|
|
21
|
|
|
22
|
elif isinstance(arg,list):
|
|
23
|
for argi in arg:
|
|
24
|
if argi in fdata[hdu].header:
|
|
25
|
arr.append(fdata[hdu].header[argi])
|
|
26
|
else:
|
|
27
|
arr.append(float('NaN'))
|
|
28
|
print fname+": Argument",argi,"not found in fits header (append NaN to array)"
|
|
29
|
else:
|
|
30
|
print "Wrong argument type in function fits_header(fname,*args,**kwargs): encountered argument",type(arg)
|
|
31
|
|
|
32
|
fdata.close()
|
|
33
|
if len(arr)==1:
|
|
34
|
return arr[0]
|
|
35
|
else:
|
|
36
|
return arr
|
|
37
|
|
|
38
|
def get_column_from_table(filename, column=0):
|
|
39
|
result = []
|
|
40
|
file = open(filename)
|
|
41
|
for line in file.readlines():
|
|
42
|
if len(line) <= 1:
|
|
43
|
continue
|
|
44
|
if line[0] == '#':
|
|
45
|
continue
|
|
46
|
if len(line.split()) > column:
|
|
47
|
result.append(line.split()[column])
|
|
48
|
else:
|
|
49
|
result.append('')
|
|
50
|
file.close()
|
|
51
|
return result
|
|
52
|
|
|
53
|
class base(object):
|
|
54
|
def __init__(self):
|
|
55
|
self.classname = self.__class__.__name__
|
|
56
|
self.errorcolor = "\033[31m"
|
|
57
|
self.infocolor = "\033[34m"
|
|
58
|
self.warningcolor = "\033[33m"
|
|
59
|
self.successcolor = "\033[32m"
|
|
60
|
self.endcolor = "\033[0m"
|
|
61
|
|
|
62
|
def error(self,message, functionname = ""):
|
|
63
|
printstring = ""
|
|
64
|
if functionname == "":
|
|
65
|
printstring = "\n"+self.errorcolor+"*** Error ["+self.classname+"]: "+message+" ***\n"+self.endcolor
|
|
66
|
else:
|
|
67
|
printstring = "\n"+self.errorcolor+"*** Error ["+self.classname+"::"+functionname+"]: "+message+" ***\n"+self.endcolor
|
|
68
|
sys.exit(printstring)
|
|
69
|
|
|
70
|
def info(self,message,newline=True):
|
|
71
|
printstring = self.infocolor+"["+self.classname+"]: "+message+self.endcolor
|
|
72
|
if newline:
|
|
73
|
print printstring
|
|
74
|
else:
|
|
75
|
print self.infocolor+message+self.endcolor,
|
|
76
|
sys.stdout.flush()
|
|
77
|
|
|
78
|
def warning(self,message,functionname = ""):
|
|
79
|
printstring = ""
|
|
80
|
if functionname == "":
|
|
81
|
printstring = self.warningcolor+"["+self.classname+"] Warning: "+message+self.endcolor
|
|
82
|
else:
|
|
83
|
printstring = self.warningcolor+"["+self.classname+"::"+functionname+"] Warning: "+message+self.endcolor
|
|
84
|
print printstring
|
|
85
|
def success(self,message):
|
|
86
|
printstring = self.successcolor+"["+self.classname+"]: "+message+self.endcolor
|
|
87
|
print printstring
|
|
88
|
|
|
89
|
class hess_obslist(base):
|
|
90
|
|
|
91
|
def __init__(self,runlist):
|
|
92
|
super(hess_obslist,self).__init__()
|
|
93
|
self.runlist = []
|
|
94
|
if isinstance(runlist,str):
|
|
95
|
self.runlistfile = os.path.abspath(runlist)
|
|
96
|
if not os.path.isfile(self.runlistfile):
|
|
97
|
self.error("Runlistfile "+self.runlistfile+" does not exist",sys._getframe().f_code.co_name)
|
|
98
|
self.runlist = get_column_from_table(self.runlistfile,0)
|
|
99
|
elif isinstance(runlist,list):
|
|
100
|
self.runlist = runlist
|
|
101
|
else:
|
|
102
|
self.error("Given runlist is neither list nor file",sys._getframe().f_code.co_name)
|
|
103
|
self.obs = GObservations()
|
|
104
|
self.emin = 0.1
|
|
105
|
self.emax = 100
|
|
106
|
self.runlist_emin = 1e3
|
|
107
|
self.runlist_emax = 1e-3
|
|
108
|
self.run_FoV = 2.5
|
|
109
|
self.datapath = os.environ["HDATA_PATH"]
|
|
110
|
self.bgpath = os.environ["HBG_PATH"]
|
|
111
|
self.verbose = True
|
|
112
|
self.nodes = 1
|
|
113
|
self.alt_bounds = [0,20,23,27,30,33,37,40,44,49,53,58,64,72,90]
|
|
114
|
self.az_bounds = [-180,-90,90,180]
|
|
115
|
|
|
116
|
def _eventsfile(self,runnr):
|
|
117
|
return self.datapath+"/events/events_"+runnr+".fits"
|
|
118
|
def _psffile(self,runnr):
|
|
119
|
return self.datapath+"/psf/psf_"+runnr+".fits"
|
|
120
|
def _Aefffile(self,runnr):
|
|
121
|
return self.datapath+"/Aeff/Aeff_"+runnr+".fits"
|
|
122
|
|
|
123
|
def _select(self,evlist,emin,emax):
|
|
124
|
|
|
125
|
[ra,dec] = fits_header(evlist,"RA_PNT","DEC_PNT")
|
|
126
|
|
|
127
|
expr = evlist+"[EVENTS][ENERGY >= "+str(emin)+" && ENERGY <= "+str(emax)+" && ANGSEP("+str(ra)+","+str(dec)+",RA,DEC) <= "+str(self.run_FoV)+"]"
|
|
128
|
fits = GFits(expr)
|
|
129
|
fits[1].card("NDSKEYS",3,"")
|
|
130
|
fits[1].card("DSTYP1","TIME","")
|
|
131
|
fits[1].card("DSUNI1",'s',"")
|
|
132
|
fits[1].card("DSVAL1","TABLE","")
|
|
133
|
fits[1].card("DSREF1",":GTI","")
|
|
134
|
fits[1].card("DSTYP2","POS(RA,DEC)","")
|
|
135
|
fits[1].card("DSUNI2",'deg',"")
|
|
136
|
fits[1].card("DSVAL2","CIRCLE("+str(ra)+","+str(dec)+","+str(self.run_FoV)+")","")
|
|
137
|
fits[1].card("DSTYP3","ENERGY","")
|
|
138
|
fits[1].card("DSUNI3",'TeV',"")
|
|
139
|
fits[1].card("DSVAL3",str(emin)+":"+str(emax),"")
|
|
140
|
|
|
141
|
obs = GCTAObservation("HESS")
|
|
142
|
obs.read(fits)
|
|
143
|
return obs
|
|
144
|
|
|
145
|
|
|
146
|
def _findbins(self,alt,az):
|
|
147
|
altbin=-1
|
|
148
|
azbin = -1
|
|
149
|
for i in range(len(self.alt_bounds)-1):
|
|
150
|
if alt >= self.alt_bounds[i] and alt <= self.alt_bounds[i+1]:
|
|
151
|
altbin=i
|
|
152
|
break
|
|
153
|
if az < -90 or az > 90:
|
|
154
|
azbin = 1
|
|
155
|
else:
|
|
156
|
azbin = 0
|
|
157
|
return altbin,azbin
|
|
158
|
|
|
159
|
def hess_inst_background(self,run,emin=0.1,emax=80):
|
|
160
|
spec = 0
|
|
161
|
if self.nodes <= 1:
|
|
162
|
spec = GModelSpectralConst()
|
|
163
|
spec["Value"].min(0.1)
|
|
164
|
spec["Value"].max(10)
|
|
165
|
if self.nodes <1:
|
|
166
|
spec["Value"].fix()
|
|
167
|
else:
|
|
168
|
spec["Value"].free()
|
|
169
|
elif self.nodes == 2:
|
|
170
|
e= GEnergy(1.0,"TeV")
|
|
171
|
spec = GModelSpectralPlaw(1.0,0.0,e)
|
|
172
|
spec[0].free()
|
|
173
|
spec[0].min(0.01)
|
|
174
|
spec[0].max(10.0)
|
|
175
|
spec[0].scale(1.0)
|
|
176
|
spec[1].free()
|
|
177
|
spec[1].min(-2.0)
|
|
178
|
spec[1].max(2.0)
|
|
179
|
spec[1].scale(1.0)
|
|
180
|
else:
|
|
181
|
spec = GModelSpectralNodes()
|
|
182
|
bounds = GEbounds(self.nodes,GEnergy(emin,"TeV"),GEnergy(emax,"TeV"),True)
|
|
183
|
for i in range(bounds.size()):
|
|
184
|
spec.append(bounds.elogmean(i),1.0)
|
|
185
|
for par in spec:
|
|
186
|
if "Energy" in par.name():
|
|
187
|
par.fix()
|
|
188
|
elif "Intensity" in par.name():
|
|
189
|
par.min(0.01)
|
|
190
|
par.max(10.0)
|
|
191
|
bck = GCTAModelIrfBackground(spec)
|
|
192
|
bck.ids(run)
|
|
193
|
bck.instruments("HESS")
|
|
194
|
bck.name("bck_"+run)
|
|
195
|
bck.tscalc(False)
|
|
196
|
return bck
|
|
197
|
|
|
198
|
|
|
199
|
def erange(self,emin,emax):
|
|
200
|
if emin>=emax:
|
|
201
|
self.error("Energy range ["+str(emin)+";"+str(emax)+"] not valid",sys._getframe().f_code.co_name)
|
|
202
|
self.emin = emin
|
|
203
|
self.emax = emax
|
|
204
|
|
|
205
|
def check_run(self,runnr):
|
|
206
|
if not os.path.exists(self._eventsfile(runnr)):
|
|
207
|
return False
|
|
208
|
if not os.path.exists(self._psffile(runnr)):
|
|
209
|
return False
|
|
210
|
if not os.path.exists(self._Aefffile(runnr)):
|
|
211
|
return False
|
|
212
|
return True
|
|
213
|
|
|
214
|
def load(self):
|
|
215
|
self.obs.clear()
|
|
216
|
models = GModels()
|
|
217
|
if self.verbose:
|
|
218
|
self.info("Loading "+str(len(self.runlist))+" H.E.S.S. runs...")
|
|
219
|
|
|
220
|
for run in self.runlist:
|
|
221
|
if not self.check_run(run):
|
|
222
|
if self.verbose:
|
|
223
|
self.warning("Run "+run+" does not exist")
|
|
224
|
continue
|
|
225
|
evlist = self._eventsfile(run)
|
|
226
|
[alt,az,ntels] = fits_header(evlist,"ALT_PNT","AZ_PNT","N_TELS")
|
|
227
|
if alt == 0.0 or az == 0.0 or ntels<=2 or ntels>5:
|
|
228
|
if self.verbose:
|
|
229
|
self.warning("Run "+run+" not valid (Alt="+str(alt)+", Az="+str(az)+", Ntels="+str(ntels)+")")
|
|
230
|
continue
|
|
231
|
altbin,azbin = self._findbins(alt, az)
|
|
232
|
bgfile = self.bgpath+"/hist_alt"+str(altbin)+"_az"+str(azbin)+".fits"
|
|
233
|
if not os.path.exists(bgfile):
|
|
234
|
if self.verbose:
|
|
235
|
self.warning("Run "+run+" has no corresponding off run (run will be skipped, Altbin="+str(altbin)+", Azbin="+str(azbin)+")")
|
|
236
|
continue
|
|
237
|
runemax = fits_header(self._Aefffile(run),"HI_THRES")
|
|
238
|
runemin = fits_header(bgfile,"E_THRES")
|
|
239
|
|
|
240
|
minenergy=0
|
|
241
|
maxenergy=0
|
|
242
|
if self.emin < runemin:
|
|
243
|
minenergy=runemin
|
|
244
|
else:
|
|
245
|
minenergy=self.emin
|
|
246
|
|
|
247
|
if self.emax > runemax:
|
|
248
|
maxenergy = runemax
|
|
249
|
else:
|
|
250
|
maxenergy = self.emax
|
|
251
|
|
|
252
|
if maxenergy <= minenergy:
|
|
253
|
if self.verbose:
|
|
254
|
self.warning("Run "+run+ " has no event after selection")
|
|
255
|
continue
|
|
256
|
|
|
257
|
if minenergy<= self.runlist_emin:
|
|
258
|
self.runlist_emin = minenergy
|
|
259
|
if maxenergy >=self.runlist_emax:
|
|
260
|
self.runlist_emax=maxenergy
|
|
261
|
|
|
262
|
response = GCTAResponseIrf()
|
|
263
|
response.load_psf(self._psffile(run))
|
|
264
|
response.load_aeff(self._Aefffile(run))
|
|
265
|
response.load_background(bgfile)
|
|
266
|
|
|
267
|
sel_hess_obs = self._select(evlist,minenergy,maxenergy)
|
|
268
|
sel_hess_obs.response(response)
|
|
269
|
sel_hess_obs.id(run)
|
|
270
|
|
|
271
|
self.obs.append(sel_hess_obs)
|
|
272
|
bck = self.hess_inst_background(run,minenergy,maxenergy)
|
|
273
|
models.append(bck)
|
|
274
|
|
|
275
|
self.obs.models(models)
|
|
276
|
|
|
277
|
if self.verbose and self.obs.size()>0:
|
|
278
|
self.success("... loaded "+str(self.obs.size())+ " H.E.S.S. runs (Energy range: "+str("%.2f" % self.runlist_emin)+" - "+ str("%.2f" % self.runlist_emax)+" TeV)")
|
|
279
|
elif self.verbose:
|
|
280
|
self.warning("no run has been loaded")
|