diff --git a/ChangeLog b/ChangeLog
index 9bdb0cf..57ccf57 100755
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,4 +1,21 @@
-2018-01-25
+2018-02-14
+
+ * Version 1.6.0 released
+ ========================
+
+ Allow specifying of input sky map in ctskymap (#2311)
+ Add ras() and decs() methods to csobsinfo (#2313)
+ Accelerate IRF and RING background methods in ctskymap (#2309)
+
+
+2018-02-16
+
+ * Bug fix version 1.5.1 released
+ ================================
+
+ Take energy dispersion in csresspec modelling into account (#2341)
+ Enable using csphagen for HESS, VERITAS and MAGIC observations (#2307)
+
* Version 1.5.0 released
========================
diff --git a/NEWS b/NEWS
index 921fdff..ae35b87 100755
--- a/NEWS
+++ b/NEWS
@@ -1,4 +1,6 @@
-25 January 2018
+New Features and Important Changes in ctools 1.6.0
+
+14 February 2018
Introduction
@@ -9,10 +11,7 @@ ctools release with respect to the precedent release.
Configuration
-------------
-In case of an error, test output is now written into log file (#2264).
-
-The cscripts::execute() method no longer creates log files, making the
-cscripts behaviour compliant with that of ctools (#2182).
+None
Python interface
@@ -27,8 +26,7 @@ None
Base classes
------------
-Add ctools::time_reference constant and remove ctools::time_reference() method
-(#2295).
+None
ctbin - CTA event binning
@@ -43,9 +41,7 @@ None
ctbutterfly - Generate butterfly diagrams
-----------------------------------------
-Add "statistic" parameter to allow setting of optimizer statistic (#2237).
-
-Add butterfly() method to return the butterfly diagram as CSV object (#1512).
+None
ctcubemask - CTA cube filtering
@@ -60,7 +56,7 @@ None
cterror - Likelihood profile error computation
----------------------------------------------
-Add "statistic" parameter to allow setting of optimizer statistic (#2237).
+None
ctexpcube - CTA exposure cube generation
@@ -70,7 +66,7 @@ None
ctlike - Maximum likelihood fitting
-----------------------------------
-Rename "stats" parameter to "statistic" (#2237).
+None
ctmapcube - CTA map cube generation
@@ -80,22 +76,12 @@ None
ctmodel - CTA counts cube modelling
-----------------------------------
-Use t-type parameters for "tmin" and "tmax" to all specifying times in various
-formats (#1864).
-
-The tool now also checks which models actually overlap with the RoI of a given
-observation and then only use these models to compute the model cube. The
-parallelization of the tool has also been optimised so that diffuse models are
-not loaded over and over again (#2189).
+None
ctobssim - CTA observation simulation
-------------------------------------
-Use t-type parameters for "tmin" and "tmax" to all specifying times in various
-formats (#1864).
-
-Set the OBS_ID keyword to the number that is used for creating the FITS file
-(#1836).
+None
ctphase - Event phase computation
@@ -110,51 +96,47 @@ None
ctpsfcube - CTA point spread function cube generation
-----------------------------------------------------
-The chatter level is now correctly forwarded to the PSF cube fill method (#2262).
+None
ctselect - CTA event selection
------------------------------
-Any .gz suffix is now correctly stripped from the event filenames when creating
-the output filenames (#2287).
-
-ctselect now informs about non-overlapping time interval; before the time interval
-was set to [0,0] which was disturbing and eventually incorrect (#2272).
-
-Use t-type parameters for "tmin" and "tmax" to all specifying times in various
-formats (#1864).
+None
ctskymap - CTA sky mapping
--------------------------
-Header keywords informing about the dates and times covered by the observations
-used to generate the skymap, as well as the event energy range, are now appended
-as header keywords to the FITS file (#2194).
-
-Add the RING background method, which is a classical analysis methods for IACT
-data analysis, for background subtraction (#2207).
+Add the "inmap", "iterations" and "threshold" parameters. The hidden
+"inmap" parameter allows specifying a sky map file containing a "COUNTS"
+and "ACCEPTANCE" extension which avoids looping over event lists in an
+observation container. The "iterations" parameter allows the iterative
+determination of an exclusion map from the resulting significance map for
+the RING background method. The exclusion map will flag all pixels with
+significances above the "threshold" value so that they will be excluded
+in the subsequent RING background determination. Since this will impact
+the pixel significances, the procedure needs to be repeated a couple of
+times to stabilise (#2311).
+
+Accelerate the IRF and RING background methods. For the RING background
+method the computation is either performed using a FFT or using a direct
+computation. The FFT computation is faster, but assumed Euclidian
+distances for the sky map grid, the direct computation is slower, but
+accurately computes the distances between sky map pixels (#2309).
cttsmap - Generation of Test Statistic map
------------------------------------------
-Add "statistic" parameter to allow setting of optimizer statistic (#2237).
+None
ctulimit - Compute upper limit
------------------------------
-The tool now also computes correct upper flux limits for diffuse map cubes
-for which the spectral component is a normalization of the spectral information
-available in the map cube (#2298).
-
-Add "statistic" parameter to allow setting of optimizer statistic (#2237).
+None
support - Support classes
-------------------------
-Change DATE_OBS, DATE_END, TIME_OBS and TIME_END FITS keywords into DATA-OBS,
-DATE-END, TIME-OBS and TIME-END.
-
-Add ctobservation.obs() method to set an observation container.
+None
cscaldb - Inspect calibration database
@@ -189,16 +171,7 @@ None
cslightcrv - Generate light curve
---------------------------------
-The script now also supports the handling of On/Off observations. For that
-purpose parameters were added that are needed for the execution of the
-csphagen script (#2273).
-
-The script now derives from the csobservation base class. A "statistic"
-parameter was added to allow setting the optimizer statistic for the light
-curve generation (#2237).
-
-Use t-type parameters for "tmin" and "tmax" to all specifying times in various
-formats (#1864).
+None
csmodelinfo - Return model information
@@ -218,8 +191,7 @@ None
csmodelsois - Puts selected sources in diffuse model cube
---------------------------------------------------------
-Add script that replaces selected models in a model container by a diffuse
-map cube (#2180).
+None
csobs2caldb - Generate calibration database entry for an IACT observation
@@ -234,47 +206,27 @@ None
csobsinfo - Return observation definition information
-----------------------------------------------------
-None
+Add ras() and decs() methods (#2313)
csobsselect - Select observations
---------------------------------
-Add tmin and tmax parameters to allow selection of observation by time (#2297).
+None
csphagen - Generate On/Off observations
---------------------------------------
-The CUSTOM background method was added to allow specifying On and Off regions
-from a region file or FITS WCS map. For multiple observations, a custom Off
-region per observation can be specified in the observation definition XML file
-using the OffRegions parameter (#2284).
-
-The script was added to allow the generation of PHA, ARF and RMF datasets for classical On/Off IACT analysis (#2208).
+None
csphasecrv - Generate phase curve
---------------------------------
-The script now also supports the handling of On/Off observations. For that
-purpose parameters were added that are needed for the execution of the
-csphagen script (#2273).
-
-The script now derives from the csobservation base class. A "statistic"
-parameter was added to allow setting the optimizer statistic for the light
-curve generation (#2237).
+None
cspull - Generate pull distributions
------------------------------------
-The script now handles On/Off observation. The parameters "onsrc" and "onrad"
-were added to allow specifying the On region test source and On region radius
-(#2260).
-
-The script now derives from the csobservation base class. A "statistic"
-parameter was added to allow setting the optimizer statistic for the pull
-distribution generation (#2237).
-
-Use t-type parameters for "tmin" and "tmax" to all specifying times in various
-formats (#1864).
+None
csresmap - Generate residual map
@@ -284,7 +236,7 @@ None
csresspec - Generate residual spectrum
--------------------------------------
-Add script to generate a residual spectrum (#2289)
+None
csroot2caldb - Generates calibration database from ROOT file
@@ -294,25 +246,12 @@ None
cssens - Computes sensitivity curve
-----------------------------------
-The script now derived from the csobservation base class. A "statistic"
-parameter was added to allow setting the optimizer statistic for the sensitivity
-computation (#2237).
+None
csspec - Generates Spectral Energy Distribution
-----------------------------------------------
-The script now derives from the csobservation base class. A "statistic"
-parameter was added to allow setting the optimizer statistic for the spectrum
-generation (#2237).
-
-The script now can also handle On/Off spectra for SED generation (#2231).
-
-The script now can also handle non-CTA data. For that purpose a "method"
-parameter has been introduced that allows to selected between slicing of the data
-(method="SLICE"), which is only supported for CTA observations, and replacing of
-the spectral model by a node function (method="NODES"). If method="AUTO" the tool
-will use "SLICE" for observations containing only CTA data and "NODES" for observations
-containing other data (#1927).
+None
cssrcdetect - Detects sources in sky map
@@ -322,18 +261,7 @@ None
cstsdist - Computes Test Statistic distribution
-----------------------------------------------
-The script now derives from the csobservation base class. A "statistic"
-parameter was added to allow setting the optimizer statistic for the Test
-Statistic distribution generation (#2237).
-
-The parameter interface was considerably revised, in particular the following
-parameters were removed: "ra", "dec", "rad", "coordsys", "proj", "emin",
-"emax", "enumbins", "tmin", "tmax", "npix" and "binsz". The "inobs" parameter
-is now always queried. The script now reads in an observation or an observation
-container and derives the TS distribution for a source based on that observation
-container. So far the script supports unbinned, binned and stacked observations.
-If an observation definition XML file is specified the script assumes that all
-observations are event lists (#1567).
+None
cstsmapmerge - Merge Test Statistic maps
@@ -343,14 +271,12 @@ None
cstsmapsplit - Prepare split Test Statistic computation
-------------------------------------------------------
-A "statistic" parameter was added to allow setting the optimizer statistic for
-the TS map generation (#2237).
+None
csviscube - Generate visibility cube
------------------------------------
-Use t-type parameters for "tmin" and "tmax" to all specifying times in various
-formats (#1864).
+None
csworkflow - Execute workflow
diff --git a/README.md b/README.md
index 86e6910..b81b128 100755
--- a/README.md
+++ b/README.md
@@ -1,6 +1,6 @@
ctools information
==================
-* Version: 1.6.0.dev1 (1 August 2017)
+* Version: 1.6.0.dev1
* Author: Juergen Knoedlseder (jurgen.knodlseder@irap.omp.eu)
* GammaLib dependency: 1.6.0.dev1
diff --git a/caldb/data/cta/prod2/bcf/North_0.5h/irf_file.fits.gz b/caldb/data/cta/prod2/bcf/North_0.5h/irf_file.fits.gz
index b4f12a7..3521fb1 100755
--- a/caldb/data/cta/prod2/bcf/North_0.5h/irf_file.fits.gz
+++ b/caldb/data/cta/prod2/bcf/North_0.5h/irf_file.fits.gz
diff --git a/caldb/data/cta/prod2/bcf/North_50h/irf_file.fits.gz b/caldb/data/cta/prod2/bcf/North_50h/irf_file.fits.gz
index 1650c4c..f15ad7a 100755
--- a/caldb/data/cta/prod2/bcf/North_50h/irf_file.fits.gz
+++ b/caldb/data/cta/prod2/bcf/North_50h/irf_file.fits.gz
diff --git a/caldb/data/cta/prod2/bcf/North_5h/irf_file.fits.gz b/caldb/data/cta/prod2/bcf/North_5h/irf_file.fits.gz
index 668517d..2ccf018 100755
--- a/caldb/data/cta/prod2/bcf/North_5h/irf_file.fits.gz
+++ b/caldb/data/cta/prod2/bcf/North_5h/irf_file.fits.gz
diff --git a/caldb/data/cta/prod2/bcf/South_0.5h/irf_file.fits.gz b/caldb/data/cta/prod2/bcf/South_0.5h/irf_file.fits.gz
index a68857d..552bb07 100755
--- a/caldb/data/cta/prod2/bcf/South_0.5h/irf_file.fits.gz
+++ b/caldb/data/cta/prod2/bcf/South_0.5h/irf_file.fits.gz
diff --git a/caldb/data/cta/prod2/bcf/South_50h/irf_file.fits.gz b/caldb/data/cta/prod2/bcf/South_50h/irf_file.fits.gz
index e9e0aa4..95b3d29 100755
--- a/caldb/data/cta/prod2/bcf/South_50h/irf_file.fits.gz
+++ b/caldb/data/cta/prod2/bcf/South_50h/irf_file.fits.gz
diff --git a/caldb/data/cta/prod2/bcf/South_5h/irf_file.fits.gz b/caldb/data/cta/prod2/bcf/South_5h/irf_file.fits.gz
index 8e7f24d..a063720 100755
--- a/caldb/data/cta/prod2/bcf/South_5h/irf_file.fits.gz
+++ b/caldb/data/cta/prod2/bcf/South_5h/irf_file.fits.gz
diff --git a/cscripts/csobsinfo.py b/cscripts/csobsinfo.py
index 36d3c39..4ae2d96 100755
--- a/cscripts/csobsinfo.py
+++ b/cscripts/csobsinfo.py
@@ -59,11 +59,11 @@ class csobsinfo(ctools.csobservation):
"""
Get parameters from parfile and setup the observation
"""
- # Get parameters
+ # Get parameters
if self.obs().size() == 0:
self._require_inobs('csobsinfo::get_parameters')
self.obs(self._get_observations(False))
-
+
# Initialise object position
self._obj_dir = gammalib.GSkyDir()
@@ -83,7 +83,7 @@ class csobsinfo(ctools.csobservation):
# Return
return
-
+
# Public methods
def run(self):
@@ -96,7 +96,7 @@ class csobsinfo(ctools.csobservation):
# Get parameters
self._get_parameters()
-
+
# Initialise arrays to store certain values for reuse
# Todo, think about using a python dictionary
self._offsets = []
@@ -107,7 +107,7 @@ class csobsinfo(ctools.csobservation):
self._ebounds = gammalib.GEbounds()
self._gti = gammalib.GGti()
obs_names = []
-
+
# Initialise output to be filled
ontime = 0.0
livetime = 0.0
@@ -115,12 +115,12 @@ class csobsinfo(ctools.csobservation):
n_eventbins = 0
n_obs_binned = 0
n_obs_unbinned = 0
-
+
# Write header
if self._logTerse():
self._log('\n')
self._log.header1(gammalib.number('Observation', self.obs().size()))
-
+
# Loop over observations
for obs in self.obs():
@@ -128,44 +128,44 @@ class csobsinfo(ctools.csobservation):
if not obs.classname() == 'GCTAObservation':
self._log('Skipping '+obs.instrument()+' observation\n')
continue
-
+
# Use observed object as observation name if name is not given
obs_name = obs.name()
if obs_name == '':
obs_name = obs.object()
-
+
# Logging
if self._logExplicit():
self._log.header2(obs_name)
-
+
# Retrieve observation name
obs_names.append(obs_name)
-
+
# Retrieve energy boundaries
obs_bounds = obs.events().ebounds()
-
+
# Retrieve time interval
obs_gti = obs.events().gti()
-
+
# Compute mean time and dead time fraction in percent
deadfrac = (1.0-obs.deadc())*100.0
-
+
# Retrieve pointing and store Ra,Dec
- pnt_dir = obs.pointing().dir()
+ pnt_dir = obs.pointing().dir()
self._pnt_ra.append(pnt_dir.ra_deg())
- self._pnt_dec.append(pnt_dir.dec_deg())
-
+ self._pnt_dec.append(pnt_dir.dec_deg())
+
# If avaliable append energy boundaries
if obs_bounds.size() > 0 :
self._ebounds.append(obs_bounds.emin(),obs_bounds.emax())
-
+
# Append time interval
self._gti.append(obs_gti.tstart(), obs_gti.tstop())
-
+
# Increment global livetime and ontime
ontime += obs.ontime()
livetime += obs.livetime()
-
+
# Bookkeeping
if obs.eventtype() == 'CountsCube':
n_eventbins += obs.events().size()
@@ -190,10 +190,10 @@ class csobsinfo(ctools.csobservation):
if self._compute_offset:
offset = pnt_dir.dist_deg(self._obj_dir)
self._offsets.append(offset)
-
+
# Optionally log details
if self._logExplicit():
-
+
# Log the observation energy range (if available)
self._log.parformat('Energy range')
if obs_bounds.size() == 0:
@@ -203,7 +203,7 @@ class csobsinfo(ctools.csobservation):
self._log(' - ')
self._log(str(obs_bounds.emax()))
self._log('\n')
-
+
# Log observation time interval
self._log.parformat('Time range (MJD)')
if obs_gti.size() == 0:
@@ -213,7 +213,7 @@ class csobsinfo(ctools.csobservation):
self._log(' - ')
self._log(str(obs_gti.tstop().mjd()))
self._log('\n')
-
+
# Log observation information
self._log_value(gammalib.EXPLICIT, 'Ontime', '%.3f s' %
obs.ontime())
@@ -233,7 +233,7 @@ class csobsinfo(ctools.csobservation):
zenith)
self._log_value(gammalib.EXPLICIT, 'Azimuth angle', '%.2f deg' %
azimuth)
-
+
# Write summary header
self._log_header1(gammalib.NORMAL, 'Summary')
@@ -258,7 +258,7 @@ class csobsinfo(ctools.csobservation):
mean_azimuth = '%.2f deg' % (sum(self._azimuths) / len(self._azimuths))
else:
mean_azimuth = 'Unknown'
-
+
# Log mean offset, azimuth and zenith angle
self._log_header3(gammalib.NORMAL, 'Pointings')
self._log_value(gammalib.NORMAL, 'Mean offset angle', mean_offset)
@@ -323,48 +323,60 @@ class csobsinfo(ctools.csobservation):
# Check if DS9 file is valid
if ds9file.url() != 'NONE':
-
+
# Write header
self._log_header1(gammalib.TERSE, 'Save pointings in DS9 file')
# Log filename
self._log_value(gammalib.NORMAL, 'DS9 filename', ds9file.url())
-
- # Open file
+
+ # Open file
f = open(ds9file.url(),'w')
-
+
# Write coordinate system
f.write('fk5\n')
-
+
# Loop over pointings
for i in range(len(self._pnt_ra)):
-
+
# Create string
line = 'point('
line += str(self._pnt_ra[i])+','+str(self._pnt_dec[i])+')'
line += ' # point=cross 20 width=3\n'
-
+
# Write to file
f.write(line)
-
+
# Close file
f.close()
-
+
# Return
- return
-
+ return
+
def zeniths(self):
"""
Return zenith angles
"""
return self._zeniths
-
+
def azimuths(self):
"""
Return azimuth angles
"""
return self._azimuths
-
+
+ def ras(self):
+ """
+ Return pointings right ascension
+ """
+ return self._pnt_ra
+
+ def decs(self):
+ """
+ Return pointings declination
+ """
+ return self._pnt_dec
+
def offsets(self):
"""
Return offset angles
@@ -376,13 +388,13 @@ class csobsinfo(ctools.csobservation):
Return energy boundaries
"""
return self._ebounds
-
+
def gti(self):
"""
Return good time intervals
"""
return self._gti
-
+
# ======================== #
# Main routine entry point #
@@ -391,6 +403,6 @@ if __name__ == '__main__':
# Create instance of application
app = csobsinfo(sys.argv)
-
+
# Execute application
app.execute()
diff --git a/cscripts/csphagen.py b/cscripts/csphagen.py
index 48cdb49..f90e78a 100755
--- a/cscripts/csphagen.py
+++ b/cscripts/csphagen.py
@@ -123,7 +123,8 @@ class csphagen(ctools.csobservation):
# Off region region the parameter 'bkgregfile' for all CTA observations
# without Off region
for obs in self.obs():
- if obs.instrument() == 'CTA':
+ #if obs.instrument() == 'CTA':
+ if obs.classname() == 'GCTAObservation':
if obs.off_regions().is_empty():
filename = self['bkgregfile'].filename()
regions = gammalib.GSkyRegions(filename)
@@ -263,8 +264,8 @@ class csphagen(ctools.csobservation):
"""
Set models in observation container
- The method replaces all "CTA" background models by "CTAOnOff"
- background models.
+ The method replaces all "CTA", "HESS", "VERITAS", and "MAGIC"
+ background models by "CTAOnOff" background models.
Parameters
----------
@@ -279,11 +280,15 @@ class csphagen(ctools.csobservation):
# Initialise model container
models = gammalib.GModels()
- # Loop over all models and replace CTA background model by CTAOnOff
- # background model
+ # Loop over all models and replace CTA/HESS/VERITAS/MAGIC background
+ # models by CTAOnOff background model
for model in self.obs().models():
- if 'GCTA' in model.classname() and 'CTA' in model.instruments():
- model.instruments('CTAOnOff')
+ if 'GCTA' in model.classname():
+ if 'CTA' in model.instruments() or \
+ 'HESS' in model.instruments() or \
+ 'VERITAS' in model.instruments() or \
+ 'MAGIC' in model.instruments():
+ model.instruments('CTAOnOff')
models.append(model)
# Append model to observation container
@@ -399,7 +404,8 @@ class csphagen(ctools.csobservation):
for obs in self.obs():
# Skip non CTA observations
- if obs.instrument() != 'CTA':
+ #if obs.instrument() != 'CTA':
+ if obs.classname() != 'GCTAObservation':
self._log_string(gammalib.NORMAL, 'Skip %s observation "%s"' % \
(obs.instrument(), obs.id()))
continue
diff --git a/cscripts/csresspec.py b/cscripts/csresspec.py
index 0c8ac79..ad95ef4 100755
--- a/cscripts/csresspec.py
+++ b/cscripts/csresspec.py
@@ -241,14 +241,14 @@ class csresspec(ctools.csobservation):
# Bin events
cntcube = ctools.ctbin(obs)
- cntcube['xref'] = ra
- cntcube['yref'] = dec
- cntcube['binsz'] = 0.02
- cntcube['nxpix'] = npix
- cntcube['nypix'] = npix
- cntcube['proj'] = 'TAN'
+ cntcube['xref'] = ra
+ cntcube['yref'] = dec
+ cntcube['binsz'] = 0.02
+ cntcube['nxpix'] = npix
+ cntcube['nypix'] = npix
+ cntcube['proj'] = 'TAN'
cntcube['coordsys'] = 'CEL'
- cntcube['ebinalg'] = self['ebinalg'].string()
+ cntcube['ebinalg'] = self['ebinalg'].string()
if self['ebinalg'].string() == 'FILE':
cntcube['ebinfile'] = self['ebinfile'].filename()
else:
@@ -356,7 +356,7 @@ class csresspec(ctools.csobservation):
def _cube_to_spectrum(self, cube, evlist_info):
"""
Derive from event cube a count spectrum
-
+
If data come from event list use only the ROI and energy range of
the original data. Apply user defined mask if requested.
@@ -550,29 +550,30 @@ class csresspec(ctools.csobservation):
# if 3D observation
if obs[0].classname() == 'GCTAObservation':
- # Prepare Observations
-
# If already binned set the evlist_info dictionary to have
# attribute was_list False
if obs[0].eventtype() == 'CountsCube':
evlist_info = {'was_list': False}
- # Otherwise bin now
+
+ # ... otherwise bin now
else:
- # we remember if we binned an event list
- # so that we can mask only the ROI for residual calculation
+ # we remember if we binned an event list so that we can
+ # mask only the ROI for residual calculation
msg = 'Setting up binned observation'
self._log_string(gammalib.NORMAL, msg)
obs, evlist_info = self._bin_evlist(obs)
- # Calculate Model and residuals
- # If model cube is provided load it
+ # Calculate Model and residuals. If model cube is provided load
+ # it
if self._use_maps:
modcube = gammalib.GCTAEventCube(self['inmodel'].filename())
- # Otherwise calculate it now
+
+ # ... otherwise calculate it now
else:
msg = 'Computing model cube'
self._log_string(gammalib.NORMAL, msg)
modelcube = ctools.ctmodel(obs)
+ modelcube['edisp'] = self['edisp'].boolean()
modelcube.run()
modcube = modelcube.cube().copy()
diff --git a/cscripts/csroot2caldb.py b/cscripts/csroot2caldb.py
index 477f4ad..0123254 100755
--- a/cscripts/csroot2caldb.py
+++ b/cscripts/csroot2caldb.py
@@ -2,7 +2,7 @@
# ==========================================================================
# Generate IRFs in CALDB format from a ROOT offaxis performance file
#
-# Copyright (C) 2016-2017 Juergen Knoedlseder
+# Copyright (C) 2016-2018 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
@@ -618,7 +618,7 @@ class csroot2caldb(ctools.cscript):
ds['HDU_EA'].card('HI_THRES', emin, '[TeV] High energy threshold')
# Create "EFFAREA" data column
- self._make_2D(etrue, ds['HDU_EA'], 'EFFAREA', 'm2')
+ self._make_2D(etrue, ds['HDU_EA'], 'EFFAREA', 'm**2')
# Return
return
@@ -745,7 +745,6 @@ class csroot2caldb(ctools.cscript):
'Generate Gaussian point spread function extension')
# Get relevant ROOT histograms
- #r68 = tfile.Get('AngRes_offaxis')
r68, _ = self._root2psf_get_r68_and_r80(tfile)
# Extract number of bins in histogram
@@ -826,8 +825,6 @@ class csroot2caldb(ctools.cscript):
'Generate King point spread function extension')
# Get relevant ROOT histograms
- #r68 = tfile.Get('AngRes_offaxis')
- #r80 = tfile.Get('AngRes80_offaxis')
r68, r80 = self._root2psf_get_r68_and_r80(tfile)
# Extract number of bins in histogram
@@ -1034,7 +1031,6 @@ class csroot2caldb(ctools.cscript):
# Get relevant ROOT histograms
matrix = self._root2edisp_get_migra(tfile)
- #matrix = tfile.Get('EestOverEtrue_offaxis')
# Write boundary keywords
self._set_cif_keywords(ds['HDU_EDISP'], irf['EDISP_NAME'],
@@ -1087,7 +1083,7 @@ class csroot2caldb(ctools.cscript):
irf['BGD_BOUNDS'], irf['BGD_DESC'], irf)
# Create "BGD" data column
- self._make_3D(array, ds['HDU_BGD'], 'BGD', '1/s/MeV/sr')
+ self._make_3D(array, ds['HDU_BGD'], 'BGD', '1/(MeV s sr)')
# Return
return
@@ -1327,7 +1323,7 @@ class csroot2caldb(ctools.cscript):
# Log power law replacement
self._log_value(gammalib.NORMAL,
'Power law replacement',
- '%e cts/s/deg2 (E=%8.4f TeV (%d) '
+ '%e cts/(s deg**2) (E=%8.4f TeV (%d) '
'Off=%.2f deg (%d)' %
(value,engs[ieng],ieng,offset,ioff))
@@ -1479,7 +1475,7 @@ class csroot2caldb(ctools.cscript):
# Log background infill
self._log_value(gammalib.NORMAL,
'Background infill',
- '%e cts/s/deg2 (E=%8.4f TeV (%d) '
+ '%e cts/(s deg**2) (E=%8.4f TeV (%d) '
'Off=%.2f deg (%d)' %
(value,engs[ieng],ieng,offset,ioff))
diff --git a/doc/source/admin/install_devel.rst b/doc/source/admin/install_devel.rst
index 29389d4..87c86ca 100644
--- a/doc/source/admin/install_devel.rst
+++ b/doc/source/admin/install_devel.rst
@@ -7,7 +7,7 @@ The current development version of the code can be downloaded as source
tarballs, Mac OS X binaries, or directly from the `Git `_
repository.
-The current ctools development release is ``ctools-1.5.0.dev1``.
+The current ctools development release is ``ctools-1.6.0.dev1``.
Source tarballs
@@ -15,8 +15,8 @@ Source tarballs
Download the source tarballs from the following links
-- `ctools `_
-- `GammaLib `_
+- `ctools `_
+- `GammaLib `_
and follow the instructions on :ref:`sec_install_source`.
@@ -26,7 +26,7 @@ Binary packages
Download the installer image from the following link
-- `Mac OS X (10.7+) `_
+- `Mac OS X (10.7+) `_
and follow the instructions on :ref:`sec_install_binary`.
diff --git a/doc/source/admin/release_history/1.5.1.rst b/doc/source/admin/release_history/1.5.1.rst
new file mode 100644
index 0000000..da8dcaf
--- /dev/null
+++ b/doc/source/admin/release_history/1.5.1.rst
@@ -0,0 +1,20 @@
+.. _1.5.1:
+
+ctools 1.5.1 release
+====================
+
+Overview
+--------
+
+ctools 1.5.1 is a bug fix release for ctools release 1.5.
+
+The following bugs were fixed:
+
+* [`2341 `_] -
+ Take energy dispersion in :ref:`csresspec` modelling into account
+
+
+The following improvements were made:
+
+* [`2307 `_] -
+ Enable using :ref:`csphagen` for HESS, VERITAS and MAGIC observations
diff --git a/doc/source/admin/release_history/1.6.rst b/doc/source/admin/release_history/1.6.rst
new file mode 100644
index 0000000..5d72ffa
--- /dev/null
+++ b/doc/source/admin/release_history/1.6.rst
@@ -0,0 +1,42 @@
+.. _1.6:
+
+ctools 1.6 release
+==================
+
+Overview
+--------
+
+ctools 1.6 is a major release that adds significant functionality.
+
+In particular, this release provides:
+
+* ...
+
+
+Bug fixes
+---------
+
+* ...
+
+
+Improvements
+------------
+
+* [`2313 `_] -
+ Add `ras()` and `decs()` methods to :ref:`csobsinfo`
+* [`2311 `_] -
+ Allow specifying of input skymap to :ref:`ctskymap` tool
+* [`2309 `_] -
+ Accelerate `IRF` and `RING` background methods in :ref:`ctskymap`
+
+
+New features
+------------
+
+* ...
+
+
+Clean up
+--------
+
+* ...
diff --git a/doc/source/admin/release_history/index.rst b/doc/source/admin/release_history/index.rst
index b969bd1..4e20099 100755
--- a/doc/source/admin/release_history/index.rst
+++ b/doc/source/admin/release_history/index.rst
@@ -6,6 +6,8 @@ Release History
.. toctree::
:maxdepth: 1
+ ctools 1.6 <1.6>
+ ctools 1.5.1 <1.5.1>
ctools 1.5 <1.5>
ctools 1.4.3 <1.4.3>
ctools 1.4.2 <1.4.2>
diff --git a/doc/source/users/reference_manual/ctskymap.rst b/doc/source/users/reference_manual/ctskymap.rst
index 80627f7..6af9c7e 100755
--- a/doc/source/users/reference_manual/ctskymap.rst
+++ b/doc/source/users/reference_manual/ctskymap.rst
@@ -41,6 +41,11 @@ General parameters
Instrument response function (only required for IRF background subtraction
if no response information is provided by ``inobs``).
+``(inmap = NONE) [file]``
+ Input sky map file containing a ``COUNTS`` and an ``ACCEPTANCE`` extension.
+ Such extensions are produce by a previous run of the :ref:`ctskymap` tool with
+ the ``IRF`` or ``RING`` background method.
+
``outmap [file]``
Output sky map file.
@@ -82,17 +87,37 @@ Background subtraction configuration parameters
Background subtraction method.
``roiradius [real]``
- Source region radius for estimating on-counts (in degrees).
+ Source region radius for ``RING`` subtraction (in degrees).
``inradius [real]``
- Inner ring radius for RING subtraction (in degrees).
+ Inner background ring radius for ``RING`` subtraction (in degrees).
``outradius [real]``
- Outer ring radius for RING subtraction (in degrees).
+ Outer background ring radius for ``RING`` subtraction (in degrees).
+
+``iterations [integer]``
+ Number of iterations for the automatic computation of exclusion regions.
+ Sky map pixels with detection significance above the value specified by the
+ ``threshold`` parameter will be defined as exclusion regions, and the ``RING``
+ background will be recomputed based on these exclusion regions. Since this
+ will likely change the significance of the pixels, the procedure can be
+ iteratively repeated, and the ``iterations`` parameter specifies how often
+ the procedure should be repeated.
+
+``threshold [real]``
+ Significance threshold above which pixels should be excluded for the
+ background computation in the ``RING`` background method.
``(inexclusion = NONE) [file]``
Exclusion region file as either a FITS map or DS9 region file.
+``(usefft = yes) [boolean]``
+ Specifies whether a FFT should be used for ring computation. If set to
+ ``no`` the ring background will be computed directly. Direct computation is
+ more accurate since the precise angular distances between pixels are taken
+ into account, but FFT is faster, but less accurate since it assumes
+ Euclidean distances between pixels.
+
Standard parameters
-------------------
diff --git a/pyext/ctobservation.i b/pyext/ctobservation.i
index 9270a21..38f221f 100755
--- a/pyext/ctobservation.i
+++ b/pyext/ctobservation.i
@@ -94,6 +94,7 @@ public:
// Make methods private in Python by prepending an underscore
%rename(_first_unbinned_observation) first_unbinned_observation;
%rename(_next_unbinned_observation) next_unbinned_observation;
+ %rename(_read_ogip_keywords) read_ogip_keywords;
%rename(_write_ogip_keywords) write_ogip_keywords;
%rename(_set_obs_statistic) set_obs_statistic;
%rename(_save_events_fits) save_events_fits;
diff --git a/src/ctskymap/ctskymap.cpp b/src/ctskymap/ctskymap.cpp
index e3bdca7..a75fe29 100755
--- a/src/ctskymap/ctskymap.cpp
+++ b/src/ctskymap/ctskymap.cpp
@@ -1,7 +1,7 @@
/***************************************************************************
* ctskymap - Sky mapping tool *
* ----------------------------------------------------------------------- *
- * copyright (C) 2011-2017 by Juergen Knoedlseder *
+ * copyright (C) 2011-2018 by Juergen Knoedlseder *
* ----------------------------------------------------------------------- *
* *
* This program is free software: you can redistribute it and/or modify *
@@ -30,14 +30,14 @@
#endif
#include
#include "ctskymap.hpp"
-#include "GSkyRegions.hpp"
/* __ Method name definitions ____________________________________________ */
#define G_GET_PARAMETERS "ctskymap::get_parameter()"
#define G_MAP_EVENTS "ctskymap::map_events(GCTAObservation*)"
#define G_MAP_BACKGROUND_IRF "ctskymap::map_background_irf(GCTAObservation*)"
-#define G_MAP_BACKGROUND_RING "ctskymap::map_background_ring("\
- "GCTAObservation*)"
+#define G_RING_BOUNDING_BOX "ctskymap::ring_bounding_box(int&, int&, int&, "\
+ "int&, int&)"
+#define G_RING_KERNEL "ctskymap::ring_kernel(double&, double&)"
/* __ Debug definitions __________________________________________________ */
@@ -222,37 +222,18 @@ void ctskymap::run(void)
// Get parameters
get_parameters();
- // Write input observation container into logger
- log_observations(NORMAL, m_obs, "Input observation");
-
- // Write header into logger
- log_header1(TERSE, gammalib::number("Map observation", m_obs.size()));
-
- // Loop over all unbinned CTA observations in the container
- for (GCTAObservation* obs = first_unbinned_observation(); obs != NULL;
- obs = next_unbinned_observation()) {
-
- // Map events into sky map
- map_events(obs);
+ // Setup maps
+ setup_maps();
- // Compute background sky map
- map_background(obs);
-
- // Dispose events to free memory
- obs->dispose_events();
-
- } // endfor: looped over observations
-
- // If background subtraction is selected then compute significance map
- if (m_bkgsubtract != "NONE") {
-
- // Compute the significance map
- map_significance();
-
- // Subtract background map from counts map
- m_skymap -= m_bkgmap;
+ // If no input map exists then fill counts and acceptance maps from
+ // observation container
+ if (!m_has_inmap) {
+ fill_maps();
}
+ // Compute maps
+ compute_maps();
+
// Optionally publish sky map
if (m_publish) {
publish();
@@ -266,8 +247,8 @@ void ctskymap::run(void)
/***********************************************************************//**
* @brief Save sky map
*
- * Saves the sky map into a FITS file. The FITS file name is specified by the
- * @p outname parameter.
+ * Saves the sky map into a FITS file. The FITS file name is specified by
+ * the @p outname parameter.
***************************************************************************/
void ctskymap::save(void)
{
@@ -284,39 +265,31 @@ void ctskymap::save(void)
GFits fits;
// Write sky map into FITS file
- GFitsHDU* hdu = m_skymap.write(fits);
-
- // Write keywords into sky map extension
- write_ogip_keywords(hdu);
- write_hdu_keywords(hdu);
+ write_map(fits, m_skymap, "SKYMAP");
- // If background subtraction is requested then write background map
- // and significance map to FITS file
+ // If background subtraction is requested then write the background map,
+ // the significance map, the counts map and the acceptance map to the
+ // FITS file
if (m_bkgsubtract != "NONE") {
-
- // Write background map into FITS file
- hdu = m_bkgmap.write(fits);
- // Set background map extension name
- if (hdu != NULL) {
- hdu->extname("BACKGROUND");
- }
+ // Write background map into FITS file
+ write_map(fits, m_bkgmap, "BACKGROUND");
- // Write keywords into background extension
- write_ogip_keywords(hdu);
- write_hdu_keywords(hdu);
// Write significance map into FITS file
- hdu = m_sigmap.write(fits);
+ write_map(fits, m_sigmap, "SIGNIFICANCE");
- // Set significance map extension name
- if (hdu != NULL) {
- hdu->extname("SIGNIFICANCE");
- }
+ // Write counts map into FITS file
+ write_map(fits, m_counts, "COUNTS");
+
+ // Write acceptance map into FITS file
+ write_map(fits, m_acceptance, "ACCEPTANCE");
- // Write keywords into significance extension
- write_ogip_keywords(hdu);
- write_hdu_keywords(hdu);
+ // If RING background subtraction is requested then write also
+ // the exclusion map to the FITS file
+ if (m_bkgsubtract == "RING") {
+ write_map(fits, m_exclmap, "EXCLUSION");
+ }
} // endif: background subtraction was requested
@@ -375,25 +348,36 @@ void ctskymap::publish(const std::string& name)
***************************************************************************/
void ctskymap::init_members(void)
{
- // Initialise members
- m_skymap.clear();
- m_bkgmap.clear();
- m_sigmap.clear();
- m_exclmap.clear();
- m_alphamap.clear();
- m_onmap.clear();
+ // Initialise User parameters
+ m_outmap.clear();
+ m_inexclusion.clear();
m_emin = 0.0;
m_emax = 0.0;
m_bkgsubtract = "NONE";
m_roiradius = 0.0;
m_inradius = 0.0;
m_outradius = 0.0;
+ m_iterations = 0;
+ m_threshold = 5.0;
+ m_usefft = true;
m_publish = false;
m_chatter = static_cast(2);
+ // Initialise members
+ m_has_inmap = false;
+ m_skymap.clear();
+ m_bkgmap.clear();
+ m_sigmap.clear();
+ m_counts.clear();
+ m_acceptance.clear();
+ m_exclmap.clear();
+
// Initialise cache
m_solidangle.clear();
m_dirs.clear();
+ m_cos_roiradius = 1.0;
+ m_cos_inradius = 1.0;
+ m_cos_outradius = 1.0;
// Return
return;
@@ -407,25 +391,36 @@ void ctskymap::init_members(void)
***************************************************************************/
void ctskymap::copy_members(const ctskymap& app)
{
- // Copy attributes
- m_skymap = app.m_skymap;
- m_bkgmap = app.m_bkgmap;
- m_sigmap = app.m_sigmap;
- m_exclmap = app.m_exclmap;
- m_alphamap = app.m_alphamap;
- m_onmap = app.m_onmap;
+ // Copy User parameters
+ m_outmap = app.m_outmap;
+ m_inexclusion = app.m_inexclusion;
m_emin = app.m_emin;
m_emax = app.m_emax;
m_bkgsubtract = app.m_bkgsubtract;
m_roiradius = app.m_roiradius;
m_inradius = app.m_inradius;
m_outradius = app.m_outradius;
+ m_iterations = app.m_iterations;
+ m_threshold = app.m_threshold;
+ m_usefft = app.m_usefft;
m_publish = app.m_publish;
m_chatter = app.m_chatter;
+ // Copy members
+ m_has_inmap = app.m_has_inmap;
+ m_skymap = app.m_skymap;
+ m_bkgmap = app.m_bkgmap;
+ m_sigmap = app.m_sigmap;
+ m_counts = app.m_counts;
+ m_acceptance = app.m_acceptance;
+ m_exclmap = app.m_exclmap;
+
// Copy cache
- m_solidangle = app.m_solidangle;
- m_dirs = app.m_dirs;
+ m_solidangle = app.m_solidangle;
+ m_dirs = app.m_dirs;
+ m_cos_roiradius = app.m_cos_roiradius;
+ m_cos_inradius = app.m_cos_inradius;
+ m_cos_outradius = app.m_cos_outradius;
// Return
return;
@@ -453,41 +448,92 @@ void ctskymap::free_members(void)
***************************************************************************/
void ctskymap::get_parameters(void)
{
- // Setup observations from "inobs" parameter. Do not request response
- // information and do not accept counts cubes.
- setup_observations(m_obs, false, true, false);
+ // So far we have no valid input map
+ m_has_inmap = false;
+
+ // If an input map was specified then extract try loading the COUNTS
+ // and ACCEPTANCE extensions from that file
+ GFilename inmap = (*this)["inmap"].filename();
+ if (inmap != "NONE") {
+ GFits fits(inmap);
+ if (fits.contains("COUNTS") && fits.contains("ACCEPTANCE")) {
+
+ // Get HDUs
+ GFitsHDU &counts = *fits["COUNTS"];
+ GFitsHDU &acceptance = *fits["ACCEPTANCE"];
+
+ // Read counts and acceptance map
+ m_counts.read(counts);
+ m_acceptance.read(acceptance);
- // Create sky map based on task parameters
- m_skymap = create_map(m_obs);
+ // Extract OGIP parameters
+ read_ogip_keywords(&counts);
- // Get further parameters
- m_emin = (*this)["emin"].real();
- m_emax = (*this)["emax"].real();
+ // Extract energy range
+ m_emin = (counts.has_card("E_MIN")) ? counts.real("E_MIN") : 0.0;
+ m_emax = (counts.has_card("E_MAX")) ? counts.real("E_MAX") : 0.0;
+
+ // Signal availability of sky map
+ m_has_inmap = true;
+ }
+ else {
+ std::string msg = "Input sky map \""+inmap.url()+"\" does not "
+ "contain \"COUNTS\" and \"ACCEPTANCE\" "
+ "extensions.";
+ throw GException::invalid_value(G_GET_PARAMETERS, msg);
+ }
+ fits.close();
+ }
+
+ // If no valid input map exists then request observations
+ if (!m_has_inmap) {
+
+ // Setup observations from "inobs" parameter. Do not request response
+ // information and do not accept counts cubes.
+ setup_observations(m_obs, false, true, false);
+
+ // Create counts and acceptance maps from the User parameters
+ m_counts = create_map(m_obs);
+ m_acceptance = m_counts;
+
+ // Get further parameters
+ m_emin = (*this)["emin"].real();
+ m_emax = (*this)["emax"].real();
+
+ } // endif: no valid input map existed
+
+ // Get background method
m_bkgsubtract = (*this)["bkgsubtract"].string();
// Get RING background parameters
if (m_bkgsubtract == "RING") {
// Get parameters
- m_roiradius = (*this)["roiradius"].real();
- m_inradius = (*this)["inradius"].real();
- m_outradius = (*this)["outradius"].real();
+ m_roiradius = (*this)["roiradius"].real();
+ m_inradius = (*this)["inradius"].real();
+ m_outradius = (*this)["outradius"].real();
+ m_iterations = (*this)["iterations"].integer();
+ if (m_iterations > 0) {
+ m_threshold = (*this)["threshold"].real();
+ }
+ m_inexclusion = (*this)["inexclusion"].filename();
+ m_usefft = (*this)["usefft"].boolean();
// Make sure that (roiradius < inradius < outradius)
if (m_roiradius > m_inradius) {
- std::string msg("'roiradius' must be smaller than 'inradius'");
+ std::string msg("\"roiradius\" must be smaller than \"inradius\"");
throw GException::invalid_value(G_GET_PARAMETERS, msg);
}
else if (m_inradius > m_outradius) {
- std::string msg("'inradius' must be smaller than 'outradius'");
+ std::string msg("\"inradius\" must be smaller than \"outradius\"");
throw GException::invalid_value(G_GET_PARAMETERS, msg);
}
} // endif: read RING background parameters
- // If IRF background subtraction is requested then make sure that the
- // CTA observations in the observation container have response information
- if (m_bkgsubtract != "NONE") {
+ // If background subtraction is requested then make sure that the CTA
+ // observations in the observation container have response information
+ if (m_bkgsubtract != "NONE" && !m_has_inmap) {
set_response(m_obs);
}
@@ -497,39 +543,209 @@ void ctskymap::get_parameters(void)
// Read ahead parameters
if (read_ahead()) {
- m_outmap = (*this)["outmap"].filename();
+ m_outmap = (*this)["outmap"].filename();
}
+ // Write parameters into logger
+ log_parameters(TERSE);
+
+ // Return
+ return;
+}
+
+
+/***********************************************************************//**
+ * @brief Setup maps
+ *
+ * Setup maps for sky map generation.
+ ***************************************************************************/
+void ctskymap::setup_maps(void)
+{
+ // Clear cache
+ m_solidangle.clear();
+ m_dirs.clear();
+
// Create background map and significance map if background subtraction
// is requested
if (m_bkgsubtract != "NONE") {
- // Create backgrond and significance maps
- m_bkgmap = create_map(m_obs);
- m_sigmap = create_map(m_obs);
-
// Setup the exclusions map
- map_exclusions((*this)["inexclusion"].filename());
+ setup_exclusion_map(m_inexclusion);
// Cache the pixel solid angles and sky directions
- m_solidangle.reserve(m_bkgmap.npix());
- m_dirs.reserve(m_bkgmap.npix());
- for (int i = 0; i < m_bkgmap.npix(); ++i) {
- m_solidangle.push_back(m_bkgmap.solidangle(i));
- m_dirs.push_back(m_bkgmap.inx2dir(i));
+ m_solidangle.reserve(m_counts.npix());
+ m_dirs.reserve(m_counts.npix());
+ for (int i = 0; i < m_counts.npix(); ++i) {
+ m_solidangle.push_back(m_counts.solidangle(i));
+ m_dirs.push_back(m_counts.inx2dir(i));
}
- // If doing a ring background subtraction, generate an alpha map and
- // make sure to cache the counts for each observation
- if (m_bkgsubtract == "RING") {
- m_alphamap = create_map(m_obs);
- m_onmap = create_map(m_obs);
+ } // endif: background subtraction selected
+
+ // Compute cosine of ring radii for RING background method
+ if (m_bkgsubtract == "RING") {
+ m_cos_roiradius = std::cos(m_roiradius * gammalib::deg2rad);
+ m_cos_inradius = std::cos(m_inradius * gammalib::deg2rad);
+ m_cos_outradius = std::cos(m_outradius * gammalib::deg2rad);
+ }
+
+ // Return
+ return;
+}
+
+
+/***********************************************************************//**
+ * @brief Generates map of pixel exclusions
+ *
+ * @param[in] filename Exclusion file name.
+ *
+ * Generates a sky map of the pixels that are to be excluded from the
+ * background estimation. Pixels with values different from 0 will be
+ * excluded.
+ ***************************************************************************/
+void ctskymap::setup_exclusion_map(const GFilename& filename)
+{
+ // Create exlusion map
+ m_exclmap = m_counts;
+
+ // Set all pixels to 0 (no pixel excluded)
+ m_exclmap = 0.0;
+
+ // Make sure the exclusions filename is valid
+ if (is_valid_filename(filename)) {
+
+ // Fill the exclusions based on the regions supplied
+ if (filename.is_fits()) {
+ setup_exclusion_map_fits(filename);
+ }
+ else {
+ setup_exclusion_map_region(filename);
}
- } // endif: background subtraction selected
+ } // endif: filename was valid
- // Write parameters into logger
- log_parameters(TERSE);
+ // Return
+ return;
+}
+
+
+/***********************************************************************//**
+ * @brief Fills exclusions map from FITS image
+ *
+ * @param[in] filename FITS image file name.
+ *
+ * Sets all exclusion map pixels to 1 that correspond to non-zero pixels in
+ * the exclusion sky map FITS file.
+ ***************************************************************************/
+void ctskymap::setup_exclusion_map_fits(const GFilename& filename)
+{
+ // Load the fits image
+ GSkyMap inmap(filename);
+
+ // Loop through the individual pixels in the exclusion map
+ for (int i = 0; i < m_exclmap.npix(); ++i) {
+
+ // Get the pixel direction
+ GSkyDir dir = m_exclmap.pix2dir(i);
+
+ // Check this sky position in the fits map
+ if (inmap.contains(dir) && (inmap(inmap.dir2inx(dir)) != 0.0)) {
+
+ // Set the pixel to 1
+ m_exclmap(i) = 1.0;
+
+ } // endif: pixel,region overlap check
+
+ } // endfor: looped over exclusion map pixels
+
+ // Return
+ return;
+}
+
+
+/***********************************************************************//**
+ * @brief Fills exclusions map from DS9 region file
+ *
+ * @param[in] filename DS9 region file name.
+ *
+ * Sets all exclusion map pixels to 1 that are contained in any of the DS9
+ * regions.
+ ***************************************************************************/
+void ctskymap::setup_exclusion_map_region(const GFilename& filename)
+{
+ // Load the exclusion regions
+ GSkyRegions regions(filename);
+
+ // Loop through the individual pixels in the exclusion map
+ for (int i = 0; i < m_exclmap.npix(); ++i) {
+
+ // Get the pixel position
+ GSkyDir dir = m_exclmap.pix2dir(i);
+
+ // If pixel position overlaps with the regions
+ if (regions.contains(dir)) {
+ m_exclmap(i) = 1.0;
+ }
+
+ } // endfor: looped over exclusion map pixels
+
+ // Return
+ return;
+}
+
+
+/***********************************************************************//**
+ * @brief Fill maps from observation container
+ *
+ * Fill counts and optionally acceptance maps from data in the observation
+ * container.
+ ***************************************************************************/
+void ctskymap::fill_maps(void)
+{
+ // Write input observation container into logger
+ log_observations(NORMAL, m_obs, "Input observation");
+
+ // Write header into logger
+ log_header1(TERSE, gammalib::number("Find unbinned observation",
+ m_obs.size()));
+
+ // Find all unbinned CTA observations in m_obs
+ std::vector obs_list(0);
+ for (GCTAObservation* obs = first_unbinned_observation(); obs != NULL;
+ obs = next_unbinned_observation()) {
+
+ // Push observation into list
+ obs_list.push_back(obs);
+
+ // Write message
+ std::string msg = " Including unbinned "+obs->instrument()+
+ " observation";
+ log_string(NORMAL, msg);
+ }
+
+ // Write header into logger
+ log_header1(TERSE, gammalib::number("Fill map from observation",
+ m_obs.size()));
+
+ // Loop over all unbinned CTA observations in the container
+ #pragma omp parallel for
+ for (int i = 0; i < obs_list.size(); ++i) {
+
+ // Get pointer to observation
+ GCTAObservation* obs = obs_list[i];
+
+ // Fill events into counts map
+ fill_maps_counts(obs);
+
+ // Optionally compute acceptance sky map
+ if (m_bkgsubtract != "NONE") {
+ fill_maps_acceptance(obs);
+ }
+
+ // Dispose events to free memory
+ obs->dispose_events();
+
+ } // endfor: looped over observations
// Return
return;
@@ -537,16 +753,17 @@ void ctskymap::get_parameters(void)
/***********************************************************************//**
- * @brief Map events into a sky map
+ * @brief Fill events into counts map
*
* @param[in] obs CTA observation.
*
* @exception GException::no_list
* No event list found in observation.
*
- * This method maps the events found in a CTA events list into a sky map.
+ * Fills the events found in a CTA events list into a sky map. The method
+ * adds the events to the m_counts member.
***************************************************************************/
-void ctskymap::map_events(GCTAObservation* obs)
+void ctskymap::fill_maps_counts(GCTAObservation* obs)
{
// Get non-const pointer on a CTA event list
GCTAEventList* events = dynamic_cast
@@ -589,11 +806,11 @@ void ctskymap::map_events(GCTAObservation* obs)
// Determine sky pixel
GCTAInstDir* inst = (GCTAInstDir*)&(event->dir());
GSkyDir dir = inst->dir();
- GSkyPixel pixel = m_skymap.dir2pix(dir);
+ GSkyPixel pixel = m_counts.dir2pix(dir);
// Skip if pixel is out of range
- if (pixel.x() < -0.5 || pixel.x() > (m_skymap.nx() - 0.5) ||
- pixel.y() < -0.5 || pixel.y() > (m_skymap.ny() - 0.5)) {
+ if (pixel.x() < -0.5 || pixel.x() > (m_counts.nx() - 0.5) ||
+ pixel.y() < -0.5 || pixel.y() > (m_counts.ny() - 0.5)) {
num_outside_map++;
continue;
}
@@ -605,49 +822,25 @@ void ctskymap::map_events(GCTAObservation* obs)
}
// Fill event in skymap
- m_skymap(pixel, 0) += 1.0;
+ #pragma omp critical(ctskymap_map_events)
+ m_counts(pixel) += 1.0;
num_in_map++;
-
+
} // endfor: looped over all events
// Log binning results
- log_value(NORMAL, "Events in list", obs->events()->size());
- log_value(NORMAL, "Events in map", num_in_map);
- log_value(NORMAL, "Events outside RoI", num_outside_roi);
- log_value(NORMAL, "Events outside map area", num_outside_map);
- log_value(NORMAL, "Events outside energies", num_outside_erange);
-
- // Write sky map into header
- log_header1(EXPLICIT, "Sky map");
- log_string(EXPLICIT, m_skymap.print(m_chatter));
-
- // Return
- return;
-}
-
-
-/***********************************************************************//**
- * @brief Estimates the background in sky map for observation
- *
- * @param[in] obs CTA observation.
- *
- * Estimates the background in the sky map for a given observation and adds
- * this estimate to the background sky map. The background estimation method
- * is specified by the "bkgsubtract" parameter which can take the following
- * values:
- *
- * NONE - No background estimation
- * IRF - Background estimation based on IRF template
- * RING - Ring background estimation
- ***************************************************************************/
-void ctskymap::map_background(GCTAObservation* obs)
-{
- // Dispatch to appropriate background estimation method
- if (m_bkgsubtract == "IRF") {
- map_background_irf(obs);
- }
- else if (m_bkgsubtract == "RING") {
- map_background_ring(obs);
+ #pragma omp critical(ctskymap_map_events)
+ {
+ log_header3(TERSE, get_obs_header(obs));
+ log_value(NORMAL, "Events in list", obs->events()->size());
+ log_value(NORMAL, "Events in map", num_in_map);
+ log_value(NORMAL, "Events outside RoI", num_outside_roi);
+ log_value(NORMAL, "Events outside map area", num_outside_map);
+ log_value(NORMAL, "Events outside energies", num_outside_erange);
+
+ // Write sky map into header
+ log_header1(EXPLICIT, "Sky map");
+ log_string(EXPLICIT, m_counts.print(m_chatter));
}
// Return
@@ -656,7 +849,7 @@ void ctskymap::map_background(GCTAObservation* obs)
/***********************************************************************//**
- * @brief Estimates the background in sky map based on IRF template
+ * @brief Compute background acceptance sky map based on IRF template
*
* @param[in] obs CTA observation.
*
@@ -664,14 +857,14 @@ void ctskymap::map_background(GCTAObservation* obs)
* No response information available for observation.
* No background template available in instrument response function.
*
- * Estimates the background in the sky map using the IRF template for a given
- * observation and adds this estimate to the background sky map. The IRF
- * template is integrated numerically in energy.
+ * Computes a background acceptance sky map using the IRF template for a
+ * given observation and adds the estimate to the acceptance sky map.
***************************************************************************/
-void ctskymap::map_background_irf(GCTAObservation* obs)
+void ctskymap::fill_maps_acceptance(GCTAObservation* obs)
{
// Get IRF response
- const GCTAResponseIrf* rsp = dynamic_cast(obs->response());
+ const GCTAResponseIrf* rsp = dynamic_cast
+ (obs->response());
// Throw an exception if observation has no instrument response function
if (rsp == NULL) {
@@ -694,60 +887,58 @@ void ctskymap::map_background_irf(GCTAObservation* obs)
throw GException::invalid_value(G_MAP_BACKGROUND_IRF, msg);
}
- // Compute natural logarithm of energy range in MeV
- double lnEmin = std::log(m_emin * 1.0e6);
- double lnEmax = std::log(m_emax * 1.0e6);
+ // Set minimum and maximum energy as GEnergy instances
+ GEnergy emin(m_emin, "TeV");
+ GEnergy emax(m_emax, "TeV");
- // Extract region of interest from observation
- GCTARoi roi = obs->roi();
+ // Extract region of interest center and radius from observation and
+ // compute the cosine of the RoI radius. This allows to mimimize as
+ // much as possible the trigonometric computations. Note that the cosine
+ // of an angle is maximal for an angle of zero. This explains later why
+ // we chose ">=" to test whether an angle is smaller than a given radius.
+ GSkyDir roi_centre = obs->roi().centre().dir();
+ double cos_roi_radius = std::cos(obs->roi().radius() * gammalib::deg2rad);
+
+ // Extract exposure
+ double exposure = obs->livetime();
// Initialise statistics
- int calls = 0;
double total = 0.0;
- // Loop over all background map pixels
- for (int i = 0; i < m_bkgmap.npix(); ++i) {
+ // Loop over all acceptance map pixels
+ for (int i = 0; i < m_acceptance.npix(); ++i) {
// Get sky direction of pixel
- GSkyDir skydir = m_bkgmap.inx2dir(i);
-
- // Convert sky direction in instrument direction
- GCTAInstDir instdir = obs->pointing().instdir(skydir);
+ GSkyDir& skydir = m_dirs[i];
- // If RoI is valid and instrument direction is not within RoI then
- // skip pixel
- if (roi.is_valid() && !roi.contains(instdir)) {
+ // Skip pixel if it is not within the RoI
+ if (roi_centre.cos_dist(skydir) < cos_roi_radius) {
continue;
}
- // Setup integration function
- ctskymap::irf_kern integrand(bkg, &instdir);
- GIntegral integral(&integrand);
-
- // Set precision (has been carefully adjusted using a test simulation
- // over the energy range 20 GeV - 120 TeV)
- integral.eps(1.0e-6);
+ // Convert sky direction to instrument direction
+ GCTAInstDir instdir = obs->pointing().instdir(skydir);
- // Do Romberg integration
- double value = integral.romberg(lnEmin, lnEmax);
+ // Compute background value
+ double value = bkg->rate_ebin(instdir, emin, emax);
- // Update number of background function calls
- calls += integral.calls();
-
// Multiply background rate with livetime and solid angle
- value *= obs->livetime() * m_solidangle[i];
+ value *= exposure * m_solidangle[i];
- // Add number of background events to background map
- m_bkgmap(i,0) += value;
+ // Add number of background events to acceptance map
+ #pragma omp critical(ctskymap_map_acceptance)
+ m_acceptance(i) += value;
// Update total number of background events
total += value;
- } // endfor: looped over background map pixels
+ } // endfor: looped over acceptance map pixels
- // Log background subtraction results
- log_value(NORMAL, "Events in background", int(total+0.5));
- log_value(NORMAL, "Background evaluations", calls);
+ // Log acceptance results
+ #pragma omp critical(ctskymap_map_events)
+ {
+ log_value(NORMAL, "Events in background", int(total+0.5));
+ }
// Return
return;
@@ -755,97 +946,99 @@ void ctskymap::map_background_irf(GCTAObservation* obs)
/***********************************************************************//**
- * @brief Estimates the background in sky map based on the ring background
- * method
+ * @brief Compute sky map, background map and significance map
*
- * @param[in] obs CTA observation.
+ * Computes the sky map, background map and significance map. The following
+ * background subtraction methods are supported:
*
- * @exception GException::invalid_value
- * No response information available for observation.
- * No background template available in instrument response function.
+ * NONE
+ * ====
+ *
+ * The sky map is simply the binned counts map, no background and
+ * significance maps are computed.
+ *
+ * IRF
+ * ===
+ *
+ * The background map is the acceptance map, and this background map is
+ * subtracted from the counts map and assigned to the sky map. The
+ * significance \f$\sigma_i\f$ is computed assuming Poisson statistic in
+ * the Gaussian limit, using
+ *
+ * \f[\sigma_i = \frac{N_i - B_i}{\sqrt{N_i}}\f]
*
- * Estimates the background in the sky map by summing the events within a
- * ring centered at a given pixel's position. The pixels in the ring are
- * weighted also by the background IRF value.
+ * where
+ * \f$N_i\f$ is the number of observed counts and
+ * \f$B_i\f$ is the estimated number of background counts.
+ *
+ * RING
+ * ====
+ *
+ * For the @c RING method the Li & Ma significance is computed for each
+ * pixel, using the @c roiradius, @c inradius and @c outradius parameters
+ * to specify the radius of the On region and the Off background ring,
+ * respectively. Depending on the @c usefft parameter, either an FFT is used
+ * to compute the number of events and the acceptance of the On and Off
+ * regions, or a direct computation is performed. While the former is faster
+ * it is less accurate since it assumes Euclidean distances in the sky map.
+ * The latter is slower but compute exact distance between pixels of the sky
+ * map. Use FFT if the sky map is close to a cartesian grid, and use the
+ * direct method if the sky map shows important distortions.
***************************************************************************/
-void ctskymap::map_background_ring(GCTAObservation* obs)
+void ctskymap::compute_maps(void)
{
- // Get IRF response (to scale background counts)
- const GCTAResponseIrf* rsp = dynamic_cast(obs->response());
+ // Write header into logger
+ log_header1(TERSE, "Compute maps");
- // Throw an exception if observation has no instrument response function
- if (rsp == NULL) {
- std::string msg = "No response information available for "+
- get_obs_header(obs)+" to compute IRF background. "
- "Please specify response information or use "
- "another background subtraction method.";
- throw GException::invalid_value(G_MAP_BACKGROUND_RING, msg);
+ // If no background subtraction is requested then simply assign the
+ // counts map to the skymap
+ if (m_bkgsubtract == "NONE") {
+ m_skymap = m_counts;
}
- // Get IRF background template
- const GCTABackground* bkg = rsp->background();
-
- // Throw an exception if observation has no IRF background template
- if (bkg == NULL) {
- std::string msg = "No IRF background template found in instrument "
- "response function for "+
- get_obs_header(obs)+". Please specify an instrument "
- "response function containing a background template.";
- throw GException::invalid_value(G_MAP_BACKGROUND_RING, msg);
+ // else if IRF background is selected then use the acceptance as background
+ else if (m_bkgsubtract == "IRF") {
+ m_skymap = m_counts - m_acceptance;
+ m_bkgmap = m_acceptance;
+ m_sigmap = (m_counts - m_acceptance) / sqrt(m_counts);
}
- // Compute natural logarithm of energy range in MeV
- double lnEmin = std::log(m_emin * 1.0e6);
- double lnEmax = std::log(m_emax * 1.0e6);
-
- // Extract region of interest from observation
- GCTARoi roi = obs->roi();
-
- // Initialise statistics
- int calls = 0;
- double exposure = obs->livetime();
-
- // Loop over all map pixels
- for (int i = 0; i < m_bkgmap.npix(); ++i) {
-
- // Get sky direction of pixel
- GSkyDir& skydir = m_dirs[i];
+ // else if RING background is selected then compute the On
+ else if (m_bkgsubtract == "RING") {
- // Convert sky direction in instrument direction
- GCTAInstDir instdir = obs->pointing().instdir(skydir);
-
- // If RoI is valid and instrument direction is not within RoI then
- // skip pixel
- if (roi.is_valid() && !roi.contains(instdir)) {
- continue;
+ // Compute initial RING background
+ if (m_usefft) {
+ compute_maps_ring_fft();
+ }
+ else {
+ compute_maps_ring_direct();
}
- // Setup integration function
- ctskymap::irf_kern integrand(bkg, &instdir);
- GIntegral integral(&integrand);
-
- // Set precision
- integral.eps(1.0e-6);
-
- // Do Romberg integration to get the sensitivity (note: this assumes
- // the background IRF is a good approximation for the radial
- // sensitivity)
- m_alphamap(i,0) += integral.romberg(lnEmin, lnEmax) *
- m_solidangle[i] * exposure;
-
- // Store the counts in this bin
- m_onmap(i,0) += m_skymap(i,0);
+ // Store copy of exclusion regions
+ GSkyMap exclmap = m_exclmap;
- // Update number of background function calls
- calls += integral.calls();
+ // Iterative computation of exclusion regions
+ for (int iter = 0; iter < m_iterations; ++iter) {
- } // endfor: Loop for caching bkg IRF sensitivity
+ // Set exclusion region from significance map
+ m_exclmap = exclmap;
+ for (int i = 0; i < m_exclmap.npix(); ++i) {
+ if (m_sigmap(i) > m_threshold) {
+ m_exclmap(i) = 1.0;
+ }
+ }
+
+ // Re-compute RING background
+ if (m_usefft) {
+ compute_maps_ring_fft();
+ }
+ else {
+ compute_maps_ring_direct();
+ }
- // Log background subtraction results
- log_value(NORMAL, "Alpha map evaluations", calls);
+ } // endfor: looped over iterations
- // Zero out the counts map to prepare it for the next observation
- m_skymap = 0.0;
+ }
// Return
return;
@@ -853,24 +1046,85 @@ void ctskymap::map_background_ring(GCTAObservation* obs)
/***********************************************************************//**
- * @brief Write keywords in FITS HDU
- *
- * @param[in,out] hdu Pointer to FITS HDU.
+ * @brief Compute the maps for RING background using a FFT
*
- * Writes keywords in FITS HDU.
+ * Computes the Li & Ma significance for each sky map pixel and replaces
+ * the sky and background maps by the On- and Off-count maps. The computation
+ * is done using a FFT.
***************************************************************************/
-void ctskymap::write_hdu_keywords(GFitsHDU* hdu) const
+void ctskymap::compute_maps_ring_fft(void)
{
- // Continue only if pointer is valid
- if (hdu != NULL) {
+ // Log message about what is being done
+ log_header3(NORMAL, "Computing ring background map (FFT method)");
- // Set keywords
- hdu->card("BKGSUB", m_bkgsubtract, "Background substraction method");
- hdu->card("E_MIN", m_emin, "[TeV] Lower energy boundary");
- hdu->card("E_MAX", m_emax, "[TeV] Upper energy boundary");
- hdu->card("EUNIT", "TeV", "Units for E_MIN and E_MAX");
+ // Initialise maps
+ m_skymap = m_counts;
+ m_bkgmap = m_counts;
+ m_sigmap = m_counts;
+ m_skymap = 0.0;
+ m_bkgmap = 0.0;
+ m_sigmap = 0.0;
- } // endif: pointer was valid
+ // Initialise statistics
+ int num_bad_alpha = 0;
+
+ // Set-up sky and acceptance maps with excluded pixels according to the
+ // exclusion map
+ GSkyMap excl_counts = m_counts;
+ GSkyMap excl_acceptance = m_acceptance;
+ for (int i = 0; i < m_counts.npix(); ++i) {
+ if (m_exclmap(i) != 0.0) {
+ excl_counts(i) = 0.0;
+ excl_acceptance(i) = 0.0;
+ }
+ }
+
+ // Convolve maps for On-region and Off-region
+ GSkyMap on_counts = ring_convolve(m_counts, 0.0, m_roiradius);
+ GSkyMap on_alpha = ring_convolve(m_acceptance, 0.0, m_roiradius);
+ GSkyMap off_counts = ring_convolve(excl_counts, m_inradius, m_outradius);
+ GSkyMap off_alpha = ring_convolve(excl_acceptance, m_inradius, m_outradius);
+
+ // Compute maps by loop over sky map pixels
+ for (int i = 0; i < m_sigmap.npix(); ++i) {
+
+ // Get On-counts, Off-counts and alpha value for this bin
+ double n_on = on_counts(i);
+ double n_off = off_counts(i);
+ double alpha;
+ if (off_alpha(i) == 0.0) {
+ alpha = 0.0;
+ n_on = 0.0;
+ n_off = 0.0;
+ }
+ else {
+ alpha = on_alpha(i) / off_alpha(i);
+ }
+
+ // Store the On and alpha-weighted Off-counts
+ m_bkgmap(i) = alpha * n_off;
+ m_skymap(i) = n_on - m_bkgmap(i);
+
+ // If alpha is zero then increment the bad alpha counter
+ if (alpha == 0.0) {
+ num_bad_alpha++;
+ }
+
+ // ... otherwise compute and store the results
+ else {
+ m_sigmap(i) = sigma_li_ma(n_on, n_off, alpha);
+ }
+
+ } // endfor: looped over all pixels
+
+ // Log the number of bad-alpha bins
+ log_value(NORMAL, "Total number of pixels", m_sigmap.npix());
+ log_value(NORMAL, "Pixels with alpha=0", num_bad_alpha);
+
+ // Now take the square root. Since some bins can be negative, first
+ // take the absolute value of the map, square root that, then multiply
+ // each bin by its sign to preserve the +/- significance.
+ m_sigmap = sign(m_sigmap) * sqrt(abs(m_sigmap));
// Return
return;
@@ -878,89 +1132,85 @@ void ctskymap::write_hdu_keywords(GFitsHDU* hdu) const
/***********************************************************************//**
- * @brief Computes the bin-by-bin significance if background subtraction
- * requested
+ * @brief Compute the pixel significance for RING background
*
- * Method computes the bin-by bin significance if a background subtraction
- * method is specified. If method is "IRF", Poisson statistics (in the
- * Gaussian limit) are assumed. If method is "RING" a Li & Ma significance
- * is computed for each bin (as well as a map of the computed alpha values).
+ * Computes the Li & Ma significance for each sky map pixel and replaces
+ * the sky and background maps by the On- and Off-count maps. The computation
+ * is done by computing the exact distances between pixels in the sky map.
***************************************************************************/
-void ctskymap::map_significance(void)
+void ctskymap::compute_maps_ring_direct(void)
{
- // Compute significance from "RING" method (Li & Ma eq. 17)
- if (m_bkgsubtract == "RING") {
-
- // Log message about what is being done
- log_header1(NORMAL, "Computing Ring background map");
- log_value(NORMAL, "Total pixels to process", m_onmap.npix());
-
- // Store the number of bins with an inappropriate alpha parameter
- int num_bad_alpha = 0;
-
- // Loop through each bin in the on-counts map
- for (int i = 0; i < m_onmap.npix(); ++i) {
-
- // Initialise the on/off counts and alpha for this bin
- double n_on = 0.0;
- double n_off = 0.0;
- double alpha = 0.0;
+ // Log message about what is being done
+ log_header3(NORMAL, "Computing ring background map (direct method)");
+ log_value(NORMAL, "Total pixels to process", m_counts.npix());
+
+ // Initialise maps
+ m_skymap = m_counts;
+ m_bkgmap = m_counts;
+ m_sigmap = m_counts;
+ m_skymap = 0.0;
+ m_bkgmap = 0.0;
+ m_sigmap = 0.0;
- // Since this can take a long time, keep the user updated on
- // the progress when another 10% of pixels is processed
- if (i % (m_onmap.npix()/10) == 0) {
- log_value(NORMAL, "Pixels remaining", m_onmap.npix()-i);
- }
+ // Initialise statistics
+ int num_bad_alpha = 0;
+
+ // Set progress logging frequency. By default the user is updated on the
+ // progress when another 10% of pixels is processed, however, if there
+ // are many pixels the user is updated each 100000 pixels.
+ int n_logpix = m_counts.npix() / 10;
+ if (n_logpix > 100000) {
+ n_logpix = 100000;
+ }
- // Get bin coordinates
- GSkyDir& skydir = m_dirs[i];
+ // Compute maps by loop over sky map pixels
+ #pragma omp parallel for
+ for (int i = 0; i < m_sigmap.npix(); ++i) {
- // Compute the alpha and counts for this bin
- compute_ring_values(m_onmap, m_alphamap, skydir,
- n_on, n_off, alpha);
+ // Initialise the on/off counts and alpha for this bin
+ double n_on = 0.0;
+ double n_off = 0.0;
+ double alpha = 0.0;
- // If alpha is zero then increment the bad alpha counter
- if (alpha == 0.0) {
- num_bad_alpha++;
- }
+ // Update user about progress (only if there are more than 10000 pixels)
+ if (i % n_logpix == 0 && n_logpix > 10000) {
+ #pragma omp critical(ctskymap_map_significance_ring)
+ log_value(NORMAL, "Pixels remaining", m_sigmap.npix()-i);
+ }
- // ... otherwise store the results
- else {
+ // Compute the alpha and counts for this bin
+ compute_ring_values(i, m_counts, m_acceptance, n_on, n_off, alpha);
- // Store the on-counts and alpha-weighted off-counts
- m_skymap(i,0) = n_on;
- m_bkgmap(i,0) = alpha * n_off;
+ // Store the On and alpha-weighted Off-counts
+ m_bkgmap(i) = alpha * n_off;
+ m_skymap(i) = n_on - m_bkgmap(i);
- // Compute and store significance (Li & Ma eq. 17)
- if (n_on == 0.0) {
- m_sigmap(i,0) = -2.0 * n_off * std::log(1.0+alpha);
- }
- else {
- m_sigmap(i,0) = (n_on < (alpha*n_off) ? -2.0 : 2.0) *
- (n_on * std::log((1.0+alpha) *
- n_on / (alpha * (n_on+n_off))) +
- n_off * std::log((1.0+alpha) *
- n_off / (n_on+n_off)));
- }
+ // If alpha is zero then increment the bad alpha counter
+ if (alpha == 0.0) {
+ num_bad_alpha++;
+ #pragma omp critical(ctskymap_compute_maps_ring_direct)
+ m_sigmap(i) = 0.0;
+ }
- } // endelse: alpha was non-zero
+ // ... otherwise store the results
+ else {
- } // endfor: looped over all pixels
+ // Compute and store significance (Li & Ma eq. 17)
+ double sigma = sigma_li_ma(n_on, n_off, alpha);
+ #pragma omp critical(ctskymap_compute_maps_ring_direct)
+ m_sigmap(i) = sigma;
- // Log the number of bad-alpha bins
- log_value(NORMAL, "Bins with alpha=0", num_bad_alpha);
+ } // endelse: alpha was non-zero
- // Now take the square root. Since some bins can be negative, first
- // take the absolute value of the map, square root that, then multiply
- // each bin by its sign to preserve the +/- significance.
- m_sigmap = sign(m_sigmap) * sqrt(abs(m_sigmap));
+ } // endfor: looped over all pixels
- } // endif: ring background selected
+ // Log the number of bad-alpha bins
+ log_value(NORMAL, "Pixels with alpha=0", num_bad_alpha);
- // Compute significance using Poisson statistics in the Gaussian limit
- else {
- m_sigmap = (m_skymap - m_bkgmap) / sqrt(m_skymap);
- }
+ // Now take the square root. Since some bins can be negative, first
+ // take the absolute value of the map, square root that, then multiply
+ // each bin by its sign to preserve the +/- significance.
+ m_sigmap = sign(m_sigmap) * sqrt(abs(m_sigmap));
// Return
return;
@@ -970,66 +1220,107 @@ void ctskymap::map_significance(void)
/***********************************************************************//**
* @brief Computes Non, Noff and alpha for a counts map and sensitivity map
*
- * @param[in] counts Counts map
- * @param[in] sensitivity Sensitivity map
- * @param[in] position Position of the ring & roi centers
- * @param[out] non Returned estimate of ON counts
- * @param[out] noff Returned estimate of OFF counts
- * @param[out] alpha Returned estimate of alpha
+ * @param[in] ipixel Sky map pixel to consider.
+ * @param[in] counts Counts map.
+ * @param[in] background Background map.
+ * @param[out] non Returned estimate of On-counts.
+ * @param[out] noff Returned estimate of Off-counts.
+ * @param[out] alpha Returned estimate of alpha.
*
- * This method computes the Non counts Noff values at a given position. It
- * also computes the alpha parameter from the passed sensitivity map.
+ * Computes the On- and Off-counts values at a given position from in counts
+ * map and the alpha parameter from the background map.
+ *
+ * To speed-up the computations the method considers only the pixels within
+ * a bounding box that comprises the outer background ring radius. The method
+ * considers wrapping around of pixel indices in the Right Ascension /
+ * Galactic longitude direction.
+ *
+ * If the alpha of the ring region is zero the values of @p non, @p noff,
+ * and @p alpha will be set to zero.
***************************************************************************/
-void ctskymap::compute_ring_values(const GSkyMap& counts,
- const GSkyMap& sensitivity,
- const GSkyDir& position,
+void ctskymap::compute_ring_values(const int& ipixel,
+ const GSkyMap& counts,
+ const GSkyMap& background,
double& non,
double& noff,
double& alpha)
{
- // Reset Non and Noff
- non = 0.0;
- noff = 0.0;
+ // Initialise return values (non, noff and alpha)
+ non = 0.0;
+ noff = 0.0;
+ alpha = 0.0;
// Initialise On and Off alpha
double alpha_on = 0.0;
double alpha_off = 0.0;
- // Define the regions necessary to do the pixel checks
- GSkyRegionCircle roi_reg(position, m_roiradius);
- GSkyRegionCircle inner_reg(position, m_inradius);
- GSkyRegionCircle outer_reg(position, m_outradius);
+ // Get sky direction of pixel
+ GSkyDir& position = m_dirs[ipixel];
- // Loop over every pixel in the observation to compute Non, Noff
- for (int j = 0; j < counts.npix(); ++j) {
+ // Get pixel bounding box
+ int ix_start;
+ int ix_stop;
+ int iy_start;
+ int iy_stop;
+ ring_bounding_box(ipixel, ix_start, ix_stop, iy_start, iy_stop);
- // Get the index and sky direction of this pixel
- GSkyDir& skydir = m_dirs[j];
+ // Get number of x pixels in sky map
+ int nx = counts.nx();
- // Check if pixel is inside the background region
- if ((m_exclmap(j,0) == 0.0) && outer_reg.contains(skydir) &&
- !inner_reg.contains(skydir)) {
+ // Loop over x pixels of bounding box
+ for (int ix_nowrap = ix_start; ix_nowrap < ix_stop; ++ix_nowrap) {
- // Update n_off
- noff += counts(j,0);
+ // Compute x pixel index by wrapping the pixel index into the
+ // interval [0,nx[
+ int ix = ix_nowrap;
+ if (ix < 0) {
+ ix += nx;
+ }
+ else if (ix >= nx) {
+ ix -= nx;
+ }
- // Update alpha_off
- alpha_off += sensitivity(j,0);
+ // Initialise pixel index
+ int i = ix + iy_start * nx;
- }
+ // Loop over y pixels of bounding box
+ for (int iy = iy_start; iy < iy_stop; ++iy, i += nx) {
+
+ // Get the index and sky direction of this pixel
+ GSkyDir& skydir = m_dirs[i];
+
+ // Only consider pixels within the outer radius of the background
+ // region
+ if (position.cos_dist(skydir) >= m_cos_outradius) {
+
+ // Check if pixel is inside the background region
+ if ((m_exclmap(i) == 0.0) &&
+ (position.cos_dist(skydir) < m_cos_inradius)) {
+
+ // Update n_off
+ noff += counts(i);
+
+ // Update alpha_off
+ alpha_off += background(i);
+
+ }
+
+ // ... otherwise check if pixel is inside source region
+ else if (position.cos_dist(skydir) >= m_cos_roiradius) {
+
+ // Update n_on
+ non += counts(i);
- // ... otherwise check if pixel is inside source region
- else if (roi_reg.contains(skydir)) {
+ // Update alpha_on
+ alpha_on += background(i);
- // Update n_on for significance computation
- non += counts(j,0);
+ } // endif: source and background region check
- // Update alpha_on
- alpha_on += sensitivity(j,0);
+ } // endif: pixel is within outer radius of background region
- } // endif: source and background region check
+ } // endfor: looped over y pixels
- } // endfor: looped over pixels
+ } // endfor: looped over x pixels
// Compute alpha. If the off region does not have any sensitivity then
// set Non = Noff = 0
@@ -1048,34 +1339,100 @@ void ctskymap::compute_ring_values(const GSkyMap& counts,
/***********************************************************************//**
- * @brief Generates map of pixel exclusions
+ * @brief Computes bounding box for RING background computation
*
- * @param[in] filename Exclusion file name.
+ * @param[in] ipixel Sky map pixel to consider.
+ * @param[out] ix1 Index of first pixel in x.
+ * @param[out] ix2 Index after last pixel in x.
+ * @param[out] iy1 Index of first pixel in y.
+ * @param[out] iy2 Index after last pixel in y.
*
- * Generates a sky map of the pixels that are to be excluded from the
- * background estimation. Pixels with values different from 0 will be
- * excluded.
+ * Computes the bounding box that contained the background ring for a
+ * specific pixel for the RING background method.
+ *
+ * The method determines the local pixel scale for the requested pixel and
+ * draws a bounding box with 1.5 times the outer background ring radius
+ * around the pixel. In the x direction the pixel indices are unconstrained,
+ * and the client has to assure that the pixel value is comprised within
+ * the validity range (this is required to handle the longitude wrap around).
+ * In the y direction the pixel indices are constrained to [0,ny], where
+ * ny is the number of y pixels in the sky map.
***************************************************************************/
-void ctskymap::map_exclusions(const GFilename& filename)
+void ctskymap::ring_bounding_box(const int& ipixel, int& ix1, int& ix2,
+ int& iy1, int& iy2) const
{
- // Create exlusion map
- m_exclmap = create_map(m_obs);
+ // Get number of pixels in x and y direction
+ int nx = m_counts.nx();
+ int ny = m_counts.ny();
+
+ // Get x and y index of pixel to consider
+ int ix0 = ipixel % nx;
+ int iy0 = ipixel / nx;
+
+ // Get pointer on WCS projection
+ const GWcs* wcs = dynamic_cast(m_counts.projection());
+ if (wcs == NULL) {
+ std::string msg = "Sky map is not a WCS projection. Method is only "
+ "valid for WCS projections.";
+ throw GException::invalid_value(G_RING_BOUNDING_BOX, msg);
+ }
- // Set all pixels to 0 (no pixel excluded)
- m_exclmap = 0.0;
+ // Get X and Y step size
+ double dx_binsz = wcs->cdelt(0);
+ double dy_binsz = wcs->cdelt(1);
- // Make sure the exclusions filename is valid
- if (is_valid_filename(filename)) {
+ // Compute pixel increment in x-direction
+ double dx = 0.0;
+ int ix = (ix0 > 0) ? ix0 - 1 : ix0 + 1;
+ if (ix < nx) {
+ dx = m_dirs[ipixel].dist_deg(m_dirs[ix+iy0*nx]);
+ }
+ if (dx < dx_binsz) {
+ dx = dx_binsz;
+ }
- // Fill the exclusions based on the regions supplied
- if (filename.is_fits()) {
- map_exclusions_fits(filename);
+ // Compute pixel increment in y-direction
+ double dy = 0.0;
+ int iy = (iy0 > 0) ? iy0 - 1 : iy0 + 1;
+ if (iy < ny) {
+ dx = m_dirs[ipixel].dist_deg(m_dirs[ix0+iy*nx]);
+ }
+ if (dy < dy_binsz) {
+ dy = dy_binsz;
+ }
+
+ // Compute bounding box half size in x and y. The outer radius is
+ // multiplied by 1.5 to have some margin in case of map distortions
+ int nx_bbox = int(1.5 * m_outradius / dx);
+ int ny_bbox = int(1.5 * m_outradius / dy);
+
+ // Compute index range for x. We allow that indices are smaller than 0
+ // or equal or larger than nx to handle wrap around, which needs to be
+ // done in the client method.
+ if (2*nx_bbox < nx) {
+ ix1 = ix0 - nx_bbox;
+ ix2 = ix0 + nx_bbox;
+ }
+ else {
+ ix1 = 0;
+ ix2 = nx;
+ }
+
+ // Compute index range for y. The index range is restricted to [0,ny].
+ if (2*ny_bbox < ny) {
+ iy1 = iy0 - ny_bbox;
+ iy2 = iy0 + ny_bbox;
+ if (iy1 < 0) {
+ iy1 = 0;
}
- else {
- map_exclusions_reg(filename);
+ if (iy2 > ny) {
+ iy2 = ny;
}
-
- } // endif: filename was valid
+ }
+ else {
+ iy1 = 0;
+ iy2 = ny;
+ }
// Return
return;
@@ -1083,64 +1440,166 @@ void ctskymap::map_exclusions(const GFilename& filename)
/***********************************************************************//**
- * @brief Fills exclusions map from FITS image
+ * @brief Return FFT kernel for background ring
*
- * @param[in] filename FITS image file name.
+ * @param[in] map Sky map to be convolved with a background ring.
+ * @param[in] rmin Minimum ring radius (degrees).
+ * @param[in] rmax Maximum ring radius (degrees).
+ * @return FFT kernel.
*
- * Sets all exclusion map pixels to 1 that correspond to non-zero pixels in
- * the exclusion sky map FITS file.
+ * Computes the FFT kernel for a backgrund ring.
***************************************************************************/
-void ctskymap::map_exclusions_fits(const GFilename& filename)
+GSkyMap ctskymap::ring_convolve(const GSkyMap& map, const double& rmin,
+ const double& rmax) const
{
- // Load the fits image
- GSkyMap inmap(filename);
+ // Copy input map
+ GSkyMap conv_map = map;
- // Loop through the individual pixels in the exclusion map
- for (int i = 0; i < m_exclmap.npix(); ++i) {
+ // Get FFT of ring kernel
+ GFft fft_kernel(ring_kernel(rmin, rmax));
- // Get the pixel direction
- GSkyDir dir = m_exclmap.pix2dir(i);
+ // Extract sky map into Ndarray
+ GNdarray array(conv_map.nx(), conv_map.ny());
+ double *dst = array.data();
+ for (int i = 0; i < conv_map.npix(); ++i) {
+ *dst++ = conv_map(i);
+ }
- // Check this sky position in the fits map
- if (inmap.contains(dir) && (inmap(inmap.dir2inx(dir)) != 0.0)) {
+ // FFT of sky map
+ GFft fft_array = GFft(array);
- // Set the pixel to 1
- m_exclmap(i) = 1.0;
+ // Multiply FFT of sky map with FFT of kernel
+ GFft fft_smooth = fft_array * fft_kernel;
- } // endif: pixel,region overlap check
+ // Backward transform sky map
+ GNdarray smooth = fft_smooth.backward();
- } // endfor: looped over exclusion map pixels
+ // Insert sky map
+ double *src = smooth.data();
+ for (int i = 0; i < conv_map.npix(); ++i) {
+ conv_map(i) = *src++;
+ }
- // Return
- return;
+ // Return convolved map
+ return conv_map;
}
/***********************************************************************//**
- * @brief Fills exclusions map from DS9 region file
+ * @brief Return FFT kernel for background ring
*
- * @param[in] filename DS9 region file name.
+ * @param[in] rmin Minimum ring radius (degrees).
+ * @param[in] rmax Maximum ring radius (degrees).
+ * @return FFT kernel.
*
- * Sets all exclusion map pixels to 1 that are contained in any of the DS9
- * regions.
+ * @exception GException::invalid_value
+ * Sky map is not a WCS projection.
+ *
+ * Computes the FFT kernel for a background ring. The computation is done
+ * assuming that the sky map represents a cartesian grid, and distances are
+ * computed in that grid using Euclidean distances. The @p rmin and @p rmax
+ * arguments refer to a ring radius in these Euclidean distances.
***************************************************************************/
-void ctskymap::map_exclusions_reg(const GFilename& filename)
+GNdarray ctskymap::ring_kernel(const double& rmin, const double& rmax) const
{
- // Load the exclusion regions
- GSkyRegions regions(filename);
+ // Get pointer on WCS projection
+ const GWcs* wcs = dynamic_cast(m_counts.projection());
+ if (wcs == NULL) {
+ std::string msg = "Sky map is not a WCS projection. Method is only "
+ "valid for WCS projections.";
+ throw GException::invalid_value(G_RING_KERNEL, msg);
+ }
- // Loop through the individual pixels in the exclusion map
- for (int i = 0; i < m_exclmap.npix(); ++i) {
+ // Get X and Y step size
+ double dx = wcs->cdelt(0);
+ double dy = wcs->cdelt(1);
+
+ // Initialise kernel
+ GNdarray kern(m_counts.nx(), m_counts.ny());
+
+ // Fill kernel
+ for (int ix1 = 0, ix2 = m_counts.nx(); ix1 < m_counts.nx(); ++ix1, --ix2) {
+ double x = ix1 * dx;
+ double xqs = x * x;
+ for (int iy1 = 0, iy2 = m_counts.ny(); iy1 < m_counts.ny(); ++iy1, --iy2) {
+ double y = iy1 * dy;
+ double r = std::sqrt(xqs + y*y);
+ if ((r >= rmin) && (r <= rmax)) {
+ kern(ix1,iy1) += 1.0;
+ if (ix2 < m_counts.nx()) {
+ kern(ix2,iy1) += 1.0;
+ }
+ if (iy2 < m_counts.ny()) {
+ kern(ix1,iy2) += 1.0;
+ }
+ if ((ix2 < m_counts.nx()) && (iy2 < m_counts.ny())) {
+ kern(ix2,iy2) += 1.0;
+ }
+ }
+ }
+ }
- // Get the pixel position
- GSkyDir dir = m_exclmap.pix2dir(i);
+ // Return kernel
+ return kern;
+}
- // If pixel position overlaps with the regions
- if (regions.contains(dir)) {
- m_exclmap(i) = 1.0;
- }
- } // endfor: looped over exclusion map pixels
+/***********************************************************************//**
+ * @brief Compute significance following Li & Ma
+ *
+ * @param[in] n_on Number of On-counts.
+ * @param[in] n_off Number of Off-counts.
+ * @param[in] alpha Alpha parameters.
+ * @return Significance.
+ *
+ * Computes the significance following Li & Ma, Equation (17).
+ ***************************************************************************/
+double ctskymap::sigma_li_ma(const double& n_on,
+ const double& n_off,
+ const double& alpha) const
+{
+ // Allocate result
+ double sigma;
+
+ // Handle special case of no On-counts
+ if (n_on == 0.0) {
+ sigma = -2.0 * n_off * std::log(1.0+alpha);
+ }
+
+ // ... otherwise handle general case
+ else {
+ sigma = (n_on < (alpha*n_off) ? -2.0 : 2.0) *
+ (n_on * std::log((1.0+alpha) * n_on / (alpha * (n_on+n_off))) +
+ n_off * std::log((1.0+alpha) * n_off / (n_on+n_off)));
+ }
+
+ // Return sigma
+ return sigma;
+}
+
+
+/***********************************************************************//**
+ * @brief Write sky map into FITS file
+ *
+ * @param[in,out] fits FITS file.
+ * @param[in] map Sky map.
+ * @param[in] extname Extension name.
+ *
+ * Write one sky map with the extension name and keywords into the FITS file.
+ ***************************************************************************/
+void ctskymap::write_map(GFits& fits, const GSkyMap& map, const std::string& extname) const
+{
+ // Write map into FITS file
+ GFitsHDU* hdu = map.write(fits);
+
+ // Set map extension name
+ if (hdu != NULL) {
+ hdu->extname(extname);
+ }
+
+ // Write keywords into map extension
+ write_ogip_keywords(hdu);
+ write_hdu_keywords(hdu);
// Return
return;
@@ -1148,21 +1607,36 @@ void ctskymap::map_exclusions_reg(const GFilename& filename)
/***********************************************************************//**
- * @brief Estimates the background in sky map based on IRF template
+ * @brief Write keywords in FITS HDU
*
- * @param[in] lnE Natural logarithm of energy in MeV.
+ * @param[in,out] hdu Pointer to FITS HDU.
+ *
+ * Writes keywords in FITS HDU.
***************************************************************************/
-double ctskymap::irf_kern::eval(const double& lnE)
+void ctskymap::write_hdu_keywords(GFitsHDU* hdu) const
{
- // Get log10 of energy in TeV
- double logE = lnE * gammalib::inv_ln10 - 6;
+ // Continue only if pointer is valid
+ if (hdu != NULL) {
+
+ // Set keywords
+ hdu->card("BKGSUB", m_bkgsubtract, "Background substraction method");
+ hdu->card("E_MIN", m_emin, "[TeV] Lower energy boundary");
+ hdu->card("E_MAX", m_emax, "[TeV] Upper energy boundary");
+ hdu->card("EUNIT", "TeV", "Units for E_MIN and E_MAX");
- // Get function value
- double value = (*m_bgd)(logE, m_dir->detx(), m_dir->dety());
+ // Write RING background keywords
+ if (m_bkgsubtract == "RING") {
+ hdu->card("USEFFT", m_usefft, "Use FFT for RING background");
+ hdu->card("ROIRAD", m_roiradius, "[deg] Source region radius");
+ hdu->card("INRAD", m_inradius, "[deg] Inner background ring radius");
+ hdu->card("OUTRAD", m_outradius, "[deg] Outer background ring radius");
+ hdu->card("ITER", m_iterations, "Exclusion map iterations");
+ hdu->card("THRES", m_threshold, "Exclusion map threshold");
+ hdu->card("EXCLMAP", m_inexclusion.url(), "Exclusion map name");
+ }
- // Correct for variable substitution
- value *= std::exp(lnE);
+ } // endif: pointer was valid
- // Return value
- return value;
+ // Return
+ return;
}
diff --git a/src/ctskymap/ctskymap.hpp b/src/ctskymap/ctskymap.hpp
index c3a1300..f23e292 100755
--- a/src/ctskymap/ctskymap.hpp
+++ b/src/ctskymap/ctskymap.hpp
@@ -1,7 +1,7 @@
/***************************************************************************
* ctskymap - Sky mapping tool *
* ----------------------------------------------------------------------- *
- * copyright (C) 2010-2017 by Juergen Knoedlseder *
+ * copyright (C) 2010-2018 by Juergen Knoedlseder *
* ----------------------------------------------------------------------- *
* *
* This program is free software: you can redistribute it and/or modify *
@@ -64,59 +64,69 @@ public:
protected:
// Protected methods
- void init_members(void);
- void copy_members(const ctskymap& app);
- void free_members(void);
- void get_parameters(void);
- void map_events(GCTAObservation* obs);
- void map_background(GCTAObservation* obs);
- void map_background_irf(GCTAObservation* obs);
- void map_background_ring(GCTAObservation* obs);
- void map_significance(void);
- void map_exclusions(const GFilename& filename);
- void map_exclusions_fits(const GFilename& filename);
- void map_exclusions_reg(const GFilename& filename);
- void compute_ring_values(const GSkyMap& counts,
- const GSkyMap& sensitivity,
- const GSkyDir& position,
- double& non, double& noff, double& alpha);
- void write_hdu_keywords(GFitsHDU* hdu) const;
-
- // Background integration kernel
- class irf_kern : public GFunction {
- public:
- irf_kern(const GCTABackground* bgd,
- const GCTAInstDir* dir) :
- m_bgd(bgd),
- m_dir(dir) { }
- double eval(const double& lnE);
- protected:
- const GCTABackground* m_bgd; //!< Pointer to background
- const GCTAInstDir* m_dir; //!< Pointer to instrument direction
- };
+ void init_members(void);
+ void copy_members(const ctskymap& app);
+ void free_members(void);
+ void get_parameters(void);
+ void setup_maps(void);
+ void setup_exclusion_map(const GFilename& filename);
+ void setup_exclusion_map_fits(const GFilename& filename);
+ void setup_exclusion_map_region(const GFilename& filename);
+ void fill_maps(void);
+ void fill_maps_counts(GCTAObservation* obs);
+ void fill_maps_acceptance(GCTAObservation* obs);
+ void compute_maps(void);
+ void compute_maps_ring_fft(void);
+ void compute_maps_ring_direct(void);
+ void compute_ring_values(const int& ipixel,
+ const GSkyMap& counts,
+ const GSkyMap& background,
+ double& non,
+ double& noff,
+ double& alpha);
+ void ring_bounding_box(const int& ipixel, int& ix1, int& ix2,
+ int& iy1, int& iy2) const;
+ GSkyMap ring_convolve(const GSkyMap& map, const double& rmin,
+ const double& rmax) const;
+ GNdarray ring_kernel(const double& rmin, const double& rmax) const;
+ double sigma_li_ma(const double& n_on,
+ const double& n_off,
+ const double& alpha) const;
+ void write_map(GFits& fits,
+ const GSkyMap& map,
+ const std::string& extname) const;
+ void write_hdu_keywords(GFitsHDU* hdu) const;
// User parameters
GFilename m_outmap; //!< Output file name
+ GFilename m_inexclusion; //!< Exclusion map file name
double m_emin; //!< Minimum energy (TeV)
double m_emax; //!< Maximum energy (TeV)
std::string m_bkgsubtract; //!< Background subtraction method
double m_roiradius; //!< Region of interest radius for RING bkg.
- double m_inradius; //!< Inner ring radius for RING bkg.
- double m_outradius; //!< Outer ring radius for RING bkg.
+ double m_inradius; //!< Inner ring radius for RING background
+ double m_outradius; //!< Outer ring radius for RING background
+ int m_iterations; //!< Number of iterations for RING background
+ double m_threshold; //!< Threshold for RING background
+ bool m_usefft; //!< Use FFT for RING background
bool m_publish; //!< Publish sky map?
GChatter m_chatter; //!< Chattiness
// Protected members
+ bool m_has_inmap; //!< Has valid input map
GSkyMap m_skymap; //!< Sky map
GSkyMap m_bkgmap; //!< Background map
GSkyMap m_sigmap; //!< Significance map
+ GSkyMap m_counts; //!< Counts map
+ GSkyMap m_acceptance; //!< Acceptance map
GSkyMap m_exclmap; //!< Exclusion map for RING background
- GSkyMap m_alphamap; //!< Alpha map values for RING background
- GSkyMap m_onmap; //!< On counts map for RING background
// Caching variables to prevent multiple computations of the same thing
- std::vector m_solidangle; //!< Cached pixel solid angles
- std::vector m_dirs; //!< Cached pixel directions
+ std::vector m_solidangle; //!< Cached pixel solid angles
+ std::vector m_dirs; //!< Cached pixel directions
+ double m_cos_roiradius; //!< Cosine of RoI radius
+ double m_cos_inradius; //!< Cosine of inner ring radius
+ double m_cos_outradius; //!< Cosine of outer ring radius
};
diff --git a/src/ctskymap/ctskymap.par b/src/ctskymap/ctskymap.par
index 86e8398..37633fb 100755
--- a/src/ctskymap/ctskymap.par
+++ b/src/ctskymap/ctskymap.par
@@ -19,6 +19,7 @@
inobs, f, a, events.fits,,, "Input event list or observation definition XML file"
caldb, s, a, prod2,,, "Calibration database"
irf, s, a, South_0.5h,,, "Instrument response function"
+inmap, f, h, NONE,,, "Input skymap file"
outmap, f, a, skymap.fits,,, "Output skymap file"
#
@@ -43,10 +44,13 @@ yref, r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC
# Background subtraction parameters
#==================================
bkgsubtract, s, a, NONE, NONE|IRF|RING,, "Background subtraction method"
-roiradius, r, a, 0.2,,, "Source region radius for estimating on-counts (degrees)"
-inradius, r, a, 0.6,,, "Inner ring radius (degrees)"
-outradius, r, a, 0.8,,, "Outer ring radius (degrees)"
+roiradius, r, a, 0.1,,, "Source region radius for estimating on-counts (degrees)"
+inradius, r, a, 0.6,,, "Inner background ring radius (degrees)"
+outradius, r, a, 0.8,,, "Outer background ring radius (degrees)"
+iterations, i, a, 0,0,100, "Number of iterations for exclusion regions computation"
+threshold, r, a, 5.0,,, "Significance threshold for exclusion regions computation"
inexclusion, f, h, NONE,,, "Exclusion map or DS9 region file name"
+usefft, b, h, yes,,, "Use FFT for ring computation?"
#
# Standard parameters
diff --git a/src/support/ctobservation.cpp b/src/support/ctobservation.cpp
index d15ee6f..e6511ec 100755
--- a/src/support/ctobservation.cpp
+++ b/src/support/ctobservation.cpp
@@ -319,6 +319,48 @@ const GCTAObservation* ctobservation::next_unbinned_observation(void) const
/***********************************************************************//**
+ * @brief Read OGIP keywords from FITS HDU
+ *
+ * @param[in,out] hdu Pointer to FITS HDU.
+ *
+ * Read OGIP keywords from FITS HDU.
+ ***************************************************************************/
+void ctobservation::read_ogip_keywords(GFitsHDU* hdu) const
+{
+ // Continue only if pointer is valid
+ if (hdu != NULL) {
+
+ // Get observation information
+ m_ogip_telescope = (hdu->has_card("TELESCOP")) ? hdu->string("TELESCOP") : "";
+
+ // Get observation time information
+ std::string date_obs = (hdu->has_card("DATE-OBS")) ? hdu->string("DATE-OBS") : "";
+ std::string time_obs = (hdu->has_card("TIME-OBS")) ? hdu->string("TIME-OBS") : "";
+ std::string date_end = (hdu->has_card("DATE-END")) ? hdu->string("DATE-END") : "";
+ std::string time_end = (hdu->has_card("TIME-END")) ? hdu->string("TIME-END") : "";
+ m_ogip_telapse = (hdu->has_card("TELAPSE")) ? hdu->real("TELAPSE") : 0.0;
+ m_ogip_ontime = (hdu->has_card("ONTIME")) ? hdu->real("ONTIME") : 0.0;
+ m_ogip_livetime = (hdu->has_card("LIVETIME")) ? hdu->real("LIVETIME") : 0.0;
+ m_ogip_exposure = (hdu->has_card("EXPOSURE")) ? hdu->real("EXPOSURE") : 0.0;
+
+ // Get OGIP start and stop time
+ std::string tstart = date_obs + "T" + time_obs;
+ std::string tstop = date_end + "T" + time_end;
+ if (tstart != "T") {
+ m_ogip_tstart.utc(tstart);
+ }
+ if (tstop != "T") {
+ m_ogip_tstop.utc(tstop);
+ }
+
+ } // endif: pointer was valid
+
+ // Return
+ return;
+}
+
+
+/***********************************************************************//**
* @brief Write OGIP keywords in FITS HDU
*
* @param[in,out] hdu Pointer to FITS HDU.
diff --git a/src/support/ctobservation.hpp b/src/support/ctobservation.hpp
index ea823e6..797ac92 100755
--- a/src/support/ctobservation.hpp
+++ b/src/support/ctobservation.hpp
@@ -1,7 +1,7 @@
/***************************************************************************
* ctobservation - Base class for observation tools *
* ----------------------------------------------------------------------- *
- * copyright (C) 2016-2017 by Juergen Knoedlseder *
+ * copyright (C) 2016-2018 by Juergen Knoedlseder *
* ----------------------------------------------------------------------- *
* *
* This program is free software: you can redistribute it and/or modify *
@@ -73,6 +73,7 @@ protected:
GCTAObservation* next_unbinned_observation(void);
const GCTAObservation* first_unbinned_observation(void) const;
const GCTAObservation* next_unbinned_observation(void) const;
+ void read_ogip_keywords(GFitsHDU* hdu) const;
void write_ogip_keywords(GFitsHDU* hdu) const;
void set_obs_statistic(const std::string& statistic);
void save_events_fits(void);
diff --git a/src/support/ctool.cpp b/src/support/ctool.cpp
index 119bead..faf9125 100755
--- a/src/support/ctool.cpp
+++ b/src/support/ctool.cpp
@@ -650,7 +650,11 @@ GEbounds ctool::create_ebounds(void)
***************************************************************************/
GSkyMap ctool::create_map(const GObservations& obs)
{
- // Read task parameters
+ // Read coordinate system and projection
+ std::string coordsys = (*this)["coordsys"].string();
+ std::string proj = (*this)["proj"].string();
+
+ // Read sky map reference
double xref = 0.0;
double yref = 0.0;
bool usepnt = (*this)["usepnt"].boolean();
@@ -658,11 +662,11 @@ GSkyMap ctool::create_map(const GObservations& obs)
xref = (*this)["xref"].real();
yref = (*this)["yref"].real();
}
- std::string proj = (*this)["proj"].string();
- std::string coordsys = (*this)["coordsys"].string();
- double binsz = (*this)["binsz"].real();
- int nxpix = (*this)["nxpix"].integer();
- int nypix = (*this)["nypix"].integer();
+
+ // Read sky map pixel size and number of bins
+ double binsz = (*this)["binsz"].real();
+ int nxpix = (*this)["nxpix"].integer();
+ int nypix = (*this)["nypix"].integer();
// If requested, get pointing from observations
if (usepnt) {
diff --git a/test/data/irf_file.fits.gz b/test/data/irf_file.fits.gz
index e9e0aa4..a0a433f 100755
--- a/test/data/irf_file.fits.gz
+++ b/test/data/irf_file.fits.gz
diff --git a/test/test_csobsinfo.py b/test/test_csobsinfo.py
index ed8a752..5391cec 100755
--- a/test/test_csobsinfo.py
+++ b/test/test_csobsinfo.py
@@ -112,7 +112,7 @@ class Test(test):
obsinfo['ra'] = 83.63
obsinfo['dec'] = 22.01
obsinfo['logfile'] = 'csobsinfo_py1.log'
- obsinfo['chatter'] = 2
+ obsinfo['chatter'] = 3
# Run csobsinfo script and save XML file
obsinfo.logFileOpen() # Make sure we get a log file
@@ -122,21 +122,20 @@ class Test(test):
# Check DS9 file
self._check_ds9_file('csobsinfo_py1.reg')
- # Set-up csobsinfo
- obsinfo = cscripts.csobsinfo()
- obsinfo['inobs'] = self._inobs
- obsinfo['outds9file'] = 'csobsinfo_py2.reg'
- obsinfo['offset'] = True
- obsinfo['ra'] = 83.63
- obsinfo['dec'] = 22.01
- obsinfo['logfile'] = 'csobsinfo_py2.log'
- obsinfo['chatter'] = 3
-
- # Execute csobsinfo script
- obsinfo.execute()
-
- # Check DS9 file
- self._check_ds9_file('csobsinfo_py2.reg')
+ # Check result access methods
+ self.test_value(len(obsinfo.zeniths()), 1, 'Test number of zenith angles')
+ self.test_value(obsinfo.zeniths()[0], 0.0, 'Test zenith angle')
+ self.test_value(len(obsinfo.azimuths()), 1, 'Test number of azimuth angles')
+ self.test_value(obsinfo.azimuths()[0], 0.0, 'Test azimuth angle')
+ self.test_value(len(obsinfo.ras()), 1, 'Test number of Right Ascensions')
+ self.test_value(obsinfo.ras()[0], 83.63, 'Test Right Ascension')
+ self.test_value(len(obsinfo.decs()), 1, 'Test number of Declinations')
+ self.test_value(obsinfo.decs()[0], 22.51, 'Test Declination')
+ self.test_value(len(obsinfo.offsets()), 1, 'Test number of offsets')
+ self.test_value(obsinfo.offsets()[0], 0.5, 'Test offset angle')
+ self.test_value(obsinfo.ebounds().size(), 1, 'Test number of energy boundaries')
+ self.test_value(obsinfo.ebounds().emin().GeV(), 100.0, 'Test minimum energy')
+ self.test_value(obsinfo.ebounds().emax().TeV(), 100.0, 'Test maximum energy')
# Return
return
@@ -151,13 +150,12 @@ class Test(test):
# Expect "fk5" in first line
line = f.readline().strip('\n')
- self.test_assert(line == 'fk5',
- 'Check for "fk5" in first line')
+ self.test_value(line, 'fk5',
+ 'Test for "fk5" in first line of DS9 file')
# Expect "fk5" in first line
line = f.readline().strip('\n')
- self.test_assert(line[0:6] == 'point(',
- 'Check for "point(" in second line')
+ self.test_value(line[0:6], 'point(', 'Test for "point(" in second line of DS9 file')
# Close file
f.close()
diff --git a/test/test_ctskymap.py b/test/test_ctskymap.py
index f3d2343..67bbc40 100755
--- a/test/test_ctskymap.py
+++ b/test/test_ctskymap.py
@@ -248,6 +248,7 @@ class Test(test):
skymap['roiradius'] = 0.1
skymap['inradius'] = 0.6
skymap['outradius'] = 0.8
+ skymap['iterations'] = 0
skymap['inexclusion'] = 'NONE'
skymap['caldb'] = self._caldb
skymap['irf'] = self._irf
diff --git a/test/test_obsutils.py b/test/test_obsutils.py
index c0111f0..fd44e89 100755
--- a/test/test_obsutils.py
+++ b/test/test_obsutils.py
@@ -2,7 +2,7 @@
# ==========================================================================
# This scripts performs unit tests for the obsutils module
#
-# Copyright (C) 2017 Juergen Knoedlseder
+# Copyright (C) 2017-2018 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
@@ -214,7 +214,7 @@ class Test(test):
# Check simulation results
self.test_value(res.size(), 1, 'Check number of observations')
self.test_value(res.models().size(), 2, 'Check number of models')
- self.test_value(res.nobserved(), 24, 'Check number of observed events')
+ self.test_value(res.nobserved(), 23, 'Check number of observed events')
self.test_value(res.npred(), 0.0, 'Check number of predicted events')
self.test_value(res[0].eventtype(), 'CountsCube', 'Check event type')
self.test_value(res[0].events().ebounds().emin().TeV(), 1.0,
@@ -223,7 +223,7 @@ class Test(test):
'Check minimum energy')
self.test_value(res[0].events().ebounds().size(), 5,
'Check number of energy bins')
- self.test_value(res[0].events().number(), 24,
+ self.test_value(res[0].events().number(), 23,
'Check number of events in cube')
# Check energy dispersion flag