GammaLib

#! /usr/bin/env python

# ==========================================================================

# This script simulates a number of CTA observations and fits a model to the

# simulated events. To make this script work, make sure that the $GAMMALIB

# environment variable points to the directory whene gammalib and ctools are

# installed. Here an example of how to configure gammalib and ctools:

#

# Suppose you have two directories, one holding a copy of gammalib and one

# holding a copy of ctools. Then install both components in a specific

# install directory by executing the following commands:

#

# cd gammalib

# ./configure --prefix=/home/knodlseder/testinstall

# make install

# cd ..

# cd ctools

# ./configure --prefix=/home/knodlseder/testinstall

# make install

#

# Now you need to configure gammalib and ctools:

#

# export GAMMALIB=/home/knodlseder/testinstall

# source $GAMMALIB/bin/gammalib-init.sh

#

# Copyright (C) 2012 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/>.

#

# ==========================================================================

from ctools import *

from gammalib import *

import obsutils

import sys

import time

# ========= #

# cta class #

# ========= #

class cta:

        """

        Implements methods for observation definition, simulations and model

        fitting.

        """

        def __init__(self, *argv):

                """

                Constructor.

                """

                # Initialise members

                self.m_obs      = GObservations()

                self.m_srcmdl   = "${GAMMALIB}/share/models/crab.xml"

                self.m_caldb    = "${GAMMALIB}/share/caldb/cta"

                self.m_irf      = "kb_E_50h_v3"

                self.m_ra       = 83.6331

                self.m_dec      = 22.0145

                self.m_emin     = 0.1

                self.m_emax     = 100.0

                self.m_duration = 1800.0

                self.m_deadc    = 0.95

                self.m_rad      = 5.0

                # Return

                return

        def __del__(self):

                """

                Destructor.

                """

                # Return

                return

        def simulate(self, number):

                """

                Simulation CTA observations.

                Parameters:

                 number - Number of observations in container.

                """

                # Set seed value

                seed = 0

                # Clear observation container

                self.m_obs.clear()

                # Set pointing direction

                pntdir = GSkyDir()

                pntdir.radec_deg(self.m_ra, self.m_dec)

                # Create CTA observation

                run = obsutils.set(pntdir, caldb=self.m_caldb, irf=self.m_irf, \

                                   duration=self.m_duration, deadc=self.m_deadc, \

                           emin=self.m_emin, emax=self.m_emax, rad=self.m_rad)

                # Append observation to container

                for i in range(number):

                        self.m_obs.append(run)

                # Append source model

                self.m_obs.models(GModels(self.m_srcmdl))

                # Simulate observations

                self.m_obs = obsutils.sim(self.m_obs, seed=seed, debug=False)

                # Return

                return

        def fit(self):

                """

                Fit model.

                """

                # Save start time

                start = time.clock()

                # Fit model

                like = obsutils.fit(self.m_obs)

                # Compute elapsed time

                elapsed = (time.clock() - start)

                print "Elapsed time: ", elapsed

                # Return

                return

# ======================== #

# Main routine entry point #

# ======================== #

if __name__ == '__main__':

        """

        This script simulates a number of CTA observations and fits a model

        to the simulated events. The command line parameter specifies the

        number of CTA observations that should be simulated. By default, 100

        observations are simulated.

        """

        # Get input arguments

        usage = "openmp [number of observations]"

        if len(sys.argv) < 1 or len(sys.argv) > 2:

                print usage

                sys.exit()

        # Set number of observations

        if len(sys.argv) > 1:

                number = int(sys.argv[1])

        else:

                number = 100

        # Create instance of test class

        obs = cta()

        # Simulate observations

        obs.simulate(number)

        # Do model fitting

        obs.fit()