Johan Bregeon / ctools

Commit 1e6cba7118f528ff798a0d7ca392e514b2a6a9e2

Authored by Michael Mayer
1 parent 3269286a

Added new csspec version and cscript to plot spectra

Inline Side-by-side
Showing 4 changed files with 552 additions and 125 deletions
scripts/Makefile.am
  View file @1e6cba7
... ... @@ -5,12 +5,12 @@ binscripts = cta_kb_bkgrates.py cta_root2irf.py cta_make_nodes.py \
5 5 tsmerge.py
6 6  
7 7 # Scripts to be distributed in pythondir
8   -pythonscripts = cspull.py cssens.py cstsdist.py cscaldb.py csresmap.py \
  8 +pythonscripts = cspull.py cssens.py cstsdist.py cscaldb.py csresmap.py csspec.py\
9 9 lsmodel.py \
10 10 obsutils.py
11 11  
12 12 # Scripts for which symbolic links will be created in bindir
13   -linkscripts = cspull cssens cstsdist cscaldb csresmap\
  13 +linkscripts = cspull cssens cstsdist cscaldb csresmap csspec csplotspec\
14 14 lsmodel
15 15  
16 16 # Distribute scripts
... ...
scripts/csplotspec.py 0 → 100755
  View file @1e6cba7
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# spectral points generation script.
  4 +#
  5 +# Copyright (C) 2014-2015 Michael Mayer
  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 gammalib
  22 +import ctools
  23 +import sys
  24 +
  25 +try:
  26 + import matplotlib.pyplot as plt
  27 +except:
  28 + sys.exit("Need matplotlib to launch this script")
  29 +
  30 +# Initialise conversion constant
  31 +MeV2erg = 1.60217e-6
  32 +
  33 +# ============== #
  34 +# cstsdist class #
  35 +# ============== #
  36 +class csplotspec(ctools.cscript):
  37 + """
  38 + This class plots a copmplete spectrum. It derives from
  39 + the ctools.cscript class which provides support for parameter files,
  40 + command line arguments, and logging. In that way the Python script
  41 + behaves just as a regular ctool.
  42 + """
  43 + def __init__(self, *argv):
  44 + """
  45 + Constructor.
  46 + """
  47 +
  48 + # Set name
  49 + self.name = "csplotspec"
  50 + self.version = "1.0.0"
  51 + self.outfile = ""
  52 +
  53 + # Make sure that parfile exists
  54 + file = self.parfile()
  55 +
  56 + # Initialise application
  57 + if len(argv) == 0:
  58 + ctools.cscript.__init__(self, self.name, self.version)
  59 + elif len(argv) ==1:
  60 + ctools.cscript.__init__(self, self.name, self.version, *argv)
  61 + else:
  62 + raise TypeError("Invalid number of arguments given.")
  63 +
  64 + # Set logger properties
  65 + self.log_header()
  66 + self.log.date(True)
  67 +
  68 + # Return
  69 + return
  70 +
  71 + def __del__(self):
  72 + """
  73 + Destructor.
  74 + """
  75 + # Write separator into logger
  76 + if self.logTerse():
  77 + self.log("\n")
  78 +
  79 + # Return
  80 + return
  81 +
  82 + def parfile(self):
  83 + """
  84 + Check if parfile exists. If parfile does not exist then create a
  85 + default parfile. This kluge avoids shipping the cscript with a parfile.
  86 + """
  87 + # Set parfile name
  88 + parfile = self.name+".par"
  89 +
  90 + try:
  91 + pars = gammalib.GApplicationPars(parfile)
  92 + except:
  93 + # Signal if parfile was not found
  94 + sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
  95 +
  96 + # Create default parfile
  97 + pars = gammalib.GApplicationPars()
  98 + pars.append(gammalib.GApplicationPar("inmodel","f","a","results.xml","","","Source model"))
  99 + pars.append(gammalib.GApplicationPar("srcname","s","a","Crab","","","Source name"))
  100 + pars.append(gammalib.GApplicationPar("inspec","f","a","spectrum.fits","","","Input spectral points"))
  101 + pars.append(gammalib.GApplicationPar("inbfly","f","a","butterfly.txt","","","Input butterfly"))
  102 + pars.append(gammalib.GApplicationPar("outfile","f","a","spectrum.png","","","Output file name"))
  103 + pars.append(gammalib.GApplicationPar("tslim","r","h","4.0","","","Minimum TS value for data point"))
  104 + pars.append(gammalib.GApplicationPar("xmin","r","h","0.1","","","Minimum energy for plotting"))
  105 + pars.append(gammalib.GApplicationPar("xmax","r","h","100","","","Maximum energy for plotting"))
  106 + pars.append(gammalib.GApplicationPar("show","b","h","yes","yes/no","","Show plot directly"))
  107 + pars.append_standard()
  108 + pars.append(gammalib.GApplicationPar("logfile","f","h","csplotspec.log","","","Log filename"))
  109 + pars.save(parfile)
  110 +
  111 + # Return
  112 + return
  113 +
  114 + def get_parameters(self):
  115 + """
  116 + Get parameters from parfile and setup the observation.
  117 + """
  118 + # Get parameters
  119 + self.models = gammalib.GModels(self["inmodel"].filename())
  120 + self.m_source = self.models[self["srcname"].string()]
  121 + self.m_bfly = self["inbfly"].filename()
  122 + self.m_inspec = self["inspec"].filename()
  123 + self.m_outfile = self["outfile"].filename()
  124 +
  125 + # Get other parameters
  126 + self.m_tslim = self["tslim"].real()
  127 + self.m_xmin = self["xmin"].real()
  128 + self.m_xmax = self["xmax"].real()
  129 + self.m_show = self["show"].boolean()
  130 +
  131 + # Set some fixed parameters
  132 + self.m_log = False # Logging in client tools
  133 + self.m_chatter = self["chatter"].integer()
  134 + self.m_clobber = self["clobber"].boolean()
  135 + self.m_debug = self["debug"].boolean()
  136 +
  137 + # Return
  138 + return
  139 +
  140 +
  141 + def get_points(self):
  142 +
  143 + # Get FITS file
  144 + fits = gammalib.GFits(self.m_inspec)
  145 + table = fits["SPECTRUM"]
  146 +
  147 + # Get FITS columns from file
  148 + c_energy = table["energy"]
  149 + c_flux = table["flux"]
  150 + c_flux_err = table["flux_err"]
  151 + c_TS = table["TS"]
  152 + c_ulim = table["ulimit"]
  153 +
  154 + # Initialise arrays
  155 + energy = []
  156 + flux = []
  157 + ed_flux = []
  158 + eu_flux = []
  159 +
  160 + # Upper limit arrays
  161 + lim_e = []
  162 + lim_flux = []
  163 + lim_edflux = []
  164 + lim_euflux = []
  165 +
  166 + # Loop over points
  167 + for i in range(c_energy.length()):
  168 +
  169 + # Check for significant point
  170 + if c_TS[i] < self.m_tslim:
  171 +
  172 + # Use upper limit if existent
  173 + if c_ulim[i] > 0.0:
  174 +
  175 + # Append value to array
  176 + lim_e.append(c_energy[i])
  177 + lim_flux.append(c_ulim[i])
  178 + lim_edflux.append(0.9*c_ulim[i])
  179 + lim_euflux.append(0.0)
  180 + else:
  181 + # If there is no upper limit available, this point is skipped
  182 + continue
  183 +
  184 + else:
  185 +
  186 + # Append significant data point to array
  187 + energy.append(c_energy[i])
  188 + flux.append(c_flux[i])
  189 + ed_flux.append(c_flux_err[i])
  190 + eu_flux.append(c_flux_err[i])
  191 +
  192 + # Return plot values
  193 + return energy,flux,ed_flux,eu_flux, lim_e, lim_flux, lim_edflux, lim_euflux
  194 +
  195 + def get_line(self, xmin = 0.1, xmax = 100.):
  196 +
  197 + # Initialise energy binning
  198 + ebounds = gammalib.GEbounds(100, gammalib.GEnergy(xmin,"TeV"), gammalib.GEnergy(xmax,"TeV"), True)
  199 +
  200 + # Initialise return array
  201 + x = []
  202 + y = []
  203 +
  204 + # Initialise dummy time
  205 + t = gammalib.GTime()
  206 +
  207 + # Loop over energy bins
  208 + for i in range(ebounds.size()):
  209 +
  210 + # Append values to array
  211 + energy = ebounds.elogmean(i)
  212 + x.append(energy.TeV())
  213 + y.append(energy.MeV() * energy.MeV() * self.m_source.spectral().eval(energy, t) * MeV2erg)
  214 +
  215 + # Return plot values
  216 + return x,y
  217 +
  218 + def get_butterfly(self):
  219 +
  220 + # Get ascii file
  221 + csv = gammalib.GCsv(self.m_bfly)
  222 +
  223 + # Initialise return arrays
  224 + x = []
  225 + y = []
  226 +
  227 + # Loop over rows
  228 + for i in range(csv.nrows()):
  229 +
  230 + #store energy
  231 + energy = csv.real(i,0) # MeV
  232 +
  233 + # store upper limit flux
  234 + flux = csv.real(i,1) + csv.real(i,2) / 2.0 #ph/cm2/s/MeV
  235 +
  236 + # append to array
  237 + x.append(energy * 1e-6) # TeV
  238 + y.append(flux * energy * energy * MeV2erg) # erg/cm2/s
  239 +
  240 + # Loop again over file from the other side
  241 + for i in range(csv.nrows()):
  242 +
  243 + # Create index counting backwards
  244 + k = csv.nrows()-1-i
  245 +
  246 + # Get values
  247 + energy = csv.real(k,0) # MeV
  248 + flux = csv.real(k,1)-csv.real(k,2) / 2.0 #ph/cm2/s/MeV
  249 +
  250 + # Append values
  251 + x.append(energy * 1e-6) # TeV
  252 + y.append(flux * energy * energy * MeV2erg) # erg/cm2/s
  253 +
  254 + # Return plot values
  255 + return x,y
  256 +
  257 +
  258 +
  259 + def execute(self):
  260 + """
  261 + Execute the script.
  262 + """
  263 + # Run the script
  264 + self.run()
  265 +
  266 + # Save image file
  267 + if not (self.m_outfile == "" or self.m_outfile == "NONE"):
  268 + plt.savefig(self.outfile)
  269 +
  270 + # Return
  271 + return
  272 +
  273 + def run(self):
  274 + """
  275 + Run the script.
  276 + """
  277 + # Switch screen logging on in debug mode
  278 + if self.logDebug():
  279 + self.log.cout(True)
  280 +
  281 + # Get parameters
  282 + self.get_parameters()
  283 +
  284 + # Write input parameters into logger
  285 + if self.logTerse():
  286 + self.log_parameters()
  287 + self.log("\n")
  288 +
  289 + plt.figure("Spectrum")
  290 + plt.loglog()
  291 + plt.xlabel("Energy [TeV]")
  292 + plt.ylabel(r'$\nu$F$\nu$ [erg cm$^{-2}$ s$^{-1}$]')
  293 +
  294 + # Write header
  295 + if self.logTerse():
  296 + self.log("\n")
  297 + self.log.header1("Get Butterfly")
  298 +
  299 + if not (self.m_bfly == "NONE" or self.m_bfly == ""):
  300 + xbf, ybf = self.get_butterfly()
  301 + plt.fill(xbf,ybf,color='green',alpha=0.5)
  302 +
  303 + # Write header
  304 + if self.logTerse():
  305 + self.log("\n")
  306 + self.log.header1("Get Spectrum Line")
  307 +
  308 + xline,yline = self.get_line(self.m_xmin, self.m_xmax)
  309 + plt.plot(xline,yline, color='black',lw=2.0)
  310 +
  311 + # Write header
  312 + if self.logTerse():
  313 + self.log("\n")
  314 + self.log.header1("Get Spectrum Points")
  315 +
  316 + if not (self.m_inspec == "NONE" or self.m_inspec == ""):
  317 + e,f,edf,euf,le,lf,ledf,leuf = self.get_points()
  318 + plt.errorbar(x=e,y=f,yerr=[edf,euf],fmt='bo',lw=2.0)
  319 + if len(le):
  320 + plt.errorbar(x=le,y=lf,yerr=[ledf,leuf],fmt='bo',lw=2.0, capsize=2.0,lowlims=True)
  321 +
  322 + if self.m_show:
  323 + plt.show()
  324 +
  325 +
  326 + # Return
  327 + return
  328 +
  329 +
  330 +# ======================== #
  331 +# Main routine entry point #
  332 +# ======================== #
  333 +if __name__ == '__main__':
  334 + """
  335 + Generates spectral points.
  336 + """
  337 + # Create instance of application
  338 + app = csplotspec(sys.argv)
  339 +
  340 + # Open logfile
  341 + app.logFileOpen()
  342 +
  343 + # Execute application
  344 + app.execute()
  345 +
0 346 \ No newline at end of file
... ...
scripts/csspec.py
  View file @1e6cba7
1 1 #! /usr/bin/env python
2 2 # ==========================================================================
3   -# Script for spectral data points generation.
  3 +# spectral points generation script.
4 4 #
5   -# Copyright (C) 2014 Michael Mayer
  5 +# Copyright (C) 2014-2015 Michael Mayer
6 6 #
7 7 # This program is free software: you can redistribute it and/or modify
8 8 # it under the terms of the GNU General Public License as published by
... ... @@ -18,33 +18,28 @@
18 18 # along with this program. If not, see <http://www.gnu.org/licenses/>.
19 19 #
20 20 # ==========================================================================
21   -import ctools
22 21 import gammalib
  22 +import ctools
23 23 import sys
24 24  
25   -import os
26   -
27   -
28 25 # ============== #
29 26 # cstsdist class #
30 27 # ============== #
31   -class csspec(ctools.ctool):
  28 +class csspec(ctools.cscript):
32 29 """
33   - This class implements the creation of spectral points. It derives
34   - from the ctools::ctool class which provides support for parameter
35   - files, command line arguments, and logging. In that way the Python
36   - script behaves just as a regular ctool.
  30 + This class implements the creation of spectral points. It derives from
  31 + the ctools.cscript class which provides support for parameter files,
  32 + command line arguments, and logging. In that way the Python script
  33 + behaves just as a regular ctool.
37 34 """
38 35 def __init__(self, *argv):
39 36 """
40 37 Constructor.
41 38 """
42   - self.name = "csspec"
43   - self.version = "0.1.0"
44   - #super(ctools.ctool).__init__(self.name,self.version,argv)
  39 +
45 40 # Set name
46   -# self.name = "csspectrum"
47   -# self.version = "0.1.0"
  41 + self.name = "csspec"
  42 + self.version = "1.0.0"
48 43  
49 44 # Initialise some members
50 45 self.obs = None
... ... @@ -55,18 +50,17 @@ class csspec(ctools.ctool):
55 50 argv = argv[1:]
56 51 else:
57 52 self.obs = gammalib.GObservations()
58   - self.obs.clear()
59   -
60   - self.outfile = ""
  53 + self.obs.clear()
  54 + self.m_outfile = ""
61 55  
62 56 # Make sure that parfile exists
63 57 file = self.parfile()
64 58  
65 59 # Initialise application
66 60 if len(argv) == 0:
67   - gammalib.GApplication.__init__(self, self.name, self.version)
  61 + ctools.cscript.__init__(self, self.name, self.version)
68 62 elif len(argv) ==1:
69   - gammalib.GApplication.__init__(self, self.name, self.version, *argv)
  63 + ctools.cscript.__init__(self, self.name, self.version, *argv)
70 64 else:
71 65 raise TypeError("Invalid number of arguments given.")
72 66  
... ... @@ -104,17 +98,24 @@ class csspec(ctools.ctool):
104 98  
105 99 # Create default parfile
106 100 pars = gammalib.GApplicationPars()
107   - pars.append(gammalib.GApplicationPar("infile","f","a","obs.xml","","","Observation definition"))
108   - pars.append(gammalib.GApplicationPar("srcmdl","f","a","crab.xml","","","model file name"))
  101 + pars.append(gammalib.GApplicationPar("inobs","f","a","events.fits","","","Event list, counts cube, or observation definition file"))
  102 + pars.append(gammalib.GApplicationPar("inmodel","f","a","$CTOOLS/share/models/crab.xml","","","Source model"))
109 103 pars.append(gammalib.GApplicationPar("srcname","s","a","Crab","","","Source name"))
110   - pars.append(gammalib.GApplicationPar("outfile","f","h","spectrum.fits","","","Output file name"))
111   - pars.append(gammalib.GApplicationPar("caldb","s","a","$GAMMALIB/share/caldb/cta","","","Calibration database"))
  104 + pars.append(gammalib.GApplicationPar("outfile","f","a","spectrum.fits","","","Output file name"))
  105 + pars.append(gammalib.GApplicationPar("expcube","f","a","NONE","","","Exposure cube file (only needed for stacked analysis)"))
  106 + pars.append(gammalib.GApplicationPar("psfcube","f","a","NONE","","","PSF cube file (only needed for stacked analysis)"))
  107 + pars.append(gammalib.GApplicationPar("caldb","s","a","dummy","","","Calibration database"))
112 108 pars.append(gammalib.GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function"))
113   - pars.append(gammalib.GApplicationPar("emin","r","a","0.1","0.0","","Lower energy limit (TeV)"))
114   - pars.append(gammalib.GApplicationPar("emax","r","a","100.0","0.0","","Upper energy limit (TeV)"))
  109 + pars.append(gammalib.GApplicationPar("emin","r","h","0.1","","","Lower energy limit (TeV)"))
  110 + pars.append(gammalib.GApplicationPar("emax","r","h","100.0","","","Upper energy limit (TeV)"))
115 111 pars.append(gammalib.GApplicationPar("enumbins","i","a","20","","","Number of energy bins"))
116 112 pars.append(gammalib.GApplicationPar("ebinalg","s","h","LOG","LIN|LOG|FILE","","Binning algorithm"))
  113 + pars.append(gammalib.GApplicationPar("calc_ts","b","h","yes","yes|no","","Compute TS value in each bin"))
  114 + pars.append(gammalib.GApplicationPar("calc_ulim","b","h","yes","yes|no","","Compute upper limit in each bin"))
  115 + pars.append(gammalib.GApplicationPar("fix_srcs","b","h","yes","yes|no","","Fix other skymodel parameters"))
  116 + pars.append(gammalib.GApplicationPar("fix_bkg","b","h","no","yes|no","","Fix background parameters"))
117 117 pars.append_standard()
  118 + pars.append(gammalib.GApplicationPar("logfile","f","h","csspec.log","","","Log filename"))
118 119 pars.save(parfile)
119 120  
120 121 # Return
... ... @@ -125,59 +126,36 @@ class csspec(ctools.ctool):
125 126 Get parameters from parfile and setup the observation.
126 127 """
127 128 # Get parameters
128   -
129 129 # Set observation if not done before
130   - if self.obs.size() == 0:
131   -
132   - obsfile = self["infile"].filename()
133   - try:
134   -
135   - self.obs = gammalib.GObservations(obsfile)
136   - except:
137   - self.obs.clear()
138   - observation = gammalib.GCTAObservation()
139   - observation.load(obsfile)
140   -
141   - self.m_irf = self["irf"].string()
142   - self.m_caldb = self["caldb"].string()
143   -
144   - caldb = gammalib.GCaldb()
145   - if os.path.isdir(self.m_caldb):
146   - caldb.rootdir(self.m_caldb)
147   - else:
148   - caldb.open("cta",self. m_caldb)
  130 + if self.obs == None or self.obs.size() == 0:
  131 + self.require_inobs("csspec::get_parameters()")
  132 + self.obs = self.get_observations()
149 133  
150   - observation.response(self.m_irf, caldb);
151   - self.obs.append(observation)
152   -
153   - # Check for models in the container
154   - # read models from file if there were none
155   - # Note that the models should be optimised for
156   - # meaningful residual map
  134 + # Set models if we have none
157 135 if self.obs.models().size() == 0:
158   - self.models = gammalib.GModels(self["srcmdl"].filename())
159   -
160   - self.m_srcname = self["srcname"].string()
161   - self.m_outfile = self["outfile"].filename()
162   -
163   -
164   - # not sure how to call get_ebounds as private member of ctool
165   - #self.ebounds = ctools.ctool.get_ebounds()
166   - self.m_emin = self["emin"].real()
167   - self.m_emax = self["emax"].real()
168   - self.m_enumbins = self["enumbins"].integer()
169   - #self.m_ebinalg = self["ebinalg"].string()
  136 + self.obs.models(self["inmodel"].filename())
  137 +
  138 + # Get source name
  139 + self.m_srcname = self["srcname"].string()
170 140  
171   - self.ebounds = gammalib.GEbounds(self.m_enumbins,gammalib.GEnergy(self.m_emin,"TeV"),gammalib.GEnergy(self.m_emax,"TeV"),True)
  141 + # Get ebounds
  142 + self.m_ebounds = self.create_ebounds()
172 143  
  144 + # Read other parameters
  145 + self.m_outfile = self["outfile"].filename()
  146 +
  147 + # Get other parameeters
  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()
173 152  
174 153 # Set some fixed parameters
175   - self.m_log = False # Logging in client tools
176   - self.m_debug = False # Debugging in client tools
  154 + self.m_log = False # Logging in client tools
  155 + self.m_chatter = self["chatter"].integer()
177 156 self.m_clobber = self["clobber"].boolean()
178   -
179   - self.parname =self.models[self.m_srcname].spectral()[0].name()
180   -
  157 + self.m_debug = self["debug"].boolean()
  158 +
181 159 # Return
182 160 return
183 161  
... ... @@ -199,49 +177,21 @@ class csspec(ctools.ctool):
199 177 self.run()
200 178  
201 179 # Save residual map
202   - self.save(self.m_outfile, self.m_clobber)
  180 + self.fits.saveto(self.m_outfile, self.m_clobber)
203 181  
204 182 # Return
205 183 return
206 184  
207   - def analyse_bin(self, emin,emax):
208   - #print "1",self.models.size()
209   - select = ctools.ctselect(self.obs)
210   - select["emin"].real(emin)
211   - select["emax"].real(emax)
212   - select["ra"].real(-1.0)
213   - select["dec"].real(-1.0)
214   - select["rad"].real(-1.0)
215   - select["tmin"].real(0.0)
216   - select["tmax"].real(0.0)
217   - select.run()
218   - #print "2",self.models.size()
219   - like = ctools.ctlike(select.obs())
220   - like.run()
221   - optmodels = like.obs().models()
222   - #print "3",self.models.size()
223   - return optmodels.copy()
224   -
225 185 def run(self):
226 186 """
227 187 Run the script.
228 188 """
229   -
230 189 # Switch screen logging on in debug mode
231 190 if self.logDebug():
232 191 self.log.cout(True)
233 192  
234 193 # Get parameters
235 194 self.get_parameters()
236   -
237   -
238   - for model in self.models:
239   - for par in model:
240   - par.fix()
241   -
242   - self.models[self.m_srcname][self.parname].free()
243   - self.models[self.m_srcname].tscalc(True)
244   - self.obs.models(self.models)
245 195  
246 196 # Write input parameters into logger
247 197 if self.logTerse():
... ... @@ -258,46 +208,151 @@ class csspec(ctools.ctool):
258 208 # Write header
259 209 if self.logTerse():
260 210 self.log("\n")
261   - self.log.header1("Generate binned map")
  211 + self.log.header1("Adjust model parameters")
  212 +
  213 + for model in self.obs.models():
  214 + model.tscalc(False)
  215 + if model.name() == self.m_srcname:
  216 + for par in model:
  217 + par.fix()
  218 + model.spectral()[0].free()
  219 + if self.m_calc_ts:
  220 + model.tscalc(True)
  221 +
  222 + elif self.m_fix_bkg and not model.classname() == "GModelData":
  223 + for par in model:
  224 + par.fix()
  225 +
  226 + elif self.m_fix_srcs and model.classname() == "GModelSky":
  227 + for par in model:
  228 + par.fix()
  229 +
  230 + # Write header
  231 + if self.logTerse():
  232 + self.log("\n")
  233 + self.log.header1("Generate spectrum")
  234 +
  235 + # Initialise FITS Table with extension "SPECTRUM"
  236 + table = gammalib.GFitsBinTable(self.m_ebounds.size())
  237 + table.extname("SPECTRUM")
262 238  
263   - binned=False
  239 + # Add Header for compatibility with gammalib.GMWLSpectrum
  240 + table.card("Instrument","CTA","Name of Instrument")
  241 + table.card("Telescope","CTA","Name of Telescope")
264 242  
  243 + # Create FITS table columns
  244 + energy = gammalib.GFitsTableDoubleCol("energy",self.m_ebounds.size())
  245 + energy.unit("TeV")
  246 + energy_err = gammalib.GFitsTableDoubleCol("energy_err",self.m_ebounds.size())
  247 + energy_err.unit("TeV")
  248 + flux = gammalib.GFitsTableDoubleCol("flux",self.m_ebounds.size())
  249 + flux.unit("erg/cm2/s")
  250 + flux_err = gammalib.GFitsTableDoubleCol("flux_err",self.m_ebounds.size())
  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())
  254 + ulim_values.unit("erg/cm2/s")
265 255  
266   - for i in range(self.ebounds.size()):
  256 + # MeV2erg
  257 + MeV2erg = 1.60217e-6
267 258  
268   - emin = self.ebounds.emin(i).TeV()
269   - emax = self.ebounds.emax(i).TeV()
270   - emean = self.ebounds.elogmean(i)
271   - optmodels = self.analyse_bin(emin,emax)
272   - #print optmodels
273   - #if optmodels.size()>0:
274   - optmodel = optmodels[self.m_srcname]
275   - ts = optmodel.ts()
  259 + for i in range(self.m_ebounds.size()):
276 260  
277   - flux = optmodel.spectral().eval(emean,gammalib.GTime())
278   - rel_error = optmodel[self.parname].error()/optmodel[self.parname].value()
  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)
  265 + elogmean2 = elogmean.MeV() * elogmean.MeV()
279 266  
280   - flux_error = flux * rel_error
  267 + # Select events
  268 + select = ctools.ctselect(self.obs)
  269 + select["emin"].real(emin.TeV())
  270 + select["emax"].real(emax.TeV())
  271 + select["tmin"].real(0.0)
  272 + select["tmax"].real(0.0)
  273 + select["rad"].real(0.0)
  274 + select["ra"].real(0.0)
  275 + select["dec"].real(-1.0)
  276 + select.run()
281 277  
282   - print emean,flux,flux_error,ts
283   -
284   - if binned:
285   - print "not implemented"
  278 + # likelihood
  279 + like = ctools.ctlike(select.obs())
  280 + like.run()
  281 +
  282 + # Get results
  283 + fitted_models = like.obs().models()
  284 + source = fitted_models[self.m_srcname]
  285 +
  286 + # Calculate Upper Limit
  287 + ulimit_value = -1.0
  288 + #if self.m_calc_ulimit:
  289 + # ulimit = ctools.ctulimit(like.obs())
  290 + # ulimit["eref"].real(elogmean.TeV())
  291 + # ulimit.run()
  292 + # ulimit_value.ulimit.diff_ulimit()
  293 +
  294 + # Get TS value
  295 + TS = -1.0
  296 + if self.m_calc_ts:
  297 + TS = source.ts()
  298 +
  299 + # Get differential flux
  300 + fitted_flux = source.spectral().eval(elogmean,gammalib.GTime())
  301 +
  302 + # Compute flux error
  303 + parvalue = source.spectral()[0].value()
  304 + rel_error = source.spectral()[0].error()/parvalue
  305 + e_flux = fitted_flux*rel_error
  306 +
  307 + # Set values for storage
  308 + TSvalues[i] = TS
  309 +
  310 + # Store energy as TeV
  311 + energy[i] = elogmean.TeV()
  312 +
  313 + # Store energy error
  314 + energy_err[i] = (emax - emean).TeV()
  315 +
  316 + # Convert fluxes to nuFnu
  317 + flux[i] = fitted_flux * elogmean2 * MeV2erg
  318 + flux_err[i] = e_flux * elogmean2 * MeV2erg
  319 + if ulimit_value > 0.0:
  320 + ulim_values[i] = ulimit_value * elogmean2 * MeV2erg
286 321  
287   - self.obs.models(self.models)
  322 +
  323 + # Log information
  324 + if self.logExplicit():
  325 + self.log("Bin "+str(i)+" ["+str(emin.TeV())+"-"+str(emax.TeV())+"] TeV: ")
  326 + self.log("Flux = "+str(flux[i]))
  327 + self.log(" +- "+str(flux_err[i])+" [erg/cm2/s]")
  328 + if self.m_calc_ts and TSvalues[i] > 0.0:
  329 + self.log(", TS = "+str(TS))
  330 + if self.m_calc_ulimit and ulim_values[i] > 0.0:
  331 + self.log(", UL = "+str(ulim_values[i])+" [erg/cm2/s]")
  332 + self.log("\n")
  333 +
  334 + # Append filles columns to fits table
  335 + table.append(energy)
  336 + table.append(energy_err)
  337 + table.append(flux)
  338 + table.append(flux_err)
  339 + table.append(TSvalues)
  340 + table.append(ulim_values)
288 341  
  342 + # Create the FITS file now
  343 + self.fits = gammalib.GFits()
  344 + self.fits.append(table)
  345 +
289 346 # Return
290 347 return
291 348  
292   - def save(self,outfile,clobber):
293   - return
294 349  
295 350 # ======================== #
296 351 # Main routine entry point #
297 352 # ======================== #
298 353 if __name__ == '__main__':
299 354 """
300   - Generates residual count map.
  355 + Generates spectral points.
301 356 """
302 357 # Create instance of application
303 358 app = csspec(sys.argv)
... ...
test/test_cscripts.sh.in
  View file @1e6cba7
... ... @@ -160,3 +160,30 @@ else
160 160 exit 1
161 161 fi
162 162 $ECHO " ok"
  163 +
  164 +#
  165 +# Test csspec
  166 +# =============
  167 +$ECHO -n "Test csspec: "
  168 +chmod +x @abs_top_srcdir@/scripts/csspec.py
  169 +@abs_top_srcdir@/scripts/csspec.py \
  170 + inobs="data/crab_events.fits" \
  171 + inmodel="data/crab.xml" \
  172 + caldb="irf" \
  173 + irf="cta_dummy_irf" \
  174 + enumbins=5 \
  175 + srcname="Crab" \
  176 + emin=0.1 \
  177 + emax=100.0 \
  178 + outfile="spectrum.fits" > csspec.log
  179 +$ECHO -n "."
  180 +if [ -s "spectrum.fits" ]
  181 +then
  182 + $ECHO -n "."
  183 +else
  184 + $ECHO " spectrum.fits file is not valid"
  185 + exit 1
  186 +fi
  187 +$ECHO " ok"
  188 +
  189 +
... ...