diff --git a/scripts/Makefile.am b/scripts/Makefile.am
index 5511cdd..9c9a44a 100644
--- a/scripts/Makefile.am
+++ b/scripts/Makefile.am
@@ -5,12 +5,12 @@ binscripts = cta_kb_bkgrates.py cta_root2irf.py cta_make_nodes.py \
tsmerge.py
# Scripts to be distributed in pythondir
-pythonscripts = cspull.py cssens.py cstsdist.py cscaldb.py csresmap.py \
+pythonscripts = cspull.py cssens.py cstsdist.py cscaldb.py csresmap.py csspec.py\
lsmodel.py \
obsutils.py
# Scripts for which symbolic links will be created in bindir
-linkscripts = cspull cssens cstsdist cscaldb csresmap\
+linkscripts = cspull cssens cstsdist cscaldb csresmap csspec csplotspec\
lsmodel
# Distribute scripts
diff --git a/scripts/csplotspec.py b/scripts/csplotspec.py
new file mode 100755
index 0000000..d74fa24
--- /dev/null
+++ b/scripts/csplotspec.py
@@ -0,0 +1,345 @@
+#! /usr/bin/env python
+# ==========================================================================
+# spectral points generation script.
+#
+# Copyright (C) 2014-2015 Michael Mayer
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+#
+# ==========================================================================
+import gammalib
+import ctools
+import sys
+
+try:
+ import matplotlib.pyplot as plt
+except:
+ sys.exit("Need matplotlib to launch this script")
+
+# Initialise conversion constant
+MeV2erg = 1.60217e-6
+
+# ============== #
+# cstsdist class #
+# ============== #
+class csplotspec(ctools.cscript):
+ """
+ This class plots a copmplete spectrum. It derives from
+ the ctools.cscript class which provides support for parameter files,
+ command line arguments, and logging. In that way the Python script
+ behaves just as a regular ctool.
+ """
+ def __init__(self, *argv):
+ """
+ Constructor.
+ """
+
+ # Set name
+ self.name = "csplotspec"
+ self.version = "1.0.0"
+ self.outfile = ""
+
+ # Make sure that parfile exists
+ file = self.parfile()
+
+ # Initialise application
+ if len(argv) == 0:
+ ctools.cscript.__init__(self, self.name, self.version)
+ elif len(argv) ==1:
+ ctools.cscript.__init__(self, self.name, self.version, *argv)
+ else:
+ raise TypeError("Invalid number of arguments given.")
+
+ # Set logger properties
+ self.log_header()
+ self.log.date(True)
+
+ # Return
+ return
+
+ def __del__(self):
+ """
+ Destructor.
+ """
+ # Write separator into logger
+ if self.logTerse():
+ self.log("\n")
+
+ # Return
+ return
+
+ def parfile(self):
+ """
+ Check if parfile exists. If parfile does not exist then create a
+ default parfile. This kluge avoids shipping the cscript with a parfile.
+ """
+ # Set parfile name
+ parfile = self.name+".par"
+
+ try:
+ pars = gammalib.GApplicationPars(parfile)
+ except:
+ # Signal if parfile was not found
+ sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
+
+ # Create default parfile
+ pars = gammalib.GApplicationPars()
+ pars.append(gammalib.GApplicationPar("inmodel","f","a","results.xml","","","Source model"))
+ pars.append(gammalib.GApplicationPar("srcname","s","a","Crab","","","Source name"))
+ pars.append(gammalib.GApplicationPar("inspec","f","a","spectrum.fits","","","Input spectral points"))
+ pars.append(gammalib.GApplicationPar("inbfly","f","a","butterfly.txt","","","Input butterfly"))
+ pars.append(gammalib.GApplicationPar("outfile","f","a","spectrum.png","","","Output file name"))
+ pars.append(gammalib.GApplicationPar("tslim","r","h","4.0","","","Minimum TS value for data point"))
+ pars.append(gammalib.GApplicationPar("xmin","r","h","0.1","","","Minimum energy for plotting"))
+ pars.append(gammalib.GApplicationPar("xmax","r","h","100","","","Maximum energy for plotting"))
+ pars.append(gammalib.GApplicationPar("show","b","h","yes","yes/no","","Show plot directly"))
+ pars.append_standard()
+ pars.append(gammalib.GApplicationPar("logfile","f","h","csplotspec.log","","","Log filename"))
+ pars.save(parfile)
+
+ # Return
+ return
+
+ def get_parameters(self):
+ """
+ Get parameters from parfile and setup the observation.
+ """
+ # Get parameters
+ self.models = gammalib.GModels(self["inmodel"].filename())
+ self.m_source = self.models[self["srcname"].string()]
+ self.m_bfly = self["inbfly"].filename()
+ self.m_inspec = self["inspec"].filename()
+ self.m_outfile = self["outfile"].filename()
+
+ # Get other parameters
+ self.m_tslim = self["tslim"].real()
+ self.m_xmin = self["xmin"].real()
+ self.m_xmax = self["xmax"].real()
+ self.m_show = self["show"].boolean()
+
+ # Set some fixed parameters
+ self.m_log = False # Logging in client tools
+ self.m_chatter = self["chatter"].integer()
+ self.m_clobber = self["clobber"].boolean()
+ self.m_debug = self["debug"].boolean()
+
+ # Return
+ return
+
+
+ def get_points(self):
+
+ # Get FITS file
+ fits = gammalib.GFits(self.m_inspec)
+ table = fits["SPECTRUM"]
+
+ # Get FITS columns from file
+ c_energy = table["energy"]
+ c_flux = table["flux"]
+ c_flux_err = table["flux_err"]
+ c_TS = table["TS"]
+ c_ulim = table["ulimit"]
+
+ # Initialise arrays
+ energy = []
+ flux = []
+ ed_flux = []
+ eu_flux = []
+
+ # Upper limit arrays
+ lim_e = []
+ lim_flux = []
+ lim_edflux = []
+ lim_euflux = []
+
+ # Loop over points
+ for i in range(c_energy.length()):
+
+ # Check for significant point
+ if c_TS[i] < self.m_tslim:
+
+ # Use upper limit if existent
+ if c_ulim[i] > 0.0:
+
+ # Append value to array
+ lim_e.append(c_energy[i])
+ lim_flux.append(c_ulim[i])
+ lim_edflux.append(0.9*c_ulim[i])
+ lim_euflux.append(0.0)
+ else:
+ # If there is no upper limit available, this point is skipped
+ continue
+
+ else:
+
+ # Append significant data point to array
+ energy.append(c_energy[i])
+ flux.append(c_flux[i])
+ ed_flux.append(c_flux_err[i])
+ eu_flux.append(c_flux_err[i])
+
+ # Return plot values
+ return energy,flux,ed_flux,eu_flux, lim_e, lim_flux, lim_edflux, lim_euflux
+
+ def get_line(self, xmin = 0.1, xmax = 100.):
+
+ # Initialise energy binning
+ ebounds = gammalib.GEbounds(100, gammalib.GEnergy(xmin,"TeV"), gammalib.GEnergy(xmax,"TeV"), True)
+
+ # Initialise return array
+ x = []
+ y = []
+
+ # Initialise dummy time
+ t = gammalib.GTime()
+
+ # Loop over energy bins
+ for i in range(ebounds.size()):
+
+ # Append values to array
+ energy = ebounds.elogmean(i)
+ x.append(energy.TeV())
+ y.append(energy.MeV() * energy.MeV() * self.m_source.spectral().eval(energy, t) * MeV2erg)
+
+ # Return plot values
+ return x,y
+
+ def get_butterfly(self):
+
+ # Get ascii file
+ csv = gammalib.GCsv(self.m_bfly)
+
+ # Initialise return arrays
+ x = []
+ y = []
+
+ # Loop over rows
+ for i in range(csv.nrows()):
+
+ #store energy
+ energy = csv.real(i,0) # MeV
+
+ # store upper limit flux
+ flux = csv.real(i,1) + csv.real(i,2) / 2.0 #ph/cm2/s/MeV
+
+ # append to array
+ x.append(energy * 1e-6) # TeV
+ y.append(flux * energy * energy * MeV2erg) # erg/cm2/s
+
+ # Loop again over file from the other side
+ for i in range(csv.nrows()):
+
+ # Create index counting backwards
+ k = csv.nrows()-1-i
+
+ # Get values
+ energy = csv.real(k,0) # MeV
+ flux = csv.real(k,1)-csv.real(k,2) / 2.0 #ph/cm2/s/MeV
+
+ # Append values
+ x.append(energy * 1e-6) # TeV
+ y.append(flux * energy * energy * MeV2erg) # erg/cm2/s
+
+ # Return plot values
+ return x,y
+
+
+
+ def execute(self):
+ """
+ Execute the script.
+ """
+ # Run the script
+ self.run()
+
+ # Save image file
+ if not (self.m_outfile == "" or self.m_outfile == "NONE"):
+ plt.savefig(self.outfile)
+
+ # Return
+ return
+
+ def run(self):
+ """
+ Run the script.
+ """
+ # Switch screen logging on in debug mode
+ if self.logDebug():
+ self.log.cout(True)
+
+ # Get parameters
+ self.get_parameters()
+
+ # Write input parameters into logger
+ if self.logTerse():
+ self.log_parameters()
+ self.log("\n")
+
+ plt.figure("Spectrum")
+ plt.loglog()
+ plt.xlabel("Energy [TeV]")
+ plt.ylabel(r'$\nu$F$\nu$ [erg cm$^{-2}$ s$^{-1}$]')
+
+ # Write header
+ if self.logTerse():
+ self.log("\n")
+ self.log.header1("Get Butterfly")
+
+ if not (self.m_bfly == "NONE" or self.m_bfly == ""):
+ xbf, ybf = self.get_butterfly()
+ plt.fill(xbf,ybf,color='green',alpha=0.5)
+
+ # Write header
+ if self.logTerse():
+ self.log("\n")
+ self.log.header1("Get Spectrum Line")
+
+ xline,yline = self.get_line(self.m_xmin, self.m_xmax)
+ plt.plot(xline,yline, color='black',lw=2.0)
+
+ # Write header
+ if self.logTerse():
+ self.log("\n")
+ self.log.header1("Get Spectrum Points")
+
+ if not (self.m_inspec == "NONE" or self.m_inspec == ""):
+ e,f,edf,euf,le,lf,ledf,leuf = self.get_points()
+ plt.errorbar(x=e,y=f,yerr=[edf,euf],fmt='bo',lw=2.0)
+ if len(le):
+ plt.errorbar(x=le,y=lf,yerr=[ledf,leuf],fmt='bo',lw=2.0, capsize=2.0,lowlims=True)
+
+ if self.m_show:
+ plt.show()
+
+
+ # Return
+ return
+
+
+# ======================== #
+# Main routine entry point #
+# ======================== #
+if __name__ == '__main__':
+ """
+ Generates spectral points.
+ """
+ # Create instance of application
+ app = csplotspec(sys.argv)
+
+ # Open logfile
+ app.logFileOpen()
+
+ # Execute application
+ app.execute()
+
\ No newline at end of file
diff --git a/scripts/csspec.py b/scripts/csspec.py
index 6fe6c3d..2be7fd8 100755
--- a/scripts/csspec.py
+++ b/scripts/csspec.py
@@ -1,8 +1,8 @@
#! /usr/bin/env python
# ==========================================================================
-# Script for spectral data points generation.
+# spectral points generation script.
#
-# Copyright (C) 2014 Michael Mayer
+# Copyright (C) 2014-2015 Michael Mayer
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
@@ -18,33 +18,28 @@
# along with this program. If not, see .
#
# ==========================================================================
-import ctools
import gammalib
+import ctools
import sys
-import os
-
-
# ============== #
# cstsdist class #
# ============== #
-class csspec(ctools.ctool):
+class csspec(ctools.cscript):
"""
- This class implements the creation of spectral points. It derives
- from the ctools::ctool class which provides support for parameter
- files, command line arguments, and logging. In that way the Python
- script behaves just as a regular ctool.
+ This class implements the creation of spectral points. It derives from
+ the ctools.cscript class which provides support for parameter files,
+ command line arguments, and logging. In that way the Python script
+ behaves just as a regular ctool.
"""
def __init__(self, *argv):
"""
Constructor.
"""
- self.name = "csspec"
- self.version = "0.1.0"
- #super(ctools.ctool).__init__(self.name,self.version,argv)
+
# Set name
-# self.name = "csspectrum"
-# self.version = "0.1.0"
+ self.name = "csspec"
+ self.version = "1.0.0"
# Initialise some members
self.obs = None
@@ -55,18 +50,17 @@ class csspec(ctools.ctool):
argv = argv[1:]
else:
self.obs = gammalib.GObservations()
- self.obs.clear()
-
- self.outfile = ""
+ self.obs.clear()
+ self.m_outfile = ""
# Make sure that parfile exists
file = self.parfile()
# Initialise application
if len(argv) == 0:
- gammalib.GApplication.__init__(self, self.name, self.version)
+ ctools.cscript.__init__(self, self.name, self.version)
elif len(argv) ==1:
- gammalib.GApplication.__init__(self, self.name, self.version, *argv)
+ ctools.cscript.__init__(self, self.name, self.version, *argv)
else:
raise TypeError("Invalid number of arguments given.")
@@ -104,17 +98,24 @@ class csspec(ctools.ctool):
# Create default parfile
pars = gammalib.GApplicationPars()
- pars.append(gammalib.GApplicationPar("infile","f","a","obs.xml","","","Observation definition"))
- pars.append(gammalib.GApplicationPar("srcmdl","f","a","crab.xml","","","model file name"))
+ pars.append(gammalib.GApplicationPar("inobs","f","a","events.fits","","","Event list, counts cube, or observation definition file"))
+ pars.append(gammalib.GApplicationPar("inmodel","f","a","$CTOOLS/share/models/crab.xml","","","Source model"))
pars.append(gammalib.GApplicationPar("srcname","s","a","Crab","","","Source name"))
- pars.append(gammalib.GApplicationPar("outfile","f","h","spectrum.fits","","","Output file name"))
- pars.append(gammalib.GApplicationPar("caldb","s","a","$GAMMALIB/share/caldb/cta","","","Calibration database"))
+ pars.append(gammalib.GApplicationPar("outfile","f","a","spectrum.fits","","","Output file name"))
+ pars.append(gammalib.GApplicationPar("expcube","f","a","NONE","","","Exposure cube file (only needed for stacked analysis)"))
+ pars.append(gammalib.GApplicationPar("psfcube","f","a","NONE","","","PSF cube file (only needed for stacked analysis)"))
+ pars.append(gammalib.GApplicationPar("caldb","s","a","dummy","","","Calibration database"))
pars.append(gammalib.GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function"))
- pars.append(gammalib.GApplicationPar("emin","r","a","0.1","0.0","","Lower energy limit (TeV)"))
- pars.append(gammalib.GApplicationPar("emax","r","a","100.0","0.0","","Upper energy limit (TeV)"))
+ pars.append(gammalib.GApplicationPar("emin","r","h","0.1","","","Lower energy limit (TeV)"))
+ pars.append(gammalib.GApplicationPar("emax","r","h","100.0","","","Upper energy limit (TeV)"))
pars.append(gammalib.GApplicationPar("enumbins","i","a","20","","","Number of energy bins"))
pars.append(gammalib.GApplicationPar("ebinalg","s","h","LOG","LIN|LOG|FILE","","Binning algorithm"))
+ pars.append(gammalib.GApplicationPar("calc_ts","b","h","yes","yes|no","","Compute TS value in each bin"))
+ pars.append(gammalib.GApplicationPar("calc_ulim","b","h","yes","yes|no","","Compute upper limit in each bin"))
+ pars.append(gammalib.GApplicationPar("fix_srcs","b","h","yes","yes|no","","Fix other skymodel parameters"))
+ pars.append(gammalib.GApplicationPar("fix_bkg","b","h","no","yes|no","","Fix background parameters"))
pars.append_standard()
+ pars.append(gammalib.GApplicationPar("logfile","f","h","csspec.log","","","Log filename"))
pars.save(parfile)
# Return
@@ -125,59 +126,36 @@ class csspec(ctools.ctool):
Get parameters from parfile and setup the observation.
"""
# Get parameters
-
# Set observation if not done before
- if self.obs.size() == 0:
-
- obsfile = self["infile"].filename()
- try:
-
- self.obs = gammalib.GObservations(obsfile)
- except:
- self.obs.clear()
- observation = gammalib.GCTAObservation()
- observation.load(obsfile)
-
- self.m_irf = self["irf"].string()
- self.m_caldb = self["caldb"].string()
-
- caldb = gammalib.GCaldb()
- if os.path.isdir(self.m_caldb):
- caldb.rootdir(self.m_caldb)
- else:
- caldb.open("cta",self. m_caldb)
+ if self.obs == None or self.obs.size() == 0:
+ self.require_inobs("csspec::get_parameters()")
+ self.obs = self.get_observations()
- observation.response(self.m_irf, caldb);
- self.obs.append(observation)
-
- # Check for models in the container
- # read models from file if there were none
- # Note that the models should be optimised for
- # meaningful residual map
+ # Set models if we have none
if self.obs.models().size() == 0:
- self.models = gammalib.GModels(self["srcmdl"].filename())
-
- self.m_srcname = self["srcname"].string()
- self.m_outfile = self["outfile"].filename()
-
-
- # not sure how to call get_ebounds as private member of ctool
- #self.ebounds = ctools.ctool.get_ebounds()
- self.m_emin = self["emin"].real()
- self.m_emax = self["emax"].real()
- self.m_enumbins = self["enumbins"].integer()
- #self.m_ebinalg = self["ebinalg"].string()
+ self.obs.models(self["inmodel"].filename())
+
+ # Get source name
+ self.m_srcname = self["srcname"].string()
- self.ebounds = gammalib.GEbounds(self.m_enumbins,gammalib.GEnergy(self.m_emin,"TeV"),gammalib.GEnergy(self.m_emax,"TeV"),True)
+ # Get ebounds
+ self.m_ebounds = self.create_ebounds()
+ # Read other parameters
+ self.m_outfile = self["outfile"].filename()
+
+ # Get other parameeters
+ self.m_calc_ulimit = self["calc_ulim"].boolean()
+ self.m_calc_ts = self["calc_ts"].boolean()
+ self.m_fix_bkg = self["fix_bkg"].boolean()
+ self.m_fix_srcs = self["fix_srcs"].boolean()
# Set some fixed parameters
- self.m_log = False # Logging in client tools
- self.m_debug = False # Debugging in client tools
+ self.m_log = False # Logging in client tools
+ self.m_chatter = self["chatter"].integer()
self.m_clobber = self["clobber"].boolean()
-
- self.parname =self.models[self.m_srcname].spectral()[0].name()
-
+ self.m_debug = self["debug"].boolean()
+
# Return
return
@@ -199,49 +177,21 @@ class csspec(ctools.ctool):
self.run()
# Save residual map
- self.save(self.m_outfile, self.m_clobber)
+ self.fits.saveto(self.m_outfile, self.m_clobber)
# Return
return
- def analyse_bin(self, emin,emax):
- #print "1",self.models.size()
- select = ctools.ctselect(self.obs)
- select["emin"].real(emin)
- select["emax"].real(emax)
- select["ra"].real(-1.0)
- select["dec"].real(-1.0)
- select["rad"].real(-1.0)
- select["tmin"].real(0.0)
- select["tmax"].real(0.0)
- select.run()
- #print "2",self.models.size()
- like = ctools.ctlike(select.obs())
- like.run()
- optmodels = like.obs().models()
- #print "3",self.models.size()
- return optmodels.copy()
-
def run(self):
"""
Run the script.
"""
-
# Switch screen logging on in debug mode
if self.logDebug():
self.log.cout(True)
# Get parameters
self.get_parameters()
-
-
- for model in self.models:
- for par in model:
- par.fix()
-
- self.models[self.m_srcname][self.parname].free()
- self.models[self.m_srcname].tscalc(True)
- self.obs.models(self.models)
# Write input parameters into logger
if self.logTerse():
@@ -258,46 +208,151 @@ class csspec(ctools.ctool):
# Write header
if self.logTerse():
self.log("\n")
- self.log.header1("Generate binned map")
+ self.log.header1("Adjust model parameters")
+
+ for model in self.obs.models():
+ model.tscalc(False)
+ if model.name() == self.m_srcname:
+ for par in model:
+ par.fix()
+ model.spectral()[0].free()
+ if self.m_calc_ts:
+ model.tscalc(True)
+
+ elif self.m_fix_bkg and not model.classname() == "GModelData":
+ for par in model:
+ par.fix()
+
+ elif self.m_fix_srcs and model.classname() == "GModelSky":
+ for par in model:
+ par.fix()
+
+ # Write header
+ if self.logTerse():
+ self.log("\n")
+ self.log.header1("Generate spectrum")
+
+ # Initialise FITS Table with extension "SPECTRUM"
+ table = gammalib.GFitsBinTable(self.m_ebounds.size())
+ table.extname("SPECTRUM")
- binned=False
+ # Add Header for compatibility with gammalib.GMWLSpectrum
+ table.card("Instrument","CTA","Name of Instrument")
+ table.card("Telescope","CTA","Name of Telescope")
+ # Create FITS table columns
+ energy = gammalib.GFitsTableDoubleCol("energy",self.m_ebounds.size())
+ energy.unit("TeV")
+ energy_err = gammalib.GFitsTableDoubleCol("energy_err",self.m_ebounds.size())
+ energy_err.unit("TeV")
+ flux = gammalib.GFitsTableDoubleCol("flux",self.m_ebounds.size())
+ flux.unit("erg/cm2/s")
+ flux_err = gammalib.GFitsTableDoubleCol("flux_err",self.m_ebounds.size())
+ flux_err.unit("erg/cm2/s")
+ TSvalues = gammalib.GFitsTableDoubleCol("TS",self.m_ebounds.size())
+ ulim_values = gammalib.GFitsTableDoubleCol("ulimit",self.m_ebounds.size())
+ ulim_values.unit("erg/cm2/s")
- for i in range(self.ebounds.size()):
+ # MeV2erg
+ MeV2erg = 1.60217e-6
- emin = self.ebounds.emin(i).TeV()
- emax = self.ebounds.emax(i).TeV()
- emean = self.ebounds.elogmean(i)
- optmodels = self.analyse_bin(emin,emax)
- #print optmodels
- #if optmodels.size()>0:
- optmodel = optmodels[self.m_srcname]
- ts = optmodel.ts()
+ for i in range(self.m_ebounds.size()):
- flux = optmodel.spectral().eval(emean,gammalib.GTime())
- rel_error = optmodel[self.parname].error()/optmodel[self.parname].value()
+ emin = self.m_ebounds.emin(i)
+ emax = self.m_ebounds.emax(i)
+ emean = self.m_ebounds.emean(i)
+ elogmean = self.m_ebounds.elogmean(i)
+ elogmean2 = elogmean.MeV() * elogmean.MeV()
- flux_error = flux * rel_error
+ # Select events
+ select = ctools.ctselect(self.obs)
+ select["emin"].real(emin.TeV())
+ select["emax"].real(emax.TeV())
+ select["tmin"].real(0.0)
+ select["tmax"].real(0.0)
+ select["rad"].real(0.0)
+ select["ra"].real(0.0)
+ select["dec"].real(-1.0)
+ select.run()
- print emean,flux,flux_error,ts
-
- if binned:
- print "not implemented"
+ # likelihood
+ like = ctools.ctlike(select.obs())
+ like.run()
+
+ # Get results
+ fitted_models = like.obs().models()
+ source = fitted_models[self.m_srcname]
+
+ # Calculate Upper Limit
+ ulimit_value = -1.0
+ #if self.m_calc_ulimit:
+ # ulimit = ctools.ctulimit(like.obs())
+ # ulimit["eref"].real(elogmean.TeV())
+ # ulimit.run()
+ # ulimit_value.ulimit.diff_ulimit()
+
+ # Get TS value
+ TS = -1.0
+ if self.m_calc_ts:
+ TS = source.ts()
+
+ # Get differential flux
+ fitted_flux = source.spectral().eval(elogmean,gammalib.GTime())
+
+ # Compute flux error
+ parvalue = source.spectral()[0].value()
+ rel_error = source.spectral()[0].error()/parvalue
+ e_flux = fitted_flux*rel_error
+
+ # Set values for storage
+ TSvalues[i] = TS
+
+ # Store energy as TeV
+ energy[i] = elogmean.TeV()
+
+ # Store energy error
+ energy_err[i] = (emax - emean).TeV()
+
+ # Convert fluxes to nuFnu
+ flux[i] = fitted_flux * elogmean2 * MeV2erg
+ flux_err[i] = e_flux * elogmean2 * MeV2erg
+ if ulimit_value > 0.0:
+ ulim_values[i] = ulimit_value * elogmean2 * MeV2erg
- self.obs.models(self.models)
+
+ # Log information
+ if self.logExplicit():
+ self.log("Bin "+str(i)+" ["+str(emin.TeV())+"-"+str(emax.TeV())+"] TeV: ")
+ self.log("Flux = "+str(flux[i]))
+ self.log(" +- "+str(flux_err[i])+" [erg/cm2/s]")
+ if self.m_calc_ts and TSvalues[i] > 0.0:
+ self.log(", TS = "+str(TS))
+ if self.m_calc_ulimit and ulim_values[i] > 0.0:
+ self.log(", UL = "+str(ulim_values[i])+" [erg/cm2/s]")
+ self.log("\n")
+
+ # Append filles columns to fits table
+ table.append(energy)
+ table.append(energy_err)
+ table.append(flux)
+ table.append(flux_err)
+ table.append(TSvalues)
+ table.append(ulim_values)
+ # Create the FITS file now
+ self.fits = gammalib.GFits()
+ self.fits.append(table)
+
# Return
return
- def save(self,outfile,clobber):
- return
# ======================== #
# Main routine entry point #
# ======================== #
if __name__ == '__main__':
"""
- Generates residual count map.
+ Generates spectral points.
"""
# Create instance of application
app = csspec(sys.argv)
diff --git a/test/test_cscripts.sh.in b/test/test_cscripts.sh.in
index c00baba..5a9cb61 100755
--- a/test/test_cscripts.sh.in
+++ b/test/test_cscripts.sh.in
@@ -160,3 +160,30 @@ else
exit 1
fi
$ECHO " ok"
+
+#
+# Test csspec
+# =============
+$ECHO -n "Test csspec: "
+chmod +x @abs_top_srcdir@/scripts/csspec.py
+@abs_top_srcdir@/scripts/csspec.py \
+ inobs="data/crab_events.fits" \
+ inmodel="data/crab.xml" \
+ caldb="irf" \
+ irf="cta_dummy_irf" \
+ enumbins=5 \
+ srcname="Crab" \
+ emin=0.1 \
+ emax=100.0 \
+ outfile="spectrum.fits" > csspec.log
+$ECHO -n "."
+if [ -s "spectrum.fits" ]
+then
+ $ECHO -n "."
+else
+ $ECHO " spectrum.fits file is not valid"
+ exit 1
+fi
+$ECHO " ok"
+
+