Veronika Vodeb / ctools

Commit 320f0d7db00c2935f5077ac5e71ddb644c08472e

Authored by Jürgen Knödlseder
1 parent 2e5b7227

Adapt to new GCaldb interface.

Inline Side-by-side
Showing 2 changed files with 298 additions and 298 deletions
scripts/obsutils.py
  View file @320f0d7
@@ -266,7 +266,7 @@ def modmap(obs, eref=0.1, proj="TAN", coord="GAL", xval=0.0, yval=0.0, \ @@ -266,7 +266,7 @@ def modmap(obs, eref=0.1, proj="TAN", coord="GAL", xval=0.0, yval=0.0, \
266 # ======================= # 266 # ======================= #
267 def set(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \ 267 def set(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \
268 emin=0.1, emax=100.0, rad=5.0, \ 268 emin=0.1, emax=100.0, rad=5.0, \
269 - irf="cta_dummy_irf", caldb=""): 269 + irf="cta_dummy_irf", caldb="$CTOOLS/share/caldb/data/cta/dummy"):
270 """ 270 """
271 Returns a single CTA observation. By looping over this function we can 271 Returns a single CTA observation. By looping over this function we can
272 add CTA observations to the observation container. 272 add CTA observations to the observation container.
test/example_survey.py
  View file @320f0d7
@@ -3,7 +3,7 @@ @@ -3,7 +3,7 @@
3 # This script simulates CTA surveys. It can be considered as the master 3 # This script simulates CTA surveys. It can be considered as the master
4 # script for CTA science simulations. Please modify as needed. 4 # script for CTA science simulations. Please modify as needed.
5 # 5 #
6 -# Copyright (C) 2011-2013 Juergen Knoedlseder 6 +# Copyright (C) 2011-2014 Juergen Knoedlseder
7 # 7 #
8 # This program is free software: you can redistribute it and/or modify 8 # This program is free software: you can redistribute it and/or modify
9 # it under the terms of the GNU General Public License as published by 9 # it under the terms of the GNU General Public License as published by
@@ -32,331 +32,331 @@ import obsutils @@ -32,331 +32,331 @@ import obsutils
32 # Plot counts # 32 # Plot counts #
33 # =========== # 33 # =========== #
34 def plot_counts(obs): 34 def plot_counts(obs):
35 - """  
36 - Plot counts.  
37 - """  
38 - # Only proceed if matplotlib is available  
39 - try:  
40 - # Import matplotlib  
41 - import matplotlib.pyplot as plt  
42 -  
43 - # Set legend fontsize  
44 - params = {'legend.fontsize': 10}  
45 - plt.rcParams.update(params)  
46 -  
47 - # Set plot styles  
48 - styles = ['b-', 'g-', 'y-', 'n-']  
49 -  
50 - # Dump header  
51 - print("")  
52 - print("Make plots (using matplotlib):")  
53 - print("==============================")  
54 -  
55 - # Create figure 1  
56 - plt.figure(1,figsize=(12,6))  
57 - plt.subplots_adjust(hspace=.7)  
58 -  
59 - # Create subplot 1  
60 - plt.subplot(121)  
61 - plt.title("Spectrum (summed over all pixels)")  
62 -  
63 - # Loop over observations  
64 - for run in obs:  
65 -  
66 - # Get event list  
67 - list = run.events()  
68 -  
69 - # Create energy axis  
70 - ebounds = GEbounds()  
71 - emin = GEnergy()  
72 - emax = GEnergy()  
73 - emin.TeV(0.1)  
74 - emax.TeV(100.0)  
75 - ebounds.setlog(emin, emax, 10)  
76 - energy = [ebounds.elogmean(i).TeV() for i in range(ebounds.size())]  
77 -  
78 - # Create spectrum  
79 - print("Extract data:")  
80 - counts = [0.0 for i in range(ebounds.size())]  
81 - print(list.size())  
82 - for atom in list:  
83 - index = ebounds.index(atom.energy())  
84 - counts[index] = counts[index] + 1.0  
85 -  
86 - # Create error bars  
87 - error = [sqrt(c) for c in counts]  
88 -  
89 - # Plot spectrum  
90 - plt.loglog(energy, counts, 'ro', label='data')  
91 - #plt.errorbar(energy, counts, error, fmt=None, ecolor='r')  
92 -  
93 - # Extract models  
94 - print("Extract models:")  
95 - sum_model = [0.0 for i in range(ebounds.size())]  
96 - for k, m in enumerate(obs.models()):  
97 - print("- "+m.name())  
98 - model = [0.0 for i in range(ebounds.size())]  
99 - for atom in list:  
100 - index = ebounds.index(atom.energy())  
101 - prob = m.eval(atom, run)  
102 - model[index] = model[index] + prob * atom.size()  
103 - for i in range(ebounds.size()):  
104 - sum_model[i] = sum_model[i] + model[i]  
105 - #plt.loglog(energy, model, styles[k], label=m.name())  
106 - #plt.loglog(energy, sum_model, 'r-', label='total')  
107 -  
108 - # Put labels  
109 - plt.xlabel("Energy (TeV)")  
110 - plt.ylabel("Counts")  
111 - plt.legend(loc="lower left")  
112 -  
113 - # Create subplot 2  
114 - plt.subplot(122)  
115 - plt.title("Offset (summed over all energies)")  
116 -  
117 - # Set Crab direction  
118 - crab = GSkyDir()  
119 - crab.radec_deg(83.63, 22.01)  
120 -  
121 - # Loop over observations  
122 - for run in obs:  
123 -  
124 - # Get event list  
125 - list = run.events()  
126 -  
127 - # Create offset histogram  
128 - print("Extract data:")  
129 - nx = 30  
130 - doffset2 = 0.01  
131 - offset2 = [(i+0.5)*doffset2 for i in range(nx)]  
132 - counts = [0.0 for i in range(nx)]  
133 - for atom in list:  
134 - off = atom.dir().dist_deg(crab)  
135 - off2 = off*off  
136 - index = int(off2/doffset2)  
137 - if index < nx:  
138 - counts[index] = counts[index] + 1.0  
139 -  
140 - # Create error bars  
141 - error = [sqrt(c) for c in counts]  
142 -  
143 - # Plot distribution  
144 - #plt.semilogy(offset2, counts, 'ro', label='data')  
145 -  
146 - # Extract models  
147 - print("Extract models:")  
148 - sum_model = [0.0 for i in range(nx)]  
149 - for k, m in enumerate(obs.models()):  
150 - print("- "+m.name())  
151 - model = [0.0 for i in range(nx)]  
152 - for atom in list:  
153 - off = atom.dir().dist_deg(crab)  
154 - off2 = off*off  
155 - index = int(off2/doffset2)  
156 - if index < nx:  
157 - prob = m.eval(atom, run)  
158 - model[index] = model[index] + prob * atom.size()  
159 - for i in range(nx):  
160 - sum_model[i] = sum_model[i] + model[i]  
161 - #plt.plot(offset2, model, styles[k], label=m.name())  
162 - #plt.plot(offset2, sum_model, 'r-', label='total')  
163 - #plt.ylim(ymin=0.1)  
164 -  
165 - # Put labels  
166 - plt.xlabel("Offset (deg^2)")  
167 - plt.ylabel("Counts")  
168 - #plt.legend(loc="upper right")  
169 -  
170 - # Show counts spectra  
171 - plt.show()  
172 -  
173 - except ImportError:  
174 - print("Matplotlib is not (correctly) installed on your system.") 35 + """
  36 + Plot counts.
  37 + """
  38 + # Only proceed if matplotlib is available
  39 + try:
  40 + # Import matplotlib
  41 + import matplotlib.pyplot as plt
  42 +
  43 + # Set legend fontsize
  44 + params = {'legend.fontsize': 10}
  45 + plt.rcParams.update(params)
  46 +
  47 + # Set plot styles
  48 + styles = ['b-', 'g-', 'y-', 'n-']
  49 +
  50 + # Dump header
  51 + print("")
  52 + print("Make plots (using matplotlib):")
  53 + print("==============================")
  54 +
  55 + # Create figure 1
  56 + plt.figure(1,figsize=(12,6))
  57 + plt.subplots_adjust(hspace=.7)
  58 +
  59 + # Create subplot 1
  60 + plt.subplot(121)
  61 + plt.title("Spectrum (summed over all pixels)")
  62 +
  63 + # Loop over observations
  64 + for run in obs:
  65 +
  66 + # Get event list
  67 + list = run.events()
  68 +
  69 + # Create energy axis
  70 + ebounds = GEbounds()
  71 + emin = GEnergy()
  72 + emax = GEnergy()
  73 + emin.TeV(0.1)
  74 + emax.TeV(100.0)
  75 + ebounds.setlog(emin, emax, 10)
  76 + energy = [ebounds.elogmean(i).TeV() for i in range(ebounds.size())]
  77 +
  78 + # Create spectrum
  79 + print("Extract data:")
  80 + counts = [0.0 for i in range(ebounds.size())]
  81 + print(list.size())
  82 + for atom in list:
  83 + index = ebounds.index(atom.energy())
  84 + counts[index] = counts[index] + 1.0
  85 +
  86 + # Create error bars
  87 + error = [sqrt(c) for c in counts]
  88 +
  89 + # Plot spectrum
  90 + plt.loglog(energy, counts, 'ro', label='data')
  91 + #plt.errorbar(energy, counts, error, fmt=None, ecolor='r')
  92 +
  93 + # Extract models
  94 + print("Extract models:")
  95 + sum_model = [0.0 for i in range(ebounds.size())]
  96 + for k, m in enumerate(obs.models()):
  97 + print("- "+m.name())
  98 + model = [0.0 for i in range(ebounds.size())]
  99 + for atom in list:
  100 + index = ebounds.index(atom.energy())
  101 + prob = m.eval(atom, run)
  102 + model[index] = model[index] + prob * atom.size()
  103 + for i in range(ebounds.size()):
  104 + sum_model[i] = sum_model[i] + model[i]
  105 + #plt.loglog(energy, model, styles[k], label=m.name())
  106 + #plt.loglog(energy, sum_model, 'r-', label='total')
  107 +
  108 + # Put labels
  109 + plt.xlabel("Energy (TeV)")
  110 + plt.ylabel("Counts")
  111 + plt.legend(loc="lower left")
  112 +
  113 + # Create subplot 2
  114 + plt.subplot(122)
  115 + plt.title("Offset (summed over all energies)")
  116 +
  117 + # Set Crab direction
  118 + crab = GSkyDir()
  119 + crab.radec_deg(83.63, 22.01)
  120 +
  121 + # Loop over observations
  122 + for run in obs:
  123 +
  124 + # Get event list
  125 + list = run.events()
  126 +
  127 + # Create offset histogram
  128 + print("Extract data:")
  129 + nx = 30
  130 + doffset2 = 0.01
  131 + offset2 = [(i+0.5)*doffset2 for i in range(nx)]
  132 + counts = [0.0 for i in range(nx)]
  133 + for atom in list:
  134 + off = atom.dir().dist_deg(crab)
  135 + off2 = off*off
  136 + index = int(off2/doffset2)
  137 + if index < nx:
  138 + counts[index] = counts[index] + 1.0
  139 +
  140 + # Create error bars
  141 + error = [sqrt(c) for c in counts]
  142 +
  143 + # Plot distribution
  144 + #plt.semilogy(offset2, counts, 'ro', label='data')
  145 +
  146 + # Extract models
  147 + print("Extract models:")
  148 + sum_model = [0.0 for i in range(nx)]
  149 + for k, m in enumerate(obs.models()):
  150 + print("- "+m.name())
  151 + model = [0.0 for i in range(nx)]
  152 + for atom in list:
  153 + off = atom.dir().dist_deg(crab)
  154 + off2 = off*off
  155 + index = int(off2/doffset2)
  156 + if index < nx:
  157 + prob = m.eval(atom, run)
  158 + model[index] = model[index] + prob * atom.size()
  159 + for i in range(nx):
  160 + sum_model[i] = sum_model[i] + model[i]
  161 + #plt.plot(offset2, model, styles[k], label=m.name())
  162 + #plt.plot(offset2, sum_model, 'r-', label='total')
  163 + #plt.ylim(ymin=0.1)
  164 +
  165 + # Put labels
  166 + plt.xlabel("Offset (deg^2)")
  167 + plt.ylabel("Counts")
  168 + #plt.legend(loc="upper right")
  169 +
  170 + # Show counts spectra
  171 + plt.show()
  172 +
  173 + except ImportError:
  174 + print("Matplotlib is not (correctly) installed on your system.")
175 175
176 - # Return  
177 - return 176 + # Return
  177 + return
178 178
179 179
180 # ======================================== # 180 # ======================================== #
181 # Add CTA background model to observations # 181 # Add CTA background model to observations #
182 # ======================================== # 182 # ======================================== #
183 def add_background_model(obs): 183 def add_background_model(obs):
184 - """  
185 - Add standard CTA background model to observations container.  
186 -  
187 - We use a simple power law here, scaled to Konrad's E configuration  
188 - performance table. The model needs still to be validated.  
189 - """  
190 - # Recover models from observation  
191 - models = obs.models()  
192 -  
193 - # Define background model  
194 - bgd_radial = GCTAModelRadialGauss(3.0)  
195 - bgd_spectrum = GModelSpectralPlaw(61.8e-6, -1.85, GEnergy(1.0, "TeV"))  
196 - bgd_model = GCTAModelRadialAcceptance(bgd_radial, bgd_spectrum)  
197 - bgd_model.name("Background")  
198 - bgd_model.instruments("CTA")  
199 -  
200 - # Add background model to container  
201 - models.append(bgd_model)  
202 -  
203 - # Put container back in observation container  
204 - obs.models(models)  
205 -  
206 - # Return observation container  
207 - return obs 184 + """
  185 + Add standard CTA background model to observations container.
  186 +
  187 + We use a simple power law here, scaled to Konrad's E configuration
  188 + performance table. The model needs still to be validated.
  189 + """
  190 + # Recover models from observation
  191 + models = obs.models()
  192 +
  193 + # Define background model
  194 + bgd_radial = GCTAModelRadialGauss(3.0)
  195 + bgd_spectrum = GModelSpectralPlaw(61.8e-6, -1.85, GEnergy(1.0, "TeV"))
  196 + bgd_model = GCTAModelRadialAcceptance(bgd_radial, bgd_spectrum)
  197 + bgd_model.name("Background")
  198 + bgd_model.instruments("CTA")
  199 +
  200 + # Add background model to container
  201 + models.append(bgd_model)
  202 +
  203 + # Put container back in observation container
  204 + obs.models(models)
  205 +
  206 + # Return observation container
  207 + return obs
208 208
209 209
210 # ================= # 210 # ================= #
211 # Set Crab spectrum # 211 # Set Crab spectrum #
212 # ================= # 212 # ================= #
213 def crab_spec(): 213 def crab_spec():
214 - """  
215 - Set Crab spectrum based on MAGIC observations  
216 - (Albert et al. 2008, ApJ, 674, 1037)  
217 - """  
218 - # Set parameters  
219 - spectrum = GModelSpectralPlaw(5.7e-16, -2.48, GEnergy(0.3, "TeV"))  
220 -  
221 - # Return spectrum  
222 - return spectrum 214 + """
  215 + Set Crab spectrum based on MAGIC observations
  216 + (Albert et al. 2008, ApJ, 674, 1037)
  217 + """
  218 + # Set parameters
  219 + spectrum = GModelSpectralPlaw(5.7e-16, -2.48, GEnergy(0.3, "TeV"))
  220 +
  221 + # Return spectrum
  222 + return spectrum
223 223
224 224
225 # =============================== # 225 # =============================== #
226 # Setup single observation survey # 226 # Setup single observation survey #
227 # =============================== # 227 # =============================== #
228 def survey_single(): 228 def survey_single():
229 - """  
230 - Creates a single observation survey for test purposes.  
231 - """  
232 - # Allocate observation container  
233 - obs = GObservations()  
234 -  
235 - # Set single pointing at galactic centre  
236 - pntdir = GSkyDir()  
237 - pntdir.lb_deg(0.0, 0.0)  
238 - #run = set_one_obs(pntdir)  
239 - run = obsutils.set(pntdir)  
240 - obs.append(run)  
241 -  
242 - # Define single point source with Crab flux at galactic centre  
243 - center = GSkyDir()  
244 - center.lb_deg(0.0, 0.0)  
245 - point_spatial = GModelSpatialPointSource(center)  
246 - point_spectrum = crab_spec()  
247 - point = GModelSky(point_spatial, point_spectrum)  
248 - point.name('GC source')  
249 -  
250 - # Create model container  
251 - models = GModels()  
252 - models.append(point)  
253 - obs.models(models)  
254 -  
255 - # Return observation container  
256 - return obs 229 + """
  230 + Creates a single observation survey for test purposes.
  231 + """
  232 + # Allocate observation container
  233 + obs = GObservations()
  234 +
  235 + # Set single pointing at galactic centre
  236 + pntdir = GSkyDir()
  237 + pntdir.lb_deg(0.0, 0.0)
  238 + #run = set_one_obs(pntdir)
  239 + run = obsutils.set(pntdir)
  240 + obs.append(run)
  241 +
  242 + # Define single point source with Crab flux at galactic centre
  243 + center = GSkyDir()
  244 + center.lb_deg(0.0, 0.0)
  245 + point_spatial = GModelSpatialPointSource(center)
  246 + point_spectrum = crab_spec()
  247 + point = GModelSky(point_spatial, point_spectrum)
  248 + point.name('GC source')
  249 +
  250 + # Create model container
  251 + models = GModels()
  252 + models.append(point)
  253 + obs.models(models)
  254 +
  255 + # Return observation container
  256 + return obs
257 257
258 258
259 # =========================== # 259 # =========================== #
260 # Setup Galactic plane survey # 260 # Setup Galactic plane survey #
261 # =========================== # 261 # =========================== #
262 def survey_gplane(lrange=10, lstep=2): 262 def survey_gplane(lrange=10, lstep=2):
263 - """  
264 - Creates a single observation survey for test purposes.  
265 -  
266 - Keywords:  
267 - lrange - Longitude range (integer deg)  
268 - lstep - Longitude step size (integer deg)  
269 - """  
270 - # Allocate observation container  
271 - obs = GObservations()  
272 -  
273 - # Loop over longitudes  
274 - for l in range(-lrange,lrange+lstep,lstep):  
275 -  
276 - # Set pointing  
277 - pntdir = GSkyDir()  
278 - pntdir.lb_deg(l, 0.0)  
279 - #run = set_one_obs(pntdir)  
280 - run = obsutils.set(pntdir)  
281 - run.id(str(l))  
282 - obs.append(run)  
283 -  
284 - # Define single point source with Crab flux at galactic centre  
285 - center = GSkyDir()  
286 - center.lb_deg(0.0, 0.0)  
287 - point_spatial = GModelSpatialPointSource(center)  
288 - point_spectrum = crab_spec()  
289 - point = GModelSky(point_spatial, point_spectrum)  
290 - point.name('GC source')  
291 -  
292 - # Create model container  
293 - models = GModels()  
294 - models.append(point)  
295 - obs.models(models)  
296 -  
297 - # Return observation container  
298 - return obs 263 + """
  264 + Creates a single observation survey for test purposes.
  265 +
  266 + Keywords:
  267 + lrange - Longitude range (integer deg)
  268 + lstep - Longitude step size (integer deg)
  269 + """
  270 + # Allocate observation container
  271 + obs = GObservations()
  272 +
  273 + # Loop over longitudes
  274 + for l in range(-lrange,lrange+lstep,lstep):
  275 +
  276 + # Set pointing
  277 + pntdir = GSkyDir()
  278 + pntdir.lb_deg(l, 0.0)
  279 + #run = set_one_obs(pntdir)
  280 + run = obsutils.set(pntdir)
  281 + run.id(str(l))
  282 + obs.append(run)
  283 +
  284 + # Define single point source with Crab flux at galactic centre
  285 + center = GSkyDir()
  286 + center.lb_deg(0.0, 0.0)
  287 + point_spatial = GModelSpatialPointSource(center)
  288 + point_spectrum = crab_spec()
  289 + point = GModelSky(point_spatial, point_spectrum)
  290 + point.name('GC source')
  291 +
  292 + # Create model container
  293 + models = GModels()
  294 + models.append(point)
  295 + obs.models(models)
  296 +
  297 + # Return observation container
  298 + return obs
299 299
300 300
301 #==========================# 301 #==========================#
302 # Main routine entry point # 302 # Main routine entry point #
303 #==========================# 303 #==========================#
304 if __name__ == '__main__': 304 if __name__ == '__main__':
305 - """  
306 - CTA survey simulation and analysis script.  
307 - """  
308 - # Initialise flags  
309 - need_help = False  
310 -  
311 - # Test for command line arguments  
312 - print(sys.argv[0])  
313 - if (len(sys.argv) > 1):  
314 - if sys.argv[1] == "-h":  
315 - need_help = True  
316 - else:  
317 - need_help = True 305 + """
  306 + CTA survey simulation and analysis script.
  307 + """
  308 + # Initialise flags
  309 + need_help = False
  310 +
  311 + # Test for command line arguments
  312 + print(sys.argv[0])
  313 + if (len(sys.argv) > 1):
  314 + if sys.argv[1] == "-h":
  315 + need_help = True
  316 + else:
  317 + need_help = True
318 318
319 - # Print help if needed and exit  
320 - if need_help:  
321 - print("Usage: example_survey.py [OPTIONS]")  
322 - print(" -h Display this usage message")  
323 - sys.exit()  
324 -  
325 - # Dump header  
326 - print("***********************")  
327 - print("* Simulate CTA survey *")  
328 - print("***********************")  
329 -  
330 - # Remove any existing result files  
331 - list = [glob.glob("*.fits"), glob.glob("*.log"), glob.glob("*.xml")]  
332 - for files in list:  
333 - for file in files:  
334 - os.remove(file)  
335 -  
336 - # Setup single observation survey  
337 - #obs = survey_single()  
338 - obs = survey_gplane()  
339 -  
340 - # Add background model  
341 - obs = add_background_model(obs)  
342 -  
343 - # Simulate events  
344 - print("Simulate events")  
345 - obs = obsutils.sim(obs)  
346 -  
347 - # Make counts map  
348 - print("Make counts map")  
349 - cntmap = obsutils.cntmap(obs)  
350 -  
351 - # Fit observations  
352 - print("Fit observations")  
353 - like = obsutils.fit(obs)  
354 - #print like.opt()  
355 - #print like.obs().models()  
356 -  
357 - # Make model map  
358 - print("Make model map (this step will take some time)")  
359 - modmap = obsutils.modmap(obs)  
360 -  
361 - # Show fit results  
362 - #plot_counts(like.obs()) 319 + # Print help if needed and exit
  320 + if need_help:
  321 + print("Usage: example_survey.py [OPTIONS]")
  322 + print(" -h Display this usage message")
  323 + sys.exit()
  324 +
  325 + # Dump header
  326 + print("***********************")
  327 + print("* Simulate CTA survey *")
  328 + print("***********************")
  329 +
  330 + # Remove any existing result files
  331 + list = [glob.glob("*.fits"), glob.glob("*.log"), glob.glob("*.xml")]
  332 + for files in list:
  333 + for file in files:
  334 + os.remove(file)
  335 +
  336 + # Setup single observation survey
  337 + #obs = survey_single()
  338 + obs = survey_gplane()
  339 +
  340 + # Add background model
  341 + obs = add_background_model(obs)
  342 +
  343 + # Simulate events
  344 + print("Simulate events")
  345 + obs = obsutils.sim(obs)
  346 +
  347 + # Make counts map
  348 + print("Make counts map")
  349 + cntmap = obsutils.cntmap(obs)
  350 +
  351 + # Fit observations
  352 + print("Fit observations")
  353 + like = obsutils.fit(obs)
  354 + #print like.opt()
  355 + #print like.obs().models()
  356 +
  357 + # Make model map
  358 + print("Make model map (this step will take some time)")
  359 + modmap = obsutils.modmap(obs)
  360 +
  361 + # Show fit results
  362 + #plot_counts(like.obs())