Iftach Sadeh / ctools

Commit b17eaf02157e070096055cd99bb7bfe95174f98d

Authored by Jürgen Knödlseder
1 parent ab0848b3

Integrate csphagen script (#2208)

Inline Side-by-side
Showing 9 changed files with 387 additions and 256 deletions
... ... @@ -15,4 +15,5 @@ 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 +Luigi Tibaldo (luigi.tibaldo@irap.omp.eu)
18 19 Alexander Ziegler (alexander.ziegler@fau.de)
... ...
ChangeLog
  View file @b17eaf0
1   -2017-08-08
  1 +2017-10-16
2 2  
3 3 * Version 1.5.0 released
4 4 ========================
5 5  
  6 + Add csphagen script (#2208)
6 7 Use t-type parameters for all "tmin" and "tmax" parameters (#1864)
7 8  
8 9  
... ...
1 1 New Features and Important Changes in ctools 1.5.0
2 2  
3   -8 August 2017
  3 +16 October 2017
4 4  
5 5  
6 6 Introduction
... ... @@ -122,6 +122,9 @@ cscripts
122 122 Use t-type parameters for "tmin" and "tmax" to all specifying times in various
123 123 formats (#1864).
124 124  
  125 +The script csphagen was added to allow the generation of PHA, ARF and RMF datasets
  126 +for classical On/Off IACT analysis (#2208).
  127 +
125 128  
126 129 Examples
127 130 --------
... ...
cscripts/csobsselect.py
  View file @b17eaf0
... ... @@ -320,7 +320,7 @@ class csobsselect(ctools.cscript):
320 320 # Otherwise log filename and save file
321 321 else:
322 322 # Log filename
323   - self._log_value(gammalib.NORMAL, 'Observation definition XML file',
  323 + self._log_value(gammalib.NORMAL, 'Obs. definition XML file',
324 324 outobs.url())
325 325  
326 326 # Save observations
... ...
cscripts/csphagen.par
  View file @b17eaf0
... ... @@ -16,11 +16,12 @@
16 16 #
17 17 # Input and output parameters
18 18 #============================
19   -inobs, f, a, obs.xml,,, "Input event list or observation definition XML file"
20   -caldb, s, a, prod2,,, "Calibration database"
21   -irf, s, a, South_0.5h,,, "Instrument response function"
22   -exclusion, f, h, NONE,,, "FITS with WCS map of exclusion regions"
23   -outroot, s, a, spectrum,,, "Root of the name for output PHA, ARF, RMF, XML, and ds9 reg files"
  19 +inobs, f, a, obs.xml,,, "Input event list or observation definition XML file"
  20 +caldb, s, a, prod2,,, "Calibration database"
  21 +irf, s, a, South_0.5h,,, "Instrument response function"
  22 +inexclusion, f, h, NONE,,, "Sky map of exclusion regions"
  23 +outobs, s, a, onoff_obs.xml,,, "Output observation definition XML file"
  24 +prefix, s, h, onoff,,, "Prefix for output files"
24 25  
25 26 #
26 27 # Script parameters
... ... @@ -30,7 +31,6 @@ emax, r, a, 100.0,,, "Upper energy limit (TeV)"
30 31 enumbins, i, a, 120,,, "Number of energy bins"
31 32 ebinalg, s, a, LOG,LIN|LOG|FILE,, "Binning algorithm"
32 33 ebinfile, f, a, NONE,,, "Name of the file containing the energy bin definition"
33   -bkgmethod, s, h, REFLECTED,REFLECTED,, "Method for background estimation"
34 34 srcshape, s, h, CIRCLE,CIRCLE,, "Shape of the source region"
35 35 coordsys, s, a, CEL,CEL|GAL,, "Coordinate system (CEL - celestial, GAL - galactic)"
36 36 ra, r, a, 83.63,0,360, "Right Ascension of source region centre (deg)"
... ... @@ -38,9 +38,10 @@ dec, r, a, 22.01,-90,90, "Declination of source region centre (deg)"
38 38 glon, r, a, 184.56,0,360, "Galactic longitude of source region centre (deg)"
39 39 glat, r, a, -5.79,-90,90, "Galactic latitude of source region centre (deg)"
40 40 rad, r, a, 0.2,0,180, "Radius of source region circle (deg)"
41   -bkgregmin, i, h, 2,1,1000, "Minimum number of background region for REFLECTED method"
42   -maxoffset, r, h, 4.,,, "Maximum offset of source from camera center to accept observation"
43   -stack, b, a, no,,, "Stack multiple observations into single On/Off PHA (if no produce run-by-run files)"
  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?"
44 45  
45 46 #
46 47 # Standard parameters
... ...
cscripts/csphagen.py
  View file @b17eaf0
... ... @@ -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 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
... ... @@ -30,7 +30,7 @@ import sys
30 30 # =============== #
31 31 class csphagen(ctools.cscript):
32 32 """
33   - Generate pha files for classical IACT spectral analysis
  33 + Generate PHA, ARF and RMF files for classical IACT spectral analysis
34 34 """
35 35  
36 36 # Constructor
... ... @@ -48,6 +48,15 @@ class csphagen(ctools.cscript):
48 48 # Initialise application by calling the appropriate class constructor
49 49 self._init_cscript(argv)
50 50  
  51 + # Initialise other variables
  52 + self._outobs = gammalib.GObservations()
  53 + self._ebounds = gammalib.GEbounds()
  54 + self._src_dir = gammalib.GSkyDir()
  55 + self._src_reg = gammalib.GSkyRegions()
  56 + self._bkg_regs = []
  57 + self._excl_reg = None
  58 + self._has_exclusion = False
  59 +
51 60 # Return
52 61 return
53 62  
... ... @@ -63,49 +72,43 @@ class csphagen(ctools.cscript):
63 72 # Set energy bounds
64 73 self._ebounds = self._create_ebounds()
65 74  
66   - # Initialise background estimation method
67   - self._bkgmethod = self["bkgmethod"].string()
68   - if self._bkgmethod == "REFLECTED":
69   - self._bkgregmin = self["bkgregmin"].integer()
70   - self._maxoffset = self["maxoffset"].real()
71   -
72 75 # Initialise source position/region querying relevant parameters
73 76 self._src_dir = gammalib.GSkyDir()
74 77 self._src_reg = gammalib.GSkyRegions()
75 78 coordsys = self['coordsys'].string()
76   - if coordsys == "CEL":
77   - ra = self['ra'].real()
  79 + if coordsys == 'CEL':
  80 + ra = self['ra'].real()
78 81 dec = self['dec'].real()
79 82 self._src_dir.radec_deg(ra, dec)
80   - elif coordsys == "GAL":
  83 + elif coordsys == 'GAL':
81 84 glon = self['glon'].real()
82 85 glat = self['glat'].real()
83 86 self._src_dir.lb_deg(glon, glat)
84 87 self._srcshape = self['srcshape'].string()
85   - if self._srcshape == "CIRCLE":
  88 + if self._srcshape == 'CIRCLE':
86 89 self._rad = self['rad'].real()
87   - self._src_reg.append(
88   - gammalib.GSkyRegionCircle(self._src_dir, self._rad))
  90 + self._src_reg.append(gammalib.GSkyRegionCircle(self._src_dir, self._rad))
  91 +
  92 + # Query background estimation method and parameters
  93 + bkgmethod = self['bkgmethod'].string()
  94 + if bkgmethod == 'REFLECTED':
  95 + self['bkgregmin'].integer()
  96 + self['maxoffset'].real()
89 97  
90   - # exclusion map
91   - if self["exclusion"].filename().is_fits():
92   - self._excl_reg = gammalib.GSkyRegionMap(
93   - self["exclusion"].filename())
  98 + # Exclusion map
  99 + if self['inexclusion'].filename().is_fits():
  100 + self._excl_reg = gammalib.GSkyRegionMap(self['inexclusion'].filename())
94 101 self._has_exclusion = True
95 102 else:
96 103 self._has_exclusion = False
97 104  
98   - # Stacking
99   - self._stack = self['stack'].boolean()
  105 + # Query remaining parametersStacking
  106 + self['stack'].boolean()
100 107  
101 108 # Query ahead output parameters
102 109 if (self._read_ahead()):
103   - self._outroot = self['outroot'].string()
104   -
105   - # Set some fixed parameters
106   - self._chatter = self["chatter"].integer()
107   - self._clobber = self["clobber"].boolean()
108   - self._debug = self["debug"].boolean()
  110 + self['outobs'].string()
  111 + self['prefix'].string()
109 112  
110 113 # If there are no observations in container then get them from the
111 114 # parameter file
... ... @@ -121,50 +124,66 @@ class csphagen(ctools.cscript):
121 124 def _reflected_regions(self, obs):
122 125 """
123 126 Calculate list of reflected regions for a single observation (pointing)
124   - :param obs: observation
125   - :return: list of reflected regions
  127 +
  128 + Parameters
  129 + ----------
  130 + obs : `~gammalib.GObservation()`
  131 + CTA observation
  132 +
  133 + Returns
  134 + -------
  135 + regions : `~gammalib.GSkyRegions`
  136 + List of reflected regions
126 137 """
127   - outregions = []
  138 + # Initialise list of reflected regions
  139 + regions = gammalib.GSkyRegions()
  140 +
  141 + # Get offset angle of source
128 142 pnt_dir = obs.pointing().dir()
129   - offset = pnt_dir.dist_deg(self._src_dir)
130   - if offset <= self._rad or offset >= self["maxoffset"].real():
131   - if self._logExplicit():
132   - msg = 'Observation {} pointed at {} deg from source\n'.format(
133   - obs.id(), offset)
134   - self._log(msg)
  143 + offset = pnt_dir.dist_deg(self._src_dir)
  144 +
  145 + # If ...
  146 + if offset <= self._rad or offset >= self['maxoffset'].real():
  147 + self._log_string(gammalib.EXPLICIT, 'Observation {} pointed at {} '
  148 + 'deg from source'.format(obs.id(), offset))
  149 +
  150 + # ... otherwise
135 151 else:
136 152 posang = pnt_dir.posang_deg(self._src_dir)
137   - if self._srcshape == "CIRCLE":
  153 + if self._srcshape == 'CIRCLE':
138 154 # angular separation of reflected regions wrt camera center
139 155 # and number
140 156 alpha = 1.05 * 2 * self._rad / offset
141   - # 1.05 ensures background regions do not overlap due to numerical precision issues
142   - N = int(2 * math.pi / alpha)
143   - if N < self._bkgregmin + 3:
144   - if self._logExplicit():
145   - msg = 'Observation {}: insufficient regions for background estimation\n'.format(
146   - obs.id(), offset)
147   - self._log = (msg)
  157 + # 1.05 ensures background regions do not overlap due to
  158 + # numerical precision issues
  159 + N = int(2.0 * math.pi / alpha)
  160 + if N < self['bkgregmin'].integer() + 3:
  161 + self._log_string(gammalib.EXPLICIT, 'Observation {}: '
  162 + 'insufficient regions for background '
  163 + 'estimation'.format(obs.id(), offset))
  164 +
  165 + # Otherwise loop over position angle to create reflected
  166 + # regions
148 167 else:
149   - alpha = 360. / N
150   - # loop to create reflected regions
  168 + alpha = 360.0 / N
151 169 for s in range(2, N - 1):
152   - dphi = s * alpha
  170 + dphi = s * alpha
153 171 ctr_dir = pnt_dir.clone()
154 172 ctr_dir.rotate_deg(posang + dphi, offset)
155 173 region = gammalib.GSkyRegionCircle(ctr_dir, self._rad)
156 174 if self._has_exclusion:
157 175 if self._excl_reg.overlaps(region):
158   - if self._logVerbose():
159   - msg = 'Observation {}: reflected region overlaps with exclusion region\n'.format(
160   - obs.id(), offset)
161   - self._log(msg)
  176 + self._log_string(gammalib.VERBOSE, 'Observation '
  177 + '{}: reflected region overlaps '
  178 + 'with exclusion region'.
  179 + format(obs.id(), offset))
162 180 else:
163   - outregions.append(region)
  181 + regions.append(region)
164 182 else:
165   - outregions.append(region)
  183 + regions.append(region)
166 184  
167   - return outregions
  185 + # Return reflected regions
  186 + return regions
168 187  
169 188 # Public methods
170 189 def run(self):
... ... @@ -181,74 +200,135 @@ class csphagen(ctools.cscript):
181 200 # Write observation into logger
182 201 self._log_observations(gammalib.NORMAL, self._obs, 'Observation')
183 202  
  203 + # Write header
184 204 self._log_header1(gammalib.NORMAL,
185 205 'Generation of source and background spectra')
  206 +
  207 + # Initialise run variables
  208 + self._outobs = gammalib.GObservations()
  209 + self._bkg_regs = []
  210 +
186 211 # Loop through observations and generate pha, arf, rmf files
187   - # Save ON/OFF regions in ds9 reg files
188   - outobs = gammalib.GObservations()
189   - etrue = self._ebounds
190   - ereco = self._ebounds
191 212 for obs in self._obs:
  213 +
  214 + # Initialise background regions for this observation
192 215 bkg_reg = gammalib.GSkyRegions()
193   - if self._bkgmethod == "REFLECTED":
194   - regions = self._reflected_regions(obs)
195   - for region in regions:
196   - bkg_reg.append(region)
197   - if bkg_reg.size() >= 1 or (
198   - self._bkgmethod == "REFLECTED" and bkg_reg.size() >= self._bkgregmin):
199   - onoff = gammalib.GCTAOnOffObservation(obs, self._src_dir, etrue,
200   - ereco, self._src_reg,
  216 +
  217 + # If reflected background is requested then created reflected
  218 + # background regions
  219 + if self['bkgmethod'].string() == 'REFLECTED':
  220 + bkg_reg = self._reflected_regions(obs)
  221 +
  222 + # If there are reflected regions then create On/Off observation
  223 + # and append it to the output container
  224 + if bkg_reg.size() >= self['bkgregmin'].integer():
  225 + onoff = gammalib.GCTAOnOffObservation(obs, self._src_dir,
  226 + self._ebounds,
  227 + self._ebounds,
  228 + self._src_reg,
201 229 bkg_reg)
202 230 onoff.id(obs.id())
203   - outobs.append(onoff)
204   - self._src_reg.save(
205   - self._outroot + '_{}_on.reg'.format(obs.id()))
206   - bkg_reg.save(self._outroot + '_{}_off.reg'.format(obs.id()))
  231 + self._outobs.append(onoff)
  232 + self._bkg_regs.append({'regions': bkg_reg, 'id': obs.id()})
207 233 else:
208   - if self._logNormal():
209   - msg = 'Observation {} not included in spectra generation\n'.format(
210   - obs.id())
211   - self._log(msg)
212   -
213   - # Save PHA, ARF and RMFs
214   - if outobs.size() < 2 or self._stack == False:
215   - for obs in outobs:
216   - obs.on_spec().save(
217   - self._outroot + '_{}_pha_on.fits'.format(obs.id()),
218   - True)
219   - obs.off_spec().save(
220   - self._outroot + '_{}_pha_off.fits'.format(obs.id()),
221   - True)
222   - obs.arf().save(self._outroot + '_{}_arf.fits'.format(obs.id()),
223   - True)
224   - obs.rmf().save(self._outroot + '_{}_rmf.fits'.format(obs.id()),
225   - True)
226   - else:
227   - if self._logVerbose():
228   - msg = 'Stacking {} observations\n'.format(
229   - outobs.size())
230   - self._log(msg)
231   - stackobs = gammalib.GCTAOnOffObservation(outobs)
232   - stackobs.on_spec().save(
233   - self._outroot + '_stacked_pha_on.fits',
234   - True)
235   - stackobs.off_spec().save(
236   - self._outroot + '_stacked_pha_off.fits',
237   - True)
238   - stackobs.arf().save(self._outroot + '_stacked_arf.fits',
239   - True)
240   - stackobs.rmf().save(self._outroot + '_stacked_rmf.fits',
241   - True)
242   - outobs = gammalib.GObservations()
243   - outobs.append(stackobs)
244   -
245   - # Save On/Off observations
246   - outname = self._outroot + '.xml'
247   - outobs.save(outname)
248   - # Log filename
249   - if self._logNormal():
250   - self._log(
251   - 'Output observation definition XML file: {}\n'.format(outname))
  234 + self._log_string(gammalib.NORMAL, 'Observation {} not included '
  235 + 'in spectra generation'.format(obs.id()))
  236 +
  237 + # Stack observations
  238 + if self._outobs.size() > 1 and self['stack'].boolean() == True:
  239 +
  240 + # Write header
  241 + self._log_header1(gammalib.NORMAL,
  242 + 'Stacking {} observations'.format(self._outobs.size()))
  243 +
  244 + # Stack observations
  245 + stacked_obs = gammalib.GCTAOnOffObservation(self._outobs)
  246 +
  247 + # Put stacked observations in output container
  248 + self._outobs = gammalib.GObservations()
  249 + self._outobs.append(stacked_obs)
  250 +
  251 + # Return
  252 + return
  253 +
  254 + def save(self):
  255 + """
  256 + Save data
  257 + """
  258 + # Write header
  259 + self._log_header1(gammalib.TERSE, 'Save data')
  260 +
  261 + # Get XMK output filename, prefix and clobber
  262 + outobs = self['outobs'].string()
  263 + prefix = self['prefix'].string()
  264 + clobber = self['clobber'].boolean()
  265 +
  266 + # Loop over all observation in container
  267 + for obs in self._outobs:
  268 +
  269 + # Set filenames
  270 + if self['stack'].boolean():
  271 + onname = prefix + '_stacked_pha_on.fits'
  272 + offname = prefix + '_stacked_pha_off.fits'
  273 + arfname = prefix + '_stacked_arf.fits'
  274 + rmfname = prefix + '_stacked_rmf.fits'
  275 + elif self._outobs.size() > 1:
  276 + onname = prefix + '_{}_pha_on.fits'.format(obs.id())
  277 + offname = prefix + '_{}_pha_off.fits'.format(obs.id())
  278 + arfname = prefix + '_{}_arf.fits'.format(obs.id())
  279 + rmfname = prefix + '_{}_rmf.fits'.format(obs.id())
  280 + else:
  281 + onname = prefix + '_pha_on.fits'
  282 + offname = prefix + '_pha_off.fits'
  283 + arfname = prefix + '_arf.fits'
  284 + rmfname = prefix + '_rmf.fits'
  285 +
  286 + # Save files
  287 + obs.on_spec().save(onname, clobber)
  288 + obs.off_spec().save(offname, clobber)
  289 + obs.arf().save(arfname, clobber)
  290 + obs.rmf().save(rmfname, clobber)
  291 +
  292 + # Log file names
  293 + self._log_value(gammalib.NORMAL, 'PHA on file', onname)
  294 + self._log_value(gammalib.NORMAL, 'PHA off file', offname)
  295 + self._log_value(gammalib.NORMAL, 'ARF file', arfname)
  296 + self._log_value(gammalib.NORMAL, 'RMF file', arfname)
  297 +
  298 + # Save observation definition XML file
  299 + self._outobs.save(outobs)
  300 +
  301 + # Log file name
  302 + self._log_value(gammalib.NORMAL, 'Obs. definition XML file', outobs)
  303 +
  304 + # Save ds9 On region file
  305 + regname = prefix + '_on.reg'
  306 + self._src_reg.save(regname)
  307 + self._log_value(gammalib.NORMAL, 'On region file', regname)
  308 +
  309 + # Save ds9 Off region files
  310 + for bkg_reg in self._bkg_regs:
  311 +
  312 + # Set filename
  313 + if len(self._bkg_regs) > 1:
  314 + regname = prefix + '_{}_off.reg'.format(bkg_reg['id'])
  315 + else:
  316 + regname = prefix + '_off.reg'
  317 +
  318 + # Save ds9 region file
  319 + bkg_reg['regions'].save(regname)
  320 +
  321 + # Log file name
  322 + self._log_value(gammalib.NORMAL, 'Off region file', regname)
  323 +
  324 + # Return
  325 + return
  326 +
  327 + def obs(self):
  328 + """
  329 + Return observation container
  330 + """
  331 + return self._outobs
252 332  
253 333 def execute(self):
254 334 """
... ... @@ -263,6 +343,9 @@ class csphagen(ctools.cscript):
263 343 # Run the script
264 344 self.run()
265 345  
  346 + # Save result
  347 + self.save()
  348 +
266 349 # Return
267 350 return
268 351  
... ... @@ -271,6 +354,7 @@ class csphagen(ctools.cscript):
271 354 # Main routine entry point #
272 355 # ======================== #
273 356 if __name__ == '__main__':
  357 +
274 358 # Create instance of application
275 359 app = csphagen(sys.argv)
276 360  
... ...
doc/source/admin/release_history/1.5.rst
  View file @b17eaf0
... ... @@ -10,6 +10,7 @@ ctools 1.5 is a major release that adds significant functionality.
10 10  
11 11 In particular, this release provides:
12 12  
  13 +* Support of classical IACT analysis
13 14 * Specifcation of times in various formats (UTC, MJD, JD, MET)
14 15  
15 16  
... ... @@ -28,5 +29,7 @@ Improvements
28 29 New features
29 30 ------------
30 31  
  32 +* [`2208 <https://cta-redmine.irap.omp.eu/issues/2208>`_] -
  33 + Add :ref:`csphagen` script
31 34 * [`1864 <https://cta-redmine.irap.omp.eu/issues/1864>`_] -
32 35 Use t-type parameters for all ``tmin`` and ``tmax`` parameters
... ...
doc/source/users/reference_manual/csphagen.rst
  View file @b17eaf0
1 1 .. _csphagen:
2 2  
3 3 csphagen
4   -=========
  4 +========
5 5  
6   -Creates files necessary to perform a classical region-based spectral analysis on IACT data, i.e., source and background PHA files, ARF and RMF files.
  6 +Creates files necessary to perform a classical region-based spectral analysis
  7 +on IACT data, i.e., source and background PHA files, ARF and RMF files.
7 8  
8 9  
9 10 Synopsis
10 11 --------
11 12  
12   -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). `See here <https://heasarc.gsfc.nasa.gov/docs/heasarc/ofwg/docs/spectra/ogip_92_007/node5.html>`__.
  13 +This script can be used to derive from an observation or a set of observations
  14 +the count spectra in a source region and in background regions, as well as the
  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).
  18 +`See here <https://heasarc.gsfc.nasa.gov/docs/heasarc/ofwg/docs/spectra/ogip_92_007/node5.html>`__.
13 19  
14   -The outputs are: 1) the PHA, ARF, RMF files, either separately for each observation, or stacked for a set of observations, 2) DS9 regions files listing the source and background regions for each observation; 3) a new observation definition XML file.
  20 +The outputs are:
  21 +1) the PHA, ARF, RMF files, either separately for each observation, or stacked
  22 + for a set of observations,
  23 +2) DS9 regions files listing the source and background regions for each
  24 + observation;
  25 +3) a new observation definition XML file.
15 26  
16 27  
17 28 General parameters
18 29 ------------------
19 30  
20 31 ``inobs [file]``
21   - Input event list or observation definition XML file
  32 + Input event list or observation definition XML file.
22 33  
23 34 ``caldb [string]``
24 35 Calibration database.
25   -
  36 +
26 37 ``irf [string]``
27 38 Instrument response function.
28 39  
29   -``(exclusion = NONE) [file]``
30   - 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)
  40 +``(inexclusion = NONE) [file]``
  41 + Optional FITS file containing a WCS map in the first hdu that defines sky
  42 + regions not to be used for background estimation (where map value != 0).
31 43  
32   -``outroot [string]``
33   - Root of the file name for output PHA, ARF, RMF, XML, and ds9 reg files
  44 +``outobs [string]``
  45 + Output observation definition XML file.
  46 +
  47 +``(prefix = onoff) [string]``
  48 + Prefix of the file name for output PHA, ARF, RMF, XML, and DS9 region files.
34 49  
35 50 ``emin [real]``
36   - Lower energy limit (TeV) if LIN or LOG binning algorithms are used
  51 + Lower energy limit (TeV) if LIN or LOG binning algorithms are used.
37 52  
38 53 ``emax [real]``
39   - Upper energy limit (TeV) if LIN or LOG binning algorithms are used
  54 + Upper energy limit (TeV) if LIN or LOG binning algorithms are used.
40 55  
41 56 ``enumbins [integer]``
42   - Number of energy bins. At least 30 bins per decade are recommended for proper evaluation of the instrument response
  57 + Number of energy bins. At least 30 bins per decade are recommended for
  58 + proper evaluation of the instrument response.
43 59  
44 60 ``ebinalg [string]’’
45   - Energy binning algorithm: LOG - logarithmically spaced energy bins; LIN - linearly spaced energy bins; FILE - energy bounds retrieved from file
  61 + Energy binning algorithm: LOG - logarithmically spaced energy bins;
  62 + LIN - linearly spaced energy bins; FILE - energy bounds retrieved from file.
46 63  
47 64 ``ebinfile [file]``
48   - Name of the file containing the energy bin definition if FILE algorithm is used
49   -
50   -``(bkgmethod = REFLECTED) [string]’’
51   - 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
  65 + Name of the file containing the energy bin definition if FILE algorithm
  66 + is used.
52 67  
53 68 ``(srcshape = CIRCLE) [string]``
54   - Shape of the source region. CIRCLE: circular region around given position
  69 + Shape of the source region.
  70 + CIRCLE: circular region around given position.
55 71  
56 72 ``coordsys [string]``
57 73 Coordinate system (CEL - celestial, GAL - galactic)
... ... @@ -71,14 +87,23 @@ General parameters
71 87 ``rad [real]``
72 88 Radius of source region circle (deg)
73 89  
  90 +``(bkgmethod = REFLECTED) [string]’’
  91 + Method for background estimation.
  92 + REFLECTED: background evaluated in regions with the same shape as source
  93 + region reflected w.r.t. pointing direction for each observation.
  94 +
74 95 ``(bkgregmin = 2) [integer]``
75   - Minimum number of background regions for REFLECTED method. If not available, observation is skipped. (Default: 2)
  96 + Minimum number of background regions that are required for an observation.
  97 + If this number of background regions is not available the observation is
  98 + skipped.
76 99  
77   -``(maxoffset = 4.) [real]``
78   - Maximum offset of source from camera center to accept observation
  100 +``(maxoffset = 4.0) [real]``
  101 + Maximum offset in degrees of source from camera center to accept the
  102 + observation.
79 103  
80 104 ``stack [boolean]``
81   - Yes: Stack multiple observations into single PHA. No: produce run-by-run PHA files
  105 + Specifies whether multiple observations should be stacked (``yes``) or
  106 + whether run-wise PHA, ARF and RMF files should be produced (``no``).
82 107  
83 108  
84 109 Standard parameters
... ... @@ -100,13 +125,13 @@ Standard parameters
100 125 Specifies whether an existing output runlist should be overwritten.
101 126  
102 127 ``(debug = no) [boolean]``
103   - Enables debug mode. In debug mode the executable will dump any log file output to the console.
  128 + Enables debug mode. In debug mode the executable will dump any log file
  129 + output to the console.
104 130  
105 131 ``(mode = ql) [string]``
106 132 Mode of automatic parameters (default is "ql", i.e. "query and learn").
107 133  
108   -``(logfile = csfindobs.log) [filename]``
  134 +``(logfile = csphagen.log) [filename]``
109 135 Log filename.
110 136  
111 137  
112   -
... ...
test/test_csphagen.py
  View file @b17eaf0
... ... @@ -41,12 +41,13 @@ class Test(test):
41 41 """
42 42 # Call base class constructor
43 43 test.__init__(self)
44   - self._myevents1 = self._datadir + '/crab_offaxis1.fits'
45   - self._myevents2 = self._datadir + '/crab_offaxis2.fits'
46   - self._exclusion = self._datadir + '/crab_exclusion.fits'
  44 + self._myevents1 = self._datadir + '/crab_offaxis1.fits'
  45 + self._myevents2 = self._datadir + '/crab_offaxis2.fits'
  46 + self._exclusion = self._datadir + '/crab_exclusion.fits'
47 47 self._nreg_with_excl = 5
48   - self._nreg_wo_excl = 8
49   - self._nreg_mul = [self._nreg_with_excl, 6]
  48 + self._nreg_wo_excl = 8
  49 + self._nreg_mul = [self._nreg_with_excl, 6]
  50 + self._nbins = 20
50 51 # number of expected background regions with/wo exclusion,
51 52 # and for two different runs
52 53  
... ... @@ -75,14 +76,17 @@ class Test(test):
75 76 """
76 77 # Set script name
77 78 csphagen = self._script('csphagen')
78   - nbins = 120
79 79  
80   - cmd = csphagen + ' inobs="' + self._myevents1 + '" caldb="' + \
81   - self._caldb + '" irf="' + self._irf + \
82   - '" ebinalg="LOG" emin=0.1 emax=100. enumbins=' + \
83   - str(nbins) + ' coordsys="CEL" ra=83.633 dec=22.0145' + \
84   - ' rad=0.2 stack="no" exclusion="' + self._exclusion + \
85   - '" outroot="phagen_cmd1" logfile="csphagen_cmd1.log" chatter=1'
  80 + # Setup csphagen command
  81 + cmd = csphagen + ' inobs="' + self._myevents1 + \
  82 + '" caldb="' + self._caldb + '" irf="' + self._irf + \
  83 + '" ebinalg="LOG" emin=0.1 emax=100. enumbins=' + \
  84 + str(self._nbins) + \
  85 + ' coordsys="CEL" ra=83.633 dec=22.0145' + \
  86 + ' rad=0.2 stack="no" inexclusion="' + \
  87 + self._exclusion + \
  88 + '" outobs="csphagen_cmd1.xml" prefix="csphagen_cmd1" ' + \
  89 + 'logfile="csphagen_cmd1.log" chatter=2'
86 90  
87 91 # Check if execution of wrong command fails
88 92 self.test_assert(self._execute('command_that_does_not_exist') != 0,
... ... @@ -93,15 +97,19 @@ class Test(test):
93 97 'Check successful execution from command line')
94 98  
95 99 # Check output files
96   - self._check_output('phagen_cmd1_', nbins, self._nreg_with_excl)
97   - self._check_outobs('phagen_cmd1', 1)
98   -
99   - cmd = csphagen + ' inobs="events_that_do_not_exist.fits" caldb="' + \
100   - self._caldb + '" irf="' + self._irf + \
101   - '" ebinalg="LOG" emin=0.1 emax=100. enumbins=' + \
102   - str(nbins) + ' coordsys="CEL" ra=83.633 dec=22.0145' + \
103   - ' rad=0.2 stack="no" exclusion="' + self._exclusion + \
104   - '" outroot="phagen_cmd1" logfile="csphagen_cmd1.log" chatter=1'
  100 + self._check_output('csphagen_cmd1', self._nbins, self._nreg_with_excl)
  101 + self._check_outobs('csphagen_cmd1', 1)
  102 +
  103 + # Setup csphagen command
  104 + cmd = csphagen + ' inobs="events_that_do_not_exist.fits" ' + \
  105 + 'caldb="' + self._caldb + '" irf="' + self._irf + \
  106 + '" ebinalg="LOG" emin=0.1 emax=100. enumbins=' + \
  107 + str(self._nbins) + \
  108 + ' coordsys="CEL" ra=83.633 dec=22.0145' + \
  109 + ' rad=0.2 stack="no" inexclusion="' + \
  110 + self._exclusion + \
  111 + '" outobs="csphagen_cmd2.xml" prefix="csphagen_cmd2" ' + \
  112 + 'logfile="csphagen_cmd2.log" chatter=1'
105 113  
106 114 # Check if execution failed
107 115 self.test_assert(self._execute(cmd) != 0,
... ... @@ -117,58 +125,58 @@ class Test(test):
117 125 """
118 126 Test csphagen from Python
119 127 """
120   - nbins = 120
121   -
122 128 # Same test as from command line
123 129 phagen = cscripts.csphagen()
124   - phagen['inobs'] = self._myevents1
125   - phagen['caldb'] = self._caldb
126   - phagen['irf'] = self._irf
127   - phagen['ebinalg'] = 'LOG'
128   - phagen['emin'] = 0.1
129   - phagen['emax'] = 100.
130   - phagen['enumbins'] = nbins
131   - phagen['coordsys'] = 'CEL'
132   - phagen['ra'] = 83.633
133   - phagen['dec'] = 22.0145
134   - phagen['rad'] = 0.2
135   - phagen['stack'] = False
136   - phagen['exclusion'] = self._exclusion
137   - phagen['outroot'] = 'phagen_py1'
138   - phagen['logfile'] = 'csphagen_py1.log'
139   - phagen['chatter'] = 1
  130 + phagen['inobs'] = self._myevents1
  131 + phagen['caldb'] = self._caldb
  132 + phagen['irf'] = self._irf
  133 + phagen['ebinalg'] = 'LOG'
  134 + phagen['emin'] = 0.1
  135 + phagen['emax'] = 100.0
  136 + phagen['enumbins'] = self._nbins
  137 + phagen['coordsys'] = 'CEL'
  138 + phagen['ra'] = 83.633
  139 + phagen['dec'] = 22.0145
  140 + phagen['rad'] = 0.2
  141 + phagen['stack'] = False
  142 + phagen['inexclusion'] = self._exclusion
  143 + phagen['outobs'] = 'csphagen_py1.xml'
  144 + phagen['prefix'] = 'csphagen_py1'
  145 + phagen['logfile'] = 'csphagen_py1.log'
  146 + phagen['chatter'] = 1
140 147  
141 148 # Run script
142 149 phagen.execute()
143 150  
144 151 # Check output
145   - self._check_output('phagen_py1_', nbins, self._nreg_with_excl)
146   - self._check_outobs('phagen_py1', 1)
  152 + self._check_output('csphagen_py1', self._nbins, self._nreg_with_excl)
  153 + self._check_outobs('csphagen_py1', 1)
147 154  
148 155 # Now test without exclusion region
149 156 phagen = cscripts.csphagen()
150   - phagen['inobs'] = self._myevents1
151   - phagen['caldb'] = self._caldb
152   - phagen['irf'] = self._irf
153   - phagen['ebinalg'] = 'LOG'
154   - phagen['emin'] = 0.1
155   - phagen['emax'] = 100.
156   - phagen['enumbins'] = nbins
  157 + phagen['inobs'] = self._myevents1
  158 + phagen['caldb'] = self._caldb
  159 + phagen['irf'] = self._irf
  160 + phagen['ebinalg'] = 'LOG'
  161 + phagen['emin'] = 0.1
  162 + phagen['emax'] = 100.0
  163 + phagen['enumbins'] = self._nbins
157 164 phagen['coordsys'] = 'CEL'
158   - phagen['ra'] = 83.633
159   - phagen['dec'] = 22.0145
160   - phagen['rad'] = 0.2
161   - phagen['stack'] = False
162   - phagen['outroot'] = 'phagen_py2'
163   - phagen['logfile'] = 'csphagen_py2.log'
164   - phagen['chatter'] = 1
  165 + phagen['ra'] = 83.633
  166 + phagen['dec'] = 22.0145
  167 + phagen['rad'] = 0.2
  168 + phagen['stack'] = False
  169 + phagen['outobs'] = 'csphagen_py2.xml'
  170 + phagen['prefix'] = 'csphagen_py2'
  171 + phagen['logfile'] = 'csphagen_py2.log'
  172 + phagen['chatter'] = 2
165 173  
166 174 # Run script
167 175 phagen.execute()
168 176  
169 177 # Check output
170   - self._check_output('phagen_py2_', nbins, self._nreg_wo_excl)
171   - self._check_outobs('phagen_py2', 1)
  178 + self._check_output('csphagen_py2', self._nbins, self._nreg_wo_excl)
  179 + self._check_outobs('csphagen_py2', 1)
172 180  
173 181 # Test with multiple input observations, no stacking
174 182 # Create observation container
... ... @@ -181,55 +189,57 @@ class Test(test):
181 189  
182 190 # Setup csphagen
183 191 phagen = cscripts.csphagen(obs)
184   - phagen['ebinalg'] = 'LOG'
185   - phagen['emin'] = 0.1
186   - phagen['emax'] = 100.
187   - phagen['enumbins'] = nbins
188   - phagen['coordsys'] = 'CEL'
189   - phagen['ra'] = 83.633
190   - phagen['dec'] = 22.0145
191   - phagen['rad'] = 0.2
192   - phagen['stack'] = False
193   - phagen['exclusion'] = self._exclusion
194   - phagen['outroot'] = 'phagen_py3'
195   - phagen['logfile'] = 'csphagen_py3.log'
196   - phagen['chatter'] = 1
  192 + phagen['ebinalg'] = 'LOG'
  193 + phagen['emin'] = 0.1
  194 + phagen['emax'] = 100.0
  195 + phagen['enumbins'] = self._nbins
  196 + phagen['coordsys'] = 'CEL'
  197 + phagen['ra'] = 83.633
  198 + phagen['dec'] = 22.0145
  199 + phagen['rad'] = 0.2
  200 + phagen['stack'] = False
  201 + phagen['inexclusion'] = self._exclusion
  202 + phagen['outobs'] = 'csphagen_py3.xml'
  203 + phagen['prefix'] = 'csphagen_py3'
  204 + phagen['logfile'] = 'csphagen_py3.log'
  205 + phagen['chatter'] = 3
197 206  
198 207 # Run script
199 208 phagen.execute()
200 209  
201 210 # Check output
202 211 for s in range(2):
203   - self._check_output('phagen_py3_' + str(s + 1), nbins,
  212 + self._check_output('csphagen_py3_' + str(s + 1), self._nbins,
204 213 self._nreg_mul[s])
205   - self._check_outobs('phagen_py3', 2)
  214 + self._check_outobs('csphagen_py3', 2)
206 215  
207 216 # Test with multiple input observations and stacking
208 217  
209 218 # Setup csphagen
210 219 phagen = cscripts.csphagen(obs)
211   - phagen['ebinalg'] = 'LOG'
212   - phagen['emin'] = 0.1
213   - phagen['emax'] = 100.
214   - phagen['enumbins'] = nbins
215   - phagen['coordsys'] = 'CEL'
216   - phagen['ra'] = 83.633
217   - phagen['dec'] = 22.0145
218   - phagen['rad'] = 0.2
219   - phagen['stack'] = True
220   - phagen['exclusion'] = self._exclusion
221   - phagen['outroot'] = 'phagen_py4'
222   - phagen['logfile'] = 'csphagen_py4.log'
223   - phagen['chatter'] = 1
  220 + phagen['ebinalg'] = 'LOG'
  221 + phagen['emin'] = 0.1
  222 + phagen['emax'] = 100.0
  223 + phagen['enumbins'] = self._nbins
  224 + phagen['coordsys'] = 'CEL'
  225 + phagen['ra'] = 83.633
  226 + phagen['dec'] = 22.0145
  227 + phagen['rad'] = 0.2
  228 + phagen['stack'] = True
  229 + phagen['inexclusion'] = self._exclusion
  230 + phagen['outobs'] = 'csphagen_py4.xml'
  231 + phagen['prefix'] = 'csphagen_py4'
  232 + phagen['logfile'] = 'csphagen_py4.log'
  233 + phagen['chatter'] = 4
224 234  
225 235 # Run script
226 236 phagen.execute()
227 237  
228 238 # Check output
229 239 for s in range(2):
230   - self._check_output('phagen_py4_stacked', nbins,
  240 + self._check_output('csphagen_py4_stacked', self._nbins,
231 241 0, check_regions=False)
232   - self._check_outobs('phagen_py4', 1)
  242 + self._check_outobs('csphagen_py4', 1)
233 243  
234 244 return
235 245  
... ... @@ -274,8 +284,7 @@ class Test(test):
274 284  
275 285 # Check FITS table structure
276 286 self.test_value(table.ncols(), len(cols),
277   - 'Check for %d columns in PHA table' % len(
278   - cols))
  287 + 'Check for %d columns in PHA table' % len(cols))
279 288 self.test_value(table.nrows(), bins,
280 289 'Check for %d rows in PHA table' % bins)
281 290 for col in cols:
... ... @@ -313,8 +322,7 @@ class Test(test):
313 322  
314 323 # Check FITS table structure
315 324 self.test_value(table.ncols(), len(cols),
316   - 'Check for %d columns in ARF table' % len(
317   - cols))
  325 + 'Check for %d columns in ARF table' % len(cols))
318 326 self.test_value(table.nrows(), bins,
319 327 'Check for %d rows in ARF table' % bins)
320 328 for col in cols:
... ... @@ -350,8 +358,7 @@ class Test(test):
350 358  
351 359 # Check FITS table structure
352 360 self.test_value(table.ncols(), len(cols),
353   - 'Check for %d columns in RMF table' % len(
354   - cols))
  361 + 'Check for %d columns in RMF table' % len(cols))
355 362 self.test_value(table.nrows(), bins,
356 363 'Check for %d rows in RMF table' % bins)
357 364 for col in cols:
... ... @@ -372,33 +379,39 @@ class Test(test):
372 379 """
373 380 Check the output from a csphagen run
374 381 """
375   -
376 382 # OGIP files
377 383 self._check_pha(filenameroot + '_pha_on.fits', bins)
378 384 self._check_pha(filenameroot + '_pha_off.fits', bins)
379 385 self._check_arf(filenameroot + '_arf.fits', bins)
380 386 self._check_rmf(filenameroot + '_rmf.fits', bins)
381 387  
382   - # Observations
  388 + # Optionally check for regions
  389 + if check_regions:
383 390  
  391 + # Check On region
  392 + onregion = (filenameroot+'_on.reg').replace('_1_', '_').replace('_2_', '_')
  393 + reg = gammalib.GSkyRegions(onregion)
  394 + self.test_value(reg.size(), 1, 'Check for 1 region in source region file')
384 395  
385   - # Regions
386   - if check_regions:
387   - reg = gammalib.GSkyRegions(filenameroot + "_on.reg")
388   - self.test_value(reg.size(), 1, 'Check for ' + str(
389   - 1) + ' region in source region file')
390   - reg = gammalib.GSkyRegions(filenameroot + "_off.reg")
391   - self.test_value(reg.size(), nreg, 'Check for ' + str(
392   - nreg) + ' region in background region file')
  396 + # Check for Off region
  397 + offregion = filenameroot + '_off.reg'
  398 + reg = gammalib.GSkyRegions(offregion)
  399 + self.test_value(reg.size(), nreg, 'Check for ' + str(nreg) +
  400 + ' region in background region file')
393 401  
  402 + # Return
394 403 return
395 404  
396 405 def _check_outobs(self, filenameroot, nout):
397 406 """
398 407 Check the output XML file containing ON/OFF observations
399 408 """
  409 + # Load observation container
400 410 obs = gammalib.GObservations(filenameroot + '.xml')
401   - self.test_value(obs.size(), nout, 'Check for ' + str(
402   - nout) + ' observations in XML file')
403 411  
  412 + # Check container size
  413 + self.test_value(obs.size(), nout, 'Check for ' + str(nout) +
  414 + ' observations in XML file')
  415 +
  416 + # Return
404 417 return
... ...