diff --git a/ChangeLog b/ChangeLog
index 7004570..6a78cb8 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,8 +1,9 @@
-2016-07-29
+2016-07-31
* Version 1.2.0 released
========================
+ Add support modules to cscripts (#1830)
Fix lacking coverage measure for cscripts execution (#1827)
Set "addbounds" parameters by default to "no" (#1819)
diff --git a/NEWS b/NEWS
index bd3dddf..0316dc6 100644
--- a/NEWS
+++ b/NEWS
@@ -1,6 +1,6 @@
New Features and Important Changes in ctools 1.2.0
-29 July 2016
+31 July 2016
Introduction
@@ -105,9 +105,11 @@ cterror - Likelihood profile error computation
None
-scripts
--------
-None
+cscripts
+--------
+Added the support modules "modutils" and "ioutils" to the cscripts module
+and removed the obsolete function "cntmap()" and "modmap()" from the
+"obsutils" module (#1830).
Examples
diff --git a/README.md b/README.md
index 8722dce..03f1f1b 100644
--- a/README.md
+++ b/README.md
@@ -1,6 +1,6 @@
ctools information
==================
-* Version: 1.2.0.dev1 (29 July 2016)
+* Version: 1.2.0.dev1 (31 July 2016)
* Author: Juergen Knoedlseder (jurgen.knodlseder@irap.omp.eu)
* GammaLib dependency: 1.2.0.dev1
diff --git a/cscripts/Makefile.am b/cscripts/Makefile.am
index dc0bece..fa8d588 100644
--- a/cscripts/Makefile.am
+++ b/cscripts/Makefile.am
@@ -42,7 +42,9 @@ pythonscripts = $(builddir)/__init__.py \
$(srcdir)/csiactcopy.py \
$(srcdir)/csiactdata.py \
$(srcdir)/csiactobs.py \
- $(srcdir)/obsutils.py
+ $(srcdir)/obsutils.py \
+ $(srcdir)/modutils.py \
+ $(srcdir)/ioutils.py
# Parfiles to be distributed in $(prefix)/syspfiles
parfiles = $(srcdir)/cscaldb.par \
@@ -203,6 +205,8 @@ clean-local:
@if test '@srcdir@' != '@builddir@'; then \
rm -rf $(builddir)/cs*.py; \
rm -rf $(builddir)/obsutils.py; \
+ rm -rf $(builddir)/modutils.py; \
+ rm -rf $(builddir)/ioutils.py; \
rm -rf $(builddir)/*.par; \
fi
rm -rf $(builddir)/*.pyc
diff --git a/cscripts/__init__.py.in b/cscripts/__init__.py.in
index 0918a44..d083943 100644
--- a/cscripts/__init__.py.in
+++ b/cscripts/__init__.py.in
@@ -41,7 +41,9 @@ __all__ = [
"cstsmapmerge",
"cstsmapsplit",
"csworkflow",
- "obsutils"
+ "obsutils",
+ "modutils",
+ "ioutils"
]
if cur_version > req_version:
__all__.extend(["csfindobs",
@@ -68,6 +70,8 @@ from cscripts.cstsmapmerge import cstsmapmerge
from cscripts.cstsmapsplit import cstsmapsplit
from cscripts.csworkflow import csworkflow
from cscripts import obsutils
+from cscripts import modutils
+from cscripts import ioutils
if cur_version > req_version:
from cscripts.csfindobs import csfindobs
from cscripts.csiactcopy import csiactcopy
diff --git a/cscripts/cspull.py b/cscripts/cspull.py
index a49e478..6f3ab54 100755
--- a/cscripts/cspull.py
+++ b/cscripts/cspull.py
@@ -19,11 +19,11 @@
#
# ==========================================================================
import sys
-import csv
import math
import gammalib
import ctools
from cscripts import obsutils
+from cscripts import ioutils
# ============ #
@@ -79,7 +79,18 @@ class cspull(ctools.cscript):
# observation parameters to set wobble pattern
self._pattern = self['pattern'].string()
if self._pattern == 'four':
- self._obs = self._set_obs()
+ self._obs = obsutils.set_observations(self['ra'].real(),
+ self['dec'].real(),
+ self['rad'].real(),
+ self['tmin'].real(),
+ self['tmax'].real(),
+ self['emin'].real(),
+ self['emax'].real(),
+ self['irf'].string(),
+ self['caldb'].string(),
+ deadc=self['deadc'].real(),
+ pattern=self['pattern'].string(),
+ offset=self['offset'].real())
# ... otherwise add response information and energy boundaries
# in case they are missing.
@@ -119,38 +130,6 @@ class cspull(ctools.cscript):
# Return
return
- def _set_obs(self):
- """
- Set CTA observation(s)
-
- Returns
- -------
- obs : `~gammalib.GObservations`
- Observation container.
- """
- # Set duration
- duration = self['tmax'].real() - self['tmin'].real()
-
- # Setup observation definition list
- obsdeflist = obsutils.set_obs_patterns(self._pattern,
- ra=self['ra'].real(),
- dec=self['dec'].real(),
- offset=self['offset'].real())
-
- # Create list of observations
- obs = obsutils.set_obs_list(obsdeflist,
- tstart=self['tmin'].real(),
- duration=duration,
- deadc=self['deadc'].real(),
- emin=self['emin'].real(),
- emax=self['emax'].real(),
- rad=self['rad'].real(),
- irf=self['irf'].string(),
- caldb=self['caldb'].string())
-
- # Return observation container
- return obs
-
def _obs_string(self, obs):
"""
Generate summary string for observation
@@ -359,22 +338,9 @@ class cspull(ctools.cscript):
# Make a trial and add initial seed
result = self._trial(seed + self._seed)
- # Set output filename
- outfile = self['outfile'].filename().url()
-
- # Write out result immediately
- if seed == 0:
- f = open(outfile, 'w')
- writer = csv.DictWriter(f, result['colnames'])
- headers = {}
- for n in result['colnames']:
- headers[n] = n
- writer.writerow(headers)
- else:
- f = open(outfile, 'a')
- writer = csv.DictWriter(f, result['colnames'])
- writer.writerow(result['values'])
- f.close()
+ # Write out trial result
+ ioutils.write_csv_row(self['outfile'].filename().url(), seed,
+ result['colnames'], result['values'])
# Return
return
diff --git a/cscripts/cssens.py b/cscripts/cssens.py
index 9828d64..6e73a20 100755
--- a/cscripts/cssens.py
+++ b/cscripts/cssens.py
@@ -1,6 +1,6 @@
#! /usr/bin/env python
# ==========================================================================
-# This script computes the CTA sensitivity.
+# Computes the array sensitivity using a test source
#
# Copyright (C) 2011-2016 Juergen Knoedlseder
#
@@ -24,6 +24,8 @@ import math
import gammalib
import ctools
from cscripts import obsutils
+from cscripts import modutils
+from cscripts import ioutils
# ============ #
@@ -190,70 +192,6 @@ class cssens(ctools.cscript):
# Return observation container
return obs
- def _set_models(self, fitspat=False, fitspec=False):
- """
- Set full model and background model
-
- Parameters
- ----------
- fitspat : bool, optional
- Fit spatial parameter?
- fitspec : bool, optional
- Fit spectral parameters?
-
- Returns
- -------
- full_model, bkg_model : `~gammalib.GModels`
- Full model and background model containers
- """
- # Retrieve full model from observation container
- full_model = self._obs.models().copy()
-
- # Get source model
- model = full_model[self._srcname]
-
- # If source model has no "Prefactor" parameter then raise an exception
- if not model.has_par('Prefactor'):
- msg = ('Model "%s" has no parameter "Prefactor". Only spectral '
- 'models with a "Prefactor" parameter are supported.' %
- self._srcname)
- raise RuntimeError(msg)
-
- # Set source position
- if self._ra != None and self._dec != None:
- if model.has_par('RA') and model.has_par('DEC'):
- model['RA'].value(self._ra)
- model['DEC'].value(self._dec)
-
- # Set possible spatial and spectral parameters
- spatial = ['RA', 'DEC', 'Sigma', 'Radius', 'Width', 'PA',
- 'MinorRadius', 'MajorRadius']
- spectral = ['Index', 'Index1', 'Index2', 'BreakEnergy', 'CutoffEnergy',
- 'InverseCutoffEnergy']
-
- # Fit or fix spatial parameters
- for par in spatial:
- if model.has_par(par):
- if fitspat:
- model[par].free()
- else:
- model[par].fix()
-
- # Fit or fix spectral parameters
- for par in spectral:
- if model.has_par(par):
- if fitspec:
- model[par].free()
- else:
- model[par].fix()
-
- # Create background model
- bkg_model = full_model.copy()
- bkg_model.remove(self._srcname)
-
- # Return models
- return full_model, bkg_model
-
def _set_obs_ebounds(self, emin, emax):
"""
Set energy boundaries for all observations in container
@@ -330,9 +268,9 @@ class cssens(ctools.cscript):
# Return photon flux
return flux
- def _get_sensitivity(self, emin, emax, bkg_model, full_model):
+ def _get_sensitivity(self, emin, emax, test_model):
"""
- Determine sensitivity for given observations.
+ Determine sensitivity for given observations
Parameters
----------
@@ -340,10 +278,8 @@ class cssens(ctools.cscript):
Minimum energy for fitting and flux computation
emax : `~gammalib.GEnergy`
Maximum energy for fitting and flux computation
- bkg_model : `~gammalib.GModels`
- Background model
- full_model : `~gammalib.GModels`
- Source plus background model
+ test_model : `~gammalib.GModels`
+ Test source model
Returns
-------
@@ -364,7 +300,7 @@ class cssens(ctools.cscript):
# Compute Crab unit (this is the factor with which the Prefactor needs
# to be multiplied to get 1 Crab
crab_flux = self._get_crab_flux(emin, emax)
- src_flux = full_model[self._srcname].spectral().flux(emin, emax)
+ src_flux = test_model[self._srcname].spectral().flux(emin, emax)
crab_unit = crab_flux/src_flux
# Write header
@@ -399,14 +335,14 @@ class cssens(ctools.cscript):
if self._logExplicit():
self._log.header2('Iteration '+str(iterations))
- # Set source model. crab_prefactor is the Prefactor that
+ # Set test source model. crab_prefactor is the Prefactor that
# corresponds to 1 Crab
- src_model = full_model.copy()
- crab_prefactor = src_model[self._srcname]['Prefactor'].value() * \
+ models = test_model.copy()
+ crab_prefactor = models[self._srcname]['Prefactor'].value() * \
crab_unit
- src_model[self._srcname]['Prefactor'].value(crab_prefactor * \
- test_crab_flux)
- self._obs.models(src_model)
+ models[self._srcname]['Prefactor'].value(crab_prefactor * \
+ test_crab_flux)
+ self._obs.models(models)
# Simulate events
sim = obsutils.sim(self._obs, nbins=self._enumbins, seed=iterations,
@@ -426,68 +362,35 @@ class cssens(ctools.cscript):
self._log(nevents)
self._log('\n')
- # Fit background only
- sim.models(bkg_model)
- like = obsutils.fit(sim, log=self._log_clients,
- debug=self._debug, edisp=self._edisp)
- result_bgm = like.obs().models().copy()
- LogL_bgm = like.opt().value()
- npred_bgm = like.obs().npred()
-
- # Assess quality based on a comparison between Npred and Nevents
- quality_bgm = npred_bgm-nevents
-
- # Write background fit results
- if self._logExplicit():
- self._log.header3('Background model fit')
- self._log.parformat('log likelihood')
- self._log(LogL_bgm)
- self._log('\n')
- self._log.parformat('Number of predicted events')
- self._log(npred_bgm)
- self._log('\n')
- self._log.parformat('Fit quality')
- self._log(quality_bgm)
- self._log('\n')
-
- # Write model fit results
- if self._logExplicit():
- for model in result_bgm:
- self._log.parformat('Model')
- self._log(model.name())
- self._log('\n')
- for par in model:
- self._log(str(par)+'\n')
+ # Fit test source
+ fit = obsutils.fit(sim, log=self._log_clients, debug=self._debug,
+ edisp=self._edisp)
- # Fit background and test source
- sim.models(src_model)
- like = obsutils.fit(sim, log=self._log_clients,
- debug=self._debug, edisp=self._edisp)
- result_all = like.obs().models().copy()
- LogL_all = like.opt().value()
- npred_all = like.obs().npred()
- ts = 2.0*(LogL_bgm-LogL_all)
+ # Get model fitting results
+ logL = fit.opt().value()
+ npred = fit.obs().npred()
+ models = fit.obs().models()
+ ts = models[self._srcname].ts()
# Assess quality based on a comparison between Npred and Nevents
- quality_all = npred_all-nevents
+ quality = npred - nevents
- # Write background and test source fit results
+ # Write test source fit results
if self._logExplicit():
- self._log.header3('Background and test source model fit')
+ self._log.header3('Test source model fit')
self._log.parformat('Test statistics')
self._log(ts)
self._log('\n')
self._log.parformat('log likelihood')
- self._log(LogL_all)
+ self._log(logL)
self._log('\n')
self._log.parformat('Number of predicted events')
- self._log(npred_all)
+ self._log(npred)
self._log('\n')
self._log.parformat('Fit quality')
- self._log(quality_all)
+ self._log(quality)
self._log('\n')
- #
- for model in result_all:
+ for model in models:
self._log.parformat('Model')
self._log(model.name())
self._log('\n')
@@ -501,17 +404,14 @@ class cssens(ctools.cscript):
continue
# Get fitted Crab, photon and energy fluxes
- crab_flux = result_all[self._srcname]['Prefactor'].value() / \
- crab_prefactor
- #crab_flux_err = result_all[self._srcname]['Prefactor'].error() / \
- # crab_prefactor
- photon_flux = result_all[self._srcname].spectral().flux(emin, emax)
- energy_flux = result_all[self._srcname].spectral().eflux(emin, emax)
+ crab_flux = models[self._srcname]['Prefactor'].value() / \
+ crab_prefactor
+ photon_flux = models[self._srcname].spectral().flux(emin, emax)
+ energy_flux = models[self._srcname].spectral().eflux(emin, emax)
# Compute differential sensitivity in unit erg/cm2/s
energy = gammalib.GEnergy(e_mean, 'TeV')
- time = gammalib.GTime()
- sensitivity = result_all[self._srcname].spectral().eval(energy, time) * \
+ sensitivity = models[self._srcname].spectral().eval(energy) * \
e_mean*erg_mean * 1.0e6
# Compute flux correction factor based on average TS
@@ -543,7 +443,7 @@ class cssens(ctools.cscript):
self._log.parformat('Differential sensitivity')
self._log(sensitivity)
self._log(' erg/cm2/s\n')
- for model in result_all:
+ for model in models:
self._log.parformat('Model')
self._log(model.name())
self._log('\n')
@@ -636,27 +536,14 @@ class cssens(ctools.cscript):
self._log.parformat('Number of simulated events')
self._log(nevents)
self._log('\n')
- self._log.header3('Background and test source model fitting')
+ self._log.header3('Test source model fitting')
self._log.parformat('log likelihood')
- self._log(LogL_all)
+ self._log(logL)
self._log('\n')
self._log.parformat('Number of predicted events')
- self._log(npred_all)
+ self._log(npred)
self._log('\n')
- for model in result_all:
- self._log.parformat('Model')
- self._log(model.name())
- self._log('\n')
- for par in model:
- self._log(str(par)+'\n')
- self._log.header3('Background model fit')
- self._log.parformat('log likelihood')
- self._log(LogL_bgm)
- self._log('\n')
- self._log.parformat('Number of predicted events')
- self._log(npred_bgm)
- self._log('\n')
- for model in result_bgm:
+ for model in models:
self._log.parformat('Model')
self._log(model.name())
self._log('\n')
@@ -698,18 +585,14 @@ class cssens(ctools.cscript):
'energy_flux', 'sensitivity']
results = []
- # Initialise models
- full_model, bkg_model = self._set_models()
+ # Set test source model for this observation
+ models = modutils.test_source(self._obs.models(), self._srcname)
# Write models into logger
if self._logTerse():
self._log('\n')
self._log.header1('Models')
- self._log.header2('Background model')
- self._log(str(bkg_model))
- self._log('\n\n')
- self._log.header2('Full model')
- self._log(str(full_model))
+ self._log(str(models))
self._log('\n')
# Write header
@@ -736,23 +619,10 @@ class cssens(ctools.cscript):
raise RuntimeError(msg)
# Determine sensitivity
- result = self._get_sensitivity(emin, emax, bkg_model, full_model)
-
- # Write results
- if ieng == 0:
- f = open(self._outfile.url(), 'w')
- writer = csv.DictWriter(f, colnames)
- headers = {}
- for n in colnames:
- headers[n] = n
- writer.writerow(headers)
- writer.writerow(result)
- f.close()
- else:
- f = open(self._outfile.url(), 'a')
- writer = csv.DictWriter(f, colnames)
- writer.writerow(result)
- f.close()
+ result = self._get_sensitivity(emin, emax, models)
+
+ # Write out trial result
+ ioutils.write_csv_row(self._outfile.url(), ieng, colnames, result)
# Append results
results.append(result)
diff --git a/cscripts/cstsdist.py b/cscripts/cstsdist.py
index 7f3e2f7..e6f6d23 100755
--- a/cscripts/cstsdist.py
+++ b/cscripts/cstsdist.py
@@ -19,10 +19,11 @@
#
# ==========================================================================
import sys
-import csv
import gammalib
import ctools
from cscripts import obsutils
+from cscripts import modutils
+from cscripts import ioutils
# ============== #
@@ -68,7 +69,18 @@ class cstsdist(ctools.cscript):
# Check for requested pattern and use above observation parameters
# to set wobble pattern
if self['pattern'].string() == 'four':
- self._obs = self._set_obs()
+ self._obs = obsutils.set_observations(self['ra'].real(),
+ self['dec'].real(),
+ self['rad'].real(),
+ self['tmin'].real(),
+ self['tmax'].real(),
+ self['emin'].real(),
+ self['emax'].real(),
+ self['irf'].string(),
+ self['caldb'].string(),
+ deadc=self['deadc'].real(),
+ pattern=self['pattern'].string(),
+ offset=self['offset'].real())
# Get source name
self._srcname = self['srcname'].string()
@@ -98,38 +110,6 @@ class cstsdist(ctools.cscript):
# Return
return
- def _set_obs(self):
- """
- Set an observation container
-
- Returns
- -------
- obs : `~gammalib.GObservations`
- Observation container
- """
- # Set duration
- duration = self['tmax'].real() - self['tmin'].real()
-
- # Setup observation definition list
- obsdeflist = obsutils.set_obs_patterns(self['pattern'].string(),
- ra = self['ra'].real(),
- dec = self['dec'].real(),
- offset = self['offset'].real())
-
- # Create list of observations
- obs = obsutils.set_obs_list(obsdeflist,
- tstart = self['tmin'].real(),
- duration = duration,
- deadc = self['deadc'].real(),
- emin = self['emin'].real(),
- emax = self['emax'].real(),
- rad = self['rad'].real(),
- irf = self['irf'].string(),
- caldb = self['caldb'].string())
-
- # Return observation container
- return obs
-
def _trial(self, seed):
"""
Create the TS for a single trial
@@ -276,8 +256,8 @@ class cstsdist(ctools.cscript):
# Get parameters
self._get_parameters()
- # Set models for this observation
- self.models(obsutils.set_ts_model(self._obs.models(), self._srcname))
+ # Set test source model for this observation
+ self.models(modutils.test_source(self._obs.models(), self._srcname))
# Write models into logger
if self._logTerse():
@@ -308,22 +288,9 @@ class cstsdist(ctools.cscript):
# Make a trial
result = self._trial(seed)
- # Set output filename
- outfile = self['outfile'].filename().url()
-
- # Write out result immediately
- if seed == 0:
- f = open(outfile, 'w')
- writer = csv.DictWriter(f, result['colnames'])
- headers = {}
- for n in result['colnames']:
- headers[n] = n
- writer.writerow(headers)
- else:
- f = open(outfile, 'a')
- writer = csv.DictWriter(f, result['colnames'])
- writer.writerow(result['values'])
- f.close()
+ # Write out trial result
+ ioutils.write_csv_row(self['outfile'].filename().url(), seed,
+ result['colnames'], result['values'])
# Return
return
diff --git a/cscripts/ioutils.py b/cscripts/ioutils.py
new file mode 100755
index 0000000..13d2646
--- /dev/null
+++ b/cscripts/ioutils.py
@@ -0,0 +1,64 @@
+# ==========================================================================
+# Utility functions for input and output
+#
+# Copyright (C) 2016 Juergen Knoedlseder
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+#
+# ==========================================================================
+import csv
+import gammalib
+import ctools
+
+
+# ========================= #
+# Write one row of CSV file #
+# ========================= #
+def write_csv_row(outfile, row, colnames, colvalues):
+ """
+ Write one row into CSV file
+
+ Parameters
+ ----------
+ outfile : str
+ Output file name
+ row : int
+ Row number
+ colnames : list of str
+ Column names
+ colvalues : list of float
+ Column values
+ """
+ # If this is the first row then open a new file and write the header
+ if row == 0:
+ f = open(outfile, 'w')
+ writer = csv.DictWriter(f, colnames)
+ headers = {}
+ for colname in colnames:
+ headers[colname] = colname
+ writer.writerow(headers)
+
+ # ... otherwise append to an existing file
+ else:
+ f = open(outfile, 'a')
+
+ # Write out row
+ writer = csv.DictWriter(f, colnames)
+ writer.writerow(colvalues)
+
+ # Close file
+ f.close()
+
+ # Return
+ return
diff --git a/cscripts/modutils.py b/cscripts/modutils.py
new file mode 100755
index 0000000..2e05cfe
--- /dev/null
+++ b/cscripts/modutils.py
@@ -0,0 +1,101 @@
+# ==========================================================================
+# Utility functions for model handling
+#
+# Copyright (C) 2016 Juergen Knoedlseder
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+#
+# ==========================================================================
+import gammalib
+import ctools
+
+
+# ============================================ #
+# Return model for TS fitting of a test source #
+# ============================================ #
+def test_source(models, srcname, ra=None, dec=None, fitspat=False,
+ fitspec=False):
+ """
+ Return model for TS fitting of a test source
+
+ Parameters
+ ----------
+ models : `~gammalib.GModels`
+ Input model container
+ srcname : str
+ Test source name
+ ra : float, optional
+ Right Ascension of test source (deg)
+ dec : float, optional
+ Declination of test source (deg)
+ fitspat : bool, optional
+ Fit spatial parameters of all models?
+ fitspec : bool, optional
+ Fit spectral parameters of all models?
+
+ Returns
+ -------
+ model : `~gammalib.GModels`
+ Model container for TS fitting
+ """
+ # Create a clone of the input model
+ outmodels = models.copy()
+
+ # Disable TS computation for all model components (will enable the
+ # test source later)
+ for model in outmodels:
+ model.tscalc(False)
+
+ # Get source model and enable TS computation
+ model = outmodels[srcname]
+ model.tscalc(True)
+
+ # If source model has no "Prefactor" parameter then raise an exception
+ if not model.has_par('Prefactor'):
+ msg = ('Model "%s" has no parameter "Prefactor". Only spectral '
+ 'models with a "Prefactor" parameter are supported.' %
+ srcname)
+ raise RuntimeError(msg)
+
+ # Set position of test source
+ if ra != None and dec != None:
+ if model.has_par('RA') and model.has_par('DEC'):
+ model['RA'].value(ra)
+ model['DEC'].value(dec)
+
+ # Set possible spatial and spectral parameters
+ spatial = ['RA', 'DEC', 'Sigma', 'Radius', 'Width', 'PA',
+ 'MinorRadius', 'MajorRadius']
+ spectral = ['Index', 'Index1', 'Index2', 'BreakEnergy', 'CutoffEnergy',
+ 'InverseCutoffEnergy']
+
+ # Fit or fix spatial parameters
+ for par in spatial:
+ if model.has_par(par):
+ if fitspat:
+ model[par].free()
+ else:
+ model[par].fix()
+
+ # Fit or fix spectral parameters
+ for par in spectral:
+ if model.has_par(par):
+ if fitspec:
+ model[par].free()
+ else:
+ model[par].fix()
+
+ # Return model container
+ return outmodels
+
diff --git a/cscripts/obsutils.py b/cscripts/obsutils.py
index d6b2052..0cd848f 100755
--- a/cscripts/obsutils.py
+++ b/cscripts/obsutils.py
@@ -1,5 +1,5 @@
# ==========================================================================
-# CTA observation handling support functions
+# Utility functions for observation handling
#
# Copyright (C) 2011-2016 Juergen Knoedlseder
#
@@ -277,152 +277,6 @@ def cterror(obs, srcname, log=False, debug=False, chatter=2):
return error
-# ================= #
-# Create counts map #
-# ================= #
-def cntmap(obs, proj='TAN', coord='GAL', xval=0.0, yval=0.0, \
- binsz=0.05, nxpix=200, nypix=200, \
- outname='cntmap.fits'):
- """
- Creates a counts map by combining the events of all observations
-
- The counts map will be a summed map over all energies
-
- Parameters
- ----------
- obs : `~gammalib.GObservations`
- Observation container without events
- proj : str, optional
- Projection for binned simulation
- coord : str, optional
- Coordinate system for binned simulation
- xval : float, optional
- Counts map centre of X direction (deg)
- yval : float, optional
- Counts map centre of Y direction (deg)
- binsz : float, optional
- Pixel size for binned simulation (deg/pixel)
- nxpix : int, optional
- Number of pixels in X direction
- nypix : int, optional
- Number of pixels in Y direction
- outname : str, optional
- Counts map file name
-
- Returns
- -------
- map : `~gammalib.GSkyMap`
- Counts map
- """
- # Allocate counts map
- map = gammalib.GSkyMap(proj, coord, xval, yval, -binsz, binsz, nxpix, nypix, 1)
-
- # Set maximum pixel number
- maxpixel = nxpix * nypix
-
- # Fill all observations
- for run in obs:
-
- # Loop over all events
- for event in run.events():
-
- # Determine sky pixel
- skydir = event.dir().dir()
- pixel = map.dir2inx(skydir)
-
- # Set pixel
- if pixel < maxpixel:
- map[pixel] += event.counts()
-
- # Save sky map. The clobber flag is set to True, so any existing FITS
- # file will be overwritten.
- map.save(outname, True)
-
- # Return counts map
- return map
-
-
-# ================ #
-# Create model map #
-# ================ #
-def modmap(obs, eref=0.1, proj='TAN', coord='GAL', xval=0.0, yval=0.0, \
- binsz=0.05, nxpix=200, nypix=200, \
- outname='modmap.fits'):
- """
- Make model map for a given reference energy by combining all observations
-
- The model map will be evaluated for a given reference energy 'eref' and will
- be given in units of [counts/(sr MeV s)].
-
- Parameters
- ----------
- obs : `~gammalib.GObservations`
- Observation container without events
- eref : float, optional
- Reference energy for which model is created (TeV)
- proj : str, optional
- Projection for binned simulation
- coord : str, optional
- Coordinate system for binned simulation
- xval : float, optional
- Counts map centre of X direction (deg)
- yval : float, optional
- Counts map centre of Y direction (deg)
- binsz : float, optional
- Pixel size for binned simulation (deg/pixel)
- nxpix : int, optional
- Number of pixels in X direction
- nypix : int, optional
- Number of pixels in Y direction
- outname : str, optional
- Model map file name
-
- Returns
- -------
- map : `~gammalib.GSkyMap`
- Model map
- """
- # Allocate model map
- map = gammalib.GSkyMap(proj, coord, xval, yval, -binsz, binsz, nxpix, nypix, 1)
-
- # Set reference energy, time and direction. The time is not initialised and is
- # in fact not used (as the IRF is assumed to be time independent for now).
- # The sky direction is set later using the pixel values.
- energy = gammalib.GEnergy()
- time = gammalib.GTime()
- instdir = gammalib.GCTAInstDir()
- energy.TeV(eref)
-
- # Loop over all map pixels
- for pixel in range(map.npix()):
-
- # Get sky direction
- skydir = map.inx2dir(pixel)
- instdir.dir(skydir)
-
- # Create event atom for map pixel
- atom = gammalib.GCTAEventAtom()
- atom.dir(instdir)
- atom.energy(energy)
- atom.time(time)
-
- # Initialise model value
- value = 0.0
-
- # Loop over all observations
- for run in obs:
- value += obs.models().eval(atom, run)
-
- # Set map value
- map[pixel] = value
-
- # Save sky map
- map.save(outname, True)
-
- # Return model map
- return map
-
-
# ======================= #
# Set one CTA observation #
# ======================= #
@@ -686,82 +540,56 @@ def set_obs_patterns(pattern, ra=83.6331, dec=22.0145, offset=1.5):
return obsdeflist
-# ======================== #
-# Set model for TS fitting #
-# ======================== #
-def set_ts_model(models, srcname, ra=None, dec=None, fitspat=False, fitspec=False):
+# ====================================================== #
+# Set observation container filled with CTA observations #
+# ====================================================== #
+def set_observations(ra, dec, rad, tstart, duration, emin, emax, irf, caldb,
+ deadc=0.95, pattern='single', offset=1.5):
"""
- Set model for TS fitting
-
+ Set an observation container filled with CTA observations
+
Parameters
----------
- models : `~gammalib.GModels`
- Input model container
- srcname : str
- Test source name
- ra : float, optional
- Right Ascension of test source (deg)
- dec : float, optional
- Declination of test source (deg)
- fitspat : bool, optional
- Fit spatial parameter?
- fitspec : bool, optional
- Fit spectral parameters?
+ ra : float
+ Right Ascension of pattern centre (deg)
+ dec : float
+ Declination of pattern centre (deg)
+ rad : float
+ ROI radius used for analysis (deg)
+ tstart : float
+ Start time of observation (seconds)
+ duration : float
+ Duration of each observation (seconds)
+ emin : float, optional
+ Minimum event energy (TeV)
+ emax : float, optional
+ Maximum event energy (TeV)
+ irf : str
+ Instrument response function
+ caldb : str
+ Calibration database path
+ deadc : float, optional
+ Deadtime correction factor
+ pattern : str, optional
+ Observation pattern ("single", "four")
+ offset : float, optional
+ Offset from pattern centre (deg)
Returns
-------
- model : `~gammalib.GModels`
- Model container for TS fitting
+ obs : `~gammalib.GObservations()
+ Observation container
"""
- # Create a clone of the input model
- outmodels = models.copy()
-
- # Disable TS computation for all model components (will enable the
- # test source later)
- for model in outmodels:
- model.tscalc(False)
-
- # Get source model and enable TS computation
- model = outmodels[srcname]
- model.tscalc(True)
-
- # If source model has no "Prefactor" parameter then raise an exception
- if not model.has_par('Prefactor'):
- msg = ('Model "%s" has no parameter "Prefactor". Only spectral '
- 'models with a "Prefactor" parameter are supported.' %
- srcname)
- raise RuntimeError(msg)
-
- # Set position of test source
- if ra != None and dec != None:
- if model.has_par('RA') and model.has_par('DEC'):
- model['RA'].value(ra)
- model['DEC'].value(dec)
-
- # Set possible spatial and spectral parameters
- spatial = ['RA', 'DEC', 'Sigma', 'Radius', 'Width', 'PA',
- 'MinorRadius', 'MajorRadius']
- spectral = ['Index', 'Index1', 'Index2', 'BreakEnergy', 'CutoffEnergy',
- 'InverseCutoffEnergy']
-
- # Fit or fix spatial parameters
- for par in spatial:
- if model.has_par(par):
- if fitspat:
- model[par].free()
- else:
- model[par].fix()
+ # Setup observation definition list
+ obsdeflist = set_obs_patterns(pattern, offset=offset, ra=ra, dec=dec)
- # Fit or fix spectral parameters
- for par in spectral:
- if model.has_par(par):
- if fitspec:
- model[par].free()
- else:
- model[par].fix()
+ # Create list of observations
+ obs = set_obs_list(obsdeflist, tstart=tstart, duration=duration,
+ emin=emin, emax=emax, rad=rad, irf=irf, caldb=caldb,
+ deadc=deadc)
- # Return model container
- return outmodels
+ # Return observation container
+ return obs
# ==================== #
@@ -833,6 +661,8 @@ def get_stacked_response(obs, xref, yref, binsz=0.05, nxpix=200, nypix=200,
expcube['emax'] = emax
expcube['coordsys'] = coordsys
expcube['proj'] = proj
+ expcube['debug'] = debug
+ expcube['chatter'] = chatter
if log:
expcube.logFileOpen()
expcube.run()
@@ -852,6 +682,8 @@ def get_stacked_response(obs, xref, yref, binsz=0.05, nxpix=200, nypix=200,
psfcube['emax'] = emax
psfcube['coordsys'] = coordsys
psfcube['proj'] = proj
+ psfcube['debug'] = debug
+ psfcube['chatter'] = chatter
if log:
psfcube.logFileOpen()
psfcube.run()
@@ -871,6 +703,8 @@ def get_stacked_response(obs, xref, yref, binsz=0.05, nxpix=200, nypix=200,
bkgcube['emax'] = emax
bkgcube['coordsys'] = coordsys
bkgcube['proj'] = proj
+ bkgcube['debug'] = debug
+ bkgcube['chatter'] = chatter
if log:
bkgcube.logFileOpen()
bkgcube.run()
@@ -891,6 +725,8 @@ def get_stacked_response(obs, xref, yref, binsz=0.05, nxpix=200, nypix=200,
edispcube['emax'] = emax
edispcube['coordsys'] = coordsys
edispcube['proj'] = proj
+ edispcube['debug'] = debug
+ edispcube['chatter'] = chatter
if log:
edispcube.logFileOpen()
edispcube.run()