Karim Sabri / ctools

Commit 0c6e2b1ade48c2cd6642a520d5e1c716550d6f2e

Authored by Jürgen Knödlseder
2 parents 2795ed02 3ed80fe9

Merge branch '2284-csphagen-load-off-regions' into integration 

# Conflicts:
#	ChangeLog
#	NEWS
#	doc/source/admin/release_history/1.5.rst
Inline Side-by-side
Showing 12 changed files with 401 additions and 83 deletions
... ... @@ -15,5 +15,6 @@ Michael Mayer (michael.mayer@physik.hu-berlin.de)
15 15 Stefan Ohm (stefan.ohm81@gmail.com)
16 16 David Sanchez (david.sanchez@lapp.in2p3.fr)
17 17 Hubert Siejkowski (h.siejkowski@cyfronet.pl)
  18 +Andreas Specovius (andreas.specovius@fau.de)
18 19 Luigi Tibaldo (luigi.tibaldo@irap.omp.eu)
19 20 Alexander Ziegler (alexander.ziegler@fau.de)
... ...
ChangeLog
  View file @0c6e2b1
1   -2018-01-09
  1 +2018-01-12
2 2  
3 3 * Version 1.5.0 released
4 4 ========================
5 5  
  6 + Add CUSTOM background method to csphagen script (#2284)
6 7 Add csresspec script (#2289)
7 8 Add ctools::time_reference constant and remove ctools::time_reference() method (#2295)
8 9 Correctly strip .gz suffix from output filenames (#2287)
... ...
1   -New Features and Important Changes in ctools 1.5.0
2   -
3   -9 January 2018
  1 +12 January 2018
4 2  
5 3  
6 4 Introduction
... ... @@ -228,6 +226,11 @@ None
228 226  
229 227 csphagen - Generate On/Off observations
230 228 ---------------------------------------
  229 +The CUSTOM background method was added to allow specifying On and Off regions
  230 +from a region file or FITS WCS map. For multiple observations, a custom Off
  231 +region per observation can be specified in the observation definition XML file
  232 +using the OffRegions parameter (#2284).
  233 +
231 234 The script was added to allow the generation of PHA, ARF and RMF datasets for classical On/Off IACT analysis (#2208).
232 235  
233 236  
... ...
cscripts/Makefile.am
  View file @0c6e2b1
1 1 #############################################################################
2 2 # Makefile for cscripts scripts #
3 3 # ------------------------------------------------------------------------- #
4   -# Copyright (C) 2016-2017 Juergen Knoedlseder #
  4 +# Copyright (C) 2016-2018 Juergen Knoedlseder #
5 5 # ------------------------------------------------------------------------- #
6 6 # #
7 7 # This program is free software: you can redistribute it and/or modify #
... ... @@ -171,6 +171,8 @@ testdata = $(top_srcdir)/test/data/crab.xml \
171 171 $(top_srcdir)/test/data/crab_offaxis1.fits.gz \
172 172 $(top_srcdir)/test/data/crab_offaxis2.fits.gz \
173 173 $(top_srcdir)/test/data/crab_exclusion.fits.gz \
  174 + $(top_srcdir)/test/data/crab_src_reg.reg \
  175 + $(top_srcdir)/test/data/crab_bkg_reg.reg \
174 176 $(top_srcdir)/test/data/tsmap_0.fits \
175 177 $(top_srcdir)/test/data/tsmap_1.fits \
176 178 $(top_srcdir)/test/data/exclusion.reg \
... ...
cscripts/csphagen.par
  View file @0c6e2b1
... ... @@ -26,25 +26,27 @@ prefix, s, h, onoff,,, "Prefix for output files"
26 26 #
27 27 # Script parameters
28 28 #==================
29   -emin, r, a, 0.1,,, "Lower energy limit (TeV)"
30   -emax, r, a, 100.0,,, "Upper energy limit (TeV)"
31   -enumbins, i, a, 120,,, "Number of energy bins"
32   -ebinalg, s, a, LOG,LIN|LOG|FILE,, "Binning algorithm"
33   -ebinfile, f, a, NONE,,, "Name of the file containing the energy bin definition"
34   -srcshape, s, h, CIRCLE,CIRCLE,, "Shape of the source region"
35   -coordsys, s, a, CEL,CEL|GAL,, "Coordinate system (CEL - celestial, GAL - galactic)"
36   -ra, r, a, 83.63,0,360, "Right Ascension of source region centre (deg)"
37   -dec, r, a, 22.01,-90,90, "Declination of source region centre (deg)"
38   -glon, r, a, 184.56,0,360, "Galactic longitude of source region centre (deg)"
39   -glat, r, a, -5.79,-90,90, "Galactic latitude of source region centre (deg)"
40   -rad, r, a, 0.2,0,180, "Radius of source region circle (deg)"
41   -bkgmethod, s, h, REFLECTED,REFLECTED,, "Method for background estimation"
42   -bkgregmin, i, h, 2,1,1000, "Minimum number of background regions"
43   -maxoffset, r, h, 4.0,,, "Maximum offset of source from camera center to accept observation"
44   -stack, b, a, no,,, "Stack multiple observations into single PHA, ARF and RMF files?"
45   -etruemin, r, h, 0.01,,, "Minimum true energy (TeV)"
46   -etruemax, r, h, 300.0,,, "Maximum true energy (TeV)"
47   -etruebins, i, h, 30,,, "Number of bins per decade for true energy bins"
  29 +emin, r, a, 0.1,,, "Lower energy limit (TeV)"
  30 +emax, r, a, 100.0,,, "Upper energy limit (TeV)"
  31 +enumbins, i, a, 120,,, "Number of energy bins"
  32 +ebinalg, s, a, LOG,LIN|LOG|FILE,, "Binning algorithm"
  33 +ebinfile, f, a, NONE,,, "Name of the file containing the energy bin definition"
  34 +srcshape, s, h, CIRCLE,CIRCLE,, "Shape of the source region"
  35 +coordsys, s, a, CEL,CEL|GAL,, "Coordinate system (CEL - celestial, GAL - galactic)"
  36 +ra, r, a, 83.63,0,360, "Right Ascension of source region centre (deg)"
  37 +dec, r, a, 22.01,-90,90, "Declination of source region centre (deg)"
  38 +glon, r, a, 184.56,0,360, "Galactic longitude of source region centre (deg)"
  39 +glat, r, a, -5.79,-90,90, "Galactic latitude of source region centre (deg)"
  40 +rad, r, a, 0.2,0,180, "Radius of source region circle (deg)"
  41 +srcregfile, f, a, on.reg,,, "Source region file (ds9 or FITS WCS)"
  42 +bkgmethod, s, a, REFLECTED,REFLECTED|CUSTOM,, "Method for background estimation"
  43 +bkgregfile, f, a, bkg.reg,,, "Background regions file (ds9 or FITS WCS)"
  44 +bkgregmin, i, h, 2,1,1000, "Minimum number of background regions"
  45 +maxoffset, r, h, 4.0,,, "Maximum offset of source from camera center to accept observation"
  46 +stack, b, a, no,,, "Stack multiple observations into single PHA, ARF and RMF files?"
  47 +etruemin, r, h, 0.01,,, "Minimum true energy (TeV)"
  48 +etruemax, r, h, 300.0,,, "Maximum true energy (TeV)"
  49 +etruebins, i, h, 30,,, "Number of bins per decade for true energy bins"
48 50  
49 51 #
50 52 # Standard parameters
... ...
cscripts/csphagen.py
  View file @0c6e2b1
... ... @@ -3,7 +3,7 @@
3 3 # Computes the PHA spectra for source/background and ARF/RMF files using the
4 4 # reflected region method
5 5 #
6   -# Copyright (C) 2017 Luigi Tibaldo
  6 +# Copyright (C) 2017-2018 Luigi Tibaldo
7 7 #
8 8 # This program is free software: you can redistribute it and/or modify
9 9 # it under the terms of the GNU General Public License as published by
... ... @@ -50,25 +50,19 @@ class csphagen(ctools.csobservation):
50 50 self._has_exclusion = False
51 51 self._srcshape = ''
52 52 self._rad = 0.0
  53 + self._bkg_reg_files = {} # dictionary: bkg region file per observation {obsid:filepath, ..}
53 54  
54 55 # Return
55 56 return
56 57  
  58 +
57 59 # Private methods
58   - def _get_parameters(self):
  60 + def _query_src_direction(self):
59 61 """
60   - Get parameters from parfile and setup observations
  62 + Set up the source direction parameter.
  63 + Query relevant parameters.
61 64 """
62   - # Setup observations (require response and allow event list, don't
63   - # allow counts cube)
64   - self._setup_observations(self.obs(), True, True, False)
65   -
66   - # Set energy bounds
67   - self._ebounds = self._create_ebounds()
68   -
69   - # Initialise source position/region querying relevant parameters
70 65 self._src_dir = gammalib.GSkyDir()
71   - self._src_reg = gammalib.GSkyRegions()
72 66 coordsys = self['coordsys'].string()
73 67 if coordsys == 'CEL':
74 68 ra = self['ra'].real()
... ... @@ -78,17 +72,139 @@ class csphagen(ctools.csobservation):
78 72 glon = self['glon'].real()
79 73 glat = self['glat'].real()
80 74 self._src_dir.lb_deg(glon, glat)
  75 +
  76 + def _get_parameters_bkgmethod_reflected(self):
  77 + """
  78 + Get parameters for REFLECTED background method
  79 + """
  80 + # Get source shape
81 81 self._srcshape = self['srcshape'].string()
  82 +
  83 + # Set source region (so far only CIRCLE is supported)
82 84 if self._srcshape == 'CIRCLE':
  85 +
  86 + # Query source direction
  87 + self._query_src_direction()
  88 +
  89 + # Query minimum number of background regions
  90 + self['bkgregmin'].integer()
  91 +
  92 + # Set circular source region
83 93 self._rad = self['rad'].real()
84 94 self._src_reg.append(gammalib.GSkyRegionCircle(self._src_dir, self._rad))
85 95  
86   - # Query background estimation method and parameters
  96 + # Return
  97 + return
  98 +
  99 + def _get_parameters_bkgmethod_custom(self):
  100 + """
  101 + Get parameters for CUSTOM background method
  102 +
  103 + Raises
  104 + ------
  105 + RuntimeError
  106 + Multiple On-Regions in on region file or missing Off-Regions files.
  107 + """
  108 + # Set up source region
  109 + src_reg_file = self['srcregfile'].filename()
  110 + self._src_reg = gammalib.GSkyRegions(src_reg_file)
  111 +
  112 + # Raise an exception if there is more than one source region
  113 + if len(self._src_reg) != 1:
  114 + raise RuntimeError('Only 1 On region is allowed.')
  115 +
  116 + # Set up source direction. Query parameters if neccessary.
  117 + if isinstance(self._src_reg[0], gammalib.GSkyRegionCircle):
  118 + self._src_dir = self._src_reg[0].centre()
  119 + self._rad = self._src_reg[0].radius()
  120 + else:
  121 + self._query_src_direction()
  122 +
  123 + # Try opening file as XML file. If this files we have probably a FITS
  124 + # events file on input, hence we query later one background region
  125 + # file
  126 + try:
  127 + # Open observation definition XML file
  128 + obs_file = self['inobs'].filename()
  129 + xml = gammalib.GXml(obs_file)
  130 +
  131 + # Go through elements in <observationlist>
  132 + for elem in xml[0]:
  133 +
  134 + # Get observation id from <observation> header
  135 + obs_id = elem.attribute('id')
  136 +
  137 + # Get region file path from <observation> body
  138 + off_regions_file = ''
  139 + for entry in elem:
  140 +
  141 + # <parameter name="OffRegions" file=".."/> contains region file
  142 + if entry.attribute('name') == 'OffRegions':
  143 + off_regions_file = entry.attribute('file')
  144 + break
  145 +
  146 + # Store file paths
  147 + self._bkg_reg_files[obs_id] = gammalib.GFilename(off_regions_file)
  148 +
  149 + # Not existing background region files are set to None (for every obs)
  150 + if not self._bkg_reg_files[obs_id].exists():
  151 + self._bkg_reg_files[obs_id] = None
  152 +
  153 + # XML file opening failed
  154 + except:
  155 + self._bkg_reg_files = {0: None}
  156 +
  157 + # Query missing background region file
  158 + if len([True for x in self._bkg_reg_files if self._bkg_reg_files[x] is None]) > 0:
  159 + if len(self._bkg_reg_files) > 1:
  160 + # Querying of multiple regions not supported
  161 + raise RuntimeError('There are missing bkg region files in the '
  162 + 'observation definition XML file. Please '
  163 + 'specify all background region files.')
  164 + else:
  165 + # Query bkg region file for single observation
  166 + self._bkg_reg_files[self._bkg_reg_files.keys()[0]] = \
  167 + self['bkgregfile'].filename()
  168 +
  169 + # Return
  170 + return
  171 +
  172 + def _get_parameters_bkgmethod(self):
  173 + """
  174 + Get background method parameters
  175 + """
  176 + # Get background method
87 177 bkgmethod = self['bkgmethod'].string()
  178 +
  179 + # Get background method dependent parameters
88 180 if bkgmethod == 'REFLECTED':
89   - self['bkgregmin'].integer()
  181 + self._get_parameters_bkgmethod_reflected()
  182 + elif bkgmethod == 'CUSTOM':
  183 + self._get_parameters_bkgmethod_custom()
  184 +
  185 + # Query parameters that are needed for all background methods
90 186 self['maxoffset'].real()
91 187  
  188 + # Return
  189 + return
  190 +
  191 + def _get_parameters(self):
  192 + """
  193 + Get parameters from parfile and setup observations
  194 + """
  195 + # Setup observations (require response and allow event list, don't
  196 + # allow counts cube)
  197 + self._setup_observations(self.obs(), True, True, False)
  198 +
  199 + # Set energy bounds
  200 + self._ebounds = self._create_ebounds()
  201 +
  202 + # Initialize empty src regions container
  203 + self._src_reg = gammalib.GSkyRegions()
  204 +
  205 + # Get background method parameters
  206 + self._get_parameters_bkgmethod()
  207 +
92 208 # Exclusion map
93 209 if self['inexclusion'].filename().is_fits():
94 210 self._excl_reg = gammalib.GSkyRegionMap(self['inexclusion'].filename())
... ... @@ -179,6 +295,77 @@ class csphagen(ctools.csobservation):
179 295 # Return reflected regions
180 296 return regions
181 297  
  298 + def _regions_from_file(self, obs, file_path):
  299 + """
  300 + Read list of reflected regions for a single observation (pointing)
  301 + from file.
  302 +
  303 + Perform basic tests for circle shaped regions.
  304 +
  305 + Parameters
  306 + ----------
  307 + obs : `~gammalib.GCTAObservation()`
  308 + CTA observation
  309 + file_path: `~gammalib.GFilename()`
  310 + File path to file containing off regions.
  311 +
  312 + Returns
  313 + -------
  314 + regions : `~gammalib.GSkyRegions`
  315 + List of reflected regions
  316 + """
  317 + # Load regions from ds9 file of FITS WCS map
  318 + regions = gammalib.GSkyRegions( file_path )
  319 +
  320 + # Do region tests only for circle shaped on-region, simply return list
  321 + # for other shapes.
  322 + on_region = self._src_reg[0]
  323 + if isinstance(on_region, gammalib.GSkyRegionCircle):
  324 +
  325 + # Prepare parameters for comparison
  326 + pnt_dir = obs.pointing().dir()
  327 + offset_on = pnt_dir.dist_deg(self._src_dir)
  328 + radius_on = on_region.radius()
  329 +
  330 + # Do some tests
  331 + for region in regions:
  332 +
  333 + # Do region tests only for circle shaped off-regions
  334 + if isinstance(region, gammalib.GSkyRegionCircle):
  335 +
  336 + # Gather parameters
  337 + offset_off = pnt_dir.dist_deg(region.centre())
  338 + radius_off = region.radius()
  339 +
  340 + # Warning: different radius
  341 + if radius_off != radius_on:
  342 + self._log_string(gammalib.EXPLICIT, 'Off region radius '
  343 + 'differs from On region radius: %.2f '
  344 + 'vs %.2f' % ((radius_off, radius_on)))
  345 +
  346 + # Warning: offset is too small or too large
  347 + if offset_off <= radius_on or offset_off >= self['maxoffset'].real():
  348 + self._log_string(gammalib.EXPLICIT, 'Region offset '
  349 + 'invalid for observation %s '
  350 + '(%.3f deg) ' % ((obs.id(), offset_off)))
  351 +
  352 + # Warning: offset of off and on region differ
  353 + if offset_off != offset_on:
  354 + self._log_string(gammalib.EXPLICIT, 'Off region offset '
  355 + 'differs %.3f deg from On region '
  356 + 'offset' % (offset_off-offset_on))
  357 +
  358 + # end: off-regions test for circle shaped regions
  359 + # end: for loop iterating over off regions
  360 + # end: do region tests only if on-region is circle shaped
  361 +
  362 + # Inform user how many regions were loaded
  363 + self._log_string(gammalib.NORMAL, 'Read %i regions from file %s' % \
  364 + (regions.size(), file_path.url()))
  365 +
  366 + # return list of regions
  367 + return regions
  368 +
182 369 def _set_models(self, obs):
183 370 """
184 371 Set models in observation container
... ... @@ -247,6 +434,45 @@ class csphagen(ctools.csobservation):
247 434 # Return energy boundaries
248 435 return ebounds
249 436  
  437 + def _set_background_regions(self, obs):
  438 + """
  439 + Set background regions for an observation
  440 +
  441 + Parameters
  442 + ----------
  443 + obs : `~gammalib.GCTAObservation()`
  444 + CTA observation
  445 +
  446 + Returns
  447 + -------
  448 + regions : `~gammalib.GSkyRegions()`
  449 + Background regions
  450 + """
  451 + # Initialise empty background regions for this observation
  452 + bkg_reg = gammalib.GSkyRegions()
  453 +
  454 + # If reflected background is requested then created reflected
  455 + # background regions
  456 + if self['bkgmethod'].string() == 'REFLECTED':
  457 + bkg_reg = self._reflected_regions(obs)
  458 +
  459 + # ... otherwise if custom background is requested the get the background
  460 + # regions from a region file
  461 + elif self['bkgmethod'].string() == 'CUSTOM':
  462 +
  463 + # Get file path of background region file from dictionary
  464 + if len(self._bkg_reg_files) > 1:
  465 + bkg_reg_file = self._bkg_reg_files[obs.id()]
  466 + else:
  467 + # for inobs==event list no obs id is set.
  468 + bkg_reg_file = self._bkg_reg_files[self._bkg_reg_files.keys()[0]]
  469 +
  470 + # Load background regions from region file for current observation
  471 + bkg_reg = self._regions_from_file(obs, bkg_reg_file)
  472 +
  473 + # Return background regions
  474 + return bkg_reg
  475 +
250 476  
251 477 # Public methods
252 478 def run(self):
... ... @@ -286,15 +512,13 @@ class csphagen(ctools.csobservation):
286 512 # Loop through observations and generate pha, arf, rmf files
287 513 for obs in self.obs():
288 514  
289   - # Initialise background regions for this observation
290   - bkg_reg = gammalib.GSkyRegions()
  515 + # Log current observation
  516 + self._log_string(gammalib.NORMAL, 'Process Observation %s' % (obs.id()))
291 517  
292   - # If reflected background is requested then created reflected
293   - # background regions
294   - if self['bkgmethod'].string() == 'REFLECTED':
295   - bkg_reg = self._reflected_regions(obs)
  518 + # Set background regions for this observation
  519 + bkg_reg = self._set_background_regions(obs)
296 520  
297   - # If there are reflected regions then create On/Off observation
  521 + # If there are background regions then create On/Off observation
298 522 # and append it to the output container
299 523 if bkg_reg.size() >= self['bkgregmin'].integer():
300 524 onoff = gammalib.GCTAOnOffObservation(obs, self._src_dir,
... ...
cscripts/obsutils.py
  View file @0c6e2b1
... ... @@ -130,6 +130,7 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0,
130 130 phagen['rad'] = onrad
131 131 phagen['stack'] = False
132 132 phagen['inexclusion'] = 'NONE'
  133 + phagen['bkgmethod'] = 'REFLECTED'
133 134  
134 135 # Optionally open the log file
135 136 if log:
... ...
doc/source/admin/release_history/1.5.rst
  View file @0c6e2b1
... ... @@ -29,6 +29,8 @@ Improvements
29 29 * [`2295 <https://cta-redmine.irap.omp.eu/issues/2295>`_] -
30 30 Add ``ctools::time_reference`` constant and remove ``ctools::time_reference()``
31 31 method
  32 +* [`2284 <https://cta-redmine.irap.omp.eu/issues/2284>`_] -
  33 + Add ``CUSTOM`` background method to :ref:`csphagen` script
32 34 * [`2273 <https://cta-redmine.irap.omp.eu/issues/2273>`_] -
33 35 Support On/Off analysis in :ref:`cslightcrv` and :ref:`csphasecrv`
34 36 * [`2267 <https://cta-redmine.irap.omp.eu/issues/2267>`_] -
... ...
doc/source/users/reference_manual/csphagen.rst
  View file @0c6e2b1
... ... @@ -13,8 +13,9 @@ Synopsis
13 13 This script can be used to derive from an observation or a set of observations
14 14 the count spectra in a source region and in background regions, as well as the
15 15 detector response (effective area, energy redistribution matrix), that can be
16   -used to perform a classical 1D spectral analysis. The output files are saved in
17   -the OGIP format normally used in X-ray astronomy (PHA, ARF, RMF).
  16 +used to perform a classical 1D spectral analysis. Regions can be placed
  17 +automatically or by hand. The output files are saved in the OGIP format normally
  18 +used in X-ray astronomy (PHA, ARF, RMF).
18 19 `See here <https://heasarc.gsfc.nasa.gov/docs/heasarc/ofwg/docs/spectra/ogip_92_007/node5.html>`__.
19 20  
20 21 The outputs are:
... ... @@ -30,50 +31,48 @@ General parameters
30 31 ------------------
31 32  
32 33 ``inobs [file]``
33   - Input event list or observation definition XML file.
  34 + Input event list or observation definition XML file
34 35  
35 36 ``caldb [string]``
36   - Calibration database.
  37 + Calibration database
37 38  
38 39 ``irf [string]``
39   - Instrument response function.
  40 + Instrument response function
40 41  
41 42 ``(inexclusion = NONE) [file]``
42 43 Optional FITS file containing a WCS map in the first hdu that defines sky
43   - regions not to be used for background estimation (where map value != 0).
  44 + regions not to be used for background estimation (where map value != 0)
44 45  
45 46 ``outobs [string]``
46   - Output observation definition XML file.
  47 + Output observation definition XML file
47 48  
48 49 ``(prefix = onoff) [string]``
49   - Prefix of the file name for output PHA, ARF, RMF, XML, and DS9 region files.
  50 + Prefix of the file name for output PHA, ARF, RMF, XML, and DS9 region files
50 51  
51 52 ``emin [real]``
52   - Lower energy limit (TeV) if ``LIN`` or ``LOG`` binning algorithms are used.
  53 + Lower energy limit (TeV) if ``LIN`` or ``LOG`` binning algorithms are used
53 54  
54 55 ``emax [real]``
55   - Upper energy limit (TeV) if ``LIN`` or ``LOG`` binning algorithms are used.
  56 + Upper energy limit (TeV) if ``LIN`` or ``LOG`` binning algorithms are used
56 57  
57 58 ``enumbins [integer]``
58 59 Number of energy bins. At least 30 bins per decade are recommended for
59 60 proper evaluation of the instrument response.
60 61  
61   -``ebinalg [string]``
62   - Energy binning algorithm:
63   - ``LOG`` : logarithmically spaced energy bins;
64   - ``LIN`` : linearly spaced energy bins;
65   - ``FILE``: energy bounds retrieved from file.
  62 +``ebinalg <FILE|LIN|LOG> [string]``
  63 + Algorithm for defining energy bins (``FILE``: energy bounds retrieved from
  64 + file, see ``ebinfile`` parameter; ``LIN``: linearly spaced energy bins;
  65 + ``LOG``: logarithmically spaced energy bins)
66 66  
67 67 ``ebinfile [file]``
68 68 Name of the file containing the energy bin definition if ``FILE`` algorithm
69   - is used.
  69 + is used
70 70  
71   -``(srcshape = CIRCLE) [string]``
72   - Shape of the source region.
73   - ``CIRCLE``: circular region around given position.
  71 +``(srcshape = CIRCLE) <CIRCLE> [string]``
  72 + Shape of the source region (``CIRCLE``: circular region around given position)
74 73  
75   -``coordsys [string]``
76   - Coordinate system (``CEL`` - celestial, ``GAL`` - galactic)
  74 +``coordsys <CEL|GAL> [string]``
  75 + Coordinate system (``CEL``: celestial; ``GAL``: galactic)
77 76  
78 77 ``ra [real]``
79 78 Right Ascension of source region centre (deg)
... ... @@ -90,10 +89,25 @@ General parameters
90 89 ``rad [real]``
91 90 Radius of source region circle (deg)
92 91  
93   -``(bkgmethod = REFLECTED) [string]``
94   - Method for background estimation.
95   - ``REFLECTED:`` background evaluated in regions with the same shape as
96   - source region reflected w.r.t. pointing direction for each observation.
  92 +``srcregfile [file]``
  93 + Source region file (ds9 or FITS WCS map)
  94 +
  95 +``bkgmethod <REFLECTED|CUSTOM> [string]``
  96 + Method for background estimation:
  97 +
  98 + - ``REFLECTED``: background is evaluated in regions with the same shape as
  99 + the source region reflected w.r.t. pointing direction for each observation
  100 +
  101 + - ``CUSTOM``: background is evaluated in regions specified by user. For an
  102 + event list or a single observation in the observation definition XML file
  103 + a region file will be queried (see ``bkgregfile`` parameter). For multiple
  104 + observations specified in the observation definition XML file the name of
  105 + the region file will be extracted from the ``OffRegions`` parameter that
  106 + needs to be specified for each observation in the observation definition
  107 + XML file. Off region files can be either ds9 region files or FITS WCS maps.
  108 +
  109 +``bkgregfile [file]``
  110 + Background regions file (ds9 or FITS WCS map)
97 111  
98 112 ``(bkgregmin = 2) [integer]``
99 113 Minimum number of background regions that are required for an observation.
... ... @@ -102,20 +116,20 @@ General parameters
102 116  
103 117 ``(maxoffset = 4.0) [real]``
104 118 Maximum offset in degrees of source from camera center to accept the
105   - observation.
  119 + observation
106 120  
107 121 ``stack [boolean]``
108 122 Specifies whether multiple observations should be stacked (``yes``) or
109   - whether run-wise PHA, ARF and RMF files should be produced (``no``).
  123 + whether run-wise PHA, ARF and RMF files should be produced (``no``)
110 124  
111 125 ``(etruemin = 0.01) [real]``
112   - Minimum true energy (TeV).
  126 + Minimum true energy (TeV)
113 127  
114 128 ``(etruemax = 0.01) [real]``
115   - Maximum true energy (TeV).
  129 + Maximum true energy (TeV)
116 130  
117 131 ``(etruebins = 30) [integer]``
118   - Number of bins per decade for true energy bins.
  132 + Number of bins per decade for true energy bins
119 133  
120 134  
121 135 Standard parameters
... ...
test/data/crab_bkg_reg.reg 0 → 100644
  View file @0c6e2b1
  1 +# Region file format: DS9 version 4.1
  2 +global color=green dashlist=8 3 width=1
  3 +fk5;circle(82.2412081162728,22.7975361111674,0.2)
  4 +fk5;circle(82.5360705305998,23.1293202741537,0.2)
  5 +fk5;circle(82.980058098929,23.260640435668,0.2)
  6 +fk5;circle(83.4304728956696,23.1492236356079,0.2)
  7 +fk5;circle(83.7424610541293,22.8309432073396,0.2)
... ...
test/data/crab_src_reg.reg 0 → 100644
  View file @0c6e2b1
  1 +# Region file format: DS9 version 4.1
  2 +global color=green dashlist=8 3 width=1
  3 +fk5;circle(83.633,22.0145,0.2)
... ...
test/test_csphagen.py
  View file @0c6e2b1
... ... @@ -48,10 +48,11 @@ class Test(test):
48 48 self._exclusion = self._datadir + '/crab_exclusion.fits'
49 49 self._nreg_with_excl = 5
50 50 self._nreg_wo_excl = 8
  51 + self._nreg_bkg_reg = 5
51 52 self._nreg_mul = [self._nreg_with_excl, 6]
52 53 self._nbins = 10
53   - # number of expected background regions with/wo exclusion,
54   - # and for two different runs
  54 + self._regfile_src = self._datadir + '/crab_src_reg.reg'
  55 + self._regfile_bkg = self._datadir + '/crab_bkg_reg.reg'
55 56  
56 57 # Return
57 58 return
... ... @@ -88,7 +89,8 @@ class Test(test):
88 89 ' coordsys="CEL" ra=83.633 dec=22.0145' + \
89 90 ' rad=0.2 stack="no" inexclusion="' + \
90 91 self._exclusion + \
91   - '" outobs="csphagen_cmd1.xml" prefix="csphagen_cmd1" ' + \
  92 + '" bkgmethod="REFLECTED" '+ \
  93 + 'outobs="csphagen_cmd1.xml" prefix="csphagen_cmd1" ' + \
92 94 'logfile="csphagen_cmd1.log" chatter=2'
93 95  
94 96 # Check if execution of wrong command fails
... ... @@ -112,7 +114,8 @@ class Test(test):
112 114 ' coordsys="CEL" ra=83.633 dec=22.0145' + \
113 115 ' rad=0.2 stack="no" inexclusion="' + \
114 116 self._exclusion + \
115   - '" outobs="csphagen_cmd2.xml" prefix="csphagen_cmd2" ' + \
  117 + '" bkgmethod="REFLECTED" '+ \
  118 + 'outobs="csphagen_cmd2.xml" prefix="csphagen_cmd2" ' + \
116 119 'logfile="csphagen_cmd2.log" debug=yes chatter=1'
117 120  
118 121 # Check if execution failed
... ... @@ -144,6 +147,7 @@ class Test(test):
144 147 phagen['rad'] = 0.2
145 148 phagen['stack'] = False
146 149 phagen['inexclusion'] = self._exclusion
  150 + phagen['bkgmethod'] = 'REFLECTED'
147 151 phagen['etruemin'] = 0.05
148 152 phagen['etruemax'] = 150.0
149 153 phagen['etruebins'] = 5
... ... @@ -152,7 +156,7 @@ class Test(test):
152 156 phagen['logfile'] = 'csphagen_py1.log'
153 157 phagen['chatter'] = 1
154 158  
155   - # Run script
  159 + # Execute script
156 160 phagen.execute()
157 161  
158 162 # Check output
... ... @@ -173,6 +177,7 @@ class Test(test):
173 177 phagen['dec'] = 22.0145
174 178 phagen['rad'] = 0.2
175 179 phagen['stack'] = False
  180 + phagen['bkgmethod'] = 'REFLECTED'
176 181 phagen['etruemin'] = 0.05
177 182 phagen['etruemax'] = 150.0
178 183 phagen['etruebins'] = 5
... ... @@ -181,7 +186,7 @@ class Test(test):
181 186 phagen['logfile'] = 'csphagen_py2.log'
182 187 phagen['chatter'] = 2
183 188  
184   - # Run script
  189 + # Execute script
185 190 phagen.execute()
186 191  
187 192 # Check output
... ... @@ -209,6 +214,7 @@ class Test(test):
209 214 phagen['rad'] = 0.2
210 215 phagen['stack'] = False
211 216 phagen['inexclusion'] = self._exclusion
  217 + phagen['bkgmethod'] = 'REFLECTED'
212 218 phagen['etruemin'] = 0.05
213 219 phagen['etruemax'] = 150.0
214 220 phagen['etruebins'] = 5
... ... @@ -239,6 +245,7 @@ class Test(test):
239 245 phagen['rad'] = 0.2
240 246 phagen['stack'] = True
241 247 phagen['inexclusion'] = self._exclusion
  248 + phagen['bkgmethod'] = 'REFLECTED'
242 249 phagen['etruemin'] = 0.05
243 250 phagen['etruemax'] = 150.0
244 251 phagen['etruebins'] = 5
... ... @@ -247,7 +254,7 @@ class Test(test):
247 254 phagen['logfile'] = 'csphagen_py4.log'
248 255 phagen['chatter'] = 4
249 256  
250   - # Run script
  257 + # Execute script
251 258 phagen.execute()
252 259  
253 260 # Check output
... ... @@ -256,6 +263,34 @@ class Test(test):
256 263 0, check_regions=False)
257 264 self._check_outobs('csphagen_py4', 1)
258 265  
  266 + # Setup csphagen for test with custom On and Off regions provided
  267 + phagen = cscripts.csphagen()
  268 + phagen['inobs'] = self._myevents1
  269 + phagen['caldb'] = self._caldb
  270 + phagen['irf'] = self._irf
  271 + phagen['ebinalg'] = 'LOG'
  272 + phagen['emin'] = 0.1
  273 + phagen['emax'] = 100.0
  274 + phagen['enumbins'] = self._nbins
  275 + phagen['bkgmethod'] = 'CUSTOM'
  276 + phagen['srcregfile'] = self._regfile_src
  277 + phagen['bkgregfile'] = self._regfile_bkg
  278 + phagen['etruemin'] = 0.05
  279 + phagen['etruemax'] = 150.0
  280 + phagen['etruebins'] = 5
  281 + phagen['stack'] = False
  282 + phagen['outobs'] = 'csphagen_py5.xml'
  283 + phagen['prefix'] = 'csphagen_py5'
  284 + phagen['logfile'] = 'csphagen_py5.log'
  285 + phagen['chatter'] = 2
  286 +
  287 + # Execute script
  288 + phagen.execute()
  289 +
  290 + # Check output
  291 + self._check_output('csphagen_py5', self._nbins, self._nreg_bkg_reg)
  292 + self._check_outobs('csphagen_py5', 1)
  293 +
259 294 # Return
260 295 return
261 296  
... ... @@ -309,7 +344,7 @@ class Test(test):
309 344 self.test_assert(table.contains(col),
310 345 'FITS file contains "' + col + '" column')
311 346  
312   - # Check SPECTRUM table
  347 + # Check EBOUNDS table
313 348 table = fits['EBOUNDS']
314 349 self._check_ebounds(table, bins)
315 350  
... ... @@ -433,3 +468,26 @@ class Test(test):
433 468  
434 469 # Return
435 470 return
  471 +
  472 + def _read_pha_counts(self, filename):
  473 + """
  474 + Read and integrate the counts in a pha file.
  475 + Pha file structure already tested in _check_pha()
  476 + """
  477 + # Open FITS file
  478 + fits = gammalib.GFits(filename)
  479 +
  480 + # Get SPECTRUM table
  481 + table = fits['SPECTRUM']
  482 +
  483 + # Integrate counts
  484 + counts_col = table['COUNTS']
  485 + counts = 0
  486 + for channel in range( counts_col.nrows() ):
  487 + counts += counts_col[channel]
  488 +
  489 + # Close FITS file
  490 + fits.close()
  491 +
  492 + # Return
  493 + return counts
... ...