Olga Sergijenko / ctools

Commit d774f1a08d403756bc8e8c118213703babafa5ca

Authored by Jürgen Knödlseder
2 parents 24451dd3 365f2934

Merge branch 'integration' into devel

Inline Side-by-side
Showing 4 changed files with 1 additions and 1001 deletions
pyext/ctools/__init__.py
  View file @d774f1a
... ... @@ -3,9 +3,6 @@ import gammalib
3 3  
4 4 # Import sub modules
5 5 from .tools import *
6   -#import obsutils (require explicit from ctools import obsutils)
7   -#import analysisutils (require explicit from ctools import obsutils)
8   -#import plotutils (require explicit from ctools import obsutils)
9 6  
10 7 # Remove swig register stuff
11 8 bad_entries = [entry for entry in list(locals())
... ...
pyext/ctools/analysisutils.py deleted
View file @24451dd
1   -# ==========================================================================
2   -# This script provides a number of functions that are useful for CTA data
3   -# analysis.
4   -#
5   -# Copyright (C) 2014 David Sanchez
6   -#
7   -# This program is free software: you can redistribute it and/or modify
8   -# it under the terms of the GNU General Public License as published by
9   -# the Free Software Foundation, either version 3 of the License, or
10   -# (at your option) any later version.
11   -#
12   -# This program is distributed in the hope that it will be useful,
13   -# but WITHOUT ANY WARRANTY; without even the implied warranty of
14   -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15   -# GNU General Public License for more details.
16   -#
17   -# You should have received a copy of the GNU General Public License
18   -# along with this program. If not, see <http://www.gnu.org/licenses/>.
19   -#
20   -# ==========================================================================
21   -import ctools as ct
22   -import gammalib as gl
23   -from math import sqrt,log10,floor
24   -import sys
25   -
26   -class base(object):
27   - def __init__(self):
28   - self.classname = self.__class__.__name__
29   - self.errorcolor = "\033[31m"#red
30   - self.infocolor = "\033[34m"#blue
31   - self.warningcolor = "\033[33m"#yellow
32   - self.successcolor = "\033[32m"#green
33   - self.endcolor = "\033[0m"#reset
34   - self.prependfunction = True
35   - self.basemembers = ["classname","errorcolor","infocolor","warningcolor","successcolor","endcolor","prependfunction"]
36   -
37   - def error(self,message, functionname = ""):
38   - printstring = ""
39   - if functionname == "":
40   - printstring = "\n"+self.errorcolor+"*** Error ["+self.classname+"]: "+message+" ***\n"+self.endcolor
41   - else:
42   - printstring = "\n"+self.errorcolor+"*** Error ["+self.classname+"::"+functionname+"]: "+message+" ***\n"+self.endcolor
43   - sys.exit(printstring)
44   -
45   - def info(self,message,newline=True):
46   - printstring = self.infocolor+"["+self.classname+"]: "+message+self.endcolor
47   - if newline:
48   - print(printstring)
49   - else:
50   - print(self.infocolor+message+self.endcolor)
51   - sys.stdout.flush()
52   -
53   - def warning(self,message,functionname = ""):
54   - printstring = ""
55   - if functionname == "":
56   - printstring = self.warningcolor+"["+self.classname+"] Warning: "+message+self.endcolor
57   - else:
58   - printstring = self.warningcolor+"["+self.classname+"::"+functionname+"] Warning: "+message+self.endcolor
59   - print(printstring)
60   -
61   - def success(self,message):
62   - printstring = self.successcolor+"["+self.classname+"]: "+message+self.endcolor
63   - print(printstring)
64   -
65   - def progress(self,message="."):
66   - string = self.infocolor+message+self.endcolor
67   - print(string)
68   -
69   -
70   -class ResultsSorage(dict,base):
71   - """class to store the results in the good format"""
72   - def __init__(self,*arg,**kw):
73   - super(ResultsSorage, self).__init__(*arg, **kw)
74   - base.__init__(self)
75   - self["dnde"] = {}
76   - self["ulim_dnde"] = []
77   - self["time"] = {}
78   - self["ener"] = {}
79   - self["TS"] = []
80   - self["Iflux"] = {}
81   -
82   - self["time"]["tmin_value"] = []
83   - self["time"]["tmax_value"] = []
84   - self["time"]["unit"] = "sec"
85   -
86   - self["ener"]["value"] = []
87   - self["ener"]["ed_value"] = []
88   - self["ener"]["eu_value"] = []
89   - self["ener"]["unit"] = "MeV"
90   -
91   - self["dnde"]["value"] = []
92   - self["dnde"]["ed_value"] = []
93   - self["dnde"]["eu_value"] = []
94   - self["dnde"]["unit"] = "ph/cm2/s/MeV"
95   -
96   - self["Iflux"]["value"] = []
97   - self["Iflux"]["ed_value"] = []
98   - self["Iflux"]["eu_value"] = []
99   - self["Iflux"]["unit"] = "ph/cm2/s"
100   -
101   - def Print(self):
102   - self.success("Results are :")
103   - for k in self.iterkeys():
104   - self.info("Key name : "+k)
105   - try :
106   - for ki in self[k].iterkeys():
107   - print("\t",ki,"\t",self[k][ki])
108   - except:
109   - print(self[k])
110   - print()
111   -
112   -class Analyser(base):
113   - def __init__(self):
114   - super(Analyser,self).__init__()
115   - self.info("Creating "+self.classname+" object")
116   - self.m_obs = None
117   - self.like = None
118   -
119   - self.m_ebinalg = "LOG"
120   - self.m_eunit = "TeV"
121   - self.m_enumbins = 10#nbins
122   - self.m_usepnt = True # Use pointing for map centre
123   - self.m_coordsys = "CEL"#coord
124   - self.m_proj = "TAN"#proj
125   - self.m_binned = False
126   - self.m_edisp = False
127   - self.m_stat = "POISSON"#stat
128   -
129   - self.m_source_set = False
130   - self.m_files_set = False
131   - self.m_time_set = False
132   - self.m_energy_set = False
133   - self.m_irfs_set = False
134   - self.m_roi_set = False
135   - self.m_xml_set = False
136   -
137   - def set_xml(self,xml):
138   - self.m_xml = xml
139   - ind = xml.find(".xml")
140   - if not ind==-1:
141   - self.m_xml_out = xml[:xml.find(".xml")]+"_results"+xml[xml.find(".xml"):]
142   - self.m_xml_set = True
143   -
144   - def set_roi(self,roi):
145   - self.m_roi = float(roi)
146   - self.m_nxpix = int(self.m_roi*100) #nxpix
147   - self.m_nypix = int(self.m_roi*100) #nypix
148   - self.m_binsz = 0.05 #binsz
149   - self.m_roi_set = True
150   -
151   - def set_irfs(self,caldb,irf):
152   - self.m_caldb = caldb
153   - self.m_irf = irf
154   - self.m_irfs_set = True
155   -
156   - def set_Fitsfiles(self,evtfile,tag="CTA_analysis"):
157   - self.m_evtfile = "event_selected_"+tag+".fits"
158   - self.m_raw_evt = evtfile
159   - self.m_cntfile = "countmap_"+tag+".fits"
160   - self.m_files_set = True
161   -
162   - def set_source(self,name,ra,dec):
163   - self.m_ra = ra
164   - self.m_dec = dec
165   - self.m_name = name
166   - self.m_source_set = True
167   -
168   - def validate(self):
169   - if not self.m_time_set :
170   - self.error("Time range of the observations not set")
171   -
172   - if not self.m_files_set :
173   - self.error("FITS files of the observations not set")
174   -
175   - if not self.m_source_set :
176   - self.error("Source not set")
177   -
178   - if not self.m_energy_set :
179   - self.error("Energy range of the observations not set")
180   -
181   - if not self.m_irfs_set :
182   - self.error("IRFs to use not set")
183   -
184   - if not self.m_roi_set :
185   - self.error("ROI to use not set")
186   -
187   - if not self.m_xml_set :
188   - self.error("XML files for the sky model not set")
189   -
190   - def copy(self,analyse,tag="CTA_analysis_copy"):
191   - self.set_xml(analyse.m_xml)
192   - self.set_irfs(analyse.m_caldb,analyse.m_irf)
193   - self.set_roi(analyse.m_roi)
194   - self.set_time_boundary(analyse.m_tmin,analyse.m_tmax)
195   - self.set_energy_boundary(analyse.m_emin,analyse.m_emax)
196   - self.set_source(analyse.m_name,analyse.m_ra,analyse.m_dec)
197   - self.set_Fitsfiles(analyse.m_evtfile,tag)
198   -
199   -
200   - def set_obs(self,obs):
201   - self.m_obs = obs #set the GObservation container
202   -
203   - def set_time_boundary(self,tmin,tmax):
204   - """Set the start and stop time of the observation"""
205   - self.m_tmin = float(tmin)
206   - self.m_tmax = float(tmax)
207   - self.m_time_set = True
208   -
209   - def set_energy_boundary(self,emin,emax):
210   - """Set the energy bounds of the observation"""
211   - self.m_emin = emin
212   - self.m_emax = emax
213   - ebounds = gl.GEbounds(gl.GEnergy(emin, self.m_eunit), \
214   - gl.GEnergy(emax, self.m_eunit))
215   -
216   - if self.m_obs:
217   - for obs in self.m_obs:
218   - obs.events().ebounds(ebounds)
219   - self.m_energy_set = True
220   -
221   - def ctselect(self):
222   - # ctselect application and set parameters
223   - self.info("Running ctselect to cut on events")
224   - self.validate()
225   - if self.m_obs:
226   - filter = ct.ctselect(self.m_obs)
227   - else:
228   - filter = ct.ctselect()
229   -
230   - filter["infile"] = self.m_raw_evt
231   - filter["outfile"] = self.m_evtfile
232   - filter["usepnt"].boolean(self.m_usepnt) # Use pointing for map centre
233   - filter["ra"] = self.m_ra
234   - filter["dec"] = self.m_dec
235   - filter["rad"] = self.m_roi
236   - filter["tmin"] = self.m_tmin
237   - filter["tmax"] = self.m_tmax
238   - filter["emin"].real(self.m_emin)
239   - filter["emax"].real(self.m_emax)
240   - filter["expr"] = ""
241   - filter.logFileOpen()
242   -
243   - filter.run()
244   - if not(self.m_obs):
245   - filter.save()
246   -
247   - if self.m_obs:
248   - # Make a deep copy of the observation that will be returned
249   - # (the ctbin object will go out of scope one the function is
250   - # left)
251   - self.m_obs = filter.obs().copy()
252   -
253   - def ctbin(self,log=False,debug=False):
254   - # ctbin application and set parameters
255   - self.validate()
256   - self.info("Running ctbin to create count map")
257   - if self.m_obs:
258   - bin = ct.ctbin(self.m_obs)
259   - else:
260   - bin = ct.ctbin()
261   - bin["evfile"] = self.m_evtfile
262   - bin["outfile"] = self.m_cntfile
263   - bin["ebinalg"] = self.m_ebinalg
264   - bin["emin"] = self.m_emin
265   - bin["emax"] = self.m_emax
266   - Nbdecade = log10(self.m_emax)-log10(self.m_emin)#Compute the number of decade
267   - bin["enumbins"] = int(self.m_enumbins*Nbdecade)
268   - bin["usepnt"] = self.m_usepnt # Use pointing for map centre
269   - bin["nxpix"] = self.m_nxpix
270   - bin["nypix"] = self.m_nypix
271   - bin["binsz"] = self.m_binsz
272   - bin["coordsys"] = self.m_coordsys
273   - bin["proj"] = self.m_proj
274   -
275   - # Optionally open the log file
276   - if log:
277   - bin.logFileOpen()
278   -
279   - # Optionally switch-on debugging model
280   - if debug:
281   - bin["debug"] = True
282   -
283   - # Run ctbin application. This will loop over all observations in
284   - # the container and bin the events in counts maps
285   - bin.execute()
286   - if self.m_obs:
287   - # Make a deep copy of the observation that will be returned
288   - # (the ctbin object will go out of scope one the function is
289   - # left)
290   - self.m_obs = bin.obs().copy()
291   -
292   - def create_fit(self,log=False,debug=False):
293   - self.validate()
294   - # create ctlike instance with given parameters
295   - self.info("Fitting Data using ctlike")
296   - if self.m_obs:
297   - self.like = ct.ctlike(self.m_obs)
298   - else:
299   - self.like = ct.ctlike()
300   - if self.m_binned:
301   - self.like["infile"] = self.m_cntfile
302   - else:
303   - self.like["infile"] = self.m_evtfile
304   -
305   - self.like["stat"] = self.m_stat
306   - self.like["caldb"] = self.m_caldb
307   - self.like["irf"] = self.m_irf
308   - self.like["srcmdl"] = self.m_xml
309   - self.like["edisp"] = self.m_edisp
310   - self.like["outmdl"] = self.m_xml_out
311   -
312   - # Optionally open the log file
313   - if log:
314   - self.like.logFileOpen()
315   - # Optionally switch-on debugging model
316   - if debug:
317   - self.like["debug"] = True
318   -
319   - def fit(self,log=False,debug=False):
320   - self.validate()
321   - if not(self.like):
322   - self.warning("ctlike object not created, creating now")
323   - self.create_fit(log,debug)
324   -
325   - # Run ctlike application.
326   - self.like.run()
327   - # Save the results in XML
328   - self.like.save()
329   - self.success("Fit performed")
330   -
331   - def PrintResults(self,srcname = ""):
332   - self.info("Results of the Fit")
333   - for m in self.like.obs().models():
334   - if srcname == m.name() or srcname=="":
335   - print("Model : "+m.name())
336   - print(m)
337   -
338   - def GetSrcResuls(self,Res,srcname,E0=None,factor = 1e6,tmin=None,tmax=None):
339   - """function to get and store the results of a fit"""
340   - results = self.GetSrcParams(srcname)
341   -
342   - if 1:
343   - # #Time
344   - # if (tmin is not None) or (tmax is not None) :
345   - # Res["time"]["tmin_value"].append(tmin)
346   - # Res["time"]["tmax_value"].append(tmax)
347   - # Res["time"]["unit"] = "sec"
348   - # else:
349   - Res["time"]["tmin_value"].append(self.m_tmin)
350   - Res["time"]["tmax_value"].append(self.m_tmax)
351   - Res["time"]["unit"] = "sec"
352   -
353   - #Energy
354   - if (E0 is not None):
355   - Res["ener"]["value"].append(E0*factor)
356   - Res["ener"]["ed_value"].append((E0-self.m_emin)*factor)
357   - Res["ener"]["eu_value"].append((self.m_emax-E0)*factor)
358   - Res["ener"]["unit"] = "MeV"
359   - else:
360   - Res["ener"]["value"].append(sqrt(self.m_emin*self.m_emax))
361   - Res["ener"]["ed_value"].append(-self.m_emin+sqrt(self.m_emin*self.m_emax))
362   - Res["ener"]["eu_value"].append(self.m_emax-sqrt(self.m_emin*self.m_emax))
363   - Res["ener"]["unit"] = self.m_eunit
364   -
365   - Res["Iflux"]["value"].append(results["flux"])
366   - Res["Iflux"]["eu_value"].append(results["eflux"])
367   - Res["Iflux"]["ed_value"].append(results["eflux"])
368   - # Res["Iflux"]["unit"] = "unit"
369   -
370   - m = self.GetModel(srcname)
371   - Res["dnde"]["value"].append(m.spectral()[0].value())
372   - Res["dnde"]["eu_value"].append(m.spectral()[0].error())
373   - Res["dnde"]["ed_value"].append(m.spectral()[0].error())
374   -
375   - Res["TS"].append(100)
376   -
377   - UL = UpperLimitComputer(self,srcname)
378   - UL.run()
379   - Res["ulim_dnde"].append(0)
380   - # del UL
381   - return Res
382   -
383   - def GetSrcParams(self,srcname):
384   - results = {}
385   - found = False
386   - m = self.GetModel(srcname)
387   - results["flux"] = m.spectral().flux(gl.GEnergy(self.m_emin,self.m_eunit),gl.GEnergy(self.m_emax,self.m_eunit))
388   - results["eflux"] = m.spectral().eflux(gl.GEnergy(self.m_emin,self.m_eunit),gl.GEnergy(self.m_emax,self.m_eunit))
389   - for par in m.spectral():
390   - if par.is_free():
391   - results[par.name()] = par.value()
392   - results["e"+par.name()] = par.error()
393   -
394   - return results
395   -
396   - def GetModel(self,srcname):
397   - found = False
398   - for m in self.like.obs().models():
399   - if m.name() == srcname:
400   - return m
401   - self.error("source "+srcname+" not in the source list")
402   -
403   -
404   -def MakeEbin(analyse,nbins,srcname,binned = False):
405   - #Store the results
406   - Res = ResultsSorage()
407   -
408   - #check the source name, exit is the srcname is false
409   - analyse.GetModel(srcname)
410   -
411   - #compute the conversion factor to have energy in MeV
412   - Energy_list = {"MeV":1,"GeV":1e3,"TeV":1e6} #TODO
413   - factor = Energy_list[analyse.m_eunit]
414   -
415   - for i in xrange(nbins):
416   - obs_copy = analyse.like.obs().copy()
417   - emin = pow(10,log10(analyse.m_emin) + (log10(analyse.m_emax)-log10(analyse.m_emin))/(nbins)*i)
418   - emax = pow(10,log10(analyse.m_emin) + (log10(analyse.m_emax)-log10(analyse.m_emin))/(nbins)*(i+1))
419   - E0 = pow(10, (log10(emin) + log10(emax)) / 2)
420   - print("Making energy bin between %2.1e"%(emin)+analyse.m_eunit+" and %2.1e"%(emax)+analyse.m_eunit)
421   -
422   - # change the PivotEnergy to E0
423   - m = obs_copy.models()[srcname]
424   - for par in m.spectral():
425   - if par.name() == "PivotEnergy":
426   - par.value(E0*factor)
427   -
428   - ana_bin = Analyser()
429   -
430   - #copy the obs object
431   - ana_bin.set_obs(obs_copy)
432   -
433   - #copy the analyser parameters
434   - ana_bin.copy(analyse)
435   -
436   - # Set energy boundaries
437   - ana_bin.set_energy_boundary(emin,emax)
438   - ana_bin.ctselect()
439   -
440   - # if binned:
441   - # ana_bin.ctbin()
442   -
443   - # Create the fit object
444   - ana_bin.create_fit(log=False,debug=False)
445   -
446   - # Run ctlike application.
447   - ana_bin.like.run()
448   -
449   - #Retrive the results in the good format
450   - Res = ana_bin.GetSrcResuls(Res,srcname,E0,factor)
451   -
452   - del ana_bin
453   - return Res
454   -
455   -def MakeLC(analyse,nbins,srcname,binned = False):
456   - #Store the results
457   - Res = ResultsSorage()
458   -
459   - #check the source name, exit is the srcname is false
460   - analyse.GetModel(srcname)
461   -
462   - for i in xrange(nbins):
463   - obs_copy = analyse.like.obs().copy()
464   - tmin = analyse.m_tmin + (analyse.m_tmax-analyse.m_tmin)/(nbins)*i
465   - tmax = analyse.m_tmin + (analyse.m_tmax-analyse.m_tmin)/(nbins)*(i+1)
466   - print("Making time bin between %2.1e"%(tmin)+" and %2.1e"%(tmax))
467   -
468   - ana_bin = Analyser()
469   -
470   - #copy the obs object
471   - ana_bin.set_obs(obs_copy)
472   -
473   - #copy the analyser parameters
474   - ana_bin.copy(analyse)
475   -
476   - # Set energy boundaries
477   - ana_bin.set_time_boundary(tmin,tmax)
478   - ana_bin.ctselect()
479   -
480   - # if binned:
481   - # ana_bin.ctbin()
482   -
483   - # Create the fit object
484   - ana_bin.create_fit(log=False,debug=False)
485   -
486   - # Run ctlike application.
487   - ana_bin.like.run()
488   -
489   - #Retrive the results in the good format
490   - Res = ana_bin.GetSrcResuls(Res,srcname,tmin=tmin,tmax=tmax)
491   -
492   - del ana_bin
493   - return Res
494   -
495   -class UpperLimitComputer(base):
496   - def __init__(self,analyser,srcname,parname = "Prefactor"):
497   - super(UpperLimitComputer,self).__init__()
498   - self.info("Creating "+self.classname+" object")
499   - # self.like = analyser.like.copy()
500   - obs = analyser.like.obs().copy()
501   - self.like = ct.ctlike(obs)
502   - self.succes = False
503   -
504   - # get spectral model
505   - self.model = self.like.obs().models()[srcname]
506   - self.spec = self.model.spectral()
507   - self.parname = parname
508   - self.bestloglike = analyser.like.opt().value()
509   - self.dloglike = 3.84/2.0 #95% CL now
510   -
511   - def _bisec(self,a,b,tol = 1e-12):
512   - """iterative function to find the roots"""
513   - if (self._func(a)==0.):
514   - return a
515   - elif (self._func(b)==0):
516   - return b
517   - elif (self._func(b)*self._func(a)>0.):
518   - self.warning("Finding root for upper limit calculation failed")
519   - return 0
520   - mid = (a+b)/2
521   - err = abs(a-b)/2
522   - fm = self._func(mid)
523   - if (err<tol or self._func(mid)==0.):
524   - self.succes = True
525   - return mid
526   - elif (fm*self._func(a)<0):
527   - b = mid
528   - else :
529   - a = mid
530   - return self._bisec(a,b)
531   -
532   - def _func(self,value):
533   - """function to be minimized"""
534   - # get parmeter value
535   - val = value*self.spec[self.parname].scale()
536   -
537   - # ensure value to be above minimum boundary
538   - if val <= self.spec[self.parname].min():
539   - return 0.0-self.self.bestloglike-self.dloglike
540   -
541   - # set value
542   - self.spec[self.parname].value(val)
543   -
544   - # Recalculate likelihood
545   - self.like.run()
546   - logl = self.like.opt().value()
547   -
548   - # Return likelihood difference
549   - return (logl-self.bestloglike-self.dloglike)
550   -
551   - def run(self):
552   -
553   - # get spectral parameters and set start value
554   - value = self.spec[self.parname].value()
555   - error = self.spec[self.parname].error()
556   - scale = self.spec[self.parname].scale()
557   -
558   - #Fix all parameters
559   - for model in self.like.obs().models():
560   - for par in model:
561   - par.fix()
562   -
563   -
564   - # set start value to 2 sigma above fit value
565   - self.startpar = value / scale + 2* error/scale
566   -
567   - # set start (a) and (b) depend on the CL asked
568   - a= value / scale + (self.dloglike/4.)*error/scale
569   - b=value / scale + (self.dloglike*2)* error/scale
570   - ul = self._bisec(a,b)
571   -
572   - self.ulimit = ul*self.spec[self.parname].scale()
573   - # Test for convergence
574   - if self.succes:
575   - self.success("Upper limit of parameter \""+self.parname +"\" determined to "+str(self.ulimit))
pyext/ctools/plotutils.py deleted
View file @24451dd
1   -# ==========================================================================
2   -# This script provides a number of functions that are useful for plotting.
3   -#
4   -# Copyright (C) 2014 David Sanchez
5   -#
6   -# This program is free software: you can redistribute it and/or modify
7   -# it under the terms of the GNU General Public License as published by
8   -# the Free Software Foundation, either version 3 of the License, or
9   -# (at your option) any later version.
10   -#
11   -# This program is distributed in the hope that it will be useful,
12   -# but WITHOUT ANY WARRANTY; without even the implied warranty of
13   -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14   -# GNU General Public License for more details.
15   -#
16   -# You should have received a copy of the GNU General Public License
17   -# along with this program. If not, see <http://www.gnu.org/licenses/>.
18   -#
19   -# ==========================================================================
20   -import ctools as ct
21   -import gammalib as gl
22   -from math import log10,pow,sqrt
23   -import os
24   -
25   -try:
26   - import matplotlib.pyplot as plt
27   - import matplotlib.gridspec as gridspec
28   - has_matplotlib = True
29   -except:
30   - has_matplotlib = False
31   -
32   -try:
33   - import aplpy
34   - has_aplpy = True
35   -except:
36   - has_aplpy = False
37   -
38   -class options(object):
39   - def __init__(self):
40   - self.fmt = 'o'
41   - self.color = BLUE
42   - self.linecolor = 'black'
43   - self.markercolor = self.color
44   - self.tslim = 4.0
45   - self.elinewidth = 2.0
46   - self.linewidth = 2.0
47   - self.linestyle = '--'
48   - self.capsize = 5.0
49   - self.markersize = 5.0
50   - self.limlength = 0.4
51   - self.label = "_nolegend_"
52   -
53   -import analysisutils
54   -
55   -class SpectralBase(object):
56   - def __init__(self,sed=True,sed_factor=1.):
57   - self.m_sed_factor = sed_factor
58   - self.m_sed = sed
59   -
60   - def _convertSED(self,tab,energy):
61   - """convert the a table storing dN/dE in sed with right units"""
62   -
63   - if len(tab) != len(energy):
64   - self.error("input table and energy table does not match")
65   - return
66   - if self.m_sed:
67   - for i in xrange(len(tab)):
68   - tab[i] *= energy[i]**2*self.m_sed_factor
69   - return tab
70   -
71   -BLUE="#337EF5"
72   -GREY='grey'
73   -GREEN = '#21851D'
74   -RED = '#BD0B0B'
75   -YELLOW = '#D4CF44'
76   -ORANGE = '#D19534'
77   -
78   -
79   -class residuals(analysisutils.base,options,SpectralBase):
80   - def __init__(self,spectrum,datapoint,sed,sed_factor,eunit="TeV"):
81   - super(residuals,self).__init__()
82   - options.__init__(self)
83   - SpectralBase.__init__(self,sed,sed_factor)
84   -
85   - # load the gammalib spectrum object
86   - self.spectrum = spectrum
87   - self.SetDataPoint(datapoint)
88   -
89   - self.xerrors = True
90   - self.refUnit = eunit # reference unit for conversion
91   -
92   - def SetDataPoint(self,datapoint):
93   - self.m_point_info = datapoint
94   - self.Npt = len(self.m_point_info["dnde"]["value"])
95   -
96   - def _plt_residuals(self):
97   - th_val = []
98   - index = []
99   - ed_y = []
100   - eu_y = []
101   - y = []
102   -
103   - Energy_list = {"MeV":1,"GeV":1e-3,"TeV":1e-6} #TODO
104   - for i in range(self.Npt):
105   - if self.m_point_info["TS"][i] >= self.tslim:
106   - th_val.append(self.spectrum.eval(gl.GEnergy(self.m_point_info["ener"]["value"][i],self.m_point_info["ener"]["unit"]),gl.GTime(0)))
107   - index.append(i)
108   - else:
109   - th_val.append(0)
110   - if self.m_sed:
111   - self._convertSED(th_val,self.m_point_info["ener"])
112   - for i in index:
113   - y.append((self.m_point_info["dnde"]["value"][i] - th_val[i]) / th_val[i])
114   - ed_y.append(self.m_point_info["dnde"]["ed_value"][i]/th_val[i])
115   - eu_y.append(self.m_point_info["dnde"]["eu_value"][i]/th_val[i])
116   -
117   - if self.xerrors:
118   - 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)
119   - else:
120   - plt.errorbar(self.m_point_info["ener"]["value"], y, yerr=[ed_y,eu_y],fmt=self.fmt,color=self.markercolor,elinewidth=self.elinewidth)
121   -
122   - plt.axhline(0.0,color=self.linecolor,lw=self.linewidth,ls=self.linestyle)
123   -
124   -
125   -class ulimgraph(analysisutils.base,options,SpectralBase):
126   - def __init__(self,datapoint,sed,sed_factor,eunit='TeV'):
127   - super(ulimgraph,self).__init__()
128   - #~ super(ulimgraph,self).__init__()
129   - options.__init__(self)
130   - SpectralBase.__init__(self,sed,sed_factor)
131   -
132   - self.xerrors=True
133   - self.nlims = -1
134   - self.refUnit = eunit
135   - self.SetDataPoint(datapoint)
136   -
137   - def SetDataPoint(self,datapoint):
138   - self.m_point_info = datapoint
139   - self.Npt = len(self.m_point_info["dnde"]["value"])
140   -
141   - def _plt_points(self):
142   - results = analysisutils.ResultsSorage()
143   -
144   - Energy_list = {"MeV":1,"GeV":1e-3,"TeV":1e-6} #TODO
145   - for i in range(self.Npt):
146   -
147   - self.m_point_info["ener"]["value"][i]*=Energy_list[self.refUnit]
148   - self.m_point_info["ener"]["ed_value"][i]*=Energy_list[self.refUnit]
149   - self.m_point_info["ener"]["eu_value"][i]*=Energy_list[self.refUnit]
150   - self.m_point_info["ener"]["unit"] = self.refUnit
151   -
152   - if self.m_point_info["TS"][i] >= self.tslim:
153   -
154   - if self.m_point_info["dnde"]["ed_value"][i]>self.m_point_info["dnde"]["value"][i]:
155   - 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]))
156   - self.warning("\tTS value is however TS="+str(self.m_point_info["TS"][i])+", weird!")
157   - else :
158   - self.m_point_info["dnde"]["ed_value"][i] = 0
159   - self.m_point_info["dnde"]["eu_value"][i] = 0
160   -
161   - if self.xerrors:
162   - 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)
163   - else:
164   - 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)
165   -
166   - def _plt_limits(self):
167   - ed_flux = []
168   - eu_flux = []
169   - flux = []
170   - ener = []
171   -
172   - drawn_lims = 0
173   - Energy_list = {"MeV":1,"GeV":1e-3,"TeV":1e-6} #TODO
174   - for i in range(self.Npt):
175   - if self.m_point_info["TS"][i] < self.tslim:
176   - self.info("Plot an upper limit for E = "+self.m_point_info["ener"]["value"][i]+" with a TS = "+self.m_point_info["TS"][i])
177   - if drawn_lims >= self.nlims and self.nlims!=-1:
178   - break
179   -
180   - ed_flux.append(self.limlength*self.m_point_info["ulim_dnde"][i])
181   - eu_flux.append(0.0)
182   - flux.append(self.m_point_info["ulim_dnde"][i])
183   - drawn_lims+=1
184   -
185   - # if self.xerrors:
186   - # 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)
187   - # else:
188   - 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)
189   -
190   -
191   -class LightCurvePlotter(analysisutils.base,options):
192   - """ Class to plot the Light curve"""
193   - def __init__(self,srcname,datapoint):
194   - super(LightCurvePlotter,self).__init__()
195   - options.__init__(self)
196   - self.m_name = srcname
197   - self.m_data = datapoint #data point in a Results Storage class
198   -
199   - def draw(self):
200   - if not(has_matplotlib):
201   - self.warning("matplotlib module not found, can draw")
202   - return
203   -
204   - time = []
205   - dtime = []
206   - flux = []
207   - dflux_ed = []
208   - dflux_eu = []
209   - for i in range(len(self.m_data["time"]["tmin_value"])):
210   - time.append((self.m_data["time"]["tmin_value"][i]+self.m_data["time"]["tmax_value"][i])/2)
211   - dtime.append((self.m_data["time"]["tmax_value"][i]-self.m_data["time"]["tmin_value"][i])/2)
212   - flux.append(self.m_data["Iflux"]["value"][i])
213   - dflux_eu.append(self.m_data["Iflux"]["eu_value"][i])
214   - dflux_ed.append(self.m_data["Iflux"]["ed_value"][i])
215   -
216   - plt.figure("results",figsize=(12,7),edgecolor="w")
217   - unit = self.m_data["time"]["unit"]
218   - plt.xlabel("Time ["+ unit+"]",fontsize=15)
219   - unit = self.m_data["Iflux"]["unit"]
220   - plt.ylabel("E$^{2}$ dN/dE ["+unit+"]",fontsize=15)
221   -
222   - print(time)
223   - print(dtime)
224   - print(flux)
225   - print(dflux_ed)
226   - print(dflux_eu)
227   - plt.errorbar(time,flux,xerr=[dtime,dtime],yerr=[dflux_ed,dflux_eu],fmt=self.fmt,color=self.markercolor,elinewidth=self.elinewidth)
228   - plt.show()
229   -
230   -class SpectrumPlotter(analysisutils.base,SpectralBase):
231   - """ Class to compute and plot the spectra (i.e best fit model and butterfly).
232   - if data points are provided, only change applied is energy unit convertion to match the butterfly one
233   - (i.e. no dN/dE to SED conversion).
234   - It is currently assumed that the unit is MeV for the energy of the data point. this will be change
235   - once the I/O container format will be defined."""
236   - def __init__(self,srcname,like,datapoint,emin=0.1,emax=100,energy = "TeV",SED = True ,npt = 50):
237   - super(SpectrumPlotter,self).__init__()
238   -
239   - self.m_ebound = [emin,emax]
240   - self.m_name = srcname
241   - self.m_like = like
242   - self.m_spectral = like.obs().models()[srcname].spectral()
243   - self.m_energy = []
244   - self.m_flux = []
245   - self.m_error = []
246   - self.m_but = []
247   - self.m_enebut = []
248   - self.m_npt = npt
249   - self.m_covar = gl.GMatrix(1,1)
250   - self.covar()
251   - #~ this part take care of the style : SED or differential flux
252   - #~ in the case of SED, the y-unit is erg/cm2/s
253   - self.m_sed = SED
254   - self.eunit = energy
255   - self.m_sed_factor = 1.
256   - self._validateUnits()
257   - SpectralBase.__init__(self,self.m_sed,self.m_sed_factor)
258   -
259   - #Assume MeV for the energy #TODO
260   - self.points = ulimgraph(datapoint,self.m_sed,gl.MeV2erg,energy)
261   - self.residuals = residuals(self.m_spectral,datapoint,self.m_sed,gl.MeV2erg,energy)
262   -
263   - def _validateUnits(self):
264   - Energy_list = {"MeV":1,"GeV":1e3,"TeV":1e6}
265   - if not(Energy_list.has_key(self.eunit)):
266   - self.warning("Change energy unit to TeV")
267   - self.eunit = "TeV"
268   -
269   - if self.m_sed:
270   - self.m_sed_factor = gl.MeV2erg*Energy_list[self.eunit]**2
271   -
272   - def covar(self):
273   - # Get covariance matrix if fit has converged
274   - try:
275   - fullcovar = self.m_like.obs().function().curvature().invert()
276   - except:
277   - self.warning("Covariance matrix not determined successfully")
278   - return
279   - #~ get the matrix element of the source of interest
280   - #~ 1) get the right element indices
281   - idx = 0
282   - par_index_map = []
283   - for m in self.m_like.obs().models() :
284   - if m.name()!= self.m_name:
285   - idx += m.size()
286   - continue
287   - idx += m.spatial().size()
288   - for par in m.spectral():
289   - if par.is_free():
290   - par_index_map.append(idx)
291   - idx += 1
292   -
293   - #~ 2) get the elemets and store them in a matrix
294   - self.m_covar = gl.GMatrix(len(par_index_map),len(par_index_map))
295   - i = 0
296   - for xpar in par_index_map:
297   - j = 0
298   - for ypar in par_index_map:
299   - self.m_covar[i,j] = fullcovar[xpar,ypar]
300   - j += 1
301   - i += 1
302   - self.success("Covariance Matrix succesfuly computed")
303   -
304   - def _makespectrum(self):
305   - for i in xrange(self.m_npt):
306   - #Compute the energy
307   - lene = log10(self.m_ebound[0])+log10(self.m_ebound[1])*i/(self.m_npt-1)
308   - self.m_energy.append(pow(10,lene))
309   - #Compute the flux
310   - self.m_flux.append(self.m_spectral.eval(gl.GEnergy(pow(10,lene),self.eunit),gl.GTime(0)))
311   -
312   - #Compute the gradients
313   - if self.m_covar.size() != 1:
314   - self.m_spectral.eval_gradients(gl.GEnergy(pow(10,lene), self.eunit),gl.GTime(0))
315   - #store them into a GVector
316   - Derivative = gl.GVector(self.m_covar.size()/2)
317   - j = 0
318   - for par in self.m_spectral:
319   - if par.is_free():
320   - Derivative[j] = par.factor_gradient()
321   - j += 1
322   - #~ computed the error
323   - self.m_error.append(sqrt(Derivative * (self.m_covar*Derivative)))
324   -
325   - #~ convert the flux in sed if asked
326   - self.m_flux = self._convertSED(self.m_flux,self.m_energy)
327   - if self.m_covar.size() != 1:
328   - self.m_error = self._convertSED(self.m_error,self.m_energy)
329   - self._makebutterfly() # make the butterfly
330   - self.success("Spectrum computed")
331   -
332   - def _makebutterfly(self):
333   - """make the butterfly by appending element in a table"""
334   - for i in xrange(self.m_npt):
335   - self.m_but.append(self.m_flux[i]+self.m_error[i])
336   - self.m_enebut.append(self.m_energy[i])
337   - for i in xrange(self.m_npt):
338   - idx = self.m_npt-i-1
339   - self.m_but.append(self.m_flux[idx]-self.m_error[idx])
340   - self.m_enebut.append(self.m_energy[idx])
341   - self.m_but.append(self.m_but[0])
342   - self.m_enebut.append(self.m_enebut[0])
343   -
344   - def write(self):
345   - self.warning("write function not yet implemented")
346   -
347   - def draw(self):
348   - """plot the results, provide a file named ptfile for the flux points"""
349   - gs = gridspec.GridSpec(2, 1,height_ratios=[3,1])
350   - self._makespectrum()
351   - self.write()
352   -
353   - if not(has_matplotlib):
354   - self.warning("matplotlib module not found, can draw")
355   - return
356   -
357   - plt.figure("results",figsize=(12,7),edgecolor="w")
358   -
359   - sb1 = plt.subplot(gs[0])
360   - plt.loglog()
361   - plt.xlabel("Energy ["+ self.eunit+"]",fontsize=15)
362   -
363   - if self.m_sed:
364   - plt.ylabel("E$^{2}$ dN/dE [erg cm$^{2}$ s$^{-1}$]",fontsize=15)
365   - else:
366   - plt.ylabel("dN/dE [cm$^{-2}$ s$^{-1}$ "+ self.eunit+"$^{-1}$]",fontsize=15)
367   -
368   - plt.plot(self.m_energy,self.m_flux)
369   - plt.plot(self.m_enebut,self.m_but)
370   -
371   - #~ plot the data points on top of the buterfly
372   - self.points._plt_points()
373   - self.points._plt_limits()
374   -
375   - #~ plot the residual points
376   - sb2 = plt.subplot(gs[1])
377   - plt.xscale('log')
378   - self.residuals._plt_residuals()
379   -
380   - plt.show()
381   -
382   -
383   -# ============= #
384   -# Plot skymaps #
385   -# ============= #
386   -class MapPlotter(analysisutils.base):
387   - def __init__(self,modelmap, countmap):
388   - super(MapPlotter,self).__init__()
389   - self.modelmap = modelmap
390   - self.countmap = countmap
391   -
392   - if not(has_aplpy) :
393   - self.error("Module aplypy needed")
394   - if not(has_matplotlib) :
395   - self.error("Module matplotlib needed")
396   -
397   - def draw(self):
398   - """ Plot the model map. """
399   - # Load model map
400   - fig = plt.figure()
401   - f1 = aplpy.FITSFigure(self.countmap, figure=fig, subplot=[0.1,0.1,0.35,0.8])
402   - f1.set_tick_labels_font(size='small')
403   - f1.set_axis_labels_font(size='small')
404   - f1.show_colorscale()
405   -
406   - f1.tick_labels.set_yformat('dd:mm:ss')
407   - f1.tick_labels.set_xformat('hh:mm:ss')
408   - f1.axis_labels.set_xtext('Right Ascension (J2000)')
409   - f1.axis_labels.set_ytext('Declination (J2000)')
410   - f1.ticks.set_length(10.5, minor_factor=0.5)
411   -
412   - f2 = aplpy.FITSFigure(self.modelmap, figure=fig, subplot=[0.5,0.1,0.35,0.8])
413   - f2.set_tick_labels_font(size='small')
414   - f2.set_axis_labels_font(size='small')
415   - f2.show_colorscale()
416   -
417   - f2.hide_yaxis_label()
418   - f2.hide_ytick_labels()
419   -
420   - fig.canvas.draw()
421   - fig.show()
pyext/setup.py.in
  View file @d774f1a
... ... @@ -121,8 +121,7 @@ setup(name=&#39;@PACKAGE_NAME@&#39;,
121 121 author='Juergen Knoedlseder',
122 122 author_email='@PACKAGE_BUGREPORT@',
123 123 url='http://gammalib.sourceforge.net/',
124   - py_modules=['ctools.tools', 'ctools.obsutils', 'ctools.analysisutils',
125   - 'ctools.plotutils'],
  124 + py_modules=['ctools.tools', 'ctools.obsutils'],
126 125 ext_modules=[Extension('ctools._tools', ['ctools/tools_wrap.cpp'],
127 126 include_dirs=include_dirs,
128 127 library_dirs=library_dirs,
... ...