1
|
import gammalib
|
2
|
import ctools
|
3
|
import cscripts
|
4
|
|
5
|
ra = 350.8524
|
6
|
dec = 58.8159
|
7
|
max_offset = 3.
|
8
|
srcrad = 0.2
|
9
|
stacking = True
|
10
|
emin = 0.1
|
11
|
emax = 50.
|
12
|
enumbins = 30
|
13
|
|
14
|
obs = 'obs.xml'
|
15
|
obsselect = cscripts.csobsselect()
|
16
|
obsselect['inobs'] = '$CTADATA/obs/obs_gps_baseline.xml'
|
17
|
obsselect['outobs'] = obs
|
18
|
obsselect['pntselect'] = 'CIRCLE'
|
19
|
obsselect['coordsys'] = 'CEL'
|
20
|
obsselect['ra'] = ra
|
21
|
obsselect['dec'] = dec
|
22
|
obsselect['rad'] = max_offset
|
23
|
obsselect['tmin'] = 'UNDEF'
|
24
|
obsselect.execute()
|
25
|
|
26
|
obs_selected = "obs_selected.xml"
|
27
|
select = ctools.ctselect()
|
28
|
select['inobs'] = obs
|
29
|
select['ra'] = 'UNDEF'
|
30
|
select['tmin'] = 'UNDEF'
|
31
|
select['emin'] = emin
|
32
|
select['emax'] = emax
|
33
|
select['expr'] = 'MC_ID==78'
|
34
|
select['outobs'] = obs_selected
|
35
|
select.execute()
|
36
|
|
37
|
|
38
|
model_sky = 'model_sky.xml'
|
39
|
modelselect = cscripts.csmodelselect()
|
40
|
modelselect['inobs'] = obs_selected
|
41
|
modelselect['inmodel'] = '$CTADATA/models/model_galactic_bright.xml'
|
42
|
modelselect['outmodel'] = model_sky
|
43
|
modelselect.execute()
|
44
|
|
45
|
|
46
|
modelfile = model_sky
|
47
|
model = gammalib.GModels(modelfile)
|
48
|
for source in model:
|
49
|
if source.name() == 'Cassiopeia A':
|
50
|
for par in source:
|
51
|
if par.name() == 'RA' or par.name() == 'DEC':
|
52
|
par.fix()
|
53
|
model['Cassiopeia A'].tscalc(True)
|
54
|
model.save(modelfile)
|
55
|
|
56
|
|
57
|
onoff_obs = 'onoff_obs.xml'
|
58
|
phatool = cscripts.csphagen()
|
59
|
phatool['inobs'] = obs_selected
|
60
|
phatool['outobs'] = onoff_obs
|
61
|
phatool['ebinalg'] = 'LOG'
|
62
|
phatool['emin'] = emin
|
63
|
phatool['emax'] = emax
|
64
|
phatool['enumbins'] = enumbins
|
65
|
phatool['coordsys'] = 'CEL'
|
66
|
phatool['ra'] = ra
|
67
|
phatool['dec'] = dec
|
68
|
phatool['rad'] = srcrad
|
69
|
phatool['bkgmethod'] = 'REFLECTED'
|
70
|
phatool['maxoffset'] = max_offset
|
71
|
phatool['stack'] = stacking
|
72
|
phatool['debug'] = True
|
73
|
phatool.execute()
|
74
|
|
75
|
|
76
|
outmodelfile_onoff = 'bestfit_model_onoff.xml'
|
77
|
like = ctools.ctlike()
|
78
|
like['inobs'] = onoff_obs
|
79
|
like['inmodel'] = modelfile
|
80
|
like['statistic'] = 'WSTAT'
|
81
|
like['outmodel'] = outmodelfile_onoff
|
82
|
like.execute()
|
83
|
|
84
|
|
85
|
sed_onoff = "spectrum_issue.fits"
|
86
|
sedmkr = cscripts.csspec()
|
87
|
sedmkr['inobs'] = onoff_obs
|
88
|
sedmkr['inmodel'] = outmodelfile_onoff
|
89
|
sedmkr['srcname'] = 'Cassiopeia A'
|
90
|
sedmkr['outfile'] = sed_onoff
|
91
|
sedmkr['statistic'] = 'WSTAT'
|
92
|
sedmkr['method'] = "AUTO"
|
93
|
sedmkr['emin'] = emin
|
94
|
sedmkr['emax'] = emax
|
95
|
sedmkr['enumbins'] = enumbins
|
96
|
sedmkr['ebinalg'] = 'LOG'
|
97
|
sedmkr['debug'] = True
|
98
|
sedmkr.execute()
|
99
|
|
100
|
|
101
|
outmodelfile_onoff_bkg = 'bestfit_model_onoff_bkg.xml'
|
102
|
bkg = gammalib.GCTAModelIrfBackground( gammalib.GModelSpectralPlaw( 0, 0, gammalib.GEnergy( 0, 'TeV' ) ) )
|
103
|
bkg.instruments( 'CTAOnOff' )
|
104
|
|
105
|
model_bkg = gammalib.GModels(outmodelfile_onoff)
|
106
|
model_bkg.append(bkg)
|
107
|
model_bkg.save(outmodelfile_onoff_bkg)
|
108
|
|
109
|
sed_onoff = "spectrum_works.fits"
|
110
|
sedmkr = cscripts.csspec()
|
111
|
sedmkr['inobs'] = onoff_obs
|
112
|
sedmkr['inmodel'] = outmodelfile_onoff_bkg
|
113
|
sedmkr['srcname'] = 'Cassiopeia A'
|
114
|
sedmkr['outfile'] = sed_onoff
|
115
|
sedmkr['statistic'] = 'WSTAT'
|
116
|
sedmkr['method'] = "AUTO"
|
117
|
sedmkr['emin'] = emin
|
118
|
sedmkr['emax'] = emax
|
119
|
sedmkr['enumbins'] = enumbins
|
120
|
sedmkr['ebinalg'] = 'LOG'
|
121
|
sedmkr['debug'] = True
|
122
|
sedmkr.execute()
|