Mathieu Servillat / ctools

Commit b31d3a498c18edab86e03026f22acf16ce85a105

Authored by Jürgen Knödlseder
2 parents eedbbc62 6d012d33

Merge branch '1527-more-bgmodels-in-cshessobs' into devel

Inline Side-by-side
Showing 1 changed file with 136 additions and 78 deletions
cscripts/cshessobs.py
  View file @b31d3a4
... ... @@ -23,6 +23,7 @@ import ctools
23 23 import sys
24 24 import os
25 25 import glob
  26 +import math
26 27  
27 28  
28 29 # =============== #
... ... @@ -43,7 +44,7 @@ class cshessobs(ctools.cscript):
43 44 """
44 45 # Set name
45 46 self.name = "cshessobs"
46   - self.version = "0.1.0"
  47 + self.version = "0.2.0"
47 48 self.datapath = ""
48 49 self.m_erange = gammalib.GEbounds()
49 50  
... ... @@ -57,7 +58,6 @@ class cshessobs(ctools.cscript):
57 58 # Initialise some members
58 59 self.inmodels = None
59 60 self.outobs = "obs.xml"
60   - self.nodes = 1
61 61  
62 62 # Make sure that parfile exists
63 63 file = self.parfile()
... ... @@ -165,8 +165,9 @@ class cshessobs(ctools.cscript):
165 165 pars.append(gammalib.GApplicationPar("outobs","f","a","obs.xml","","","Observation XML outfile"))
166 166 pars.append(gammalib.GApplicationPar("outmodel","f","a","bgmodels.xml","","","Output model XML file"))
167 167 pars.append(gammalib.GApplicationPar("datapath","s","h",self.datapath,"","","Data directoy"))
168   - pars.append(gammalib.GApplicationPar("nodes","i","h","1","","","Number of free parameters per background model"))
  168 + pars.append(gammalib.GApplicationPar("bkgpars","i","a","1","","","Number of free parameters per background model"))
169 169 pars.append(gammalib.GApplicationPar("usetrig","b","h","yes","yes|no","","Pre-estimate background normalisation? Recommended for binned analysis"))
  170 + pars.append(gammalib.GApplicationPar("bkgtype","s","a","irf","gauss|irf|aeff","","What background model should be used"))
170 171 pars.append_standard()
171 172 pars.save(parfile)
172 173  
... ... @@ -192,7 +193,8 @@ class cshessobs(ctools.cscript):
192 193 self.outmodel = self["outmodel"].filename()
193 194  
194 195 # Read hidden parameters
195   - self.nodes = self["nodes"].integer()
  196 + self.bkgpars = self["bkgpars"].integer()
  197 + self.bkgtype = self["bkgtype"].string()
196 198 self.usetrig = self["usetrig"].boolean()
197 199 self.datapath = self["datapath"].string()
198 200  
... ... @@ -245,82 +247,133 @@ class cshessobs(ctools.cscript):
245 247  
246 248 # Return valid configs
247 249 return newfolders
248   -
249   - # Create instrument background model
250   - def hess_inst_background(self,run,aefffile,trgrate=0.0):
251   -
252   - # Initialise return spectrum
253   - spec = 0
254   -
255   - # Initialise default value
256   - value = 1.0
257   -
258   - if self.usetrig and trgrate > 0.0:
259   - # Use value from lookup
260   - value = 0.00427233352611 * trgrate + 0.181244590076
261   -
262   - if self.nodes <= 1:
263   -
264   - # Use a constant scaling function
  250 +
  251 + def background_spectrum(self, run, prefactor, index, aefffile = ""):
  252 +
  253 + scale = pow(10,math.floor(math.log10(abs(prefactor))))
  254 +
  255 + if index == 0.0 and self.bkgpars <= 1:
265 256 spec = gammalib.GModelSpectralConst()
266   - spec["Value"].min(0.1)
267   - spec["Value"].max(10)
268   - spec["Value"].value(value)
269   - if self.nodes < 1:
  257 + spec["Value"].min(0.1*scale)
  258 + spec["Value"].max(10*scale)
  259 + spec["Value"].value(prefactor)
  260 + spec["Value"].scale(scale)
  261 + if self.bkgpars == 0:
270 262 spec["Value"].fix()
271 263 else:
272 264 spec["Value"].free()
  265 +
  266 + else:
  267 +
  268 + # Create power law if number of free parameters are less 2 or less
  269 + if self.bkgpars <= 2:
  270 + e = gammalib.GEnergy(1.0,"TeV")
  271 + spec = gammalib.GModelSpectralPlaw(prefactor, index, e)
  272 +
  273 + # Set parameter ranges
  274 + spec[0].scale(scale)
  275 + spec[0].min(0.01 * scale)
  276 + spec[0].max(100.0 * scale)
  277 + spec[1].scale(1)
  278 + spec[1].min(-5.0)
  279 + spec[1].max(5.0)
  280 +
  281 + # Set number of free parameters
  282 + if self.bkgpars == 0:
  283 + spec[0].fix()
  284 + spec[1].fix()
  285 + elif self.bkgpars == 1:
  286 + spec[0].free()
  287 + spec[1].fix()
  288 + else:
  289 + spec[0].free()
  290 + spec[1].free()
  291 +
  292 + else:
  293 +
  294 + # Use several energy-dependent scaling factors
  295 + fits = gammalib.GFits(aefffile)
  296 + emin = fits["EFFECTIVE AREA"].real("LO_THRES")
  297 + emax = fits["EFFECTIVE AREA"].real("HI_THRES")
  298 + fits.close()
  299 +
  300 + # Create reference powerlaw
  301 + plaw = gammalib.GModelSpectralPlaw(prefactor, index, gammalib.GEnergy(1.0,"TeV"))
  302 +
  303 + # Create spectral model and energy values
  304 + spec = gammalib.GModelSpectralNodes()
  305 + bounds = gammalib.GEbounds(self.bkgpars,gammalib.GEnergy(emin,"TeV"),gammalib.GEnergy(emax,"TeV"),True)
  306 + for i in range(bounds.size()):
  307 + energy = bounds.elogmean(i)
  308 + spec.append(energy, plaw.eval(energy, gammalib.GTime()))
  309 + for par in spec:
  310 +
  311 + if "Energy" in par.name():
  312 + par.fix()
  313 + elif "Intensity" in par.name():
  314 + parscale = pow(10,math.floor(math.log10(abs(par.value()))))
  315 + par.scale(parscale)
  316 + par.min(0.01*parscale)
  317 + par.max(100.0*parscale)
  318 +
  319 +
  320 + return spec
273 321  
274   - elif self.nodes == 2:
275   -
276   - # User a power law scaling function
277   - e = gammalib.GEnergy(1.0,"TeV")
278   - spec = gammalib.GModelSpectralPlaw(1.0,0.0,e)
279   - spec[0].free()
280   - spec[0].min(0.01)
281   - spec[0].max(10.0)
282   - spec[0].scale(1.0)
283   - spec[0].value(value)
284   - spec[1].free()
285   - spec[1].min(-2.0)
286   - spec[1].max(2.0)
287   - spec[1].scale(1.0)
288   - else:
289   -
290   - # Use several energy-dependent scaling factors
291   - fits = gammalib.GFits(aefffile)
292   - emin = fits["EFFECTIVE AREA"].real("LO_THRES")
293   - emax = fits["EFFECTIVE AREA"].real("HI_THRES")
294   - fits.close()
295   - spec = gammalib.GModelSpectralNodes()
296   - bounds = gammalib.GEbounds(self.nodes,gammalib.GEnergy(emin,"TeV"),gammalib.GEnergy(emax,"TeV"),True)
297   - for i in range(bounds.size()):
298   - spec.append(bounds.elogmean(i),value)
299   - for par in spec:
300   - if "Energy" in par.name():
301   - par.fix()
302   - elif "Intensity" in par.name():
303   - par.min(0.01)
304   - par.max(10.0)
305   -
306   - # Create background model instance
307   - bck = gammalib.GCTAModelIrfBackground(spec)
308   -
  322 + def hess_background(self, run, aefffile, trgrate, zenith):
  323 +
  324 + # Set IrfBackground
  325 + if self.bkgtype == "irf":
  326 + prefactor = 1.0
  327 + index = 0.0
  328 + if self.usetrig and trgrate > 0.0:
  329 + # Use value from lookup
  330 + prefactor = 0.00427233352611 * trgrate + 0.181244590076
  331 + spec = self.background_spectrum(run, prefactor, index, aefffile)
  332 + # Create background model instance
  333 + bck = gammalib.GCTAModelIrfBackground(spec)
  334 +
  335 + # Set AeffBackground
  336 + elif self.bkgtype == "aeff":
  337 + prefactor = 5e-14
  338 + index = -2.5
  339 + # Use values from lookup if possible
  340 + if self.usetrig and trgrate > 0.0:
  341 + prefactor = 1.52976410372e-16 * trgrate + 2.71805547108e-14
  342 + if zenith > 0.0:
  343 + index = 0.000271928721732 * zenith * zenith -0.0029602286453 * zenith -2.76404600416
  344 +
  345 + spec = self.background_spectrum(run, prefactor, index, aefffile)
  346 +
  347 + # Create background model instance
  348 + bck = gammalib.GCTAModelAeffBackground(spec)
  349 +
  350 + # Set Gaussian Background
  351 + elif self.bkgtype == "gauss":
  352 + spec = self.background_spectrum(run, 1e-4, -1.8, aefffile)
  353 + radial = gammalib.GCTAModelRadialGauss(2.5)
  354 + bck = gammalib.GCTAModelRadialAcceptance(radial, spec)
  355 +
  356 + else:
  357 + sys.exit("Background type \""+self.bkgtype+"\" unsupported")
  358 +
  359 + # Copy model
  360 + model = bck.clone()
  361 +
309 362 # Assign specific run id
310   - bck.ids(str(int(run)))
311   -
  363 + model.ids(str(int(run)))
  364 +
312 365 # Assign instrument
313   - bck.instruments("HESS")
314   -
  366 + model.instruments("HESS")
  367 +
315 368 # Set name (arbitrary)
316   - bck.name("bkg_"+run)
317   -
  369 + model.name("bkg_"+run)
  370 +
318 371 # Turn off TS calculation for background model
319   - bck.tscalc(False)
320   -
  372 + model.tscalc(False)
  373 +
321 374 # Return model
322   - return bck
323   -
  375 + return model
  376 +
324 377  
325 378 def erange(self):
326 379 """
... ... @@ -441,11 +494,13 @@ class cshessobs(ctools.cscript):
441 494 msg = "Run "+str(int(run))+" has no PSF - Run is skipped"
442 495 skip = True
443 496 elif not os.path.isfile(edispfile):
444   - msg = "Run "+str(int(run))+" has no energy dispersion - Run is skipped"
445   - skip = True
  497 + msg = "Warning: Run "+str(int(run))+" has no energy dispersion - usage of energy dispersion not possible"
  498 + skip = False
  499 + edispfile = ""
446 500 elif not os.path.isfile(bgfile):
447   - msg = "Run "+str(int(run))+" has no background - Run is skipped"
448   - skip = True
  501 + msg = "Warning: Run "+str(int(run))+" has no background file - use of IRF background not possible"
  502 + skip = False
  503 + bgfile = ""
449 504 else:
450 505 msg = "Adding run "+str(int(run))
451 506  
... ... @@ -465,10 +520,13 @@ class cshessobs(ctools.cscript):
465 520 fits = gammalib.GFits(eventfile)
466 521 object_name = fits["EVENTS"].string("OBJECT")
467 522  
468   - # Retrieve trigger rate if available
  523 + # Retrieve trigger rate and zenith if available
469 524 trgrate = 0.0
  525 + zenith = 0.0
470 526 if fits["EVENTS"].has_card("ZTRGRATE"):
471 527 trgrate = fits["EVENTS"].real("ZTRGRATE")
  528 + if fits["EVENTS"].has_card("ALT_PNT"):
  529 + zenith = 90.-fits["EVENTS"].real("ALT_PNT")
472 530  
473 531 # Close FITS file
474 532 fits.close()
... ... @@ -515,10 +573,10 @@ class cshessobs(ctools.cscript):
515 573 obs.append(psf)
516 574 obs.append(edisp)
517 575 obs.append(bck)
518   -
  576 +
519 577 # Append instrumental background model
520   - self.models.append(self.hess_inst_background(run,aefffile,trgrate))
521   -
  578 + self.models.append(self.hess_background(run, aefffile, trgrate, zenith))
  579 +
522 580 # Continue only if there are observations available
523 581 if lib.size():
524 582  
... ...