Alexander Ziegler / ctools

Commit e285d731fa0bab30ce5ebc04237ff920687edb1d

Authored by Jürgen Knödlseder
1 parent 5003e285

Use obsutils module and adapt to Python 3.

Inline Side-by-side
Showing 1 changed file with 35 additions and 220 deletions
test/example_survey.py
  View file @e285d73
... ... @@ -3,7 +3,7 @@
3 3 # This script simulates CTA surveys. It can be considered as the master
4 4 # script for CTA science simulations. Please modify as needed.
5 5 #
6   -# Copyright (C) 2011 Jurgen Knodlseder
  6 +# Copyright (C) 2011-2012 Jurgen Knodlseder
7 7 #
8 8 # This program is free software: you can redistribute it and/or modify
9 9 # it under the terms of the GNU General Public License as published by
... ... @@ -25,6 +25,7 @@ from math import *
25 25 import os
26 26 import glob
27 27 import sys
  28 +import obsutils
28 29  
29 30  
30 31 # =========== #
... ... @@ -47,9 +48,9 @@ def plot_counts(obs):
47 48 styles = ['b-', 'g-', 'y-', 'n-']
48 49  
49 50 # Dump header
50   - print ""
51   - print "Make plots (using matplotlib):"
52   - print "=============================="
  51 + print("")
  52 + print("Make plots (using matplotlib):")
  53 + print("==============================")
53 54  
54 55 # Create figure 1
55 56 plt.figure(1,figsize=(12,6))
... ... @@ -63,7 +64,7 @@ def plot_counts(obs):
63 64 for run in obs:
64 65  
65 66 # Get event list
66   - list = cast_GCTAEventList(run.events())
  67 + list = run.events()
67 68  
68 69 # Create energy axis
69 70 ebounds = GEbounds()
... ... @@ -75,9 +76,9 @@ def plot_counts(obs):
75 76 energy = [ebounds.elogmean(i).TeV() for i in range(ebounds.size())]
76 77  
77 78 # Create spectrum
78   - print "Extract data:"
  79 + print("Extract data:")
79 80 counts = [0.0 for i in range(ebounds.size())]
80   - print list.size()
  81 + print(list.size())
81 82 for atom in list:
82 83 index = ebounds.index(atom.energy())
83 84 counts[index] = counts[index] + 1.0
... ... @@ -90,10 +91,10 @@ def plot_counts(obs):
90 91 #plt.errorbar(energy, counts, error, fmt=None, ecolor='r')
91 92  
92 93 # Extract models
93   - print "Extract models:"
  94 + print("Extract models:")
94 95 sum_model = [0.0 for i in range(ebounds.size())]
95 96 for k, m in enumerate(obs.models()):
96   - print "- "+m.name()
  97 + print("- "+m.name())
97 98 model = [0.0 for i in range(ebounds.size())]
98 99 for atom in list:
99 100 index = ebounds.index(atom.energy())
... ... @@ -120,11 +121,11 @@ def plot_counts(obs):
120 121 # Loop over observations
121 122 for run in obs:
122 123  
123   - # Get event cube
124   - list = cast_GCTAEventList(run.events())
  124 + # Get event list
  125 + list = run.events()
125 126  
126 127 # Create offset histogram
127   - print "Extract data:"
  128 + print("Extract data:")
128 129 nx = 30
129 130 doffset2 = 0.01
130 131 offset2 = [(i+0.5)*doffset2 for i in range(nx)]
... ... @@ -143,10 +144,10 @@ def plot_counts(obs):
143 144 #plt.semilogy(offset2, counts, 'ro', label='data')
144 145  
145 146 # Extract models
146   - print "Extract models:"
  147 + print("Extract models:")
147 148 sum_model = [0.0 for i in range(nx)]
148 149 for k, m in enumerate(obs.models()):
149   - print "- "+m.name()
  150 + print("- "+m.name())
150 151 model = [0.0 for i in range(nx)]
151 152 for atom in list:
152 153 off = atom.dir().dist_deg(crab)
... ... @@ -170,200 +171,12 @@ def plot_counts(obs):
170 171 plt.show()
171 172  
172 173 except ImportError:
173   - print "Matplotlib is not (correctly) installed on your system."
  174 + print("Matplotlib is not (correctly) installed on your system.")
174 175  
175 176 # Return
176 177 return
177 178  
178 179  
179   -# ===================== #
180   -# Simulate observations #
181   -# ===================== #
182   -def sim_obs(obs):
183   - """
184   - Simulate events for all observations in the container.
185   - """
186   - # Simulate events
187   - sim = ctobssim(obs)
188   - sim.logFileOpen() # We need this to explicitely open the log file in Python mode
189   - sim.run()
190   -
191   - # Extract copy of observations (we need a copy here as sim goes out of
192   - # scope once we leave the function)
193   - obs = sim.obs().copy()
194   -
195   - # Return observations
196   - return obs
197   -
198   -
199   -# ================ #
200   -# Fit observations #
201   -# ================ #
202   -def fit_obs(obs):
203   - """
204   - Perform maximum likelihood fitting of observations.
205   - """
206   - # Perform maximum likelihood fitting
207   - like = ctlike(obs)
208   - like.logFileOpen() # We need this to explicitely open the log file in Python mode
209   - like.run()
210   -
211   - # Extract copy of observations (we need a copy here as like goes out of
212   - # scope once we leave the function)
213   - #obs = like.obs().copy()
214   -
215   - #print like.opt()
216   -
217   - # Return observations
218   - return like
219   -
220   -
221   -# ================= #
222   -# Create counts map #
223   -# ================= #
224   -def make_counts_map(obs):
225   - """
226   - Make counts map by combining the events of all observations.
227   - """
228   - # Allocate counts map
229   - nxpix = 600
230   - nypix = 200
231   - cntmap = GSkymap("CAR", "GAL", 0.0, 0.0, -0.05, 0.05, nxpix, nypix, 1)
232   -
233   - # Fill all observations
234   - for run in obs:
235   -
236   - # Loop over all events
237   - for event in run.events():
238   -
239   - # Cast to CTA event
240   - event = cast_GCTAEventAtom(event)
241   -
242   - # Determine sky pixel
243   - skydir = event.dir().skydir()
244   - pixel = cntmap.dir2pix(skydir)
245   -
246   - # Set pixel
247   - cntmap[pixel] += 1.0
248   -
249   - # Save sky map
250   - cntmap.save("survey_cntmap.fits", True)
251   -
252   - # Return counts map
253   - return cntmap
254   -
255   -
256   -# ================ #
257   -# Create model map #
258   -# ================ #
259   -def make_model_map(obs):
260   - """
261   - Make model map by combining all observations. The map will
262   - be in units of [counts/(sr MeV s)].
263   - """
264   - # Allocate model map
265   - nxpix = 600
266   - nypix = 200
267   - modmap = GSkymap("CAR", "GAL", 0.0, 0.0, -0.05, 0.05, nxpix, nypix, 1)
268   -
269   - # Set energy and time
270   - energy = GEnergy()
271   - time = GTime()
272   - instdir = GCTAInstDir()
273   - energy.TeV(0.1)
274   -
275   - # Loop over all map pixels
276   - for pixel in range(modmap.npix()):
277   -
278   - # Get sky direction
279   - skydir = modmap.pix2dir(pixel)
280   - instdir.skydir(skydir)
281   -
282   - # Create event atom for map pixel
283   - atom = GCTAEventAtom()
284   - atom.dir(instdir)
285   - atom.energy(energy)
286   - atom.time(time)
287   -
288   - # Initialise model value
289   - value = 0.0
290   -
291   - # Loop over all observations
292   - for run in obs:
293   - value += obs.models().eval(atom, run)
294   -
295   - # Set map value
296   - modmap[pixel] = value
297   -
298   - # Save sky map
299   - modmap.save("survey_modmap.fits", True)
300   -
301   - # Return model map
302   - return modmap
303   -
304   -
305   -# ======================= #
306   -# Set one CTA observation #
307   -# ======================= #
308   -def set_one_obs(pntdir, tstart=0.0, duration=1800.0, emin=0.1, emax=100.0, rad=3.0, \
309   - irf="kb_E_50h_v3", caldb="irf"):
310   - """
311   - Parameters:
312   - pntdir - Pointing direction
313   - Keywords:
314   - tstart - Start time (seconds)
315   - duration - Duration of observation (seconds)
316   - emin - Minimum event energy (TeV)
317   - emax - Maximum event energy (TeV)
318   - rad - ROI radius (deg)
319   - irf - Instrument response function
320   - caldb - Calibration database path
321   - """
322   - # Allocate CTA observation
323   - obs = GCTAObservation()
324   -
325   - # Set pointing direction
326   - pnt = GCTAPointing()
327   - pnt.dir(pntdir)
328   - obs.pointing(pnt)
329   -
330   - # Set ROI
331   - roi = GCTARoi()
332   - instdir = GCTAInstDir()
333   - instdir.skydir(pntdir)
334   - roi.centre(instdir)
335   - roi.radius(rad)
336   -
337   - # Set GTI
338   - gti = GGti()
339   - start = GTime()
340   - stop = GTime()
341   - start.met(tstart)
342   - stop.met(tstart+duration)
343   - gti.append(start, stop)
344   -
345   - # Set energy boundaries
346   - ebounds = GEbounds()
347   - e_min = GEnergy()
348   - e_max = GEnergy()
349   - e_min.TeV(emin)
350   - e_max.TeV(emax)
351   - ebounds.append(e_min, e_max)
352   -
353   - # Allocate event list
354   - events = GCTAEventList()
355   - events.roi(roi)
356   - events.gti(gti)
357   - events.ebounds(ebounds)
358   - obs.events(events)
359   -
360   - # Set instrument response
361   - obs.response(irf, caldb)
362   -
363   - # Return observation
364   - return obs
365   -
366   -
367 180 # ======================================== #
368 181 # Add CTA background model to observations #
369 182 # ======================================== #
... ... @@ -428,7 +241,8 @@ def survey_single():
428 241 # Set single pointing at galactic centre
429 242 pntdir = GSkyDir()
430 243 pntdir.lb_deg(0.0, 0.0)
431   - run = set_one_obs(pntdir)
  244 + #run = set_one_obs(pntdir)
  245 + run = obsutils.set(pntdir)
432 246 obs.append(run)
433 247  
434 248 # Define single point source with Crab flux at galactic centre
... ... @@ -468,7 +282,8 @@ def survey_gplane(lrange=10, lstep=2):
468 282 # Set pointing
469 283 pntdir = GSkyDir()
470 284 pntdir.lb_deg(l, 0.0)
471   - run = set_one_obs(pntdir)
  285 + #run = set_one_obs(pntdir)
  286 + run = obsutils.set(pntdir)
472 287 obs.append(run)
473 288  
474 289 # Define single point source with Crab flux at galactic centre
... ... @@ -499,7 +314,7 @@ if __name__ == '__main__':
499 314 need_help = False
500 315  
501 316 # Test for command line arguments
502   - print sys.argv[0]
  317 + print(sys.argv[0])
503 318 if (len(sys.argv) > 1):
504 319 if sys.argv[1] == "-h":
505 320 need_help = True
... ... @@ -508,14 +323,14 @@ if __name__ == '__main__':
508 323  
509 324 # Print help if needed and exit
510 325 if need_help:
511   - print "Usage: example_survey.py [OPTIONS]"
512   - print " -h Display this usage message"
  326 + print("Usage: example_survey.py [OPTIONS]")
  327 + print(" -h Display this usage message")
513 328 sys.exit()
514 329  
515 330 # Dump header
516   - print "**********************"
517   - print "* Simulate CTA survey *"
518   - print "**********************"
  331 + print("***********************")
  332 + print("* Simulate CTA survey *")
  333 + print("***********************")
519 334  
520 335 # Remove any existing result files
521 336 list = [glob.glob("*.fits"), glob.glob("*.log"), glob.glob("*.xml")]
... ... @@ -524,29 +339,29 @@ if __name__ == '__main__':
524 339 os.remove(file)
525 340  
526 341 # Setup single observation survey
527   - obs = survey_single()
  342 + #obs = survey_single()
528 343 obs = survey_gplane()
529 344  
530 345 # Add background model
531 346 obs = add_background_model(obs)
532 347  
533 348 # Simulate events
534   - print "Simulate events"
535   - obs = sim_obs(obs)
  349 + print("Simulate events")
  350 + obs = obsutils.sim(obs)
536 351  
537 352 # Make counts map
538   - print "Make counts map"
539   - cntmap = make_counts_map(obs)
  353 + print("Make counts map")
  354 + cntmap = obsutils.cntmap(obs)
540 355  
541 356 # Fit observations
542   - print "Fit observatins"
543   - like = fit_obs(obs)
  357 + print("Fit observatins")
  358 + like = obsutils.fit(obs)
544 359 #print like.opt()
545 360 #print like.obs().models()
546 361  
547 362 # Make model map
548   - print "Make model map (this step will take some time)"
549   - modmap = make_model_map(obs)
  363 + print("Make model map (this step will take some time)")
  364 + modmap = obsutils.modmap(obs)
550 365  
551 366 # Show fit results
552 367 #plot_counts(like.obs())
... ...