diff --git a/pyext/ctools/analysisutils.py b/pyext/ctools/analysisutils.py
deleted file mode 100644
index be634c7..0000000
--- a/pyext/ctools/analysisutils.py
+++ /dev/null
@@ -1,575 +0,0 @@
-# ==========================================================================
-# This script provides a number of functions that are useful for CTA data
-# analysis.
-#
-# Copyright (C) 2014 David Sanchez
-#
-# 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 sqrt,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["time"]      = {}
-        self["ener"]      = {}
-        self["TS"]        = []
-        self["Iflux"]     = {}
-        
-        self["time"]["tmin_value"] = []
-        self["time"]["tmax_value"] = []
-        self["time"]["unit"]       = "sec"
-        
-        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_ebinalg  = "LOG"
-        self.m_eunit    = "TeV"
-        self.m_enumbins = 10#nbins
-        self.m_usepnt   = True # Use pointing for map centre
-        self.m_coordsys = "CEL"#coord
-        self.m_proj     = "TAN"#proj
-        self.m_binned   = False
-        self.m_edisp    = False
-        self.m_stat     = "POISSON"#stat
-
-        self.m_source_set = False
-        self.m_files_set  = False
-        self.m_time_set   = False
-        self.m_energy_set = False
-        self.m_irfs_set   = False
-        self.m_roi_set    = False
-        self.m_xml_set    = False
-        
-    def set_xml(self,xml):
-        self.m_xml     = xml
-        ind = xml.find(".xml")
-        if not ind==-1:
-            self.m_xml_out = xml[:xml.find(".xml")]+"_results"+xml[xml.find(".xml"):]
-        self.m_xml_set = True
-        
-    def set_roi(self,roi):
-        self.m_roi      = float(roi)
-        self.m_nxpix    = int(self.m_roi*100) #nxpix
-        self.m_nypix    = int(self.m_roi*100) #nypix
-        self.m_binsz    = 0.05                #binsz
-        self.m_roi_set  = True
-
-    def set_irfs(self,caldb,irf):
-        self.m_caldb    = caldb
-        self.m_irf      = irf
-        self.m_irfs_set = True
-
-    def set_Fitsfiles(self,evtfile,tag="CTA_analysis"):
-        self.m_evtfile   = "event_selected_"+tag+".fits"    
-        self.m_raw_evt   = evtfile
-        self.m_cntfile   = "countmap_"+tag+".fits"
-        self.m_files_set = True
-
-    def set_source(self,name,ra,dec):
-        self.m_ra         = ra
-        self.m_dec        = dec
-        self.m_name       = name
-        self.m_source_set = True
-
-    def validate(self):
-        if not self.m_time_set :
-            self.error("Time range of the observations not set")
-
-        if not self.m_files_set :
-            self.error("FITS files of the observations not set")
-
-        if not self.m_source_set :
-            self.error("Source not set")
-
-        if not self.m_energy_set :
-            self.error("Energy range of the observations not set")
-            
-        if not self.m_irfs_set :
-            self.error("IRFs to use not set")
-
-        if not self.m_roi_set :
-            self.error("ROI to use not set")
-
-        if not self.m_xml_set :
-            self.error("XML files for the sky model not set")
-
-    def copy(self,analyse,tag="CTA_analysis_copy"):
-        self.set_xml(analyse.m_xml)
-        self.set_irfs(analyse.m_caldb,analyse.m_irf)
-        self.set_roi(analyse.m_roi)
-        self.set_time_boundary(analyse.m_tmin,analyse.m_tmax)
-        self.set_energy_boundary(analyse.m_emin,analyse.m_emax)
-        self.set_source(analyse.m_name,analyse.m_ra,analyse.m_dec)
-        self.set_Fitsfiles(analyse.m_evtfile,tag)
-
-
-    def set_obs(self,obs):
-        self.m_obs = obs #set the GObservation container
-
-    def set_time_boundary(self,tmin,tmax):
-        """Set the start and stop time of the observation"""
-        self.m_tmin        = float(tmin)
-        self.m_tmax        = float(tmax)
-        self.m_time_set    = True
-
-    def set_energy_boundary(self,emin,emax):
-        """Set the energy bounds of the observation"""
-        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)
-        self.m_energy_set = True
-        
-    def ctselect(self):
-        # ctselect application and set parameters
-        self.info("Running ctselect to cut on events")
-        self.validate()
-        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.validate()
-        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"] = self.m_ebinalg
-        bin["emin"] = self.m_emin
-        bin["emax"] = self.m_emax
-        Nbdecade = log10(self.m_emax)-log10(self.m_emin)#Compute the number of decade
-        bin["enumbins"] = int(self.m_enumbins*Nbdecade)
-        bin["usepnt"] = self.m_usepnt # Use pointing for map centre
-        bin["nxpix"] = self.m_nxpix
-        bin["nypix"] = self.m_nypix
-        bin["binsz"] = self.m_binsz
-        bin["coordsys"] = self.m_coordsys
-        bin["proj"] = self.m_proj
-        
-        # Optionally open the log file
-        if log:
-            bin.logFileOpen()
-
-        # Optionally switch-on debugging model
-        if debug:
-            bin["debug"] = True
-
-        # Run ctbin application. This will loop over all observations in
-        # the container and bin the events in counts maps
-        bin.execute()
-        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):
-        self.validate()
-        # 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"] = self.m_stat
-            self.like["caldb"] = self.m_caldb
-            self.like["irf"] = self.m_irf
-        self.like["srcmdl"] = self.m_xml 
-        self.like["edisp"] = 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"] = True
-            
-    def fit(self,log=False,debug=False):
-        self.validate()
-        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=None,factor = 1e6,tmin=None,tmax=None):
-        """function to get and store the results of a fit"""
-        results = self.GetSrcParams(srcname)
-
-        if 1:
-        # #Time
-        # if (tmin is not None) or (tmax is not None) :
-            # Res["time"]["tmin_value"].append(tmin)
-            # Res["time"]["tmax_value"].append(tmax)
-            # Res["time"]["unit"] = "sec"
-        # else:
-            Res["time"]["tmin_value"].append(self.m_tmin)
-            Res["time"]["tmax_value"].append(self.m_tmax)
-            Res["time"]["unit"] = "sec"
-
-        #Energy
-        if (E0 is not None):
-            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"
-        else:
-            Res["ener"]["value"].append(sqrt(self.m_emin*self.m_emax))
-            Res["ener"]["ed_value"].append(-self.m_emin+sqrt(self.m_emin*self.m_emax))
-            Res["ener"]["eu_value"].append(self.m_emax-sqrt(self.m_emin*self.m_emax))
-            Res["ener"]["unit"] = self.m_eunit
-
-        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["TS"].append(100)
-
-        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
-        m = obs_copy.models()[srcname]
-        for par in m.spectral():
-            if par.name() == "PivotEnergy":
-                par.value(E0*factor)
-                            
-        ana_bin = Analyser()
-
-        #copy the obs object
-        ana_bin.set_obs(obs_copy)
-
-        #copy the analyser parameters
-        ana_bin.copy(analyse)
-
-        # 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
-    
-def MakeLC(analyse,nbins,srcname,binned = False):
-    #Store the results
-    Res = ResultsSorage()
-
-    #check the source name, exit is the srcname is false
-    analyse.GetModel(srcname)
-    
-    for i in xrange(nbins):
-        obs_copy = analyse.like.obs().copy()
-        tmin = analyse.m_tmin + (analyse.m_tmax-analyse.m_tmin)/(nbins)*i
-        tmax = analyse.m_tmin + (analyse.m_tmax-analyse.m_tmin)/(nbins)*(i+1)
-        print("Making time bin between %2.1e"%(tmin)+" and %2.1e"%(tmax))
-                            
-        ana_bin = Analyser()
-
-        #copy the obs object
-        ana_bin.set_obs(obs_copy)
-        
-        #copy the analyser parameters
-        ana_bin.copy(analyse)
-
-        # Set energy boundaries
-        ana_bin.set_time_boundary(tmin,tmax)
-        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,tmin=tmin,tmax=tmax)
-        
-        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()
-
-        #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
-        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))
diff --git a/pyext/ctools/plotutils.py b/pyext/ctools/plotutils.py
deleted file mode 100644
index f2ada6b..0000000
--- a/pyext/ctools/plotutils.py
+++ /dev/null
@@ -1,421 +0,0 @@
-# ==========================================================================
-# This script provides a number of functions that are useful for plotting.
-#
-# Copyright (C) 2014 David Sanchez
-#
-# 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,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 LightCurvePlotter(analysisutils.base,options):
-    """ Class to plot the Light curve"""
-    def __init__(self,srcname,datapoint):
-        super(LightCurvePlotter,self).__init__()
-        options.__init__(self)
-        self.m_name = srcname
-        self.m_data = datapoint #data point in a Results Storage class
-
-    def draw(self):
-        if not(has_matplotlib):
-            self.warning("matplotlib module not found, can draw")
-            return
-
-        time     = []
-        dtime    = []
-        flux     = []
-        dflux_ed = []
-        dflux_eu = []
-        for i in range(len(self.m_data["time"]["tmin_value"])):
-            time.append((self.m_data["time"]["tmin_value"][i]+self.m_data["time"]["tmax_value"][i])/2)
-            dtime.append((self.m_data["time"]["tmax_value"][i]-self.m_data["time"]["tmin_value"][i])/2)
-            flux.append(self.m_data["Iflux"]["value"][i])
-            dflux_eu.append(self.m_data["Iflux"]["eu_value"][i])
-            dflux_ed.append(self.m_data["Iflux"]["ed_value"][i])
-
-        plt.figure("results",figsize=(12,7),edgecolor="w")
-        unit = self.m_data["time"]["unit"]
-        plt.xlabel("Time ["+ unit+"]",fontsize=15)
-        unit = self.m_data["Iflux"]["unit"]
-        plt.ylabel("E$^{2}$ dN/dE ["+unit+"]",fontsize=15)
-
-        print(time)
-        print(dtime)
-        print(flux)
-        print(dflux_ed)
-        print(dflux_eu)
-        plt.errorbar(time,flux,xerr=[dtime,dtime],yerr=[dflux_ed,dflux_eu],fmt=self.fmt,color=self.markercolor,elinewidth=self.elinewidth) 
-        plt.show()
-
-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  #
-# ============= #
-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()