Dmitry Malyshev / ctools

Commit 003f01ccbc5dac3a9ebb7b86a8a4b71733159390

Authored by Jürgen Knödlseder
1 parent a90f5fd4

Improve quality of cssens script (#1828)

Inline Side-by-side
Showing 5 changed files with 252 additions and 249 deletions
cscripts/csfindobs.py
  View file @003f01c
@@ -40,7 +40,7 @@ class csfindobs(ctools.cscript): @@ -40,7 +40,7 @@ class csfindobs(ctools.cscript):
40 """ 40 """
41 # Set name 41 # Set name
42 self._name = 'csfindobs' 42 self._name = 'csfindobs'
43 - self._version = '1.1.0' 43 + self._version = '1.2.0'
44 self._datapath = os.getenv('VHEFITS','') 44 self._datapath = os.getenv('VHEFITS','')
45 self._prodname = '' 45 self._prodname = ''
46 self._select_radec = True 46 self._select_radec = True
cscripts/cslightcrv.py
  View file @003f01c
@@ -68,7 +68,7 @@ class cslightcrv(ctools.cscript): @@ -68,7 +68,7 @@ class cslightcrv(ctools.cscript):
68 """ 68 """
69 # Set name 69 # Set name
70 self._name = 'cslightcrv' 70 self._name = 'cslightcrv'
71 - self._version = '1.1.0' 71 + self._version = '1.2.0'
72 72
73 # Initialise some members 73 # Initialise some members
74 self._srcname = '' 74 self._srcname = ''
cscripts/csresmap.py
  View file @003f01c
@@ -38,7 +38,7 @@ class csresmap(ctools.cscript): @@ -38,7 +38,7 @@ class csresmap(ctools.cscript):
38 """ 38 """
39 # Set name 39 # Set name
40 self._name = 'csresmap' 40 self._name = 'csresmap'
41 - self._version = '1.1.0' 41 + self._version = '1.2.0'
42 42
43 # Initialise class members 43 # Initialise class members
44 self._resmap = None 44 self._resmap = None
cscripts/cssens.py
  View file @003f01c
@@ -31,7 +31,7 @@ from cscripts import obsutils @@ -31,7 +31,7 @@ from cscripts import obsutils
31 # ============ # 31 # ============ #
32 class cssens(ctools.cscript): 32 class cssens(ctools.cscript):
33 """ 33 """
34 - Computes the CTA sensitivity. 34 + Computes the CTA sensitivity
35 35
36 This class computes the CTA sensitivity for a number of energy bins using 36 This class computes the CTA sensitivity for a number of energy bins using
37 ctlike. Spectra are fitted in narrow energy bins to simulated data, 37 ctlike. Spectra are fitted in narrow energy bins to simulated data,
@@ -47,14 +47,14 @@ class cssens(ctools.cscript): @@ -47,14 +47,14 @@ class cssens(ctools.cscript):
47 Constructor. 47 Constructor.
48 """ 48 """
49 # Set name 49 # Set name
50 - self._name = "cssens"  
51 - self._version = "1.1.0" 50 + self._name = 'cssens'
  51 + self._version = '1.2.0'
52 52
53 # Initialise class members 53 # Initialise class members
54 self._obs = gammalib.GObservations() 54 self._obs = gammalib.GObservations()
55 self._ebounds = gammalib.GEbounds() 55 self._ebounds = gammalib.GEbounds()
56 self._obs_ebounds = [] 56 self._obs_ebounds = []
57 - self._srcname = "" 57 + self._srcname = ''
58 self._outfile = gammalib.GFilename() 58 self._outfile = gammalib.GFilename()
59 self._ra = None 59 self._ra = None
60 self._dec = None 60 self._dec = None
@@ -65,7 +65,7 @@ class cssens(ctools.cscript): @@ -65,7 +65,7 @@ class cssens(ctools.cscript):
65 self._enumbins = 0 65 self._enumbins = 0
66 self._npix = 200 66 self._npix = 200
67 self._binsz = 0.05 67 self._binsz = 0.05
68 - self._type = "Differential" 68 + self._type = 'Differential'
69 self._ts_thres = 25.0 69 self._ts_thres = 25.0
70 self._max_iter = 50 70 self._max_iter = 50
71 self._num_avg = 3 71 self._num_avg = 3
@@ -83,7 +83,7 @@ class cssens(ctools.cscript): @@ -83,7 +83,7 @@ class cssens(ctools.cscript):
83 # Private methods 83 # Private methods
84 def _get_parameters(self): 84 def _get_parameters(self):
85 """ 85 """
86 - Get user parameters from parfile. 86 + Get user parameters from parfile
87 """ 87 """
88 # Set observation if not done before 88 # Set observation if not done before
89 if self._obs == None or self._obs.size() == 0: 89 if self._obs == None or self._obs.size() == 0:
@@ -91,74 +91,80 @@ class cssens(ctools.cscript): @@ -91,74 +91,80 @@ class cssens(ctools.cscript):
91 91
92 # Set models if we have none 92 # Set models if we have none
93 if self._obs.models().size() == 0: 93 if self._obs.models().size() == 0:
94 - self._obs.models(self["inmodel"].filename()) 94 + self._obs.models(self['inmodel'].filename())
95 95
96 # Get source name 96 # Get source name
97 - self._srcname = self["srcname"].string() 97 + self._srcname = self['srcname'].string()
98 98
99 # Read further parameters 99 # Read further parameters
100 - self._outfile = self["outfile"].filename()  
101 - self._emin = self["emin"].real()  
102 - self._emax = self["emax"].real()  
103 - self._bins = self["bins"].integer() 100 + self._outfile = self['outfile'].filename()
  101 + self._emin = self['emin'].real()
  102 + self._emax = self['emax'].real()
  103 + self._bins = self['bins'].integer()
104 104
105 # Read parameters for binned if requested 105 # Read parameters for binned if requested
106 - self._enumbins = self["enumbins"].integer() 106 + self._enumbins = self['enumbins'].integer()
107 if not self._enumbins == 0: 107 if not self._enumbins == 0:
108 - self._npix = self["npix"].integer()  
109 - self._binsz = self["binsz"].real() 108 + self._npix = self['npix'].integer()
  109 + self._binsz = self['binsz'].real()
110 110
111 # Read remaining parameters 111 # Read remaining parameters
112 - self._edisp = self["edisp"].boolean()  
113 - self._ts_thres = self["sigma"].real()*self["sigma"].real()  
114 - self._max_iter = self["max_iter"].integer()  
115 - self._num_avg = self["num_avg"].integer()  
116 - self._type = self["type"].string() 112 + self._edisp = self['edisp'].boolean()
  113 + self._ts_thres = self['sigma'].real() * self['sigma'].real()
  114 + self._max_iter = self['max_iter'].integer()
  115 + self._num_avg = self['num_avg'].integer()
  116 + self._type = self['type'].string()
117 117
118 # Set some fixed parameters 118 # Set some fixed parameters
119 - self._debug = self["debug"].boolean() # Debugging in client tools 119 + self._debug = self['debug'].boolean() # Debugging in client tools
120 120
121 # Derive some parameters 121 # Derive some parameters
122 self._ebounds = gammalib.GEbounds(self._bins, 122 self._ebounds = gammalib.GEbounds(self._bins,
123 - gammalib.GEnergy(self._emin, "TeV"),  
124 - gammalib.GEnergy(self._emax, "TeV")) 123 + gammalib.GEnergy(self._emin, 'TeV'),
  124 + gammalib.GEnergy(self._emax, 'TeV'))
125 125
126 # Write input parameters into logger 126 # Write input parameters into logger
127 if self._logTerse(): 127 if self._logTerse():
128 self._log_parameters() 128 self._log_parameters()
129 - self._log("\n") 129 + self._log('\n')
130 130
131 # Return 131 # Return
132 return 132 return
133 133
134 def _set_obs(self, lpnt=0.0, bpnt=0.0, emin=0.1, emax=100.0): 134 def _set_obs(self, lpnt=0.0, bpnt=0.0, emin=0.1, emax=100.0):
135 """ 135 """
136 - Set an observation container.  
137 -  
138 - Kwargs:  
139 - lpnt: Galactic longitude of pointing [deg] (default: 0.0)  
140 - bpnt: Galactic latitude of pointing [deg] (default: 0.0)  
141 - emin: Minimum energy [TeV] (default: 0.1)  
142 - emax: Maximum energy [TeV] (default: 100.0)  
143 -  
144 - Returns:  
145 - Observation container. 136 + Set an observation container
  137 +
  138 + Parameters
  139 + ----------
  140 + lpnt : float, optional
  141 + Galactic longitude of pointing (deg)
  142 + bpnt : float, optional
  143 + Galactic latitude of pointing (deg)
  144 + emin : float, optional
  145 + Minimum energy (TeV)
  146 + emax : float, optional
  147 + Maximum energy (TeV)
  148 +
  149 + Returns
  150 + -------
  151 + obs : `~gammalib.GObservations`
  152 + Observation container
146 """ 153 """
147 - # If an observation was provided on input then load it from XML  
148 - # file  
149 - filename = self["inobs"].filename()  
150 - if filename != "NONE" and filename != "": 154 + # If an observation was provided on input then load it from XML file
  155 + filename = self['inobs'].filename()
  156 + if filename != 'NONE' and filename != '':
151 obs = self._get_observations() 157 obs = self._get_observations()
152 158
153 # ... otherwise allocate a single observation 159 # ... otherwise allocate a single observation
154 else: 160 else:
155 161
156 # Read relevant user parameters 162 # Read relevant user parameters
157 - caldb = self["caldb"].string()  
158 - irf = self["irf"].string()  
159 - deadc = self["deadc"].real()  
160 - duration = self["duration"].real()  
161 - rad = self["rad"].real() 163 + caldb = self['caldb'].string()
  164 + irf = self['irf'].string()
  165 + deadc = self['deadc'].real()
  166 + duration = self['duration'].real()
  167 + rad = self['rad'].real()
162 168
163 # Allocate observation container 169 # Allocate observation container
164 obs = gammalib.GObservations() 170 obs = gammalib.GObservations()
@@ -176,7 +182,7 @@ class cssens(ctools.cscript): @@ -176,7 +182,7 @@ class cssens(ctools.cscript):
176 obs.append(run) 182 obs.append(run)
177 183
178 # Set source position 184 # Set source position
179 - offset = self["offset"].real() 185 + offset = self['offset'].real()
180 pntdir.lb_deg(lpnt, bpnt+offset) 186 pntdir.lb_deg(lpnt, bpnt+offset)
181 self._ra = pntdir.ra_deg() 187 self._ra = pntdir.ra_deg()
182 self._dec = pntdir.dec_deg() 188 self._dec = pntdir.dec_deg()
@@ -186,14 +192,19 @@ class cssens(ctools.cscript): @@ -186,14 +192,19 @@ class cssens(ctools.cscript):
186 192
187 def _set_models(self, fitspat=False, fitspec=False): 193 def _set_models(self, fitspat=False, fitspec=False):
188 """ 194 """
189 - Set full model and background model.  
190 -  
191 - Kwargs:  
192 - fitspat: Fit spatial parameter (default: False).  
193 - fitspec: Fit spectral parameters (default: False).  
194 -  
195 - Returns:  
196 - Tuple containing full model and background model. 195 + Set full model and background model
  196 +
  197 + Parameters
  198 + ----------
  199 + fitspat : bool, optional
  200 + Fit spatial parameter?
  201 + fitspec : bool, optional
  202 + Fit spectral parameters?
  203 +
  204 + Returns
  205 + -------
  206 + full_model, bkg_model : `~gammalib.GModels`
  207 + Full model and background model containers
197 """ 208 """
198 # Retrieve full model from observation container 209 # Retrieve full model from observation container
199 full_model = self._obs.models().copy() 210 full_model = self._obs.models().copy()
@@ -201,54 +212,40 @@ class cssens(ctools.cscript): @@ -201,54 +212,40 @@ class cssens(ctools.cscript):
201 # Get source model 212 # Get source model
202 model = full_model[self._srcname] 213 model = full_model[self._srcname]
203 214
204 - # Check that model has a Prefactor  
205 - if not model.has_par("Prefactor"):  
206 - msg = "Model \""+self._srcname+"\" has no parameter "+\  
207 - "\"Prefactor\". Only spectral models with a "+\  
208 - "\"Prefactor\" parameter are supported." 215 + # If source model has no "Prefactor" parameter then raise an exception
  216 + if not model.has_par('Prefactor'):
  217 + msg = ('Model "%s" has no parameter "Prefactor". Only spectral '
  218 + 'models with a "Prefactor" parameter are supported.' %
  219 + self._srcname)
209 raise RuntimeError(msg) 220 raise RuntimeError(msg)
210 221
211 # Set source position 222 # Set source position
212 if self._ra != None and self._dec != None: 223 if self._ra != None and self._dec != None:
213 - if model.has_par("RA") and model.has_par("DEC"):  
214 - model["RA"].value(self._ra)  
215 - model["DEC"].value(self._dec) 224 + if model.has_par('RA') and model.has_par('DEC'):
  225 + model['RA'].value(self._ra)
  226 + model['DEC'].value(self._dec)
  227 +
  228 + # Set possible spatial and spectral parameters
  229 + spatial = ['RA', 'DEC', 'Sigma', 'Radius', 'Width', 'PA',
  230 + 'MinorRadius', 'MajorRadius']
  231 + spectral = ['Index', 'Index1', 'Index2', 'BreakEnergy', 'CutoffEnergy',
  232 + 'InverseCutoffEnergy']
216 233
217 # Fit or fix spatial parameters 234 # Fit or fix spatial parameters
218 - if fitspat:  
219 - if model.has_par("RA"):  
220 - model["RA"].free()  
221 - if model.has_par("DEC"):  
222 - model["DEC"].free()  
223 - if model.has_par("Sigma"):  
224 - model["Sigma"].free()  
225 - if model.has_par("Radius"):  
226 - model["Radius"].free()  
227 - if model.has_par("Width"):  
228 - model["Width"].free()  
229 - else:  
230 - if model.has_par("RA"):  
231 - model["RA"].fix()  
232 - if model.has_par("DEC"):  
233 - model["DEC"].fix()  
234 - if model.has_par("Sigma"):  
235 - model["Sigma"].fix()  
236 - if model.has_par("Radius"):  
237 - model["Radius"].fix()  
238 - if model.has_par("Width"):  
239 - model["Width"].fix() 235 + for par in spatial:
  236 + if model.has_par(par):
  237 + if fitspat:
  238 + model[par].free()
  239 + else:
  240 + model[par].fix()
240 241
241 # Fit or fix spectral parameters 242 # Fit or fix spectral parameters
242 - if fitspec:  
243 - if model.has_par("Index"):  
244 - model["Index"].free()  
245 - if model.has_par("Cutoff"):  
246 - model["Cutoff"].free()  
247 - else:  
248 - if model.has_par("Index"):  
249 - model["Index"].fix()  
250 - if model.has_par("Cutoff"):  
251 - model["Cutoff"].fix() 243 + for par in spectral:
  244 + if model.has_par(par):
  245 + if fitspec:
  246 + model[par].free()
  247 + else:
  248 + model[par].fix()
252 249
253 # Create background model 250 # Create background model
254 bkg_model = full_model.copy() 251 bkg_model = full_model.copy()
@@ -259,14 +256,14 @@ class cssens(ctools.cscript): @@ -259,14 +256,14 @@ class cssens(ctools.cscript):
259 256
260 def _set_obs_ebounds(self, emin, emax): 257 def _set_obs_ebounds(self, emin, emax):
261 """ 258 """
262 - Set energy boundaries for observation container.  
263 -  
264 - Sets the energy boundaries for all observations in the observation  
265 - container.  
266 -  
267 - Args:  
268 - emin: Minimum energy  
269 - emax: Maximum energy 259 + Set energy boundaries for all observations in container
  260 +
  261 + Parameters
  262 + ----------
  263 + emin : `~gammalib.GEnergy`
  264 + Minimum energy
  265 + emax : `~gammalib.GEnergy`
  266 + Maximum energy
270 """ 267 """
271 # Loop over all observations in container 268 # Loop over all observations in container
272 for obs in self._obs: 269 for obs in self._obs:
@@ -274,38 +271,34 @@ class cssens(ctools.cscript): @@ -274,38 +271,34 @@ class cssens(ctools.cscript):
274 # Get observation energy boundaries 271 # Get observation energy boundaries
275 obs_ebounds = obs.events().ebounds() 272 obs_ebounds = obs.events().ebounds()
276 273
277 - # Get minimum and maximum energy 274 + # Get minimum and maximum energy of the observation
278 obs_emin = obs_ebounds.emin() 275 obs_emin = obs_ebounds.emin()
279 obs_emax = obs_ebounds.emax() 276 obs_emax = obs_ebounds.emax()
280 277
281 - # Case A: bin fully contained in observation ebounds 278 + # If [emin,emax] is fully contained in the observation energy range
  279 + # the use [emin,emax] as energy boundaries
282 if obs_emin <= emin and obs_emax >= emax: 280 if obs_emin <= emin and obs_emax >= emax:
283 ebounds = gammalib.GEbounds(emin, emax) 281 ebounds = gammalib.GEbounds(emin, emax)
284 282
285 - # Case B: bin fully outside of obs ebounds 283 + # ... otherwise, if [emin,emax] is completely outside the
  284 + # observation energy range then set the energy boundaries to the
  285 + # zero-width interval [0,0]
286 elif emax < obs_emin or emin > obs_emax: 286 elif emax < obs_emin or emin > obs_emax:
287 -  
288 - # Set zero range (inspired by ctselect)  
289 - e0 = gammalib.GEnergy(0.0, "TeV") 287 + e0 = gammalib.GEnergy(0.0, 'TeV')
290 ebounds = gammalib.GEbounds(e0, e0) 288 ebounds = gammalib.GEbounds(e0, e0)
291 289
292 - # Case C: bin partly overlapping with observation ebounds 290 + # ... otherwise, if [emin,emax] overlaps partially with the
  291 + # observation energy range then set the energy boundaries to the
  292 + # overlapping part
293 else: 293 else:
  294 + # Set overlapping energy range
  295 + set_emin = max(emin, obs_emin)
  296 + set_emax = min(emax, obs_emax)
294 297
295 - # Set energy range as obs ebounds were fully contained inside energy bin  
296 - set_emin = emin  
297 - set_emax = emax  
298 -  
299 - # Adjust energy bin to respect observation energy boundary  
300 - if emin < obs_emin:  
301 - set_emin = obs_emin  
302 - if emax > obs_emax:  
303 - set_emax = obs_emax  
304 -  
305 - #Set energy boundaries 298 + # Set energy boundaries
306 ebounds = gammalib.GEbounds(set_emin, set_emax) 299 ebounds = gammalib.GEbounds(set_emin, set_emax)
307 300
308 - # Set energy boundaries 301 + # Set the energy boundaries as the boundaries of the observation
309 obs.events().ebounds(ebounds) 302 obs.events().ebounds(ebounds)
310 303
311 # Return 304 # Return
@@ -313,36 +306,48 @@ class cssens(ctools.cscript): @@ -313,36 +306,48 @@ class cssens(ctools.cscript):
313 306
314 def _get_crab_flux(self, emin, emax): 307 def _get_crab_flux(self, emin, emax):
315 """ 308 """
316 - Return Crab photon flux in a given energy interval.  
317 -  
318 - Args:  
319 - emin: Minimum energy  
320 - emax: Maximum energy  
321 -  
322 - Returns:  
323 - Crab photon flux in specified energy interval. 309 + Return Crab photon flux in a given energy interval
  310 +
  311 + Parameters
  312 + ----------
  313 + emin : `~gammalib.GEnergy`
  314 + Minimum energy
  315 + emax : `~gammalib.GEnergy`
  316 + Maximum energy
  317 +
  318 + Returns
  319 + -------
  320 + flux : float
  321 + Crab photon flux in specified energy interval (ph/cm2/s)
324 """ 322 """
325 # Set Crab TeV spectral model based on a power law 323 # Set Crab TeV spectral model based on a power law
326 crab = gammalib.GModelSpectralPlaw(5.7e-16, -2.48, 324 crab = gammalib.GModelSpectralPlaw(5.7e-16, -2.48,
327 - gammalib.GEnergy(0.3, "TeV")) 325 + gammalib.GEnergy(0.3, 'TeV'))
328 326
329 - # Determine flux 327 + # Determine photon flux
330 flux = crab.flux(emin, emax) 328 flux = crab.flux(emin, emax)
331 329
332 - # Return flux 330 + # Return photon flux
333 return flux 331 return flux
334 332
335 def _get_sensitivity(self, emin, emax, bkg_model, full_model): 333 def _get_sensitivity(self, emin, emax, bkg_model, full_model):
336 """ 334 """
337 Determine sensitivity for given observations. 335 Determine sensitivity for given observations.
338 336
339 - Args:  
340 - emin: Minimum energy for fitting and flux computation  
341 - emax: Maximum energy for fitting and flux computation  
342 - bkg_model: Background model  
343 - full_model: Source model  
344 -  
345 - Returns: 337 + Parameters
  338 + ----------
  339 + emin : `~gammalib.GEnergy`
  340 + Minimum energy for fitting and flux computation
  341 + emax : `~gammalib.GEnergy`
  342 + Maximum energy for fitting and flux computation
  343 + bkg_model : `~gammalib.GModels`
  344 + Background model
  345 + full_model : `~gammalib.GModels`
  346 + Source plus background model
  347 +
  348 + Returns
  349 + -------
  350 + result : dict
346 Result dictionary 351 Result dictionary
347 """ 352 """
348 # Set TeV->erg conversion factor 353 # Set TeV->erg conversion factor
@@ -364,17 +369,17 @@ class cssens(ctools.cscript): @@ -364,17 +369,17 @@ class cssens(ctools.cscript):
364 369
365 # Write header 370 # Write header
366 if self._logTerse(): 371 if self._logTerse():
367 - self._log("\n")  
368 - self._log.header2("Energies: "+str(emin)+" - "+str(emax))  
369 - self._log.parformat("Crab flux") 372 + self._log('\n')
  373 + self._log.header2('Energies: '+str(emin)+' - '+str(emax))
  374 + self._log.parformat('Crab flux')
370 self._log(crab_flux) 375 self._log(crab_flux)
371 - self._log(" ph/cm2/s\n")  
372 - self._log.parformat("Source model flux") 376 + self._log(' ph/cm2/s\n')
  377 + self._log.parformat('Source model flux')
373 self._log(src_flux) 378 self._log(src_flux)
374 - self._log(" ph/cm2/s\n")  
375 - self._log.parformat("Crab unit factor") 379 + self._log(' ph/cm2/s\n')
  380 + self._log.parformat('Crab unit factor')
376 self._log(crab_unit) 381 self._log(crab_unit)
377 - self._log("\n") 382 + self._log('\n')
378 383
379 # Initialise loop 384 # Initialise loop
380 crab_flux_value = [] 385 crab_flux_value = []
@@ -392,7 +397,7 @@ class cssens(ctools.cscript): @@ -392,7 +397,7 @@ class cssens(ctools.cscript):
392 397
393 # Write header 398 # Write header
394 if self._logExplicit(): 399 if self._logExplicit():
395 - self._log.header2("Iteration "+str(iterations)) 400 + self._log.header2('Iteration '+str(iterations))
396 401
397 # Set source model. crab_prefactor is the Prefactor that 402 # Set source model. crab_prefactor is the Prefactor that
398 # corresponds to 1 Crab 403 # corresponds to 1 Crab
@@ -416,10 +421,10 @@ class cssens(ctools.cscript): @@ -416,10 +421,10 @@ class cssens(ctools.cscript):
416 421
417 # Write simulation results 422 # Write simulation results
418 if self._logExplicit(): 423 if self._logExplicit():
419 - self._log.header3("Simulation")  
420 - self._log.parformat("Number of simulated events") 424 + self._log.header3('Simulation')
  425 + self._log.parformat('Number of simulated events')
421 self._log(nevents) 426 self._log(nevents)
422 - self._log("\n") 427 + self._log('\n')
423 428
424 # Fit background only 429 # Fit background only
425 sim.models(bkg_model) 430 sim.models(bkg_model)
@@ -434,25 +439,25 @@ class cssens(ctools.cscript): @@ -434,25 +439,25 @@ class cssens(ctools.cscript):
434 439
435 # Write background fit results 440 # Write background fit results
436 if self._logExplicit(): 441 if self._logExplicit():
437 - self._log.header3("Background model fit")  
438 - self._log.parformat("log likelihood") 442 + self._log.header3('Background model fit')
  443 + self._log.parformat('log likelihood')
439 self._log(LogL_bgm) 444 self._log(LogL_bgm)
440 - self._log("\n")  
441 - self._log.parformat("Number of predicted events") 445 + self._log('\n')
  446 + self._log.parformat('Number of predicted events')
442 self._log(npred_bgm) 447 self._log(npred_bgm)
443 - self._log("\n")  
444 - self._log.parformat("Fit quality") 448 + self._log('\n')
  449 + self._log.parformat('Fit quality')
445 self._log(quality_bgm) 450 self._log(quality_bgm)
446 - self._log("\n") 451 + self._log('\n')
447 452
448 # Write model fit results 453 # Write model fit results
449 if self._logExplicit(): 454 if self._logExplicit():
450 for model in result_bgm: 455 for model in result_bgm:
451 - self._log.parformat("Model") 456 + self._log.parformat('Model')
452 self._log(model.name()) 457 self._log(model.name())
453 - self._log("\n") 458 + self._log('\n')
454 for par in model: 459 for par in model:
455 - self._log(str(par)+"\n") 460 + self._log(str(par)+'\n')
456 461
457 # Fit background and test source 462 # Fit background and test source
458 sim.models(src_model) 463 sim.models(src_model)
@@ -468,31 +473,31 @@ class cssens(ctools.cscript): @@ -468,31 +473,31 @@ class cssens(ctools.cscript):
468 473
469 # Write background and test source fit results 474 # Write background and test source fit results
470 if self._logExplicit(): 475 if self._logExplicit():
471 - self._log.header3("Background and test source model fit")  
472 - self._log.parformat("Test statistics") 476 + self._log.header3('Background and test source model fit')
  477 + self._log.parformat('Test statistics')
473 self._log(ts) 478 self._log(ts)
474 - self._log("\n")  
475 - self._log.parformat("log likelihood") 479 + self._log('\n')
  480 + self._log.parformat('log likelihood')
476 self._log(LogL_all) 481 self._log(LogL_all)
477 - self._log("\n")  
478 - self._log.parformat("Number of predicted events") 482 + self._log('\n')
  483 + self._log.parformat('Number of predicted events')
479 self._log(npred_all) 484 self._log(npred_all)
480 - self._log("\n")  
481 - self._log.parformat("Fit quality") 485 + self._log('\n')
  486 + self._log.parformat('Fit quality')
482 self._log(quality_all) 487 self._log(quality_all)
483 - self._log("\n") 488 + self._log('\n')
484 # 489 #
485 for model in result_all: 490 for model in result_all:
486 - self._log.parformat("Model") 491 + self._log.parformat('Model')
487 self._log(model.name()) 492 self._log(model.name())
488 - self._log("\n") 493 + self._log('\n')
489 for par in model: 494 for par in model:
490 - self._log(str(par)+"\n") 495 + self._log(str(par)+'\n')
491 496
492 # Start over if TS was non-positive 497 # Start over if TS was non-positive
493 if ts <= 0.0: 498 if ts <= 0.0:
494 if self._logExplicit(): 499 if self._logExplicit():
495 - self._log("Non positive TS. Start over.\n") 500 + self._log('Non positive TS. Start over.\n')
496 continue 501 continue
497 502
498 # Get fitted Crab, photon and energy fluxes 503 # Get fitted Crab, photon and energy fluxes
@@ -504,7 +509,7 @@ class cssens(ctools.cscript): @@ -504,7 +509,7 @@ class cssens(ctools.cscript):
504 energy_flux = result_all[self._srcname].spectral().eflux(emin, emax) 509 energy_flux = result_all[self._srcname].spectral().eflux(emin, emax)
505 510
506 # Compute differential sensitivity in unit erg/cm2/s 511 # Compute differential sensitivity in unit erg/cm2/s
507 - energy = gammalib.GEnergy(e_mean, "TeV") 512 + energy = gammalib.GEnergy(e_mean, 'TeV')
508 time = gammalib.GTime() 513 time = gammalib.GTime()
509 sensitivity = result_all[self._srcname].spectral().eval(energy, time) * \ 514 sensitivity = result_all[self._srcname].spectral().eval(energy, time) * \
510 e_mean*erg_mean * 1.0e6 515 e_mean*erg_mean * 1.0e6
@@ -526,40 +531,40 @@ class cssens(ctools.cscript): @@ -526,40 +531,40 @@ class cssens(ctools.cscript):
526 531
527 # Write background and test source fit results 532 # Write background and test source fit results
528 if self._logExplicit(): 533 if self._logExplicit():
529 - self._log.parformat("Photon flux") 534 + self._log.parformat('Photon flux')
530 self._log(photon_flux) 535 self._log(photon_flux)
531 - self._log(" ph/cm2/s\n")  
532 - self._log.parformat("Energy flux") 536 + self._log(' ph/cm2/s\n')
  537 + self._log.parformat('Energy flux')
533 self._log(energy_flux) 538 self._log(energy_flux)
534 - self._log(" erg/cm2/s\n")  
535 - self._log.parformat("Crab flux") 539 + self._log(' erg/cm2/s\n')
  540 + self._log.parformat('Crab flux')
536 self._log(crab_flux*1000.0) 541 self._log(crab_flux*1000.0)
537 - self._log(" mCrab\n")  
538 - self._log.parformat("Differential sensitivity") 542 + self._log(' mCrab\n')
  543 + self._log.parformat('Differential sensitivity')
539 self._log(sensitivity) 544 self._log(sensitivity)
540 - self._log(" erg/cm2/s\n") 545 + self._log(' erg/cm2/s\n')
541 for model in result_all: 546 for model in result_all:
542 - self._log.parformat("Model") 547 + self._log.parformat('Model')
543 self._log(model.name()) 548 self._log(model.name())
544 - self._log("\n") 549 + self._log('\n')
545 for par in model: 550 for par in model:
546 - self._log(str(par)+"\n") 551 + self._log(str(par)+'\n')
547 elif self._logTerse(): 552 elif self._logTerse():
548 - self._log.parformat("Iteration "+str(iterations))  
549 - self._log("TS=") 553 + self._log.parformat('Iteration '+str(iterations))
  554 + self._log('TS=')
550 self._log(ts) 555 self._log(ts)
551 - self._log(" ")  
552 - self._log("corr=") 556 + self._log(' ')
  557 + self._log('corr=')
553 self._log(correct) 558 self._log(correct)
554 - self._log(" ") 559 + self._log(' ')
555 self._log(photon_flux) 560 self._log(photon_flux)
556 - self._log(" ph/cm2/s = ") 561 + self._log(' ph/cm2/s = ')
557 self._log(energy_flux) 562 self._log(energy_flux)
558 - self._log(" erg/cm2/s = ") 563 + self._log(' erg/cm2/s = ')
559 self._log(crab_flux*1000.0) 564 self._log(crab_flux*1000.0)
560 - self._log(" mCrab = ") 565 + self._log(' mCrab = ')
561 self._log(sensitivity) 566 self._log(sensitivity)
562 - self._log(" erg/cm2/s\n") 567 + self._log(' erg/cm2/s\n')
563 568
564 # Compute sliding average of extrapolated fitted prefactor, 569 # Compute sliding average of extrapolated fitted prefactor,
565 # photon and energy flux. This damps out fluctuations and 570 # photon and energy flux. This damps out fluctuations and
@@ -593,11 +598,11 @@ class cssens(ctools.cscript): @@ -593,11 +598,11 @@ class cssens(ctools.cscript):
593 if ratio >= 0.99 and ratio <= 1.01 and \ 598 if ratio >= 0.99 and ratio <= 1.01 and \
594 correct >= 0.9 and correct <= 1.1: 599 correct >= 0.9 and correct <= 1.1:
595 if self._logTerse(): 600 if self._logTerse():
596 - self._log(" Converged ("+str(ratio)+")\n") 601 + self._log(' Converged ('+str(ratio)+')\n')
597 break 602 break
598 else: 603 else:
599 if self._logTerse(): 604 if self._logTerse():
600 - self._log(" Flux is zero.\n") 605 + self._log(' Flux is zero.\n')
601 break 606 break
602 607
603 # Use average for next iteration 608 # Use average for next iteration
@@ -606,57 +611,57 @@ class cssens(ctools.cscript): @@ -606,57 +611,57 @@ class cssens(ctools.cscript):
606 # Exit loop if number of trials exhausted 611 # Exit loop if number of trials exhausted
607 if (iterations >= self._max_iter): 612 if (iterations >= self._max_iter):
608 if self._logTerse(): 613 if self._logTerse():
609 - self._log(" Test ended after "+str(self._max_iter)+  
610 - " iterations.\n") 614 + self._log(' Test ended after '+str(self._max_iter)+
  615 + ' iterations.\n')
611 break 616 break
612 617
613 # Write fit results 618 # Write fit results
614 if self._logTerse(): 619 if self._logTerse():
615 - self._log.header3("Fit results")  
616 - self._log.parformat("Test statistics") 620 + self._log.header3('Fit results')
  621 + self._log.parformat('Test statistics')
617 self._log(ts) 622 self._log(ts)
618 - self._log("\n")  
619 - self._log.parformat("Photon flux") 623 + self._log('\n')
  624 + self._log.parformat('Photon flux')
620 self._log(photon_flux) 625 self._log(photon_flux)
621 - self._log(" ph/cm2/s\n")  
622 - self._log.parformat("Energy flux") 626 + self._log(' ph/cm2/s\n')
  627 + self._log.parformat('Energy flux')
623 self._log(energy_flux) 628 self._log(energy_flux)
624 - self._log(" erg/cm2/s\n")  
625 - self._log.parformat("Crab flux") 629 + self._log(' erg/cm2/s\n')
  630 + self._log.parformat('Crab flux')
626 self._log(crab_flux*1000.0) 631 self._log(crab_flux*1000.0)
627 - self._log(" mCrab\n")  
628 - self._log.parformat("Differential sensitivity") 632 + self._log(' mCrab\n')
  633 + self._log.parformat('Differential sensitivity')
629 self._log(sensitivity) 634 self._log(sensitivity)
630 - self._log(" erg/cm2/s\n")  
631 - self._log.parformat("Number of simulated events") 635 + self._log(' erg/cm2/s\n')
  636 + self._log.parformat('Number of simulated events')
632 self._log(nevents) 637 self._log(nevents)
633 - self._log("\n")  
634 - self._log.header3("Background and test source model fitting")  
635 - self._log.parformat("log likelihood") 638 + self._log('\n')
  639 + self._log.header3('Background and test source model fitting')
  640 + self._log.parformat('log likelihood')
636 self._log(LogL_all) 641 self._log(LogL_all)
637 - self._log("\n")  
638 - self._log.parformat("Number of predicted events") 642 + self._log('\n')
  643 + self._log.parformat('Number of predicted events')
639 self._log(npred_all) 644 self._log(npred_all)
640 - self._log("\n") 645 + self._log('\n')
641 for model in result_all: 646 for model in result_all:
642 - self._log.parformat("Model") 647 + self._log.parformat('Model')
643 self._log(model.name()) 648 self._log(model.name())
644 - self._log("\n") 649 + self._log('\n')
645 for par in model: 650 for par in model:
646 - self._log(str(par)+"\n")  
647 - self._log.header3("Background model fit")  
648 - self._log.parformat("log likelihood") 651 + self._log(str(par)+'\n')
  652 + self._log.header3('Background model fit')
  653 + self._log.parformat('log likelihood')
649 self._log(LogL_bgm) 654 self._log(LogL_bgm)
650 - self._log("\n")  
651 - self._log.parformat("Number of predicted events") 655 + self._log('\n')
  656 + self._log.parformat('Number of predicted events')
652 self._log(npred_bgm) 657 self._log(npred_bgm)
653 - self._log("\n") 658 + self._log('\n')
654 for model in result_bgm: 659 for model in result_bgm:
655 - self._log.parformat("Model") 660 + self._log.parformat('Model')
656 self._log(model.name()) 661 self._log(model.name())
657 - self._log("\n") 662 + self._log('\n')
658 for par in model: 663 for par in model:
659 - self._log(str(par)+"\n") 664 + self._log(str(par)+'\n')
660 665
661 # Restore energy boundaries of observation container 666 # Restore energy boundaries of observation container
662 for i, obs in enumerate(self._obs): 667 for i, obs in enumerate(self._obs):
@@ -675,7 +680,7 @@ class cssens(ctools.cscript): @@ -675,7 +680,7 @@ class cssens(ctools.cscript):
675 # Public methods 680 # Public methods
676 def run(self): 681 def run(self):
677 """ 682 """
678 - Run the script. 683 + Run the script
679 """ 684 """
680 # Switch screen logging on in debug mode 685 # Switch screen logging on in debug mode
681 if self._logDebug(): 686 if self._logDebug():
@@ -698,42 +703,40 @@ class cssens(ctools.cscript): @@ -698,42 +703,40 @@ class cssens(ctools.cscript):
698 703
699 # Write models into logger 704 # Write models into logger
700 if self._logTerse(): 705 if self._logTerse():
701 - self._log("\n")  
702 - self._log.header1("Models")  
703 - self._log.header2("Background model") 706 + self._log('\n')
  707 + self._log.header1('Models')
  708 + self._log.header2('Background model')
704 self._log(str(bkg_model)) 709 self._log(str(bkg_model))
705 - self._log("\n\n")  
706 - self._log.header2("Full model") 710 + self._log('\n\n')
  711 + self._log.header2('Full model')
707 self._log(str(full_model)) 712 self._log(str(full_model))
708 - self._log("\n") 713 + self._log('\n')
709 714
710 # Write header 715 # Write header
711 if self._logTerse(): 716 if self._logTerse():
712 - self._log("\n")  
713 - self._log.header1("Sensitivity determination")  
714 - self._log.parformat("Type") 717 + self._log('\n')
  718 + self._log.header1('Sensitivity determination')
  719 + self._log.parformat('Type')
715 self._log(self._type) 720 self._log(self._type)
716 - self._log("\n") 721 + self._log('\n')
717 722
718 # Loop over energy bins 723 # Loop over energy bins
719 for ieng in range(self._ebounds.size()): 724 for ieng in range(self._ebounds.size()):
720 725
721 # Set energies 726 # Set energies
722 - if self._type == "Differential": 727 + if self._type == 'Differential':
723 emin = self._ebounds.emin(ieng) 728 emin = self._ebounds.emin(ieng)
724 emax = self._ebounds.emax(ieng) 729 emax = self._ebounds.emax(ieng)
725 - elif self._type == "Integral": 730 + elif self._type == 'Integral':
726 emin = self._ebounds.emin(ieng) 731 emin = self._ebounds.emin(ieng)
727 emax = self._ebounds.emax() 732 emax = self._ebounds.emax()
728 else: 733 else:
729 - msg = "Invalid sensitivity type \""+self._type+"\" "+\  
730 - "encountered. Either specify \"Differential\" or "+\  
731 - "\"Integral\"." 734 + msg = ('Invalid sensitivity type "%s" encountered. Either '
  735 + 'specify "Differential" or "Integral".' % self._type)
732 raise RuntimeError(msg) 736 raise RuntimeError(msg)
733 737
734 # Determine sensitivity 738 # Determine sensitivity
735 - result = self._get_sensitivity(emin, emax,  
736 - bkg_model, full_model) 739 + result = self._get_sensitivity(emin, emax, bkg_model, full_model)
737 740
738 # Write results 741 # Write results
739 if ieng == 0: 742 if ieng == 0:
@@ -759,7 +762,7 @@ class cssens(ctools.cscript): @@ -759,7 +762,7 @@ class cssens(ctools.cscript):
759 762
760 def execute(self): 763 def execute(self):
761 """ 764 """
762 - Execute the script. 765 + Execute the script
763 """ 766 """
764 # Open logfile 767 # Open logfile
765 self.logFileOpen() 768 self.logFileOpen()
test/test_cssens.py
  View file @003f01c
@@ -121,7 +121,7 @@ class Test(test): @@ -121,7 +121,7 @@ class Test(test):
121 sens['emax'] = 10.0 121 sens['emax'] = 10.0
122 sens['bins'] = 1 122 sens['bins'] = 1
123 sens['logfile'] = 'cssens_py1.log' 123 sens['logfile'] = 'cssens_py1.log'
124 - sens['chatter'] = 2 124 + sens['chatter'] = 4
125 125
126 # Run cssens script 126 # Run cssens script
127 sens.logFileOpen() # Make sure we get a log file 127 sens.logFileOpen() # Make sure we get a log file