Commit 7955dfc9afa801c898f446c43b68b0089e0be98d

Authored by Jürgen Knödlseder
1 parent 50d6ad88

Integrate code

scripts/csspec.py
@@ -23,9 +23,9 @@ import ctools @@ -23,9 +23,9 @@ import ctools
23 import sys 23 import sys
24 import tempfile 24 import tempfile
25 25
26 -# ============== #  
27 -# cstsdist class #  
28 -# ============== # 26 +# ============ #
  27 +# csspec class #
  28 +# ============ #
29 class csspec(ctools.cscript): 29 class csspec(ctools.cscript):
30 """ 30 """
31 This class implements the creation of spectral points. It derives from 31 This class implements the creation of spectral points. It derives from
@@ -43,14 +43,14 @@ class csspec(ctools.cscript): @@ -43,14 +43,14 @@ class csspec(ctools.cscript):
43 self.version = "1.0.0" 43 self.version = "1.0.0"
44 44
45 # Initialise some members 45 # Initialise some members
46 - self.obs = None 46 + self.obs = None
47 47
48 # Initialise some members 48 # Initialise some members
49 if isinstance(argv[0],gammalib.GObservations): 49 if isinstance(argv[0],gammalib.GObservations):
50 self.obs = argv[0] 50 self.obs = argv[0]
51 - argv = argv[1:] 51 + argv = argv[1:]
52 else: 52 else:
53 - self.obs = gammalib.GObservations() 53 + self.obs = gammalib.GObservations()
54 self.obs.clear() 54 self.obs.clear()
55 self.m_outfile = "" 55 self.m_outfile = ""
56 56
@@ -105,6 +105,7 @@ class csspec(ctools.cscript): @@ -105,6 +105,7 @@ class csspec(ctools.cscript):
105 pars.append(gammalib.GApplicationPar("outfile","f","a","spectrum.fits","","","Output file name")) 105 pars.append(gammalib.GApplicationPar("outfile","f","a","spectrum.fits","","","Output file name"))
106 pars.append(gammalib.GApplicationPar("expcube","f","a","NONE","","","Exposure cube file (only needed for stacked analysis)")) 106 pars.append(gammalib.GApplicationPar("expcube","f","a","NONE","","","Exposure cube file (only needed for stacked analysis)"))
107 pars.append(gammalib.GApplicationPar("psfcube","f","a","NONE","","","PSF cube file (only needed for stacked analysis)")) 107 pars.append(gammalib.GApplicationPar("psfcube","f","a","NONE","","","PSF cube file (only needed for stacked analysis)"))
  108 + pars.append(gammalib.GApplicationPar("bkgcube","s","a","NONE","","","Background cube file (only needed for stacked analysis)"))
108 pars.append(gammalib.GApplicationPar("caldb","s","a","dummy","","","Calibration database")) 109 pars.append(gammalib.GApplicationPar("caldb","s","a","dummy","","","Calibration database"))
109 pars.append(gammalib.GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function")) 110 pars.append(gammalib.GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function"))
110 pars.append(gammalib.GApplicationPar("emin","r","h","0.1","","","Lower energy limit for spectral points(TeV)")) 111 pars.append(gammalib.GApplicationPar("emin","r","h","0.1","","","Lower energy limit for spectral points(TeV)"))
@@ -135,7 +136,6 @@ class csspec(ctools.cscript): @@ -135,7 +136,6 @@ class csspec(ctools.cscript):
135 """ 136 """
136 Get parameters from parfile and setup the observation. 137 Get parameters from parfile and setup the observation.
137 """ 138 """
138 - # Get parameters  
139 # Set observation if not done before 139 # Set observation if not done before
140 if self.obs == None or self.obs.size() == 0: 140 if self.obs == None or self.obs.size() == 0:
141 self.require_inobs("csspec::get_parameters()") 141 self.require_inobs("csspec::get_parameters()")
@@ -154,14 +154,14 @@ class csspec(ctools.cscript): @@ -154,14 +154,14 @@ class csspec(ctools.cscript):
154 # Get binning flag 154 # Get binning flag
155 self.m_binned = self["binned"].boolean() 155 self.m_binned = self["binned"].boolean()
156 if self.m_binned: 156 if self.m_binned:
157 - self.m_xref = self["xref"].real()  
158 - self.m_yref = self["yref"].real()  
159 - self.m_nxpix = self["nxpix"].integer()  
160 - self.m_nypix = self["nypix"].integer()  
161 - self.m_binsz = self["binsz"].real() 157 + self.m_xref = self["xref"].real()
  158 + self.m_yref = self["yref"].real()
  159 + self.m_nxpix = self["nxpix"].integer()
  160 + self.m_nypix = self["nypix"].integer()
  161 + self.m_binsz = self["binsz"].real()
162 self.m_coordsys = self["coordsys"].string() 162 self.m_coordsys = self["coordsys"].string()
163 - self.m_proj = self["proj"].string()  
164 - self.m_ebins = self["nebins"].integer() 163 + self.m_proj = self["proj"].string()
  164 + self.m_ebins = self["nebins"].integer()
165 165
166 # Read other parameters 166 # Read other parameters
167 self.m_outfile = self["outfile"].filename() 167 self.m_outfile = self["outfile"].filename()
@@ -277,8 +277,8 @@ class csspec(ctools.cscript): @@ -277,8 +277,8 @@ class csspec(ctools.cscript):
277 table.extname("SPECTRUM") 277 table.extname("SPECTRUM")
278 278
279 # Add Header for compatibility with gammalib.GMWLSpectrum 279 # Add Header for compatibility with gammalib.GMWLSpectrum
280 - table.card("Instrument", "CTA", "Name of Instrument")  
281 - table.card("Telescope", "CTA", "Name of Telescope") 280 + table.card("INSTRUME", "CTA", "Name of Instrument")
  281 + table.card("TELESCOP", "CTA", "Name of Telescope")
282 282
283 # Create FITS table columns 283 # Create FITS table columns
284 energy = gammalib.GFitsTableDoubleCol("Energy", self.m_ebounds.size()) 284 energy = gammalib.GFitsTableDoubleCol("Energy", self.m_ebounds.size())
@@ -313,7 +313,7 @@ class csspec(ctools.cscript): @@ -313,7 +313,7 @@ class csspec(ctools.cscript):
313 energy[i] = elogmean.TeV() 313 energy[i] = elogmean.TeV()
314 314
315 # Store energy errors 315 # Store energy errors
316 - energy_low[i] = (elogmean - emin).TeV() 316 + energy_low[i] = (elogmean - emin).TeV()
317 energy_high[i] = (emax - elogmean).TeV() 317 energy_high[i] = (emax - elogmean).TeV()
318 318
319 # Log information 319 # Log information
@@ -330,15 +330,18 @@ class csspec(ctools.cscript): @@ -330,15 +330,18 @@ class csspec(ctools.cscript):
330 select["ra"].value("UNDEFINED") 330 select["ra"].value("UNDEFINED")
331 select["dec"].value("UNDEFINED") 331 select["dec"].value("UNDEFINED")
332 select.run() 332 select.run()
333 - 333 +
  334 + # Retrieve observation
334 obs = select.obs() 335 obs = select.obs()
335 - 336 +
  337 + # Binned analysis
336 if self.m_binned: 338 if self.m_binned:
337 - 339 +
  340 + # Header
338 if self.logTerse(): 341 if self.logTerse():
339 self.log.header3("Binning events") 342 self.log.header3("Binning events")
340 343
341 - # bin events 344 + # Bin events
342 bin = ctools.ctbin(select.obs()) 345 bin = ctools.ctbin(select.obs())
343 bin["usepnt"].boolean(False) 346 bin["usepnt"].boolean(False)
344 bin["ebinalg"].string("LOG") 347 bin["ebinalg"].string("LOG")
@@ -354,10 +357,11 @@ class csspec(ctools.cscript): @@ -354,10 +357,11 @@ class csspec(ctools.cscript):
354 bin["proj"].string(self.m_proj) 357 bin["proj"].string(self.m_proj)
355 bin.run() 358 bin.run()
356 359
  360 + # Header
357 if self.logTerse(): 361 if self.logTerse():
358 self.log.header3("Creating exposure cube") 362 self.log.header3("Creating exposure cube")
359 363
360 - # create exposure cube 364 + # Create exposure cube
361 expcube = ctools.ctexpcube(select.obs()) 365 expcube = ctools.ctexpcube(select.obs())
362 expcube["incube"].filename("NONE") 366 expcube["incube"].filename("NONE")
363 expcube["usepnt"].boolean(False) 367 expcube["usepnt"].boolean(False)
@@ -374,10 +378,11 @@ class csspec(ctools.cscript): @@ -374,10 +378,11 @@ class csspec(ctools.cscript):
374 expcube["proj"].string(self.m_proj) 378 expcube["proj"].string(self.m_proj)
375 expcube.run() 379 expcube.run()
376 380
  381 + # Header
377 if self.logTerse(): 382 if self.logTerse():
378 self.log.header3("Creating PSF cube") 383 self.log.header3("Creating PSF cube")
379 384
380 - # create psf cube 385 + # Create psf cube
381 psfcube = ctools.ctpsfcube(select.obs()) 386 psfcube = ctools.ctpsfcube(select.obs())
382 psfcube["incube"].filename("NONE") 387 psfcube["incube"].filename("NONE")
383 psfcube["usepnt"].boolean(False) 388 psfcube["usepnt"].boolean(False)
@@ -394,10 +399,11 @@ class csspec(ctools.cscript): @@ -394,10 +399,11 @@ class csspec(ctools.cscript):
394 psfcube["proj"].string(self.m_proj) 399 psfcube["proj"].string(self.m_proj)
395 psfcube.run() 400 psfcube.run()
396 401
  402 + # Header
397 if self.logTerse(): 403 if self.logTerse():
398 self.log.header3("Creating background cube") 404 self.log.header3("Creating background cube")
399 405
400 - # create background cube 406 + # Create background cube
401 bkgcube = ctools.ctbkgcube(select.obs()) 407 bkgcube = ctools.ctbkgcube(select.obs())
402 bkgcube["incube"].filename("NONE") 408 bkgcube["incube"].filename("NONE")
403 bkgcube["usepnt"].boolean(False) 409 bkgcube["usepnt"].boolean(False)
@@ -433,10 +439,11 @@ class csspec(ctools.cscript): @@ -433,10 +439,11 @@ class csspec(ctools.cscript):
433 # Set new models to binned observation 439 # Set new models to binned observation
434 obs.models(models) 440 obs.models(models)
435 441
  442 + # Header
436 if self.logTerse(): 443 if self.logTerse():
437 self.log.header3("Performing fit") 444 self.log.header3("Performing fit")
438 445
439 - # likelihood 446 + # Likelihood
440 like = ctools.ctlike(obs) 447 like = ctools.ctlike(obs)
441 like.run() 448 like.run()
442 449
@@ -446,10 +453,11 @@ class csspec(ctools.cscript): @@ -446,10 +453,11 @@ class csspec(ctools.cscript):
446 # Log information 453 # Log information
447 if self.logTerse(): 454 if self.logTerse():
448 self.log("No event in this bin. Bin is skipped\n") 455 self.log("No event in this bin. Bin is skipped\n")
449 -  
450 - flux[i] = 0.0  
451 - flux_err[i] = 0.0  
452 - TSvalues[i] = 0.0 456 +
  457 + # Set all values to 0
  458 + flux[i] = 0.0
  459 + flux_err[i] = 0.0
  460 + TSvalues[i] = 0.0
453 ulim_values[i] = 0.0 461 ulim_values[i] = 0.0
454 continue 462 continue
455 463
@@ -475,7 +483,6 @@ class csspec(ctools.cscript): @@ -475,7 +483,6 @@ class csspec(ctools.cscript):
475 ulimit.run() 483 ulimit.run()
476 ulimit_value = ulimit.diff_ulimit() 484 ulimit_value = ulimit.diff_ulimit()
477 except: 485 except:
478 -  
479 if self.logTerse(): 486 if self.logTerse():
480 self.log("Upper limit calculation failed\n") 487 self.log("Upper limit calculation failed\n")
481 ulimit_value = -1.0 488 ulimit_value = -1.0
test/test_cscripts.sh.in
@@ -191,7 +191,7 @@ $ECHO " ok" @@ -191,7 +191,7 @@ $ECHO " ok"
191 191
192 # 192 #
193 # Test csspec 193 # Test csspec
194 -# ============= 194 +# ===========
195 $ECHO -n "Test csspec: " 195 $ECHO -n "Test csspec: "
196 chmod +x @abs_top_srcdir@/scripts/csspec.py 196 chmod +x @abs_top_srcdir@/scripts/csspec.py
197 @abs_top_srcdir@/scripts/csspec.py \ 197 @abs_top_srcdir@/scripts/csspec.py \
@@ -217,7 +217,7 @@ $ECHO " ok" @@ -217,7 +217,7 @@ $ECHO " ok"
217 217
218 # 218 #
219 # Test csspec Run2 219 # Test csspec Run2
220 -# ============= 220 +# ================
221 $ECHO -n "Test csspec: " 221 $ECHO -n "Test csspec: "
222 chmod +x @abs_top_srcdir@/scripts/csspec.py 222 chmod +x @abs_top_srcdir@/scripts/csspec.py
223 @abs_top_srcdir@/scripts/csspec.py \ 223 @abs_top_srcdir@/scripts/csspec.py \