openmp.py
1 |
#! /usr/bin/env python
|
---|---|
2 |
# ==========================================================================
|
3 |
# This script simulates a number of CTA observations and fits a model to the
|
4 |
# simulated events. To make this script work, make sure that the $GAMMALIB
|
5 |
# environment variable points to the directory whene gammalib and ctools are
|
6 |
# installed. Here an example of how to configure gammalib and ctools:
|
7 |
#
|
8 |
# Suppose you have two directories, one holding a copy of gammalib and one
|
9 |
# holding a copy of ctools. Then install both components in a specific
|
10 |
# install directory by executing the following commands:
|
11 |
#
|
12 |
# cd gammalib
|
13 |
# ./configure --prefix=/home/knodlseder/testinstall
|
14 |
# make install
|
15 |
# cd ..
|
16 |
# cd ctools
|
17 |
# ./configure --prefix=/home/knodlseder/testinstall
|
18 |
# make install
|
19 |
#
|
20 |
# Now you need to configure gammalib and ctools:
|
21 |
#
|
22 |
# export GAMMALIB=/home/knodlseder/testinstall
|
23 |
# source $GAMMALIB/bin/gammalib-init.sh
|
24 |
#
|
25 |
# Copyright (C) 2012 Juergen Knoedlseder
|
26 |
#
|
27 |
# This program is free software: you can redistribute it and/or modify
|
28 |
# it under the terms of the GNU General Public License as published by
|
29 |
# the Free Software Foundation, either version 3 of the License, or
|
30 |
# (at your option) any later version.
|
31 |
#
|
32 |
# This program is distributed in the hope that it will be useful,
|
33 |
# but WITHOUT ANY WARRANTY; without even the implied warranty of
|
34 |
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
35 |
# GNU General Public License for more details.
|
36 |
#
|
37 |
# You should have received a copy of the GNU General Public License
|
38 |
# along with this program. If not, see <http://www.gnu.org/licenses/>.
|
39 |
#
|
40 |
# ==========================================================================
|
41 |
from ctools import * |
42 |
from gammalib import * |
43 |
import obsutils |
44 |
import sys |
45 |
import time |
46 |
|
47 |
|
48 |
# ========= #
|
49 |
# cta class #
|
50 |
# ========= #
|
51 |
class cta: |
52 |
"""
|
53 |
Implements methods for observation definition, simulations and model
|
54 |
fitting.
|
55 |
"""
|
56 |
def __init__(self, *argv): |
57 |
"""
|
58 |
Constructor.
|
59 |
"""
|
60 |
# Initialise members
|
61 |
self.m_obs = GObservations()
|
62 |
self.m_srcmdl = "crab.xml" |
63 |
self.m_caldb = "${GAMMALIB}/share/caldb/cta" |
64 |
self.m_irf = "kb_E_50h_v3" |
65 |
self.m_ra = 83.6331 |
66 |
self.m_dec = 22.0145 |
67 |
self.m_emin = 0.1 |
68 |
self.m_emax = 100.0 |
69 |
self.m_duration = 1800.0 |
70 |
self.m_deadc = 0.95 |
71 |
self.m_rad = 5.0 |
72 |
|
73 |
# Return
|
74 |
return
|
75 |
|
76 |
def __del__(self): |
77 |
"""
|
78 |
Destructor.
|
79 |
"""
|
80 |
# Return
|
81 |
return
|
82 |
|
83 |
def load(self, filename): |
84 |
"""
|
85 |
"""
|
86 |
# Load data
|
87 |
self.m_obs.load(filename)
|
88 |
|
89 |
# Load model
|
90 |
self.m_obs.models(GModels(self.m_srcmdl)) |
91 |
print self.m_obs |
92 |
|
93 |
# Return
|
94 |
return
|
95 |
|
96 |
def simulate(self, number): |
97 |
"""
|
98 |
Simulation CTA observations.
|
99 |
|
100 |
Parameters:
|
101 |
number - Number of observations in container.
|
102 |
"""
|
103 |
# Set seed value
|
104 |
seed = 0
|
105 |
|
106 |
# Clear observation container
|
107 |
self.m_obs.clear()
|
108 |
|
109 |
# Set pointing direction
|
110 |
pntdir = GSkyDir() |
111 |
pntdir.radec_deg(self.m_ra, self.m_dec) |
112 |
|
113 |
# Create CTA observation
|
114 |
#run = obsutils.set(pntdir, caldb=self.m_caldb, irf=self.m_irf, \
|
115 |
# duration=self.m_duration, deadc=self.m_deadc, \
|
116 |
# emin=self.m_emin, emax=self.m_emax, rad=self.m_rad)
|
117 |
|
118 |
# Append observation to container
|
119 |
#for i in range(number):
|
120 |
# self.m_obs.append(run)
|
121 |
|
122 |
# Append source model
|
123 |
self.m_obs.models(GModels(self.m_srcmdl)) |
124 |
|
125 |
# Simulate observations
|
126 |
#self.m_obs = obsutils.sim(self.m_obs, seed=seed, debug=True)
|
127 |
|
128 |
# Return
|
129 |
return
|
130 |
|
131 |
def fit(self): |
132 |
"""
|
133 |
Fit model.
|
134 |
"""
|
135 |
# Save start time
|
136 |
start = time.clock() |
137 |
|
138 |
# Fit model
|
139 |
like = obsutils.fit(self.m_obs, debug=True) |
140 |
|
141 |
# Compute elapsed time
|
142 |
elapsed = (time.clock() - start) |
143 |
print "Elapsed time: ", elapsed |
144 |
|
145 |
# Return
|
146 |
return
|
147 |
|
148 |
|
149 |
# ======================== #
|
150 |
# Main routine entry point #
|
151 |
# ======================== #
|
152 |
if __name__ == '__main__': |
153 |
"""
|
154 |
This script simulates a number of CTA observations and fits a model
|
155 |
to the simulated events. The command line parameter specifies the
|
156 |
number of CTA observations that should be simulated. By default, 100
|
157 |
observations are simulated.
|
158 |
"""
|
159 |
# Get input arguments
|
160 |
usage = "openmp [number of observations]"
|
161 |
if len(sys.argv) < 1 or len(sys.argv) > 2: |
162 |
print usage
|
163 |
sys.exit() |
164 |
|
165 |
# Set number of observations
|
166 |
if len(sys.argv) > 1: |
167 |
number = int(sys.argv[1]) |
168 |
else:
|
169 |
number = 100
|
170 |
|
171 |
# Create instance of test class
|
172 |
obs = cta() |
173 |
|
174 |
# Simulate observations
|
175 |
print "1" |
176 |
#obs.simulate(number)
|
177 |
|
178 |
# Load data
|
179 |
print "2" |
180 |
obs.load("events.xml")
|
181 |
|
182 |
# Do model fitting
|
183 |
print "3" |
184 |
obs.fit() |
185 |
|