first commit of the high level analysis script : analyse, Ebin and UL

sanchez
1 parent 148a99f9
Showing 4 changed files with 814 additions and 38 deletions
pyext/ctools/__init__.py
pyext/ctools/analysisutils.py
pyext/ctools/plotutils.py
pyext/setup.py.in
 from ctools import *
 import obsutils
+import analysisutils
+import plotutils
+# ==========================================================================
+# This script provides a number of functions that are useful for handling
+# CTA observations.
+#
+# Copyright (C) 2011-2014 David Sanchez Michael Mayer Rolf Buehler 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 ctools as ct
+import gammalib as gl
+from math import log10,floor
+import sys
+
+class base(object):
+    def __init__(self):
+        self.classname = self.__class__.__name__
+        self.errorcolor = "\033[31m"#red
+        self.infocolor = "\033[34m"#blue
+        self.warningcolor = "\033[33m"#yellow
+        self.successcolor = "\033[32m"#green
+        self.endcolor = "\033[0m"#reset
+        self.prependfunction = True
+        self.basemembers = ["classname","errorcolor","infocolor","warningcolor","successcolor","endcolor","prependfunction"]
+        
+    def error(self,message, functionname = ""):
+        printstring = ""
+        if functionname == "":
+            printstring = "\n"+self.errorcolor+"*** Error ["+self.classname+"]: "+message+" ***\n"+self.endcolor
+        else:
+            printstring = "\n"+self.errorcolor+"*** Error ["+self.classname+"::"+functionname+"]: "+message+" ***\n"+self.endcolor
+        sys.exit(printstring)
+    
+    def info(self,message,newline=True):
+        printstring = self.infocolor+"["+self.classname+"]: "+message+self.endcolor     
+        if newline:
+            print printstring
+        else:
+            print self.infocolor+message+self.endcolor,
+            sys.stdout.flush()  
+        
+    def warning(self,message,functionname = ""):
+        printstring = ""
+        if functionname == "":
+            printstring = self.warningcolor+"["+self.classname+"] Warning: "+message+self.endcolor
+        else:
+            printstring = self.warningcolor+"["+self.classname+"::"+functionname+"] Warning: "+message+self.endcolor
+        print printstring
+        
+    def success(self,message):
+        printstring = self.successcolor+"["+self.classname+"]: "+message+self.endcolor
+        print printstring
+        
+    def progress(self,message="."):
+        string = self.infocolor+message+self.endcolor
+        print string
+        
+        
+class ResultsSorage(dict,base):
+    """class to store the results in the good format"""
+    def __init__(self,*arg,**kw):
+        super(ResultsSorage, self).__init__(*arg, **kw)
+        base.__init__(self)
+        self["dnde"]      = {}
+        self["ulim_dnde"] = []
+        self["ener"]      = {}
+        self["TS"]        = []
+        self["Iflux"]     = {}
+        
+        self["ener"]["value"]      = []
+        self["ener"]["ed_value"]   = []
+        self["ener"]["eu_value"]   = []
+        self["ener"]["unit"]       = "MeV"
+        
+        self["dnde"]["value"]      = []
+        self["dnde"]["ed_value"]   = []
+        self["dnde"]["eu_value"]   = []
+        self["dnde"]["unit"]       = "ph/cm2/s/MeV"
+        
+        self["Iflux"]["value"]     = []
+        self["Iflux"]["ed_value"]  = []
+        self["Iflux"]["eu_value"]  = []
+        self["Iflux"]["unit"]      = "ph/cm2/s"
+        
+    def Print(self):
+        self.success("Results are :")
+        for k in self.iterkeys():
+            self.info("Key name : "+k)
+            try :
+                for ki in self[k].iterkeys():
+                    print "\t",ki,"\t",self[k][ki]
+            except:
+                print self[k]
+        print 
+        
+class Analyser(base):
+    def __init__(self):
+        super(Analyser,self).__init__()
+        self.info("Creating "+self.classname+" object")
+        self.m_obs = None
+        self.like  = None
+        
+        self.m_evtfile  = "events_selected.fits"        
+        self.m_raw_evt  = "events.fits"
+        self.m_cntfile  = "cntmap.fits"
+        self.m_ebinalg  = "LOG"
+        self.m_eunit    = "TeV"
+        self.m_ra       = 83.63
+        self.m_dec      = 22.01
+        self.m_roi      = 2.
+        self.m_tmin     = 0.
+        self.m_tmax     = 0.
+        self.m_emin     = 0.1#emin
+        self.m_emax     = 10.#emax
+        self.m_enumbins = 10#nbins
+        self.m_usepnt   = True # Use pointing for map centre
+        self.m_nxpix    = 200#npix
+        self.m_nypix    = 200#npix
+        self.m_binsz    = 0.05#binsz
+        self.m_coordsys = "GAL"#coord
+        self.m_proj     = "TAN"#proj
+        self.m_binned   = False
+        
+        self.m_stat     = "POISSON"#stat
+        self.m_caldb    = "dummy"
+        self.m_irf      = "cta_dummy_irf"
+        self.m_edisp    = False
+        self.m_xml      = "$PWD/crab.xml"
+        self.m_xml_out  = "$PWD/crab_results.xml"
+
+
+    def set_obs(self,obs):
+        self.m_obs = obs #set the GObservation container
+
+    def set_energy_boundary(self,emin,emax):
+        self.m_emin = emin
+        self.m_emax = emax
+        ebounds = gl.GEbounds(gl.GEnergy(emin, self.m_eunit), \
+                                gl.GEnergy(emax, self.m_eunit))
+
+        if self.m_obs:
+            for obs in self.m_obs:
+                obs.events().ebounds(ebounds)
+        
+    def ctselect(self):
+        # ctselect application and set parameters
+        self.info("Running ctselect to cut on events")
+        if self.m_obs:
+            filter = ct.ctselect(self.m_obs)
+        else:
+            filter = ct.ctselect()
+        
+        filter["infile"] = self.m_raw_evt               
+        filter["outfile"] = self.m_evtfile  
+        filter["usepnt"].boolean(self.m_usepnt) # Use pointing for map centre
+        filter["ra"]  = self.m_ra
+        filter["dec"]  = self.m_dec
+        filter["rad"]  = self.m_roi
+        filter["tmin"] = self.m_tmin
+        filter["tmax"] = self.m_tmax
+        filter["emin"].real(self.m_emin)
+        filter["emax"].real(self.m_emax)
+        filter["expr"] = ""
+        filter.logFileOpen()
+
+        filter.run()
+        if not(self.m_obs):
+            filter.save()
+
+        if self.m_obs:
+            # Make a deep copy of the observation that will be returned
+            # (the ctbin object will go out of scope one the function is
+            # left)
+            self.m_obs = filter.obs().copy()
+        
+    def ctbin(self,log=False,debug=False):
+        # ctbin application and set parameters
+        self.info("Running ctbin to create count map")
+        if self.m_obs:
+            bin = ct.ctbin(self.m_obs)
+        else:
+            bin = ct.ctbin()
+            bin["evfile"] = self.m_evtfile
+        bin["outfile"] = self.m_cntfile
+        bin["ebinalg"].string(self.m_ebinalg)
+        bin["emin"].real(self.m_emin)
+        bin["emax"].real(self.m_emax)
+        Nbdecade = log10(self.m_emax)-log10(self.m_emin)#Compute the number of decade
+        bin["enumbins"].integer(int(self.m_enumbins*Nbdecade))
+        bin["usepnt"].boolean(self.m_usepnt) # Use pointing for map centre
+        bin["nxpix"].integer(self.m_nxpix)
+        bin["nypix"].integer(self.m_nypix)
+        bin["binsz"].real(self.m_binsz)
+        bin["coordsys"].string(self.m_coordsys)
+        bin["proj"].string(self.m_proj)
+        
+        # Optionally open the log file
+        if log:
+            bin.logFileOpen()
+
+        # Optionally switch-on debugging model
+        if debug:
+            bin["debug"].boolean(True)
+
+        # Run ctbin application. This will loop over all observations in
+        # the container and bin the events in counts maps
+        bin.excute()
+        if self.m_obs:
+            # Make a deep copy of the observation that will be returned
+            # (the ctbin object will go out of scope one the function is
+            # left)
+            self.m_obs = bin.obs().copy()
+            
+    def create_fit(self,log=False,debug=False):
+        # create ctlike instance with given parameters
+        self.info("Fitting Data using ctlike")
+        if self.m_obs:
+            self.like = ct.ctlike(self.m_obs)
+        else:
+            self.like = ct.ctlike()
+            if self.m_binned:
+                self.like["infile"] = self.m_cntfile
+            else:
+                self.like["infile"] = self.m_evtfile
+                
+            self.like["stat"].string(self.m_stat)
+            self.like["caldb"].string(self.m_caldb)
+            self.like["irf"].string(self.m_irf)
+            self.like["srcmdl"] = self.m_xml 
+        self.like["edisp"].boolean(self.m_edisp)
+        self.like["outmdl"] = self.m_xml_out
+
+        # Optionally open the log file
+        if log:
+            self.like.logFileOpen()
+        # Optionally switch-on debugging model
+        if debug:
+            self.like["debug"].boolean(True)
+            
+    def fit(self,log=False,debug=False):
+        if not(self.like):
+            self.warning("ctlike object not created, creating now")
+            self.create_fit(log,debug)
+
+        # Run ctlike application.
+        self.like.run()
+        # Save the results in XML
+        self.like.save()
+        self.success("Fit performed")
+        
+    def PrintResults(self,srcname = ""):
+        self.info("Results of the Fit")
+        for m in self.like.obs().models():
+            if srcname == m.name() or srcname=="":
+                print "Model : "+m.name()
+                print m
+                
+    def GetSrcResuls(self,Res,srcname,E0=1,factor = 1e6):
+
+        results = self.GetSrcParams(srcname)
+
+        Res["ener"]["value"].append(E0*factor)
+        Res["ener"]["ed_value"].append((E0-self.m_emin)*factor)
+        Res["ener"]["eu_value"].append((self.m_emax-E0)*factor)
+        Res["ener"]["unit"] = "MeV"
+        
+        Res["Iflux"]["value"].append(results["flux"])
+        Res["Iflux"]["eu_value"].append(results["eflux"])
+        Res["Iflux"]["ed_value"].append(results["eflux"])
+        # Res["Iflux"]["unit"] = "unit"
+        
+        m = self.GetModel(srcname)
+        Res["dnde"]["value"].append(m.spectral()[0].value())
+        Res["dnde"]["eu_value"].append(m.spectral()[0].error())
+        Res["dnde"]["ed_value"].append(m.spectral()[0].error())
+        # Res["dnde"]["value"].append(m.spectral().eval(gl.GEnergy(E0,self.m_eunit),gl.GTime(0)))
+        # m.spectral().eval_gradients(gl.GEnergy(E0,self.m_eunit),gl.GTime(0))
+        # Assume that par 0 is the norm
+        # Res["dnde"]["eu_value"].append(m.spectral()[0].factor_gradient())
+        # Res["dnde"]["ed_value"].append(m.spectral()[0].factor_gradient())
+        # Res["dnde"]["unit"] = "unit"
+
+        Res["TS"].append(100)
+        # Res["TS"].append(self.like.obs().models()[srcname].ts())
+
+        UL = UpperLimitComputer(self,srcname) 
+        UL.run()
+        Res["ulim_dnde"].append(0)
+        # del UL
+        return Res
+
+    def GetSrcParams(self,srcname):
+        results = {}
+        found = False
+        m = self.GetModel(srcname)
+        results["flux"] = m.spectral().flux(gl.GEnergy(self.m_emin,self.m_eunit),gl.GEnergy(self.m_emax,self.m_eunit))
+        results["eflux"] = m.spectral().eflux(gl.GEnergy(self.m_emin,self.m_eunit),gl.GEnergy(self.m_emax,self.m_eunit))
+        for par in m.spectral():
+            if par.is_free():
+                results[par.name()] = par.value()
+                results["e"+par.name()] = par.error()
+            
+        return results
+        
+    def GetModel(self,srcname):
+        found = False
+        for m in self.like.obs().models():
+            if m.name() == srcname:
+                return m
+        self.error("source "+srcname+" not in the source list")
+            
+    
+def MakeEbin(analyse,nbins,srcname,binned = False):
+    #Store the results
+    Res = ResultsSorage()
+    
+    #check the source name, exit is the srcname is false
+    analyse.GetModel(srcname)
+    
+    #compute the conversion factor to have energy in MeV
+    Energy_list = {"MeV":1,"GeV":1e3,"TeV":1e6} #TODO
+    factor = Energy_list[analyse.m_eunit]
+    
+    for i in xrange(nbins):
+        obs_copy = analyse.like.obs().copy()
+        emin = pow(10,log10(analyse.m_emin) + (log10(analyse.m_emax)-log10(analyse.m_emin))/(nbins)*i)
+        emax = pow(10,log10(analyse.m_emin) + (log10(analyse.m_emax)-log10(analyse.m_emin))/(nbins)*(i+1))
+        E0 = pow(10, (log10(emin) + log10(emax)) / 2)
+        print "Making energy bin between %2.1e"%(emin)+analyse.m_eunit+" and %2.1e"%(emax)+analyse.m_eunit
+
+        # change the PivotEnergy to E0, froze all spectral parameters but the prefactor
+        m = obs_copy.models()[srcname]
+        for par in m.spectral():
+            if par.name() == "PivotEnergy":
+                par.value(E0*factor)
+            # if not(par.name() == "Prefactor") and par.is_free(): # assume that all model have a prefactor
+                # par.fix()
+            # if par.name() == "Prefactor":
+                # Set the norm to a better guess
+                # new_val = m.spectral().eval(gl.GEnergy(E0,analyse.m_eunit),gl.GTime(0))
+                # m.spectral()["Prefactor"].value(new_val)
+                            
+        ana_bin = Analyser()
+
+        #copy the obs object
+        ana_bin.set_obs(obs_copy)
+
+        # Set energy boundaries
+        ana_bin.set_energy_boundary(emin,emax)
+        ana_bin.ctselect()
+
+        # if binned:
+            # ana_bin.ctbin()
+
+        # Create the fit object
+        ana_bin.create_fit(log=False,debug=False)
+        
+        # Run ctlike application.
+        ana_bin.like.run()
+
+        #Retrive the results in the good format
+        Res = ana_bin.GetSrcResuls(Res,srcname,E0,factor)
+        
+        del ana_bin
+    return Res
+        
+class UpperLimitComputer(base):
+    def __init__(self,analyser,srcname,parname = "Prefactor"):
+        super(UpperLimitComputer,self).__init__()
+        self.info("Creating "+self.classname+" object")
+        # self.like = analyser.like.copy()
+        obs = analyser.like.obs().copy()
+        self.like = ct.ctlike(obs)
+        self.succes = False
+
+        # get spectral model
+        self.model = self.like.obs().models()[srcname]
+        self.spec = self.model.spectral()
+        self.parname = parname
+        self.bestloglike = analyser.like.opt().value()
+        self.dloglike = 3.84/2.0 #95% CL now
+
+    def _bisec(self,a,b,tol = 1e-12):
+        """iterative function to find the roots"""
+        if (self._func(a)==0.):
+            return a
+        elif (self._func(b)==0):
+            return b
+        elif (self._func(b)*self._func(a)>0.):
+            self.warning("Finding root for upper limit calculation failed")
+            return 0
+        mid = (a+b)/2
+        err = abs(a-b)/2
+        fm = self._func(mid)
+        if (err<tol or self._func(mid)==0.):
+            self.succes = True
+            return mid
+        elif (fm*self._func(a)<0):
+            b = mid
+        else :
+            a = mid
+        return self._bisec(a,b)
+           
+    def _func(self,value): 
+        """function to be minimized"""
+        # get parmeter value
+        val = value*self.spec[self.parname].scale()
+        
+        # ensure value to be above minimum boundary
+        if val <= self.spec[self.parname].min():
+            return 0.0-self.self.bestloglike-self.dloglike
+        
+        # set value
+        self.spec[self.parname].value(val)
+        
+        # Recalculate likelihood
+        self.like.run()
+        logl = self.like.opt().value()   
+        
+        # Return likelihood difference
+        return (logl-self.bestloglike-self.dloglike)
+    
+    def run(self):
+        
+        # get spectral parameters and set start value
+        value = self.spec[self.parname].value()
+        error = self.spec[self.parname].error()
+        scale = self.spec[self.parname].scale()
+        print value," ",error," ",scale
+
+        #Fix all parameters
+        for model in self.like.obs().models():
+            for par in model:
+                par.fix()
+                
+        
+        # set start value to 2 sigma above fit value
+        self.startpar = value / scale + 2* error/scale
+
+        # set start (a) and (b) depend on the CL asked
+        a= value / scale +  (self.dloglike/4.)*error/scale
+        b=value / scale + (self.dloglike*2)* error/scale
+        print a," ",b," ",self._func(a)," ",self._func(b)
+        ul = self._bisec(a,b)
+        
+        self.ulimit = ul*self.spec[self.parname].scale()
+        # Test for convergence
+        if self.succes:
+            self.success("Upper limit of parameter \""+self.parname +"\"  determined to "+str(self.ulimit))
@@ -2,7 +2,7 @@
 # This script provides a number of functions that are useful for handling
 # CTA observations.
 #
-# Copyright (C) 2014 Rolf Buehler
+# Copyright (C) 2011-2014 David Sanchez Michael Mayer Rolf Buehler 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
@@ -18,54 +18,368 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
 # ==========================================================================
-import ctools
-import gammalib
+import ctools as ct
+import gammalib as gl
+from math import log10,pow,sqrt
+import os
+
 try:
         import matplotlib.pyplot as plt
+        import matplotlib.gridspec as gridspec
         has_matplotlib = True
 except:
         has_matplotlib = False
+
 try:
         import aplpy
         has_aplpy = True
 except:
         has_aplpy = False
+class options(object):
+    def __init__(self):
+        self.fmt         = 'o'
+        self.color       = BLUE
+        self.linecolor   = 'black'
+        self.markercolor = self.color
+        self.tslim       = 4.0
+        self.elinewidth  = 2.0
+        self.linewidth   = 2.0
+        self.linestyle   = '--'
+        self.capsize     = 5.0
+        self.markersize  = 5.0
+        self.limlength   = 0.4 
+        self.label       = "_nolegend_"
+
+import analysisutils 
+
+class SpectralBase(object):
+    def __init__(self,sed=True,sed_factor=1.):
+        self.m_sed_factor = sed_factor
+        self.m_sed        = sed
+        
+    def _convertSED(self,tab,energy):
+        """convert the a table storing dN/dE in sed with right units"""
+
+        if len(tab) != len(energy):
+            self.error("input table and energy table does not match")
+            return
+        if self.m_sed:
+            for i in xrange(len(tab)):
+                tab[i] *= energy[i]**2*self.m_sed_factor
+        return tab
+
+BLUE="#337EF5"
+GREY='grey'
+GREEN = '#21851D'
+RED = '#BD0B0B'
+YELLOW = '#D4CF44'
+ORANGE = '#D19534'
+
+
+class residuals(analysisutils.base,options,SpectralBase):
+    def __init__(self,spectrum,datapoint,sed,sed_factor,eunit="TeV"):
+        super(residuals,self).__init__()
+        options.__init__(self)
+        SpectralBase.__init__(self,sed,sed_factor)
+
+        # load the gammalib spectrum object
+        self.spectrum  = spectrum
+        self.SetDataPoint(datapoint)
+        
+        self.xerrors = True
+        self.refUnit = eunit     # reference unit for conversion
+        
+    def SetDataPoint(self,datapoint):
+        self.m_point_info = datapoint
+        self.Npt = len(self.m_point_info["dnde"]["value"])
+
+    def _plt_residuals(self):
+        th_val = []
+        index  = []
+        ed_y   = []
+        eu_y   = []
+        y      = []
+        
+        Energy_list = {"MeV":1,"GeV":1e-3,"TeV":1e-6} #TODO
+        for i in range(self.Npt):
+            if self.m_point_info["TS"][i] >= self.tslim:
+                th_val.append(self.spectrum.eval(gl.GEnergy(self.m_point_info["ener"]["value"][i],self.m_point_info["ener"]["unit"]),gl.GTime(0)))
+                index.append(i)
+            else:
+                th_val.append(0)
+        if self.m_sed:
+            self._convertSED(th_val,self.m_point_info["ener"])
+        for i in index:
+                y.append((self.m_point_info["dnde"]["value"][i] - th_val[i]) / th_val[i])
+                ed_y.append(self.m_point_info["dnde"]["ed_value"][i]/th_val[i])
+                eu_y.append(self.m_point_info["dnde"]["eu_value"][i]/th_val[i])
+                
+        if self.xerrors:
+            plt.errorbar(self.m_point_info["ener"]["value"], y, xerr=[self.m_point_info["ener"]["ed_value"],self.m_point_info["ener"]["ed_value"]],yerr=[ed_y,eu_y],fmt=self.fmt,color=self.markercolor,elinewidth=self.elinewidth) 
+        else:
+            plt.errorbar(self.m_point_info["ener"]["value"], y, yerr=[ed_y,eu_y],fmt=self.fmt,color=self.markercolor,elinewidth=self.elinewidth) 
+        
+        plt.axhline(0.0,color=self.linecolor,lw=self.linewidth,ls=self.linestyle)
+
+
+class ulimgraph(analysisutils.base,options,SpectralBase):
+    def __init__(self,datapoint,sed,sed_factor,eunit='TeV'):
+        super(ulimgraph,self).__init__()
+        #~ super(ulimgraph,self).__init__()
+        options.__init__(self)
+        SpectralBase.__init__(self,sed,sed_factor)
+        
+        self.xerrors=True  
+        self.nlims = -1
+        self.refUnit = eunit
+        self.SetDataPoint(datapoint)
+    
+    def SetDataPoint(self,datapoint):
+        self.m_point_info = datapoint
+        self.Npt = len(self.m_point_info["dnde"]["value"])
+       
+    def _plt_points(self):
+        results = analysisutils.ResultsSorage()
+        
+        Energy_list = {"MeV":1,"GeV":1e-3,"TeV":1e-6} #TODO
+        for i in range(self.Npt):
+
+            self.m_point_info["ener"]["value"][i]*=Energy_list[self.refUnit]
+            self.m_point_info["ener"]["ed_value"][i]*=Energy_list[self.refUnit]
+            self.m_point_info["ener"]["eu_value"][i]*=Energy_list[self.refUnit]
+            self.m_point_info["ener"]["unit"] = self.refUnit
+
+            if self.m_point_info["TS"][i] >= self.tslim:
+
+                if self.m_point_info["dnde"]["ed_value"][i]>self.m_point_info["dnde"]["value"][i]:
+                    self.warning("\tflux error: "+str(self.m_point_info["dnde"]["ed_value"][i])+" is larger than flux value "+str(self.m_point_info["dnde"]["ed_value"][i]))
+                    self.warning("\tTS value is however TS="+str(self.m_point_info["TS"][i])+", weird!")
+            else :
+                self.m_point_info["dnde"]["ed_value"][i] = 0
+                self.m_point_info["dnde"]["eu_value"][i] = 0
+            
+        if self.xerrors:
+            plt.errorbar(self.m_point_info["ener"]["value"], self.m_point_info["dnde"]["value"], xerr=[self.m_point_info["ener"]["ed_value"], self.m_point_info["ener"]["eu_value"]],yerr=[self.m_point_info["dnde"]["ed_value"], self.m_point_info["dnde"]["eu_value"]],fmt=self.fmt,color=self.markercolor,elinewidth=self.elinewidth) 
+        else:
+            plt.errorbar(self.m_point_info["ener"]["value"], self.m_point_info["dnde"]["value"],yerr=[results["dnde"]["ed_value"],results["dnde"]["eu_value"]],fmt=self.fmt,color=self.markercolor,elinewidth=self.elinewidth) 
+        
+    def _plt_limits(self):
+        ed_flux = []
+        eu_flux = []
+        flux    = []
+        ener    = []
+        
+        drawn_lims = 0
+        Energy_list = {"MeV":1,"GeV":1e-3,"TeV":1e-6} #TODO
+        for i in range(self.Npt):
+            if self.m_point_info["TS"][i] < self.tslim:
+                self.info("Plot an upper limit for E = "+self.m_point_info["ener"]["value"][i]+" with a TS = "+self.m_point_info["TS"][i])
+                if drawn_lims >= self.nlims and self.nlims!=-1:
+                    break
+
+                ed_flux.append(self.limlength*self.m_point_info["ulim_dnde"][i])
+                eu_flux.append(0.0)
+                flux.append(self.m_point_info["ulim_dnde"][i])
+                drawn_lims+=1
+                
+        # if self.xerrors:
+            # plt.errorbar(self.m_point_info["ener"]["value"], self.m_point_info["dnde"]["value"], xerr=[self.m_point_info["ener"]["ed_value"], self.m_point_info["ener"]["eu_value"]],yerr=[self.m_point_info["dnde"]["ed_value"],self.m_point_info["dnde"]["eu_value"]],fmt=self.fmt,markersize=0.0,elinewidth=self.elinewidth,lolims=True,capsize=self.capsize,label='_nolegend_',mfc = self.color,ecolor=self.color,ms =self.markersize)      
+        # else:
+        plt.errorbar(ener, flux, yerr=[ed_flux, eu_flux],fmt=self.fmt,markersize=0.0,elinewidth=self.elinewidth,lolims=True,capsize=self.capsize,label='_nolegend_',mfc = self.color,ecolor=self.color,ms =self.markersize)      
+            
+
+class SpectrumPlotter(analysisutils.base,SpectralBase):
+    """ Class to compute and plot the spectra (i.e best fit model and butterfly).
+        if data points are provided, only change applied  is energy unit convertion to match the butterfly one
+        (i.e. no dN/dE to SED conversion).
+        It is currently assumed that the unit is MeV for the energy of the data point. this will be change 
+        once the I/O container format will be defined."""
+    def __init__(self,srcname,like,datapoint,emin=0.1,emax=100,energy = "TeV",SED = True ,npt = 50):
+        super(SpectrumPlotter,self).__init__()
+
+        self.m_ebound   = [emin,emax]
+        self.m_name     = srcname
+        self.m_like     = like
+        self.m_spectral = like.obs().models()[srcname].spectral()
+        self.m_energy   = []
+        self.m_flux     = []
+        self.m_error    = []
+        self.m_but      = []
+        self.m_enebut   = []
+        self.m_npt      = npt
+        self.m_covar    = gl.GMatrix(1,1)
+        self.covar()
+        #~ this part take care of the style : SED or differential flux
+        #~ in the case of SED, the y-unit is erg/cm2/s
+        self.m_sed = SED
+        self.eunit = energy
+        self.m_sed_factor = 1.
+        self._validateUnits()
+        SpectralBase.__init__(self,self.m_sed,self.m_sed_factor)
+
+        #Assume MeV for the energy #TODO
+        self.points = ulimgraph(datapoint,self.m_sed,gl.MeV2erg,energy)
+        self.residuals = residuals(self.m_spectral,datapoint,self.m_sed,gl.MeV2erg,energy)
+
+
+    def _validateUnits(self):
+        Energy_list = {"MeV":1,"GeV":1e3,"TeV":1e6}
+        if not(Energy_list.has_key(self.eunit)):
+            self.warning("Change energy unit to TeV")
+            self.eunit = "TeV"
+        
+        if self.m_sed:
+            self.m_sed_factor = gl.MeV2erg*Energy_list[self.eunit]**2
+
+    def covar(self):
+        # Get covariance matrix if fit has converged
+        try:
+            fullcovar = self.m_like.obs().function().curvature().invert()
+        except:
+            self.warning("Covariance matrix not determined successfully")
+            return
+        #~ get the matrix element of the source of interest
+        #~ 1) get the right element indices
+        idx = 0
+        par_index_map = []
+        for m in self.m_like.obs().models() :
+            if m.name()!= self.m_name:
+                idx += m.size()
+                continue
+            idx += m.spatial().size()
+            for par in m.spectral():
+                if par.is_free():
+                    par_index_map.append(idx)
+                    idx += 1
+
+        #~ 2) get the elemets and store them in a matrix
+        self.m_covar = gl.GMatrix(len(par_index_map),len(par_index_map))
+        i = 0
+        for xpar in par_index_map:
+            j = 0
+            for ypar in par_index_map:
+                self.m_covar[i,j] = fullcovar[xpar,ypar]
+                j += 1
+            i += 1
+        self.success("Covariance Matrix succesfuly computed")
+
+    def _makespectrum(self):
+        for i in xrange(self.m_npt):
+            #Compute the energy
+            lene = log10(self.m_ebound[0])+log10(self.m_ebound[1])*i/(self.m_npt-1)
+            self.m_energy.append(pow(10,lene))
+            #Compute the flux
+            self.m_flux.append(self.m_spectral.eval(gl.GEnergy(pow(10,lene),self.eunit),gl.GTime(0)))
+
+            #Compute the gradients
+            if self.m_covar.size() != 1:
+                self.m_spectral.eval_gradients(gl.GEnergy(pow(10,lene), self.eunit),gl.GTime(0))
+                #store them into a GVector
+                Derivative = gl.GVector(self.m_covar.size()/2)
+                j = 0
+                for par in self.m_spectral:
+                    if par.is_free():
+                        Derivative[j] = par.factor_gradient()
+                    j += 1
+                #~ computed the error
+                self.m_error.append(sqrt(Derivative * (self.m_covar*Derivative)))
+
+        #~ convert the flux in sed if asked
+        self.m_flux = self._convertSED(self.m_flux,self.m_energy)
+        if self.m_covar.size() != 1:
+            self.m_error = self._convertSED(self.m_error,self.m_energy)
+            self._makebutterfly() # make the butterfly
+        self.success("Spectrum computed")
+
+        
+    def _makebutterfly(self):
+        """make the butterfly by appending element in a table"""
+        for i in xrange(self.m_npt):
+            self.m_but.append(self.m_flux[i]+self.m_error[i])
+            self.m_enebut.append(self.m_energy[i])
+        for i in xrange(self.m_npt):
+            idx = self.m_npt-i-1
+            self.m_but.append(self.m_flux[idx]-self.m_error[idx])
+            self.m_enebut.append(self.m_energy[idx])
+        self.m_but.append(self.m_but[0])
+        self.m_enebut.append(self.m_enebut[0])
+
+    def write(self):
+        self.warning("write function not yet implemented")
+
+    def draw(self):
+        """plot the results, provide a file named ptfile for the flux points"""
+        gs = gridspec.GridSpec(2, 1,height_ratios=[3,1])
+        self._makespectrum()
+        self.write()
+
+        if not(has_matplotlib):
+            self.warning("matplotlib module not found, can draw")
+            return
+
+        plt.figure("results",figsize=(12,7),edgecolor="w")
+
+        sb1 = plt.subplot(gs[0])
+        plt.loglog()
+        plt.xlabel("Energy ["+ self.eunit+"]",fontsize=15)
+        
+        if self.m_sed:
+            plt.ylabel("E$^{2}$ dN/dE [erg cm$^{2}$ s$^{-1}$]",fontsize=15)
+        else:
+            plt.ylabel("dN/dE [cm$^{-2}$ s$^{-1}$ "+ self.eunit+"$^{-1}$]",fontsize=15)
+
+        plt.plot(self.m_energy,self.m_flux)
+        plt.plot(self.m_enebut,self.m_but)
+        
+        #~ plot the data points on top of the buterfly
+        self.points._plt_points()
+        self.points._plt_limits()
+
+        #~ plot the residual points
+        sb2 = plt.subplot(gs[1])     
+        plt.xscale('log')     
+        self.residuals._plt_residuals()
+
+        plt.show()
+
 # ============= #
 # Plot skymaps  #
 # ============= #
-def plotskymaps(modelmap, countmap):
-    """
-    Plot the model map.
-
-    Parameters:    
-     modname - model input name
-     evtname - events input name
-    Keywords:
-     None
-    """
-    # Load model map
-    fig = plt.figure()
-    f1 = aplpy.FITSFigure(countmap, figure=fig, subplot=[0.1,0.1,0.35,0.8])
-    f1.set_tick_labels_font(size='small')
-    f1.set_axis_labels_font(size='small')
-    f1.show_colorscale()
-
-    f1.tick_labels.set_yformat('dd:mm:ss')
-    f1.tick_labels.set_xformat('hh:mm:ss')
-    f1.axis_labels.set_xtext('Right Ascension (J2000)')
-    f1.axis_labels.set_ytext('Declination (J2000)')
-    f1.ticks.set_length(10.5, minor_factor=0.5)
-
-    f2 = aplpy.FITSFigure(modelmap, figure=fig, subplot=[0.5,0.1,0.35,0.8])
-    f2.set_tick_labels_font(size='small')
-    f2.set_axis_labels_font(size='small')
-    f2.show_colorscale()
-
-    f2.hide_yaxis_label()
-    f2.hide_ytick_labels()
-
-    fig.canvas.draw()
-    fig.show()
-    
+class MapPlotter(analysisutils.base):
+    def __init__(self,modelmap, countmap):
+        super(MapPlotter,self).__init__()
+        self.modelmap = modelmap
+        self.countmap = countmap
+        
+        if not(has_aplpy) :
+            self.error("Module aplypy needed")
+        if not(has_matplotlib) :
+            self.error("Module matplotlib needed")
+            
+    def draw(self):
+        """ Plot the model map.    """
+        # Load model map
+        fig = plt.figure()
+        f1 = aplpy.FITSFigure(self.countmap, figure=fig, subplot=[0.1,0.1,0.35,0.8])
+        f1.set_tick_labels_font(size='small')
+        f1.set_axis_labels_font(size='small')
+        f1.show_colorscale()
+
+        f1.tick_labels.set_yformat('dd:mm:ss')
+        f1.tick_labels.set_xformat('hh:mm:ss')
+        f1.axis_labels.set_xtext('Right Ascension (J2000)')
+        f1.axis_labels.set_ytext('Declination (J2000)')
+        f1.ticks.set_length(10.5, minor_factor=0.5)
+
+        f2 = aplpy.FITSFigure(self.modelmap, figure=fig, subplot=[0.5,0.1,0.35,0.8])
+        f2.set_tick_labels_font(size='small')
+        f2.set_axis_labels_font(size='small')
+        f2.show_colorscale()
+
+        f2.hide_yaxis_label()
+        f2.hide_ytick_labels()
+
+        fig.canvas.draw()
+        fig.show()
@@ -121,7 +121,7 @@ setup(name=&#39;@PACKAGE_NAME@&#39;,
       author='Juergen Knoedlseder',
       author_email='@PACKAGE_BUGREPORT@',
       url='http://gammalib.sourceforge.net/',
-      py_modules=['ctools.ctools', 'ctools.obsutils'],
+      py_modules=['ctools.ctools', 'ctools.obsutils','ctools.plotutils','ctools.analysisutils'],
       ext_modules=[Extension('ctools._ctools', ['ctools/ctools_wrap.cpp'],
                              include_dirs=include_dirs,
                              library_dirs=library_dirs,