Fabio Iocco / ctools

Browse Code »

Commit 4e5640586d55059a016cb66b8fcf4b77b866e6f6

Authored by Jürgen Knödlseder 2015-06-30 23:07:39 +0200

1 parent 996feeb5

Remove analysisutils and plotutils from distribution.

Inline Side-by-side

Showing 2 changed files with 0 additions and 996 deletions

pyext/ctools/analysisutils.py deleted

View file @996feeb

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.dloglike2) 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 @996feeb

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]2self.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.MeV2ergEnergy_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()