diff --git a/AUTHORS b/AUTHORS index 42793a0..3a059d3 100644 --- a/AUTHORS +++ b/AUTHORS @@ -15,5 +15,6 @@ Michael Mayer (michael.mayer@physik.hu-berlin.de) Stefan Ohm (stefan.ohm81@gmail.com) David Sanchez (david.sanchez@lapp.in2p3.fr) Hubert Siejkowski (h.siejkowski@cyfronet.pl) +Andreas Specovius (andreas.specovius@fau.de) Luigi Tibaldo (luigi.tibaldo@irap.omp.eu) Alexander Ziegler (alexander.ziegler@fau.de) diff --git a/ChangeLog b/ChangeLog index 1c73b6a..8dd8a12 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,9 @@ -2018-01-09 +2018-01-12 * Version 1.5.0 released ======================== + Add CUSTOM background method to csphagen script (#2284) Add csresspec script (#2289) Add ctools::time_reference constant and remove ctools::time_reference() method (#2295) Correctly strip .gz suffix from output filenames (#2287) diff --git a/NEWS b/NEWS index 904941b..5f1f60d 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,4 @@ -New Features and Important Changes in ctools 1.5.0 - -9 January 2018 +12 January 2018 Introduction @@ -228,6 +226,11 @@ 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). diff --git a/cscripts/Makefile.am b/cscripts/Makefile.am index 0049fdc..94fc226 100644 --- a/cscripts/Makefile.am +++ b/cscripts/Makefile.am @@ -1,7 +1,7 @@ ############################################################################# # Makefile for cscripts scripts # # ------------------------------------------------------------------------- # -# Copyright (C) 2016-2017 Juergen Knoedlseder # +# Copyright (C) 2016-2018 Juergen Knoedlseder # # ------------------------------------------------------------------------- # # # # This program is free software: you can redistribute it and/or modify # @@ -172,6 +172,8 @@ testdata = $(top_srcdir)/test/data/crab.xml \ $(top_srcdir)/test/data/crab_offaxis1.fits.gz \ $(top_srcdir)/test/data/crab_offaxis2.fits.gz \ $(top_srcdir)/test/data/crab_exclusion.fits.gz \ + $(top_srcdir)/test/data/crab_src_reg.reg \ + $(top_srcdir)/test/data/crab_bkg_reg.reg \ $(top_srcdir)/test/data/tsmap_0.fits \ $(top_srcdir)/test/data/tsmap_1.fits \ $(top_srcdir)/test/data/exclusion.reg \ diff --git a/cscripts/csphagen.par b/cscripts/csphagen.par index 45fc441..9366518 100644 --- a/cscripts/csphagen.par +++ b/cscripts/csphagen.par @@ -26,25 +26,27 @@ prefix, s, h, onoff,,, "Prefix for output files" # # Script parameters #================== -emin, r, a, 0.1,,, "Lower energy limit (TeV)" -emax, r, a, 100.0,,, "Upper energy limit (TeV)" -enumbins, i, a, 120,,, "Number of energy bins" -ebinalg, s, a, LOG,LIN|LOG|FILE,, "Binning algorithm" -ebinfile, f, a, NONE,,, "Name of the file containing the energy bin definition" -srcshape, s, h, CIRCLE,CIRCLE,, "Shape of the source region" -coordsys, s, a, CEL,CEL|GAL,, "Coordinate system (CEL - celestial, GAL - galactic)" -ra, r, a, 83.63,0,360, "Right Ascension of source region centre (deg)" -dec, r, a, 22.01,-90,90, "Declination of source region centre (deg)" -glon, r, a, 184.56,0,360, "Galactic longitude of source region centre (deg)" -glat, r, a, -5.79,-90,90, "Galactic latitude of source region centre (deg)" -rad, r, a, 0.2,0,180, "Radius of source region circle (deg)" -bkgmethod, s, h, REFLECTED,REFLECTED,, "Method for background estimation" -bkgregmin, i, h, 2,1,1000, "Minimum number of background regions" -maxoffset, r, h, 4.0,,, "Maximum offset of source from camera center to accept observation" -stack, b, a, no,,, "Stack multiple observations into single PHA, ARF and RMF files?" -etruemin, r, h, 0.01,,, "Minimum true energy (TeV)" -etruemax, r, h, 300.0,,, "Maximum true energy (TeV)" -etruebins, i, h, 30,,, "Number of bins per decade for true energy bins" +emin, r, a, 0.1,,, "Lower energy limit (TeV)" +emax, r, a, 100.0,,, "Upper energy limit (TeV)" +enumbins, i, a, 120,,, "Number of energy bins" +ebinalg, s, a, LOG,LIN|LOG|FILE,, "Binning algorithm" +ebinfile, f, a, NONE,,, "Name of the file containing the energy bin definition" +srcshape, s, h, CIRCLE,CIRCLE,, "Shape of the source region" +coordsys, s, a, CEL,CEL|GAL,, "Coordinate system (CEL - celestial, GAL - galactic)" +ra, r, a, 83.63,0,360, "Right Ascension of source region centre (deg)" +dec, r, a, 22.01,-90,90, "Declination of source region centre (deg)" +glon, r, a, 184.56,0,360, "Galactic longitude of source region centre (deg)" +glat, r, a, -5.79,-90,90, "Galactic latitude of source region centre (deg)" +rad, r, a, 0.2,0,180, "Radius of source region circle (deg)" +srcregfile, f, a, on.reg,,, "Source region file (ds9 or FITS WCS)" +bkgmethod, s, a, REFLECTED,REFLECTED|CUSTOM,, "Method for background estimation" +bkgregfile, f, a, bkg.reg,,, "Background regions file (ds9 or FITS WCS)" +bkgregmin, i, h, 2,1,1000, "Minimum number of background regions" +maxoffset, r, h, 4.0,,, "Maximum offset of source from camera center to accept observation" +stack, b, a, no,,, "Stack multiple observations into single PHA, ARF and RMF files?" +etruemin, r, h, 0.01,,, "Minimum true energy (TeV)" +etruemax, r, h, 300.0,,, "Maximum true energy (TeV)" +etruebins, i, h, 30,,, "Number of bins per decade for true energy bins" # # Standard parameters diff --git a/cscripts/csphagen.py b/cscripts/csphagen.py index a99a515..48cdb49 100755 --- a/cscripts/csphagen.py +++ b/cscripts/csphagen.py @@ -3,7 +3,7 @@ # Computes the PHA spectra for source/background and ARF/RMF files using the # reflected region method # -# Copyright (C) 2017 Luigi Tibaldo +# Copyright (C) 2017-2018 Luigi Tibaldo # # 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 @@ -54,21 +54,14 @@ class csphagen(ctools.csobservation): # Return return + # Private methods - def _get_parameters(self): + def _query_src_direction(self): """ - Get parameters from parfile and setup observations + Set up the source direction parameter. + Query relevant parameters. """ - # Setup observations (require response and allow event list, don't - # allow counts cube) - self._setup_observations(self.obs(), True, True, False) - - # Set energy bounds - self._ebounds = self._create_ebounds() - - # Initialise source position/region querying relevant parameters self._src_dir = gammalib.GSkyDir() - self._src_reg = gammalib.GSkyRegions() coordsys = self['coordsys'].string() if coordsys == 'CEL': ra = self['ra'].real() @@ -78,17 +71,100 @@ class csphagen(ctools.csobservation): glon = self['glon'].real() glat = self['glat'].real() self._src_dir.lb_deg(glon, glat) + + def _get_parameters_bkgmethod_reflected(self): + """ + Get parameters for REFLECTED background method + """ + # Get source shape self._srcshape = self['srcshape'].string() + + # Set source region (so far only CIRCLE is supported) if self._srcshape == 'CIRCLE': + + # Query source direction + self._query_src_direction() + + # Query minimum number of background regions + self['bkgregmin'].integer() + + # Set circular source region self._rad = self['rad'].real() self._src_reg.append(gammalib.GSkyRegionCircle(self._src_dir, self._rad)) - # Query background estimation method and parameters + # Return + return + + def _get_parameters_bkgmethod_custom(self): + """ + Get parameters for CUSTOM background method + + Raises + ------ + RuntimeError + Only one On region is allowed + """ + # Set up source region + filename = self['srcregfile'].filename() + self._src_reg = gammalib.GSkyRegions(filename) + + # Raise an exception if there is more than one source region + if len(self._src_reg) != 1: + raise RuntimeError('Only one On region is allowed') + + # Set up source direction. Query parameters if neccessary. + if isinstance(self._src_reg[0], gammalib.GSkyRegionCircle): + self._src_dir = self._src_reg[0].centre() + self._rad = self._src_reg[0].radius() + else: + self._query_src_direction() + + # Make sure that all CTA observations have an Off region by loading the + # Off region region the parameter 'bkgregfile' for all CTA observations + # without Off region + for obs in self.obs(): + if obs.instrument() == 'CTA': + if obs.off_regions().is_empty(): + filename = self['bkgregfile'].filename() + regions = gammalib.GSkyRegions(filename) + obs.off_regions(regions) + + # Return + return + + def _get_parameters_bkgmethod(self): + """ + Get background method parameters + """ + # Get background method bkgmethod = self['bkgmethod'].string() + + # Get background method dependent parameters if bkgmethod == 'REFLECTED': - self['bkgregmin'].integer() + self._get_parameters_bkgmethod_reflected() + elif bkgmethod == 'CUSTOM': + self._get_parameters_bkgmethod_custom() + + # Query parameters that are needed for all background methods self['maxoffset'].real() + # Return + return + + def _get_parameters(self): + """ + Get parameters from parfile and setup observations + """ + # Setup observations (require response and allow event list, don't + # allow counts cube) + self._setup_observations(self.obs(), True, True, False) + + # Set energy bounds + self._ebounds = self._create_ebounds() + + # Initialize empty src regions container + self._src_reg = gammalib.GSkyRegions() + # Exclusion map if self['inexclusion'].filename().is_fits(): self._excl_reg = gammalib.GSkyRegionMap(self['inexclusion'].filename()) @@ -109,6 +185,10 @@ class csphagen(ctools.csobservation): if self.obs().is_empty(): self.obs(self._get_observations(False)) + # Get background method parameters (have to come after setting up of + # observations) + self._get_parameters_bkgmethod() + # Write input parameters into logger self._log_parameters(gammalib.TERSE) @@ -247,6 +327,38 @@ class csphagen(ctools.csobservation): # Return energy boundaries return ebounds + def _set_background_regions(self, obs): + """ + Set background regions for an observation + + Parameters + ---------- + obs : `~gammalib.GCTAObservation()` + CTA observation + + Returns + ------- + regions : `~gammalib.GSkyRegions()` + Background regions + """ + # Initialise empty background regions for this observation + bkg_reg = gammalib.GSkyRegions() + + # If reflected background is requested then created reflected + # background regions + if self['bkgmethod'].string() == 'REFLECTED': + bkg_reg = self._reflected_regions(obs) + + # ... otherwise if custom background is requested then get the + # background regions from the observation. We use a copy here since + # otherwise the background regions go out of scope once the observations + # are replaced by the On/Off observations. + elif self['bkgmethod'].string() == 'CUSTOM': + bkg_reg = obs.off_regions().copy() + + # Return background regions + return bkg_reg + # Public methods def run(self): @@ -286,15 +398,20 @@ class csphagen(ctools.csobservation): # Loop through observations and generate pha, arf, rmf files for obs in self.obs(): - # Initialise background regions for this observation - bkg_reg = gammalib.GSkyRegions() + # Skip non CTA observations + if obs.instrument() != 'CTA': + self._log_string(gammalib.NORMAL, 'Skip %s observation "%s"' % \ + (obs.instrument(), obs.id())) + continue + + # Log current observation + self._log_string(gammalib.NORMAL, 'Process observation "%s"' % \ + (obs.id())) - # If reflected background is requested then created reflected - # background regions - if self['bkgmethod'].string() == 'REFLECTED': - bkg_reg = self._reflected_regions(obs) + # Set background regions for this observation + bkg_reg = self._set_background_regions(obs) - # If there are reflected regions then create On/Off observation + # If there are background regions then create On/Off observation # and append it to the output container if bkg_reg.size() >= self['bkgregmin'].integer(): onoff = gammalib.GCTAOnOffObservation(obs, self._src_dir, diff --git a/cscripts/obsutils.py b/cscripts/obsutils.py index 5138794..9d80c96 100755 --- a/cscripts/obsutils.py +++ b/cscripts/obsutils.py @@ -130,6 +130,7 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0, phagen['rad'] = onrad phagen['stack'] = False phagen['inexclusion'] = 'NONE' + phagen['bkgmethod'] = 'REFLECTED' # Optionally open the log file if log: diff --git a/doc/source/admin/release_history/1.5.rst b/doc/source/admin/release_history/1.5.rst index b11d7ea..d3977cc 100644 --- a/doc/source/admin/release_history/1.5.rst +++ b/doc/source/admin/release_history/1.5.rst @@ -29,6 +29,8 @@ Improvements * [`2295 `_] - Add ``ctools::time_reference`` constant and remove ``ctools::time_reference()`` method +* [`2284 `_] - + Add ``CUSTOM`` background method to :ref:`csphagen` script * [`2273 `_] - Support On/Off analysis in :ref:`cslightcrv` and :ref:`csphasecrv` * [`2267 `_] - diff --git a/doc/source/users/reference_manual/csphagen.rst b/doc/source/users/reference_manual/csphagen.rst index 61176ae..9bcb5fe 100644 --- a/doc/source/users/reference_manual/csphagen.rst +++ b/doc/source/users/reference_manual/csphagen.rst @@ -13,8 +13,9 @@ Synopsis This script can be used to derive from an observation or a set of observations the count spectra in a source region and in background regions, as well as the detector response (effective area, energy redistribution matrix), that can be -used to perform a classical 1D spectral analysis. The output files are saved in -the OGIP format normally used in X-ray astronomy (PHA, ARF, RMF). +used to perform a classical 1D spectral analysis. Regions can be placed +automatically or by hand. The output files are saved in the OGIP format normally +used in X-ray astronomy (PHA, ARF, RMF). `See here `__. The outputs are: @@ -30,50 +31,48 @@ General parameters ------------------ ``inobs [file]`` - Input event list or observation definition XML file. + Input event list or observation definition XML file ``caldb [string]`` - Calibration database. + Calibration database ``irf [string]`` - Instrument response function. + Instrument response function ``(inexclusion = NONE) [file]`` Optional FITS file containing a WCS map in the first hdu that defines sky - regions not to be used for background estimation (where map value != 0). + regions not to be used for background estimation (where map value != 0) ``outobs [string]`` - Output observation definition XML file. + Output observation definition XML file ``(prefix = onoff) [string]`` - Prefix of the file name for output PHA, ARF, RMF, XML, and DS9 region files. + Prefix of the file name for output PHA, ARF, RMF, XML, and DS9 region files ``emin [real]`` - Lower energy limit (TeV) if ``LIN`` or ``LOG`` binning algorithms are used. + Lower energy limit (TeV) if ``LIN`` or ``LOG`` binning algorithms are used ``emax [real]`` - Upper energy limit (TeV) if ``LIN`` or ``LOG`` binning algorithms are used. + Upper energy limit (TeV) if ``LIN`` or ``LOG`` binning algorithms are used ``enumbins [integer]`` Number of energy bins. At least 30 bins per decade are recommended for proper evaluation of the instrument response. -``ebinalg [string]`` - Energy binning algorithm: - ``LOG`` : logarithmically spaced energy bins; - ``LIN`` : linearly spaced energy bins; - ``FILE``: energy bounds retrieved from file. +``ebinalg [string]`` + Algorithm for defining energy bins (``FILE``: energy bounds retrieved from + file, see ``ebinfile`` parameter; ``LIN``: linearly spaced energy bins; + ``LOG``: logarithmically spaced energy bins) ``ebinfile [file]`` Name of the file containing the energy bin definition if ``FILE`` algorithm - is used. + is used -``(srcshape = CIRCLE) [string]`` - Shape of the source region. - ``CIRCLE``: circular region around given position. +``(srcshape = CIRCLE) [string]`` + Shape of the source region (``CIRCLE``: circular region around given position) -``coordsys [string]`` - Coordinate system (``CEL`` - celestial, ``GAL`` - galactic) +``coordsys [string]`` + Coordinate system (``CEL``: celestial; ``GAL``: galactic) ``ra [real]`` Right Ascension of source region centre (deg) @@ -90,10 +89,25 @@ General parameters ``rad [real]`` Radius of source region circle (deg) -``(bkgmethod = REFLECTED) [string]`` - Method for background estimation. - ``REFLECTED:`` background evaluated in regions with the same shape as - source region reflected w.r.t. pointing direction for each observation. +``srcregfile [file]`` + Source region file (ds9 or FITS WCS map) + +``bkgmethod [string]`` + Method for background estimation: + + - ``REFLECTED``: background is evaluated in regions with the same shape as + the source region reflected w.r.t. pointing direction for each observation + + - ``CUSTOM``: background is evaluated in regions specified by user. For an + event list or a single observation in the observation definition XML file + a region file will be queried (see ``bkgregfile`` parameter). For multiple + observations specified in the observation definition XML file the name of + the region file will be extracted from the ``OffRegions`` parameter that + needs to be specified for each observation in the observation definition + XML file. Off region files can be either ds9 region files or FITS WCS maps. + +``bkgregfile [file]`` + Background regions file (ds9 or FITS WCS map) ``(bkgregmin = 2) [integer]`` Minimum number of background regions that are required for an observation. @@ -102,20 +116,20 @@ General parameters ``(maxoffset = 4.0) [real]`` Maximum offset in degrees of source from camera center to accept the - observation. + observation ``stack [boolean]`` Specifies whether multiple observations should be stacked (``yes``) or - whether run-wise PHA, ARF and RMF files should be produced (``no``). + whether run-wise PHA, ARF and RMF files should be produced (``no``) ``(etruemin = 0.01) [real]`` - Minimum true energy (TeV). + Minimum true energy (TeV) ``(etruemax = 0.01) [real]`` - Maximum true energy (TeV). + Maximum true energy (TeV) ``(etruebins = 30) [integer]`` - Number of bins per decade for true energy bins. + Number of bins per decade for true energy bins Standard parameters diff --git a/test/data/crab_bkg_reg.reg b/test/data/crab_bkg_reg.reg new file mode 100644 index 0000000..5c65fc2 --- /dev/null +++ b/test/data/crab_bkg_reg.reg @@ -0,0 +1,7 @@ +# Region file format: DS9 version 4.1 +global color=green dashlist=8 3 width=1 +fk5;circle(82.2412081162728,22.7975361111674,0.2) +fk5;circle(82.5360705305998,23.1293202741537,0.2) +fk5;circle(82.980058098929,23.260640435668,0.2) +fk5;circle(83.4304728956696,23.1492236356079,0.2) +fk5;circle(83.7424610541293,22.8309432073396,0.2) diff --git a/test/data/crab_src_reg.reg b/test/data/crab_src_reg.reg new file mode 100644 index 0000000..148f7a0 --- /dev/null +++ b/test/data/crab_src_reg.reg @@ -0,0 +1,3 @@ +# Region file format: DS9 version 4.1 +global color=green dashlist=8 3 width=1 +fk5;circle(83.633,22.0145,0.2) diff --git a/test/test_csphagen.py b/test/test_csphagen.py index cc47967..df6d099 100644 --- a/test/test_csphagen.py +++ b/test/test_csphagen.py @@ -2,7 +2,7 @@ # ========================================================================== # This scripts performs unit tests for the csphagen script. # -# Copyright (C) 2017 Luigi Tibaldo +# Copyright (C) 2017-2018 Luigi Tibaldo # # 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 @@ -48,10 +48,11 @@ class Test(test): self._exclusion = self._datadir + '/crab_exclusion.fits' self._nreg_with_excl = 5 self._nreg_wo_excl = 8 + self._nreg_bkg_reg = 5 self._nreg_mul = [self._nreg_with_excl, 6] self._nbins = 10 - # number of expected background regions with/wo exclusion, - # and for two different runs + self._regfile_src = self._datadir + '/crab_src_reg.reg' + self._regfile_bkg = self._datadir + '/crab_bkg_reg.reg' # Return return @@ -88,7 +89,8 @@ class Test(test): ' coordsys="CEL" ra=83.633 dec=22.0145' + \ ' rad=0.2 stack="no" inexclusion="' + \ self._exclusion + \ - '" outobs="csphagen_cmd1.xml" prefix="csphagen_cmd1" ' + \ + '" bkgmethod="REFLECTED" '+ \ + 'outobs="csphagen_cmd1.xml" prefix="csphagen_cmd1" ' + \ 'logfile="csphagen_cmd1.log" chatter=2' # Check if execution of wrong command fails @@ -112,7 +114,8 @@ class Test(test): ' coordsys="CEL" ra=83.633 dec=22.0145' + \ ' rad=0.2 stack="no" inexclusion="' + \ self._exclusion + \ - '" outobs="csphagen_cmd2.xml" prefix="csphagen_cmd2" ' + \ + '" bkgmethod="REFLECTED" '+ \ + 'outobs="csphagen_cmd2.xml" prefix="csphagen_cmd2" ' + \ 'logfile="csphagen_cmd2.log" debug=yes chatter=1' # Check if execution failed @@ -144,6 +147,7 @@ class Test(test): phagen['rad'] = 0.2 phagen['stack'] = False phagen['inexclusion'] = self._exclusion + phagen['bkgmethod'] = 'REFLECTED' phagen['etruemin'] = 0.05 phagen['etruemax'] = 150.0 phagen['etruebins'] = 5 @@ -152,7 +156,7 @@ class Test(test): phagen['logfile'] = 'csphagen_py1.log' phagen['chatter'] = 1 - # Run script + # Execute script phagen.execute() # Check output @@ -173,6 +177,7 @@ class Test(test): phagen['dec'] = 22.0145 phagen['rad'] = 0.2 phagen['stack'] = False + phagen['bkgmethod'] = 'REFLECTED' phagen['etruemin'] = 0.05 phagen['etruemax'] = 150.0 phagen['etruebins'] = 5 @@ -181,7 +186,7 @@ class Test(test): phagen['logfile'] = 'csphagen_py2.log' phagen['chatter'] = 2 - # Run script + # Execute script phagen.execute() # Check output @@ -209,6 +214,7 @@ class Test(test): phagen['rad'] = 0.2 phagen['stack'] = False phagen['inexclusion'] = self._exclusion + phagen['bkgmethod'] = 'REFLECTED' phagen['etruemin'] = 0.05 phagen['etruemax'] = 150.0 phagen['etruebins'] = 5 @@ -226,8 +232,7 @@ class Test(test): self._nreg_mul[s]) self._check_outobs('csphagen_py3', 2) - # Setup csphagen for test with multiple input observations and - # stacking + # Setup csphagen for test with multiple input observations and stacking phagen = cscripts.csphagen(obs) phagen['ebinalg'] = 'LOG' phagen['emin'] = 0.1 @@ -239,6 +244,7 @@ class Test(test): phagen['rad'] = 0.2 phagen['stack'] = True phagen['inexclusion'] = self._exclusion + phagen['bkgmethod'] = 'REFLECTED' phagen['etruemin'] = 0.05 phagen['etruemax'] = 150.0 phagen['etruebins'] = 5 @@ -247,7 +253,7 @@ class Test(test): phagen['logfile'] = 'csphagen_py4.log' phagen['chatter'] = 4 - # Run script + # Execute script phagen.execute() # Check output @@ -256,6 +262,66 @@ class Test(test): 0, check_regions=False) self._check_outobs('csphagen_py4', 1) + # Setup csphagen for test with custom On and Off regions provided + phagen = cscripts.csphagen() + phagen['inobs'] = self._myevents1 + phagen['caldb'] = self._caldb + phagen['irf'] = self._irf + phagen['ebinalg'] = 'LOG' + phagen['emin'] = 0.1 + phagen['emax'] = 100.0 + phagen['enumbins'] = self._nbins + phagen['bkgmethod'] = 'CUSTOM' + phagen['srcregfile'] = self._regfile_src + phagen['bkgregfile'] = self._regfile_bkg + phagen['etruemin'] = 0.05 + phagen['etruemax'] = 150.0 + phagen['etruebins'] = 5 + phagen['stack'] = False + phagen['outobs'] = 'csphagen_py5.xml' + phagen['prefix'] = 'csphagen_py5' + phagen['logfile'] = 'csphagen_py5.log' + phagen['chatter'] = 2 + + # Execute script + phagen.execute() + + # Check output + self._check_output('csphagen_py5', self._nbins, self._nreg_bkg_reg) + self._check_outobs('csphagen_py5', 1) + + # Append off regions to observation container + for run in obs: + run.off_regions(gammalib.GSkyRegions(self._regfile_bkg)) + + # Setup csphagen for test with multiple input observations and stacking + phagen = cscripts.csphagen(obs) + phagen['ebinalg'] = 'LOG' + phagen['emin'] = 0.1 + phagen['emax'] = 100.0 + phagen['enumbins'] = self._nbins + phagen['coordsys'] = 'CEL' + phagen['stack'] = True + phagen['inexclusion'] = self._exclusion + phagen['bkgmethod'] = 'CUSTOM' + phagen['srcregfile'] = self._regfile_src + phagen['etruemin'] = 0.05 + phagen['etruemax'] = 150.0 + phagen['etruebins'] = 5 + phagen['outobs'] = 'csphagen_py6.xml' + phagen['prefix'] = 'csphagen_py6' + phagen['logfile'] = 'csphagen_py6.log' + phagen['chatter'] = 4 + + # Execute script + phagen.execute() + + # Check output + for s in range(2): + self._check_output('csphagen_py6_stacked', self._nbins, + 0, check_regions=False) + self._check_outobs('csphagen_py6', 1) + # Return return @@ -309,7 +375,7 @@ class Test(test): self.test_assert(table.contains(col), 'FITS file contains "' + col + '" column') - # Check SPECTRUM table + # Check EBOUNDS table table = fits['EBOUNDS'] self._check_ebounds(table, bins) @@ -433,3 +499,26 @@ class Test(test): # Return return + + def _read_pha_counts(self, filename): + """ + Read and integrate the counts in a pha file. + Pha file structure already tested in _check_pha() + """ + # Open FITS file + fits = gammalib.GFits(filename) + + # Get SPECTRUM table + table = fits['SPECTRUM'] + + # Integrate counts + counts_col = table['COUNTS'] + counts = 0 + for channel in range( counts_col.nrows() ): + counts += counts_col[channel] + + # Close FITS file + fits.close() + + # Return + return counts