From d11303337da842bf84a741ca82bff8d1645690d6 Mon Sep 17 00:00:00 2001 From: sanchez Date: Mon, 1 Dec 2014 16:33:25 +0100 Subject: [PATCH] use the new ctbutterfly + some doc (need more) --- pyext/ctools/analysisutils.py | 197 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------ pyext/ctools/plotutils.py | 131 ++++++++++++++++++++++++++++++++++++++++++++++++++--------------------------------------------------------------------------------- 2 files changed, 187 insertions(+), 141 deletions(-) diff --git a/pyext/ctools/analysisutils.py b/pyext/ctools/analysisutils.py index 748e900..72fb572 100644 --- a/pyext/ctools/analysisutils.py +++ b/pyext/ctools/analysisutils.py @@ -21,9 +21,10 @@ import ctools as ct import gammalib as gl from math import sqrt,log10,floor -import sys +import sys,string class base(object): + '''This class is used to defined the colors of the print''' def __init__(self): self.classname = self.__class__.__name__ self.errorcolor = "\033[31m"#red @@ -35,6 +36,7 @@ class base(object): self.basemembers = ["classname","errorcolor","infocolor","warningcolor","successcolor","endcolor","prependfunction"] def error(self,message, functionname = ""): + '''Print the error message''' printstring = "" if functionname == "": printstring = "\n"+self.errorcolor+"*** Error ["+self.classname+"]: "+message+" ***\n"+self.endcolor @@ -43,6 +45,7 @@ class base(object): sys.exit(printstring) def info(self,message,newline=True): + '''Print the information message''' printstring = self.infocolor+"["+self.classname+"]: "+message+self.endcolor if newline: print(printstring) @@ -50,7 +53,8 @@ class base(object): print(self.infocolor+message+self.endcolor) sys.stdout.flush() - def warning(self,message,functionname = ""): + def warning(self,message,functionname = ""):* + '''Print the warning message''' printstring = "" if functionname == "": printstring = self.warningcolor+"["+self.classname+"] Warning: "+message+self.endcolor @@ -59,19 +63,21 @@ class base(object): print(printstring) def success(self,message): + '''Print the success message''' printstring = self.successcolor+"["+self.classname+"]: "+message+self.endcolor print(printstring) def progress(self,message="."): + '''Print the progess message''' string = self.infocolor+message+self.endcolor print(string) - -class ResultsSorage(dict,base): +class ResultsStorage(dict,base): """class to store the results in the good format""" def __init__(self,*arg,**kw): - super(ResultsSorage, self).__init__(*arg, **kw) + super(ResultsStorage, self).__init__(*arg, **kw) base.__init__(self) + # Can store the differential flux, upper limit, time, energy and so one self["dnde"] = {} self["ulim_dnde"] = [] self["time"] = {} @@ -79,26 +85,31 @@ class ResultsSorage(dict,base): self["TS"] = [] self["Iflux"] = {} + #store the time : Min, Max and units of the variable self["time"]["tmin_value"] = [] self["time"]["tmax_value"] = [] self["time"]["unit"] = "sec" + #store the energy : Value min,max (ed_* and eu_*) and units of the variable self["ener"]["value"] = [] self["ener"]["ed_value"] = [] self["ener"]["eu_value"] = [] self["ener"]["unit"] = "MeV" + #store the differential flux : Value min,max (ed_* and eu_*) and units of the variable self["dnde"]["value"] = [] self["dnde"]["ed_value"] = [] self["dnde"]["eu_value"] = [] self["dnde"]["unit"] = "ph/cm2/s/MeV" + #store the integral flux : Value min,max (ed_* and eu_*) and units of the variable self["Iflux"]["value"] = [] self["Iflux"]["ed_value"] = [] self["Iflux"]["eu_value"] = [] self["Iflux"]["unit"] = "ph/cm2/s" def Print(self): + """ Print the value stored in the class""" self.success("Results are :") for k in self.iterkeys(): self.info("Key name : "+k) @@ -110,12 +121,17 @@ class ResultsSorage(dict,base): print() class Analyser(base): + '''Class to run the analysis of the data. Currently the parameters are hard coded but this will change once + the method to feed the class will be defined + The classe can work with GObservation as input (simulation) or fits files (simulation+real data)''' def __init__(self): super(Analyser,self).__init__() self.info("Creating "+self.classname+" object") + #Observation and ctlike object self.m_obs = None self.like = None + #parameters self.m_ebinalg = "LOG" self.m_eunit = "TeV" self.m_enumbins = 10#nbins @@ -126,6 +142,8 @@ class Analyser(base): self.m_edisp = False self.m_stat = "POISSON"#stat + # booleans to ensure that input files and parameters have been defined + # set to false at the begining and turn to true while a file/parameter is defined. self.m_source_set = False self.m_files_set = False self.m_time_set = False @@ -135,6 +153,7 @@ class Analyser(base): self.m_xml_set = False def set_xml(self,xml): + '''Set the model XML file''' self.m_xml = xml ind = xml.find(".xml") if not ind==-1: @@ -142,6 +161,7 @@ class Analyser(base): self.m_xml_set = True def set_roi(self,roi): + '''Set the ROI parameter and compute the number of pixel, bin size of the map''' self.m_roi = float(roi) self.m_nxpix = int(self.m_roi*100) #nxpix self.m_nypix = int(self.m_roi*100) #nypix @@ -149,23 +169,45 @@ class Analyser(base): self.m_roi_set = True def set_irfs(self,caldb,irf): + '''Set the IRFs definition''' self.m_caldb = caldb self.m_irf = irf self.m_irfs_set = True def set_Fitsfiles(self,evtfile,tag="CTA_analysis"): + '''Set the input and output fits files''' 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): + '''Set the source name and the RA/DEC position of the source''' self.m_ra = ra self.m_dec = dec self.m_name = name self.m_source_set = True + 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 validate(self): + '''Check is all the parameters has been set''' if not self.m_time_set : self.error("Time range of the observations not set") @@ -188,6 +230,7 @@ class Analyser(base): self.error("XML files for the sky model not set") def copy(self,analyse,tag="CTA_analysis_copy"): + '''Make a copy of the class''' self.set_xml(analyse.m_xml) self.set_irfs(analyse.m_caldb,analyse.m_irf) self.set_roi(analyse.m_roi) @@ -198,35 +241,20 @@ class Analyser(base): def set_obs(self,obs): - self.m_obs = obs #set the GObservation container + '''set the GObservation container''' + self.m_obs = obs - 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 + '''ctselect application and set parameters''' self.info("Running ctselect to cut on events") - self.validate() + self.validate()# check the parameters + if self.m_obs: filter = ct.ctselect(self.m_obs) else: filter = ct.ctselect() + #set the parameters for running the tool filter["infile"] = self.m_raw_evt filter["outfile"] = self.m_evtfile filter["usepnt"].boolean(self.m_usepnt) # Use pointing for map centre @@ -241,6 +269,7 @@ class Analyser(base): filter.logFileOpen() filter.run() + # if no Gobservation have been defined, then save the fits file. if not(self.m_obs): filter.save() @@ -251,14 +280,17 @@ class Analyser(base): self.m_obs = filter.obs().copy() def ctbin(self,log=False,debug=False): - # ctbin application and set parameters - self.validate() + '''ctbin application and set parameters''' self.info("Running ctbin to create count map") + self.validate()# check the parameters + if self.m_obs: bin = ct.ctbin(self.m_obs) else: bin = ct.ctbin() bin["evfile"] = self.m_evtfile + + #set the parameters for running the tool bin["outfile"] = self.m_cntfile bin["ebinalg"].string(self.m_ebinalg) bin["emin"].real(self.m_emin) @@ -290,14 +322,17 @@ class Analyser(base): self.m_obs = bin.obs().copy() def create_fit(self,log=False,debug=False): - self.validate() + '''Create the like object from ctlike''' # create ctlike instance with given parameters self.info("Fitting Data using ctlike") + self.validate() + + #set the parameters of the tools if self.m_obs: self.like = ct.ctlike(self.m_obs) else: self.like = ct.ctlike() - if self.m_binned: + if self.m_binned#binned or unbinned analysis self.like["infile"] = self.m_cntfile else: self.like["infile"] = self.m_evtfile @@ -312,10 +347,50 @@ class Analyser(base): # Optionally open the log file if log: self.like.logFileOpen() + # Optionally switch-on debugging model if debug: self.like["debug"].boolean(True) + def compute_butterfly(self,name,log=False,debug=False): + '''create ctlike instance with given parameters''' + self.info("Computing the data with ctbutterfly") + self.validate() + + if self.m_obs: + butterfly = ct.ctbutterfly(self.m_obs) + else:#parameters of the tools + butterfly = ct.ctbutterfly() + butterfly["infile"] = self.m_evtfile + butterfly['srcname'] = name + butterfly["emin"].real(self.m_emin) + butterfly["emax"].real(self.m_emax) + butterfly["caldb"].string(self.m_caldb) + butterfly["irf"].string(self.m_irf) + butterfly["srcmdl"] = self.m_xml + filename = name+"_butterfly.txt" + butterfly["outfile"] = filename + + # Optionally open the log file + if log: + butterfly.logFileOpen() + # Optionally switch-on debugging model + if debug: + butterfly["debug"].boolean(True) + + butterfly.run() #Compute the butterfly + butterfly.save() # save the txt files + + #Store the results in the Storage class + self.m_butterfly = ResultsStorage() + data = open(filename,'r').readlines() + for line in data: + words = string.split(line) + self.m_butterfly["ener"]["value"].append(float(words[0])) + self.m_butterfly["dnde"]["value"].append(float(words[1])) + self.m_butterfly["dnde"]["ed_value"].append(float(words[1]) - float(words[2])) + self.m_butterfly["dnde"]["eu_value"].append(float(words[1]) + float(words[2])) + def fit(self,log=False,debug=False): self.validate() if not(self.like): @@ -329,81 +404,82 @@ class Analyser(base): self.success("Fit performed") def PrintResults(self,srcname = ""): + ''' Print results of the fit ''' 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""" + def GetSrcResuls(self,Res,srcname,E0=None,factor = 1e6): + """function to get and store the results of a fit. Store the results in the Storage class element""" 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): + #Store the time + Res["time"]["tmin_value"].append(self.m_tmin) + Res["time"]["tmax_value"].append(self.m_tmax) + Res["time"]["unit"] = "sec" + + #Store the energy + if (E0 is not None):#if a 'pivot' energy is given 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: + else:#if not, the geometric mean energy is used 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 + #save the flux values 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()) + #save the differential flux in a src model independent maner + srcmodel = self.GetModel(srcname) + Res["dnde"]["value"].append(srcmodel.spectral()[0].value()) + Res["dnde"]["eu_value"].append(srcmodel.spectral()[0].error()) + Res["dnde"]["ed_value"].append(srcmodel.spectral()[0].error()) - Res["TS"].append(100) + Res["TS"].append(100) ## TODO use the real TS value - UL = UpperLimitComputer(self,srcname) + UL = UpperLimitComputer(self,srcname) #Compute the UL. To be updated UL.run() Res["ulim_dnde"].append(0) - # del UL + return Res def GetSrcParams(self,srcname): + '''Retrive the source parameters from the like object and store it in a dictonnary''' results = {} found = False m = self.GetModel(srcname) + + #get the integral flux 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)) + + # Get free parameters 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): + ''' Get the spectral parameters of the sources stored in the like object''' 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() + Res = ResultsStorage() #check the source name, exit is the srcname is false analyse.GetModel(srcname) @@ -454,7 +530,7 @@ def MakeEbin(analyse,nbins,srcname,binned = False): def MakeLC(analyse,nbins,srcname,binned = False): #Store the results - Res = ResultsSorage() + Res = ResultsStorage() #check the source name, exit is the srcname is false analyse.GetModel(srcname) @@ -493,10 +569,11 @@ def MakeLC(analyse,nbins,srcname,binned = False): return Res class UpperLimitComputer(base): + '''This class compute the ul and will be removed''' 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 diff --git a/pyext/ctools/plotutils.py b/pyext/ctools/plotutils.py index f2ada6b..fa5c7bc 100644 --- a/pyext/ctools/plotutils.py +++ b/pyext/ctools/plotutils.py @@ -56,7 +56,7 @@ 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""" @@ -68,6 +68,26 @@ class SpectralBase(object): tab[i] *= energy[i]**2*self.m_sed_factor return tab + def _convertE(self,energy,unit,newunit): + """convert the a energy which has units 'unit' to "newunit'""" + factor = 1 + if unit == "MeV" and newunit == "GeV" : + factor = 1e-3 + if unit == "MeV" and newunit == "TeV" : + factor = 1e-6 + if unit == "GeV" and newunit == "TeV" : + factor = 1e-3 + if unit == "GeV" and newunit == "MeV" : + factor = 1e3 + if unit == "TeV" and newunit == "GeV" : + factor = 1e3 + if unit == "TeV" and newunit == "MeV" : + factor = 1e6 + for i in xrange(len(energy)): + energy[i] *= factor + return energy + + BLUE="#337EF5" GREY='grey' GREEN = '#21851D' @@ -103,6 +123,7 @@ class residuals(analysisutils.base,options,SpectralBase): 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: + #TODO remove the like 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: @@ -139,7 +160,7 @@ class ulimgraph(analysisutils.base,options,SpectralBase): self.Npt = len(self.m_point_info["dnde"]["value"]) def _plt_points(self): - results = analysisutils.ResultsSorage() + results = analysisutils.ResultsStorage() Energy_list = {"MeV":1,"GeV":1e-3,"TeV":1e-6} #TODO for i in range(self.Npt): @@ -232,22 +253,19 @@ class SpectrumPlotter(analysisutils.base,SpectralBase): 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): + once the I/O container format will be defined. + like is used to compute the residual. if like = None, no residual will be computed TODO + butterfly is a results Storage object read from the file produced by ctbutterfly + datapoint is the results Storage object for the datapoints.""" + def __init__(self,srcname,like, butterfly,datapoint,emin=0.1,emax=100,energy = "TeV",SED = True): super(SpectrumPlotter,self).__init__() + self.m_ebound = [emin,emax] + self.m_name = srcname + self.m_butterfly = butterfly + self.m_spectral = like.obs().models()[srcname].spectral() + self.m_but = [] + self.m_enebut = [] - 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 @@ -260,6 +278,7 @@ class SpectrumPlotter(analysisutils.base,SpectralBase): 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)): @@ -269,74 +288,24 @@ class SpectrumPlotter(analysisutils.base,SpectralBase): 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") - + #convert the flux in sed if asked + self.m_energy = self._convertE(self.m_butterfly['ener']["value"],"MeV","TeV") #gtbutterfly gives MeV + self.m_flux = self._convertSED(self.m_butterfly['dnde']["value"],self.m_butterfly['ener']["value"]) + self.m_error_d = self._convertSED(self.m_butterfly['dnde']["ed_value"],self.m_butterfly['ener']["value"]) + self.m_error_u = self._convertSED(self.m_butterfly['dnde']["eu_value"],self.m_butterfly['ener']["value"]) + self._makebutterfly() # make the butterfly + + self.success("Butterfly ready to be plotted") +# 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]) + """make the butterfly by reading the result and making a closed table""" + for i in xrange(len(self.m_energy)): + self.m_but.append(self.m_error_u[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]) + for i in xrange(len(self.m_energy)): + idx = len(self.m_energy)-i-1 + self.m_but.append(self.m_error_d[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]) -- libgit2 0.21.4