Johan Bregeon / ctools

Commit 0a5b1abee12549fdc2271f704d7df471e1cc4d34

Authored by Jürgen Knödlseder
1 parent 63fe865c

Add support modules to cscripts (#1830) 

To avoid code duplication in cscripts the support modules "modutils' and "ioutils" have been added. The obsolete functions cntmap() and modmap() have been removed from the "obsutils" module.
Inline Side-by-side
Showing 11 changed files with 315 additions and 500 deletions
ChangeLog
  View file @0a5b1ab
1   -2016-07-29
  1 +2016-07-31
2 2  
3 3 * Version 1.2.0 released
4 4 ========================
5 5  
  6 + Add support modules to cscripts (#1830)
6 7 Fix lacking coverage measure for cscripts execution (#1827)
7 8 Set "addbounds" parameters by default to "no" (#1819)
8 9  
... ...
1 1 New Features and Important Changes in ctools 1.2.0
2 2  
3   -29 July 2016
  3 +31 July 2016
4 4  
5 5  
6 6 Introduction
... ... @@ -105,9 +105,11 @@ cterror - Likelihood profile error computation
105 105 None
106 106  
107 107  
108   -scripts
109   --------
110   -None
  108 +cscripts
  109 +--------
  110 +Added the support modules "modutils" and "ioutils" to the cscripts module
  111 +and removed the obsolete function "cntmap()" and "modmap()" from the
  112 +"obsutils" module (#1830).
111 113  
112 114  
113 115 Examples
... ...
README.md
  View file @0a5b1ab
1 1 ctools information
2 2 ==================
3   -* Version: 1.2.0.dev1 (29 July 2016)
  3 +* Version: 1.2.0.dev1 (31 July 2016)
4 4 * Author: Juergen Knoedlseder (jurgen.knodlseder@irap.omp.eu)
5 5 * GammaLib dependency: 1.2.0.dev1
6 6  
... ...
cscripts/Makefile.am
  View file @0a5b1ab
... ... @@ -42,7 +42,9 @@ pythonscripts = $(builddir)/__init__.py \
42 42 $(srcdir)/csiactcopy.py \
43 43 $(srcdir)/csiactdata.py \
44 44 $(srcdir)/csiactobs.py \
45   - $(srcdir)/obsutils.py
  45 + $(srcdir)/obsutils.py \
  46 + $(srcdir)/modutils.py \
  47 + $(srcdir)/ioutils.py
46 48  
47 49 # Parfiles to be distributed in $(prefix)/syspfiles
48 50 parfiles = $(srcdir)/cscaldb.par \
... ... @@ -203,6 +205,8 @@ clean-local:
203 205 @if test '@srcdir@' != '@builddir@'; then \
204 206 rm -rf $(builddir)/cs*.py; \
205 207 rm -rf $(builddir)/obsutils.py; \
  208 + rm -rf $(builddir)/modutils.py; \
  209 + rm -rf $(builddir)/ioutils.py; \
206 210 rm -rf $(builddir)/*.par; \
207 211 fi
208 212 rm -rf $(builddir)/*.pyc
... ...
cscripts/__init__.py.in
  View file @0a5b1ab
... ... @@ -41,7 +41,9 @@ __all__ = [
41 41 "cstsmapmerge",
42 42 "cstsmapsplit",
43 43 "csworkflow",
44   - "obsutils"
  44 + "obsutils",
  45 + "modutils",
  46 + "ioutils"
45 47 ]
46 48 if cur_version > req_version:
47 49 __all__.extend(["csfindobs",
... ... @@ -68,6 +70,8 @@ from cscripts.cstsmapmerge import cstsmapmerge
68 70 from cscripts.cstsmapsplit import cstsmapsplit
69 71 from cscripts.csworkflow import csworkflow
70 72 from cscripts import obsutils
  73 +from cscripts import modutils
  74 +from cscripts import ioutils
71 75 if cur_version > req_version:
72 76 from cscripts.csfindobs import csfindobs
73 77 from cscripts.csiactcopy import csiactcopy
... ...
cscripts/cspull.py
  View file @0a5b1ab
... ... @@ -19,11 +19,11 @@
19 19 #
20 20 # ==========================================================================
21 21 import sys
22   -import csv
23 22 import math
24 23 import gammalib
25 24 import ctools
26 25 from cscripts import obsutils
  26 +from cscripts import ioutils
27 27  
28 28  
29 29 # ============ #
... ... @@ -79,7 +79,18 @@ class cspull(ctools.cscript):
79 79 # observation parameters to set wobble pattern
80 80 self._pattern = self['pattern'].string()
81 81 if self._pattern == 'four':
82   - self._obs = self._set_obs()
  82 + self._obs = obsutils.set_observations(self['ra'].real(),
  83 + self['dec'].real(),
  84 + self['rad'].real(),
  85 + self['tmin'].real(),
  86 + self['tmax'].real(),
  87 + self['emin'].real(),
  88 + self['emax'].real(),
  89 + self['irf'].string(),
  90 + self['caldb'].string(),
  91 + deadc=self['deadc'].real(),
  92 + pattern=self['pattern'].string(),
  93 + offset=self['offset'].real())
83 94  
84 95 # ... otherwise add response information and energy boundaries
85 96 # in case they are missing.
... ... @@ -119,38 +130,6 @@ class cspull(ctools.cscript):
119 130 # Return
120 131 return
121 132  
122   - def _set_obs(self):
123   - """
124   - Set CTA observation(s)
125   -
126   - Returns
127   - -------
128   - obs : `~gammalib.GObservations`
129   - Observation container.
130   - """
131   - # Set duration
132   - duration = self['tmax'].real() - self['tmin'].real()
133   -
134   - # Setup observation definition list
135   - obsdeflist = obsutils.set_obs_patterns(self._pattern,
136   - ra=self['ra'].real(),
137   - dec=self['dec'].real(),
138   - offset=self['offset'].real())
139   -
140   - # Create list of observations
141   - obs = obsutils.set_obs_list(obsdeflist,
142   - tstart=self['tmin'].real(),
143   - duration=duration,
144   - deadc=self['deadc'].real(),
145   - emin=self['emin'].real(),
146   - emax=self['emax'].real(),
147   - rad=self['rad'].real(),
148   - irf=self['irf'].string(),
149   - caldb=self['caldb'].string())
150   -
151   - # Return observation container
152   - return obs
153   -
154 133 def _obs_string(self, obs):
155 134 """
156 135 Generate summary string for observation
... ... @@ -359,22 +338,9 @@ class cspull(ctools.cscript):
359 338 # Make a trial and add initial seed
360 339 result = self._trial(seed + self._seed)
361 340  
362   - # Set output filename
363   - outfile = self['outfile'].filename().url()
364   -
365   - # Write out result immediately
366   - if seed == 0:
367   - f = open(outfile, 'w')
368   - writer = csv.DictWriter(f, result['colnames'])
369   - headers = {}
370   - for n in result['colnames']:
371   - headers[n] = n
372   - writer.writerow(headers)
373   - else:
374   - f = open(outfile, 'a')
375   - writer = csv.DictWriter(f, result['colnames'])
376   - writer.writerow(result['values'])
377   - f.close()
  341 + # Write out trial result
  342 + ioutils.write_csv_row(self['outfile'].filename().url(), seed,
  343 + result['colnames'], result['values'])
378 344  
379 345 # Return
380 346 return
... ...
cscripts/cssens.py
  View file @0a5b1ab
1 1 #! /usr/bin/env python
2 2 # ==========================================================================
3   -# This script computes the CTA sensitivity.
  3 +# Computes the array sensitivity using a test source
4 4 #
5 5 # Copyright (C) 2011-2016 Juergen Knoedlseder
6 6 #
... ... @@ -24,6 +24,8 @@ import math
24 24 import gammalib
25 25 import ctools
26 26 from cscripts import obsutils
  27 +from cscripts import modutils
  28 +from cscripts import ioutils
27 29  
28 30  
29 31 # ============ #
... ... @@ -190,70 +192,6 @@ class cssens(ctools.cscript):
190 192 # Return observation container
191 193 return obs
192 194  
193   - def _set_models(self, fitspat=False, fitspec=False):
194   - """
195   - Set full model and background model
196   -
197   - Parameters
198   - ----------
199   - fitspat : bool, optional
200   - Fit spatial parameter?
201   - fitspec : bool, optional
202   - Fit spectral parameters?
203   -
204   - Returns
205   - -------
206   - full_model, bkg_model : `~gammalib.GModels`
207   - Full model and background model containers
208   - """
209   - # Retrieve full model from observation container
210   - full_model = self._obs.models().copy()
211   -
212   - # Get source model
213   - model = full_model[self._srcname]
214   -
215   - # If source model has no "Prefactor" parameter then raise an exception
216   - if not model.has_par('Prefactor'):
217   - msg = ('Model "%s" has no parameter "Prefactor". Only spectral '
218   - 'models with a "Prefactor" parameter are supported.' %
219   - self._srcname)
220   - raise RuntimeError(msg)
221   -
222   - # Set source position
223   - if self._ra != None and self._dec != None:
224   - if model.has_par('RA') and model.has_par('DEC'):
225   - model['RA'].value(self._ra)
226   - model['DEC'].value(self._dec)
227   -
228   - # Set possible spatial and spectral parameters
229   - spatial = ['RA', 'DEC', 'Sigma', 'Radius', 'Width', 'PA',
230   - 'MinorRadius', 'MajorRadius']
231   - spectral = ['Index', 'Index1', 'Index2', 'BreakEnergy', 'CutoffEnergy',
232   - 'InverseCutoffEnergy']
233   -
234   - # Fit or fix spatial parameters
235   - for par in spatial:
236   - if model.has_par(par):
237   - if fitspat:
238   - model[par].free()
239   - else:
240   - model[par].fix()
241   -
242   - # Fit or fix spectral parameters
243   - for par in spectral:
244   - if model.has_par(par):
245   - if fitspec:
246   - model[par].free()
247   - else:
248   - model[par].fix()
249   -
250   - # Create background model
251   - bkg_model = full_model.copy()
252   - bkg_model.remove(self._srcname)
253   -
254   - # Return models
255   - return full_model, bkg_model
256   -
257 195 def _set_obs_ebounds(self, emin, emax):
258 196 """
259 197 Set energy boundaries for all observations in container
... ... @@ -330,9 +268,9 @@ class cssens(ctools.cscript):
330 268 # Return photon flux
331 269 return flux
332 270  
333   - def _get_sensitivity(self, emin, emax, bkg_model, full_model):
  271 + def _get_sensitivity(self, emin, emax, test_model):
334 272 """
335   - Determine sensitivity for given observations.
  273 + Determine sensitivity for given observations
336 274  
337 275 Parameters
338 276 ----------
... ... @@ -340,10 +278,8 @@ class cssens(ctools.cscript):
340 278 Minimum energy for fitting and flux computation
341 279 emax : `~gammalib.GEnergy`
342 280 Maximum energy for fitting and flux computation
343   - bkg_model : `~gammalib.GModels`
344   - Background model
345   - full_model : `~gammalib.GModels`
346   - Source plus background model
  281 + test_model : `~gammalib.GModels`
  282 + Test source model
347 283  
348 284 Returns
349 285 -------
... ... @@ -364,7 +300,7 @@ class cssens(ctools.cscript):
364 300 # Compute Crab unit (this is the factor with which the Prefactor needs
365 301 # to be multiplied to get 1 Crab
366 302 crab_flux = self._get_crab_flux(emin, emax)
367   - src_flux = full_model[self._srcname].spectral().flux(emin, emax)
  303 + src_flux = test_model[self._srcname].spectral().flux(emin, emax)
368 304 crab_unit = crab_flux/src_flux
369 305  
370 306 # Write header
... ... @@ -399,14 +335,14 @@ class cssens(ctools.cscript):
399 335 if self._logExplicit():
400 336 self._log.header2('Iteration '+str(iterations))
401 337  
402   - # Set source model. crab_prefactor is the Prefactor that
  338 + # Set test source model. crab_prefactor is the Prefactor that
403 339 # corresponds to 1 Crab
404   - src_model = full_model.copy()
405   - crab_prefactor = src_model[self._srcname]['Prefactor'].value() * \
  340 + models = test_model.copy()
  341 + crab_prefactor = models[self._srcname]['Prefactor'].value() * \
406 342 crab_unit
407   - src_model[self._srcname]['Prefactor'].value(crab_prefactor * \
408   - test_crab_flux)
409   - self._obs.models(src_model)
  343 + models[self._srcname]['Prefactor'].value(crab_prefactor * \
  344 + test_crab_flux)
  345 + self._obs.models(models)
410 346  
411 347 # Simulate events
412 348 sim = obsutils.sim(self._obs, nbins=self._enumbins, seed=iterations,
... ... @@ -426,68 +362,35 @@ class cssens(ctools.cscript):
426 362 self._log(nevents)
427 363 self._log('\n')
428 364  
429   - # Fit background only
430   - sim.models(bkg_model)
431   - like = obsutils.fit(sim, log=self._log_clients,
432   - debug=self._debug, edisp=self._edisp)
433   - result_bgm = like.obs().models().copy()
434   - LogL_bgm = like.opt().value()
435   - npred_bgm = like.obs().npred()
436   -
437   - # Assess quality based on a comparison between Npred and Nevents
438   - quality_bgm = npred_bgm-nevents
439   -
440   - # Write background fit results
441   - if self._logExplicit():
442   - self._log.header3('Background model fit')
443   - self._log.parformat('log likelihood')
444   - self._log(LogL_bgm)
445   - self._log('\n')
446   - self._log.parformat('Number of predicted events')
447   - self._log(npred_bgm)
448   - self._log('\n')
449   - self._log.parformat('Fit quality')
450   - self._log(quality_bgm)
451   - self._log('\n')
452   -
453   - # Write model fit results
454   - if self._logExplicit():
455   - for model in result_bgm:
456   - self._log.parformat('Model')
457   - self._log(model.name())
458   - self._log('\n')
459   - for par in model:
460   - self._log(str(par)+'\n')
  365 + # Fit test source
  366 + fit = obsutils.fit(sim, log=self._log_clients, debug=self._debug,
  367 + edisp=self._edisp)
461 368  
462   - # Fit background and test source
463   - sim.models(src_model)
464   - like = obsutils.fit(sim, log=self._log_clients,
465   - debug=self._debug, edisp=self._edisp)
466   - result_all = like.obs().models().copy()
467   - LogL_all = like.opt().value()
468   - npred_all = like.obs().npred()
469   - ts = 2.0*(LogL_bgm-LogL_all)
  369 + # Get model fitting results
  370 + logL = fit.opt().value()
  371 + npred = fit.obs().npred()
  372 + models = fit.obs().models()
  373 + ts = models[self._srcname].ts()
470 374  
471 375 # Assess quality based on a comparison between Npred and Nevents
472   - quality_all = npred_all-nevents
  376 + quality = npred - nevents
473 377  
474   - # Write background and test source fit results
  378 + # Write test source fit results
475 379 if self._logExplicit():
476   - self._log.header3('Background and test source model fit')
  380 + self._log.header3('Test source model fit')
477 381 self._log.parformat('Test statistics')
478 382 self._log(ts)
479 383 self._log('\n')
480 384 self._log.parformat('log likelihood')
481   - self._log(LogL_all)
  385 + self._log(logL)
482 386 self._log('\n')
483 387 self._log.parformat('Number of predicted events')
484   - self._log(npred_all)
  388 + self._log(npred)
485 389 self._log('\n')
486 390 self._log.parformat('Fit quality')
487   - self._log(quality_all)
  391 + self._log(quality)
488 392 self._log('\n')
489   - #
490   - for model in result_all:
  393 + for model in models:
491 394 self._log.parformat('Model')
492 395 self._log(model.name())
493 396 self._log('\n')
... ... @@ -501,17 +404,14 @@ class cssens(ctools.cscript):
501 404 continue
502 405  
503 406 # Get fitted Crab, photon and energy fluxes
504   - crab_flux = result_all[self._srcname]['Prefactor'].value() / \
505   - crab_prefactor
506   - #crab_flux_err = result_all[self._srcname]['Prefactor'].error() / \
507   - # crab_prefactor
508   - photon_flux = result_all[self._srcname].spectral().flux(emin, emax)
509   - energy_flux = result_all[self._srcname].spectral().eflux(emin, emax)
  407 + crab_flux = models[self._srcname]['Prefactor'].value() / \
  408 + crab_prefactor
  409 + photon_flux = models[self._srcname].spectral().flux(emin, emax)
  410 + energy_flux = models[self._srcname].spectral().eflux(emin, emax)
510 411  
511 412 # Compute differential sensitivity in unit erg/cm2/s
512 413 energy = gammalib.GEnergy(e_mean, 'TeV')
513   - time = gammalib.GTime()
514   - sensitivity = result_all[self._srcname].spectral().eval(energy, time) * \
  414 + sensitivity = models[self._srcname].spectral().eval(energy) * \
515 415 e_mean*erg_mean * 1.0e6
516 416  
517 417 # Compute flux correction factor based on average TS
... ... @@ -543,7 +443,7 @@ class cssens(ctools.cscript):
543 443 self._log.parformat('Differential sensitivity')
544 444 self._log(sensitivity)
545 445 self._log(' erg/cm2/s\n')
546   - for model in result_all:
  446 + for model in models:
547 447 self._log.parformat('Model')
548 448 self._log(model.name())
549 449 self._log('\n')
... ... @@ -636,27 +536,14 @@ class cssens(ctools.cscript):
636 536 self._log.parformat('Number of simulated events')
637 537 self._log(nevents)
638 538 self._log('\n')
639   - self._log.header3('Background and test source model fitting')
  539 + self._log.header3('Test source model fitting')
640 540 self._log.parformat('log likelihood')
641   - self._log(LogL_all)
  541 + self._log(logL)
642 542 self._log('\n')
643 543 self._log.parformat('Number of predicted events')
644   - self._log(npred_all)
  544 + self._log(npred)
645 545 self._log('\n')
646   - for model in result_all:
647   - self._log.parformat('Model')
648   - self._log(model.name())
649   - self._log('\n')
650   - for par in model:
651   - self._log(str(par)+'\n')
652   - self._log.header3('Background model fit')
653   - self._log.parformat('log likelihood')
654   - self._log(LogL_bgm)
655   - self._log('\n')
656   - self._log.parformat('Number of predicted events')
657   - self._log(npred_bgm)
658   - self._log('\n')
659   - for model in result_bgm:
  546 + for model in models:
660 547 self._log.parformat('Model')
661 548 self._log(model.name())
662 549 self._log('\n')
... ... @@ -698,18 +585,14 @@ class cssens(ctools.cscript):
698 585 'energy_flux', 'sensitivity']
699 586 results = []
700 587  
701   - # Initialise models
702   - full_model, bkg_model = self._set_models()
  588 + # Set test source model for this observation
  589 + models = modutils.test_source(self._obs.models(), self._srcname)
703 590  
704 591 # Write models into logger
705 592 if self._logTerse():
706 593 self._log('\n')
707 594 self._log.header1('Models')
708   - self._log.header2('Background model')
709   - self._log(str(bkg_model))
710   - self._log('\n\n')
711   - self._log.header2('Full model')
712   - self._log(str(full_model))
  595 + self._log(str(models))
713 596 self._log('\n')
714 597  
715 598 # Write header
... ... @@ -736,23 +619,10 @@ class cssens(ctools.cscript):
736 619 raise RuntimeError(msg)
737 620  
738 621 # Determine sensitivity
739   - result = self._get_sensitivity(emin, emax, bkg_model, full_model)
740   -
741   - # Write results
742   - if ieng == 0:
743   - f = open(self._outfile.url(), 'w')
744   - writer = csv.DictWriter(f, colnames)
745   - headers = {}
746   - for n in colnames:
747   - headers[n] = n
748   - writer.writerow(headers)
749   - writer.writerow(result)
750   - f.close()
751   - else:
752   - f = open(self._outfile.url(), 'a')
753   - writer = csv.DictWriter(f, colnames)
754   - writer.writerow(result)
755   - f.close()
  622 + result = self._get_sensitivity(emin, emax, models)
  623 +
  624 + # Write out trial result
  625 + ioutils.write_csv_row(self._outfile.url(), ieng, colnames, result)
756 626  
757 627 # Append results
758 628 results.append(result)
... ...
cscripts/cstsdist.py
  View file @0a5b1ab
... ... @@ -19,10 +19,11 @@
19 19 #
20 20 # ==========================================================================
21 21 import sys
22   -import csv
23 22 import gammalib
24 23 import ctools
25 24 from cscripts import obsutils
  25 +from cscripts import modutils
  26 +from cscripts import ioutils
26 27  
27 28  
28 29 # ============== #
... ... @@ -68,7 +69,18 @@ class cstsdist(ctools.cscript):
68 69 # Check for requested pattern and use above observation parameters
69 70 # to set wobble pattern
70 71 if self['pattern'].string() == 'four':
71   - self._obs = self._set_obs()
  72 + self._obs = obsutils.set_observations(self['ra'].real(),
  73 + self['dec'].real(),
  74 + self['rad'].real(),
  75 + self['tmin'].real(),
  76 + self['tmax'].real(),
  77 + self['emin'].real(),
  78 + self['emax'].real(),
  79 + self['irf'].string(),
  80 + self['caldb'].string(),
  81 + deadc=self['deadc'].real(),
  82 + pattern=self['pattern'].string(),
  83 + offset=self['offset'].real())
72 84  
73 85 # Get source name
74 86 self._srcname = self['srcname'].string()
... ... @@ -98,38 +110,6 @@ class cstsdist(ctools.cscript):
98 110 # Return
99 111 return
100 112  
101   - def _set_obs(self):
102   - """
103   - Set an observation container
104   -
105   - Returns
106   - -------
107   - obs : `~gammalib.GObservations`
108   - Observation container
109   - """
110   - # Set duration
111   - duration = self['tmax'].real() - self['tmin'].real()
112   -
113   - # Setup observation definition list
114   - obsdeflist = obsutils.set_obs_patterns(self['pattern'].string(),
115   - ra = self['ra'].real(),
116   - dec = self['dec'].real(),
117   - offset = self['offset'].real())
118   -
119   - # Create list of observations
120   - obs = obsutils.set_obs_list(obsdeflist,
121   - tstart = self['tmin'].real(),
122   - duration = duration,
123   - deadc = self['deadc'].real(),
124   - emin = self['emin'].real(),
125   - emax = self['emax'].real(),
126   - rad = self['rad'].real(),
127   - irf = self['irf'].string(),
128   - caldb = self['caldb'].string())
129   -
130   - # Return observation container
131   - return obs
132   -
133 113 def _trial(self, seed):
134 114 """
135 115 Create the TS for a single trial
... ... @@ -276,8 +256,8 @@ class cstsdist(ctools.cscript):
276 256 # Get parameters
277 257 self._get_parameters()
278 258  
279   - # Set models for this observation
280   - self.models(obsutils.set_ts_model(self._obs.models(), self._srcname))
  259 + # Set test source model for this observation
  260 + self.models(modutils.test_source(self._obs.models(), self._srcname))
281 261  
282 262 # Write models into logger
283 263 if self._logTerse():
... ... @@ -308,22 +288,9 @@ class cstsdist(ctools.cscript):
308 288 # Make a trial
309 289 result = self._trial(seed)
310 290  
311   - # Set output filename
312   - outfile = self['outfile'].filename().url()
313   -
314   - # Write out result immediately
315   - if seed == 0:
316   - f = open(outfile, 'w')
317   - writer = csv.DictWriter(f, result['colnames'])
318   - headers = {}
319   - for n in result['colnames']:
320   - headers[n] = n
321   - writer.writerow(headers)
322   - else:
323   - f = open(outfile, 'a')
324   - writer = csv.DictWriter(f, result['colnames'])
325   - writer.writerow(result['values'])
326   - f.close()
  291 + # Write out trial result
  292 + ioutils.write_csv_row(self['outfile'].filename().url(), seed,
  293 + result['colnames'], result['values'])
327 294  
328 295 # Return
329 296 return
... ...
cscripts/ioutils.py 0 → 100755
  View file @0a5b1ab
  1 +# ==========================================================================
  2 +# Utility functions for input and output
  3 +#
  4 +# Copyright (C) 2016 Juergen Knoedlseder
  5 +#
  6 +# This program is free software: you can redistribute it and/or modify
  7 +# it under the terms of the GNU General Public License as published by
  8 +# the Free Software Foundation, either version 3 of the License, or
  9 +# (at your option) any later version.
  10 +#
  11 +# This program is distributed in the hope that it will be useful,
  12 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  14 +# GNU General Public License for more details.
  15 +#
  16 +# You should have received a copy of the GNU General Public License
  17 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  18 +#
  19 +# ==========================================================================
  20 +import csv
  21 +import gammalib
  22 +import ctools
  23 +
  24 +
  25 +# ========================= #
  26 +# Write one row of CSV file #
  27 +# ========================= #
  28 +def write_csv_row(outfile, row, colnames, colvalues):
  29 + """
  30 + Write one row into CSV file
  31 +
  32 + Parameters
  33 + ----------
  34 + outfile : str
  35 + Output file name
  36 + row : int
  37 + Row number
  38 + colnames : list of str
  39 + Column names
  40 + colvalues : list of float
  41 + Column values
  42 + """
  43 + # If this is the first row then open a new file and write the header
  44 + if row == 0:
  45 + f = open(outfile, 'w')
  46 + writer = csv.DictWriter(f, colnames)
  47 + headers = {}
  48 + for colname in colnames:
  49 + headers[colname] = colname
  50 + writer.writerow(headers)
  51 +
  52 + # ... otherwise append to an existing file
  53 + else:
  54 + f = open(outfile, 'a')
  55 +
  56 + # Write out row
  57 + writer = csv.DictWriter(f, colnames)
  58 + writer.writerow(colvalues)
  59 +
  60 + # Close file
  61 + f.close()
  62 +
  63 + # Return
  64 + return
... ...
cscripts/modutils.py 0 → 100755
  View file @0a5b1ab
  1 +# ==========================================================================
  2 +# Utility functions for model handling
  3 +#
  4 +# Copyright (C) 2016 Juergen Knoedlseder
  5 +#
  6 +# This program is free software: you can redistribute it and/or modify
  7 +# it under the terms of the GNU General Public License as published by
  8 +# the Free Software Foundation, either version 3 of the License, or
  9 +# (at your option) any later version.
  10 +#
  11 +# This program is distributed in the hope that it will be useful,
  12 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  14 +# GNU General Public License for more details.
  15 +#
  16 +# You should have received a copy of the GNU General Public License
  17 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  18 +#
  19 +# ==========================================================================
  20 +import gammalib
  21 +import ctools
  22 +
  23 +
  24 +# ============================================ #
  25 +# Return model for TS fitting of a test source #
  26 +# ============================================ #
  27 +def test_source(models, srcname, ra=None, dec=None, fitspat=False,
  28 + fitspec=False):
  29 + """
  30 + Return model for TS fitting of a test source
  31 +
  32 + Parameters
  33 + ----------
  34 + models : `~gammalib.GModels`
  35 + Input model container
  36 + srcname : str
  37 + Test source name
  38 + ra : float, optional
  39 + Right Ascension of test source (deg)
  40 + dec : float, optional
  41 + Declination of test source (deg)
  42 + fitspat : bool, optional
  43 + Fit spatial parameters of all models?
  44 + fitspec : bool, optional
  45 + Fit spectral parameters of all models?
  46 +
  47 + Returns
  48 + -------
  49 + model : `~gammalib.GModels`
  50 + Model container for TS fitting
  51 + """
  52 + # Create a clone of the input model
  53 + outmodels = models.copy()
  54 +
  55 + # Disable TS computation for all model components (will enable the
  56 + # test source later)
  57 + for model in outmodels:
  58 + model.tscalc(False)
  59 +
  60 + # Get source model and enable TS computation
  61 + model = outmodels[srcname]
  62 + model.tscalc(True)
  63 +
  64 + # If source model has no "Prefactor" parameter then raise an exception
  65 + if not model.has_par('Prefactor'):
  66 + msg = ('Model "%s" has no parameter "Prefactor". Only spectral '
  67 + 'models with a "Prefactor" parameter are supported.' %
  68 + srcname)
  69 + raise RuntimeError(msg)
  70 +
  71 + # Set position of test source
  72 + if ra != None and dec != None:
  73 + if model.has_par('RA') and model.has_par('DEC'):
  74 + model['RA'].value(ra)
  75 + model['DEC'].value(dec)
  76 +
  77 + # Set possible spatial and spectral parameters
  78 + spatial = ['RA', 'DEC', 'Sigma', 'Radius', 'Width', 'PA',
  79 + 'MinorRadius', 'MajorRadius']
  80 + spectral = ['Index', 'Index1', 'Index2', 'BreakEnergy', 'CutoffEnergy',
  81 + 'InverseCutoffEnergy']
  82 +
  83 + # Fit or fix spatial parameters
  84 + for par in spatial:
  85 + if model.has_par(par):
  86 + if fitspat:
  87 + model[par].free()
  88 + else:
  89 + model[par].fix()
  90 +
  91 + # Fit or fix spectral parameters
  92 + for par in spectral:
  93 + if model.has_par(par):
  94 + if fitspec:
  95 + model[par].free()
  96 + else:
  97 + model[par].fix()
  98 +
  99 + # Return model container
  100 + return outmodels
  101 +
... ...
cscripts/obsutils.py
  View file @0a5b1ab
1 1 # ==========================================================================
2   -# CTA observation handling support functions
  2 +# Utility functions for observation handling
3 3 #
4 4 # Copyright (C) 2011-2016 Juergen Knoedlseder
5 5 #
... ... @@ -277,152 +277,6 @@ def cterror(obs, srcname, log=False, debug=False, chatter=2):
277 277 return error
278 278  
279 279  
280   -# ================= #
281   -# Create counts map #
282   -# ================= #
283   -def cntmap(obs, proj='TAN', coord='GAL', xval=0.0, yval=0.0, \
284   - binsz=0.05, nxpix=200, nypix=200, \
285   - outname='cntmap.fits'):
286   - """
287   - Creates a counts map by combining the events of all observations
288   -
289   - The counts map will be a summed map over all energies
290   -
291   - Parameters
292   - ----------
293   - obs : `~gammalib.GObservations`
294   - Observation container without events
295   - proj : str, optional
296   - Projection for binned simulation
297   - coord : str, optional
298   - Coordinate system for binned simulation
299   - xval : float, optional
300   - Counts map centre of X direction (deg)
301   - yval : float, optional
302   - Counts map centre of Y direction (deg)
303   - binsz : float, optional
304   - Pixel size for binned simulation (deg/pixel)
305   - nxpix : int, optional
306   - Number of pixels in X direction
307   - nypix : int, optional
308   - Number of pixels in Y direction
309   - outname : str, optional
310   - Counts map file name
311   -
312   - Returns
313   - -------
314   - map : `~gammalib.GSkyMap`
315   - Counts map
316   - """
317   - # Allocate counts map
318   - map = gammalib.GSkyMap(proj, coord, xval, yval, -binsz, binsz, nxpix, nypix, 1)
319   -
320   - # Set maximum pixel number
321   - maxpixel = nxpix * nypix
322   -
323   - # Fill all observations
324   - for run in obs:
325   -
326   - # Loop over all events
327   - for event in run.events():
328   -
329   - # Determine sky pixel
330   - skydir = event.dir().dir()
331   - pixel = map.dir2inx(skydir)
332   -
333   - # Set pixel
334   - if pixel < maxpixel:
335   - map[pixel] += event.counts()
336   -
337   - # Save sky map. The clobber flag is set to True, so any existing FITS
338   - # file will be overwritten.
339   - map.save(outname, True)
340   -
341   - # Return counts map
342   - return map
343   -
344   -
345   -# ================ #
346   -# Create model map #
347   -# ================ #
348   -def modmap(obs, eref=0.1, proj='TAN', coord='GAL', xval=0.0, yval=0.0, \
349   - binsz=0.05, nxpix=200, nypix=200, \
350   - outname='modmap.fits'):
351   - """
352   - Make model map for a given reference energy by combining all observations
353   -
354   - The model map will be evaluated for a given reference energy 'eref' and will
355   - be given in units of [counts/(sr MeV s)].
356   -
357   - Parameters
358   - ----------
359   - obs : `~gammalib.GObservations`
360   - Observation container without events
361   - eref : float, optional
362   - Reference energy for which model is created (TeV)
363   - proj : str, optional
364   - Projection for binned simulation
365   - coord : str, optional
366   - Coordinate system for binned simulation
367   - xval : float, optional
368   - Counts map centre of X direction (deg)
369   - yval : float, optional
370   - Counts map centre of Y direction (deg)
371   - binsz : float, optional
372   - Pixel size for binned simulation (deg/pixel)
373   - nxpix : int, optional
374   - Number of pixels in X direction
375   - nypix : int, optional
376   - Number of pixels in Y direction
377   - outname : str, optional
378   - Model map file name
379   -
380   - Returns
381   - -------
382   - map : `~gammalib.GSkyMap`
383   - Model map
384   - """
385   - # Allocate model map
386   - map = gammalib.GSkyMap(proj, coord, xval, yval, -binsz, binsz, nxpix, nypix, 1)
387   -
388   - # Set reference energy, time and direction. The time is not initialised and is
389   - # in fact not used (as the IRF is assumed to be time independent for now).
390   - # The sky direction is set later using the pixel values.
391   - energy = gammalib.GEnergy()
392   - time = gammalib.GTime()
393   - instdir = gammalib.GCTAInstDir()
394   - energy.TeV(eref)
395   -
396   - # Loop over all map pixels
397   - for pixel in range(map.npix()):
398   -
399   - # Get sky direction
400   - skydir = map.inx2dir(pixel)
401   - instdir.dir(skydir)
402   -
403   - # Create event atom for map pixel
404   - atom = gammalib.GCTAEventAtom()
405   - atom.dir(instdir)
406   - atom.energy(energy)
407   - atom.time(time)
408   -
409   - # Initialise model value
410   - value = 0.0
411   -
412   - # Loop over all observations
413   - for run in obs:
414   - value += obs.models().eval(atom, run)
415   -
416   - # Set map value
417   - map[pixel] = value
418   -
419   - # Save sky map
420   - map.save(outname, True)
421   -
422   - # Return model map
423   - return map
424   -
425   -
426 280 # ======================= #
427 281 # Set one CTA observation #
428 282 # ======================= #
... ... @@ -686,82 +540,56 @@ def set_obs_patterns(pattern, ra=83.6331, dec=22.0145, offset=1.5):
686 540 return obsdeflist
687 541  
688 542  
689   -# ======================== #
690   -# Set model for TS fitting #
691   -# ======================== #
692   -def set_ts_model(models, srcname, ra=None, dec=None, fitspat=False, fitspec=False):
  543 +# ====================================================== #
  544 +# Set observation container filled with CTA observations #
  545 +# ====================================================== #
  546 +def set_observations(ra, dec, rad, tstart, duration, emin, emax, irf, caldb,
  547 + deadc=0.95, pattern='single', offset=1.5):
693 548 """
694   - Set model for TS fitting
695   -
  549 + Set an observation container filled with CTA observations
  550 +
696 551 Parameters
697 552 ----------
698   - models : `~gammalib.GModels`
699   - Input model container
700   - srcname : str
701   - Test source name
702   - ra : float, optional
703   - Right Ascension of test source (deg)
704   - dec : float, optional
705   - Declination of test source (deg)
706   - fitspat : bool, optional
707   - Fit spatial parameter?
708   - fitspec : bool, optional
709   - Fit spectral parameters?
  553 + ra : float
  554 + Right Ascension of pattern centre (deg)
  555 + dec : float
  556 + Declination of pattern centre (deg)
  557 + rad : float
  558 + ROI radius used for analysis (deg)
  559 + tstart : float
  560 + Start time of observation (seconds)
  561 + duration : float
  562 + Duration of each observation (seconds)
  563 + emin : float, optional
  564 + Minimum event energy (TeV)
  565 + emax : float, optional
  566 + Maximum event energy (TeV)
  567 + irf : str
  568 + Instrument response function
  569 + caldb : str
  570 + Calibration database path
  571 + deadc : float, optional
  572 + Deadtime correction factor
  573 + pattern : str, optional
  574 + Observation pattern ("single", "four")
  575 + offset : float, optional
  576 + Offset from pattern centre (deg)
710 577  
711 578 Returns
712 579 -------
713   - model : `~gammalib.GModels`
714   - Model container for TS fitting
  580 + obs : `~gammalib.GObservations()
  581 + Observation container
715 582 """
716   - # Create a clone of the input model
717   - outmodels = models.copy()
718   -
719   - # Disable TS computation for all model components (will enable the
720   - # test source later)
721   - for model in outmodels:
722   - model.tscalc(False)
723   -
724   - # Get source model and enable TS computation
725   - model = outmodels[srcname]
726   - model.tscalc(True)
727   -
728   - # If source model has no "Prefactor" parameter then raise an exception
729   - if not model.has_par('Prefactor'):
730   - msg = ('Model "%s" has no parameter "Prefactor". Only spectral '
731   - 'models with a "Prefactor" parameter are supported.' %
732   - srcname)
733   - raise RuntimeError(msg)
734   -
735   - # Set position of test source
736   - if ra != None and dec != None:
737   - if model.has_par('RA') and model.has_par('DEC'):
738   - model['RA'].value(ra)
739   - model['DEC'].value(dec)
740   -
741   - # Set possible spatial and spectral parameters
742   - spatial = ['RA', 'DEC', 'Sigma', 'Radius', 'Width', 'PA',
743   - 'MinorRadius', 'MajorRadius']
744   - spectral = ['Index', 'Index1', 'Index2', 'BreakEnergy', 'CutoffEnergy',
745   - 'InverseCutoffEnergy']
746   -
747   - # Fit or fix spatial parameters
748   - for par in spatial:
749   - if model.has_par(par):
750   - if fitspat:
751   - model[par].free()
752   - else:
753   - model[par].fix()
  583 + # Setup observation definition list
  584 + obsdeflist = set_obs_patterns(pattern, offset=offset, ra=ra, dec=dec)
754 585  
755   - # Fit or fix spectral parameters
756   - for par in spectral:
757   - if model.has_par(par):
758   - if fitspec:
759   - model[par].free()
760   - else:
761   - model[par].fix()
  586 + # Create list of observations
  587 + obs = set_obs_list(obsdeflist, tstart=tstart, duration=duration,
  588 + emin=emin, emax=emax, rad=rad, irf=irf, caldb=caldb,
  589 + deadc=deadc)
762 590  
763   - # Return model container
764   - return outmodels
  591 + # Return observation container
  592 + return obs
765 593  
766 594  
767 595 # ==================== #
... ... @@ -833,6 +661,8 @@ def get_stacked_response(obs, xref, yref, binsz=0.05, nxpix=200, nypix=200,
833 661 expcube['emax'] = emax
834 662 expcube['coordsys'] = coordsys
835 663 expcube['proj'] = proj
  664 + expcube['debug'] = debug
  665 + expcube['chatter'] = chatter
836 666 if log:
837 667 expcube.logFileOpen()
838 668 expcube.run()
... ... @@ -852,6 +682,8 @@ def get_stacked_response(obs, xref, yref, binsz=0.05, nxpix=200, nypix=200,
852 682 psfcube['emax'] = emax
853 683 psfcube['coordsys'] = coordsys
854 684 psfcube['proj'] = proj
  685 + psfcube['debug'] = debug
  686 + psfcube['chatter'] = chatter
855 687 if log:
856 688 psfcube.logFileOpen()
857 689 psfcube.run()
... ... @@ -871,6 +703,8 @@ def get_stacked_response(obs, xref, yref, binsz=0.05, nxpix=200, nypix=200,
871 703 bkgcube['emax'] = emax
872 704 bkgcube['coordsys'] = coordsys
873 705 bkgcube['proj'] = proj
  706 + bkgcube['debug'] = debug
  707 + bkgcube['chatter'] = chatter
874 708 if log:
875 709 bkgcube.logFileOpen()
876 710 bkgcube.run()
... ... @@ -891,6 +725,8 @@ def get_stacked_response(obs, xref, yref, binsz=0.05, nxpix=200, nypix=200,
891 725 edispcube['emax'] = emax
892 726 edispcube['coordsys'] = coordsys
893 727 edispcube['proj'] = proj
  728 + edispcube['debug'] = debug
  729 + edispcube['chatter'] = chatter
894 730 if log:
895 731 edispcube.logFileOpen()
896 732 edispcube.run()
... ...