diff --git a/test/example_survey.py b/test/example_survey.py index 54056b6..dfa6fa4 100755 --- a/test/example_survey.py +++ b/test/example_survey.py @@ -3,7 +3,7 @@ # This script simulates CTA surveys. It can be considered as the master # script for CTA science simulations. Please modify as needed. # -# Copyright (C) 2011 Jurgen Knodlseder +# Copyright (C) 2011-2012 Jurgen Knodlseder # # 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 @@ -25,6 +25,7 @@ from math import * import os import glob import sys +import obsutils # =========== # @@ -47,9 +48,9 @@ def plot_counts(obs): styles = ['b-', 'g-', 'y-', 'n-'] # Dump header - print "" - print "Make plots (using matplotlib):" - print "==============================" + print("") + print("Make plots (using matplotlib):") + print("==============================") # Create figure 1 plt.figure(1,figsize=(12,6)) @@ -63,7 +64,7 @@ def plot_counts(obs): for run in obs: # Get event list - list = cast_GCTAEventList(run.events()) + list = run.events() # Create energy axis ebounds = GEbounds() @@ -75,9 +76,9 @@ def plot_counts(obs): energy = [ebounds.elogmean(i).TeV() for i in range(ebounds.size())] # Create spectrum - print "Extract data:" + print("Extract data:") counts = [0.0 for i in range(ebounds.size())] - print list.size() + print(list.size()) for atom in list: index = ebounds.index(atom.energy()) counts[index] = counts[index] + 1.0 @@ -90,10 +91,10 @@ def plot_counts(obs): #plt.errorbar(energy, counts, error, fmt=None, ecolor='r') # Extract models - print "Extract models:" + print("Extract models:") sum_model = [0.0 for i in range(ebounds.size())] for k, m in enumerate(obs.models()): - print "- "+m.name() + print("- "+m.name()) model = [0.0 for i in range(ebounds.size())] for atom in list: index = ebounds.index(atom.energy()) @@ -120,11 +121,11 @@ def plot_counts(obs): # Loop over observations for run in obs: - # Get event cube - list = cast_GCTAEventList(run.events()) + # Get event list + list = run.events() # Create offset histogram - print "Extract data:" + print("Extract data:") nx = 30 doffset2 = 0.01 offset2 = [(i+0.5)*doffset2 for i in range(nx)] @@ -143,10 +144,10 @@ def plot_counts(obs): #plt.semilogy(offset2, counts, 'ro', label='data') # Extract models - print "Extract models:" + print("Extract models:") sum_model = [0.0 for i in range(nx)] for k, m in enumerate(obs.models()): - print "- "+m.name() + print("- "+m.name()) model = [0.0 for i in range(nx)] for atom in list: off = atom.dir().dist_deg(crab) @@ -170,200 +171,12 @@ def plot_counts(obs): plt.show() except ImportError: - print "Matplotlib is not (correctly) installed on your system." + print("Matplotlib is not (correctly) installed on your system.") # Return return -# ===================== # -# Simulate observations # -# ===================== # -def sim_obs(obs): - """ - Simulate events for all observations in the container. - """ - # Simulate events - sim = ctobssim(obs) - sim.logFileOpen() # We need this to explicitely open the log file in Python mode - sim.run() - - # Extract copy of observations (we need a copy here as sim goes out of - # scope once we leave the function) - obs = sim.obs().copy() - - # Return observations - return obs - - -# ================ # -# Fit observations # -# ================ # -def fit_obs(obs): - """ - Perform maximum likelihood fitting of observations. - """ - # Perform maximum likelihood fitting - like = ctlike(obs) - like.logFileOpen() # We need this to explicitely open the log file in Python mode - like.run() - - # Extract copy of observations (we need a copy here as like goes out of - # scope once we leave the function) - #obs = like.obs().copy() - - #print like.opt() - - # Return observations - return like - - -# ================= # -# Create counts map # -# ================= # -def make_counts_map(obs): - """ - Make counts map by combining the events of all observations. - """ - # Allocate counts map - nxpix = 600 - nypix = 200 - cntmap = GSkymap("CAR", "GAL", 0.0, 0.0, -0.05, 0.05, nxpix, nypix, 1) - - # Fill all observations - for run in obs: - - # Loop over all events - for event in run.events(): - - # Cast to CTA event - event = cast_GCTAEventAtom(event) - - # Determine sky pixel - skydir = event.dir().skydir() - pixel = cntmap.dir2pix(skydir) - - # Set pixel - cntmap[pixel] += 1.0 - - # Save sky map - cntmap.save("survey_cntmap.fits", True) - - # Return counts map - return cntmap - - -# ================ # -# Create model map # -# ================ # -def make_model_map(obs): - """ - Make model map by combining all observations. The map will - be in units of [counts/(sr MeV s)]. - """ - # Allocate model map - nxpix = 600 - nypix = 200 - modmap = GSkymap("CAR", "GAL", 0.0, 0.0, -0.05, 0.05, nxpix, nypix, 1) - - # Set energy and time - energy = GEnergy() - time = GTime() - instdir = GCTAInstDir() - energy.TeV(0.1) - - # Loop over all map pixels - for pixel in range(modmap.npix()): - - # Get sky direction - skydir = modmap.pix2dir(pixel) - instdir.skydir(skydir) - - # Create event atom for map pixel - atom = GCTAEventAtom() - atom.dir(instdir) - atom.energy(energy) - atom.time(time) - - # Initialise model value - value = 0.0 - - # Loop over all observations - for run in obs: - value += obs.models().eval(atom, run) - - # Set map value - modmap[pixel] = value - - # Save sky map - modmap.save("survey_modmap.fits", True) - - # Return model map - return modmap - - -# ======================= # -# Set one CTA observation # -# ======================= # -def set_one_obs(pntdir, tstart=0.0, duration=1800.0, emin=0.1, emax=100.0, rad=3.0, \ - irf="kb_E_50h_v3", caldb="irf"): - """ - Parameters: - pntdir - Pointing direction - Keywords: - tstart - Start time (seconds) - duration - Duration of observation (seconds) - emin - Minimum event energy (TeV) - emax - Maximum event energy (TeV) - rad - ROI radius (deg) - irf - Instrument response function - caldb - Calibration database path - """ - # Allocate CTA observation - obs = GCTAObservation() - - # Set pointing direction - pnt = GCTAPointing() - pnt.dir(pntdir) - obs.pointing(pnt) - - # Set ROI - roi = GCTARoi() - instdir = GCTAInstDir() - instdir.skydir(pntdir) - roi.centre(instdir) - roi.radius(rad) - - # Set GTI - gti = GGti() - start = GTime() - stop = GTime() - start.met(tstart) - stop.met(tstart+duration) - gti.append(start, stop) - - # Set energy boundaries - ebounds = GEbounds() - e_min = GEnergy() - e_max = GEnergy() - e_min.TeV(emin) - e_max.TeV(emax) - ebounds.append(e_min, e_max) - - # Allocate event list - events = GCTAEventList() - events.roi(roi) - events.gti(gti) - events.ebounds(ebounds) - obs.events(events) - - # Set instrument response - obs.response(irf, caldb) - - # Return observation - return obs - - # ======================================== # # Add CTA background model to observations # # ======================================== # @@ -428,7 +241,8 @@ def survey_single(): # Set single pointing at galactic centre pntdir = GSkyDir() pntdir.lb_deg(0.0, 0.0) - run = set_one_obs(pntdir) + #run = set_one_obs(pntdir) + run = obsutils.set(pntdir) obs.append(run) # Define single point source with Crab flux at galactic centre @@ -468,7 +282,8 @@ def survey_gplane(lrange=10, lstep=2): # Set pointing pntdir = GSkyDir() pntdir.lb_deg(l, 0.0) - run = set_one_obs(pntdir) + #run = set_one_obs(pntdir) + run = obsutils.set(pntdir) obs.append(run) # Define single point source with Crab flux at galactic centre @@ -499,7 +314,7 @@ if __name__ == '__main__': need_help = False # Test for command line arguments - print sys.argv[0] + print(sys.argv[0]) if (len(sys.argv) > 1): if sys.argv[1] == "-h": need_help = True @@ -508,14 +323,14 @@ if __name__ == '__main__': # Print help if needed and exit if need_help: - print "Usage: example_survey.py [OPTIONS]" - print " -h Display this usage message" + print("Usage: example_survey.py [OPTIONS]") + print(" -h Display this usage message") sys.exit() # Dump header - print "**********************" - print "* Simulate CTA survey *" - print "**********************" + print("***********************") + print("* Simulate CTA survey *") + print("***********************") # Remove any existing result files list = [glob.glob("*.fits"), glob.glob("*.log"), glob.glob("*.xml")] @@ -524,29 +339,29 @@ if __name__ == '__main__': os.remove(file) # Setup single observation survey - obs = survey_single() + #obs = survey_single() obs = survey_gplane() # Add background model obs = add_background_model(obs) # Simulate events - print "Simulate events" - obs = sim_obs(obs) + print("Simulate events") + obs = obsutils.sim(obs) # Make counts map - print "Make counts map" - cntmap = make_counts_map(obs) + print("Make counts map") + cntmap = obsutils.cntmap(obs) # Fit observations - print "Fit observatins" - like = fit_obs(obs) + print("Fit observatins") + like = obsutils.fit(obs) #print like.opt() #print like.obs().models() # Make model map - print "Make model map (this step will take some time)" - modmap = make_model_map(obs) + print("Make model map (this step will take some time)") + modmap = obsutils.modmap(obs) # Show fit results #plot_counts(like.obs())