From 4e5640586d55059a016cb66b8fcf4b77b866e6f6 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Tue, 30 Jun 2015 23:07:39 +0200 Subject: [PATCH] Remove analysisutils and plotutils from distribution. --- pyext/ctools/analysisutils.py | 575 ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- pyext/ctools/plotutils.py | 421 ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- 2 files changed, 0 insertions(+), 996 deletions(-) delete mode 100644 pyext/ctools/analysisutils.py delete mode 100644 pyext/ctools/plotutils.py 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 . -# -# ========================================================================== -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. -# -# ========================================================================== -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() -- libgit2 0.21.4