Alexander Ziegler / ctools

Commit 64e71e40f4c82d94c3671198249f47fb204e4d40

Authored by Jürgen Knödlseder
2 parents 45396196 7c914de9

Merge branch '00-07-02' into integration 

Conflicts:
	ChangeLog
	NEWS
	scripts/cspull.py
Inline Side-by-side
Showing 8 changed files with 1596 additions and 1517 deletions
ChangeLog
  View file @64e71e4
... ... @@ -6,6 +6,14 @@
6 6 Support energy dispersion in ctobssim, ctlike and ctmodel
7 7  
8 8  
  9 +2014-02-10 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
  10 +
  11 + * Bug fix Version 0.7.2 released
  12 + ================================
  13 +
  14 + Add logL and event statistics to cspull output file
  15 +
  16 +
9 17 2014-02-05 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
10 18  
11 19 * Bug fix Version 0.7.1 released
... ...
... ... @@ -2,7 +2,7 @@ New Features and Important Changes in ctools 00-08-00
2 2  
3 3 Author(s) : Juergen Knoedlseder
4 4  
5   -5 February 2014
  5 +10 February 2014
6 6  
7 7  
8 8 1. Introduction
... ... @@ -13,13 +13,7 @@ important changes made since the last release of ctools.
13 13  
14 14 2. Configuration
15 15 -----------------
16   -Some hyphens were missing in the configuration script when synchronising
17   -the user and system parameter files. This led to an error message on the
18   -console, but did not impact the proper configuration. Hyphens have now
19   -been added as needed.
20   -
21   -The configuration test script testconf.sh has been added to allow testing
22   -configuration script during continous integration.
  16 +None
23 17  
24 18  
25 19 3. Python interface
... ... @@ -29,13 +23,7 @@ None
29 23  
30 24 4. ctobssim - CTA observation simulation
31 25 ----------------------------------------
32   -The instrument and id attributes provided by an observation were not
33   -considering when simulating source and background events. This has been
34   -changed now that only events will be simulated when the instrument
35   -and id are valid for a given model.
36   -
37   -ctobssim now has full access to the calibration database, allowing
38   -the use for 2D response functions in FITS files.
  26 +None
39 27  
40 28 ctobssim now also supports the handling of energy dispersion. Energy
41 29 dispersion can be switched on or off using the hidden edisp parameter.
... ... @@ -81,8 +69,9 @@ None
81 69  
82 70 11. Scripts
83 71 -----------
84   -The obsutils.py script has been adapted to the new calibration database
85   -interface.
  72 +Add logL and event statistics to cspull output file.
  73 +
  74 +Convert all tabs into spaces.
86 75  
87 76  
88 77 12. caldb
... ...
scripts/cspull.py
  View file @64e71e4
... ... @@ -30,54 +30,55 @@ import math
30 30 # cspull class #
31 31 # ============ #
32 32 class cspull(GApplication):
33   - """
34   - This class implements the pull distribution generation script. It derives
35   - from the GammaLib::GApplication class which provides support for parameter
36   - files, command line arguments, and logging. In that way the Python
37   - script behaves just as a regular ctool.
38   - """
39   - def __init__(self, *argv):
40   - """
41   - Constructor.
42   - """
43   - # Set name
44   - self.name = "cspull"
45   - self.version = "0.2.0"
46   -
47   - # Initialise some members
48   - self.obs = None
49   - self.model = None
50   - self.m_srcmdl = None
51   -
52   - # Make sure that parfile exists
53   - file = self.parfile()
  33 + """
  34 + This class implements the pull distribution generation script. It derives
  35 + from the GammaLib::GApplication class which provides support for parameter
  36 + files, command line arguments, and logging. In that way the Python
  37 + script behaves just as a regular ctool.
  38 + """
  39 + def __init__(self, *argv):
  40 + """
  41 + Constructor.
  42 + """
  43 + # Set name
  44 + self.name = "cspull"
  45 + self.version = "0.2.0"
  46 +
  47 + # Initialise some members
  48 + self.obs = None
  49 + self.model = None
  50 + self.m_srcmdl = None
  51 +
  52 + # Make sure that parfile exists
  53 + file = self.parfile()
54 54  
55   - # Initialise application
56   - if len(argv) == 0:
57   - GApplication.__init__(self, self.name, self.version)
58   - elif len(argv) ==1:
59   - GApplication.__init__(self, self.name, self.version, *argv)
60   - else:
61   - raise TypeError("Invalid number of arguments given.")
  55 + # Initialise application
  56 + if len(argv) == 0:
  57 + GApplication.__init__(self, self.name, self.version)
  58 + elif len(argv) ==1:
  59 + GApplication.__init__(self, self.name, self.version, *argv)
  60 + else:
  61 + raise TypeError("Invalid number of arguments given.")
62 62  
63   - # Set logger properties
64   - self.log_header()
65   - self.log.date(True)
  63 + # Set logger properties
  64 + self.log_header()
  65 + self.log.date(True)
66 66  
67   - # Return
68   - return
69   -
70   - def __del__(self):
71   - """
72   - Destructor.
73   - """
74   - # Write separator into logger
75   - if self.logTerse():
76   - self.log("\n")
77   -
78   - # Return
79   - return
  67 + # Return
  68 + return
  69 +
  70 + def __del__(self):
  71 + """
  72 + Destructor.
  73 + """
  74 + # Write separator into logger
  75 + if self.logTerse():
  76 + self.log("\n")
  77 +
  78 + # Return
  79 + return
80 80  
  81 +<<<<<<< HEAD
81 82 def parfile(self):
82 83 """
83 84 Check if parfile exists. If parfile does not exist then create a
... ... @@ -138,132 +139,193 @@ class cspull(GApplication):
138 139 self.m_rad = self["rad"].real()
139 140 self.m_npix = self["npix"].integer()
140 141 self.m_binsz = self["binsz"].real()
  142 +=======
  143 + def parfile(self):
  144 + """
  145 + Check if parfile exists. If parfile does not exist then create a
  146 + default parfile. This kluge avoids shipping the cscript with a parfile.
  147 + """
  148 + # Set parfile name
  149 + parfile = self.name+".par"
  150 +
  151 + try:
  152 + pars = GApplicationPars(parfile)
  153 + except:
  154 + # Signal if parfile was not found
  155 + sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
  156 +
  157 + # Create default parfile
  158 + pars = GApplicationPars()
  159 + pars.append(GApplicationPar("srcmdl","f","a","$GAMMALIB/share/models/crab.xml","","","Source model"))
  160 + pars.append(GApplicationPar("outfile","f","a","pull.dat","","","Output file name"))
  161 + pars.append(GApplicationPar("ntrials","i","a","10","","","Number of trials"))
  162 + pars.append(GApplicationPar("caldb","s","h","","","","Calibration database"))
  163 + pars.append(GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function"))
  164 + pars.append(GApplicationPar("ra","r","a","83.6331","0","360","RA of pointing (deg)"))
  165 + pars.append(GApplicationPar("dec","r","a","22.0145","-90","90","Dec of pointing (deg)"))
  166 + pars.append(GApplicationPar("emin","r","a","0.1","0.0","","Lower energy limit (TeV)"))
  167 + pars.append(GApplicationPar("emax","r","a","100.0","0.0","","Upper energy limit (TeV)"))
  168 + pars.append(GApplicationPar("enumbins","i","a","0","","","Number of energy bins (0=unbinned)"))
  169 + pars.append(GApplicationPar("duration","r","a","180000.0","","","Effective exposure time (s)"))
  170 + pars.append(GApplicationPar("deadc","r","h","0.95","","","Deadtime correction factor"))
  171 + pars.append(GApplicationPar("rad","r","h","5.0","","","Radius of ROI (deg)"))
  172 + pars.append(GApplicationPar("npix","i","h","200","","","Number of pixels for binned"))
  173 + pars.append(GApplicationPar("binsz","r","h","0.05","","","Pixel size for binned (deg/pixel)"))
  174 + pars.append_standard()
  175 + pars.save(parfile)
  176 +
  177 + # Return
  178 + return
  179 +
  180 + def get_parameters(self):
  181 + """
  182 + Get parameters from parfile and setup the observation.
  183 + """
  184 + # Get parameters
  185 + if self.model == None:
  186 + self.m_srcmdl = self["srcmdl"].filename()
  187 + self.m_outfile = self["outfile"].filename()
  188 + self.m_ntrials = self["ntrials"].integer()
  189 + self.m_caldb = self["caldb"].string()
  190 + self.m_irf = self["irf"].string()
  191 + self.m_ra = self["ra"].real()
  192 + self.m_dec = self["dec"].real()
  193 + self.m_emin = self["emin"].real()
  194 + self.m_emax = self["emax"].real()
  195 + self.m_enumbins = self["enumbins"].integer()
  196 + self.m_duration = self["duration"].real()
  197 + self.m_deadc = self["deadc"].real()
  198 + self.m_rad = self["rad"].real()
  199 + self.m_npix = self["npix"].integer()
  200 + self.m_binsz = self["binsz"].real()
  201 +>>>>>>> 00-07-02
141 202  
142   - # Set some fixed parameters
143   - self.m_log = False # Logging in client tools
144   - self.m_debug = False # Debugging in client tools
145   -
146   - # Setup observation
147   - self.obs = self.set_obs(ra=self.m_ra, dec=self.m_dec, \
148   - emin=self.m_emin, emax=self.m_emax)
149   -
150   - # Load source model
151   - if self.m_srcmdl != None:
152   - self.model = GModels(self.m_srcmdl)
153   -
154   - # Append source model to observation
155   - self.obs.models(self.model)
  203 + # Set some fixed parameters
  204 + self.m_log = False # Logging in client tools
  205 + self.m_debug = False # Debugging in client tools
  206 +
  207 + # Setup observation
  208 + self.obs = self.set_obs(ra=self.m_ra, dec=self.m_dec, \
  209 + emin=self.m_emin, emax=self.m_emax)
  210 +
  211 + # Load source model
  212 + if self.m_srcmdl != None:
  213 + self.model = GModels(self.m_srcmdl)
  214 +
  215 + # Append source model to observation
  216 + self.obs.models(self.model)
156 217  
157   - # Return
158   - return
159   -
160   - def models(self, models):
161   - """
162   - Set model.
163   - """
164   - # Copy models
165   - self.model = models.copy()
166   -
167   - # Return
168   - return
169   -
170   - def execute(self):
171   - """
172   - Execute the script.
173   - """
174   - # Run the script
175   - self.run()
176   -
177   - # Return
178   - return
  218 + # Return
  219 + return
  220 +
  221 + def models(self, models):
  222 + """
  223 + Set model.
  224 + """
  225 + # Copy models
  226 + self.model = models.copy()
  227 +
  228 + # Return
  229 + return
  230 +
  231 + def execute(self):
  232 + """
  233 + Execute the script.
  234 + """
  235 + # Run the script
  236 + self.run()
  237 +
  238 + # Return
  239 + return
179 240  
180   - def run(self):
181   - """
182   - Run the script.
183   - """
184   - # Switch screen logging on in debug mode
185   - if self.logDebug():
186   - self.log.cout(True)
  241 + def run(self):
  242 + """
  243 + Run the script.
  244 + """
  245 + # Switch screen logging on in debug mode
  246 + if self.logDebug():
  247 + self.log.cout(True)
187 248  
188   - # Get parameters
189   - self.get_parameters()
190   -
191   - # Write input parameters into logger
192   - if self.logTerse():
193   - self.log_parameters()
194   - self.log("\n")
195   -
196   - # Write observation into logger
197   - if self.logTerse():
198   - self.log("\n")
199   - self.log.header1("Observation")
200   - self.log(str(self.obs))
201   - self.log("\n")
  249 + # Get parameters
  250 + self.get_parameters()
  251 +
  252 + # Write input parameters into logger
  253 + if self.logTerse():
  254 + self.log_parameters()
  255 + self.log("\n")
  256 +
  257 + # Write observation into logger
  258 + if self.logTerse():
  259 + self.log("\n")
  260 + self.log.header1("Observation")
  261 + self.log(str(self.obs))
  262 + self.log("\n")
202 263  
203   - # Write header
204   - if self.logTerse():
205   - self.log("\n")
206   - self.log.header1("Generate pull distribution")
  264 + # Write header
  265 + if self.logTerse():
  266 + self.log("\n")
  267 + self.log.header1("Generate pull distribution")
207 268  
208   - # Loop over trials
209   - for seed in range(self.m_ntrials):
210   -
211   - # Make a trial
212   - result = self.trial(seed)
213   -
214   - # Write out result immediately
215   - if seed == 0:
216   - file = open(self.m_outfile, 'w')
217   - writer = csv.DictWriter(file, result['colnames'])
218   - writer.writerow(dict((_,_) for _ in result['colnames']))
219   - else:
220   - file = open(self.m_outfile, 'a')
221   - writer = csv.DictWriter(file, result['colnames'])
222   - writer.writerow(result['values'])
223   - file.close()
224   -
225   - # Return
226   - return
227   -
228   - def set_obs(self, ra=83.6331, dec=22.0145, emin=0.1, emax=100.0):
229   - """
230   - Returns an observation container with a single CTA observation.
231   -
232   - Keywords:
233   - ra - Right Ascension of pointing [deg] (default: 83.6331)
234   - dec - Declination of pointing [deg] (default: 22.0145)
235   - emin - Minimum energy [TeV] (default: 0.1)
236   - emax - Maximum energy [TeV] (default: 100.0)
237   - """
238   - # Allocate observation container
239   - obs = GObservations()
240   -
241   - # Set single pointing
242   - pntdir = GSkyDir()
243   - pntdir.radec_deg(ra, dec)
244   -
245   - # Create CTA observation
246   - run = obsutils.set(pntdir, caldb=self.m_caldb, irf=self.m_irf, \
247   - duration=self.m_duration, deadc=self.m_deadc, \
  269 + # Loop over trials
  270 + for seed in range(self.m_ntrials):
  271 +
  272 + # Make a trial
  273 + result = self.trial(seed)
  274 +
  275 + # Write out result immediately
  276 + if seed == 0:
  277 + file = open(self.m_outfile, 'w')
  278 + writer = csv.DictWriter(file, result['colnames'])
  279 + writer.writerow(dict((_,_) for _ in result['colnames']))
  280 + else:
  281 + file = open(self.m_outfile, 'a')
  282 + writer = csv.DictWriter(file, result['colnames'])
  283 + writer.writerow(result['values'])
  284 + file.close()
  285 +
  286 + # Return
  287 + return
  288 +
  289 + def set_obs(self, ra=83.6331, dec=22.0145, emin=0.1, emax=100.0):
  290 + """
  291 + Returns an observation container with a single CTA observation.
  292 +
  293 + Keywords:
  294 + ra - Right Ascension of pointing [deg] (default: 83.6331)
  295 + dec - Declination of pointing [deg] (default: 22.0145)
  296 + emin - Minimum energy [TeV] (default: 0.1)
  297 + emax - Maximum energy [TeV] (default: 100.0)
  298 + """
  299 + # Allocate observation container
  300 + obs = GObservations()
  301 +
  302 + # Set single pointing
  303 + pntdir = GSkyDir()
  304 + pntdir.radec_deg(ra, dec)
  305 +
  306 + # Create CTA observation
  307 + run = obsutils.set(pntdir, caldb=self.m_caldb, irf=self.m_irf, \
  308 + duration=self.m_duration, deadc=self.m_deadc, \
248 309 emin=emin, emax=emax, rad=self.m_rad)
249   -
250   - # Append observation to container
251   - obs.append(run)
252   -
253   - # Return observation container
254   - return obs
255   -
256   - def trial(self, seed):
257   - """
258   - Create the pull for a single trial.
259   -
260   - Parameters:
261   - seed - Random number generator seed
262   - """
263   - # Write header
264   - if self.logExplicit():
265   - self.log.header2("Trial "+str(seed+1))
  310 +
  311 + # Append observation to container
  312 + obs.append(run)
  313 +
  314 + # Return observation container
  315 + return obs
  316 +
  317 + def trial(self, seed):
  318 + """
  319 + Create the pull for a single trial.
  320 +
  321 + Parameters:
  322 + seed - Random number generator seed
  323 + """
  324 + # Write header
  325 + if self.logExplicit():
  326 + self.log.header2("Trial "+str(seed+1))
266 327  
  328 +<<<<<<< HEAD
267 329 # Simulate events
268 330 obs = obsutils.sim(self.obs, \
269 331 nbins=self.m_enumbins, \
... ... @@ -272,94 +334,114 @@ class cspull(GApplication):
272 334 npix=self.m_npix, \
273 335 edisp=self.m_edisp, \
274 336 log=self.m_log, debug=self.m_debug)
  337 +=======
  338 + # Simulate events
  339 + obs = obsutils.sim(self.obs, \
  340 + nbins=self.m_enumbins, \
  341 + seed=seed, \
  342 + binsz=self.m_binsz, \
  343 + npix=self.m_npix, \
  344 + log=self.m_log, debug=self.m_debug)
  345 +>>>>>>> 00-07-02
275 346  
276   - # Determine number of events in simulation
277   - nevents = 0.0
278   - for run in obs:
279   - nevents += run.events().number()
  347 + # Determine number of events in simulation
  348 + nevents = 0.0
  349 + for run in obs:
  350 + nevents += run.events().number()
280 351  
281   - # Write simulation results
282   - if self.logExplicit():
283   - self.log.header3("Simulation")
284   - self.log.parformat("Number of simulated events")
285   - self.log(nevents)
286   - self.log("\n")
  352 + # Write simulation results
  353 + if self.logExplicit():
  354 + self.log.header3("Simulation")
  355 + self.log.parformat("Number of simulated events")
  356 + self.log(nevents)
  357 + self.log("\n")
287 358  
  359 +<<<<<<< HEAD
288 360 # Fit model
289 361 like = obsutils.fit(obs, edisp=self.m_edisp, \
290 362 log=self.m_log, debug=self.m_debug)
  363 +=======
  364 + # Fit model
  365 + like = obsutils.fit(obs, log=self.m_log, debug=self.m_debug)
  366 +>>>>>>> 00-07-02
291 367  
292   - # Store results
293   - logL = like.opt().value()
294   - npred = like.obs().npred()
295   - models = like.obs().models()
  368 + # Store results
  369 + logL = like.opt().value()
  370 + npred = like.obs().npred()
  371 + models = like.obs().models()
296 372  
297   - # Write result header
298   - if self.logExplicit():
299   - self.log.header3("Pulls")
300   -
301   - # Gather results
302   - colnames = []
303   - values = {}
304   - for i in range(models.size()):
305   - model = models[i]
306   - model_name = model.name()
307   - for k in range(model.size()):
308   - par = model[k]
309   - if par.is_free():
310   -
311   - # Set parameter name
312   - name = model_name+"_"+par.name()
313   -
314   - # Append parameter, Pull_parameter and Unc_parameter
315   - colnames.append(name)
316   - colnames.append("Pull_"+name)
317   - colnames.append("Unc_"+name)
318   -
319   - # Compute pull
320   - fitted_value = par.value()
321   - real_value = self.model[i][k].value()
322   - error = par.error()
323   - if error != 0.0:
324   - pull = (fitted_value - real_value) / error
325   - else:
326   - pull = 99.0
327   -
328   - # Store results
329   - values[name] = fitted_value
330   - values["Pull_"+name] = pull
331   - values["Unc_"+name] = error
  373 + # Write result header
  374 + if self.logExplicit():
  375 + self.log.header3("Pulls")
  376 +
  377 + # Gather results
  378 + colnames = []
  379 + values = {}
  380 + colnames.append("LogL")
  381 + colnames.append("Sim_Events")
  382 + colnames.append("Npred_Events")
  383 + values["LogL"] = logL
  384 + values["Sim_Events"] = nevents
  385 + values["Npred_Events"] = npred
  386 + for i in range(models.size()):
  387 + model = models[i]
  388 + model_name = model.name()
  389 + for k in range(model.size()):
  390 + par = model[k]
  391 + if par.is_free():
  392 +
  393 + # Set parameter name
  394 + name = model_name+"_"+par.name()
  395 +
  396 + # Append parameter, Pull_parameter and Unc_parameter
  397 + colnames.append(name)
  398 + colnames.append("Pull_"+name)
  399 + colnames.append("Unc_"+name)
  400 +
  401 + # Compute pull
  402 + fitted_value = par.value()
  403 + real_value = self.model[i][k].value()
  404 + error = par.error()
  405 + if error != 0.0:
  406 + pull = (fitted_value - real_value) / error
  407 + else:
  408 + pull = 99.0
  409 +
  410 + # Store results
  411 + values[name] = fitted_value
  412 + values["Pull_"+name] = pull
  413 + values["Unc_"+name] = error
332 414  
333   - # Write result
334   - if self.logExplicit():
335   - self.log.parformat(name)
336   - self.log(pull)
337   - self.log(" (")
338   - self.log(fitted_value)
339   - self.log(" +/- ")
340   - self.log(error)
341   - self.log(")\n")
342   -
343   - # Bundle together results
344   - result = {'colnames': colnames, 'values': values}
345   -
346   - # Return
347   - return result
  415 + # Write result
  416 + if self.logExplicit():
  417 + self.log.parformat(name)
  418 + self.log(pull)
  419 + self.log(" (")
  420 + self.log(fitted_value)
  421 + self.log(" +/- ")
  422 + self.log(error)
  423 + self.log(")\n")
  424 +
  425 + # Bundle together results
  426 + result = {'colnames': colnames, 'values': values}
  427 +
  428 + # Return
  429 + return result
348 430  
349 431  
350 432 # ======================== #
351 433 # Main routine entry point #
352 434 # ======================== #
353 435 if __name__ == '__main__':
354   - """
355   - Generates pull distribution for a source model.
356   - """
357   - # Create instance of application
358   - app = cspull(sys.argv)
359   -
360   - # Open logfile
361   - app.logFileOpen()
362   -
363   - # Execute application
364   - app.execute()
365   -
366 436 \ No newline at end of file
  437 + """
  438 + Generates pull distribution for a source model.
  439 + """
  440 + # Create instance of application
  441 + app = cspull(sys.argv)
  442 +
  443 + # Open logfile
  444 + app.logFileOpen()
  445 +
  446 + # Execute application
  447 + app.execute()
  448 +
... ...
scripts/cssens.py
  View file @64e71e4
... ... @@ -17,7 +17,7 @@
17 17 #
18 18 # The source is modelled as a point source
19 19 #
20   -# Copyright (C) 2011 Jurgen Knodlseder
  20 +# Copyright (C) 2011-2014 Juergen Knoedlseder
21 21 #
22 22 # This program is free software: you can redistribute it and/or modify
23 23 # it under the terms of the GNU General Public License as published by
... ... @@ -45,682 +45,682 @@ import math
45 45 # cssens class #
46 46 # ============ #
47 47 class cssens(GApplication):
48   - """
49   - This class implements the sensitivity computation script. It derives
50   - from the GammaLib::GApplication class which provides support for parameter
51   - files, command line arguments, and logging. In that way the Python
52   - script behaves just as a regular ctool.
53   - """
54   - def __init__(self, *argv):
55   - """
56   - Constructor.
57   - """
58   - # Set name
59   - self.name = "cssens"
60   - self.version = "0.3.0"
61   -
62   - # Make sure that parfile exists
63   - file = self.parfile()
64   -
65   - # Initialise application
66   - if len(argv) == 0:
67   - GApplication.__init__(self, self.name, self.version)
68   - elif len(argv) ==1:
69   - GApplication.__init__(self, self.name, self.version, *argv)
70   - else:
71   - raise TypeError("Invalid number of arguments given.")
72   -
73   - # Set logger properties
74   - self.logFileOpen()
75   - self.log_header()
76   - self.log.date(True)
77   -
78   - # Return
79   - return
80   -
81   - def __del__(self):
82   - """
83   - Destructor.
84   - """
85   - # Write separator into logger
86   - if self.logTerse():
87   - self.log("\n")
88   -
89   - # Return
90   - return
91   -
92   - def parfile(self):
93   - """
94   - Check if parfile exists. If parfile does not exist then create a
95   - default parfile. This kluge avoids shipping the cscript with a parfile.
96   - """
97   - # Set parfile name
98   - parfile = self.name+".par"
99   -
100   - try:
101   - pars = GApplicationPars(parfile)
102   - except:
103   - # Signal if parfile was not found
104   - sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
105   -
106   - # Create default parfile
107   - pars = GApplicationPars()
108   - pars.append(GApplicationPar("outfile","f","h","sensitivity.dat","","","Output file name"))
109   - pars.append(GApplicationPar("caldb","s","a","$GAMMALIB/share/caldb/cta","","","Calibration database"))
110   - pars.append(GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function"))
111   - pars.append(GApplicationPar("type","s","a","point","","","Source model type (point/gauss/shell/disk)"))
112   - pars.append(GApplicationPar("offset","r","a","0.0","0.0","","Source offset angle (deg)"))
113   - pars.append(GApplicationPar("bkg","s","a","$GAMMALIB/share/models/bkg_dummy.txt","","","Background model file function (none=power law for E)"))
114   - pars.append(GApplicationPar("duration","r","a","180000.0","","","Effective exposure time (s)"))
115   - pars.append(GApplicationPar("rad","r","a","5.0","","","Radius of ROI (deg)"))
116   - pars.append(GApplicationPar("enumbins","i","h","0","","","Number of energy bins (0=unbinned)"))
117   - pars.append(GApplicationPar("npix","i","h","200","","","Number of pixels for binned"))
118   - pars.append(GApplicationPar("binsz","r","h","0.05","","","Pixel size for binned (deg/pixel)"))
119   - pars.append(GApplicationPar("sigma","r","h","5.0","","","Significance threshold"))
120   - pars.append(GApplicationPar("ts_use","b","h","yes","","","Use TS?"))
121   - pars.append(GApplicationPar("index","r","h","-2.48","","","Assumed spectral index"))
122   - pars.append(GApplicationPar("radius","r","h","0.1","","","Extended source model radius"))
123   - pars.append(GApplicationPar("width","r","h","0.05","","","Extended source model width"))
124   - pars.append(GApplicationPar("max_iter","i","h","50","","","Maximum number of iterations"))
125   - pars.append(GApplicationPar("num_avg","i","h","3","","","Number of iterations for sliding average"))
126   - pars.append_standard()
127   - pars.save(parfile)
128   -
129   - # Return
130   - return
131   -
132   - def get_parameters(self):
133   - """
134   - Get parameters from parfile.
135   - """
136   - # Get parameters
137   - self.m_outfile = self["outfile"].filename()
138   - self.m_caldb = self["caldb"].string()
139   - self.m_irf = self["irf"].string()
140   - self.m_type = self["type"].string()
141   - self.m_offset = self["offset"].real()
142   - self.m_bkg = self["bkg"].string()
143   - self.m_duration = self["duration"].real()
144   - self.m_roi = self["rad"].real()
145   - self.m_nbins = self["enumbins"].integer()
146   - self.m_npix = self["npix"].integer()
147   - self.m_binsz = self["binsz"].real()
148   - self.m_ts_thres = self["sigma"].real()*self["sigma"].real()
149   - self.m_use_ts = self["ts_use"].boolean()
150   - self.m_index = self["index"].real()
151   - self.m_radius = self["radius"].real()
152   - self.m_width = self["width"].real()
153   - self.m_max_iter = self["max_iter"].integer()
154   - self.m_num_avg = self["num_avg"].integer()
155   -
156   - # Set some fixed parameters
157   - self.m_log = False # Logging in client tools
158   - self.m_debug = False # Debugging in client tools
159   -
160   - # Return
161   - return
162   -
163   - def execute(self):
164   - """
165   - Execute the script.
166   - """
167   - # Run the script
168   - self.run()
169   -
170   - # Return
171   - return
172   -
173   - def run(self):
174   - """
175   - Run the script.
176   - """
177   - # Switch screen logging on in debug mode
178   - if self.logDebug():
179   - self.log.cout(True)
180   -
181   - # Get parameters
182   - self.get_parameters()
183   -
184   - # Write input parameters into logger
185   - if self.logTerse():
186   - self.log_parameters()
187   - self.log("\n")
188   -
189   - # Initialise script
190   - colnames = ['loge', 'emin', 'emax', 'crab', 'pflux', 'eflux', 'diffSens']
191   - results = []
192   -
193   - # Initialise models. Note that we centre the point source at the Galactic
194   - # center as our observation is also centred at the Galactic centre, so
195   - # we're onaxis.
196   - bkg_model = GModels()
197   - full_model = GModels()
198   - bkg_model.append(self.set_bkg_model())
199   - full_model.append(self.set_bkg_model())
200   - full_model.append(self.set_src_model(0.0, \
201   - self.m_offset, \
202   - type=self.m_type, \
203   - index=self.m_index))
204   -
205   - # Write models into logger
206   - if self.logTerse():
207   - self.log("\n")
208   - self.log.header1("Models")
209   - self.log.header2("Background model")
210   - self.log(str(bkg_model))
211   - self.log("\n\n")
212   - self.log.header2("Full model")
213   - self.log(str(full_model))
214   - self.log("\n")
215   -
216   - # Write heder
217   - if self.logTerse():
218   - self.log("\n")
219   - self.log.header1("Sensitivity determination")
220   -
221   - # Loop over energy bands. The number of energy band is still
222   - # fixed ...
223   - for ieng in range(21):
224   -
225   - # Set energies
226   - loge = -1.7 + ieng * 0.2
227   - emean = pow(10.0, loge)
228   - emin = pow(10.0, loge-0.1)
229   - emax = pow(10.0, loge+0.1)
230   -
231   - # Setup observation(s)
232   - obs = self.set_obs(emin=emin, emax=emax)
233   -
234   - # Determine sensitivity
235   - result = self.get_sensitivity(obs, bkg_model, full_model)
236   -
237   - # Write results
238   - if ieng == 0:
239   - f = open(self.m_outfile, 'w')
240   - writer = csv.DictWriter(f, colnames)
241   - writer.writerow(dict((_,_) for _ in colnames))
242   - writer.writerow(result)
243   - f.close()
244   - else:
245   - f = open(self.m_outfile, 'a')
246   - writer = csv.DictWriter(f, colnames)
247   - writer.writerow(result)
248   - f.close()
249   -
250   - # Append results
251   - results.append(result)
252   -
253   - # Return
254   - return
255   -
256   - def set_obs(self, lpnt=0.0, bpnt=0.0, emin=0.1, emax=100.0):
257   - """
258   - Returns an observation container with a single CTA observation.
259   -
260   - Keywords:
261   - lpnt - Galactic longitude of pointing [deg] (default: 0.0)
262   - bpnt - Galactic latitude of pointing [deg] (default: 0.0)
263   - emin - Minimum energy [TeV] (default: 0.1)
264   - emax - Maximum energy [TeV] (default: 100.0)
265   - """
266   - # Allocate observation container
267   - obs = GObservations()
268   -
269   - # Set single pointing
270   - pntdir = GSkyDir()
271   - pntdir.lb_deg(lpnt, bpnt)
272   -
273   - # Create CTA observation
274   - run = obsutils.set(pntdir, caldb=self.m_caldb, irf=self.m_irf, \
275   - duration=self.m_duration, \
  48 + """
  49 + This class implements the sensitivity computation script. It derives
  50 + from the GammaLib::GApplication class which provides support for parameter
  51 + files, command line arguments, and logging. In that way the Python
  52 + script behaves just as a regular ctool.
  53 + """
  54 + def __init__(self, *argv):
  55 + """
  56 + Constructor.
  57 + """
  58 + # Set name
  59 + self.name = "cssens"
  60 + self.version = "0.3.0"
  61 +
  62 + # Make sure that parfile exists
  63 + file = self.parfile()
  64 +
  65 + # Initialise application
  66 + if len(argv) == 0:
  67 + GApplication.__init__(self, self.name, self.version)
  68 + elif len(argv) ==1:
  69 + GApplication.__init__(self, self.name, self.version, *argv)
  70 + else:
  71 + raise TypeError("Invalid number of arguments given.")
  72 +
  73 + # Set logger properties
  74 + self.logFileOpen()
  75 + self.log_header()
  76 + self.log.date(True)
  77 +
  78 + # Return
  79 + return
  80 +
  81 + def __del__(self):
  82 + """
  83 + Destructor.
  84 + """
  85 + # Write separator into logger
  86 + if self.logTerse():
  87 + self.log("\n")
  88 +
  89 + # Return
  90 + return
  91 +
  92 + def parfile(self):
  93 + """
  94 + Check if parfile exists. If parfile does not exist then create a
  95 + default parfile. This kluge avoids shipping the cscript with a parfile.
  96 + """
  97 + # Set parfile name
  98 + parfile = self.name+".par"
  99 +
  100 + try:
  101 + pars = GApplicationPars(parfile)
  102 + except:
  103 + # Signal if parfile was not found
  104 + sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
  105 +
  106 + # Create default parfile
  107 + pars = GApplicationPars()
  108 + pars.append(GApplicationPar("outfile","f","h","sensitivity.dat","","","Output file name"))
  109 + pars.append(GApplicationPar("caldb","s","a","$GAMMALIB/share/caldb/cta","","","Calibration database"))
  110 + pars.append(GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function"))
  111 + pars.append(GApplicationPar("type","s","a","point","","","Source model type (point/gauss/shell/disk)"))
  112 + pars.append(GApplicationPar("offset","r","a","0.0","0.0","","Source offset angle (deg)"))
  113 + pars.append(GApplicationPar("bkg","s","a","$GAMMALIB/share/models/bkg_dummy.txt","","","Background model file function (none=power law for E)"))
  114 + pars.append(GApplicationPar("duration","r","a","180000.0","","","Effective exposure time (s)"))
  115 + pars.append(GApplicationPar("rad","r","a","5.0","","","Radius of ROI (deg)"))
  116 + pars.append(GApplicationPar("enumbins","i","h","0","","","Number of energy bins (0=unbinned)"))
  117 + pars.append(GApplicationPar("npix","i","h","200","","","Number of pixels for binned"))
  118 + pars.append(GApplicationPar("binsz","r","h","0.05","","","Pixel size for binned (deg/pixel)"))
  119 + pars.append(GApplicationPar("sigma","r","h","5.0","","","Significance threshold"))
  120 + pars.append(GApplicationPar("ts_use","b","h","yes","","","Use TS?"))
  121 + pars.append(GApplicationPar("index","r","h","-2.48","","","Assumed spectral index"))
  122 + pars.append(GApplicationPar("radius","r","h","0.1","","","Extended source model radius"))
  123 + pars.append(GApplicationPar("width","r","h","0.05","","","Extended source model width"))
  124 + pars.append(GApplicationPar("max_iter","i","h","50","","","Maximum number of iterations"))
  125 + pars.append(GApplicationPar("num_avg","i","h","3","","","Number of iterations for sliding average"))
  126 + pars.append_standard()
  127 + pars.save(parfile)
  128 +
  129 + # Return
  130 + return
  131 +
  132 + def get_parameters(self):
  133 + """
  134 + Get parameters from parfile.
  135 + """
  136 + # Get parameters
  137 + self.m_outfile = self["outfile"].filename()
  138 + self.m_caldb = self["caldb"].string()
  139 + self.m_irf = self["irf"].string()
  140 + self.m_type = self["type"].string()
  141 + self.m_offset = self["offset"].real()
  142 + self.m_bkg = self["bkg"].string()
  143 + self.m_duration = self["duration"].real()
  144 + self.m_roi = self["rad"].real()
  145 + self.m_nbins = self["enumbins"].integer()
  146 + self.m_npix = self["npix"].integer()
  147 + self.m_binsz = self["binsz"].real()
  148 + self.m_ts_thres = self["sigma"].real()*self["sigma"].real()
  149 + self.m_use_ts = self["ts_use"].boolean()
  150 + self.m_index = self["index"].real()
  151 + self.m_radius = self["radius"].real()
  152 + self.m_width = self["width"].real()
  153 + self.m_max_iter = self["max_iter"].integer()
  154 + self.m_num_avg = self["num_avg"].integer()
  155 +
  156 + # Set some fixed parameters
  157 + self.m_log = False # Logging in client tools
  158 + self.m_debug = False # Debugging in client tools
  159 +
  160 + # Return
  161 + return
  162 +
  163 + def execute(self):
  164 + """
  165 + Execute the script.
  166 + """
  167 + # Run the script
  168 + self.run()
  169 +
  170 + # Return
  171 + return
  172 +
  173 + def run(self):
  174 + """
  175 + Run the script.
  176 + """
  177 + # Switch screen logging on in debug mode
  178 + if self.logDebug():
  179 + self.log.cout(True)
  180 +
  181 + # Get parameters
  182 + self.get_parameters()
  183 +
  184 + # Write input parameters into logger
  185 + if self.logTerse():
  186 + self.log_parameters()
  187 + self.log("\n")
  188 +
  189 + # Initialise script
  190 + colnames = ['loge', 'emin', 'emax', 'crab', 'pflux', 'eflux', 'diffSens']
  191 + results = []
  192 +
  193 + # Initialise models. Note that we centre the point source at the Galactic
  194 + # center as our observation is also centred at the Galactic centre, so
  195 + # we're onaxis.
  196 + bkg_model = GModels()
  197 + full_model = GModels()
  198 + bkg_model.append(self.set_bkg_model())
  199 + full_model.append(self.set_bkg_model())
  200 + full_model.append(self.set_src_model(0.0, \
  201 + self.m_offset, \
  202 + type=self.m_type, \
  203 + index=self.m_index))
  204 +
  205 + # Write models into logger
  206 + if self.logTerse():
  207 + self.log("\n")
  208 + self.log.header1("Models")
  209 + self.log.header2("Background model")
  210 + self.log(str(bkg_model))
  211 + self.log("\n\n")
  212 + self.log.header2("Full model")
  213 + self.log(str(full_model))
  214 + self.log("\n")
  215 +
  216 + # Write heder
  217 + if self.logTerse():
  218 + self.log("\n")
  219 + self.log.header1("Sensitivity determination")
  220 +
  221 + # Loop over energy bands. The number of energy band is still
  222 + # fixed ...
  223 + for ieng in range(21):
  224 +
  225 + # Set energies
  226 + loge = -1.7 + ieng * 0.2
  227 + emean = pow(10.0, loge)
  228 + emin = pow(10.0, loge-0.1)
  229 + emax = pow(10.0, loge+0.1)
  230 +
  231 + # Setup observation(s)
  232 + obs = self.set_obs(emin=emin, emax=emax)
  233 +
  234 + # Determine sensitivity
  235 + result = self.get_sensitivity(obs, bkg_model, full_model)
  236 +
  237 + # Write results
  238 + if ieng == 0:
  239 + f = open(self.m_outfile, 'w')
  240 + writer = csv.DictWriter(f, colnames)
  241 + writer.writerow(dict((_,_) for _ in colnames))
  242 + writer.writerow(result)
  243 + f.close()
  244 + else:
  245 + f = open(self.m_outfile, 'a')
  246 + writer = csv.DictWriter(f, colnames)
  247 + writer.writerow(result)
  248 + f.close()
  249 +
  250 + # Append results
  251 + results.append(result)
  252 +
  253 + # Return
  254 + return
  255 +
  256 + def set_obs(self, lpnt=0.0, bpnt=0.0, emin=0.1, emax=100.0):
  257 + """
  258 + Returns an observation container with a single CTA observation.
  259 +
  260 + Keywords:
  261 + lpnt - Galactic longitude of pointing [deg] (default: 0.0)
  262 + bpnt - Galactic latitude of pointing [deg] (default: 0.0)
  263 + emin - Minimum energy [TeV] (default: 0.1)
  264 + emax - Maximum energy [TeV] (default: 100.0)
  265 + """
  266 + # Allocate observation container
  267 + obs = GObservations()
  268 +
  269 + # Set single pointing
  270 + pntdir = GSkyDir()
  271 + pntdir.lb_deg(lpnt, bpnt)
  272 +
  273 + # Create CTA observation
  274 + run = obsutils.set(pntdir, caldb=self.m_caldb, irf=self.m_irf, \
  275 + duration=self.m_duration, \
276 276 emin=emin, emax=emax, rad=self.m_roi)
277   -
278   - # Append observation to container
279   - obs.append(run)
280   -
281   - # Return observation container
282   - return obs
283   -
284   - def set_bkg_model(self, fitidx=False):
285   - """
286   - Setup standard CTA background model. We use a simple power law here,
287   - scaled to Konrad's E configuration performance table.
288   -
289   - Keywords:
290   - fitidx - Fit spectral index (default: False)
291   - """
292   - # Define radial component
293   - radial = GCTAModelRadialGauss(3.0)
294   -
295   - # Define spectral component. If a background model is given then
296   - # use a file function. Otherwise use a static power law.
297   - if len(self.m_bkg) > 0:
298   - spectrum = GModelSpectralFunc(self.m_bkg, 1.0)
299   - else:
300   - spectrum = GModelSpectralPlaw(61.8e-6, -1.85, GEnergy(1.0, "TeV"))
301   - if fitidx:
302   - spectrum["Index"].free()
303   - else:
304   - spectrum["Index"].fix()
305   -
306   - # Create background model
307   - model = GCTAModelRadialAcceptance(radial, spectrum)
308   - model.name("Background")
309   - model.instruments("CTA")
310   -
311   - # Return background model
312   - return model
313   -
314   - def set_src_model(self, l, b, flux=1.0, index=-2.48, \
  277 +
  278 + # Append observation to container
  279 + obs.append(run)
  280 +
  281 + # Return observation container
  282 + return obs
  283 +
  284 + def set_bkg_model(self, fitidx=False):
  285 + """
  286 + Setup standard CTA background model. We use a simple power law here,
  287 + scaled to Konrad's E configuration performance table.
  288 +
  289 + Keywords:
  290 + fitidx - Fit spectral index (default: False)
  291 + """
  292 + # Define radial component
  293 + radial = GCTAModelRadialGauss(3.0)
  294 +
  295 + # Define spectral component. If a background model is given then
  296 + # use a file function. Otherwise use a static power law.
  297 + if len(self.m_bkg) > 0:
  298 + spectrum = GModelSpectralFunc(self.m_bkg, 1.0)
  299 + else:
  300 + spectrum = GModelSpectralPlaw(61.8e-6, -1.85, GEnergy(1.0, "TeV"))
  301 + if fitidx:
  302 + spectrum["Index"].free()
  303 + else:
  304 + spectrum["Index"].fix()
  305 +
  306 + # Create background model
  307 + model = GCTAModelRadialAcceptance(radial, spectrum)
  308 + model.name("Background")
  309 + model.instruments("CTA")
  310 +
  311 + # Return background model
  312 + return model
  313 +
  314 + def set_src_model(self, l, b, flux=1.0, index=-2.48, \
315 315 type="point", fitpos=False, fitidx=False):
316   - """
317   - Returns a single source with Crab-like spectrum. The source flux
318   - can be scaled in Crab units. The Crab spectrum is based on MAGIC
319   - observations (Albert et al. 2008, ApJ, 674, 1037).
320   -
321   - Parameters:
322   - l - Galactic longitude of source location [deg]
323   - b - Galactic latitude of source location [deg]
324   - Keywords:
325   - flux - Source flux [Crabs]
326   - index - Spectral index
327   - type - Source type ("point", "gauss", "disk", "shell")
328   - fitpos - Fit position and size? (default: False)
329   - fitidx - Fit index? (default: False)
330   - """
331   - # Set source location
332   - location = GSkyDir()
333   - location.lb_deg(l, b)
334   -
335   - # Set source spectrum
336   - spectrum = GModelSpectralPlaw(flux*5.7e-16, index, GEnergy(0.3, "TeV")
337   - if fitidx:
338   - spectrum["Index"].free()
339   - else:
340   - spectrum["Index"].fix()
341   -
342   - # Set source
343   - if type == "point":
344   - spatial = GModelSpatialPointSource(location)
345   - if fitpos:
346   - spatial["RA"].free()
347   - spatial["DEC"].free()
348   - else:
349   - spatial["RA"].fix()
350   - spatial["DEC"].fix()
351   - elif type == "gauss":
352   - spatial = GModelSpatialRadialGauss(location, self.m_radius)
353   - if fitpos:
354   - spatial["RA"].free()
355   - spatial["DEC"].free()
356   - spatial["Sigma"].free()
357   - else:
358   - spatial["RA"].fix()
359   - spatial["DEC"].fix()
360   - spatial["Sigma"].fix()
361   - elif type == "disk":
362   - spatial = GModelSpatialRadialDisk(location, self.m_radius)
363   - if fitpos:
364   - spatial["RA"].free()
365   - spatial["DEC"].free()
366   - spatial["Radius"].free()
367   - else:
368   - spatial["RA"].fix()
369   - spatial["DEC"].fix()
370   - spatial["Radius"].fix()
371   - elif type == "shell":
372   - spatial = GModelSpatialRadialShell(location, self.m_radius, self.m_width)
373   - if fitpos:
374   - spatial["RA"].free()
375   - spatial["DEC"].free()
376   - spatial["Radius"].free()
377   - spatial["Width"].free()
378   - else:
379   - spatial["RA"].fix()
380   - spatial["DEC"].fix()
381   - spatial["Radius"].fix()
382   - spatial["Width"].fix()
383   - else:
384   - self.log("ERROR: Unknown source type '"+type+"'.\n")
385   - return None
  316 + """
  317 + Returns a single source with Crab-like spectrum. The source flux
  318 + can be scaled in Crab units. The Crab spectrum is based on MAGIC
  319 + observations (Albert et al. 2008, ApJ, 674, 1037).
  320 +
  321 + Parameters:
  322 + l - Galactic longitude of source location [deg]
  323 + b - Galactic latitude of source location [deg]
  324 + Keywords:
  325 + flux - Source flux [Crabs]
  326 + index - Spectral index
  327 + type - Source type ("point", "gauss", "disk", "shell")
  328 + fitpos - Fit position and size? (default: False)
  329 + fitidx - Fit index? (default: False)
  330 + """
  331 + # Set source location
  332 + location = GSkyDir()
  333 + location.lb_deg(l, b)
  334 +
  335 + # Set source spectrum
  336 + spectrum = GModelSpectralPlaw(flux*5.7e-16, index, GEnergy(0.3, "TeV")
  337 + if fitidx:
  338 + spectrum["Index"].free()
  339 + else:
  340 + spectrum["Index"].fix()
  341 +
  342 + # Set source
  343 + if type == "point":
  344 + spatial = GModelSpatialPointSource(location)
  345 + if fitpos:
  346 + spatial["RA"].free()
  347 + spatial["DEC"].free()
  348 + else:
  349 + spatial["RA"].fix()
  350 + spatial["DEC"].fix()
  351 + elif type == "gauss":
  352 + spatial = GModelSpatialRadialGauss(location, self.m_radius)
  353 + if fitpos:
  354 + spatial["RA"].free()
  355 + spatial["DEC"].free()
  356 + spatial["Sigma"].free()
  357 + else:
  358 + spatial["RA"].fix()
  359 + spatial["DEC"].fix()
  360 + spatial["Sigma"].fix()
  361 + elif type == "disk":
  362 + spatial = GModelSpatialRadialDisk(location, self.m_radius)
  363 + if fitpos:
  364 + spatial["RA"].free()
  365 + spatial["DEC"].free()
  366 + spatial["Radius"].free()
  367 + else:
  368 + spatial["RA"].fix()
  369 + spatial["DEC"].fix()
  370 + spatial["Radius"].fix()
  371 + elif type == "shell":
  372 + spatial = GModelSpatialRadialShell(location, self.m_radius, self.m_width)
  373 + if fitpos:
  374 + spatial["RA"].free()
  375 + spatial["DEC"].free()
  376 + spatial["Radius"].free()
  377 + spatial["Width"].free()
  378 + else:
  379 + spatial["RA"].fix()
  380 + spatial["DEC"].fix()
  381 + spatial["Radius"].fix()
  382 + spatial["Width"].fix()
  383 + else:
  384 + self.log("ERROR: Unknown source type '"+type+"'.\n")
  385 + return None
386 386 source = GModelSky(spatial, spectrum)
387   -
388   - # Set source name
389   - source.name("Test")
390   -
391   - # Return source
392   - return source
393   -
394   - def get_sensitivity(self, obs, bkg_model, full_model):
395   - """
396   - Determine sensitivity for a given observations.
397   -
398   - Parameters:
399   - obs - Observation container
400   - bkg_model - Background model
401   - full_model - Source model
402   - """
403   - # Set TeV->erg conversion factor
404   - tev2erg = 1.6021764
405   -
406   - # Determine energy boundaries from first observation in
407   - # the container
408   - for run in obs:
409   - emin = run.events().ebounds().emin()
410   - emax = run.events().ebounds().emax()
411   - loge = math.log10(math.sqrt(emin.TeV()*emax.TeV()))
412   - erg_mean = math.pow(10.0, loge) * tev2erg
413   - erg_width = (emax.TeV()-emin.TeV()) * tev2erg
414   - break
415   -
416   - # Write header
417   - if self.logTerse():
418   - self.log("\n")
419   - self.log.header2("Energies: "+str(emin)+" - "+str(emax))
420   -
421   - # Initialise loop
422   - flux_value = []
423   - pflux_value = []
424   - eflux_value = []
425   - diffSens_value = []
426   - iter = 0
427   - test_flux = 0.1 # This is the initial test flux in Crab units
428   -
429   - # Loop until we break
430   - while True:
431   -
432   - # Update iteration counter
433   - iter += 1
434   -
435   - # Write header
436   - if self.logExplicit():
437   - self.log.header2("Iteration "+str(iter))
438   -
439   - # Set test flux
440   - full_model["Test"]['Prefactor'].value(test_flux)
441   - obs.models(full_model)
442   -
443   - # Simulate events
444   - sim = obsutils.sim(obs, nbins=self.m_nbins, seed=iter, \
445   - binsz=self.m_binsz, npix=self.m_npix, \
446   - log=self.m_log, debug=self.m_debug)
447   -
448   - # Determine number of events in simulation
449   - nevents = 0.0
450   - for run in sim:
451   - nevents += run.events().number()
452   -
453   - # Write simulation results
454   - if self.logExplicit():
455   - self.log.header3("Simulation")
456   - self.log.parformat("Number of simulated events")
457   - self.log(nevents)
458   - self.log("\n")
459   -
460   - # Fit background only
461   - sim.models(bkg_model)
462   - like = obsutils.fit(sim, log=self.m_log, debug=self.m_debug)
463   - result_bgm = like.obs().models().copy()
464   - LogL_bgm = like.opt().value()
465   - npred_bgm = like.obs().npred()
466   -
467   - # Assess quality based on a comparison between Npred and Nevents
468   - quality_bgm = npred_bgm-nevents
469   -
470   - # Write background fit results
471   - if self.logExplicit():
472   - self.log.header3("Background model fit")
473   - self.log.parformat("log likelihood")
474   - self.log(LogL_bgm)
475   - self.log("\n")
476   - self.log.parformat("Number of predicted events")
477   - self.log(npred_bgm)
478   - self.log("\n")
479   - self.log.parformat("Fit quality")
480   - self.log(quality_bgm)
481   - self.log("\n")
482   -
483   - # Start over if the fit quality was bad
484   - if abs(quality_bgm) > 3.0:
485   - if self.logExplicit():
486   - self.log("Fit quality outside required range. Start over.\n")
487   - continue
488   -
489   - # Write model fit results
490   - if self.logExplicit():
491   - for model in result_bgm:
492   - self.log.parformat("Model")
493   - self.log(model.name())
494   - self.log("\n")
495   - for par in model:
496   - self.log(str(par)+"\n")
497   -
498   - # Fit background and test source
499   - sim.models(full_model)
500   - like = obsutils.fit(sim, log=self.m_log, debug=self.m_debug)
501   - result_all = like.obs().models().copy()
502   - LogL_all = like.opt().value()
503   - npred_all = like.obs().npred()
504   - ts = 2.0*(LogL_bgm-LogL_all)
505   -
506   - # Assess quality based on a comparison between Npred and Nevents
507   - quality_all = npred_all-nevents
508   -
509   - # Write background and test source fit results
510   - if self.logExplicit():
511   - self.log.header3("Background and test source model fit")
512   - self.log.parformat("Test statistics")
513   - self.log(ts)
514   - self.log("\n")
515   - self.log.parformat("log likelihood")
516   - self.log(LogL_all)
517   - self.log("\n")
518   - self.log.parformat("Number of predicted events")
519   - self.log(npred_all)
520   - self.log("\n")
521   - self.log.parformat("Fit quality")
522   - self.log(quality_all)
523   - self.log("\n")
524   - #
525   - for model in result_all:
526   - self.log.parformat("Model")
527   - self.log(model.name())
528   - self.log("\n")
529   - for par in model:
530   - self.log(str(par)+"\n")
531   -
532   - # Start over if the fit quality was bad
533   - if abs(quality_all) > 3.0:
534   - if self.logExplicit():
535   - self.log("Fit quality outside required range. Start over.\n")
536   - continue
537   -
538   - # Start over if TS was non-positive
539   - if ts <= 0.0:
540   - if self.logExplicit():
541   - self.log("Non positive TS. Start over.\n")
542   - continue
543   -
544   - # Get fitted Crab, photon and energy fluxes
545   - cflux = result_all["Test"]['Prefactor'].value()
546   - cflux_err = result_all["Test"]['Prefactor'].error()
547   - pflux = result_all["Test"].spectral().flux(emin, emax)
548   - eflux = result_all["Test"].spectral().eflux(emin, emax)
549   -
550   - # Compute differential sensitivity
551   - diffSens = pflux / erg_width * erg_mean*erg_mean
552   -
553   - # Compute flux correction factor based on average TS
554   - correct = 1.0
555   - if self.m_use_ts:
556   - if ts > 0:
557   - correct = math.sqrt(self.m_ts_thres/ts)
558   - else:
559   - if cflux > 0 and cflux_err > 0:
560   - correct = math.sqrt(self.m_ts_thres)/(cflux/cflux_err)
561   -
562   - # Compute extrapolated fluxes
563   - flux = correct * cflux
564   - pflux = correct * pflux
565   - eflux = correct * eflux
566   - diffSens = correct * diffSens
567   - flux_value.append(flux)
568   - pflux_value.append(pflux)
569   - eflux_value.append(eflux)
570   - diffSens_value.append(diffSens)
571   -
572   - # Write background and test source fit results
573   - if self.logExplicit():
574   - self.log.parformat("Photon flux")
575   - self.log(pflux)
576   - self.log(" ph/cm2/s\n")
577   - self.log.parformat("Energy flux")
578   - self.log(eflux)
579   - self.log(" erg/cm2/s\n")
580   - self.log.parformat("Crab flux")
581   - self.log(cflux*1000.0)
582   - self.log(" mCrab\n")
583   - self.log.parformat("Differential sensitivity")
584   - self.log(diffSens)
585   - self.log(" erg/cm2/s\n")
586   - for model in result_all:
587   - self.log.parformat("Model")
588   - self.log(model.name())
589   - self.log("\n")
590   - for par in model:
591   - self.log(str(par)+"\n")
592   - elif self.logTerse():
593   - self.log.parformat("Iteration "+str(iter))
594   - self.log("TS=")
595   - self.log(ts)
596   - self.log(" ")
597   - self.log("corr=")
598   - self.log(correct)
599   - self.log(" ")
600   - self.log(pflux)
601   - self.log(" ph/cm2/s = ")
602   - self.log(eflux)
603   - self.log(" erg/cm2/s = ")
604   - self.log(cflux*1000.0)
605   - self.log(" mCrab = ")
606   - self.log(diffSens)
607   - self.log(" erg/cm2/s\n")
608   -
609   - # Compute sliding average of extrapolated fitted prefactor,
610   - # photon and energy flux. This damps out fluctuations and
611   - # improves convergence
612   - flux = 0.0
613   - pflux = 0.0
614   - eflux = 0.0
615   - diffSens = 0.0
616   - num = 0.0
617   - for k in range(self.m_num_avg):
618   - inx = len(flux_value) - k - 1
619   - if inx >= 0:
620   - flux += flux_value[inx]
621   - pflux += pflux_value[inx]
622   - eflux += eflux_value[inx]
623   - diffSens += diffSens_value[inx]
624   - num += 1.0
625   - flux /= num
626   - pflux /= num
627   - eflux /= num
628   - diffSens /= num
629   -
630   - # Compare average flux to last average
631   - if iter > self.m_num_avg:
632   - if test_flux > 0:
633   - ratio = flux/test_flux
634   -
635   - # We have 2 convergence criteria:
636   - # 1. The average flux does not change
637   - # 2. The flux correction factor is small
638   - if ratio >= 0.99 and ratio <= 1.01 and \
639   - correct >= 0.9 and correct <= 1.1:
640   - if self.logTerse():
641   - self.log(" Converged ("+str(ratio)+")\n")
642   - break
643   - else:
644   - if self.logTerse():
645   - self.log(" Flux is zero.\n")
646   - break
647   -
648   - # Use average for next iteration
649   - test_flux = flux
650   -
651   - # Exit loop if number of trials exhausted
652   - if (iter >= self.m_max_iter):
653   - if self.logTerse():
654   - self.log(" Test ended after "+str(self.m_max_iter)+" iterations.\n")
655   - break
656   -
657   - # Write fit results
658   - if self.logTerse():
659   - self.log.header3("Fit results")
660   - self.log.parformat("Test statistics")
661   - self.log(ts)
662   - self.log("\n")
663   - self.log.parformat("Photon flux")
664   - self.log(pflux)
665   - self.log(" ph/cm2/s\n")
666   - self.log.parformat("Energy flux")
667   - self.log(eflux)
668   - self.log(" erg/cm2/s\n")
669   - self.log.parformat("Crab flux")
670   - self.log(cflux*1000.0)
671   - self.log(" mCrab\n")
672   - self.log.parformat("Differential sensitivity")
673   - self.log(diffSens)
674   - self.log(" erg/cm2/s\n")
675   - self.log.parformat("Number of simulated events")
676   - self.log(nevents)
677   - self.log("\n")
678   - self.log.header3("Background and test source model fitting")
679   - self.log.parformat("log likelihood")
680   - self.log(LogL_all)
681   - self.log("\n")
682   - self.log.parformat("Number of predicted events")
683   - self.log(npred_all)
684   - self.log("\n")
685   - for model in result_all:
686   - self.log.parformat("Model")
687   - self.log(model.name())
688   - self.log("\n")
689   - for par in model:
690   - self.log(str(par)+"\n")
691   - self.log.header3("Background model fit")
692   - self.log.parformat("log likelihood")
693   - self.log(LogL_bgm)
694   - self.log("\n")
695   - self.log.parformat("Number of predicted events")
696   - self.log(npred_bgm)
697   - self.log("\n")
698   - for model in result_bgm:
699   - self.log.parformat("Model")
700   - self.log(model.name())
701   - self.log("\n")
702   - for par in model:
703   - self.log(str(par)+"\n")
704   -
705   - # Store result
706   - result = {'loge': loge, 'emin': emin.TeV(), 'emax': emax.TeV(), \
707   - 'crab': flux, 'pflux': pflux, 'eflux': eflux, \
708   - 'diffSens': diffSens}
709   -
710   - # Return result
711   - return result
  387 +
  388 + # Set source name
  389 + source.name("Test")
  390 +
  391 + # Return source
  392 + return source
  393 +
  394 + def get_sensitivity(self, obs, bkg_model, full_model):
  395 + """
  396 + Determine sensitivity for a given observations.
  397 +
  398 + Parameters:
  399 + obs - Observation container
  400 + bkg_model - Background model
  401 + full_model - Source model
  402 + """
  403 + # Set TeV->erg conversion factor
  404 + tev2erg = 1.6021764
  405 +
  406 + # Determine energy boundaries from first observation in
  407 + # the container
  408 + for run in obs:
  409 + emin = run.events().ebounds().emin()
  410 + emax = run.events().ebounds().emax()
  411 + loge = math.log10(math.sqrt(emin.TeV()*emax.TeV()))
  412 + erg_mean = math.pow(10.0, loge) * tev2erg
  413 + erg_width = (emax.TeV()-emin.TeV()) * tev2erg
  414 + break
  415 +
  416 + # Write header
  417 + if self.logTerse():
  418 + self.log("\n")
  419 + self.log.header2("Energies: "+str(emin)+" - "+str(emax))
  420 +
  421 + # Initialise loop
  422 + flux_value = []
  423 + pflux_value = []
  424 + eflux_value = []
  425 + diffSens_value = []
  426 + iter = 0
  427 + test_flux = 0.1 # This is the initial test flux in Crab units
  428 +
  429 + # Loop until we break
  430 + while True:
  431 +
  432 + # Update iteration counter
  433 + iter += 1
  434 +
  435 + # Write header
  436 + if self.logExplicit():
  437 + self.log.header2("Iteration "+str(iter))
  438 +
  439 + # Set test flux
  440 + full_model["Test"]['Prefactor'].value(test_flux)
  441 + obs.models(full_model)
  442 +
  443 + # Simulate events
  444 + sim = obsutils.sim(obs, nbins=self.m_nbins, seed=iter, \
  445 + binsz=self.m_binsz, npix=self.m_npix, \
  446 + log=self.m_log, debug=self.m_debug)
  447 +
  448 + # Determine number of events in simulation
  449 + nevents = 0.0
  450 + for run in sim:
  451 + nevents += run.events().number()
  452 +
  453 + # Write simulation results
  454 + if self.logExplicit():
  455 + self.log.header3("Simulation")
  456 + self.log.parformat("Number of simulated events")
  457 + self.log(nevents)
  458 + self.log("\n")
  459 +
  460 + # Fit background only
  461 + sim.models(bkg_model)
  462 + like = obsutils.fit(sim, log=self.m_log, debug=self.m_debug)
  463 + result_bgm = like.obs().models().copy()
  464 + LogL_bgm = like.opt().value()
  465 + npred_bgm = like.obs().npred()
  466 +
  467 + # Assess quality based on a comparison between Npred and Nevents
  468 + quality_bgm = npred_bgm-nevents
  469 +
  470 + # Write background fit results
  471 + if self.logExplicit():
  472 + self.log.header3("Background model fit")
  473 + self.log.parformat("log likelihood")
  474 + self.log(LogL_bgm)
  475 + self.log("\n")
  476 + self.log.parformat("Number of predicted events")
  477 + self.log(npred_bgm)
  478 + self.log("\n")
  479 + self.log.parformat("Fit quality")
  480 + self.log(quality_bgm)
  481 + self.log("\n")
  482 +
  483 + # Start over if the fit quality was bad
  484 + if abs(quality_bgm) > 3.0:
  485 + if self.logExplicit():
  486 + self.log("Fit quality outside required range. Start over.\n")
  487 + continue
  488 +
  489 + # Write model fit results
  490 + if self.logExplicit():
  491 + for model in result_bgm:
  492 + self.log.parformat("Model")
  493 + self.log(model.name())
  494 + self.log("\n")
  495 + for par in model:
  496 + self.log(str(par)+"\n")
  497 +
  498 + # Fit background and test source
  499 + sim.models(full_model)
  500 + like = obsutils.fit(sim, log=self.m_log, debug=self.m_debug)
  501 + result_all = like.obs().models().copy()
  502 + LogL_all = like.opt().value()
  503 + npred_all = like.obs().npred()
  504 + ts = 2.0*(LogL_bgm-LogL_all)
  505 +
  506 + # Assess quality based on a comparison between Npred and Nevents
  507 + quality_all = npred_all-nevents
  508 +
  509 + # Write background and test source fit results
  510 + if self.logExplicit():
  511 + self.log.header3("Background and test source model fit")
  512 + self.log.parformat("Test statistics")
  513 + self.log(ts)
  514 + self.log("\n")
  515 + self.log.parformat("log likelihood")
  516 + self.log(LogL_all)
  517 + self.log("\n")
  518 + self.log.parformat("Number of predicted events")
  519 + self.log(npred_all)
  520 + self.log("\n")
  521 + self.log.parformat("Fit quality")
  522 + self.log(quality_all)
  523 + self.log("\n")
  524 + #
  525 + for model in result_all:
  526 + self.log.parformat("Model")
  527 + self.log(model.name())
  528 + self.log("\n")
  529 + for par in model:
  530 + self.log(str(par)+"\n")
  531 +
  532 + # Start over if the fit quality was bad
  533 + if abs(quality_all) > 3.0:
  534 + if self.logExplicit():
  535 + self.log("Fit quality outside required range. Start over.\n")
  536 + continue
  537 +
  538 + # Start over if TS was non-positive
  539 + if ts <= 0.0:
  540 + if self.logExplicit():
  541 + self.log("Non positive TS. Start over.\n")
  542 + continue
  543 +
  544 + # Get fitted Crab, photon and energy fluxes
  545 + cflux = result_all["Test"]['Prefactor'].value()
  546 + cflux_err = result_all["Test"]['Prefactor'].error()
  547 + pflux = result_all["Test"].spectral().flux(emin, emax)
  548 + eflux = result_all["Test"].spectral().eflux(emin, emax)
  549 +
  550 + # Compute differential sensitivity
  551 + diffSens = pflux / erg_width * erg_mean*erg_mean
  552 +
  553 + # Compute flux correction factor based on average TS
  554 + correct = 1.0
  555 + if self.m_use_ts:
  556 + if ts > 0:
  557 + correct = math.sqrt(self.m_ts_thres/ts)
  558 + else:
  559 + if cflux > 0 and cflux_err > 0:
  560 + correct = math.sqrt(self.m_ts_thres)/(cflux/cflux_err)
  561 +
  562 + # Compute extrapolated fluxes
  563 + flux = correct * cflux
  564 + pflux = correct * pflux
  565 + eflux = correct * eflux
  566 + diffSens = correct * diffSens
  567 + flux_value.append(flux)
  568 + pflux_value.append(pflux)
  569 + eflux_value.append(eflux)
  570 + diffSens_value.append(diffSens)
  571 +
  572 + # Write background and test source fit results
  573 + if self.logExplicit():
  574 + self.log.parformat("Photon flux")
  575 + self.log(pflux)
  576 + self.log(" ph/cm2/s\n")
  577 + self.log.parformat("Energy flux")
  578 + self.log(eflux)
  579 + self.log(" erg/cm2/s\n")
  580 + self.log.parformat("Crab flux")
  581 + self.log(cflux*1000.0)
  582 + self.log(" mCrab\n")
  583 + self.log.parformat("Differential sensitivity")
  584 + self.log(diffSens)
  585 + self.log(" erg/cm2/s\n")
  586 + for model in result_all:
  587 + self.log.parformat("Model")
  588 + self.log(model.name())
  589 + self.log("\n")
  590 + for par in model:
  591 + self.log(str(par)+"\n")
  592 + elif self.logTerse():
  593 + self.log.parformat("Iteration "+str(iter))
  594 + self.log("TS=")
  595 + self.log(ts)
  596 + self.log(" ")
  597 + self.log("corr=")
  598 + self.log(correct)
  599 + self.log(" ")
  600 + self.log(pflux)
  601 + self.log(" ph/cm2/s = ")
  602 + self.log(eflux)
  603 + self.log(" erg/cm2/s = ")
  604 + self.log(cflux*1000.0)
  605 + self.log(" mCrab = ")
  606 + self.log(diffSens)
  607 + self.log(" erg/cm2/s\n")
  608 +
  609 + # Compute sliding average of extrapolated fitted prefactor,
  610 + # photon and energy flux. This damps out fluctuations and
  611 + # improves convergence
  612 + flux = 0.0
  613 + pflux = 0.0
  614 + eflux = 0.0
  615 + diffSens = 0.0
  616 + num = 0.0
  617 + for k in range(self.m_num_avg):
  618 + inx = len(flux_value) - k - 1
  619 + if inx >= 0:
  620 + flux += flux_value[inx]
  621 + pflux += pflux_value[inx]
  622 + eflux += eflux_value[inx]
  623 + diffSens += diffSens_value[inx]
  624 + num += 1.0
  625 + flux /= num
  626 + pflux /= num
  627 + eflux /= num
  628 + diffSens /= num
  629 +
  630 + # Compare average flux to last average
  631 + if iter > self.m_num_avg:
  632 + if test_flux > 0:
  633 + ratio = flux/test_flux
  634 +
  635 + # We have 2 convergence criteria:
  636 + # 1. The average flux does not change
  637 + # 2. The flux correction factor is small
  638 + if ratio >= 0.99 and ratio <= 1.01 and \
  639 + correct >= 0.9 and correct <= 1.1:
  640 + if self.logTerse():
  641 + self.log(" Converged ("+str(ratio)+")\n")
  642 + break
  643 + else:
  644 + if self.logTerse():
  645 + self.log(" Flux is zero.\n")
  646 + break
  647 +
  648 + # Use average for next iteration
  649 + test_flux = flux
  650 +
  651 + # Exit loop if number of trials exhausted
  652 + if (iter >= self.m_max_iter):
  653 + if self.logTerse():
  654 + self.log(" Test ended after "+str(self.m_max_iter)+" iterations.\n")
  655 + break
  656 +
  657 + # Write fit results
  658 + if self.logTerse():
  659 + self.log.header3("Fit results")
  660 + self.log.parformat("Test statistics")
  661 + self.log(ts)
  662 + self.log("\n")
  663 + self.log.parformat("Photon flux")
  664 + self.log(pflux)
  665 + self.log(" ph/cm2/s\n")
  666 + self.log.parformat("Energy flux")
  667 + self.log(eflux)
  668 + self.log(" erg/cm2/s\n")
  669 + self.log.parformat("Crab flux")
  670 + self.log(cflux*1000.0)
  671 + self.log(" mCrab\n")
  672 + self.log.parformat("Differential sensitivity")
  673 + self.log(diffSens)
  674 + self.log(" erg/cm2/s\n")
  675 + self.log.parformat("Number of simulated events")
  676 + self.log(nevents)
  677 + self.log("\n")
  678 + self.log.header3("Background and test source model fitting")
  679 + self.log.parformat("log likelihood")
  680 + self.log(LogL_all)
  681 + self.log("\n")
  682 + self.log.parformat("Number of predicted events")
  683 + self.log(npred_all)
  684 + self.log("\n")
  685 + for model in result_all:
  686 + self.log.parformat("Model")
  687 + self.log(model.name())
  688 + self.log("\n")
  689 + for par in model:
  690 + self.log(str(par)+"\n")
  691 + self.log.header3("Background model fit")
  692 + self.log.parformat("log likelihood")
  693 + self.log(LogL_bgm)
  694 + self.log("\n")
  695 + self.log.parformat("Number of predicted events")
  696 + self.log(npred_bgm)
  697 + self.log("\n")
  698 + for model in result_bgm:
  699 + self.log.parformat("Model")
  700 + self.log(model.name())
  701 + self.log("\n")
  702 + for par in model:
  703 + self.log(str(par)+"\n")
  704 +
  705 + # Store result
  706 + result = {'loge': loge, 'emin': emin.TeV(), 'emax': emax.TeV(), \
  707 + 'crab': flux, 'pflux': pflux, 'eflux': eflux, \
  708 + 'diffSens': diffSens}
  709 +
  710 + # Return result
  711 + return result
712 712  
713 713  
714 714 # ======================== #
715 715 # Main routine entry point #
716 716 # ======================== #
717 717 if __name__ == '__main__':
718   - """
719   - CTA ctlike-based sensitivity estimator.
720   - """
721   - # Create instance of application
722   - app = cssens(sys.argv)
723   -
724   - # Run application
725   - app.execute()
726   -
727 718 \ No newline at end of file
  719 + """
  720 + CTA ctlike-based sensitivity estimator.
  721 + """
  722 + # Create instance of application
  723 + app = cssens(sys.argv)
  724 +
  725 + # Run application
  726 + app.execute()
  727 +
728 728 \ No newline at end of file
... ...
scripts/cstsdist.py
  View file @64e71e4
... ... @@ -3,7 +3,7 @@
3 3 # This script generates the TS distribution for a particular model based
4 4 # on Monte-Carlo simulations.
5 5 #
6   -# Copyright (C) 2011-2013 Jurgen Knodlseder
  6 +# Copyright (C) 2011-2014 Juergen Knoedlseder
7 7 #
8 8 # This program is free software: you can redistribute it and/or modify
9 9 # it under the terms of the GNU General Public License as published by
... ... @@ -31,486 +31,486 @@ import math
31 31 # cstsdist class #
32 32 # ============== #
33 33 class cstsdist(GApplication):
34   - """
35   - This class implements the TS distribution generation script. It derives
36   - from the GammaLib::GApplication class which provides support for parameter
37   - files, command line arguments, and logging. In that way the Python
38   - script behaves just as a regular ctool.
39   - """
40   - def __init__(self, *argv):
41   - """
42   - Constructor.
43   - """
44   - # Set name
45   - self.name = "cstsdist"
46   - self.version = "0.1.0"
47   -
48   - # Initialise some members
49   - self.obs = None
50   - self.bkg_model = None
51   - self.full_model = None
52   -
53   - # Make sure that parfile exists
54   - file = self.parfile()
  34 + """
  35 + This class implements the TS distribution generation script. It derives
  36 + from the GammaLib::GApplication class which provides support for parameter
  37 + files, command line arguments, and logging. In that way the Python
  38 + script behaves just as a regular ctool.
  39 + """
  40 + def __init__(self, *argv):
  41 + """
  42 + Constructor.
  43 + """
  44 + # Set name
  45 + self.name = "cstsdist"
  46 + self.version = "0.1.0"
  47 +
  48 + # Initialise some members
  49 + self.obs = None
  50 + self.bkg_model = None
  51 + self.full_model = None
  52 +
  53 + # Make sure that parfile exists
  54 + file = self.parfile()
55 55  
56   - # Initialise application
57   - if len(argv) == 0:
58   - GApplication.__init__(self, self.name, self.version)
59   - elif len(argv) ==1:
60   - GApplication.__init__(self, self.name, self.version, *argv)
61   - else:
62   - raise TypeError("Invalid number of arguments given.")
  56 + # Initialise application
  57 + if len(argv) == 0:
  58 + GApplication.__init__(self, self.name, self.version)
  59 + elif len(argv) ==1:
  60 + GApplication.__init__(self, self.name, self.version, *argv)
  61 + else:
  62 + raise TypeError("Invalid number of arguments given.")
63 63  
64   - # Set logger properties
65   - self.log_header()
66   - self.log.date(True)
  64 + # Set logger properties
  65 + self.log_header()
  66 + self.log.date(True)
67 67  
68   - # Return
69   - return
70   -
71   - def __del__(self):
72   - """
73   - Destructor.
74   - """
75   - # Write separator into logger
76   - if self.logTerse():
77   - self.log("\n")
78   -
79   - # Return
80   - return
  68 + # Return
  69 + return
  70 +
  71 + def __del__(self):
  72 + """
  73 + Destructor.
  74 + """
  75 + # Write separator into logger
  76 + if self.logTerse():
  77 + self.log("\n")
  78 +
  79 + # Return
  80 + return
81 81  
82   - def parfile(self):
83   - """
84   - Check if parfile exists. If parfile does not exist then create a
85   - default parfile. This kluge avoids shipping the cscript with a parfile.
86   - """
87   - # Set parfile name
88   - parfile = self.name+".par"
89   -
90   - try:
91   - pars = GApplicationPars(parfile)
92   - except:
93   - # Signal if parfile was not found
94   - sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
95   -
96   - # Create default parfile
97   - pars = GApplicationPars()
98   - pars.append(GApplicationPar("outfile","f","h","ts.dat","","","Output file name"))
99   - pars.append(GApplicationPar("ntrials","i","a","10","","","Number of trials"))
100   - pars.append(GApplicationPar("caldb","s","a","$GAMMALIB/share/caldb/cta","","","Calibration database"))
101   - pars.append(GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function"))
102   - pars.append(GApplicationPar("type","s","a","point","","","Source model type (point/gauss/shell/disk)"))
103   - pars.append(GApplicationPar("index","r","h","-2.48","","","Spectral index"))
104   - pars.append(GApplicationPar("offset","r","a","0.0","0.0","","Source offset angle (deg)"))
105   - pars.append(GApplicationPar("bkg","s","a","$GAMMALIB/share/models/bkg_dummy.txt","","","Background model file function (none=power law for E)"))
106   - pars.append(GApplicationPar("emin","r","a","0.1","0.0","","Lower energy limit (TeV)"))
107   - pars.append(GApplicationPar("emax","r","a","100.0","0.0","","Upper energy limit (TeV)"))
108   - pars.append(GApplicationPar("enumbins","i","a","0","","","Number of energy bins (0=unbinned)"))
109   - pars.append(GApplicationPar("duration","r","a","180000.0","","","Effective exposure time (s)"))
110   - pars.append(GApplicationPar("rad","r","h","5.0","","","Radius of ROI (deg)"))
111   - pars.append(GApplicationPar("npix","i","h","200","","","Number of pixels for binned"))
112   - pars.append(GApplicationPar("binsz","r","h","0.05","","","Pixel size for binned (deg/pixel)"))
113   - pars.append_standard()
114   - pars.save(parfile)
115   -
116   - # Return
117   - return
118   -
119   - def get_parameters(self):
120   - """
121   - Get parameters from parfile and setup the observation.
122   - """
123   - # Get parameters
124   - self.m_outfile = self["outfile"].filename()
125   - self.m_ntrials = self["ntrials"].integer()
126   - self.m_caldb = self["caldb"].string()
127   - self.m_irf = self["irf"].string()
128   - self.m_type = self["type"].string()
129   - self.m_index = self["index"].real()
130   - self.m_offset = self["offset"].real()
131   - self.m_bkg = self["bkg"].string()
132   - self.m_emin = self["emin"].real()
133   - self.m_emax = self["emax"].real()
134   - self.m_enumbins = self["enumbins"].integer()
135   - self.m_duration = self["duration"].real()
136   - self.m_rad = self["rad"].real()
137   - self.m_npix = self["npix"].integer()
138   - self.m_binsz = self["binsz"].real()
  82 + def parfile(self):
  83 + """
  84 + Check if parfile exists. If parfile does not exist then create a
  85 + default parfile. This kluge avoids shipping the cscript with a parfile.
  86 + """
  87 + # Set parfile name
  88 + parfile = self.name+".par"
  89 +
  90 + try:
  91 + pars = GApplicationPars(parfile)
  92 + except:
  93 + # Signal if parfile was not found
  94 + sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
  95 +
  96 + # Create default parfile
  97 + pars = GApplicationPars()
  98 + pars.append(GApplicationPar("outfile","f","h","ts.dat","","","Output file name"))
  99 + pars.append(GApplicationPar("ntrials","i","a","10","","","Number of trials"))
  100 + pars.append(GApplicationPar("caldb","s","a","$GAMMALIB/share/caldb/cta","","","Calibration database"))
  101 + pars.append(GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function"))
  102 + pars.append(GApplicationPar("type","s","a","point","","","Source model type (point/gauss/shell/disk)"))
  103 + pars.append(GApplicationPar("index","r","h","-2.48","","","Spectral index"))
  104 + pars.append(GApplicationPar("offset","r","a","0.0","0.0","","Source offset angle (deg)"))
  105 + pars.append(GApplicationPar("bkg","s","a","$GAMMALIB/share/models/bkg_dummy.txt","","","Background model file function (none=power law for E)"))
  106 + pars.append(GApplicationPar("emin","r","a","0.1","0.0","","Lower energy limit (TeV)"))
  107 + pars.append(GApplicationPar("emax","r","a","100.0","0.0","","Upper energy limit (TeV)"))
  108 + pars.append(GApplicationPar("enumbins","i","a","0","","","Number of energy bins (0=unbinned)"))
  109 + pars.append(GApplicationPar("duration","r","a","180000.0","","","Effective exposure time (s)"))
  110 + pars.append(GApplicationPar("rad","r","h","5.0","","","Radius of ROI (deg)"))
  111 + pars.append(GApplicationPar("npix","i","h","200","","","Number of pixels for binned"))
  112 + pars.append(GApplicationPar("binsz","r","h","0.05","","","Pixel size for binned (deg/pixel)"))
  113 + pars.append_standard()
  114 + pars.save(parfile)
  115 +
  116 + # Return
  117 + return
  118 +
  119 + def get_parameters(self):
  120 + """
  121 + Get parameters from parfile and setup the observation.
  122 + """
  123 + # Get parameters
  124 + self.m_outfile = self["outfile"].filename()
  125 + self.m_ntrials = self["ntrials"].integer()
  126 + self.m_caldb = self["caldb"].string()
  127 + self.m_irf = self["irf"].string()
  128 + self.m_type = self["type"].string()
  129 + self.m_index = self["index"].real()
  130 + self.m_offset = self["offset"].real()
  131 + self.m_bkg = self["bkg"].string()
  132 + self.m_emin = self["emin"].real()
  133 + self.m_emax = self["emax"].real()
  134 + self.m_enumbins = self["enumbins"].integer()
  135 + self.m_duration = self["duration"].real()
  136 + self.m_rad = self["rad"].real()
  137 + self.m_npix = self["npix"].integer()
  138 + self.m_binsz = self["binsz"].real()
139 139  
140   - # Set some fixed parameters
141   - self.m_log = False # Logging in client tools
142   - self.m_debug = False # Debugging in client tools
143   -
144   - # Setup observation
145   - self.obs = self.set_obs(emin=self.m_emin, emax=self.m_emax)
146   -
147   - # Initialise models. Note that we centre the point source at the Galactic
148   - # center as our observation is also centred at the Galactic centre, so
149   - # we're onaxis.
150   - self.bkg_model = GModels()
151   - self.full_model = GModels()
152   - self.bkg_model.append(self.set_bkg_model())
153   - self.full_model.append(self.set_bkg_model())
154   - self.full_model.append(self.set_src_model(0.0, self.m_offset, \
155   - flux=0.010, \
156   - type=self.m_type, \
157   - index=self.m_index))
  140 + # Set some fixed parameters
  141 + self.m_log = False # Logging in client tools
  142 + self.m_debug = False # Debugging in client tools
  143 +
  144 + # Setup observation
  145 + self.obs = self.set_obs(emin=self.m_emin, emax=self.m_emax)
  146 +
  147 + # Initialise models. Note that we centre the point source at the Galactic
  148 + # center as our observation is also centred at the Galactic centre, so
  149 + # we're onaxis.
  150 + self.bkg_model = GModels()
  151 + self.full_model = GModels()
  152 + self.bkg_model.append(self.set_bkg_model())
  153 + self.full_model.append(self.set_bkg_model())
  154 + self.full_model.append(self.set_src_model(0.0, self.m_offset, \
  155 + flux=0.010, \
  156 + type=self.m_type, \
  157 + index=self.m_index))
158 158  
159   - # Attach background model to observation container
160   - self.obs.models(self.bkg_model)
161   -
162   - # Return
163   - return
164   -
165   - def models(self, models):
166   - """
167   - Set model.
168   - """
169   - # Copy models
170   - self.model = models.copy()
171   -
172   - # Return
173   - return
174   -
175   - def execute(self):
176   - """
177   - Execute the script.
178   - """
179   - # Run the script
180   - self.run()
181   -
182   - # Return
183   - return
  159 + # Attach background model to observation container
  160 + self.obs.models(self.bkg_model)
  161 +
  162 + # Return
  163 + return
  164 +
  165 + def models(self, models):
  166 + """
  167 + Set model.
  168 + """
  169 + # Copy models
  170 + self.model = models.copy()
  171 +
  172 + # Return
  173 + return
  174 +
  175 + def execute(self):
  176 + """
  177 + Execute the script.
  178 + """
  179 + # Run the script
  180 + self.run()
  181 +
  182 + # Return
  183 + return
184 184  
185   - def run(self):
186   - """
187   - Run the script.
188   - """
189   - # Switch screen logging on in debug mode
190   - if self.logDebug():
191   - self.log.cout(True)
  185 + def run(self):
  186 + """
  187 + Run the script.
  188 + """
  189 + # Switch screen logging on in debug mode
  190 + if self.logDebug():
  191 + self.log.cout(True)
192 192  
193   - # Get parameters
194   - self.get_parameters()
195   -
196   - # Write input parameters into logger
197   - if self.logTerse():
198   - self.log_parameters()
199   - self.log("\n")
200   -
201   - # Write observation into logger
202   - if self.logTerse():
203   - self.log("\n")
204   - self.log.header1("Observation")
205   - self.log(str(self.obs))
206   - self.log("\n")
  193 + # Get parameters
  194 + self.get_parameters()
  195 +
  196 + # Write input parameters into logger
  197 + if self.logTerse():
  198 + self.log_parameters()
  199 + self.log("\n")
  200 +
  201 + # Write observation into logger
  202 + if self.logTerse():
  203 + self.log("\n")
  204 + self.log.header1("Observation")
  205 + self.log(str(self.obs))
  206 + self.log("\n")
207 207  
208   - # Write models into logger
209   - if self.logTerse():
210   - self.log("\n")
211   - self.log.header1("Test model")
212   - self.log(str(self.full_model))
213   - self.log("\n")
  208 + # Write models into logger
  209 + if self.logTerse():
  210 + self.log("\n")
  211 + self.log.header1("Test model")
  212 + self.log(str(self.full_model))
  213 + self.log("\n")
214 214  
215   - # Write header
216   - if self.logTerse():
217   - self.log("\n")
218   - self.log.header1("Generate TS distribution")
  215 + # Write header
  216 + if self.logTerse():
  217 + self.log("\n")
  218 + self.log.header1("Generate TS distribution")
219 219  
220   - # Loop over trials
221   - for seed in range(self.m_ntrials):
222   -
223   - # Make a trial
224   - result = self.trial(seed)
225   -
226   - # Write out result immediately
227   - if seed == 0:
228   - file = open(self.m_outfile, 'w')
229   - writer = csv.DictWriter(file, result['colnames'])
230   - writer.writerow(dict((_,_) for _ in result['colnames']))
231   - else:
232   - file = open(self.m_outfile, 'a')
233   - writer = csv.DictWriter(file, result['colnames'])
234   - writer.writerow(result['values'])
235   - file.close()
236   -
237   - # Return
238   - return
239   -
240   - def set_obs(self, lpnt=0.0, bpnt=0.0, emin=0.1, emax=100.0):
241   - """
242   - Returns an observation container with a single CTA observation.
243   -
244   - Keywords:
245   - lpnt - Galactic longitude of pointing [deg] (default: 0.0)
246   - bpnt - Galactic latitude of pointing [deg] (default: 0.0)
247   - emin - Minimum energy [TeV] (default: 0.1)
248   - emax - Maximum energy [TeV] (default: 100.0)
249   - """
250   - # Allocate observation container
251   - obs = GObservations()
252   -
253   - # Set single pointing
254   - pntdir = GSkyDir()
255   - pntdir.lb_deg(lpnt, bpnt)
256   -
257   - # Create CTA observation
258   - run = obsutils.set(pntdir, caldb=self.m_caldb, irf=self.m_irf, \
259   - duration=self.m_duration, \
  220 + # Loop over trials
  221 + for seed in range(self.m_ntrials):
  222 +
  223 + # Make a trial
  224 + result = self.trial(seed)
  225 +
  226 + # Write out result immediately
  227 + if seed == 0:
  228 + file = open(self.m_outfile, 'w')
  229 + writer = csv.DictWriter(file, result['colnames'])
  230 + writer.writerow(dict((_,_) for _ in result['colnames']))
  231 + else:
  232 + file = open(self.m_outfile, 'a')
  233 + writer = csv.DictWriter(file, result['colnames'])
  234 + writer.writerow(result['values'])
  235 + file.close()
  236 +
  237 + # Return
  238 + return
  239 +
  240 + def set_obs(self, lpnt=0.0, bpnt=0.0, emin=0.1, emax=100.0):
  241 + """
  242 + Returns an observation container with a single CTA observation.
  243 +
  244 + Keywords:
  245 + lpnt - Galactic longitude of pointing [deg] (default: 0.0)
  246 + bpnt - Galactic latitude of pointing [deg] (default: 0.0)
  247 + emin - Minimum energy [TeV] (default: 0.1)
  248 + emax - Maximum energy [TeV] (default: 100.0)
  249 + """
  250 + # Allocate observation container
  251 + obs = GObservations()
  252 +
  253 + # Set single pointing
  254 + pntdir = GSkyDir()
  255 + pntdir.lb_deg(lpnt, bpnt)
  256 +
  257 + # Create CTA observation
  258 + run = obsutils.set(pntdir, caldb=self.m_caldb, irf=self.m_irf, \
  259 + duration=self.m_duration, \
260 260 emin=emin, emax=emax, rad=self.m_rad)
261   -
262   - # Append observation to container
263   - obs.append(run)
264   -
265   - # Return observation container
266   - return obs
267   -
268   - def set_bkg_model(self, fitsigma=False):
269   - """
270   - Setup CTA background model.
271   - """
272   - # Define radial component
273   - radial = GCTAModelRadialGauss(3.0)
274   - if fitsigma:
275   - radial["Sigma"].free()
276   - else:
277   - radial["Sigma"].fix()
278   -
279   - # Define spectral component
280   - spectrum = GModelSpectralFunc(self.m_bkg, 1.0)
281   -
282   - # Create background model
283   - model = GCTAModelRadialAcceptance(radial, spectrum)
284   - model.name("Background")
285   - model.instruments("CTA")
286   -
287   - # Return background model
288   - return model
289   -
290   - def set_src_model(self, l, b, flux=1.0, index=-2.48, \
  261 +
  262 + # Append observation to container
  263 + obs.append(run)
  264 +
  265 + # Return observation container
  266 + return obs
  267 +
  268 + def set_bkg_model(self, fitsigma=False):
  269 + """
  270 + Setup CTA background model.
  271 + """
  272 + # Define radial component
  273 + radial = GCTAModelRadialGauss(3.0)
  274 + if fitsigma:
  275 + radial["Sigma"].free()
  276 + else:
  277 + radial["Sigma"].fix()
  278 +
  279 + # Define spectral component
  280 + spectrum = GModelSpectralFunc(self.m_bkg, 1.0)
  281 +
  282 + # Create background model
  283 + model = GCTAModelRadialAcceptance(radial, spectrum)
  284 + model.name("Background")
  285 + model.instruments("CTA")
  286 +
  287 + # Return background model
  288 + return model
  289 +
  290 + def set_src_model(self, l, b, flux=1.0, index=-2.48, \
291 291 type="point", sigma=1.0, radius=1.0, width=0.1, \
292 292 fitpos=False, fitidx=False):
293   - """
294   - Returns a single source with Crab-like spectrum. The source flux
295   - can be scaled in Crab units. The Crab spectrum is based on MAGIC
296   - observations (Albert et al. 2008, ApJ, 674, 1037).
  293 + """
  294 + Returns a single source with Crab-like spectrum. The source flux
  295 + can be scaled in Crab units. The Crab spectrum is based on MAGIC
  296 + observations (Albert et al. 2008, ApJ, 674, 1037).
297 297  
298   - Parameters:
299   - l - Galactic longitude of source location [deg]
300   - b - Galactic latitude of source location [deg]
301   - Keywords:
302   - flux - Source flux [Crabs]
303   - index - Spectral index
304   - type - Source type ("point", "gauss", "disk", "shell")
305   - sigma - Gaussian sigma (for type="gauss")
306   - radius - Disk or shell inner radius [deg] (for type="disk" and type="shell")
307   - width - Shell width [deg] (for type="shell")
308   - fitpos - Fit position? (default: True)
309   - fitidx - Fit index? (default: True)
310   - """
311   - # Set source location
312   - location = GSkyDir()
313   - location.lb_deg(l, b)
314   -
315   - # Set source spectrum
316   - spectrum = GModelSpectralPlaw(flux*5.7e-16, index, GEnergy(0.3, "TeV"))
317   - if fitidx:
318   - spectrum["Index"].free()
319   - else:
320   - spectrum["Index"].fix()
  298 + Parameters:
  299 + l - Galactic longitude of source location [deg]
  300 + b - Galactic latitude of source location [deg]
  301 + Keywords:
  302 + flux - Source flux [Crabs]
  303 + index - Spectral index
  304 + type - Source type ("point", "gauss", "disk", "shell")
  305 + sigma - Gaussian sigma (for type="gauss")
  306 + radius - Disk or shell inner radius [deg] (for type="disk" and type="shell")
  307 + width - Shell width [deg] (for type="shell")
  308 + fitpos - Fit position? (default: True)
  309 + fitidx - Fit index? (default: True)
  310 + """
  311 + # Set source location
  312 + location = GSkyDir()
  313 + location.lb_deg(l, b)
  314 +
  315 + # Set source spectrum
  316 + spectrum = GModelSpectralPlaw(flux*5.7e-16, index, GEnergy(0.3, "TeV"))
  317 + if fitidx:
  318 + spectrum["Index"].free()
  319 + else:
  320 + spectrum["Index"].fix()
321 321  
322   - # Set source
323   - if type == "point":
324   - spatial = GModelSpatialPointSource(location)
325   - if fitpos:
326   - spatial[0].free()
327   - spatial[1].free()
328   - elif type == "gauss":
329   - spatial = GModelSpatialRadialGauss(location, sigma)
330   - if fitpos:
331   - spatial[0].free()
332   - spatial[1].free()
333   - elif type == "disk":
334   - spatial = GModelSpatialRadialDisk(location, radius)
335   - if fitpos:
336   - spatial[0].free()
337   - spatial[1].free()
338   - elif type == "shell":
339   - spatial = GModelSpatialRadialShell(location, radius, width)
340   - if fitpos:
341   - spatial[0].free()
342   - spatial[1].free()
343   - else:
344   - self.log("ERROR: Unknown source type '"+type+"'.\n")
345   - return None
346   - source = GModelSky(spatial, spectrum)
347   -
348   - # Set source name
349   - source.name("Test")
350   -
351   - # Return source
352   - return source
353   -
354   - def trial(self, seed):
355   - """
356   - Create the TS for a single trial.
357   -
358   - Parameters:
359   - seed - Random number generator seed
360   - """
361   - # Write header
362   - if self.logExplicit():
363   - self.log.header2("Trial "+str(seed+1))
  322 + # Set source
  323 + if type == "point":
  324 + spatial = GModelSpatialPointSource(location)
  325 + if fitpos:
  326 + spatial[0].free()
  327 + spatial[1].free()
  328 + elif type == "gauss":
  329 + spatial = GModelSpatialRadialGauss(location, sigma)
  330 + if fitpos:
  331 + spatial[0].free()
  332 + spatial[1].free()
  333 + elif type == "disk":
  334 + spatial = GModelSpatialRadialDisk(location, radius)
  335 + if fitpos:
  336 + spatial[0].free()
  337 + spatial[1].free()
  338 + elif type == "shell":
  339 + spatial = GModelSpatialRadialShell(location, radius, width)
  340 + if fitpos:
  341 + spatial[0].free()
  342 + spatial[1].free()
  343 + else:
  344 + self.log("ERROR: Unknown source type '"+type+"'.\n")
  345 + return None
  346 + source = GModelSky(spatial, spectrum)
  347 +
  348 + # Set source name
  349 + source.name("Test")
  350 +
  351 + # Return source
  352 + return source
  353 +
  354 + def trial(self, seed):
  355 + """
  356 + Create the TS for a single trial.
  357 +
  358 + Parameters:
  359 + seed - Random number generator seed
  360 + """
  361 + # Write header
  362 + if self.logExplicit():
  363 + self.log.header2("Trial "+str(seed+1))
364 364  
365   - # Simulate events
366   - sim = obsutils.sim(self.obs, \
367   - nbins=self.m_enumbins, \
368   - seed=seed, \
369   - binsz=self.m_binsz, \
370   - npix=self.m_npix, \
371   - log=self.m_log, debug=self.m_debug)
  365 + # Simulate events
  366 + sim = obsutils.sim(self.obs, \
  367 + nbins=self.m_enumbins, \
  368 + seed=seed, \
  369 + binsz=self.m_binsz, \
  370 + npix=self.m_npix, \
  371 + log=self.m_log, debug=self.m_debug)
372 372  
373   - # Determine number of events in simulation
374   - nevents = 0.0
375   - for run in sim:
376   - nevents += run.events().number()
  373 + # Determine number of events in simulation
  374 + nevents = 0.0
  375 + for run in sim:
  376 + nevents += run.events().number()
377 377  
378   - # Write simulation results
379   - if self.logExplicit():
380   - self.log.header3("Simulation")
381   - self.log.parformat("Number of simulated events")
382   - self.log(nevents)
383   - self.log("\n")
  378 + # Write simulation results
  379 + if self.logExplicit():
  380 + self.log.header3("Simulation")
  381 + self.log.parformat("Number of simulated events")
  382 + self.log(nevents)
  383 + self.log("\n")
384 384  
385   - # Fit background only
386   - sim.models(self.bkg_model)
387   - like_bgm = obsutils.fit(sim, log=self.m_log, debug=self.m_debug)
388   - result_bgm = like_bgm.obs().models()
389   - LogL_bgm = like_bgm.opt().value()
390   - npred_bgm = like_bgm.obs().npred()
  385 + # Fit background only
  386 + sim.models(self.bkg_model)
  387 + like_bgm = obsutils.fit(sim, log=self.m_log, debug=self.m_debug)
  388 + result_bgm = like_bgm.obs().models()
  389 + LogL_bgm = like_bgm.opt().value()
  390 + npred_bgm = like_bgm.obs().npred()
391 391  
392   - # Write background fit results
393   - if self.logExplicit():
394   - self.log.header3("Background model fit")
395   - self.log.parformat("log likelihood")
396   - self.log(LogL_bgm)
397   - self.log("\n")
398   - self.log.parformat("Number of predicted events")
399   - self.log(npred_bgm)
400   - self.log("\n")
401   - for model in result_bgm:
402   - self.log.parformat("Model")
403   - self.log(model.name())
404   - self.log("\n")
405   - for par in model:
406   - self.log(str(par)+"\n")
  392 + # Write background fit results
  393 + if self.logExplicit():
  394 + self.log.header3("Background model fit")
  395 + self.log.parformat("log likelihood")
  396 + self.log(LogL_bgm)
  397 + self.log("\n")
  398 + self.log.parformat("Number of predicted events")
  399 + self.log(npred_bgm)
  400 + self.log("\n")
  401 + for model in result_bgm:
  402 + self.log.parformat("Model")
  403 + self.log(model.name())
  404 + self.log("\n")
  405 + for par in model:
  406 + self.log(str(par)+"\n")
407 407  
408   - # Fit background and test source
409   - sim.models(self.full_model)
410   - like_all = obsutils.fit(sim, log=self.m_log, debug=self.m_debug)
411   - result_all = like_all.obs().models()
412   - LogL_all = like_all.opt().value()
413   - npred_all = like_all.obs().npred()
414   - ts = 2.0*(LogL_bgm-LogL_all)
415   -
416   - # Write background and test source fit results
417   - if self.logExplicit():
418   - self.log.header3("Background and test source model fit")
419   - self.log.parformat("Test statistics")
420   - self.log(ts)
421   - self.log("\n")
422   - self.log.parformat("log likelihood")
423   - self.log(LogL_all)
424   - self.log("\n")
425   - self.log.parformat("Number of predicted events")
426   - self.log(npred_all)
427   - self.log("\n")
428   - for model in result_all:
429   - self.log.parformat("Model")
430   - self.log(model.name())
431   - self.log("\n")
432   - for par in model:
433   - self.log(str(par)+"\n")
434   -
435   - # Write result
436   - elif self.logTerse():
437   - self.log.parformat("Trial "+str(seed))
438   - self.log("TS=")
439   - self.log(ts)
440   - self.log(" Prefactor=")
441   - self.log(result_all["Test"]["Prefactor"].value())
442   - self.log("+/-")
443   - self.log(result_all["Test"]["Prefactor"].error())
444   - self.log("\n")
445   -
446   - # Initialise results
447   - colnames = []
448   - values = {}
449   -
450   - # Set TS value
451   - colnames.append("TS")
452   - values["TS"] = ts
  408 + # Fit background and test source
  409 + sim.models(self.full_model)
  410 + like_all = obsutils.fit(sim, log=self.m_log, debug=self.m_debug)
  411 + result_all = like_all.obs().models()
  412 + LogL_all = like_all.opt().value()
  413 + npred_all = like_all.obs().npred()
  414 + ts = 2.0*(LogL_bgm-LogL_all)
  415 +
  416 + # Write background and test source fit results
  417 + if self.logExplicit():
  418 + self.log.header3("Background and test source model fit")
  419 + self.log.parformat("Test statistics")
  420 + self.log(ts)
  421 + self.log("\n")
  422 + self.log.parformat("log likelihood")
  423 + self.log(LogL_all)
  424 + self.log("\n")
  425 + self.log.parformat("Number of predicted events")
  426 + self.log(npred_all)
  427 + self.log("\n")
  428 + for model in result_all:
  429 + self.log.parformat("Model")
  430 + self.log(model.name())
  431 + self.log("\n")
  432 + for par in model:
  433 + self.log(str(par)+"\n")
  434 +
  435 + # Write result
  436 + elif self.logTerse():
  437 + self.log.parformat("Trial "+str(seed))
  438 + self.log("TS=")
  439 + self.log(ts)
  440 + self.log(" Prefactor=")
  441 + self.log(result_all["Test"]["Prefactor"].value())
  442 + self.log("+/-")
  443 + self.log(result_all["Test"]["Prefactor"].error())
  444 + self.log("\n")
  445 +
  446 + # Initialise results
  447 + colnames = []
  448 + values = {}
  449 +
  450 + # Set TS value
  451 + colnames.append("TS")
  452 + values["TS"] = ts
453 453  
454   - # Set logL for background fit
455   - colnames.append("LogL_bgm")
456   - values["LogL_bgm"] = LogL_bgm
  454 + # Set logL for background fit
  455 + colnames.append("LogL_bgm")
  456 + values["LogL_bgm"] = LogL_bgm
457 457  
458   - # Set logL for full fit
459   - colnames.append("LogL_all")
460   - values["LogL_all"] = LogL_all
  458 + # Set logL for full fit
  459 + colnames.append("LogL_all")
  460 + values["LogL_all"] = LogL_all
461 461  
462   - # Set Nevents
463   - colnames.append("Nevents")
464   - values["Nevents"] = nevents
  462 + # Set Nevents
  463 + colnames.append("Nevents")
  464 + values["Nevents"] = nevents
465 465  
466   - # Set Npred for background fit
467   - colnames.append("Npred_bkg")
468   - values["Npred_bkg"] = npred_bgm
  466 + # Set Npred for background fit
  467 + colnames.append("Npred_bkg")
  468 + values["Npred_bkg"] = npred_bgm
469 469  
470   - # Set Npred for full fit
471   - colnames.append("Npred_all")
472   - values["Npred_all"] = npred_all
473   -
474   - # Gather free full fit parameters
475   - for i in range(result_all.size()):
476   - model = result_all[i]
477   - model_name = model.name()
478   - for k in range(model.size()):
479   - par = model[k]
480   - if par.is_free():
481   -
482   - # Set parameter name
483   - name = model_name+"_"+par.name()
484   -
485   - # Append value
486   - colnames.append(name)
487   - values[name] = par.value()
488   -
489   - # Append error
490   - name = "Unc_"+name
491   - colnames.append(name)
492   - values[name] = par.error()
493   -
494   - # Bundle together results
495   - result = {'colnames': colnames, 'values': values}
496   -
497   - # Return
498   - return result
  470 + # Set Npred for full fit
  471 + colnames.append("Npred_all")
  472 + values["Npred_all"] = npred_all
  473 +
  474 + # Gather free full fit parameters
  475 + for i in range(result_all.size()):
  476 + model = result_all[i]
  477 + model_name = model.name()
  478 + for k in range(model.size()):
  479 + par = model[k]
  480 + if par.is_free():
  481 +
  482 + # Set parameter name
  483 + name = model_name+"_"+par.name()
  484 +
  485 + # Append value
  486 + colnames.append(name)
  487 + values[name] = par.value()
  488 +
  489 + # Append error
  490 + name = "Unc_"+name
  491 + colnames.append(name)
  492 + values[name] = par.error()
  493 +
  494 + # Bundle together results
  495 + result = {'colnames': colnames, 'values': values}
  496 +
  497 + # Return
  498 + return result
499 499  
500 500  
501 501 # ======================== #
502 502 # Main routine entry point #
503 503 # ======================== #
504 504 if __name__ == '__main__':
505   - """
506   - Generates TS distribution.
507   - """
508   - # Create instance of application
509   - app = cstsdist(sys.argv)
510   -
511   - # Open logfile
512   - app.logFileOpen()
513   -
514   - # Execute application
515   - app.execute()
516   -
517 505 \ No newline at end of file
  506 + """
  507 + Generates TS distribution.
  508 + """
  509 + # Create instance of application
  510 + app = cstsdist(sys.argv)
  511 +
  512 + # Open logfile
  513 + app.logFileOpen()
  514 +
  515 + # Execute application
  516 + app.execute()
  517 +
518 518 \ No newline at end of file
... ...
scripts/cta_kb_bkgrates.py
  View file @64e71e4
... ... @@ -4,7 +4,7 @@
4 4 # into a column separated file function table. The table is automatically
5 5 # named "bkg_XXX.txt" where "XXX" is for example "kb_E_50h_v3".
6 6 #
7   -# Copyright (C) 2011 Jurgen Knodlseder
  7 +# Copyright (C) 2011-2014 Juergen Knoedlseder
8 8 #
9 9 # This program is free software: you can redistribute it and/or modify
10 10 # it under the terms of the GNU General Public License as published by
... ... @@ -29,78 +29,78 @@ import sys
29 29 # Generate file function from model #
30 30 # ================================= #
31 31 def make_file_function(irfname, filename):
32   - """
33   - Generate file function from model.
34   -
35   - Parameters:
36   - irfname - IRF filename
37   - filename - File function filename
38   - """
39   - # Open performance file and file function
40   - irf = open(irfname, "r")
41   - file = open(filename, "w")
42   -
43   - # Read data
44   - for row in irf:
45   -
46   - # Skip header and trailer
47   - if row.find("log(E)") != -1:
48   - continue
49   - if row.find("--------") != -1:
50   - break
51   -
52   - # Extract data
53   - fields = row.replace(" ",",").replace(",,",",").rstrip("\n").split(",")
54   - logE = float(fields[0])
55   - aeff = float(fields[1])
56   - r68 = float(fields[2])
57   - r80 = float(fields[3])
58   - eres = float(fields[4])
59   - bgd = float(fields[5])
60   - sens = float(fields[6])
61   -
62   - # Compute energy (in MeV)
63   - energy = math.pow(10.0, logE)*1.0e6
64   - emin = math.pow(10.0, logE-0.1)*1.0e6
65   - emax = math.pow(10.0, logE+0.1)*1.0e6
66   - ewidth = emax-emin
67   -
68   - # Compute solid angle of r80 region in sr
69   - omega = 2.0 * math.pi * (1.0 - math.cos(math.radians(r80)))
70   -
71   - # Compute background rate per steradian and MeV
72   - bkg_rate = bgd/omega/ewidth
  32 + """
  33 + Generate file function from model.
  34 +
  35 + Parameters:
  36 + irfname - IRF filename
  37 + filename - File function filename
  38 + """
  39 + # Open performance file and file function
  40 + irf = open(irfname, "r")
  41 + file = open(filename, "w")
  42 +
  43 + # Read data
  44 + for row in irf:
  45 +
  46 + # Skip header and trailer
  47 + if row.find("log(E)") != -1:
  48 + continue
  49 + if row.find("--------") != -1:
  50 + break
  51 +
  52 + # Extract data
  53 + fields = row.replace(" ",",").replace(",,",",").rstrip("\n").split(",")
  54 + logE = float(fields[0])
  55 + aeff = float(fields[1])
  56 + r68 = float(fields[2])
  57 + r80 = float(fields[3])
  58 + eres = float(fields[4])
  59 + bgd = float(fields[5])
  60 + sens = float(fields[6])
  61 +
  62 + # Compute energy (in MeV)
  63 + energy = math.pow(10.0, logE)*1.0e6
  64 + emin = math.pow(10.0, logE-0.1)*1.0e6
  65 + emax = math.pow(10.0, logE+0.1)*1.0e6
  66 + ewidth = emax-emin
  67 +
  68 + # Compute solid angle of r80 region in sr
  69 + omega = 2.0 * math.pi * (1.0 - math.cos(math.radians(r80)))
  70 +
  71 + # Compute background rate per steradian and MeV
  72 + bkg_rate = bgd/omega/ewidth
73 73  
74   - # Write
75   - file.write(str(energy)+" "+str(bkg_rate)+"\n")
76   -
77   - # Close files
78   - irf.close()
79   - file.close()
80   -
81   - # Return
82   - return
  74 + # Write
  75 + file.write(str(energy)+" "+str(bkg_rate)+"\n")
  76 +
  77 + # Close files
  78 + irf.close()
  79 + file.close()
  80 +
  81 + # Return
  82 + return
83 83  
84 84  
85 85 # ================ #
86 86 # Main entry point #
87 87 # ================ #
88 88 if __name__ == '__main__':
89   - """
90   - Main entry point
91   - """
92   - # Get IRF from command line argument
93   - if (len(sys.argv) > 1):
94   - irf = sys.argv[1]
95   - else:
96   - irf = "kb_E_50h_v3"
  89 + """
  90 + Main entry point
  91 + """
  92 + # Get IRF from command line argument
  93 + if (len(sys.argv) > 1):
  94 + irf = sys.argv[1]
  95 + else:
  96 + irf = "kb_E_50h_v3"
97 97  
98   - # Set parameters
99   - caldb = "/usr/local/gamma/share/caldb/cta"
100   -
101   - # Set filenames
102   - irfname = caldb+"/"+irf+".dat"
103   - filename = "bkg_"+irf+".txt"
104   -
105   - # Read
106   - make_file_function(irfname, filename)
  98 + # Set parameters
  99 + caldb = "/usr/local/gamma/share/caldb/cta"
  100 +
  101 + # Set filenames
  102 + irfname = caldb+"/"+irf+".dat"
  103 + filename = "bkg_"+irf+".txt"
  104 +
  105 + # Read
  106 + make_file_function(irfname, filename)
... ...
scripts/cta_make_nodes.py
  View file @64e71e4
... ... @@ -4,7 +4,7 @@
4 4 # logarithmically spaced energy bins. The intensity scale is set to the
5 5 # HESS Crab intensity (assuming a power law).
6 6 #
7   -# Copyright (C) 2012 Juergen Knoedlseder
  7 +# Copyright (C) 2012-2014 Juergen Knoedlseder
8 8 #
9 9 # This program is free software: you can redistribute it and/or modify
10 10 # it under the terms of the GNU General Public License as published by
... ... @@ -29,91 +29,91 @@ import math
29 29 # Write one node #
30 30 # ============== #
31 31 def write_node(file, energy, scale):
32   - """
33   - Writes one node to XML file.
34   - """
35   - # Convert to strings
36   - str_energy = str(energy)
37   - str_scale = str(scale)
  32 + """
  33 + Writes one node to XML file.
  34 + """
  35 + # Convert to strings
  36 + str_energy = str(energy)
  37 + str_scale = str(scale)
38 38  
39   - # Write start tag
40   - file.write(' <node>\n')
  39 + # Write start tag
  40 + file.write(' <node>\n')
41 41  
42   - # Write energy
43   - file.write(' <parameter scale="1e6" name="Energy"')
44   - file.write(' min="'+str_energy+'" max="'+str_energy+'"')
45   - file.write(' value="'+str_energy+'"')
46   - file.write(' free="0"/>\n')
  42 + # Write energy
  43 + file.write(' <parameter scale="1e6" name="Energy"')
  44 + file.write(' min="'+str_energy+'" max="'+str_energy+'"')
  45 + file.write(' value="'+str_energy+'"')
  46 + file.write(' free="0"/>\n')
47 47  
48   - # Write intensity
49   - file.write(' <parameter scale="'+str_scale+'" name="Intensity"')
50   - file.write(' min="1e-5" max="1e5"')
51   - file.write(' value="1.0"')
52   - file.write(' free="1"/>\n')
  48 + # Write intensity
  49 + file.write(' <parameter scale="'+str_scale+'" name="Intensity"')
  50 + file.write(' min="1e-5" max="1e5"')
  51 + file.write(' value="1.0"')
  52 + file.write(' free="1"/>\n')
53 53  
54   - # Write end tag
55   - file.write(' </node>\n')
  54 + # Write end tag
  55 + file.write(' </node>\n')
56 56  
57   - # Return
58   - return
59   -
  57 + # Return
  58 + return
  59 +
60 60  
61 61 # ============ #
62 62 # Create nodes #
63 63 # ============ #
64 64 def create_nodes(emin, emax, enumbins):
65   - """
66   - Create nodes (energies in TeV).
67   - """
68   - # Open file
69   - file = open("nodes.xml", "w")
  65 + """
  66 + Create nodes (energies in TeV).
  67 + """
  68 + # Open file
  69 + file = open("nodes.xml", "w")
70 70  
71   - # Set node boundaries
72   - elogmin = math.log10(float(emin))
73   - elogmax = math.log10(float(emax))
74   - elogbin = (elogmax - elogmin)/(float(enumbins)-1.0)
75   -
76   - # Fill arrays
77   - for i in range(int(enumbins)):
78   -
79   - # Compute energy
80   - energy = math.pow(10.0, i*elogbin+elogmin)
81   -
82   - # Compute scale (HESS Crab spectrum)
83   - scale = 3.45e-17 * math.pow(energy, -2.63)
84   -
85   - # Write node
86   - write_node(file, energy, scale)
87   -
88   - # Debug
89   - #sys.stdout.write(energy)
  71 + # Set node boundaries
  72 + elogmin = math.log10(float(emin))
  73 + elogmax = math.log10(float(emax))
  74 + elogbin = (elogmax - elogmin)/(float(enumbins)-1.0)
  75 +
  76 + # Fill arrays
  77 + for i in range(int(enumbins)):
  78 +
  79 + # Compute energy
  80 + energy = math.pow(10.0, i*elogbin+elogmin)
  81 +
  82 + # Compute scale (HESS Crab spectrum)
  83 + scale = 3.45e-17 * math.pow(energy, -2.63)
  84 +
  85 + # Write node
  86 + write_node(file, energy, scale)
  87 +
  88 + # Debug
  89 + #sys.stdout.write(energy)
90 90 #sys.stdout.write(math.pow(10.0, int(math.log10(scale))-1.0)±"\n")
91   -
92   - # Close file
93   - file.close()
94   -
95   - # Return
96   - return
97   -
  91 +
  92 + # Close file
  93 + file.close()
  94 +
  95 + # Return
  96 + return
  97 +
98 98  
99 99 # ======================= #
100 100 # Main script entry point #
101 101 # ======================= #
102 102 if __name__ == '__main__':
103   - """
104   - Compute nodes.
105   - """
106   - # Check command line
107   - usage = "Usage: cta_make_nodes emin emax enumbins"
108   - if len(sys.argv) < 3:
109   - sys.stdout.write(usage+"\n")
110   - sys.exit()
111   -
112   - # Extract parameters
113   - emin = sys.argv[1]
114   - emax = sys.argv[2]
115   - enumbins = sys.argv[3]
  103 + """
  104 + Compute nodes.
  105 + """
  106 + # Check command line
  107 + usage = "Usage: cta_make_nodes emin emax enumbins"
  108 + if len(sys.argv) < 3:
  109 + sys.stdout.write(usage+"\n")
  110 + sys.exit()
  111 +
  112 + # Extract parameters
  113 + emin = sys.argv[1]
  114 + emax = sys.argv[2]
  115 + enumbins = sys.argv[3]
116 116  
117   - # Create nodes
118   - create_nodes(emin, emax, enumbins)
119   -
  117 + # Create nodes
  118 + create_nodes(emin, emax, enumbins)
  119 +
... ...
scripts/cta_root2irf.py
  View file @64e71e4
... ... @@ -11,7 +11,7 @@
11 11 # area (i.e. the area before applying any theta cut). This is required by
12 12 # GammaLib as the maximum likelihood analysis does not imply any theta cut.
13 13 #
14   -# Copyright (C) 2011 Jurgen Knodlseder
  14 +# Copyright (C) 2011-2014 Juergen Knoedlseder
15 15 #
16 16 # This program is free software: you can redistribute it and/or modify
17 17 # it under the terms of the GNU General Public License as published by
... ...