Johan Bregeon / ctools

Commit 29bbc390ffb2cc94e8931e741400c9d70d5625e5

Authored by Jürgen Knödlseder
1 parent 15960584

Review cslightcrv script (#1828)

Inline Side-by-side
Showing 1 changed file with 230 additions and 224 deletions
cscripts/cslightcrv.py
  View file @29bbc39
@@ -28,7 +28,7 @@ import ctools @@ -28,7 +28,7 @@ import ctools
28 # ================ # 28 # ================ #
29 class cslightcrv(ctools.cscript): 29 class cslightcrv(ctools.cscript):
30 """ 30 """
31 - Generates a lightcurve. 31 + Generates a lightcurve
32 32
33 The cslightcrv class generates a light curve for Imaging Air Cherenkov 33 The cslightcrv class generates a light curve for Imaging Air Cherenkov
34 Telescope event data by performing a maximum likelihood fit using 34 Telescope event data by performing a maximum likelihood fit using
@@ -64,14 +64,14 @@ class cslightcrv(ctools.cscript): @@ -64,14 +64,14 @@ class cslightcrv(ctools.cscript):
64 # Constructor 64 # Constructor
65 def __init__(self, *argv): 65 def __init__(self, *argv):
66 """ 66 """
67 - Constructor. 67 + Constructor
68 """ 68 """
69 # Set name 69 # Set name
70 - self._name = "cslightcrv"  
71 - self._version = "1.1.0" 70 + self._name = 'cslightcrv'
  71 + self._version = '1.1.0'
72 72
73 # Initialise some members 73 # Initialise some members
74 - self._srcname = "" 74 + self._srcname = ''
75 self._tbins = gammalib.GGti() 75 self._tbins = gammalib.GGti()
76 self._stacked = False 76 self._stacked = False
77 self._fits = gammalib.GFits() 77 self._fits = gammalib.GFits()
@@ -93,17 +93,17 @@ class cslightcrv(ctools.cscript): @@ -93,17 +93,17 @@ class cslightcrv(ctools.cscript):
93 # Private methods 93 # Private methods
94 def _get_parameters(self): 94 def _get_parameters(self):
95 """ 95 """
96 - Get parameters from parfile. 96 + Get parameters from parfile
97 """ 97 """
98 # Setup observations 98 # Setup observations
99 self._setup_observations(self._obs) 99 self._setup_observations(self._obs)
100 100
101 # Set models if there are none in the container 101 # Set models if there are none in the container
102 if self._obs.models().size() == 0: 102 if self._obs.models().size() == 0:
103 - self._obs.models(self["inmodel"].filename()) 103 + self._obs.models(self['inmodel'].filename())
104 104
105 # Get source name 105 # Get source name
106 - self._srcname = self["srcname"].string() 106 + self._srcname = self['srcname'].string()
107 107
108 # Get time boundaries 108 # Get time boundaries
109 self._tbins = self._create_tbounds() 109 self._tbins = self._create_tbounds()
@@ -111,7 +111,7 @@ class cslightcrv(ctools.cscript): @@ -111,7 +111,7 @@ class cslightcrv(ctools.cscript):
111 # Set stacked analysis flag to True if the requested number of 111 # Set stacked analysis flag to True if the requested number of
112 # energy bins is positive. Otherwise an unbinned analysis will 112 # energy bins is positive. Otherwise an unbinned analysis will
113 # be done and the stacked analysis flag will be set to False. 113 # be done and the stacked analysis flag will be set to False.
114 - if self["enumbins"].integer() > 0: 114 + if self['enumbins'].integer() > 0:
115 self._stacked = True 115 self._stacked = True
116 else: 116 else:
117 self._stacked = False 117 self._stacked = False
@@ -120,61 +120,56 @@ class cslightcrv(ctools.cscript): @@ -120,61 +120,56 @@ class cslightcrv(ctools.cscript):
120 # do not store the actual parameter values as we do not want 120 # do not store the actual parameter values as we do not want
121 # too many instance attributes with enhances the maintenance 121 # too many instance attributes with enhances the maintenance
122 # costs. 122 # costs.
123 - self["emin"].real()  
124 - self["emax"].real() 123 + self['emin'].real()
  124 + self['emax'].real()
125 if self._stacked: 125 if self._stacked:
126 - self["coordsys"].string()  
127 - self["proj"].string()  
128 - self["xref"].real()  
129 - self["yref"].real()  
130 - self["nxpix"].integer()  
131 - self["nypix"].integer()  
132 - self["binsz"].real() 126 + self['coordsys'].string()
  127 + self['proj'].string()
  128 + self['xref'].real()
  129 + self['yref'].real()
  130 + self['nxpix'].integer()
  131 + self['nypix'].integer()
  132 + self['binsz'].real()
133 133
134 # Do the same for the hidden parameters, just in case 134 # Do the same for the hidden parameters, just in case
135 - self["edisp"].boolean()  
136 - self["calc_ulim"].boolean()  
137 - self["calc_ts"].boolean()  
138 - self["fix_bkg"].boolean()  
139 - self["fix_srcs"].boolean() 135 + self['edisp'].boolean()
  136 + self['calc_ulim'].boolean()
  137 + self['calc_ts'].boolean()
  138 + self['fix_bkg'].boolean()
  139 + self['fix_srcs'].boolean()
140 140
141 # Read ahead output parameters 141 # Read ahead output parameters
142 if self._read_ahead(): 142 if self._read_ahead():
143 - self["outfile"].filename() 143 + self['outfile'].filename()
144 144
145 # Write input parameters into logger 145 # Write input parameters into logger
146 if self._logTerse(): 146 if self._logTerse():
147 self._log_parameters() 147 self._log_parameters()
148 - self._log("\n") 148 + self._log('\n')
149 149
150 # Return 150 # Return
151 return 151 return
152 152
153 def _create_tbounds(self): 153 def _create_tbounds(self):
154 """ 154 """
155 - Creates light curve time bins.  
156 -  
157 - The method reads the following user parameters:  
158 - tbinalg: Time binning algorithm  
159 - tbinfile: Time binning file (FITS or ASCII)  
160 - tmin: Start time (MJD)  
161 - tmax: Stop time (MJD)  
162 - tbins: Number of time bins 155 + Creates light curve time bins
163 156
164 - Returns:  
165 - Light curve bins in form of a GTI. 157 + Returns
  158 + -------
  159 + gti : `~gammalib.GGti`
  160 + Light curve bins in form of Good Time Intervals
166 """ 161 """
167 # Initialise Good Time Intervals 162 # Initialise Good Time Intervals
168 gti = gammalib.GGti() 163 gti = gammalib.GGti()
169 164
170 # Get algorithm to use for defining the time intervals 165 # Get algorithm to use for defining the time intervals
171 - algorithm = self["tbinalg"].string() 166 + algorithm = self['tbinalg'].string()
172 167
173 # Handle a FITS or a ASCII file for time bin definition 168 # Handle a FITS or a ASCII file for time bin definition
174 - if algorithm == "FILE": 169 + if algorithm == 'FILE':
175 170
176 # Get the filename 171 # Get the filename
177 - filename = self["tbinfile"].filename() 172 + filename = self['tbinfile'].filename()
178 173
179 # If the file a FITS file then load GTIs 174 # If the file a FITS file then load GTIs
180 if filename.is_fits(): 175 if filename.is_fits():
@@ -190,12 +185,12 @@ class cslightcrv(ctools.cscript): @@ -190,12 +185,12 @@ class cslightcrv(ctools.cscript):
190 gti.append(tmin,tmax) 185 gti.append(tmin,tmax)
191 186
192 # Handle linear time binning 187 # Handle linear time binning
193 - elif algorithm == "LIN": 188 + elif algorithm == 'LIN':
194 189
195 # Get start and stop time and number of time bins 190 # Get start and stop time and number of time bins
196 - time_min = self["tmin"].real()  
197 - time_max = self["tmax"].real()  
198 - nbins = self["tbins"].integer() 191 + time_min = self['tmin'].real()
  192 + time_max = self['tmax'].real()
  193 + nbins = self['tbins'].integer()
199 194
200 # Compute time step and setup time intervals 195 # Compute time step and setup time intervals
201 time_step = (time_max - time_min) / float(nbins) 196 time_step = (time_max - time_min) / float(nbins)
@@ -207,7 +202,7 @@ class cslightcrv(ctools.cscript): @@ -207,7 +202,7 @@ class cslightcrv(ctools.cscript):
207 gti.append(tmin,tmax) 202 gti.append(tmin,tmax)
208 203
209 # Handle usage of observation GTIs 204 # Handle usage of observation GTIs
210 - elif algorithm == "GTI": 205 + elif algorithm == 'GTI':
211 206
212 # Append the GTIs of all observations 207 # Append the GTIs of all observations
213 for obs in self._obs: 208 for obs in self._obs:
@@ -225,9 +220,11 @@ class cslightcrv(ctools.cscript): @@ -225,9 +220,11 @@ class cslightcrv(ctools.cscript):
225 220
226 def _get_free_par_names(self): 221 def _get_free_par_names(self):
227 """ 222 """
228 - Return list of free parameter names. 223 + Return list of free parameter names
229 224
230 - Returns: 225 + Returns
  226 + -------
  227 + names : list of str
231 List of free parameter names. 228 List of free parameter names.
232 """ 229 """
233 # Initialise list of free parameter names 230 # Initialise list of free parameter names
@@ -243,75 +240,78 @@ class cslightcrv(ctools.cscript): @@ -243,75 +240,78 @@ class cslightcrv(ctools.cscript):
243 240
244 def _adjust_model_pars(self): 241 def _adjust_model_pars(self):
245 """ 242 """
246 - Adjust model parameters dependent on user parameters. 243 + Adjust model parameters dependent on user parameters
247 """ 244 """
248 # Write header 245 # Write header
249 if self._logTerse(): 246 if self._logTerse():
250 - self._log("\n")  
251 - self._log.header1("Adjust model parameters") 247 + self._log('\n')
  248 + self._log.header1('Adjust model parameters')
252 249
253 # Adjust model parameters dependent on input user parameters 250 # Adjust model parameters dependent on input user parameters
254 for model in self._obs.models(): 251 for model in self._obs.models():
255 252
256 - # Set TS flag for all models to false. The source of interest  
257 - # will be set to true later 253 + # Set TS flag for all models to false. The source of interest will
  254 + # be set to true later
258 model.tscalc(False) 255 model.tscalc(False)
259 256
  257 + # Get the model classname to distinguish later sky from background
  258 + # models
  259 + classname = model.classname()
  260 +
260 # Log model name 261 # Log model name
261 if self._logNormal(): 262 if self._logNormal():
262 self._log.header3(model.name()) 263 self._log.header3(model.name())
263 264
264 - # Deal with the source of interest 265 + # If the model is the source of interest and the 'calc_ts' parameter
  266 + # is true then enable the TS computation for the source
265 if model.name() == self._srcname: 267 if model.name() == self._srcname:
266 - if self["calc_ts"].boolean(): 268 + if self['calc_ts'].boolean():
267 model.tscalc(True) 269 model.tscalc(True)
268 270
269 - elif (self["fix_bkg"].boolean() and  
270 - model.classname() != "GModelSky"):  
271 - for par in model:  
272 - if par.is_free():  
273 - par.fix()  
274 - if self._logNormal():  
275 - self._log(gammalib.parformat(par.name()))  
276 - self._log("fixed\n")  
277 -  
278 - elif (self["fix_srcs"].boolean() and  
279 - model.classname() == "GModelSky"): 271 + # ... otherwise, if the model is not a sky model and the 'fix_bkg'
  272 + # parameter is true or if the model is a sky model and the
  273 + # 'fix_srcs' parameter is true then fix all parameters of the model
  274 + elif ((self['fix_bkg'].boolean() and classname != 'GModelSky') or
  275 + (self['fix_srcs'].boolean() and classname == 'GModelSky')):
280 for par in model: 276 for par in model:
281 if par.is_free(): 277 if par.is_free():
282 par.fix() 278 par.fix()
283 if self._logNormal(): 279 if self._logNormal():
284 self._log(gammalib.parformat(par.name())) 280 self._log(gammalib.parformat(par.name()))
285 - self._log("fixed\n") 281 + self._log('fixed\n')
286 282
287 # Return 283 # Return
288 return 284 return
289 285
290 def _create_fits_table(self, results): 286 def _create_fits_table(self, results):
291 """ 287 """
292 - Creates FITS binary table containing light curve results. 288 + Creates FITS binary table containing light curve results
293 289
294 - Args:  
295 - results: List of result dictionaries. 290 + Parameters
  291 + ----------
  292 + results : list of dict
  293 + List of result dictionaries
296 294
297 - Returns:  
298 - FITS binary table containing light curve. 295 + Returns
  296 + -------
  297 + table : `~gammalib.GFitsBinTable`
  298 + FITS binary table containing light curve
299 """ 299 """
300 # Determine number of rows in FITS table 300 # Determine number of rows in FITS table
301 nrows = len(results) 301 nrows = len(results)
302 302
303 # Create FITS Table with extension "LIGHTCURVE" 303 # Create FITS Table with extension "LIGHTCURVE"
304 table = gammalib.GFitsBinTable(nrows) 304 table = gammalib.GFitsBinTable(nrows)
305 - table.extname("LIGHTCURVE") 305 + table.extname('LIGHTCURVE')
306 306
307 # Append time columns 307 # Append time columns
308 - mjd = gammalib.GFitsTableDoubleCol("MJD", nrows)  
309 - e_mjd = gammalib.GFitsTableDoubleCol("e_MJD", nrows)  
310 - mjd.unit("days")  
311 - e_mjd.unit("days") 308 + mjd = gammalib.GFitsTableDoubleCol('MJD', nrows)
  309 + e_mjd = gammalib.GFitsTableDoubleCol('e_MJD', nrows)
  310 + mjd.unit('days')
  311 + e_mjd.unit('days')
312 for i, result in enumerate(results): 312 for i, result in enumerate(results):
313 - mjd[i] = result["mjd"]  
314 - e_mjd[i] = result["e_mjd"] 313 + mjd[i] = result['mjd']
  314 + e_mjd[i] = result['e_mjd']
315 table.append(mjd) 315 table.append(mjd)
316 table.append(e_mjd) 316 table.append(e_mjd)
317 317
@@ -319,7 +319,7 @@ class cslightcrv(ctools.cscript): @@ -319,7 +319,7 @@ class cslightcrv(ctools.cscript):
319 for par in self._obs.models()[self._srcname]: 319 for par in self._obs.models()[self._srcname]:
320 if par.is_free(): 320 if par.is_free():
321 name = par.name() 321 name = par.name()
322 - e_name = "e_"+par.name() 322 + e_name = 'e_'+par.name()
323 col = gammalib.GFitsTableDoubleCol(name, nrows) 323 col = gammalib.GFitsTableDoubleCol(name, nrows)
324 e_col = gammalib.GFitsTableDoubleCol(e_name, nrows) 324 e_col = gammalib.GFitsTableDoubleCol(e_name, nrows)
325 col.unit(par.unit()) 325 col.unit(par.unit())
@@ -331,14 +331,14 @@ class cslightcrv(ctools.cscript): @@ -331,14 +331,14 @@ class cslightcrv(ctools.cscript):
331 table.append(e_col) 331 table.append(e_col)
332 332
333 # Append Test Statistic column 333 # Append Test Statistic column
334 - ts = gammalib.GFitsTableDoubleCol("TS", nrows) 334 + ts = gammalib.GFitsTableDoubleCol('TS', nrows)
335 for i, result in enumerate(results): 335 for i, result in enumerate(results):
336 - ts[i] = result["ts"] 336 + ts[i] = result['ts']
337 table.append(ts) 337 table.append(ts)
338 338
339 # Append upper limit column 339 # Append upper limit column
340 - ulimit = gammalib.GFitsTableDoubleCol("UpperLimit", nrows)  
341 - ulimit.unit("ph/cm2/s") 340 + ulimit = gammalib.GFitsTableDoubleCol('UpperLimit', nrows)
  341 + ulimit.unit('ph/cm2/s')
342 for i, result in enumerate(results): 342 for i, result in enumerate(results):
343 ulimit[i] = result['ulimit'] 343 ulimit[i] = result['ulimit']
344 table.append(ulimit) 344 table.append(ulimit)
@@ -348,28 +348,32 @@ class cslightcrv(ctools.cscript): @@ -348,28 +348,32 @@ class cslightcrv(ctools.cscript):
348 348
349 def _compute_ulimit(self, obs): 349 def _compute_ulimit(self, obs):
350 """ 350 """
351 - Computes upper flux limit. 351 + Computes upper flux limit
352 352
353 - Args:  
354 - obs: Observation container. 353 + Parameters
  354 + ----------
  355 + obs : `~gammalib.GObservations`
  356 + Observation container
355 357
356 - Returns:  
357 - Upper flux limit (-1 of not computed). 358 + Returns
  359 + -------
  360 + ulimit : float
  361 + Upper flux limit (-1.0 of not computed)
358 """ 362 """
359 # Initialise upper flux limit 363 # Initialise upper flux limit
360 ulimit_value = -1.0 364 ulimit_value = -1.0
361 365
362 # Perform computation only if requested 366 # Perform computation only if requested
363 - if self["calc_ulim"].boolean(): 367 + if self['calc_ulim'].boolean():
364 368
365 # Write header in logger 369 # Write header in logger
366 if self._logExplicit(): 370 if self._logExplicit():
367 - self._log.header3("Computing upper flux limit") 371 + self._log.header3('Computing upper flux limit')
368 372
369 # Create upper limit object 373 # Create upper limit object
370 ulimit = ctools.ctulimit(obs) 374 ulimit = ctools.ctulimit(obs)
371 - ulimit["srcname"] = self._srcname  
372 - ulimit["eref"] = 1.0 375 + ulimit['srcname'] = self._srcname
  376 + ulimit['eref'] = 1.0
373 377
374 # Try to run upper limit and catch exceptions 378 # Try to run upper limit and catch exceptions
375 try: 379 try:
@@ -377,7 +381,7 @@ class cslightcrv(ctools.cscript): @@ -377,7 +381,7 @@ class cslightcrv(ctools.cscript):
377 ulimit_value = ulimit.flux_ulimit() 381 ulimit_value = ulimit.flux_ulimit()
378 except: 382 except:
379 if self._logTerse(): 383 if self._logTerse():
380 - self._log("Upper limit flux calculation failed.s\n") 384 + self._log('Upper limit flux calculation failed.\n')
381 ulimit_value = -1.0 385 ulimit_value = -1.0
382 386
383 # Return upper limit 387 # Return upper limit
@@ -385,129 +389,129 @@ class cslightcrv(ctools.cscript): @@ -385,129 +389,129 @@ class cslightcrv(ctools.cscript):
385 389
386 def _bin_observation(self, obs): 390 def _bin_observation(self, obs):
387 """ 391 """
388 - Bin an observation is a binned analysis was requested. 392 + Bin an observation if a binned analysis was requested
389 393
390 - Args:  
391 - obs: Observation container. 394 + Parameters
  395 + ----------
  396 + obs : `~gammalib.GObservations`
  397 + Observation container
392 398
393 - Returns:  
394 - Observation container with a binned. 399 + Returns
  400 + -------
  401 + obs : `~gammalib.GObservations`
  402 + Observation container where the first observation is a binned observation
395 """ 403 """
396 # Header 404 # Header
397 if self._logExplicit(): 405 if self._logExplicit():
398 - self._log.header3("Binning events") 406 + self._log.header3('Binning events')
399 407
400 # Bin events 408 # Bin events
401 cntcube = ctools.ctbin(obs) 409 cntcube = ctools.ctbin(obs)
402 - cntcube["usepnt"] = False  
403 - cntcube["ebinalg"] = "LOG"  
404 - cntcube["xref"] = self["xref"].real()  
405 - cntcube["yref"] = self["yref"].real()  
406 - cntcube["binsz"] = self["binsz"].real()  
407 - cntcube["nxpix"] = self["nxpix"].integer()  
408 - cntcube["nypix"] = self["nypix"].integer()  
409 - cntcube["enumbins"] = self["enumbins"].integer()  
410 - cntcube["emin"] = self["emin"].real()  
411 - cntcube["emax"] = self["emax"].real()  
412 - cntcube["coordsys"] = self["coordsys"].string()  
413 - cntcube["proj"] = self["proj"].string() 410 + cntcube['usepnt'] = False
  411 + cntcube['ebinalg'] = 'LOG'
  412 + cntcube['xref'] = self['xref'].real()
  413 + cntcube['yref'] = self['yref'].real()
  414 + cntcube['binsz'] = self['binsz'].real()
  415 + cntcube['nxpix'] = self['nxpix'].integer()
  416 + cntcube['nypix'] = self['nypix'].integer()
  417 + cntcube['enumbins'] = self['enumbins'].integer()
  418 + cntcube['emin'] = self['emin'].real()
  419 + cntcube['emax'] = self['emax'].real()
  420 + cntcube['coordsys'] = self['coordsys'].string()
  421 + cntcube['proj'] = self['proj'].string()
414 cntcube.run() 422 cntcube.run()
415 423
416 # Header 424 # Header
417 if self._logExplicit(): 425 if self._logExplicit():
418 - self._log.header3("Creating exposure cube") 426 + self._log.header3('Creating exposure cube')
419 427
420 # Create exposure cube 428 # Create exposure cube
421 expcube = ctools.ctexpcube(obs) 429 expcube = ctools.ctexpcube(obs)
422 - expcube["incube"] = "NONE"  
423 - expcube["usepnt"] = False  
424 - expcube["ebinalg"] = "LOG"  
425 - expcube["xref"] = self["xref"].real()  
426 - expcube["yref"] = self["yref"].real()  
427 - expcube["binsz"] = self["binsz"].real()  
428 - expcube["nxpix"] = self["nxpix"].integer()  
429 - expcube["nypix"] = self["nypix"].integer()  
430 - expcube["enumbins"] = self["enumbins"].integer()  
431 - expcube["emin"] = self["emin"].real()  
432 - expcube["emax"] = self["emax"].real()  
433 - expcube["coordsys"] = self["coordsys"].string()  
434 - expcube["proj"] = self["proj"].string() 430 + expcube['incube'] = 'NONE'
  431 + expcube['usepnt'] = False
  432 + expcube['ebinalg'] = 'LOG'
  433 + expcube['xref'] = self['xref'].real()
  434 + expcube['yref'] = self['yref'].real()
  435 + expcube['binsz'] = self['binsz'].real()
  436 + expcube['nxpix'] = self['nxpix'].integer()
  437 + expcube['nypix'] = self['nypix'].integer()
  438 + expcube['enumbins'] = self['enumbins'].integer()
  439 + expcube['emin'] = self['emin'].real()
  440 + expcube['emax'] = self['emax'].real()
  441 + expcube['coordsys'] = self['coordsys'].string()
  442 + expcube['proj'] = self['proj'].string()
435 expcube.run() 443 expcube.run()
436 444
437 # Header 445 # Header
438 if self._logExplicit(): 446 if self._logExplicit():
439 - self._log.header3("Creating point spread function cube") 447 + self._log.header3('Creating point spread function cube')
440 448
441 # Compute spatial binning for point spread function and 449 # Compute spatial binning for point spread function and
442 # energy dispersion cubes 450 # energy dispersion cubes
443 - binsz = 10.0 * self["binsz"].real()  
444 - nxpix = self["nxpix"].integer() // 10 # Make sure result is int  
445 - nypix = self["nypix"].integer() // 10 # Make sure result is int  
446 - if nxpix < 2:  
447 - nxpix = 2  
448 - if nypix < 2:  
449 - nypix = 2 451 + binsz = 10.0 * self['binsz'].real()
  452 + nxpix = min(self['nxpix'].integer() // 10, 2) # Make sure result is int
  453 + nypix = min(self['nypix'].integer() // 10, 2) # Make sure result is int
450 454
451 # Create point spread function cube 455 # Create point spread function cube
452 psfcube = ctools.ctpsfcube(obs) 456 psfcube = ctools.ctpsfcube(obs)
453 - psfcube["incube"] = "NONE"  
454 - psfcube["usepnt"] = False  
455 - psfcube["ebinalg"] = "LOG"  
456 - psfcube["xref"] = self["xref"].real()  
457 - psfcube["yref"] = self["yref"].real()  
458 - psfcube["binsz"] = binsz  
459 - psfcube["nxpix"] = nxpix  
460 - psfcube["nypix"] = nypix  
461 - psfcube["enumbins"] = self["enumbins"].integer()  
462 - psfcube["emin"] = self["emin"].real()  
463 - psfcube["emax"] = self["emax"].real()  
464 - psfcube["coordsys"] = self["coordsys"].string()  
465 - psfcube["proj"] = self["proj"].string() 457 + psfcube['incube'] = 'NONE'
  458 + psfcube['usepnt'] = False
  459 + psfcube['ebinalg'] = 'LOG'
  460 + psfcube['xref'] = self['xref'].real()
  461 + psfcube['yref'] = self['yref'].real()
  462 + psfcube['binsz'] = binsz
  463 + psfcube['nxpix'] = nxpix
  464 + psfcube['nypix'] = nypix
  465 + psfcube['enumbins'] = self['enumbins'].integer()
  466 + psfcube['emin'] = self['emin'].real()
  467 + psfcube['emax'] = self['emax'].real()
  468 + psfcube['coordsys'] = self['coordsys'].string()
  469 + psfcube['proj'] = self['proj'].string()
466 psfcube.run() 470 psfcube.run()
467 471
468 # Check if we need to include energy dispersion 472 # Check if we need to include energy dispersion
469 - if self["edisp"].boolean(): 473 + if self['edisp'].boolean():
470 474
471 # Header 475 # Header
472 if self._logExplicit(): 476 if self._logExplicit():
473 - self._log.header3("Creating energy dispersion cube") 477 + self._log.header3('Creating energy dispersion cube')
474 478
475 # Create energy dispersion cube 479 # Create energy dispersion cube
476 edispcube = ctools.ctedispcube(obs) 480 edispcube = ctools.ctedispcube(obs)
477 - edispcube["incube"] = "NONE"  
478 - edispcube["usepnt"] = False  
479 - edispcube["ebinalg"] = "LOG"  
480 - edispcube["xref"] = self["xref"].real()  
481 - edispcube["yref"] = self["yref"].real()  
482 - edispcube["binsz"] = binsz  
483 - edispcube["nxpix"] = nxpix  
484 - edispcube["nypix"] = nypix  
485 - edispcube["enumbins"] = self["enumbins"].integer()  
486 - edispcube["emin"] = self["emin"].real()  
487 - edispcube["emax"] = self["emax"].real()  
488 - edispcube["coordsys"] = self["coordsys"].string()  
489 - edispcube["proj"] = self["proj"].string() 481 + edispcube['incube'] = 'NONE'
  482 + edispcube['usepnt'] = False
  483 + edispcube['ebinalg'] = 'LOG'
  484 + edispcube['xref'] = self['xref'].real()
  485 + edispcube['yref'] = self['yref'].real()
  486 + edispcube['binsz'] = binsz
  487 + edispcube['nxpix'] = nxpix
  488 + edispcube['nypix'] = nypix
  489 + edispcube['enumbins'] = self['enumbins'].integer()
  490 + edispcube['emin'] = self['emin'].real()
  491 + edispcube['emax'] = self['emax'].real()
  492 + edispcube['coordsys'] = self['coordsys'].string()
  493 + edispcube['proj'] = self['proj'].string()
490 edispcube.run() 494 edispcube.run()
491 495
492 # Header 496 # Header
493 if self._logExplicit(): 497 if self._logExplicit():
494 - self._log.header3("Creating background cube") 498 + self._log.header3('Creating background cube')
495 499
496 # Create background cube 500 # Create background cube
497 bkgcube = ctools.ctbkgcube(obs) 501 bkgcube = ctools.ctbkgcube(obs)
498 - bkgcube["incube"] = "NONE"  
499 - bkgcube["usepnt"] = False  
500 - bkgcube["ebinalg"] = "LOG"  
501 - bkgcube["xref"] = self["xref"].real()  
502 - bkgcube["yref"] = self["yref"].real()  
503 - bkgcube["binsz"] = self["binsz"].real()  
504 - bkgcube["nxpix"] = self["nxpix"].integer()  
505 - bkgcube["nypix"] = self["nypix"].integer()  
506 - bkgcube["enumbins"] = self["enumbins"].integer()  
507 - bkgcube["emin"] = self["emin"].real()  
508 - bkgcube["emax"] = self["emax"].real()  
509 - bkgcube["coordsys"] = self["coordsys"].string()  
510 - bkgcube["proj"] = self["proj"].string() 502 + bkgcube['incube'] = 'NONE'
  503 + bkgcube['usepnt'] = False
  504 + bkgcube['ebinalg'] = 'LOG'
  505 + bkgcube['xref'] = self['xref'].real()
  506 + bkgcube['yref'] = self['yref'].real()
  507 + bkgcube['binsz'] = self['binsz'].real()
  508 + bkgcube['nxpix'] = self['nxpix'].integer()
  509 + bkgcube['nypix'] = self['nypix'].integer()
  510 + bkgcube['enumbins'] = self['enumbins'].integer()
  511 + bkgcube['emin'] = self['emin'].real()
  512 + bkgcube['emax'] = self['emax'].real()
  513 + bkgcube['coordsys'] = self['coordsys'].string()
  514 + bkgcube['proj'] = self['proj'].string()
511 bkgcube.run() 515 bkgcube.run()
512 516
513 # Retrieve a new oberservation container 517 # Retrieve a new oberservation container
@@ -517,20 +521,17 @@ class cslightcrv(ctools.cscript): @@ -517,20 +521,17 @@ class cslightcrv(ctools.cscript):
517 models = bkgcube.models() 521 models = bkgcube.models()
518 522
519 # Set stacked response 523 # Set stacked response
520 - if self["edisp"].boolean():  
521 - new_obs[0].response(expcube.expcube(),  
522 - psfcube.psfcube(),  
523 - edispcube.edispcube(),  
524 - bkgcube.bkgcube()) 524 + if self['edisp'].boolean():
  525 + new_obs[0].response(expcube.expcube(), psfcube.psfcube(),
  526 + edispcube.edispcube(), bkgcube.bkgcube())
525 else: 527 else:
526 - new_obs[0].response(expcube.expcube(),  
527 - psfcube.psfcube(), 528 + new_obs[0].response(expcube.expcube(), psfcube.psfcube(),
528 bkgcube.bkgcube()) 529 bkgcube.bkgcube())
529 530
530 # Fix background models if required 531 # Fix background models if required
531 - if self["fix_bkg"].boolean(): 532 + if self['fix_bkg'].boolean():
532 for model in models: 533 for model in models:
533 - if model.classname() != "GModelSky": 534 + if model.classname() != 'GModelSky':
534 for par in model: 535 for par in model:
535 par.fix() 536 par.fix()
536 537
@@ -544,7 +545,7 @@ class cslightcrv(ctools.cscript): @@ -544,7 +545,7 @@ class cslightcrv(ctools.cscript):
544 # Public methods 545 # Public methods
545 def run(self): 546 def run(self):
546 """ 547 """
547 - Run the script. 548 + Run the script
548 """ 549 """
549 # Switch screen logging on in debug mode 550 # Switch screen logging on in debug mode
550 if self._logDebug(): 551 if self._logDebug():
@@ -555,18 +556,18 @@ class cslightcrv(ctools.cscript): @@ -555,18 +556,18 @@ class cslightcrv(ctools.cscript):
555 556
556 # Write observation into logger 557 # Write observation into logger
557 if self._logTerse(): 558 if self._logTerse():
558 - self._log("\n")  
559 - self._log.header1("Observation") 559 + self._log('\n')
  560 + self._log.header1(gammalib.number('Observation',len(self._obs)))
560 self._log(str(self._obs)) 561 self._log(str(self._obs))
561 - self._log("\n") 562 + self._log('\n')
562 563
563 # Adjust model parameters dependent on user parameters 564 # Adjust model parameters dependent on user parameters
564 self._adjust_model_pars() 565 self._adjust_model_pars()
565 566
566 # Write header 567 # Write header
567 if self._logTerse(): 568 if self._logTerse():
568 - self._log("\n")  
569 - self._log.header1("Generate lightcurve") 569 + self._log('\n')
  570 + self._log.header1('Generate lightcurve')
570 571
571 # Initialise list of result dictionaries 572 # Initialise list of result dictionaries
572 results = [] 573 results = []
@@ -583,9 +584,7 @@ class cslightcrv(ctools.cscript): @@ -583,9 +584,7 @@ class cslightcrv(ctools.cscript):
583 584
584 # Write time bin into header 585 # Write time bin into header
585 if self._logTerse(): 586 if self._logTerse():
586 - self._log.header2("MJD "+  
587 - str(tmin.mjd())+"-"+  
588 - str(tmax.mjd())) 587 + self._log.header2('MJD '+str(tmin.mjd())+'-'+str(tmax.mjd()))
589 588
590 # Compute time bin center and time width 589 # Compute time bin center and time width
591 twidth = 0.5 * (tmax - tmin) # in seconds 590 twidth = 0.5 * (tmax - tmin) # in seconds
@@ -601,17 +600,17 @@ class cslightcrv(ctools.cscript): @@ -601,17 +600,17 @@ class cslightcrv(ctools.cscript):
601 600
602 # Log information 601 # Log information
603 if self._logExplicit(): 602 if self._logExplicit():
604 - self._log.header3("Selecting events") 603 + self._log.header3('Selecting events')
605 604
606 # Select events 605 # Select events
607 select = ctools.ctselect(self._obs) 606 select = ctools.ctselect(self._obs)
608 - select["emin"] = self["emin"].real()  
609 - select["emax"] = self["emax"].real()  
610 - select["tmin"] = tmin.convert(self._time_reference())  
611 - select["tmax"] = tmax.convert(self._time_reference())  
612 - select["rad"] = "UNDEFINED"  
613 - select["ra"] = "UNDEFINED"  
614 - select["dec"] = "UNDEFINED" 607 + select['emin'] = self['emin'].real()
  608 + select['emax'] = self['emax'].real()
  609 + select['tmin'] = tmin.convert(self._time_reference())
  610 + select['tmax'] = tmax.convert(self._time_reference())
  611 + select['rad'] = 'UNDEFINED'
  612 + select['ra'] = 'UNDEFINED'
  613 + select['dec'] = 'UNDEFINED'
615 select.run() 614 select.run()
616 615
617 # Retrieve observation 616 # Retrieve observation
@@ -625,11 +624,11 @@ class cslightcrv(ctools.cscript): @@ -625,11 +624,11 @@ class cslightcrv(ctools.cscript):
625 624
626 # Header 625 # Header
627 if self._logExplicit(): 626 if self._logExplicit():
628 - self._log.header3("Fitting the data") 627 + self._log.header3('Fitting the data')
629 628
630 # Do maximum likelihood model fitting 629 # Do maximum likelihood model fitting
631 like = ctools.ctlike(obs) 630 like = ctools.ctlike(obs)
632 - like["edisp"] = self["edisp"].boolean() 631 + like['edisp'] = self['edisp'].boolean()
633 like.run() 632 like.run()
634 633
635 # Skip bin if no event was present 634 # Skip bin if no event was present
@@ -637,13 +636,13 @@ class cslightcrv(ctools.cscript): @@ -637,13 +636,13 @@ class cslightcrv(ctools.cscript):
637 636
638 # Signal skipping of bin 637 # Signal skipping of bin
639 if self._logTerse(): 638 if self._logTerse():
640 - self._log(gammalib.parformat("Warning"))  
641 - self._log("No event in this time bin, skip bin.\n") 639 + self._log(gammalib.parformat('Warning'))
  640 + self._log('No event in this time bin, skip bin.\n')
642 641
643 # Set all results to 0 642 # Set all results to 0
644 for par in pars: 643 for par in pars:
645 result['values'][par] = 0.0 644 result['values'][par] = 0.0
646 - result['values']["e_"+par] = 0.0 645 + result['values']['e_'+par] = 0.0
647 646
648 # Append result 647 # Append result
649 results.append(result) 648 results.append(result)
@@ -657,7 +656,7 @@ class cslightcrv(ctools.cscript): @@ -657,7 +656,7 @@ class cslightcrv(ctools.cscript):
657 # Extract parameter values 656 # Extract parameter values
658 for par in pars: 657 for par in pars:
659 result['values'][par] = source.spectral()[par].value() 658 result['values'][par] = source.spectral()[par].value()
660 - result['values']["e_"+par] = source.spectral()[par].error() 659 + result['values']['e_'+par] = source.spectral()[par].error()
661 660
662 # Calculate upper limit (-1 if not computed) 661 # Calculate upper limit (-1 if not computed)
663 ulimit_value = self._compute_ulimit(like.obs()) 662 ulimit_value = self._compute_ulimit(like.obs())
@@ -665,7 +664,7 @@ class cslightcrv(ctools.cscript): @@ -665,7 +664,7 @@ class cslightcrv(ctools.cscript):
665 result['ulimit'] = ulimit_value 664 result['ulimit'] = ulimit_value
666 665
667 # Extract Test Statistic value 666 # Extract Test Statistic value
668 - if self["calc_ts"].boolean(): 667 + if self['calc_ts'].boolean():
669 result['ts'] = source.ts() 668 result['ts'] = source.ts()
670 669
671 # Append result to list of dictionaries 670 # Append result to list of dictionaries
@@ -673,7 +672,7 @@ class cslightcrv(ctools.cscript): @@ -673,7 +672,7 @@ class cslightcrv(ctools.cscript):
673 672
674 # Log results for this time bin 673 # Log results for this time bin
675 if self._logNormal(): 674 if self._logNormal():
676 - self._log.header3("Results") 675 + self._log.header3('Results')
677 pars = self._get_free_par_names() 676 pars = self._get_free_par_names()
678 for par in pars: 677 for par in pars:
679 value = source.spectral()[par].value() 678 value = source.spectral()[par].value()
@@ -681,17 +680,17 @@ class cslightcrv(ctools.cscript): @@ -681,17 +680,17 @@ class cslightcrv(ctools.cscript):
681 unit = source.spectral()[par].unit() 680 unit = source.spectral()[par].unit()
682 self._log(gammalib.parformat(par)) 681 self._log(gammalib.parformat(par))
683 self._log(str(value)) 682 self._log(str(value))
684 - self._log(" +/- ") 683 + self._log(' +/- ')
685 self._log(str(error)) 684 self._log(str(error))
686 - self._log(" ") 685 + self._log(' ')
687 self._log(unit) 686 self._log(unit)
688 - self._log("\n") 687 + self._log('\n')
689 if result['ulimit'] > 0.0: 688 if result['ulimit'] > 0.0:
690 - self._log(gammalib.parformat("Upper flux limit"))  
691 - self._log(str(result['ulimit'])+" ph/cm2/s\n")  
692 - if self["calc_ts"].boolean():  
693 - self._log(gammalib.parformat("Test Statistic"))  
694 - self._log(str(result['ts'])+"\n") 689 + self._log(gammalib.parformat('Upper flux limit'))
  690 + self._log(str(result['ulimit'])+' ph/cm2/s\n')
  691 + if self['calc_ts'].boolean():
  692 + self._log(gammalib.parformat('Test Statistic'))
  693 + self._log(str(result['ts'])+'\n')
695 694
696 # Create FITS table from results 695 # Create FITS table from results
697 table = self._create_fits_table(results) 696 table = self._create_fits_table(results)
@@ -705,7 +704,7 @@ class cslightcrv(ctools.cscript): @@ -705,7 +704,7 @@ class cslightcrv(ctools.cscript):
705 704
706 def execute(self): 705 def execute(self):
707 """ 706 """
708 - Execute the script. 707 + Execute the script
709 """ 708 """
710 # Open logfile 709 # Open logfile
711 self.logFileOpen() 710 self.logFileOpen()
@@ -724,7 +723,7 @@ class cslightcrv(ctools.cscript): @@ -724,7 +723,7 @@ class cslightcrv(ctools.cscript):
724 723
725 def save(self): 724 def save(self):
726 """ 725 """
727 - Save light curve. 726 + Save light curve
728 """ 727 """
729 # Write header 728 # Write header
730 if self._logTerse(): 729 if self._logTerse():
@@ -732,16 +731,16 @@ class cslightcrv(ctools.cscript): @@ -732,16 +731,16 @@ class cslightcrv(ctools.cscript):
732 self._log.header1('Save light curve') 731 self._log.header1('Save light curve')
733 732
734 # Get light curve filename 733 # Get light curve filename
735 - outfile = self["outfile"].filename() 734 + outfile = self['outfile'].filename()
736 735
737 # Continue only filename and residual map are valid 736 # Continue only filename and residual map are valid
738 if self._fits != None: 737 if self._fits != None:
739 738
740 # Log file name 739 # Log file name
741 if self._logTerse(): 740 if self._logTerse():
742 - self._log(gammalib.parformat("Light curve file")) 741 + self._log(gammalib.parformat('Light curve file'))
743 self._log(outfile.url()) 742 self._log(outfile.url())
744 - self._log("\n") 743 + self._log('\n')
745 744
746 # Save spectrum 745 # Save spectrum
747 self._fits.saveto(outfile, self._clobber()) 746 self._fits.saveto(outfile, self._clobber())
@@ -753,18 +752,25 @@ class cslightcrv(ctools.cscript): @@ -753,18 +752,25 @@ class cslightcrv(ctools.cscript):
753 """ 752 """
754 Return light curve FITS file. 753 Return light curve FITS file.
755 754
756 - Returns:  
757 - FITS file containing light curve. 755 + Returns
  756 + -------
  757 + fits : `~gammalib.GFits()`
  758 + FITS file containing light curve
758 """ 759 """
759 # Return 760 # Return
760 return self._fits 761 return self._fits
761 762
762 def models(self, models): 763 def models(self, models):
763 """ 764 """
764 - Set model. 765 + Set model
  766 +
  767 + Parameters
  768 + ----------
  769 + models : `~gammalib.GModels()`
  770 + Set model container
765 """ 771 """
766 # Copy models 772 # Copy models
767 - self._obs.models(models.clone()) 773 + self._obs.models(models)
768 774
769 # Return 775 # Return
770 return 776 return