Francois Brun / ctools

Commit d11303337da842bf84a741ca82bff8d1645690d6

Authored by sanchez
1 parent 70701752

use the new ctbutterfly + some doc (need more)

Inline Side-by-side
Showing 2 changed files with 187 additions and 141 deletions
pyext/ctools/analysisutils.py
  View file @d113033
... ... @@ -21,9 +21,10 @@
21 21 import ctools as ct
22 22 import gammalib as gl
23 23 from math import sqrt,log10,floor
24   -import sys
  24 +import sys,string
25 25  
26 26 class base(object):
  27 + '''This class is used to defined the colors of the print'''
27 28 def __init__(self):
28 29 self.classname = self.__class__.__name__
29 30 self.errorcolor = "\033[31m"#red
... ... @@ -35,6 +36,7 @@ class base(object):
35 36 self.basemembers = ["classname","errorcolor","infocolor","warningcolor","successcolor","endcolor","prependfunction"]
36 37  
37 38 def error(self,message, functionname = ""):
  39 + '''Print the error message'''
38 40 printstring = ""
39 41 if functionname == "":
40 42 printstring = "\n"+self.errorcolor+"*** Error ["+self.classname+"]: "+message+" ***\n"+self.endcolor
... ... @@ -43,6 +45,7 @@ class base(object):
43 45 sys.exit(printstring)
44 46  
45 47 def info(self,message,newline=True):
  48 + '''Print the information message'''
46 49 printstring = self.infocolor+"["+self.classname+"]: "+message+self.endcolor
47 50 if newline:
48 51 print(printstring)
... ... @@ -50,7 +53,8 @@ class base(object):
50 53 print(self.infocolor+message+self.endcolor)
51 54 sys.stdout.flush()
52 55  
53   - def warning(self,message,functionname = ""):
  56 + def warning(self,message,functionname = ""):*
  57 + '''Print the warning message'''
54 58 printstring = ""
55 59 if functionname == "":
56 60 printstring = self.warningcolor+"["+self.classname+"] Warning: "+message+self.endcolor
... ... @@ -59,19 +63,21 @@ class base(object):
59 63 print(printstring)
60 64  
61 65 def success(self,message):
  66 + '''Print the success message'''
62 67 printstring = self.successcolor+"["+self.classname+"]: "+message+self.endcolor
63 68 print(printstring)
64 69  
65 70 def progress(self,message="."):
  71 + '''Print the progess message'''
66 72 string = self.infocolor+message+self.endcolor
67 73 print(string)
68 74  
69   -
70   -class ResultsSorage(dict,base):
  75 +class ResultsStorage(dict,base):
71 76 """class to store the results in the good format"""
72 77 def __init__(self,*arg,**kw):
73   - super(ResultsSorage, self).__init__(*arg, **kw)
  78 + super(ResultsStorage, self).__init__(*arg, **kw)
74 79 base.__init__(self)
  80 + # Can store the differential flux, upper limit, time, energy and so one
75 81 self["dnde"] = {}
76 82 self["ulim_dnde"] = []
77 83 self["time"] = {}
... ... @@ -79,26 +85,31 @@ class ResultsSorage(dict,base):
79 85 self["TS"] = []
80 86 self["Iflux"] = {}
81 87  
  88 + #store the time : Min, Max and units of the variable
82 89 self["time"]["tmin_value"] = []
83 90 self["time"]["tmax_value"] = []
84 91 self["time"]["unit"] = "sec"
85 92  
  93 + #store the energy : Value min,max (ed_* and eu_*) and units of the variable
86 94 self["ener"]["value"] = []
87 95 self["ener"]["ed_value"] = []
88 96 self["ener"]["eu_value"] = []
89 97 self["ener"]["unit"] = "MeV"
90 98  
  99 + #store the differential flux : Value min,max (ed_* and eu_*) and units of the variable
91 100 self["dnde"]["value"] = []
92 101 self["dnde"]["ed_value"] = []
93 102 self["dnde"]["eu_value"] = []
94 103 self["dnde"]["unit"] = "ph/cm2/s/MeV"
95 104  
  105 + #store the integral flux : Value min,max (ed_* and eu_*) and units of the variable
96 106 self["Iflux"]["value"] = []
97 107 self["Iflux"]["ed_value"] = []
98 108 self["Iflux"]["eu_value"] = []
99 109 self["Iflux"]["unit"] = "ph/cm2/s"
100 110  
101 111 def Print(self):
  112 + """ Print the value stored in the class"""
102 113 self.success("Results are :")
103 114 for k in self.iterkeys():
104 115 self.info("Key name : "+k)
... ... @@ -110,12 +121,17 @@ class ResultsSorage(dict,base):
110 121 print()
111 122  
112 123 class Analyser(base):
  124 + '''Class to run the analysis of the data. Currently the parameters are hard coded but this will change once
  125 + the method to feed the class will be defined
  126 + The classe can work with GObservation as input (simulation) or fits files (simulation+real data)'''
113 127 def __init__(self):
114 128 super(Analyser,self).__init__()
115 129 self.info("Creating "+self.classname+" object")
  130 + #Observation and ctlike object
116 131 self.m_obs = None
117 132 self.like = None
118 133  
  134 + #parameters
119 135 self.m_ebinalg = "LOG"
120 136 self.m_eunit = "TeV"
121 137 self.m_enumbins = 10#nbins
... ... @@ -126,6 +142,8 @@ class Analyser(base):
126 142 self.m_edisp = False
127 143 self.m_stat = "POISSON"#stat
128 144  
  145 + # booleans to ensure that input files and parameters have been defined
  146 + # set to false at the begining and turn to true while a file/parameter is defined.
129 147 self.m_source_set = False
130 148 self.m_files_set = False
131 149 self.m_time_set = False
... ... @@ -135,6 +153,7 @@ class Analyser(base):
135 153 self.m_xml_set = False
136 154  
137 155 def set_xml(self,xml):
  156 + '''Set the model XML file'''
138 157 self.m_xml = xml
139 158 ind = xml.find(".xml")
140 159 if not ind==-1:
... ... @@ -142,6 +161,7 @@ class Analyser(base):
142 161 self.m_xml_set = True
143 162  
144 163 def set_roi(self,roi):
  164 + '''Set the ROI parameter and compute the number of pixel, bin size of the map'''
145 165 self.m_roi = float(roi)
146 166 self.m_nxpix = int(self.m_roi*100) #nxpix
147 167 self.m_nypix = int(self.m_roi*100) #nypix
... ... @@ -149,23 +169,45 @@ class Analyser(base):
149 169 self.m_roi_set = True
150 170  
151 171 def set_irfs(self,caldb,irf):
  172 + '''Set the IRFs definition'''
152 173 self.m_caldb = caldb
153 174 self.m_irf = irf
154 175 self.m_irfs_set = True
155 176  
156 177 def set_Fitsfiles(self,evtfile,tag="CTA_analysis"):
  178 + '''Set the input and output fits files'''
157 179 self.m_evtfile = "event_selected_"+tag+".fits"
158 180 self.m_raw_evt = evtfile
159 181 self.m_cntfile = "countmap_"+tag+".fits"
160 182 self.m_files_set = True
161 183  
162 184 def set_source(self,name,ra,dec):
  185 + '''Set the source name and the RA/DEC position of the source'''
163 186 self.m_ra = ra
164 187 self.m_dec = dec
165 188 self.m_name = name
166 189 self.m_source_set = True
167 190  
  191 + def set_time_boundary(self,tmin,tmax):
  192 + """Set the start and stop time of the observation"""
  193 + self.m_tmin = float(tmin)
  194 + self.m_tmax = float(tmax)
  195 + self.m_time_set = True
  196 +
  197 + def set_energy_boundary(self,emin,emax):
  198 + """Set the energy bounds of the observation"""
  199 + self.m_emin = emin
  200 + self.m_emax = emax
  201 + ebounds = gl.GEbounds(gl.GEnergy(emin, self.m_eunit), \
  202 + gl.GEnergy(emax, self.m_eunit))
  203 +
  204 + if self.m_obs:
  205 + for obs in self.m_obs:
  206 + obs.events().ebounds(ebounds)
  207 + self.m_energy_set = True
  208 +
168 209 def validate(self):
  210 + '''Check is all the parameters has been set'''
169 211 if not self.m_time_set :
170 212 self.error("Time range of the observations not set")
171 213  
... ... @@ -188,6 +230,7 @@ class Analyser(base):
188 230 self.error("XML files for the sky model not set")
189 231  
190 232 def copy(self,analyse,tag="CTA_analysis_copy"):
  233 + '''Make a copy of the class'''
191 234 self.set_xml(analyse.m_xml)
192 235 self.set_irfs(analyse.m_caldb,analyse.m_irf)
193 236 self.set_roi(analyse.m_roi)
... ... @@ -198,35 +241,20 @@ class Analyser(base):
198 241  
199 242  
200 243 def set_obs(self,obs):
201   - self.m_obs = obs #set the GObservation container
  244 + '''set the GObservation container'''
  245 + self.m_obs = obs
202 246  
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 247 def ctselect(self):
222   - # ctselect application and set parameters
  248 + '''ctselect application and set parameters'''
223 249 self.info("Running ctselect to cut on events")
224   - self.validate()
  250 + self.validate()# check the parameters
  251 +
225 252 if self.m_obs:
226 253 filter = ct.ctselect(self.m_obs)
227 254 else:
228 255 filter = ct.ctselect()
229 256  
  257 + #set the parameters for running the tool
230 258 filter["infile"] = self.m_raw_evt
231 259 filter["outfile"] = self.m_evtfile
232 260 filter["usepnt"].boolean(self.m_usepnt) # Use pointing for map centre
... ... @@ -241,6 +269,7 @@ class Analyser(base):
241 269 filter.logFileOpen()
242 270  
243 271 filter.run()
  272 + # if no Gobservation have been defined, then save the fits file.
244 273 if not(self.m_obs):
245 274 filter.save()
246 275  
... ... @@ -251,14 +280,17 @@ class Analyser(base):
251 280 self.m_obs = filter.obs().copy()
252 281  
253 282 def ctbin(self,log=False,debug=False):
254   - # ctbin application and set parameters
255   - self.validate()
  283 + '''ctbin application and set parameters'''
256 284 self.info("Running ctbin to create count map")
  285 + self.validate()# check the parameters
  286 +
257 287 if self.m_obs:
258 288 bin = ct.ctbin(self.m_obs)
259 289 else:
260 290 bin = ct.ctbin()
261 291 bin["evfile"] = self.m_evtfile
  292 +
  293 + #set the parameters for running the tool
262 294 bin["outfile"] = self.m_cntfile
263 295 bin["ebinalg"].string(self.m_ebinalg)
264 296 bin["emin"].real(self.m_emin)
... ... @@ -290,14 +322,17 @@ class Analyser(base):
290 322 self.m_obs = bin.obs().copy()
291 323  
292 324 def create_fit(self,log=False,debug=False):
293   - self.validate()
  325 + '''Create the like object from ctlike'''
294 326 # create ctlike instance with given parameters
295 327 self.info("Fitting Data using ctlike")
  328 + self.validate()
  329 +
  330 + #set the parameters of the tools
296 331 if self.m_obs:
297 332 self.like = ct.ctlike(self.m_obs)
298 333 else:
299 334 self.like = ct.ctlike()
300   - if self.m_binned:
  335 + if self.m_binned#binned or unbinned analysis
301 336 self.like["infile"] = self.m_cntfile
302 337 else:
303 338 self.like["infile"] = self.m_evtfile
... ... @@ -312,10 +347,50 @@ class Analyser(base):
312 347 # Optionally open the log file
313 348 if log:
314 349 self.like.logFileOpen()
  350 +
315 351 # Optionally switch-on debugging model
316 352 if debug:
317 353 self.like["debug"].boolean(True)
318 354  
  355 + def compute_butterfly(self,name,log=False,debug=False):
  356 + '''create ctlike instance with given parameters'''
  357 + self.info("Computing the data with ctbutterfly")
  358 + self.validate()
  359 +
  360 + if self.m_obs:
  361 + butterfly = ct.ctbutterfly(self.m_obs)
  362 + else:#parameters of the tools
  363 + butterfly = ct.ctbutterfly()
  364 + butterfly["infile"] = self.m_evtfile
  365 + butterfly['srcname'] = name
  366 + butterfly["emin"].real(self.m_emin)
  367 + butterfly["emax"].real(self.m_emax)
  368 + butterfly["caldb"].string(self.m_caldb)
  369 + butterfly["irf"].string(self.m_irf)
  370 + butterfly["srcmdl"] = self.m_xml
  371 + filename = name+"_butterfly.txt"
  372 + butterfly["outfile"] = filename
  373 +
  374 + # Optionally open the log file
  375 + if log:
  376 + butterfly.logFileOpen()
  377 + # Optionally switch-on debugging model
  378 + if debug:
  379 + butterfly["debug"].boolean(True)
  380 +
  381 + butterfly.run() #Compute the butterfly
  382 + butterfly.save() # save the txt files
  383 +
  384 + #Store the results in the Storage class
  385 + self.m_butterfly = ResultsStorage()
  386 + data = open(filename,'r').readlines()
  387 + for line in data:
  388 + words = string.split(line)
  389 + self.m_butterfly["ener"]["value"].append(float(words[0]))
  390 + self.m_butterfly["dnde"]["value"].append(float(words[1]))
  391 + self.m_butterfly["dnde"]["ed_value"].append(float(words[1]) - float(words[2]))
  392 + self.m_butterfly["dnde"]["eu_value"].append(float(words[1]) + float(words[2]))
  393 +
319 394 def fit(self,log=False,debug=False):
320 395 self.validate()
321 396 if not(self.like):
... ... @@ -329,81 +404,82 @@ class Analyser(base):
329 404 self.success("Fit performed")
330 405  
331 406 def PrintResults(self,srcname = ""):
  407 + ''' Print results of the fit '''
332 408 self.info("Results of the Fit")
333 409 for m in self.like.obs().models():
334 410 if srcname == m.name() or srcname=="":
335 411 print("Model : "+m.name())
336 412 print(m)
337 413  
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"""
  414 + def GetSrcResuls(self,Res,srcname,E0=None,factor = 1e6):
  415 + """function to get and store the results of a fit. Store the results in the Storage class element"""
340 416 results = self.GetSrcParams(srcname)
341 417  
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):
  418 + #Store the time
  419 + Res["time"]["tmin_value"].append(self.m_tmin)
  420 + Res["time"]["tmax_value"].append(self.m_tmax)
  421 + Res["time"]["unit"] = "sec"
  422 +
  423 + #Store the energy
  424 + if (E0 is not None):#if a 'pivot' energy is given
355 425 Res["ener"]["value"].append(E0*factor)
356 426 Res["ener"]["ed_value"].append((E0-self.m_emin)*factor)
357 427 Res["ener"]["eu_value"].append((self.m_emax-E0)*factor)
358 428 Res["ener"]["unit"] = "MeV"
359   - else:
  429 + else:#if not, the geometric mean energy is used
360 430 Res["ener"]["value"].append(sqrt(self.m_emin*self.m_emax))
361 431 Res["ener"]["ed_value"].append(-self.m_emin+sqrt(self.m_emin*self.m_emax))
362 432 Res["ener"]["eu_value"].append(self.m_emax-sqrt(self.m_emin*self.m_emax))
363 433 Res["ener"]["unit"] = self.m_eunit
364 434  
  435 + #save the flux values
365 436 Res["Iflux"]["value"].append(results["flux"])
366 437 Res["Iflux"]["eu_value"].append(results["eflux"])
367 438 Res["Iflux"]["ed_value"].append(results["eflux"])
368   - # Res["Iflux"]["unit"] = "unit"
369 439  
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())
  440 + #save the differential flux in a src model independent maner
  441 + srcmodel = self.GetModel(srcname)
  442 + Res["dnde"]["value"].append(srcmodel.spectral()[0].value())
  443 + Res["dnde"]["eu_value"].append(srcmodel.spectral()[0].error())
  444 + Res["dnde"]["ed_value"].append(srcmodel.spectral()[0].error())
374 445  
375   - Res["TS"].append(100)
  446 + Res["TS"].append(100) ## TODO use the real TS value
376 447  
377   - UL = UpperLimitComputer(self,srcname)
  448 + UL = UpperLimitComputer(self,srcname) #Compute the UL. To be updated
378 449 UL.run()
379 450 Res["ulim_dnde"].append(0)
380   - # del UL
  451 +
381 452 return Res
382 453  
383 454 def GetSrcParams(self,srcname):
  455 + '''Retrive the source parameters from the like object and store it in a dictonnary'''
384 456 results = {}
385 457 found = False
386 458 m = self.GetModel(srcname)
  459 +
  460 + #get the integral flux
387 461 results["flux"] = m.spectral().flux(gl.GEnergy(self.m_emin,self.m_eunit),gl.GEnergy(self.m_emax,self.m_eunit))
388 462 results["eflux"] = m.spectral().eflux(gl.GEnergy(self.m_emin,self.m_eunit),gl.GEnergy(self.m_emax,self.m_eunit))
  463 +
  464 + # Get free parameters
389 465 for par in m.spectral():
390 466 if par.is_free():
391 467 results[par.name()] = par.value()
392 468 results["e"+par.name()] = par.error()
393   -
  469 +
394 470 return results
395   -
  471 +
396 472 def GetModel(self,srcname):
  473 + ''' Get the spectral parameters of the sources stored in the like object'''
397 474 found = False
398 475 for m in self.like.obs().models():
399 476 if m.name() == srcname:
400 477 return m
401 478 self.error("source "+srcname+" not in the source list")
402   -
403 479  
404 480 def MakeEbin(analyse,nbins,srcname,binned = False):
405 481 #Store the results
406   - Res = ResultsSorage()
  482 + Res = ResultsStorage()
407 483  
408 484 #check the source name, exit is the srcname is false
409 485 analyse.GetModel(srcname)
... ... @@ -454,7 +530,7 @@ def MakeEbin(analyse,nbins,srcname,binned = False):
454 530  
455 531 def MakeLC(analyse,nbins,srcname,binned = False):
456 532 #Store the results
457   - Res = ResultsSorage()
  533 + Res = ResultsStorage()
458 534  
459 535 #check the source name, exit is the srcname is false
460 536 analyse.GetModel(srcname)
... ... @@ -493,10 +569,11 @@ def MakeLC(analyse,nbins,srcname,binned = False):
493 569 return Res
494 570  
495 571 class UpperLimitComputer(base):
  572 + '''This class compute the ul and will be removed'''
496 573 def __init__(self,analyser,srcname,parname = "Prefactor"):
497 574 super(UpperLimitComputer,self).__init__()
498 575 self.info("Creating "+self.classname+" object")
499   - # self.like = analyser.like.copy()
  576 +
500 577 obs = analyser.like.obs().copy()
501 578 self.like = ct.ctlike(obs)
502 579 self.succes = False
... ...
pyext/ctools/plotutils.py
  View file @d113033
... ... @@ -56,7 +56,7 @@ class SpectralBase(object):
56 56 def __init__(self,sed=True,sed_factor=1.):
57 57 self.m_sed_factor = sed_factor
58 58 self.m_sed = sed
59   -
  59 +
60 60 def _convertSED(self,tab,energy):
61 61 """convert the a table storing dN/dE in sed with right units"""
62 62  
... ... @@ -68,6 +68,26 @@ class SpectralBase(object):
68 68 tab[i] *= energy[i]**2*self.m_sed_factor
69 69 return tab
70 70  
  71 + def _convertE(self,energy,unit,newunit):
  72 + """convert the a energy which has units 'unit' to "newunit'"""
  73 + factor = 1
  74 + if unit == "MeV" and newunit == "GeV" :
  75 + factor = 1e-3
  76 + if unit == "MeV" and newunit == "TeV" :
  77 + factor = 1e-6
  78 + if unit == "GeV" and newunit == "TeV" :
  79 + factor = 1e-3
  80 + if unit == "GeV" and newunit == "MeV" :
  81 + factor = 1e3
  82 + if unit == "TeV" and newunit == "GeV" :
  83 + factor = 1e3
  84 + if unit == "TeV" and newunit == "MeV" :
  85 + factor = 1e6
  86 + for i in xrange(len(energy)):
  87 + energy[i] *= factor
  88 + return energy
  89 +
  90 +
71 91 BLUE="#337EF5"
72 92 GREY='grey'
73 93 GREEN = '#21851D'
... ... @@ -103,6 +123,7 @@ class residuals(analysisutils.base,options,SpectralBase):
103 123 Energy_list = {"MeV":1,"GeV":1e-3,"TeV":1e-6} #TODO
104 124 for i in range(self.Npt):
105 125 if self.m_point_info["TS"][i] >= self.tslim:
  126 + #TODO remove the like
106 127 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 128 index.append(i)
108 129 else:
... ... @@ -139,7 +160,7 @@ class ulimgraph(analysisutils.base,options,SpectralBase):
139 160 self.Npt = len(self.m_point_info["dnde"]["value"])
140 161  
141 162 def _plt_points(self):
142   - results = analysisutils.ResultsSorage()
  163 + results = analysisutils.ResultsStorage()
143 164  
144 165 Energy_list = {"MeV":1,"GeV":1e-3,"TeV":1e-6} #TODO
145 166 for i in range(self.Npt):
... ... @@ -232,22 +253,19 @@ class SpectrumPlotter(analysisutils.base,SpectralBase):
232 253 if data points are provided, only change applied is energy unit convertion to match the butterfly one
233 254 (i.e. no dN/dE to SED conversion).
234 255 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):
  256 + once the I/O container format will be defined.
  257 + like is used to compute the residual. if like = None, no residual will be computed TODO
  258 + butterfly is a results Storage object read from the file produced by ctbutterfly
  259 + datapoint is the results Storage object for the datapoints."""
  260 + def __init__(self,srcname,like, butterfly,datapoint,emin=0.1,emax=100,energy = "TeV",SED = True):
237 261 super(SpectrumPlotter,self).__init__()
  262 + self.m_ebound = [emin,emax]
  263 + self.m_name = srcname
  264 + self.m_butterfly = butterfly
  265 + self.m_spectral = like.obs().models()[srcname].spectral()
  266 + self.m_but = []
  267 + self.m_enebut = []
238 268  
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 269 #~ this part take care of the style : SED or differential flux
252 270 #~ in the case of SED, the y-unit is erg/cm2/s
253 271 self.m_sed = SED
... ... @@ -260,6 +278,7 @@ class SpectrumPlotter(analysisutils.base,SpectralBase):
260 278 self.points = ulimgraph(datapoint,self.m_sed,gl.MeV2erg,energy)
261 279 self.residuals = residuals(self.m_spectral,datapoint,self.m_sed,gl.MeV2erg,energy)
262 280  
  281 +
263 282 def _validateUnits(self):
264 283 Energy_list = {"MeV":1,"GeV":1e3,"TeV":1e6}
265 284 if not(Energy_list.has_key(self.eunit)):
... ... @@ -269,74 +288,24 @@ class SpectrumPlotter(analysisutils.base,SpectralBase):
269 288 if self.m_sed:
270 289 self.m_sed_factor = gl.MeV2erg*Energy_list[self.eunit]**2
271 290  
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 291 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   -
  292 + #convert the flux in sed if asked
  293 + self.m_energy = self._convertE(self.m_butterfly['ener']["value"],"MeV","TeV") #gtbutterfly gives MeV
  294 + self.m_flux = self._convertSED(self.m_butterfly['dnde']["value"],self.m_butterfly['ener']["value"])
  295 + self.m_error_d = self._convertSED(self.m_butterfly['dnde']["ed_value"],self.m_butterfly['ener']["value"])
  296 + self.m_error_u = self._convertSED(self.m_butterfly['dnde']["eu_value"],self.m_butterfly['ener']["value"])
  297 + self._makebutterfly() # make the butterfly
  298 +
  299 + self.success("Butterfly ready to be plotted")
  300 +#
332 301 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])
  302 + """make the butterfly by reading the result and making a closed table"""
  303 + for i in xrange(len(self.m_energy)):
  304 + self.m_but.append(self.m_error_u[i])
336 305 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])
  306 + for i in xrange(len(self.m_energy)):
  307 + idx = len(self.m_energy)-i-1
  308 + self.m_but.append(self.m_error_d[idx])
340 309 self.m_enebut.append(self.m_energy[idx])
341 310 self.m_but.append(self.m_but[0])
342 311 self.m_enebut.append(self.m_enebut[0])
... ...