Pierrick Martin / ctools

Commit d0742fe74ee7a903675937093c10a0b23cba4e9c

Authored by Jürgen Knödlseder
1 parent 6011bf73

Change implementation of stacked cube handling in csspec.

Inline Side-by-side
Showing 1 changed file with 114 additions and 65 deletions
cscripts/csspec.py
  View file @d0742fe
@@ -110,10 +110,11 @@ class csspec(ctools.cscript): @@ -110,10 +110,11 @@ class csspec(ctools.cscript):
110 pars.append(gammalib.GApplicationPar("caldb","s","a","prod2","","","Calibration database")) 110 pars.append(gammalib.GApplicationPar("caldb","s","a","prod2","","","Calibration database"))
111 pars.append(gammalib.GApplicationPar("irf","s","a","South_50h","","","Instrument response function")) 111 pars.append(gammalib.GApplicationPar("irf","s","a","South_50h","","","Instrument response function"))
112 pars.append(gammalib.GApplicationPar("edisp","b","h","no","","","Apply energy dispersion?")) 112 pars.append(gammalib.GApplicationPar("edisp","b","h","no","","","Apply energy dispersion?"))
113 - pars.append(gammalib.GApplicationPar("emin","r","h","0.1","","","Lower energy limit for spectral points(TeV)"))  
114 - pars.append(gammalib.GApplicationPar("emax","r","h","100.0","","","Upper energy limit for spectral points(TeV)"))  
115 - pars.append(gammalib.GApplicationPar("enumbins","i","a","20","","","Number of spectral points"))  
116 - pars.append(gammalib.GApplicationPar("ebinalg","s", "h","LOG","FILE|LIN|LOG","", "Algorithm for defining energy bins")) 113 + pars.append(gammalib.GApplicationPar("emin","r","a","0.1","","","Lower energy limit for spectral points (TeV)"))
  114 + pars.append(gammalib.GApplicationPar("emax","r","a","100.0","","","Upper energy limit for spectral points (TeV)"))
  115 + pars.append(gammalib.GApplicationPar("enumbins","i","a","20","1","10000","Number of spectral points"))
  116 + pars.append(gammalib.GApplicationPar("ebinalg","s","a","LOG","FILE|LIN|LOG","", "Algorithm for defining energy bins"))
  117 + pars.append(gammalib.GApplicationPar("ebinfile","f","a","NONE","","", "Name of the file containing the energy bin definition"))
117 pars.append(gammalib.GApplicationPar("calc_ts","b","h","yes","yes|no","","Compute TS value in each bin")) 118 pars.append(gammalib.GApplicationPar("calc_ts","b","h","yes","yes|no","","Compute TS value in each bin"))
118 pars.append(gammalib.GApplicationPar("calc_ulim","b","h","yes","yes|no","","Compute upper limit in each bin")) 119 pars.append(gammalib.GApplicationPar("calc_ulim","b","h","yes","yes|no","","Compute upper limit in each bin"))
119 pars.append(gammalib.GApplicationPar("fix_srcs","b","h","yes","yes|no","","Fix other skymodel parameters")) 120 pars.append(gammalib.GApplicationPar("fix_srcs","b","h","yes","yes|no","","Fix other skymodel parameters"))
@@ -129,7 +130,6 @@ class csspec(ctools.cscript): @@ -129,7 +130,6 @@ class csspec(ctools.cscript):
129 """ 130 """
130 Get parameters from parfile and setup the observation. 131 Get parameters from parfile and setup the observation.
131 """ 132 """
132 -  
133 # Set observation if not done before 133 # Set observation if not done before
134 if self.obs == None or self.obs.size() == 0: 134 if self.obs == None or self.obs.size() == 0:
135 self.require_inobs("csspec::get_parameters()") 135 self.require_inobs("csspec::get_parameters()")
@@ -149,41 +149,8 @@ class csspec(ctools.cscript): @@ -149,41 +149,8 @@ class csspec(ctools.cscript):
149 # Get source name 149 # Get source name
150 self.m_srcname = self["srcname"].string() 150 self.m_srcname = self["srcname"].string()
151 151
152 - # Get ebounds  
153 - if self.m_binned_mode:  
154 -  
155 - # Extract cube ebounds  
156 - cube_ebounds = self.obs[0].events().ebounds()  
157 -  
158 - # Read user parameters  
159 - self.m_enumbins = self["enumbins"].integer()  
160 - self.m_emin = self["emin"].real() - 1e-6 #rounding tolerance  
161 - self.m_emax = self["emax"].real() + 1e-6 #rounding tolerance  
162 -  
163 - # Compute number of energy layers used for one spectral bin  
164 - n_layers = cube_ebounds.size() / self.m_enumbins  
165 -  
166 - # Don't allow zero layers  
167 - if n_layers == 0:  
168 - n_layers = 1  
169 -  
170 - # Create new ebounds object and fill it  
171 - self.m_ebounds = gammalib.GEbounds()  
172 - for i in range(cube_ebounds.size())[::n_layers]:  
173 - if i + n_layers > cube_ebounds.size():  
174 - break  
175 -  
176 - # Get new emin and emax of bin  
177 - emin = cube_ebounds.emin(i)  
178 - emax = cube_ebounds.emax(i + n_layers - 1)  
179 -  
180 - # Append energy bin if within user range  
181 - if self.m_emin <= emin.TeV() and self.m_emax >= emax.TeV():  
182 - self.m_ebounds.append(emin, emax)  
183 -  
184 - # Unbinned mode  
185 - else:  
186 - self.m_ebounds = self.create_ebounds() 152 + # Set ebounds
  153 + self.set_ebounds()
187 154
188 # Get edisp flag 155 # Get edisp flag
189 self.m_edisp = self["edisp"].boolean() 156 self.m_edisp = self["edisp"].boolean()
@@ -206,6 +173,66 @@ class csspec(ctools.cscript): @@ -206,6 +173,66 @@ class csspec(ctools.cscript):
206 # Return 173 # Return
207 return 174 return
208 175
  176 + def set_ebounds(self):
  177 + """
  178 + Set energy boundaries.
  179 + """
  180 + # Get ebounds
  181 + if self.m_binned_mode:
  182 +
  183 + # Extract cube ebounds
  184 + cube_ebounds = self.obs[0].events().ebounds()
  185 +
  186 + # Read user parameters
  187 + self.m_enumbins = self["enumbins"].integer()
  188 + self.m_emin = self["emin"].real() - 1e-6 # Rounding tolerance
  189 + self.m_emax = self["emax"].real() + 1e-6 # Rounding tolerance
  190 +
  191 + # Set energy range
  192 + emin = gammalib.GEnergy(self.m_emin, "TeV")
  193 + emax = gammalib.GEnergy(self.m_emax, "TeV")
  194 +
  195 + # Determine all cube layers that overlap with the [emin,emax]
  196 + # energy range.
  197 + use_layers = []
  198 + for i in range(cube_ebounds.size()):
  199 + if cube_ebounds.emin(i) < emax and cube_ebounds.emax(i) > emin:
  200 + use_layers.append(i)
  201 +
  202 + # If no layers overlap then throw an exception
  203 + if (len(use_layers) == 0):
  204 + msg = "Energy range ["+str(cube_ebounds.emin(i))+ \
  205 + ", "+str(cube_ebounds.emax(i))+"] of counts "+ \
  206 + "cube does not overlap with specified energy range ["+ \
  207 + str(emin)+", "+str(emax)+"]. Specify "+ \
  208 + "overlapping energy range."
  209 + raise gammalib.GException.invalid_argument("csspec", msg)
  210 +
  211 + # Determine number of layers to use for each spectral bin
  212 + n_layers = int(len(use_layers) / self.m_enumbins)
  213 +
  214 + # Don't allow zero layers
  215 + if n_layers == 0:
  216 + n_layers = 1
  217 +
  218 + # Create energy boundaries
  219 + self.m_ebounds = gammalib.GEbounds()
  220 + for i in range(len(use_layers))[::n_layers]:
  221 + i_start = i
  222 + i_stop = i + n_layers - 1
  223 + if i_stop >= len(use_layers):
  224 + i_stop = len(use_layers) - 1
  225 + ebinmin = cube_ebounds.emin(use_layers[i_start])
  226 + ebinmax = cube_ebounds.emax(use_layers[i_stop])
  227 + self.m_ebounds.append(ebinmin, ebinmax)
  228 +
  229 + # Unbinned mode
  230 + else:
  231 + self.m_ebounds = self.create_ebounds()
  232 +
  233 + # Return
  234 + return
  235 +
209 def models(self, models): 236 def models(self, models):
210 """ 237 """
211 Set model. 238 Set model.
@@ -252,6 +279,24 @@ class csspec(ctools.cscript): @@ -252,6 +279,24 @@ class csspec(ctools.cscript):
252 self.log_parameters() 279 self.log_parameters()
253 self.log("\n") 280 self.log("\n")
254 281
  282 + # Write spectral binning into header
  283 + if self.logTerse():
  284 + self.log("\n")
  285 + self.log.header1("Spectral binning")
  286 + if self.m_binned_mode:
  287 + cube_ebounds = self.obs[0].events().ebounds()
  288 + self.log.parformat("Counts cube energy range")
  289 + self.log(str(cube_ebounds.emin()))
  290 + self.log(" - ")
  291 + self.log(str(cube_ebounds.emax()))
  292 + self.log("\n")
  293 + for i in range(self.m_ebounds.size()):
  294 + self.log.parformat("Bin "+str(i+1))
  295 + self.log(str(self.m_ebounds.emin(i)))
  296 + self.log(" - ")
  297 + self.log(str(self.m_ebounds.emax(i)))
  298 + self.log("\n")
  299 +
255 # Write observation into logger 300 # Write observation into logger
256 if self.logTerse(): 301 if self.logTerse():
257 self.log("\n") 302 self.log("\n")
@@ -316,26 +361,25 @@ class csspec(ctools.cscript): @@ -316,26 +361,25 @@ class csspec(ctools.cscript):
316 361
317 # Create FITS table columns 362 # Create FITS table columns
318 energy = gammalib.GFitsTableDoubleCol("Energy", self.m_ebounds.size()) 363 energy = gammalib.GFitsTableDoubleCol("Energy", self.m_ebounds.size())
319 - energy.unit("TeV")  
320 energy_low = gammalib.GFitsTableDoubleCol("ed_Energy", self.m_ebounds.size()) 364 energy_low = gammalib.GFitsTableDoubleCol("ed_Energy", self.m_ebounds.size())
321 - energy_low.unit("TeV")  
322 energy_high = gammalib.GFitsTableDoubleCol("eu_Energy", self.m_ebounds.size()) 365 energy_high = gammalib.GFitsTableDoubleCol("eu_Energy", self.m_ebounds.size())
323 - energy_high.unit("TeV")  
324 flux = gammalib.GFitsTableDoubleCol("Flux", self.m_ebounds.size()) 366 flux = gammalib.GFitsTableDoubleCol("Flux", self.m_ebounds.size())
325 - flux.unit("erg/cm2/s")  
326 flux_err = gammalib.GFitsTableDoubleCol("e_Flux", self.m_ebounds.size()) 367 flux_err = gammalib.GFitsTableDoubleCol("e_Flux", self.m_ebounds.size())
327 - flux_err.unit("erg/cm2/s")  
328 TSvalues = gammalib.GFitsTableDoubleCol("TS", self.m_ebounds.size()) 368 TSvalues = gammalib.GFitsTableDoubleCol("TS", self.m_ebounds.size())
329 ulim_values = gammalib.GFitsTableDoubleCol("UpperLimit", self.m_ebounds.size()) 369 ulim_values = gammalib.GFitsTableDoubleCol("UpperLimit", self.m_ebounds.size())
330 - ulim_values.unit("erg/cm2/s")  
331 Npred_values = gammalib.GFitsTableDoubleCol("Npred", self.m_ebounds.size()) 370 Npred_values = gammalib.GFitsTableDoubleCol("Npred", self.m_ebounds.size())
332 - 371 + energy.unit("TeV")
  372 + energy_low.unit("TeV")
  373 + energy_high.unit("TeV")
  374 + flux.unit("erg/cm2/s")
  375 + flux_err.unit("erg/cm2/s")
  376 + ulim_values.unit("erg/cm2/s")
333 377
334 # Loop over energy bins 378 # Loop over energy bins
335 for i in range(self.m_ebounds.size()): 379 for i in range(self.m_ebounds.size()):
336 380
337 # Log information 381 # Log information
338 - if self.logTerse(): 382 + if self.logExplicit():
339 self.log("\n") 383 self.log("\n")
340 self.log.header2("Energy bin "+str(i+1)) 384 self.log.header2("Energy bin "+str(i+1))
341 385
@@ -377,9 +421,10 @@ class csspec(ctools.cscript): @@ -377,9 +421,10 @@ class csspec(ctools.cscript):
377 else: 421 else:
378 422
379 # Header 423 # Header
380 - if self.logTerse():  
381 - self.log.header3("Selecting events") 424 + if self.logExplicit():
  425 + self.log.header3("Filtering cube")
382 426
  427 + # Select layers
383 cubemask = ctools.ctcubemask(self.obs) 428 cubemask = ctools.ctcubemask(self.obs)
384 cubemask["regfile"] = "NONE" 429 cubemask["regfile"] = "NONE"
385 cubemask["ra"] = "UNDEFINED" 430 cubemask["ra"] = "UNDEFINED"
@@ -393,7 +438,7 @@ class csspec(ctools.cscript): @@ -393,7 +438,7 @@ class csspec(ctools.cscript):
393 obs = cubemask.obs() 438 obs = cubemask.obs()
394 439
395 # Header 440 # Header
396 - if self.logTerse(): 441 + if self.logExplicit():
397 self.log.header3("Performing fit") 442 self.log.header3("Performing fit")
398 443
399 # Likelihood 444 # Likelihood
@@ -405,8 +450,10 @@ class csspec(ctools.cscript): @@ -405,8 +450,10 @@ class csspec(ctools.cscript):
405 if like.obs().logL() == 0.0: 450 if like.obs().logL() == 0.0:
406 451
407 # Log information 452 # Log information
408 - if self.logTerse():  
409 - self.log("No event in this bin. Likelihood is zero. Bin is skipped\n") 453 + if self.logExplicit():
  454 + self.log("No event in this bin. ")
  455 + self.log("Likelihood is zero. ")
  456 + self.log("Bin is skipped.")
410 457
411 # Set all values to 0 458 # Set all values to 0
412 flux[i] = 0.0 459 flux[i] = 0.0
@@ -425,21 +472,21 @@ class csspec(ctools.cscript): @@ -425,21 +472,21 @@ class csspec(ctools.cscript):
425 if self.m_calc_ulimit: 472 if self.m_calc_ulimit:
426 473
427 # Logging information 474 # Logging information
428 - if self.logTerse(): 475 + if self.logExplicit():
429 self.log.header3("Computing upper limit") 476 self.log.header3("Computing upper limit")
430 477
431 # Create upper limit object 478 # Create upper limit object
432 ulimit = ctools.ctulimit(like.obs()) 479 ulimit = ctools.ctulimit(like.obs())
433 ulimit["srcname"] = self.m_srcname 480 ulimit["srcname"] = self.m_srcname
434 - ulimit["eref"] = elogmean.TeV() 481 + ulimit["eref"] = elogmean.TeV()
435 482
436 # Try to run upper limit and catch exceptions 483 # Try to run upper limit and catch exceptions
437 try: 484 try:
438 ulimit.run() 485 ulimit.run()
439 ulimit_value = ulimit.diff_ulimit() 486 ulimit_value = ulimit.diff_ulimit()
440 except: 487 except:
441 - if self.logTerse():  
442 - self.log("Upper limit calculation failed\n") 488 + if self.logExplicit():
  489 + self.log("Upper limit calculation failed.")
443 ulimit_value = -1.0 490 ulimit_value = -1.0
444 491
445 # Get TS value 492 # Get TS value
@@ -474,15 +521,17 @@ class csspec(ctools.cscript): @@ -474,15 +521,17 @@ class csspec(ctools.cscript):
474 ulim_values[i] = ulimit_value * elogmean2 * gammalib.MeV2erg 521 ulim_values[i] = ulimit_value * elogmean2 * gammalib.MeV2erg
475 522
476 # Log information 523 # Log information
477 - if self.logExplicit():  
478 - self.log("Bin "+str(i)+" ["+str(emin.TeV())+"-"+str(emax.TeV())+"] TeV: ")  
479 - self.log("Flux = "+str(flux[i]))  
480 - self.log(" +- "+str(flux_err[i])+" [erg/cm2/s]")  
481 - if self.m_calc_ts and TSvalues[i] > 0.0:  
482 - self.log(", TS = "+str(TS))  
483 - if self.m_calc_ulimit and ulim_values[i] > 0.0:  
484 - self.log(", UL = "+str(ulim_values[i])+" [erg/cm2/s]") 524 + if self.logTerse():
485 self.log("\n") 525 self.log("\n")
  526 + self.log.parformat("Bin "+str(i+1))
  527 + self.log(str(flux[i]))
  528 + self.log(" +/- ")
  529 + self.log(str(flux_err[i]))
  530 + if self.m_calc_ulimit and ulim_values[i] > 0.0:
  531 + self.log(" [< "+str(ulim_values[i])+"]")
  532 + self.log(" erg/cm2/s")
  533 + if self.m_calc_ts and TSvalues[i] > 0.0:
  534 + self.log(" (TS = "+str(TS)+")")
486 535
487 # Append filled columns to fits table 536 # Append filled columns to fits table
488 table.append(energy) 537 table.append(energy)