ctools

#! /usr/bin/env python

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

# This script generates the pull distribution for all model parameters.

#

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

#

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

import gammalib as gl

import ctools as ct

import sys

import math

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

# csbutterfly class #

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

class csbutterfly(gl.GApplication):

    """

    This class implements an the computation of a butterfly. 

    It derives from the GammaLib::GApplication class which 

    provides support for parameter files, command line arguments,

    and logging. In that way the Python script behaves just as a 

    regular ctool.

    """

    def __init__(self, *argv):

        """

        Constructor.

        """

        # Set name

        self.name    = "csbutterfly"

        self.version = "0.0.1"

        # Initialise some members

        if isinstance(argv[0],gl.GObservations):

            self.obs = argv[0]

        else:      

            self.obs      = gl.GObservations()

        self.source = "CrabNebula"

        self.m_sigma = 1 

        self.m_emin = 0.0 #TeV

        self.m_emax = 0.0 #TeV 

        self.m_nbins = 100

        # Make sure that parfile exists

        file = self.parfile()

        # Initialise application

        if len(argv) == 1:

            gl.GApplication.__init__(self, self.name, self.version)

        elif len(argv) ==2:

            gl.GApplication.__init__(self, self.name, self.version, argv[1:])

        else:

            raise TypeError("Invalid number of arguments given.")

        # Set logger properties

        self.log_header()

        self.log.date(True)

        # Return

        return

    def __del__(self):

        """

        Destructor.

        """

        #  Write separator into logger

        if self.logTerse():

            self.log("\n")

        # Return

        return

    def parfile(self):

        """

        Check if parfile exists. If parfile does not exist then create a

        default parfile. This kluge avoids shipping the cscript with a parfile.

        """

        # Set parfile name

        parfile = self.name+".par"

        try:

            pars = gl.GApplicationPars(parfile)

        except:

            # Signal if parfile was not found

            sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")

            # Create default parfile

            pars = gl.GApplicationPars()

            pars.append(gl.GApplicationPar("infile", "f", "a", "events.fits","","","observation definition file"))

            pars.append(gl.GApplicationPar("stat","s","h","POISSON","","","Optimization statistics"))

            pars.append(gl.GApplicationPar("refit",  "b", "h", "no","","", "Do refitting?"))

            pars.append(gl.GApplicationPar("srcmdl","f","a","$CTOOLS/share/models/crab.xml","","","Source model"))

            pars.append(gl.GApplicationPar("caldb","s","h","","","","Calibration database"))

            pars.append(gl.GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function"))

            pars.append(gl.GApplicationPar("edisp","b","h","no","","","Apply energy dispersion?"))

            pars.append(gl.GApplicationPar("logfile", "f", "h", "csbutterfly.log","","", "Log filename"))

            pars.append(gl.GApplicationPar("outfile", "f", "a", "CrabNebula.bf","","", "Outfile name"))

            pars.append(gl.GApplicationPar("Emin", "r", "a", "0.1","","TeV", "Minimum Energy"))

            pars.append(gl.GApplicationPar("Emax", "r", "a", "100","","TeV", "Maximum Energy"))

            pars.append(gl.GApplicationPar("nbins", "i", "h", "100","","", "Butterfly spacing"))

            pars.append(gl.GApplicationPar("sigma", "f", "h", "1","","", "Butterfly spacing"))

            pars.append_standard()

            pars.save(parfile)

        # Return

        return

    def get_parameters(self):

        """

        Get parameters from parfile and setup the observation.

        """

        # Setup obs if it isn't already set

        if self.obs.size() == 0:

            infile = self["infile"].filename()

            try:      

                self.obs = gl.GObservations(infile)

            except:

                obs0 = gl.GCTAObservation(infile)

                self.obs.append(obs0)

            mdlfile = self["srcmdl"].filename()

            self.obs.models(mdlfile)

            self.m_caldb    = self["caldb"].string()

            self.m_irf      = self["irf"].string()

        else:

            self.m_caldb = ""

            self.m_irf = ""

        # get hidden parameters

        self.m_stat = self["stat"].string()

        self.m_refit = self["refit"].boolean()

        self.m_edisp = self["edisp"].boolean()

        self.m_chatter = self["chatter"].integer()

        self.m_debug = self["debug"].boolean()

        self.m_clobber = self["clobber"].boolean()

        self.m_mode = self["mode"].string()

        self.m_logfile = self["logfile"].filename()

        if self.m_emin == 0.0 or self.m_emax == 0.0:

            self.m_emin = self["Emin"].real()

            self.m_emax = self["Emax"].real()

        self.m_nbins = self["nbins"].integer()

    def execute(self):

        """

        Execute the script.

        """

        # Run the script

        self.run()

        # Return

        return

    def obs(self):

        """

        Return observations

        """

        return self.obs.clone()    

    def run(self):

        """

        Run the fit.

        """

        # Switch screen logging on in debug mode

        if self.logDebug():

            self.log.cout(True)

        # Get parameters

        self.get_parameters()

        #  Write input parameters into logger

        if self.logTerse():

            self.log_parameters()

            self.log("\n")

        # Write observation into logger

        if self.logTerse():

            self.log("\n")

            self.log.header1("Observation")

            self.log(str(self.obs))

            self.log("\n")

        # Write header

        if self.logTerse():

            self.log("\n")

            self.log.header1("Generate butterfly")

        # Check if observation container is already optimised

        # might add a flag in the future for that

        if self.obs.function().value() == 0.0:

            print "Optimising observation container"

            like = ct.ctlike(self.obs)

            like.run()

            self.covar = like.obs().function().curvature().invert()

            self.obs.models(like.obs().models().clone())

        else:

            self.covar = self.obs.function().curvature().invert()

        try:

            m = self.obs.models()[self.source]

        except RuntimeError:

            sys.exit("Source \""+self.source+"\" not found in model container",sys._getframe().f_code.co_name)

        # Initialise energies in log binning

        energies = []

        step = (math.log10(self.m_emax)-math.log10(self.m_emin) )/ (self.m_nbins-1)

        for ibin in range(self.m_nbins):

            loge = math.log10(self.m_emin)+ibin*step

            energies.append(pow(10,loge))

        fluxes = []

        errors = []

        t=gl.GTime()

        self.fluxes = []

        self.errors = []

        for energy in energies:

            npars = self.obs.models().npars()

            grad = gl.GVector(npars)

            e = gl.GEnergy(energy,"TeV") 

            k=0

            value=0            

            for m in self.obs.models():                

                if m.name() == self.source: 

                    k+=m.spatial().size()

                    spec = m.spectral()

                    value = spec.eval_gradients(e,t)

                    for i in range(spec.size()):  

                        grad[k] = spec[i].gradient()

                        k+=1

                    k+=1

                else:

                    k+=m.size()

            fluxes.append(value)

            vect = self.covar * grad

            err2 = grad * vect

            err = math.sqrt(err2)            

            errors.append(err)

#             if self.verbose:

#                 self.success("E = "+str(e)+", Flux = "+str("%.3e" % value)+" +- "+str("%.3e" % err)+" 1/cm2/s/MeV")

        energies_MeV = []

        for e in energies:

            energies_MeV.append(e*1e6)

        self.save(self.outfile,energies_MeV,fluxes,errors)

        print "Finished: output written to "+self.outfile

    def save(self,outfile,energies,fluxes,errors):

        f=open(outfile,"w")

        f.write("# 1 energy [MeV] \n")

        f.write("# 2 value [1/cm2/s/MeV] \n")

        for i in range(len(energies)):

            y= fluxes[i]+errors[i]

            f.write(str(energies[i])+" "+str(y)+" \n")

        for i in range(len(energies)):

            y = fluxes[-(i+1)]-errors[-(i+1)]

            if y <1e-35:

                y=1e-35

            f.write(str(energies[-(i+1)])+" "+str(y)+" \n")            

        f.write(str(energies[0])+" "+str((fluxes[0]+errors[0]))+" \n")

        f.close()   

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

# Main routine entry point #

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

if __name__ == '__main__':

    """

    Generates butterfly for a source model.

    """

    # Create instance of application

    app = csbutterfly(sys.argv)

    # Open logfile

    app.logFileOpen()

    # Execute application

    app.execute()