Pierrick Martin / ctools

Commit e446f62c5d1a674194987f0880017ab61963fb66

Authored by Jürgen Knödlseder
1 parent ab3bcbff

Refactor csspec script (#2231)

Inline Side-by-side
Showing 2 changed files with 356 additions and 257 deletions
cscripts/csobsselect.py
  View file @e446f62
@@ -283,7 +283,7 @@ class csobsselect(ctools.cscript): @@ -283,7 +283,7 @@ class csobsselect(ctools.cscript):
283 self._log_observations(gammalib.NORMAL, self._obs, 'Input observation') 283 self._log_observations(gammalib.NORMAL, self._obs, 'Input observation')
284 284
285 # Write header 285 # Write header
286 - self._log_header1(gammalib.NORMAL, 'Observation selection') 286 + self._log_header1(gammalib.TERSE, 'Observation selection')
287 287
288 # Loop over observations 288 # Loop over observations
289 for obs in self._obs: 289 for obs in self._obs:
cscripts/csspec.py
  View file @e446f62
@@ -28,42 +28,31 @@ import ctools @@ -28,42 +28,31 @@ import ctools
28 # ============ # 28 # ============ #
29 class csspec(ctools.cscript): 29 class csspec(ctools.cscript):
30 """ 30 """
31 - Generates a spectrum. 31 + Generates a spectrum
32 32
33 - This class implements the creation of spectral points. It derives from  
34 - the ctools.cscript class which provides support for parameter files,  
35 - command line arguments, and logging. In that way the Python script  
36 - behaves just as a regular ctool. 33 + This class implements the generation of a spectral energy distribution
  34 + from data.
37 """ 35 """
38 36
39 # Constructors and destructors 37 # Constructors and destructors
40 def __init__(self, *argv): 38 def __init__(self, *argv):
41 """ 39 """
42 - Constructor. 40 + Constructor
43 """ 41 """
44 # Set name 42 # Set name
45 self._name = 'csspec' 43 self._name = 'csspec'
46 self._version = ctools.__version__ 44 self._version = ctools.__version__
47 45
48 - # Initialise some members 46 + # Initialise data members
49 self._ebounds = gammalib.GEbounds() 47 self._ebounds = gammalib.GEbounds()
50 self._fits = None 48 self._fits = None
51 self._binned_mode = False 49 self._binned_mode = False
52 - self._srcname = ""  
53 - self._edisp = False  
54 - self._calc_ulimit = True  
55 - self._calc_ts = True  
56 - self._fix_bkg = False  
57 - self._fix_srcs = True  
58 - self._chatter = 2  
59 - self._clobber = True  
60 - self._debug = False  
61 -  
62 - # Initialise observation container from constructor arguments. 50 + self._onoff_mode = False
  51 +
  52 + # Initialise observation container from constructor arguments
63 self._obs, argv = self._set_input_obs(argv) 53 self._obs, argv = self._set_input_obs(argv)
64 54
65 - # Initialise application by calling the appropriate class  
66 - # constructor. 55 + # Initialise application by calling the appropriate class constructor
67 self._init_cscript(argv) 56 self._init_cscript(argv)
68 57
69 # Return 58 # Return
@@ -71,7 +60,7 @@ class csspec(ctools.cscript): @@ -71,7 +60,7 @@ class csspec(ctools.cscript):
71 60
72 def __del__(self): 61 def __del__(self):
73 """ 62 """
74 - Destructor. 63 + Destructor
75 """ 64 """
76 # Return 65 # Return
77 return 66 return
@@ -80,47 +69,44 @@ class csspec(ctools.cscript): @@ -80,47 +69,44 @@ class csspec(ctools.cscript):
80 # Private methods 69 # Private methods
81 def _get_parameters(self): 70 def _get_parameters(self):
82 """ 71 """
83 - Get parameters from parfile and setup the observation. 72 + Get parameters from parfile and setup the observation
84 """ 73 """
85 # Set observation if not done before 74 # Set observation if not done before
86 if self._obs == None or self._obs.size() == 0: 75 if self._obs == None or self._obs.size() == 0:
87 - self._require_inobs("csspec::get_parameters()") 76 + self._require_inobs('csspec::get_parameters()')
88 self._obs = self._get_observations() 77 self._obs = self._get_observations()
89 -  
90 - # Check if we have one binned cta observation, i.e. if we are in binned mode  
91 - self._binned_mode = False 78 +
  79 + # Check if we have one binned CTA observation, i.e. if we are in
  80 + # binned mode
  81 + self._binned_mode = False
  82 + self._onoff_mode = False
92 if self._obs.size() == 1: 83 if self._obs.size() == 1:
93 - if self._obs[0].classname() == "GCTAObservation":  
94 - if self._obs[0].eventtype() == "CountsCube": 84 + if self._obs[0].classname() == 'GCTAObservation':
  85 + if self._obs[0].eventtype() == 'CountsCube':
95 self._binned_mode = True 86 self._binned_mode = True
  87 + elif self._obs[0].classname() == 'GCTAOnOffObservation':
  88 + self._onoff_mode = True
96 89
97 # Set models if we have none 90 # Set models if we have none
98 if self._obs.models().size() == 0: 91 if self._obs.models().size() == 0:
99 - self._obs.models(self["inmodel"].filename()) 92 + self._obs.models(self['inmodel'].filename())
100 93
101 - # Get source name  
102 - self._srcname = self["srcname"].string() 94 + # Query source name
  95 + self['srcname'].string()
103 96
104 # Set ebounds 97 # Set ebounds
105 self._set_ebounds() 98 self._set_ebounds()
106 -  
107 - # Get edisp flag  
108 - self._edisp = self["edisp"].boolean()  
109 99
110 - # Get other parameeters  
111 - self._calc_ulimit = self["calc_ulim"].boolean()  
112 - self._calc_ts = self["calc_ts"].boolean()  
113 - self._fix_bkg = self["fix_bkg"].boolean()  
114 - self._fix_srcs = self["fix_srcs"].boolean()  
115 -  
116 - # Set some fixed parameters  
117 - self._chatter = self["chatter"].integer()  
118 - self._clobber = self["clobber"].boolean()  
119 - self._debug = self["debug"].boolean() 100 + # Query other parameeters
  101 + self['edisp'].boolean()
  102 + self['calc_ulim'].boolean()
  103 + self['calc_ts'].boolean()
  104 + self['fix_bkg'].boolean()
  105 + self['fix_srcs'].boolean()
120 106
121 # Read ahead output parameters 107 # Read ahead output parameters
122 if self._read_ahead(): 108 if self._read_ahead():
123 - self["outfile"].filename() 109 + self['outfile'].filename()
124 110
125 # Write input parameters into logger 111 # Write input parameters into logger
126 self._log_parameters(gammalib.TERSE) 112 self._log_parameters(gammalib.TERSE)
@@ -130,10 +116,11 @@ class csspec(ctools.cscript): @@ -130,10 +116,11 @@ class csspec(ctools.cscript):
130 116
131 def _set_ebounds(self): 117 def _set_ebounds(self):
132 """ 118 """
133 - Set energy boundaries. 119 + Set energy boundaries
134 """ 120 """
135 - # Get ebounds  
136 - if self._binned_mode: 121 + # If we are in binned or On/Off mode then align the energy boundaries
  122 + # with the cube of RMF spectrum
  123 + if self._binned_mode or self._onoff_mode:
137 124
138 # Initialise energy boundaries for spectrum 125 # Initialise energy boundaries for spectrum
139 self._ebounds = gammalib.GEbounds() 126 self._ebounds = gammalib.GEbounds()
@@ -142,9 +129,12 @@ class csspec(ctools.cscript): @@ -142,9 +129,12 @@ class csspec(ctools.cscript):
142 ebounds = self._create_ebounds() 129 ebounds = self._create_ebounds()
143 130
144 # Extract cube ebounds 131 # Extract cube ebounds
145 - cube_ebounds = self._obs[0].events().ebounds() 132 + if self._binned_mode:
  133 + cube_ebounds = self._obs[0].events().ebounds()
  134 + else:
  135 + cube_ebounds = self._obs[0].rmf().etrue()
146 136
147 - # Loop over user energy boundaries and collect all layers 137 + # Loop over user energy boundaries and collect all energy bins
148 # that overlap 138 # that overlap
149 for i in range(ebounds.size()): 139 for i in range(ebounds.size()):
150 140
@@ -193,270 +183,345 @@ class csspec(ctools.cscript): @@ -193,270 +183,345 @@ class csspec(ctools.cscript):
193 # Return 183 # Return
194 return 184 return
195 185
196 -  
197 - # Public methods  
198 - def run(self): 186 + def _log_spectral_binning(self):
199 """ 187 """
200 - Run the script. 188 + Log spectral binning
201 """ 189 """
202 - # Switch screen logging on in debug mode  
203 - if self._logDebug():  
204 - self._log.cout(True)  
205 -  
206 - # Get parameters  
207 - self._get_parameters() 190 + # Write header
  191 + self._log_header1(gammalib.TERSE, 'Spectral binning')
208 192
209 - # Write spectral binning into header  
210 - self._log("\n")  
211 - self._log.header1("Spectral binning") 193 + # Log counts cube energy range for binned mode
212 if self._binned_mode: 194 if self._binned_mode:
213 cube_ebounds = self._obs[0].events().ebounds() 195 cube_ebounds = self._obs[0].events().ebounds()
214 value = '%s - %s' % (str(cube_ebounds.emin()), 196 value = '%s - %s' % (str(cube_ebounds.emin()),
215 str(cube_ebounds.emax())) 197 str(cube_ebounds.emax()))
216 self._log_value(gammalib.TERSE, 'Counts cube energy range', value) 198 self._log_value(gammalib.TERSE, 'Counts cube energy range', value)
  199 +
  200 + # Log energy bins
217 for i in range(self._ebounds.size()): 201 for i in range(self._ebounds.size()):
218 value = '%s - %s' % (str(self._ebounds.emin(i)), 202 value = '%s - %s' % (str(self._ebounds.emin(i)),
219 str(self._ebounds.emax(i))) 203 str(self._ebounds.emax(i)))
220 self._log_value(gammalib.TERSE, 'Bin %d' % (i+1), value) 204 self._log_value(gammalib.TERSE, 'Bin %d' % (i+1), value)
221 205
222 - # Write observation into logger  
223 - if self._logTerse():  
224 - self._log("\n")  
225 - self._log.header1("Observation")  
226 - self._log(str(self._obs))  
227 - self._log("\n") 206 + # Return
  207 + return
228 208
  209 + def _adjust_model_pars(self):
  210 + """
  211 + Adjust model parameters
  212 + """
229 # Write header 213 # Write header
230 - if self._logTerse():  
231 - self._log("\n")  
232 - self._log.header1("Adjust model parameters") 214 + self._log_header1(gammalib.TERSE, 'Adjust model parameters')
233 215
234 # Adjust model parameters dependent on input user parameters 216 # Adjust model parameters dependent on input user parameters
235 for model in self._obs.models(): 217 for model in self._obs.models():
236 218
237 - # Set TS flag for all models to false.  
238 - # Source of interest will be set to true later 219 + # Initialise TS flag for all models to false
239 model.tscalc(False) 220 model.tscalc(False)
240 221
241 # Log model name 222 # Log model name
242 - if self._logExplicit():  
243 - self._log.header3(model.name()) 223 + self._log_header3(gammalib.EXPLICIT, model.name())
244 224
245 # Deal with the source of interest 225 # Deal with the source of interest
246 - if model.name() == self._srcname: 226 + if model.name() == self['srcname'].string():
  227 +
  228 + # Fix all model parameters
247 for par in model: 229 for par in model:
248 - if par.is_free() and self._logExplicit():  
249 - self._log(" Fixing \""+par.name()+"\"\n") 230 + if par.is_free():
  231 + self._log_string(gammalib.EXPLICIT,
  232 + ' Fixing "'+par.name()+'"')
250 par.fix() 233 par.fix()
  234 +
  235 + # Free the normalisation parameter which is assumed to be
  236 + # the first spectral parameter
251 normpar = model.spectral()[0] 237 normpar = model.spectral()[0]
252 - if normpar.is_fixed() and self._logExplicit():  
253 - self._log(" Freeing \""+normpar.name()+"\"\n") 238 + if normpar.is_fixed():
  239 + self._log_string(gammalib.EXPLICIT,
  240 + ' Freeing "'+normpar.name()+'"')
254 normpar.free() 241 normpar.free()
255 - if self._calc_ts: 242 +
  243 + # Optionally compute Test Statistic value
  244 + if self['calc_ts'].boolean():
256 model.tscalc(True) 245 model.tscalc(True)
257 246
258 - elif self._fix_bkg and not model.classname() == "GModelSky": 247 + # Deal with background models
  248 + elif self['fix_bkg'].boolean() and \
  249 + not model.classname() == 'GModelSky':
259 for par in model: 250 for par in model:
260 - if par.is_free() and self._logExplicit():  
261 - self._log(" Fixing \""+par.name()+"\"\n") 251 + if par.is_free():
  252 + self._log_string(gammalib.EXPLICIT,
  253 + ' Fixing "'+par.name()+'"')
262 par.fix() 254 par.fix()
263 255
264 - elif self._fix_srcs and model.classname() == "GModelSky": 256 + # Deal with source models
  257 + elif self['fix_srcs'].boolean() and \
  258 + model.classname() == 'GModelSky':
265 for par in model: 259 for par in model:
266 - if par.is_free() and self._logExplicit():  
267 - self._log(" Fixing \""+par.name()+"\"\n") 260 + if par.is_free():
  261 + self._log_string(gammalib.EXPLICIT,
  262 + ' Fixing "'+par.name()+'"')
268 par.fix() 263 par.fix()
269 264
270 - # Write header  
271 - self._log_header1(gammalib.TERSE, 'Generate spectrum')  
272 - if self._logTerse():  
273 - self._log(str(self._ebounds)) 265 + # Return
  266 + return
274 267
275 - # Initialise FITS Table with extension "SPECTRUM"  
276 - table = gammalib.GFitsBinTable(self._ebounds.size())  
277 - table.extname('SPECTRUM') 268 + def _select_obs(self, emin, emax):
  269 + """
  270 + Select observations for energy interval
278 271
279 - # Add Header for compatibility with gammalib.GMWLSpectrum  
280 - table.card('INSTRUME', 'CTA', 'Name of Instrument')  
281 - table.card('TELESCOP', 'CTA', 'Name of Telescope') 272 + Parameters
  273 + ----------
  274 + emin : `~gammalib.GEnergy()`
  275 + Minimum energy
  276 + emax : `~gammalib.GEnergy()`
  277 + Maximum energy
  278 +
  279 + Returns
  280 + -------
  281 + obs : `~gammalib.GObservations()`
  282 + Observation container
  283 + """
  284 + # Use ctcubemask for binned analysis
  285 + if self._binned_mode:
282 286
283 - # Create FITS table columns  
284 - nrows = self._ebounds.size()  
285 - energy = gammalib.GFitsTableDoubleCol("Energy", nrows)  
286 - energy_low = gammalib.GFitsTableDoubleCol("ed_Energy", nrows)  
287 - energy_high = gammalib.GFitsTableDoubleCol("eu_Energy", nrows)  
288 - flux = gammalib.GFitsTableDoubleCol("Flux", nrows)  
289 - flux_err = gammalib.GFitsTableDoubleCol("e_Flux", nrows)  
290 - TSvalues = gammalib.GFitsTableDoubleCol("TS", nrows)  
291 - ulim_values = gammalib.GFitsTableDoubleCol("UpperLimit", nrows)  
292 - Npred_values = gammalib.GFitsTableDoubleCol("Npred", nrows)  
293 - energy.unit("TeV")  
294 - energy_low.unit("TeV")  
295 - energy_high.unit("TeV")  
296 - flux.unit("erg/cm2/s")  
297 - flux_err.unit("erg/cm2/s")  
298 - ulim_values.unit("erg/cm2/s") 287 + # Write header
  288 + self._log_header3(gammalib.EXPLICIT, 'Filtering cube')
299 289
300 - # Loop over energy bins  
301 - for i in range(nrows): 290 + # Select layers
  291 + cubemask = ctools.ctcubemask(self._obs)
  292 + cubemask['regfile'] = 'NONE'
  293 + cubemask['ra'] = 'UNDEFINED'
  294 + cubemask['dec'] = 'UNDEFINED'
  295 + cubemask['rad'] = 'UNDEFINED'
  296 + cubemask['emin'] = emin.TeV()
  297 + cubemask['emax'] = emax.TeV()
  298 + cubemask.run()
  299 +
  300 + # Set new binned observation
  301 + obs = cubemask.obs().copy()
302 302
303 - # Log information  
304 - if self._logExplicit():  
305 - self._log("\n")  
306 - self._log.header2("Energy bin "+str(i+1)) 303 + # Use ...
  304 + elif self._onoff_mode:
307 305
308 - # Get energy boundaries  
309 - emin = self._ebounds.emin(i)  
310 - emax = self._ebounds.emax(i)  
311 - elogmean = self._ebounds.elogmean(i)  
312 - elogmean2 = elogmean.MeV() * elogmean.MeV() 306 + # Write header
  307 + self._log_header3(gammalib.EXPLICIT, 'Filtering PHA, ARF and RMF')
313 308
314 - # Store energy as TeV  
315 - energy[i] = elogmean.TeV() 309 + # DUMMY
  310 + obs = self._obs.copy()
316 311
317 - # Store energy errors  
318 - energy_low[i] = (elogmean - emin).TeV()  
319 - energy_high[i] = (emax - elogmean).TeV()  
320 312
321 - # Use ctselect for unbinned analysis  
322 - if not self._binned_mode:  
323 -  
324 - # Log information  
325 - if self._logExplicit():  
326 - self._log.header3("Selecting events") 313 + # Use ctselect for unbinned analysis
  314 + else:
  315 +
  316 + # Write header
  317 + self._log_header3(gammalib.EXPLICIT, 'Selecting events')
327 318
328 - # Select events  
329 - select = ctools.ctselect(self._obs)  
330 - select["emin"] = emin.TeV()  
331 - select["emax"] = emax.TeV()  
332 - select["tmin"] = "UNDEFINED"  
333 - select["tmax"] = "UNDEFINED"  
334 - select["rad"] = "UNDEFINED"  
335 - select["ra"] = "UNDEFINED"  
336 - select["dec"] = "UNDEFINED"  
337 - select.run() 319 + # Select events
  320 + select = ctools.ctselect(self._obs)
  321 + select['emin'] = emin.TeV()
  322 + select['emax'] = emax.TeV()
  323 + select['tmin'] = 'UNDEFINED'
  324 + select['tmax'] = 'UNDEFINED'
  325 + select['rad'] = 'UNDEFINED'
  326 + select['ra'] = 'UNDEFINED'
  327 + select['dec'] = 'UNDEFINED'
  328 + select.run()
338 329
339 - # Retrieve observation  
340 - obs = select.obs() 330 + # Retrieve observation
  331 + obs = select.obs().copy()
341 332
342 - # Use ctcubemask for binned analysis  
343 - else: 333 + # Return observation container
  334 + return obs
  335 +
  336 + def _fit_energy_bin(self, i):
  337 + """
  338 + Fit data for one energy bin
  339 +
  340 + Parameters
  341 + ----------
  342 + i : int
  343 + Energy bin index
344 344
345 - # Header  
346 - if self._logExplicit():  
347 - self._log.header3("Filtering cube")  
348 -  
349 - # Select layers  
350 - cubemask = ctools.ctcubemask(self._obs)  
351 - cubemask["regfile"] = "NONE"  
352 - cubemask["ra"] = "UNDEFINED"  
353 - cubemask["dec"] = "UNDEFINED"  
354 - cubemask["rad"] = "UNDEFINED"  
355 - cubemask["emin"] = emin.TeV()  
356 - cubemask["emax"] = emax.TeV()  
357 - cubemask.run()  
358 -  
359 - # Set new binned observation  
360 - obs = cubemask.obs()  
361 -  
362 - # Header  
363 - if self._logExplicit():  
364 - self._log.header3("Performing fit")  
365 -  
366 - # Likelihood  
367 - like = ctools.ctlike(obs)  
368 - like["edisp"] = self._edisp  
369 - like.run()  
370 -  
371 - # Skip bin if no event was present  
372 - if like.obs().logL() == 0.0:  
373 -  
374 - # Log information  
375 - if self._logExplicit():  
376 - self._log("No event in this bin. ")  
377 - self._log("Likelihood is zero. ")  
378 - self._log("Bin is skipped.")  
379 -  
380 - # Set all values to 0  
381 - flux[i] = 0.0  
382 - flux_err[i] = 0.0  
383 - TSvalues[i] = 0.0  
384 - ulim_values[i] = 0.0  
385 - Npred_values[i] = 0.0  
386 - continue 345 + Returns
  346 + -------
  347 + result : dict
  348 + Dictionary with fit results
  349 + """
  350 + # Get energy boundaries
  351 + emin = self._ebounds.emin(i)
  352 + emax = self._ebounds.emax(i)
  353 + elogmean = self._ebounds.elogmean(i)
  354 +
  355 + # Select observations for energy bin
  356 + obs = self._select_obs(emin, emax)
  357 +
  358 + # Initialise dictionary
  359 + result = {'energy': elogmean.TeV(),
  360 + 'energy_low': (elogmean - emin).TeV(),
  361 + 'energy_high': (emax - elogmean).TeV(),
  362 + 'flux': 0.0,
  363 + 'flux_err': 0.0,
  364 + 'TS': 0.0,
  365 + 'ulimit': 0.0,
  366 + 'Npred': 0.0}
  367 +
  368 + # Write header for fitting
  369 + self._log_header3(gammalib.EXPLICIT, 'Performing fit')
  370 +
  371 + # Perform maximum likelihood fit
  372 + like = ctools.ctlike(obs)
  373 + like['edisp'] = self['edisp'].boolean()
  374 + like.run()
  375 +
  376 + # Continue only if log-likelihood is non-zero
  377 + if like.obs().logL() != 0.0:
387 378
388 # Get results 379 # Get results
389 fitted_models = like.obs().models() 380 fitted_models = like.obs().models()
390 - source = fitted_models[self._srcname] 381 + source = fitted_models[self['srcname'].string()]
  382 +
  383 + # Extract Test Statistic value
  384 + if self['calc_ts'].boolean():
  385 + result['TS'] = source.ts()
391 386
392 - # Calculate Upper Limit 387 + # Compute Npred value (only works for unbinned analysis)
  388 + if not self._binned_mode:
  389 + for observation in like.obs():
  390 + result['Npred'] += observation.npred(source)
  391 +
  392 + # Compute upper flux limit
393 ulimit_value = -1.0 393 ulimit_value = -1.0
394 - if self._calc_ulimit: 394 + if self['calc_ulim'].boolean():
395 395
396 # Logging information 396 # Logging information
397 - if self._logExplicit():  
398 - self._log.header3("Computing upper limit") 397 + self._log_header3(gammalib.EXPLICIT, 'Computing upper limit')
399 398
400 # Create upper limit object 399 # Create upper limit object
401 ulimit = ctools.ctulimit(like.obs()) 400 ulimit = ctools.ctulimit(like.obs())
402 - ulimit["srcname"] = self._srcname  
403 - ulimit["eref"] = elogmean.TeV() 401 + ulimit['srcname'] = self['srcname'].string()
  402 + ulimit['eref'] = elogmean.TeV()
404 403
405 # Try to run upper limit and catch exceptions 404 # Try to run upper limit and catch exceptions
406 try: 405 try:
407 ulimit.run() 406 ulimit.run()
408 ulimit_value = ulimit.diff_ulimit() 407 ulimit_value = ulimit.diff_ulimit()
409 except: 408 except:
410 - if self._logExplicit():  
411 - self._log("Upper limit calculation failed.") 409 + self._log_string(gammalib.EXPLICIT, 'Upper limit '
  410 + 'calculation failed.')
412 ulimit_value = -1.0 411 ulimit_value = -1.0
413 412
414 - # Get TS value  
415 - TS = -1.0  
416 - if self._calc_ts:  
417 - TS = source.ts() 413 + # Compute upper limit
  414 + if ulimit_value > 0.0:
  415 + result['ulimit'] = ulimit_value * elogmean.MeV() * \
  416 + elogmean.MeV() * gammalib.MeV2erg
418 417
419 - # Compute Npred value (only works for unbinned analysis)  
420 - Npred = 0.0  
421 - if not self._binned_mode:  
422 - for observation in like.obs():  
423 - Npred += observation.npred(source)  
424 -  
425 - # Get differential flux 418 + # Compute differential flux and flux error
426 fitted_flux = source.spectral().eval(elogmean) 419 fitted_flux = source.spectral().eval(elogmean)
427 -  
428 - # Compute flux error  
429 - parvalue = source.spectral()[0].value() 420 + parvalue = source.spectral()[0].value()
430 if parvalue != 0.0: 421 if parvalue != 0.0:
431 rel_error = source.spectral()[0].error() / parvalue 422 rel_error = source.spectral()[0].error() / parvalue
432 e_flux = fitted_flux * rel_error 423 e_flux = fitted_flux * rel_error
433 else: 424 else:
434 e_flux = 0.0 425 e_flux = 0.0
435 426
436 - # Set values for storage  
437 - TSvalues[i] = TS 427 + # Convert differential flux and flux error to nuFnu
  428 + elogmean2 = elogmean.MeV() * elogmean.MeV()
  429 + result['flux'] = fitted_flux * elogmean2 * gammalib.MeV2erg
  430 + result['flux_err'] = e_flux * elogmean2 * gammalib.MeV2erg
  431 +
  432 + # Log information
  433 + value = '%e +/- %e' % (result['flux'], result['flux_err'])
  434 + if self['calc_ulim'].boolean() and result['ulimit'] > 0.0:
  435 + value += ' [< %e]' % (result['ulimit'])
  436 + value += ' erg/cm2/s'
  437 + if self['calc_ts'].boolean() and result['TS'] > 0.0:
  438 + value += ' (TS = %.3f)' % (result['TS'])
  439 + self._log_value(gammalib.TERSE, 'Bin '+str(i+1), value)
  440 +
  441 + # ... otherwise if logL is zero then signal that bin is
  442 + # skipped
  443 + else:
  444 + value = 'No event in this bin. Likelihood is zero. Bin is skipped.'
  445 + self._log_value(gammalib.TERSE, 'Bin '+str(i+1), value)
438 446
439 - # Set npred values  
440 - Npred_values[i] = Npred 447 + # Return result
  448 + return result
441 449
442 - # Convert fluxes to nuFnu  
443 - flux[i] = fitted_flux * elogmean2 * gammalib.MeV2erg  
444 - flux_err[i] = e_flux * elogmean2 * gammalib.MeV2erg  
445 - if ulimit_value > 0.0:  
446 - ulim_values[i] = ulimit_value * elogmean2 * gammalib.MeV2erg 450 + def _fit_energy_bins(self):
  451 + """
  452 + Fit model to energy bins
447 453
448 - # Log information  
449 - if self._logTerse():  
450 - self._log("\n")  
451 - self._log.parformat("Bin "+str(i+1))  
452 - self._log(str(flux[i]))  
453 - self._log(" +/- ")  
454 - self._log(str(flux_err[i]))  
455 - if self._calc_ulimit and ulim_values[i] > 0.0:  
456 - self._log(" [< "+str(ulim_values[i])+"]")  
457 - self._log(" erg/cm2/s")  
458 - if self._calc_ts and TSvalues[i] > 0.0:  
459 - self._log(" (TS = "+str(TS)+")") 454 + Returns
  455 + -------
  456 + results : list of dict
  457 + List of dictionaries with fit results
  458 + """
  459 + # Write header
  460 + self._log_header1(gammalib.TERSE, 'Generate spectrum')
  461 + self._log_string(gammalib.TERSE, str(self._ebounds))
  462 +
  463 + # Initialise results
  464 + results = []
  465 +
  466 + # Loop over energy bins
  467 + for i in range(self._ebounds.size()):
  468 +
  469 + # Write header for energy bin
  470 + self._log_header2(gammalib.EXPLICIT, 'Energy bin '+str(i+1))
  471 +
  472 + # Fit energy bin
  473 + result = self._fit_energy_bin(i)
  474 +
  475 + # Append results
  476 + results.append(result)
  477 +
  478 + # Return results
  479 + return results
  480 +
  481 + def _create_fits(self, results):
  482 + """
  483 + Create FITS file
  484 +
  485 + Parameters
  486 + ----------
  487 + results : list of dict
  488 + List of result dictionaries
  489 + """
  490 + # Create FITS table columns
  491 + nrows = self._ebounds.size()
  492 + energy = gammalib.GFitsTableDoubleCol('Energy', nrows)
  493 + energy_low = gammalib.GFitsTableDoubleCol('ed_Energy', nrows)
  494 + energy_high = gammalib.GFitsTableDoubleCol('eu_Energy', nrows)
  495 + flux = gammalib.GFitsTableDoubleCol('Flux', nrows)
  496 + flux_err = gammalib.GFitsTableDoubleCol('e_Flux', nrows)
  497 + TSvalues = gammalib.GFitsTableDoubleCol('TS', nrows)
  498 + ulim_values = gammalib.GFitsTableDoubleCol('UpperLimit', nrows)
  499 + Npred_values = gammalib.GFitsTableDoubleCol('Npred', nrows)
  500 + energy.unit('TeV')
  501 + energy_low.unit('TeV')
  502 + energy_high.unit('TeV')
  503 + flux.unit('erg/cm2/s')
  504 + flux_err.unit('erg/cm2/s')
  505 + ulim_values.unit('erg/cm2/s')
  506 +
  507 + # File FITS table columns
  508 + for i, result in enumerate(results):
  509 + energy[i] = result['energy']
  510 + energy_low[i] = result['energy_low']
  511 + energy_high[i] = result['energy_high']
  512 + flux[i] = result['flux']
  513 + flux_err[i] = result['flux_err']
  514 + TSvalues[i] = result['TS']
  515 + ulim_values[i] = result['ulimit']
  516 + Npred_values[i] = result['Npred']
  517 +
  518 + # Initialise FITS Table with extension "SPECTRUM"
  519 + table = gammalib.GFitsBinTable(nrows)
  520 + table.extname('SPECTRUM')
  521 +
  522 + # Add Header for compatibility with gammalib.GMWLSpectrum
  523 + table.card('INSTRUME', 'CTA', 'Name of Instrument')
  524 + table.card('TELESCOP', 'CTA', 'Name of Telescope')
460 525
461 # Append filled columns to fits table 526 # Append filled columns to fits table
462 table.append(energy) 527 table.append(energy)
@@ -472,7 +537,38 @@ class csspec(ctools.cscript): @@ -472,7 +537,38 @@ class csspec(ctools.cscript):
472 self._fits = gammalib.GFits() 537 self._fits = gammalib.GFits()
473 self._fits.append(table) 538 self._fits.append(table)
474 539
475 - # Optionally publish light curve 540 + # Return
  541 + return
  542 +
  543 +
  544 + # Public methods
  545 + def run(self):
  546 + """
  547 + Run the script
  548 + """
  549 + # Switch screen logging on in debug mode
  550 + if self._logDebug():
  551 + self._log.cout(True)
  552 +
  553 + # Get parameters
  554 + self._get_parameters()
  555 +
  556 + # Write input observation container into logger
  557 + self._log_observations(gammalib.NORMAL, self._obs, 'Input observation')
  558 +
  559 + # Adjust model parameters dependent on input user parameters
  560 + self._adjust_model_pars()
  561 +
  562 + # Write spectral binning into logger
  563 + self._log_spectral_binning()
  564 +
  565 + # Fit energy bins
  566 + results = self._fit_energy_bins()
  567 +
  568 + # Create FITS file
  569 + self._create_fits(results)
  570 +
  571 + # Optionally publish spectrum
476 if self['publish'].boolean(): 572 if self['publish'].boolean():
477 self.publish() 573 self.publish()
478 574
@@ -505,17 +601,17 @@ class csspec(ctools.cscript): @@ -505,17 +601,17 @@ class csspec(ctools.cscript):
505 # Write header 601 # Write header
506 self._log_header1(gammalib.TERSE, 'Save spectrum') 602 self._log_header1(gammalib.TERSE, 'Save spectrum')
507 603
508 - # Get outmap parameter  
509 - outfile = self['outfile'].filename()  
510 -  
511 - # Continue only filename and residual map are valid 604 + # Continue only if FITS file is valid
512 if self._fits != None: 605 if self._fits != None:
513 606
  607 + # Get outmap parameter
  608 + outfile = self['outfile'].filename()
  609 +
514 # Log file name 610 # Log file name
515 self._log_value(gammalib.NORMAL, 'Spectrum file', outfile.url()) 611 self._log_value(gammalib.NORMAL, 'Spectrum file', outfile.url())
516 612
517 # Save spectrum 613 # Save spectrum
518 - self._fits.saveto(outfile, self._clobber) 614 + self._fits.saveto(outfile, self['clobber'].boolean())
519 615
520 # Return 616 # Return
521 return 617 return
@@ -532,37 +628,40 @@ class csspec(ctools.cscript): @@ -532,37 +628,40 @@ class csspec(ctools.cscript):
532 # Write header 628 # Write header
533 self._log_header1(gammalib.TERSE, 'Publish spectrum') 629 self._log_header1(gammalib.TERSE, 'Publish spectrum')
534 630
535 - # Continue only if spectrum is valid  
536 - if self._fits.contains('SPECTRUM'): 631 + # Continue only if FITS file is valid
  632 + if self._fits != None:
537 633
538 - # Set default name is user name is empty  
539 - if not name:  
540 - user_name = self._name  
541 - else:  
542 - user_name = name 634 + # Continue only if spectrum is valid
  635 + if self._fits.contains('SPECTRUM'):
543 636
544 - # Log file name  
545 - self._log_value(gammalib.NORMAL, 'Spectrum name', user_name) 637 + # Set default name is user name is empty
  638 + if not name:
  639 + user_name = self._name
  640 + else:
  641 + user_name = name
  642 +
  643 + # Log file name
  644 + self._log_value(gammalib.NORMAL, 'Spectrum name', user_name)
546 645
547 - # Publish spectrum  
548 - self._fits.publish('SPECTRUM', user_name) 646 + # Publish spectrum
  647 + self._fits.publish('SPECTRUM', user_name)
549 648
550 # Return 649 # Return
551 return 650 return
552 651
553 def spectrum(self): 652 def spectrum(self):
554 """ 653 """
555 - Return spectrum FITS file. 654 + Return spectrum FITS file
556 655
557 Returns: 656 Returns:
558 - FITS file containing spectrum. 657 + FITS file containing spectrum
559 """ 658 """
560 # Return 659 # Return
561 return self._fits 660 return self._fits
562 661
563 def models(self, models): 662 def models(self, models):
564 """ 663 """
565 - Set model. 664 + Set model
566 """ 665 """
567 # Copy models 666 # Copy models
568 self._obs.models(models.clone()) 667 self._obs.models(models.clone())