Johan Bregeon / ctools

Commit c2b4bfbe98d4798c04ccfe72804b0de60057790d

Authored by Michael Mayer
1 parent 4d7dc233

added bgtype parameter to cshessobs

Inline Side-by-side
Showing 1 changed file with 128 additions and 75 deletions
cscripts/cshessobs.py
  View file @c2b4bfb
... ... @@ -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,130 @@ 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()
273   -
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   -
  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 + # Use several energy-dependent scaling factors
  294 + fits = gammalib.GFits(aefffile)
  295 + emin = fits["EFFECTIVE AREA"].real("LO_THRES")
  296 + emax = fits["EFFECTIVE AREA"].real("HI_THRES")
  297 + fits.close()
  298 +
  299 + # Create reference powerlaw
  300 + plaw = gammalib.GModelSpectralPlaw(prefactor, index, gammalib.GEnergy(1.0,"TeV"))
  301 +
  302 + # Create spectral model and energy values
  303 + spec = gammalib.GModelSpectralNodes()
  304 + bounds = gammalib.GEbounds(self.bkgpars,gammalib.GEnergy(emin,"TeV"),gammalib.GEnergy(emax,"TeV"),True)
  305 + for i in range(bounds.size()):
  306 + energy = bounds.elogmean(i)
  307 + spec.append(energy, plaw.eval(energy, gammalib.GTime()))
  308 + for par in spec:
  309 + if "Energy" in par.name():
  310 + par.fix()
  311 + elif "Intensity" in par.name():
  312 + par.scale(scale)
  313 + par.min(0.01*par.scale())
  314 + par.max(100.0*par.scale())
  315 +# print prefactor, index
  316 +# print spec
  317 + return spec
  318 +
  319 + def hess_background(self, run, aefffile, trgrate, zenith):
  320 +
  321 + # Set IrfBackground
  322 + if self.bkgtype == "irf":
  323 + prefactor = 1.0
  324 + index = 0.0
  325 + if self.usetrig and trgrate > 0.0:
  326 + # Use value from lookup
  327 + prefactor = 0.00427233352611 * trgrate + 0.181244590076
  328 + spec = self.background_spectrum(run, prefactor, index, aefffile)
  329 + # Create background model instance
  330 + bck = gammalib.GCTAModelIrfBackground(spec)
  331 +
  332 + # Set AeffBackground
  333 + elif self.bkgtype == "aeff":
  334 + prefactor = 5e-14
  335 + index = -2.5
  336 + # Use values from lookup if possible
  337 + if self.usetrig and trgrate > 0.0:
  338 + prefactor = 1.52976410372e-16 * trgrate + 2.71805547108e-14
  339 + if zenith > 0.0:
  340 + index = 0.000271928721732 * zenith * zenith -0.0029602286453 * zenith -2.76404600416
  341 +
  342 + spec = self.background_spectrum(run, prefactor, index, aefffile)
  343 +
  344 + # Create background model instance
  345 + bck = gammalib.GCTAModelAeffBackground(spec)
  346 +
  347 + # Set Gaussian Background
  348 + elif self.bkgtype == "gauss":
  349 + spec = self.background_spectrum(run, 1e-4, -1.8, aefffile)
  350 + radial = gammalib.GCTAModelRadialGauss(2.5)
  351 + bck = gammalib.GCTAModelRadialAcceptance(radial, spec)
  352 +
  353 + else:
  354 + sys.exit("Background type \""+self.bkgtype+"\" unsupported")
  355 +
  356 + # Copy model
  357 + model = bck.clone()
  358 +
309 359 # Assign specific run id
310   - bck.ids(str(int(run)))
311   -
  360 + model.ids(str(int(run)))
  361 +
312 362 # Assign instrument
313   - bck.instruments("HESS")
314   -
  363 + model.instruments("HESS")
  364 +
315 365 # Set name (arbitrary)
316   - bck.name("bkg_"+run)
317   -
  366 + model.name("bkg_"+run)
  367 +
318 368 # Turn off TS calculation for background model
319   - bck.tscalc(False)
320   -
  369 + model.tscalc(False)
  370 +
321 371 # Return model
322   - return bck
323   -
  372 + return model
  373 +
324 374  
325 375 def erange(self):
326 376 """
... ... @@ -465,10 +515,13 @@ class cshessobs(ctools.cscript):
465 515 fits = gammalib.GFits(eventfile)
466 516 object_name = fits["EVENTS"].string("OBJECT")
467 517  
468   - # Retrieve trigger rate if available
  518 + # Retrieve trigger rate and zenith if available
469 519 trgrate = 0.0
  520 + zenith = 0.0
470 521 if fits["EVENTS"].has_card("ZTRGRATE"):
471 522 trgrate = fits["EVENTS"].real("ZTRGRATE")
  523 + if fits["EVENTS"].has_card("ALT_PNT"):
  524 + zenith = 90.-fits["EVENTS"].real("ALT_PNT")
472 525  
473 526 # Close FITS file
474 527 fits.close()
... ... @@ -515,10 +568,10 @@ class cshessobs(ctools.cscript):
515 568 obs.append(psf)
516 569 obs.append(edisp)
517 570 obs.append(bck)
518   -
  571 +
519 572 # Append instrumental background model
520   - self.models.append(self.hess_inst_background(run,aefffile,trgrate))
521   -
  573 + self.models.append(self.hess_background(run, aefffile, trgrate, zenith))
  574 +
522 575 # Continue only if there are observations available
523 576 if lib.size():
524 577  
... ...