check_models.py
3.62 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
136
137
#! /usr/bin/env python
# ==========================================================================
# This script performs an unbinned maximum likelihood analysis for a
# variety of models. This allows checking the models.
#
# Usage:
# ./check_models.py
#
# ==========================================================================
#
# Copyright (C) 2011-2015 Juergen Knoedlseder
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# ==========================================================================
import gammalib
import ctools
import sys
# ================= #
# Analysis pipeline #
# ================= #
def pipeline(model_name):
"""
Unbinned analysis pipeline - keep intermediate results in memory.
This function implements an analysis pipeline that successively calls
ctobssim, ctselect and ctlike without saving the intermediate results as
FITS files on disk. All data is only hold in memory.
"""
# Set script parameters
caldb = "prod2"
irf = "South_50h"
ra = 83.63
dec = 22.01
rad_sim = 10.0
tstart = 0.0
tstop = 1800.0
emin = 0.1
emax = 100.0
rad_select = 3.0
# Write model name
print("*** Model: "+model_name+" ************************************")
# Initialise timing
wall_seconds = 0.0
cpu_seconds = 0.0
# Simulate events
sim = ctools.ctobssim()
sim["inmodel"].filename(model_name)
sim["caldb"].string(caldb)
sim["irf"].string(irf)
sim["ra"].real(ra)
sim["dec"].real(dec)
sim["rad"].real(rad_sim)
sim["tmin"].real(tstart)
sim["tmax"].real(tstop)
sim["emin"].real(emin)
sim["emax"].real(emax)
sim.run()
# Select events
select = ctools.ctselect(sim.obs())
select["ra"].real(ra)
select["dec"].real(dec)
select["rad"].real(rad_select)
select["tmin"].real(tstart)
select["tmax"].real(tstop)
select["emin"].real(emin)
select["emax"].real(emax)
select.run()
# Perform maximum likelihood fitting
like = ctools.ctlike(select.obs())
like.run()
# Show model fitting results
print(like.obs().models()[0])
# Return
return
#==========================#
# Main routine entry point #
#==========================#
if __name__ == '__main__':
"""
Perform unbinned analyses for various models.
"""
# Initialise flags
need_help = False
# Test for command line arguments
print(sys.argv[0])
if (len(sys.argv) > 1):
if sys.argv[1] == "-h":
need_help = True
else:
need_help = True
# Print help if needed and exit
if need_help:
print("Usage: check_models.py [OPTIONS]")
print(" -h Display this usage message")
sys.exit()
# Dump header
print("*******************************************")
print("* Check models using an unbinned analysis *")
print("*******************************************")
# Perform analysis for Crab model
pipeline("$CTOOLS/share/models/crab.xml")
# Perform analysis for disk model
pipeline("$CTOOLS/share/models/disk.xml")
# Perform analysis for Gaussian model
pipeline("$CTOOLS/share/models/gauss.xml")
# Perform analysis for shell model
pipeline("$CTOOLS/share/models/shell.xml")