-
Notifications
You must be signed in to change notification settings - Fork 2
/
VOTResponse.py
405 lines (294 loc) · 17.4 KB
/
VOTResponse.py
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
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
#!/usr/bin/env python
import numpy as np
from VOTUtils import initLogger, checkParams
'''This reads in a VHE IRF and calculates the response to a inputted
spectrum. This can be from the extrapolatedSpectrum class or you can
input a file with energy and dNdE as colums. Energy is in GeV and
dNdE is in photons*cm^-2*GeV^-1.'''
class VOTResponse:
def __init__(self,instrument="VERITAS", **pars):
self.logger = initLogger('VOTSpectrum')
self.loadEA(instrument, **pars)
def loadEA(self, instrument = "VERITAS", **pars):
''' - Instruments -
These are the currently avaialble instruments.
'VERITAS': this is an array of IACTs located in Southern
Arizona and sensitive from about 100 GeV to tens of TeV. You
can choose a variety of zenith angles (0 is pointing straight
up) and noise levels (ie. night sky background), and expected
source spectrum. The default noise level is adequate for most
fields and the default expected source spectrum (-4) is
adequate for most sources at VHE.
'VERITAS2': same as VERITAS but with multivariate
interpolation in energy and zenith angle. You can pick any
zenith angle from 20 to 90.
'HESS': This is an array of IACTs located in Nambia and
sensitive from about 100 GeV to tens of TeV. You can choose
three zenith angles (20, 45 and 60). The effective areas are
from A&A 457, 899 (2006), Figure 13. The threshold is set to
150 GeV.
'HAWC': This is a water Cherenkov experiment located in
Mexico. The effective areas are from APh 35, 6412 (2012),
Figure 2. You can choose from a high trigger rate (nHit > 70)
or a low trigger rate (nHit > 30). You can also choose from
four zenith angle rages.
'CTA': this is a future IACT observatory that will be an order
of magnitude more sensitive then the current generation. The
Effective Area curve is from "Monte Carlo design studies for
the Cherenkov Telescope Array", Bernlohr et
al. arXiv:1210.3503, Figure 15 and assume the baseline/MPIK
curve. Note that this is for demonstration purposes only and
there are no modifiers for zenith angle, noise level or
anything else. The safe energy range is hardcoded from 100
GeV to 100 TeV and the threshold is set to 100 GeV.
'''
EAParamTemplates = {"VERITAS" : {"zenith":20, "azimuth":0, "noise":4.07},
"VERITAS2" : {"zenith":20, "azimuth":0, "noise":4.07},
"CTA" : {"zenith":0, "azimuth":0, "noise":0},
"HESS" : {"zenith":20, "azimuth":0, "noise":0},
"HESS2" : {"zenith":0, "azimuth":0, "noise":0},
"HAWC" : {"zenith":20}}
EffectiveAreas = {"VERITAS": self.loadVERITASEA,
"VERITAS2": self.loadVERITASEA2,
"CTA": self.loadCTAEA,
"HESS": self.loadHESSEA,
"HESS2": self.loadHESS2EA,
"HAWC": self.loadHAWCEA}
try:
checkParams(EAParamTemplates[instrument],pars, self.logger)
except KeyError:
self.logger.critical("Unsupported instrument template.")
return
try:
self.EASummary,self.EACurve,self.SensCurve,self.EAFile,self.EATable,self.crabRate = EffectiveAreas[instrument](**pars)
except KeyError:
self.logger.critical("Unsupported instrument effective area.")
return
def interpolateEA(self, E, EA):
'''Calculates EA for specific energies [GeV]'''
return np.interp(E, 10**EA[0:,0], EA[0:,1])
def convolveSpectrum(self, EnergyBins, dNdE, EA, minE, maxE):
'''Integrates a spectrum over an effective area. dNdE and EA must
have the same energy axis (ie. the same x axis points). minE
and maxE are assumed to be inside the limits of E. '''
if(EnergyBins[0] > minE or EnergyBins[-1] <maxE):
self.logger.warn("Not integrating over full VHE response")
counts_per_second = 0
it = np.nditer(EnergyBins, flags=['f_index'])
while not it.finished:
if(it[0] >= minE and it[0] <= maxE):
counts_per_second = counts_per_second + (EnergyBins[it.index] - EnergyBins[it.index-1])*dNdE[it.index]*EA[it.index]
it.iternext()
return counts_per_second*60.
def loadSensitivity(self, instrument='VERITAS'):
'''This loads a sensitivity curve for a specific instrument. Note
that I don't have a sensitivity curve for CTA yet.'''
if instrument == 'CTA':
self.logger.warning("Assuming 10 times better sensitivity than HESS for CTA!")
filename = "Effective_Areas/HESS/sensitivity.csv"
sensitivity = np.genfromtxt(filename,delimiter=",")
sensitivity[:,1] = 0.1*sensitivity[:,1]
return sensitivity
elif instrument == 'HAWC':
self.logger.warning("Assuming 1000 times worse sensitivity than HESS for HAWC!")
filename = "Effective_Areas/HESS/sensitivity.csv"
sensitivity = np.genfromtxt(filename,delimiter=",")
sensitivity[:,1] = 1000.*sensitivity[:,1]
return sensitivity
elif instrument == 'HESS2':
self.logger.warning("Assuming HESS sensitivity but using HESS2 effective areas!")
filename = "Effective_Areas/HESS/sensitivity.csv"
return np.genfromtxt(filename,delimiter=",")
else:
self.logger.warning("Time to detection assumes a Crab Nebula Spectrum.")
filename = "Effective_Areas/"+instrument+"/sensitivity.csv"
return np.genfromtxt(filename,delimiter=",")
def loadVERITASEA(self, cuts = 'soft', **pars):
'''This loads in the effective areas for VERITAS. All units returned
in GeV and cm^2.'''
name = {'soft': "Effective_Areas/VERITAS/ea_Nov2010_na_ATM21_vegasv240rc1_7sam_050off_soft-1",
'medium': "Effective_Areas/VERITAS/ea_Nov2010_na_ATM21_vegasv240rc1_7sam_050off_med-1",
'hard': "Effective_Areas/VERITAS/ea_Nov2010_na_ATM21_vegasv240rc1_7sam_050off_hard-1"}
crabRate = 659.0
self.logger.warning("Using {} counts/hr as the rate from the Crab Nebula.".format(crabRate))
filename = name[cuts] + ".summary.csv"
summary_data = np.genfromtxt(filename, delimiter=",", unpack=True,
dtype=[('eanames', 'S50'),
('minSafeE', 'float'),
('maxSafeE', 'float'),
('peakArea', 'float'),
('threshold', 'float')])
EAName = "EffectiveArea_Azimuth_" + str(pars["azimuth"]) + "_Zenith_" + str(pars["zenith"]) + "_Noise_" + str(pars["noise"])
array_mask = summary_data['eanames'] == EAName
try:
EASummary = {'eaname': EAName,
'minSafeE': ((summary_data['minSafeE'])[array_mask])[0] + 3,
'maxSafeE': ((summary_data['maxSafeE'])[array_mask])[0] + 3,
'peakArea': ((summary_data['peakArea'])[array_mask])[0] * 10000.,
'threshold':((summary_data['threshold'])[array_mask])[0] * 1000.}
except IndexError:
self.logger.critical('Could not find that EA curve.')
return 0,0
EACurveFileName = name[cuts]+"/"+EAName+".csv"
EACurve_data = np.genfromtxt(EACurveFileName,delimiter=",")
EACurve_data = EACurve_data + [3.,0.]
EACurve_data = EACurve_data * [1.,10000.]
Sensitivity_data = self.loadSensitivity('VERITAS')
return EASummary, EACurve_data, Sensitivity_data, filename, EAName, crabRate
def loadVERITASEA2(self, cuts = 'soft', **pars):
name = {'soft': "Effective_Areas/VERITAS/ea_Nov2010_na_ATM21_vegasv240rc1_7sam_050off_soft-1",
'medium': "Effective_Areas/VERITAS/ea_Nov2010_na_ATM21_vegasv240rc1_7sam_050off_med-1",
'hard': "Effective_Areas/VERITAS/ea_Nov2010_na_ATM21_vegasv240rc1_7sam_050off_hard-1"}
crabRate = 659.0
self.logger.warning("Using {} counts/hr as the rate from the Crab Nebula.".format(crabRate))
filename = name[cuts] + ".summary.csv"
summary_data = np.genfromtxt(filename, delimiter=",", unpack=True,
dtype=[('eanames', 'S50'),
('minSafeE', 'float'),
('maxSafeE', 'float'),
('peakArea', 'float'),
('threshold', 'float')])
zeniths = list({eaname.split('_')[4] for eaname in summary_data['eanames']})
zeniths.sort()
idx = np.searchsorted(np.float64(zeniths),pars['zenith'])
print idx
EANames = ("EffectiveArea_Azimuth_" + str(pars["azimuth"]) + "_Zenith_" + zeniths[idx-1] + "_Noise_" + str(pars["noise"]),
"EffectiveArea_Azimuth_" + str(pars["azimuth"]) + "_Zenith_" + zeniths[idx] + "_Noise_" + str(pars["noise"]))
array_mask = summary_data['eanames'] == EANames[0]
try:
EASummary = {'eaname': EANames[0],
'minSafeE': ((summary_data['minSafeE'])[array_mask])[0] + 3,
'maxSafeE': ((summary_data['maxSafeE'])[array_mask])[0] + 3,
'peakArea': ((summary_data['peakArea'])[array_mask])[0] * 10000.,
'threshold':((summary_data['threshold'])[array_mask])[0] * 1000.}
except IndexError:
self.logger.critical('Could not find that EA curve.')
return 0,0
EACurveFileNames = (name[cuts]+"/"+EANames[0]+".csv",name[cuts]+"/"+EANames[1]+".csv")
EACurves_data = [np.genfromtxt(EACurveFileNames[0],delimiter=","),np.genfromtxt(EACurveFileNames[1],delimiter=",")]
EACurves_data[0] = EACurves_data[0] + [3.,0.]
EACurves_data[0] = EACurves_data[0] * [1.,10000.]
EACurves_data[1] = EACurves_data[1] + [3.,0.]
EACurves_data[1] = EACurves_data[1] * [1.,10000.]
# Sometimes the length of one of the EA's is smaller than the other.
if(np.shape(EACurves_data[0]) != np.shape(EACurves_data[1])):
small_idx = np.argmin([np.shape(EACurve_data)[0] for EACurve_data in EACurves_data])
start_idx = np.argmax(EACurves_data[small_idx-1][:,0] == EACurves_data[small_idx][0,0])
if start_idx > 0:
EACurves_data[small_idx] = np.insert(EACurves_data[small_idx],0,
zip(EACurves_data[small_idx-1][:start_idx,0],np.zeros(start_idx)),axis=0)
else:
end_idx = np.shape(EACurves_data[small_idx-1][:,0])[0] - np.shape(EACurves_data[small_idx][:,0])[0]
EACurves_data[small_idx] = np.insert(EACurves_data[small_idx],np.shape(EACurves_data[small_idx])[0],
zip(EACurves_data[small_idx-1][-end_idx:,0],np.zeros(end_idx)),axis=0)
interpEA = np.array([[EACurve_data[0],np.interp(pars['zenith'],[zeniths[idx-1],zeniths[idx]],[EACurve_data[1],EACurves_data[1][i,1]])]
for i,EACurve_data in enumerate(EACurves_data[0])])
print "hello"
Sensitivity_data = self.loadSensitivity('VERITAS')
return EASummary, interpEA, Sensitivity_data, filename, EANames[0], crabRate
def loadHESSEA(self, zenith = 20, **pars):
'''This loads in the effective areas for HESS. All units returned in
GeV and cm^2. There are three zenith angles (20, 45 and 60)
which each have their own safe energy range.'''
crabRate = 1040.0
self.logger.warning("Using {} counts/hr as the rate from the Crab Nebula.".format(crabRate))
minEnergies = {"20": 0.174399764328,
"45": 0.395019338585,
"60": 1.14280205981}
maxEnergies = {"20": 25.0,
"45": 45.0,
"60": 70.0}
try:
minE = minEnergies[str(zenith)]
maxE = maxEnergies[str(zenith)]
except IndexError:
self.logger.critical('Could not find that EA curve.')
return 0,0
EAName = "EA_AA457_Fig13_True_Zenith_"+str(zenith)
EASummary = {'eaname': EAName,
'minSafeE': np.log10(minE*1000.),
'maxSafeE': np.log10(maxE*1000.),
'peakArea': 0,
'threshold': 150.}
EACurveFileName = "Effective_Areas/HESS/" + EAName + ".csv"
EACurve_data = np.genfromtxt(EACurveFileName,delimiter=",")
EASummary['peakArea'] = np.max(EACurve_data[:,0]*10000.)
EACurve_data[:,0] = np.log10(EACurve_data[:,0]*1000.)
EACurve_data = EACurve_data * [1.,10000.]
Sensitivity_data = self.loadSensitivity('HESS')
return EASummary, EACurve_data, Sensitivity_data, EACurveFileName, EAName, crabRate
def loadHESS2EA(self, analysisChain=2,**pars):
'''This loads in the effective areas for HESS2 (mono). All units returned in
GeV and cm^2. There is only one (unknown) zenith.'''
crabRate = 2150.0
self.logger.warning('This EA is for demonstration purposes only. All paramters (zenith etc.) are ignored.')
self.logger.warning("Using {} counts/hr as the rate from the Crab Nebula.".format(crabRate))
if(analysisChain == 2):
EAName = 'EA_1307.6003v1_Fig2_AC2'
else:
EAName = 'EA_1307.6003v1_Fig2_AC1'
EACurveFileName="Effective_Areas/HESS2/"+EAName+".csv"
EASummary = {'eaname': EAName,
'minSafeE': np.log10(50.),
'maxSafeE': np.log10(990.),
'peakArea': 92000*10000.,
'threshold': 50.}
EACurve_data = np.genfromtxt(EACurveFileName,delimiter=",")
EACurve_data[:,0] = EACurve_data[:,0] + 3.0
EACurve_data = EACurve_data * [1.,10000.]
Sensitivity_data = self.loadSensitivity('HESS2')
return EASummary, EACurve_data, Sensitivity_data, EACurveFileName, "EA_1307.6003v1_Fig2_AC2", crabRate
def loadCTAEA(self, EACurveFileName="Effective_Areas/CTA/EA_1210.3503_Fig15_MPIK.csv", **pars):
'''This loads in the effective areas for CTA. All units returned in
GeV and cm^2. Note that this is for demonstration purposes
only and there are no modifiers for zenith angle, noise level
or anything else. The safe energy range is hardcoded from 100
GeV to 100 TeV and the threshold is set to 100 GeV.'''
crabRate = 1760.0
self.logger.warning('This EA is for demonstration purposes only. All paramters (zenith etc.) are ignored.')
self.logger.warning("Using {} counts/hr as the rate from the Crab Nebula.".format(crabRate))
EASummary = {'eaname': 'EA_1210.3503_Fig15_MPIK',
'minSafeE': 2.0,
'maxSafeE': 5.0,
'peakArea': 3.37550925e10,
'threshold': 100.}
EACurve_data = np.genfromtxt(EACurveFileName,delimiter=",")
EACurve_data[:,0] = np.log10(EACurve_data[:,0]*1000.)
EACurve_data = EACurve_data * [1.,10000.]
Sensitivity_data = self.loadSensitivity('CTA')
return EASummary, EACurve_data, Sensitivity_data, EACurveFileName, "EA_1210.3503_Fig15_MPIK", crabRate
def loadHAWCEA(self, zenith = 20, trigger_rate = "high", **pars):
'''This loads the effective area for HAWC. All units returned in GeV
and cm^2.'''
crabRate = 19.
theta_bins = ['0607','0708','0809','0910']
theta_bounds = [0.6,0.7,0.8,0.9,1.0]
if trigger_rate == 'high':
nHit = '70'
minEnergies = {"0607": 200.,
"0708": 100.,
"0809": 63.1,
"0910": 31.6}
else:
nHit = '30'
minEnergies = {"0607": 100.,
"0708": 50.,
"0809": 25.,
"0910": 15.}
cos_zenith = np.cos(zenith*np.pi/180)
if cos_zenith <= 0.6:
theta_bin = theta_bins[0]
else:
theta_bin = theta_bins[np.searchsorted(theta_bounds,cos_zenith)-2]
EAName = "EA_APh35_Fig2_nHit_"+nHit+"_cosTh_"+theta_bin
EACurveFileName = "Effective_Areas/HAWC/" + EAName + ".csv"
EASummary = {'eaname': EAName,
'minSafeE': np.log10(minEnergies[theta_bin]),
'maxSafeE': np.log10(99000.),
'peakArea': 1.0,
'threshold': 100.}
EACurve_data = np.genfromtxt(EACurveFileName,delimiter=",")
#EACurve_data[:,0] = np.log10(EACurve_data[:,0])+3.0
EACurve_data = EACurve_data * [1.,10000.]
Sensitivity_data = self.loadSensitivity('HAWC')
return EASummary, EACurve_data, Sensitivity_data, EACurveFileName, EAName, crabRate