jean perrot / ctools

Commit 7c45f419e3fa08cfa09b6d29fe9343a726a4e061

Authored by Jürgen Knödlseder
1 parent e3a962c3

Integrate script

Inline Side-by-side
Showing 1 changed file with 26 additions and 27 deletions
scripts/csspec.py
  View file @7c45f41
1 #! /usr/bin/env python 1 #! /usr/bin/env python
2 # ========================================================================== 2 # ==========================================================================
3 -# spectral points generation script. 3 +# Spectral points generation script.
4 # 4 #
5 # Copyright (C) 2014-2015 Michael Mayer 5 # Copyright (C) 2014-2015 Michael Mayer
6 # 6 #
@@ -142,13 +142,13 @@ class csspec(ctools.cscript): @@ -142,13 +142,13 @@ class csspec(ctools.cscript):
142 self.m_ebounds = self.create_ebounds() 142 self.m_ebounds = self.create_ebounds()
143 143
144 # Read other parameters 144 # Read other parameters
145 - self.m_outfile = self["outfile"].filename() 145 + self.m_outfile = self["outfile"].filename()
146 146
147 # Get other parameeters 147 # Get other parameeters
148 self.m_calc_ulimit = self["calc_ulim"].boolean() 148 self.m_calc_ulimit = self["calc_ulim"].boolean()
149 - self.m_calc_ts = self["calc_ts"].boolean()  
150 - self.m_fix_bkg = self["fix_bkg"].boolean()  
151 - self.m_fix_srcs = self["fix_srcs"].boolean() 149 + self.m_calc_ts = self["calc_ts"].boolean()
  150 + self.m_fix_bkg = self["fix_bkg"].boolean()
  151 + self.m_fix_srcs = self["fix_srcs"].boolean()
152 152
153 # Set some fixed parameters 153 # Set some fixed parameters
154 self.m_log = False # Logging in client tools 154 self.m_log = False # Logging in client tools
@@ -237,31 +237,30 @@ class csspec(ctools.cscript): @@ -237,31 +237,30 @@ class csspec(ctools.cscript):
237 table.extname("SPECTRUM") 237 table.extname("SPECTRUM")
238 238
239 # Add Header for compatibility with gammalib.GMWLSpectrum 239 # Add Header for compatibility with gammalib.GMWLSpectrum
240 - table.card("Instrument","CTA","Name of Instrument")  
241 - table.card("Telescope","CTA","Name of Telescope") 240 + table.card("Instrument", "CTA", "Name of Instrument")
  241 + table.card("Telescope", "CTA", "Name of Telescope")
242 242
243 # Create FITS table columns 243 # Create FITS table columns
244 - energy = gammalib.GFitsTableDoubleCol("energy",self.m_ebounds.size()) 244 + energy = gammalib.GFitsTableDoubleCol("Energy", self.m_ebounds.size())
245 energy.unit("TeV") 245 energy.unit("TeV")
246 - energy_err = gammalib.GFitsTableDoubleCol("energy_err",self.m_ebounds.size()) 246 + energy_err = gammalib.GFitsTableDoubleCol("e_Energy", self.m_ebounds.size())
247 energy_err.unit("TeV") 247 energy_err.unit("TeV")
248 - flux = gammalib.GFitsTableDoubleCol("flux",self.m_ebounds.size()) 248 + flux = gammalib.GFitsTableDoubleCol("Flux", self.m_ebounds.size())
249 flux.unit("erg/cm2/s") 249 flux.unit("erg/cm2/s")
250 - flux_err = gammalib.GFitsTableDoubleCol("flux_err",self.m_ebounds.size()) 250 + flux_err = gammalib.GFitsTableDoubleCol("e_Flux", self.m_ebounds.size())
251 flux_err.unit("erg/cm2/s") 251 flux_err.unit("erg/cm2/s")
252 - TSvalues = gammalib.GFitsTableDoubleCol("TS",self.m_ebounds.size())  
253 - ulim_values = gammalib.GFitsTableDoubleCol("ulimit",self.m_ebounds.size()) 252 + TSvalues = gammalib.GFitsTableDoubleCol("TS", self.m_ebounds.size())
  253 + ulim_values = gammalib.GFitsTableDoubleCol("UpperLimit", self.m_ebounds.size())
254 ulim_values.unit("erg/cm2/s") 254 ulim_values.unit("erg/cm2/s")
255 255
256 - # MeV2erg  
257 - MeV2erg = 1.60217e-6  
258 - 256 + # Loop over energy bins
259 for i in range(self.m_ebounds.size()): 257 for i in range(self.m_ebounds.size()):
260 -  
261 - emin = self.m_ebounds.emin(i)  
262 - emax = self.m_ebounds.emax(i)  
263 - emean = self.m_ebounds.emean(i)  
264 - elogmean = self.m_ebounds.elogmean(i) 258 +
  259 + # Get energy boundaries
  260 + emin = self.m_ebounds.emin(i)
  261 + emax = self.m_ebounds.emax(i)
  262 + emean = self.m_ebounds.emean(i)
  263 + elogmean = self.m_ebounds.elogmean(i)
265 elogmean2 = elogmean.MeV() * elogmean.MeV() 264 elogmean2 = elogmean.MeV() * elogmean.MeV()
266 265
267 # Select events 266 # Select events
@@ -281,7 +280,7 @@ class csspec(ctools.cscript): @@ -281,7 +280,7 @@ class csspec(ctools.cscript):
281 280
282 # Get results 281 # Get results
283 fitted_models = like.obs().models() 282 fitted_models = like.obs().models()
284 - source = fitted_models[self.m_srcname] 283 + source = fitted_models[self.m_srcname]
285 284
286 # Calculate Upper Limit 285 # Calculate Upper Limit
287 ulimit_value = -1.0 286 ulimit_value = -1.0
@@ -300,9 +299,9 @@ class csspec(ctools.cscript): @@ -300,9 +299,9 @@ class csspec(ctools.cscript):
300 fitted_flux = source.spectral().eval(elogmean,gammalib.GTime()) 299 fitted_flux = source.spectral().eval(elogmean,gammalib.GTime())
301 300
302 # Compute flux error 301 # Compute flux error
303 - parvalue = source.spectral()[0].value() 302 + parvalue = source.spectral()[0].value()
304 rel_error = source.spectral()[0].error()/parvalue 303 rel_error = source.spectral()[0].error()/parvalue
305 - e_flux = fitted_flux*rel_error 304 + e_flux = fitted_flux*rel_error
306 305
307 # Set values for storage 306 # Set values for storage
308 TSvalues[i] = TS 307 TSvalues[i] = TS
@@ -314,10 +313,10 @@ class csspec(ctools.cscript): @@ -314,10 +313,10 @@ class csspec(ctools.cscript):
314 energy_err[i] = (emax - emean).TeV() 313 energy_err[i] = (emax - emean).TeV()
315 314
316 # Convert fluxes to nuFnu 315 # Convert fluxes to nuFnu
317 - flux[i] = fitted_flux * elogmean2 * MeV2erg  
318 - flux_err[i] = e_flux * elogmean2 * MeV2erg 316 + flux[i] = fitted_flux * elogmean2 * gammalib.MeV2erg
  317 + flux_err[i] = e_flux * elogmean2 * gammalib.MeV2erg
319 if ulimit_value > 0.0: 318 if ulimit_value > 0.0:
320 - ulim_values[i] = ulimit_value * elogmean2 * MeV2erg 319 + ulim_values[i] = ulimit_value * elogmean2 * gammalib.MeV2erg
321 320
322 321
323 # Log information 322 # Log information