Leonardo Di Venere / ctools

Commit 5a7b1baff23b7df20e5e0991c7c27ee0055ea048

Authored by sanchez
1 parent 121c7b21

add the LC

Inline Side-by-side
Showing 2 changed files with 200 additions and 48 deletions
pyext/ctools/analysisutils.py
  View file @5a7b1ba
... ... @@ -20,7 +20,7 @@
20 20 # ==========================================================================
21 21 import ctools as ct
22 22 import gammalib as gl
23   -from math import log10,floor
  23 +from math import sqrt,log10,floor
24 24 import sys
25 25  
26 26 class base(object):
... ... @@ -74,10 +74,15 @@ class ResultsSorage(dict,base):
74 74 base.__init__(self)
75 75 self["dnde"] = {}
76 76 self["ulim_dnde"] = []
  77 + self["time"] = {}
77 78 self["ener"] = {}
78 79 self["TS"] = []
79 80 self["Iflux"] = {}
80 81  
  82 + self["time"]["tmin_value"] = []
  83 + self["time"]["tmax_value"] = []
  84 + self["time"]["unit"] = "sec"
  85 +
81 86 self["ener"]["value"] = []
82 87 self["ener"]["ed_value"] = []
83 88 self["ener"]["eu_value"] = []
... ... @@ -111,39 +116,98 @@ class Analyser(base):
111 116 self.m_obs = None
112 117 self.like = None
113 118  
114   - self.m_evtfile = "events_selected.fits"
115   - self.m_raw_evt = "events.fits"
116   - self.m_cntfile = "cntmap.fits"
117 119 self.m_ebinalg = "LOG"
118 120 self.m_eunit = "TeV"
119   - self.m_ra = 83.63
120   - self.m_dec = 22.01
121   - self.m_roi = 2.
122   - self.m_tmin = 0.
123   - self.m_tmax = 0.
124   - self.m_emin = 0.1#emin
125   - self.m_emax = 10.#emax
126 121 self.m_enumbins = 10#nbins
127 122 self.m_usepnt = True # Use pointing for map centre
128   - self.m_nxpix = 200#npix
129   - self.m_nypix = 200#npix
130   - self.m_binsz = 0.05#binsz
131   - self.m_coordsys = "GAL"#coord
  123 + self.m_coordsys = "CEL"#coord
132 124 self.m_proj = "TAN"#proj
133 125 self.m_binned = False
134   -
135   - self.m_stat = "POISSON"#stat
136   - self.m_caldb = "dummy"
137   - self.m_irf = "cta_dummy_irf"
138 126 self.m_edisp = False
139   - self.m_xml = "$PWD/crab.xml"
140   - self.m_xml_out = "$PWD/crab_results.xml"
  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)#npix
  147 + self.m_nypix = int(self.m_roi*100)#npix
  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)
141 198  
142 199  
143 200 def set_obs(self,obs):
144 201 self.m_obs = obs #set the GObservation container
145 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 +
146 209 def set_energy_boundary(self,emin,emax):
  210 + """Set the energy bounds of the observation"""
147 211 self.m_emin = emin
148 212 self.m_emax = emax
149 213 ebounds = gl.GEbounds(gl.GEnergy(emin, self.m_eunit), \
... ... @@ -152,10 +216,12 @@ class Analyser(base):
152 216 if self.m_obs:
153 217 for obs in self.m_obs:
154 218 obs.events().ebounds(ebounds)
  219 + self.m_energy_set = True
155 220  
156 221 def ctselect(self):
157 222 # ctselect application and set parameters
158 223 self.info("Running ctselect to cut on events")
  224 + self.validate()
159 225 if self.m_obs:
160 226 filter = ct.ctselect(self.m_obs)
161 227 else:
... ... @@ -186,6 +252,7 @@ class Analyser(base):
186 252  
187 253 def ctbin(self,log=False,debug=False):
188 254 # ctbin application and set parameters
  255 + self.validate()
189 256 self.info("Running ctbin to create count map")
190 257 if self.m_obs:
191 258 bin = ct.ctbin(self.m_obs)
... ... @@ -215,7 +282,7 @@ class Analyser(base):
215 282  
216 283 # Run ctbin application. This will loop over all observations in
217 284 # the container and bin the events in counts maps
218   - bin.excute()
  285 + bin.execute()
219 286 if self.m_obs:
220 287 # Make a deep copy of the observation that will be returned
221 288 # (the ctbin object will go out of scope one the function is
... ... @@ -223,6 +290,7 @@ class Analyser(base):
223 290 self.m_obs = bin.obs().copy()
224 291  
225 292 def create_fit(self,log=False,debug=False):
  293 + self.validate()
226 294 # create ctlike instance with given parameters
227 295 self.info("Fitting Data using ctlike")
228 296 if self.m_obs:
... ... @@ -237,7 +305,7 @@ class Analyser(base):
237 305 self.like["stat"].string(self.m_stat)
238 306 self.like["caldb"].string(self.m_caldb)
239 307 self.like["irf"].string(self.m_irf)
240   - self.like["srcmdl"] = self.m_xml
  308 + self.like["srcmdl"] = self.m_xml
241 309 self.like["edisp"].boolean(self.m_edisp)
242 310 self.like["outmdl"] = self.m_xml_out
243 311  
... ... @@ -249,6 +317,7 @@ class Analyser(base):
249 317 self.like["debug"].boolean(True)
250 318  
251 319 def fit(self,log=False,debug=False):
  320 + self.validate()
252 321 if not(self.like):
253 322 self.warning("ctlike object not created, creating now")
254 323 self.create_fit(log,debug)
... ... @@ -266,15 +335,33 @@ class Analyser(base):
266 335 print "Model : "+m.name()
267 336 print m
268 337  
269   - def GetSrcResuls(self,Res,srcname,E0=1,factor = 1e6):
270   -
  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"""
271 340 results = self.GetSrcParams(srcname)
272 341  
273   - Res["ener"]["value"].append(E0*factor)
274   - Res["ener"]["ed_value"].append((E0-self.m_emin)*factor)
275   - Res["ener"]["eu_value"].append((self.m_emax-E0)*factor)
276   - Res["ener"]["unit"] = "MeV"
277   -
  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 +
278 365 Res["Iflux"]["value"].append(results["flux"])
279 366 Res["Iflux"]["eu_value"].append(results["eflux"])
280 367 Res["Iflux"]["ed_value"].append(results["eflux"])
... ... @@ -284,15 +371,8 @@ class Analyser(base):
284 371 Res["dnde"]["value"].append(m.spectral()[0].value())
285 372 Res["dnde"]["eu_value"].append(m.spectral()[0].error())
286 373 Res["dnde"]["ed_value"].append(m.spectral()[0].error())
287   - # Res["dnde"]["value"].append(m.spectral().eval(gl.GEnergy(E0,self.m_eunit),gl.GTime(0)))
288   - # m.spectral().eval_gradients(gl.GEnergy(E0,self.m_eunit),gl.GTime(0))
289   - # Assume that par 0 is the norm
290   - # Res["dnde"]["eu_value"].append(m.spectral()[0].factor_gradient())
291   - # Res["dnde"]["ed_value"].append(m.spectral()[0].factor_gradient())
292   - # Res["dnde"]["unit"] = "unit"
293 374  
294 375 Res["TS"].append(100)
295   - # Res["TS"].append(self.like.obs().models()[srcname].ts())
296 376  
297 377 UL = UpperLimitComputer(self,srcname)
298 378 UL.run()
... ... @@ -339,23 +419,20 @@ def MakeEbin(analyse,nbins,srcname,binned = False):
339 419 E0 = pow(10, (log10(emin) + log10(emax)) / 2)
340 420 print "Making energy bin between %2.1e"%(emin)+analyse.m_eunit+" and %2.1e"%(emax)+analyse.m_eunit
341 421  
342   - # change the PivotEnergy to E0, froze all spectral parameters but the prefactor
  422 + # change the PivotEnergy to E0
343 423 m = obs_copy.models()[srcname]
344 424 for par in m.spectral():
345 425 if par.name() == "PivotEnergy":
346 426 par.value(E0*factor)
347   - # if not(par.name() == "Prefactor") and par.is_free(): # assume that all model have a prefactor
348   - # par.fix()
349   - # if par.name() == "Prefactor":
350   - # Set the norm to a better guess
351   - # new_val = m.spectral().eval(gl.GEnergy(E0,analyse.m_eunit),gl.GTime(0))
352   - # m.spectral()["Prefactor"].value(new_val)
353 427  
354 428 ana_bin = Analyser()
355 429  
356 430 #copy the obs object
357 431 ana_bin.set_obs(obs_copy)
358 432  
  433 + #copy the analyser parameters
  434 + ana_bin.copy(analyse)
  435 +
359 436 # Set energy boundaries
360 437 ana_bin.set_energy_boundary(emin,emax)
361 438 ana_bin.ctselect()
... ... @@ -374,6 +451,46 @@ def MakeEbin(analyse,nbins,srcname,binned = False):
374 451  
375 452 del ana_bin
376 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
377 494  
378 495 class UpperLimitComputer(base):
379 496 def __init__(self,analyser,srcname,parname = "Prefactor"):
... ... @@ -437,7 +554,6 @@ class UpperLimitComputer(base):
437 554 value = self.spec[self.parname].value()
438 555 error = self.spec[self.parname].error()
439 556 scale = self.spec[self.parname].scale()
440   - print value," ",error," ",scale
441 557  
442 558 #Fix all parameters
443 559 for model in self.like.obs().models():
... ... @@ -451,7 +567,6 @@ class UpperLimitComputer(base):
451 567 # set start (a) and (b) depend on the CL asked
452 568 a= value / scale + (self.dloglike/4.)*error/scale
453 569 b=value / scale + (self.dloglike*2)* error/scale
454   - print a," ",b," ",self._func(a)," ",self._func(b)
455 570 ul = self._bisec(a,b)
456 571  
457 572 self.ulimit = ul*self.spec[self.parname].scale()
... ...
pyext/ctools/plotutils.py
  View file @5a7b1ba
... ... @@ -189,6 +189,45 @@ class ulimgraph(analysisutils.base,options,SpectralBase):
189 189 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)
190 190  
191 191  
  192 +class LightCurvePlotter(analysisutils.base,options):
  193 + """ Class to plot the Light curve"""
  194 + def __init__(self,srcname,datapoint):
  195 + super(LightCurvePlotter,self).__init__()
  196 + options.__init__(self)
  197 + self.m_name = srcname
  198 + self.m_data = datapoint #data point in a Results Storage class
  199 +
  200 + def draw(self):
  201 + if not(has_matplotlib):
  202 + self.warning("matplotlib module not found, can draw")
  203 + return
  204 +
  205 + time = []
  206 + dtime = []
  207 + flux = []
  208 + dflux_ed = []
  209 + dflux_eu = []
  210 + for i in range(len(self.m_data["time"]["tmin_value"])):
  211 + time.append((self.m_data["time"]["tmin_value"][i]+self.m_data["time"]["tmax_value"][i])/2)
  212 + dtime.append((self.m_data["time"]["tmax_value"][i]-self.m_data["time"]["tmin_value"][i])/2)
  213 + flux.append(self.m_data["Iflux"]["value"][i])
  214 + dflux_eu.append(self.m_data["Iflux"]["eu_value"][i])
  215 + dflux_ed.append(self.m_data["Iflux"]["ed_value"][i])
  216 +
  217 + plt.figure("results",figsize=(12,7),edgecolor="w")
  218 + unit = self.m_data["time"]["unit"]
  219 + plt.xlabel("Time ["+ unit+"]",fontsize=15)
  220 + unit = self.m_data["Iflux"]["unit"]
  221 + plt.ylabel("E$^{2}$ dN/dE ["+unit+"]",fontsize=15)
  222 +
  223 + print time
  224 + print dtime
  225 + print flux
  226 + print dflux_ed
  227 + print dflux_eu
  228 + plt.errorbar(time,flux,xerr=[dtime,dtime],yerr=[dflux_ed,dflux_eu],fmt=self.fmt,color=self.markercolor,elinewidth=self.elinewidth)
  229 + plt.show()
  230 +
192 231 class SpectrumPlotter(analysisutils.base,SpectralBase):
193 232 """ Class to compute and plot the spectra (i.e best fit model and butterfly).
194 233 if data points are provided, only change applied is energy unit convertion to match the butterfly one
... ... @@ -222,7 +261,6 @@ class SpectrumPlotter(analysisutils.base,SpectralBase):
222 261 self.points = ulimgraph(datapoint,self.m_sed,gl.MeV2erg,energy)
223 262 self.residuals = residuals(self.m_spectral,datapoint,self.m_sed,gl.MeV2erg,energy)
224 263  
225   -
226 264 def _validateUnits(self):
227 265 Energy_list = {"MeV":1,"GeV":1e3,"TeV":1e6}
228 266 if not(Energy_list.has_key(self.eunit)):
... ... @@ -292,7 +330,6 @@ class SpectrumPlotter(analysisutils.base,SpectralBase):
292 330 self._makebutterfly() # make the butterfly
293 331 self.success("Spectrum computed")
294 332  
295   -
296 333 def _makebutterfly(self):
297 334 """make the butterfly by appending element in a table"""
298 335 for i in xrange(self.m_npt):
... ...