-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathrun_as.py
108 lines (91 loc) · 3.18 KB
/
run_as.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
import numpy as np
import os
import astropy
import regions
import sherpa
import matplotlib.pyplot as plt
import astropy.units as u
from astropy.coordinates import SkyCoord, Angle
from astropy.table import vstack as vstack_table
from regions import CircleSkyRegion
from gammapy.data import DataStore, ObservationList
from gammapy.data import ObservationStats, ObservationSummary
from gammapy.background.reflected import ReflectedRegionsBackgroundEstimator
from gammapy.utils.energy import EnergyBounds
from gammapy.spectrum import SpectrumExtraction, SpectrumObservation, SpectrumFit, SpectrumResult
#from gammapy.spectrum import SimulationRealBkg
from gammapy.spectrum.models import PowerLaw, ExponentialCutoffPowerLaw, LogParabola
from gammapy.spectrum import FluxPoints, SpectrumEnergyGroupMaker, FluxPointEstimator
from gammapy.image import SkyImage
from gammapy.spectrum import models
from simulation_realbkg import *
#from gammapy.irf import EnergyDispersion, EnergyDispersion2D, EffectiveAreaTable, EffectiveAreaTable2D
# Setup the logger
import logging
logging.basicConfig()
log = logging.getLogger('gammapy.spectrum')
log.setLevel(logging.WARNING)
#choose the obsevartion
name="Crab"
datastore = DataStore.from_dir("$HESS_DATA")
src=SkyCoord.from_name(name)
sep=SkyCoord.separation(src,datastore.obs_table.pointing_radec)
srcruns=(datastore.obs_table[sep<2.0*u.deg]) #2 deg because of HESS field of view.. choose accordingly
obsid=srcruns['OBS_ID'].data
# Define obs parameters
lo_threshold = 0.1 * u.TeV
hi_threshold = 60 * u.TeV
# Define spectral model
index = 3.0 * u.Unit('')
amplitude = 2.5 * 1e-11 * u.Unit('cm-2 s-1 TeV-1')
reference = 1 * u.TeV
model = PowerLaw(index=index, amplitude=amplitude, reference=reference)
# for one obs only
mylist=datastore.obs_list((obsid[0],))
n_obs=1
seeds = np.arange(n_obs)
sim = SimulationRealBkg(source_model=model, obsrun=mylist[0], obspos=src)
sim.run(seeds)
obs=sim.result[0]
fit = SpectrumFit(obs, model=model, stat='wstat')
fit.model.parameters['index'].value = 2
fit.run()
fit.est_errors()
print fit.result[0]
#now run in a loop
n_obs=200
mylist=datastore.obs_list(obsid[0:n_obs])
seeds = np.arange(n_obs)
sims=[]
for i in range(n_obs):
sim = SimulationRealBkg(source_model=model, obsrun=mylist[i], obspos=src)
sim.run(seeds)
sims.append(sim.result[0])
n_on = [obs.total_stats.n_on for obs in sims]
n_off = [obs.total_stats.n_off for obs in sims]
excess = [obs.total_stats.excess for obs in sims]
"""
fix, axes = plt.subplots(1,3, figsize=(12, 4))
axes[0].hist(n_on)
axes[0].set_xlabel('n_on')
axes[1].hist(n_off)
axes[1].set_xlabel('n_off')
axes[2].hist(excess)
axes[2].set_xlabel('excess')
"""
best_fit_index = []
best_fit_flux = []
for obs in sims:
fit = SpectrumFit(obs, model=model, stat='wstat')
fit.model.parameters['index'].value = 2.0
fit.run()
best_fit_index.append(fit.result[0].model.parameters['index'].value)
best_fit_flux.append(fit.result[0].model.parameters['amplitude'].value)
fix, axes = plt.subplots(1,3)
axes[0].hist(best_fit_index)
axes[0].set_xlabel('index')
axes[1].hist(best_fit_flux)
axes[1].set_xlabel('amplitude')
axes[2].plot(best_fit_index,best_fit_flux,"ro")
#print('best_fit_index:', best_fit_index)
plt.show()