ctools  2.0.0
 All Classes Namespaces Files Functions Variables Macros Pages
csphagen.py
Go to the documentation of this file.
1 #! /usr/bin/env python
2 # ==========================================================================
3 # Computes the PHA spectra for source/background and ARF/RMF files using the
4 # reflected region method
5 #
6 # Copyright (C) 2017-2022 Luigi Tibaldo
7 #
8 # This program is free software: you can redistribute it and/or modify
9 # it under the terms of the GNU General Public License as published by
10 # the Free Software Foundation, either version 3 of the License, or
11 # (at your option) any later version.
12 #
13 # This program is distributed in the hope that it will be useful,
14 # but WITHOUT ANY WARRANTY; without even the implied warranty of
15 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 # GNU General Public License for more details.
17 #
18 # You should have received a copy of the GNU General Public License
19 # along with this program. If not, see <http://www.gnu.org/licenses/>.
20 #
21 # ==========================================================================
22 import gammalib
23 import ctools
24 import math
25 import sys
26 from cscripts import mputils
27 
28 
29 # =============== #
30 # csfindobs class #
31 # =============== #
32 class csphagen(ctools.csobservation):
33  """
34  Generate PHA, ARF and RMF files for classical IACT spectral analysis
35  """
36 
37  # Constructor
38  def __init__(self, *argv):
39  """
40  Constructor
41  """
42  # Initialise application by calling the appropriate class constructor
43  self._init_csobservation(self.__class__.__name__, ctools.__version__, argv)
44 
45  # Initialise other variables
46  self._obs_off = gammalib.GObservations()
47  self._ebounds = gammalib.GEbounds()
48  self._etruebounds = gammalib.GEbounds()
49  self._src_dir = gammalib.GSkyDir()
50  self._src_reg = gammalib.GSkyRegions()
51  self._models = gammalib.GModels()
52  self._srcname = ''
53  self._bkg_regs = []
54  self._excl_reg = None
55  self._has_exclusion = False
56  self._srcshape = ''
57  self._rad = 0.0
58  self._reg_width = 0.0
59  self._reg_height = 0.0
60  self._reg_posang = 0.0
61  self._nthreads = 0
62 
63  # Return
64  return
65 
66  # State methods por pickling
67  def __getstate__(self):
68  """
69  Extend ctools.csobservation getstate method to include some members
70 
71  Returns
72  -------
73  state : dict
74  Pickled instance
75  """
76  # Set pickled dictionary
77  state = {'base' : ctools.csobservation.__getstate__(self),
78  'obs_off' : self._obs_off,
79  'ebounds' : self._ebounds,
80  'etruebounds' : self._etruebounds,
81  'src_dir' : self._src_dir,
82  'src_reg' : self._src_reg,
83  'models' : self._models,
84  'srcname' : self._srcname,
85  'bkg_regs' : self._bkg_regs,
86  'excl_reg' : self._excl_reg,
87  'has_exclusion' : self._has_exclusion,
88  'srcshape' : self._srcshape,
89  'rad' : self._rad,
90  'reg_width' : self._reg_width,
91  'reg_height' : self._reg_height,
92  'reg_posang' : self._reg_posang,
93  'nthreads' : self._nthreads}
94 
95  # Return pickled dictionary
96  return state
97 
98  def __setstate__(self, state):
99  """
100  Extend ctools.csobservation setstate method to include some members
101 
102  Parameters
103  ----------
104  state : dict
105  Pickled instance
106  """
107  ctools.csobservation.__setstate__(self, state['base'])
108  self._obs_off = state['obs_off']
109  self._ebounds = state['ebounds']
110  self._etruebounds = state['etruebounds']
111  self._src_dir = state['src_dir']
112  self._src_reg = state['src_reg']
113  self._models = state['models']
114  self._srcname = state['srcname']
115  self._bkg_regs = state['bkg_regs']
116  self._excl_reg = state['excl_reg']
117  self._has_exclusion = state['has_exclusion']
118  self._srcshape = state['srcshape']
119  self._rad = state['rad']
120  self._reg_width = state['reg_width']
121  self._reg_height = state['reg_height']
122  self._reg_posang = state['reg_posang']
123  self._nthreads = state['nthreads']
124 
125  # Return
126  return
127 
128  # Private methods
130  """
131  Set up the source direction parameter
132  """
133  # Initialise source direction
134  self._src_dir = gammalib.GSkyDir()
135 
136  # Get coordinate systel
137  coordsys = self['coordsys'].string()
138 
139  # If coordinate system is celestial then query "ra" and "dec"
140  if coordsys == 'CEL':
141  ra = self['ra'].real()
142  dec = self['dec'].real()
143  self._src_dir.radec_deg(ra, dec)
144 
145  # ... otherwise, if coordinate system is galactic then query "glon"
146  # and "glat"
147  elif coordsys == 'GAL':
148  glon = self['glon'].real()
149  glat = self['glat'].real()
150  self._src_dir.lb_deg(glon, glat)
151 
152  # Return
153  return
154 
155  def _compute_posang(self, pnt_dir, a, b):
156  """
157  Compute the difference in position angle wrt the pointing in degrees
158 
159  Parameters
160  ----------
161  pnt_dir : `~gammalib.GSkyDir`
162  Pointing direction
163  a : `~gammalib.GSkyDir`
164  First sky direction
165  a : `~gammalib.GSkyDir`
166  Second sky direction
167 
168  Returns
169  -------
170  posang : float
171  Position angle (degrees)
172  """
173  # Compute position angles
174  posang_a = pnt_dir.posang_deg(a) % 360
175  posang_b = pnt_dir.posang_deg(b) % 360
176 
177  # Compute difference
178  posang = abs(posang_a - posang_b)
179 
180  # Return position angle
181  return posang
182 
183  def _get_regions(self, filename):
184  """
185  Get regions from DS9 file or FITS file
186 
187  Parameters
188  ----------
189  filename : `~gammalib.GFilename`
190  Filename
191 
192  Returns
193  -------
194  regs : `~gammalib.GSkyRegions`
195  Region container
196  """
197  # If filename is a FITS file then load region map and append to
198  # list of regions
199  if filename.is_fits():
200  map = gammalib.GSkyRegionMap(filename)
201  regs = gammalib.GSkyRegions()
202  regs.append(map)
203 
204  # ... otherwise load DS9 file
205  else:
206  regs = gammalib.GSkyRegions(filename)
207 
208  # Return region container
209  return regs
210 
212  """
213  Get parameters to define source/On region
214  """
215 
216  # Get source shape
217  self._srcshape = self['srcshape'].string()
218 
219  # Query source direction
220  self._query_src_direction()
221 
222  # If source shape is a circle the append GSkyRegionCircle
223  if self._srcshape == 'CIRCLE':
224 
225  # Set circular source region
226  self._rad = self['rad'].real()
227  self._src_reg.append(gammalib.GSkyRegionCircle(self._src_dir, self._rad))
228 
229  # ... otherwise if source shape is a rectangle then append
230  # GSkyRegionRectangle
231  elif self._srcshape == 'RECT':
232 
233  # Set rectangular source region
234  self._reg_width = self['width'].real()
235  self._reg_height = self['height'].real()
236  self._reg_posang = self['posang'].real()
237  self._src_reg.append(gammalib.GSkyRegionRectangle(self._src_dir,
238  self._reg_width,
239  self._reg_height,
240  self._reg_posang))
241 
242  # Return
243  return
244 
246  """
247  Get parameters for REFLECTED background method
248  """
249 
250  # Query parameters for source/On region definition
252 
253  # Query minimum number of background regions and
254  # number of background regions to skip next to On region
255  self['bkgregmin'].integer()
256  self['bkgregskip'].integer()
257 
258  # Return
259  return
260 
262  """
263  Get parameters for CUSTOM background method
264 
265  Raises
266  ------
267  RuntimeError
268  Only one On region is allowed
269  """
270  # Set up source region
271  filename = self['srcregfile'].filename()
272  self._src_reg = self._get_regions(filename)
273 
274  # Raise an exception if there is more than one source region
275  if len(self._src_reg) != 1:
276  raise RuntimeError('Only one On region is allowed')
277 
278  # Set up source direction. Query parameters if neccessary.
279  if self._models.is_empty():
280  if isinstance(self._src_reg[0], gammalib.GSkyRegionCircle):
281  self._src_dir = self._src_reg[0].centre()
282  self._rad = self._src_reg[0].radius()
283  else:
284  self._query_src_direction()
285 
286  # Make sure that all CTA observations have an Off region by loading the
287  # Off region region the parameter 'bkgregfile' for all CTA observations
288  # without Off region
289  for obs in self.obs():
290  if obs.classname() == 'GCTAObservation':
291  if obs.off_regions().is_empty():
292  filename = self['bkgregfile'].filename()
293  regions = self._get_regions(filename)
294  obs.off_regions(regions)
295 
296  # Return
297  return
298 
300  """
301  Get parameters for OFF background method
302 
303  Raises
304  ------
305  RuntimeError
306  On and Off observations must have same size
307  RuntimeError
308  Off observations must be event lists
309  """
310 
311  # Set up Off observations. If there are no Off observations in the
312  # container then load them via user parameters
313  if self.obs_off().is_empty():
314 
315  # Get Off observation file name
316  filename = self['inobsoff'].filename()
317 
318  # If Off observation is a FITS file then load observation and
319  # append it to the Off observation container
320  if gammalib.GFilename(filename).is_fits():
321  self._obs_off.append(gammalib.GCTAObservation(filename))
322 
323  # ... otherwise load XML file into Off observation container
324  else:
325  self._obs_off.load(filename)
326 
327  # Check that size of On and Off observations are the same, otherwise
328  # raise error
329  if self.obs().size() != self.obs_off().size():
330  raise RuntimeError('On and Off observations must have the same size')
331 
332  # Loop through observations
333  for obs in self.obs_off():
334 
335  # Check that observation is event list, otherwise throw error
336  if obs.eventtype() != "EventList":
337  raise RuntimeError('Off observations must be event lists')
338 
339  # Check that they have response, otherwise assign based on user parameter
340  if obs.has_response() == False:
341 
342  # Get database and IRF
343  database = self["caldb"].string()
344  irf = self["irf"].string()
345 
346  # Create an XML element for response
347  parameter = "parameter name=\"Calibration\"" +\
348  " database=\"" + database + "\"" +\
349  " response=\"" + irf + "\""
350  xml = gammalib.GXmlElement()
351  xml.append(parameter)
352 
353  # Create CTA response
354  response = gammalib.GCTAResponseIrf(xml)
355 
356  # Attach response to observation
357  obs.response(response)
358 
359  # Add models from Off observations to model container
360  for model in self.obs_off().models():
361  self._models.append(model)
362 
363  # Query parameters for source/On region definition
365 
366  # Return
367  return
368 
370  """
371  Get background method parameters
372  """
373  # Get background method
374  bkgmethod = self['bkgmethod'].string()
375 
376  # Get background method dependent parameters
377  if bkgmethod == 'REFLECTED':
379  elif bkgmethod == 'CUSTOM':
381  elif bkgmethod == 'OFF':
383 
384  # Query parameters that are needed for all background methods
385  self['maxoffset'].real()
386  self['use_model_bkg'].boolean()
387 
388  # Return
389  return
390 
391  def _get_parameters(self):
392  """
393  Get parameters from parfile and setup observations
394  """
395  # Clear source models
396  self._models.clear()
397 
398  # Setup observations (require response and allow event list, don't
399  # allow counts cube)
400  self._setup_observations(self.obs(), True, True, False)
401 
402  # Get source model and source name. First try to extract models from
403  # observation container. If this does not work then try creating
404  # model from the inmodel parameter
405  if self.obs().models().size() > 0:
406  self._models = self.obs().models().clone()
407  self._srcname = self['srcname'].string()
408  elif self['inmodel'].is_valid():
409  inmodel = self['inmodel'].filename()
410  self._models = gammalib.GModels(inmodel)
411  self._srcname = self['srcname'].string()
412 
413  # Set energy bounds
414  self._ebounds = self._create_ebounds()
415 
416  # Initialize empty src regions container
417  self._src_reg = gammalib.GSkyRegions()
418 
419  # Exclusion map
420  if (self._excl_reg is not None) and (self._excl_reg.map().npix() > 0):
421  # Exclusion map set and is not empty
422  self._has_exclusion = True
423  elif self['inexclusion'].is_valid():
424  inexclusion = self['inexclusion'].filename()
425  # If the user has not specified the extension to use
426  # and there is an extension called 'EXCLUSION' ...
427  if not inexclusion.has_extname()\
428  and not inexclusion.has_extno()\
429  and gammalib.GFits(inexclusion).contains('EXCLUSION'):
430  # ... choose it for the exclusion map
431  extname = inexclusion.url() + '[EXCLUSION]'
432  inexclusion = gammalib.GFilename(extname)
433  # Otherwise will pick the default (primary) HDU
434  self._excl_reg = gammalib.GSkyRegionMap(inexclusion)
435  self._has_exclusion = True
436  else:
437  self._has_exclusion = False
438 
439  # Get background method parameters (have to come after setting up of
440  # observations and models)
442 
443  # If there are multiple observations query whether to stack them
444  if self.obs().size() > 1:
445  self['stack'].boolean()
446 
447  # Query ahead output parameters
448  if (self._read_ahead()):
449  self['outobs'].filename()
450  self['outmodel'].filename()
451  self['prefix'].string()
452 
453  # Write input parameters into logger
454  self._log_parameters(gammalib.TERSE)
455 
456  # Set number of processes for multiprocessing
457  self._nthreads = mputils.nthreads(self)
458 
459  # If we have no model then create now a dummy model
460  if self._models.is_empty():
461  spatial = gammalib.GModelSpatialPointSource(self._src_dir)
462  spectral = gammalib.GModelSpectralPlaw(1.0e-18, -2.0,
463  gammalib.GEnergy(1.0, 'TeV'))
464  model = gammalib.GModelSky(spatial, spectral)
465  model.name('Dummy')
466  self._models.append(model)
467  self._srcname = 'Dummy'
468  self['use_model_bkg'].boolean(False)
469 
470  # Return
471  return
472 
473  def _compute_region_separation(self, pnt_dir):
474  """
475  Compute the separation angle for reflected off regions in radians
476 
477  Returns
478  -------
479  angle : float
480  Separation angle of two off regions (radians)
481  """
482  # Initialise the result
483  separation = -1.0
484 
485  # Compute offset of reflected regions to pointing position
486  offset = pnt_dir.dist_deg(self._src_dir)
487 
488  # If shape is a circle then compute apparent diameter of the circle
489  # as separation
490  if self._srcshape == 'CIRCLE':
491  separation = 2.0 * self._rad / offset
492 
493  # ... otherwise if shape is a rectangle then compute the opening angle
494  # towards combinations of rectangle corners. This method overestimates
495  # the real need of space between the ectangles, so the method may be
496  # optimised to gain more off regions! Anyway, it is assured that the
497  # off regions will never overlap.
498  elif self._srcshape == 'RECT':
499 
500  # Get the sky directions of the corners of the rectangle
501  cs = [self._src_reg[0].corner(icorner) for icorner in range(4)]
502 
503  # Compute the 6 opening angles
504  combinations = [[0,1], [0,2], [0,3], [1,2], [1,3], [2,3]]
505  angles = [self._compute_posang(pnt_dir, cs[i], cs[j]) \
506  for i,j in combinations]
507 
508  # The desired separation is the maximum opening angle
509  separation = max(angles) * gammalib.deg2rad
510 
511  # Return
512  return separation
513 
514  def _reflected_regions(self, obs):
515  """
516  Calculate list of reflected regions for a single observation (pointing)
517 
518  Parameters
519  ----------
520  obs : `~gammalib.GCTAObservation()`
521  CTA observation
522 
523  Returns
524  -------
525  regions : `~gammalib.GSkyRegions`
526  List of reflected regions
527  """
528  # Initialise list of reflected regions
529  regions = gammalib.GSkyRegions()
530 
531  # Get offset angle of source
532  pnt_dir = obs.pointing().dir()
533  offset = pnt_dir.dist_deg(self._src_dir)
534 
535  # Skip observation if it is too close to source
536  if self._src_reg.contains(pnt_dir):
537  msg = ' Skip because observation is pointed at %.3f deg from source'\
538  % (offset)
539  if self._srcshape == 'CIRCLE':
540  msg += ' (circle rad=%.3f).' % (self._rad)
541  self._log_string(gammalib.NORMAL, msg)
542  # ... otherwise
543  else:
544  posang = pnt_dir.posang_deg(self._src_dir)
545  if (self._srcshape == 'CIRCLE') or (self._srcshape == 'RECT'):
546 
547  # Determine number of background regions to skip
548  N_skip = self['bkgregskip'].integer()
549  N_lim = 1 + 2 * N_skip
550 
551  # Compute the angular separation of reflected regions wrt
552  # camera center. The factor 1.05 ensures background regions
553  # do not overlap due to numerical precision issues
554  alpha = 1.05 * self._compute_region_separation(pnt_dir)
555 
556  # Compute number of reflected regions by dividing the angular
557  # separation by 2 pi.
558  N = int(2.0 * math.pi / alpha)
559 
560  # If there are not enough reflected regions then skip the
561  # observation ...
562  if N < self['bkgregmin'].integer() + N_lim:
563  msg = ' Skip because the number %d of reflected regions '\
564  'for background estimation is smaller than '\
565  '"bkgregmin"=%d.' % (N-N_lim, self['bkgregmin'].integer())
566  self._log_string(gammalib.NORMAL, msg)
567 
568  # ... otherwise loop over position angle to create reflected
569  # regions
570  else:
571 
572  # Append reflected regions
573  alpha = 360.0 / N
574  dphi_max = 360.0 - alpha * (1 + N_skip)
575  dphi = alpha * (1 + N_skip)
576  while dphi <= dphi_max:
577  ctr_dir = pnt_dir.clone()
578  ctr_dir.rotate_deg(posang + dphi, offset)
579  if self._srcshape == 'CIRCLE':
580  region = gammalib.GSkyRegionCircle(ctr_dir, self._rad)
581  elif self._srcshape == 'RECT':
582  # Adjust the posang of the rectangle correspondingly
583  region = gammalib.GSkyRegionRectangle(ctr_dir,
584  self._reg_width,
585  self._reg_height,
586  self._reg_posang + dphi)
587  if self._has_exclusion:
588  if self._excl_reg.overlaps(region):
589 
590  # Signal region overlap
591  msg = ' Reflected region overlaps with '\
592  'exclusion region.'
593  self._log_string(gammalib.EXPLICIT, msg)
594 
595  # If region overlaps with exclusion region
596  # try to increment by 10% of angular step
597  dphi += 0.1 * alpha
598 
599  else:
600  regions.append(region)
601  dphi += alpha
602  else:
603  regions.append(region)
604  dphi += alpha
605 
606  # Check again that we have enough background regions
607  # now that we have checked for overlap with exclusion region
608  if regions.size() >= self['bkgregmin'].integer():
609  # Log number of reflected regions
610  msg = ' Use %d reflected regions.' % (regions.size())
611  self._log_string(gammalib.NORMAL, msg)
612  # Otherwise log observation skipped and return empty region container
613  else:
614  msg = ' Skip because the number %d of regions' \
615  'for background estimation not overlapping ' \
616  'with the exclusion region is smaller than ' \
617  '"bkgregmin"=%d.' % \
618  (regions.size(), self['bkgregmin'].integer())
619  self._log_string(gammalib.NORMAL, msg)
620  regions = gammalib.GSkyRegions()
621 
622  # Return reflected regions
623  return regions
624 
625  def _instrument_regions(self, obs, obs_off):
626  """
627  Compute background regions for Off observation
628 
629  Calculate background region in Off observation that corresponds to the
630  source region in the On observation in instrument coordinates
631 
632  Parameters
633  ----------
634  obs : `~gammalib.GCTAObservation()`
635  On CTA observation
636  obs_off : `~gammalib.GCTAObservation()`
637  Off CTA observation
638 
639  Returns
640  -------
641  regions : `~gammalib.GSkyRegions`
642  Container with background region
643  """
644  # Initialise region container
645  regions = gammalib.GSkyRegions()
646 
647  # Convert source position in On observation to instrument coordinates
648  instdir = obs.pointing().instdir(self._src_dir)
649 
650  # Convert instrument position to sky direction for Off observation
651  off_dir = obs_off.pointing().skydir(instdir)
652 
653  # Build region according to shape specified by user
654  # If circle
655  if self._srcshape == 'CIRCLE':
656  region = gammalib.GSkyRegionCircle(off_dir, self._rad)
657 
658  # ... otherwise if rectangle
659  elif self._srcshape == 'RECT':
660  # Instrument coordinates take sky direction as reference
661  # so no need to change the position angle
662  region = gammalib.GSkyRegionRectangle(off_dir,
663  self._reg_width,
664  self._reg_height,
665  self._reg_posang)
666 
667  # Check if background region overlaps with exclusion region
668  is_valid = True
669  if self._has_exclusion:
670  if self._excl_reg.overlaps(region):
671  # Signal region overlap
672  msg = ' Background region overlaps with exclusion region.'
673  self._log_string(gammalib.EXPLICIT, msg)
674  is_valid = False
675 
676  # If region is valid then append it to container
677  if is_valid:
678  regions.append(region)
679 
680  # Return
681  return regions
682 
683  def _set_models(self, results):
684  """
685  Set models for On/Off fitting
686 
687  The method does the following
688  - append "OnOff" to the instrument name of all background models
689  - fix all spatial and temporal parameters
690 
691  Parameters
692  ----------
693  results : list of dict
694  Result dictionaries
695 
696  Returns
697  -------
698  models : `~gammalib.GModels()`
699  Model container
700  """
701  # Write header
702  self._log_header1(gammalib.NORMAL, 'Set models')
703 
704  # Initialise model container
705  models = gammalib.GModels()
706 
707  # Initialies stacked model flag
708  has_stacked_model = False
709 
710  # Loop over all models in observation and append "OnOff" to instrument
711  # names
712  for model in self._models:
713 
714  # Initialise model usage
715  use_model = False
716 
717  # If model is a background model then check if it will be
718  # used
719  if 'GCTA' in model.classname():
720 
721  # Skip model if background model should not be used
722  if not self['use_model_bkg'].boolean():
723  self._log_string(gammalib.NORMAL, ' Skip "%s" model "%s" (%s)' % \
724  (model.instruments(), model.name(), model.ids()))
725  continue
726 
727  # Check if model corresponds to one of the relevant
728  # observations
729  for result in results:
730  if model.is_valid(result['instrument'], result['id']):
731  if result['bkg_reg'].size() > 0:
732  use_model = True
733  break
734 
735  # If stacked analysis is requested then just use for model
736  # and remove instrument ID
737  if self['stack'].boolean():
738 
739  # If there is already a model for stacked analysis then
740  # skip this one
741  if has_stacked_model:
742  msg = ' Skip "%s" model "%s" (%s). There is already ' \
743  'a model for stacked analysis.' % \
744  (model.instruments(), model.name(), model.ids())
745  self._log_string(gammalib.NORMAL, msg)
746  continue
747 
748  # ... otherwise use model for stacked analysis
749  else:
750  has_stacked_model = True
751  use_model = True
752  model.ids('')
753 
754  # Append "OnOff" to instrument name
755  model.instruments(model.instruments()+'OnOff')
756 
757  # ... otherwise, if model is not a background model then use it
758  else:
759  use_model = True
760 
761  # If model is relevant then append it now to the model
762  # container
763  if use_model:
764 
765  # Log model usage
766  self._log_string(gammalib.NORMAL, ' Use "%s" model "%s" (%s)' % \
767  (model.instruments(), model.name(), model.ids()))
768 
769  # Append model to container
770  models.append(model)
771 
772  # ... otherwise signal that model is skipped
773  else:
774  self._log_string(gammalib.NORMAL, ' Skip "%s" model "%s" (%s)' % \
775  (model.instruments(), model.name(), model.ids()))
776 
777  # Return model container
778  return models
779 
780  def _set_statistic(self, obs):
781  """
782  Set statistic for observation
783 
784  If the "use_model_bkg" is true then set statistic to "cstat",
785  otherwise set it to "wstat"
786 
787  Parameters
788  ----------
789  obs : `~gammalib.GObservation()`
790  Observation
791 
792  Returns
793  -------
794  obs : `~gammalib.GObservation()`
795  Observation
796  """
797  # Set statistic according to background model usage
798  if self['use_model_bkg'].boolean():
799  obs.statistic('cstat')
800  else:
801  obs.statistic('wstat')
802 
803  # Return observation
804  return obs
805 
806  def _etrue_ebounds(self):
807  """
808  Set true energy boundaries
809 
810  Returns
811  -------
812  ebounds : `~gammalib.GEbounds()`
813  True energy boundaries
814  """
815  # Determine minimum and maximum energies
816  emin = self._ebounds.emin() * 0.5
817  emax = self._ebounds.emax() * 1.2
818  if emin.TeV() < self['etruemin'].real():
819  emin = gammalib.GEnergy(self['etruemin'].real(), 'TeV')
820  if emax.TeV() > self['etruemax'].real():
821  emax = gammalib.GEnergy(self['etruemax'].real(), 'TeV')
822 
823  # Determine number of energy bins
824  n_decades = (emax.log10TeV() - emin.log10TeV())
825  n_bins = int(n_decades * float(self['etruebins'].integer()) + 0.5)
826  if n_bins < 1:
827  n_bins = 1
828 
829  # Set energy boundaries
830  ebounds = gammalib.GEbounds(n_bins, emin, emax)
831 
832  # Write header
833  self._log_header1(gammalib.TERSE, 'True energy binning')
834 
835  # Log true energy bins
836  for i in range(ebounds.size()):
837  value = '%s - %s' % (str(ebounds.emin(i)), str(ebounds.emax(i)))
838  self._log_value(gammalib.TERSE, 'Bin %d' % (i+1), value)
839 
840  # Return energy boundaries
841  return ebounds
842 
843  def _set_background_regions(self, obs, obs_off=None):
844  """
845  Set background regions for an observation
846 
847  Parameters
848  ----------
849  obs : `~gammalib.GCTAObservation()`
850  CTA observation
851 
852  Returns
853  -------
854  regions : `~gammalib.GSkyRegions()`
855  Background regions
856  """
857  # Initialise empty background regions for this observation
858  bkg_reg = gammalib.GSkyRegions()
859 
860  # If reflected background is requested then create reflected
861  # background regions
862  if self['bkgmethod'].string() == 'REFLECTED':
863  bkg_reg = self._reflected_regions(obs)
864 
865  # ... otherwise if custom background is requested then get the
866  # background regions from the observation. We use a copy here since
867  # otherwise the background regions go out of scope once the observations
868  # are replaced by the On/Off observations.
869  elif self['bkgmethod'].string() == 'CUSTOM':
870  bkg_reg = obs.off_regions().copy()
871 
872  # ... otherwise if dedicated Off runs are use then
873  # use background region that correspond to the same instrument coordinates
874  if self['bkgmethod'].string() == 'OFF':
875  bkg_reg = self._instrument_regions(obs,obs_off)
876 
877  # Return background regions
878  return bkg_reg
879 
880  def _process_observation(self,i):
881  """
882  Generate On/Off spectra for individual observation
883 
884  Parameters
885  ----------
886  i : int
887  Observation number
888 
889  Returns
890  -------
891  result : dict
892  On/Off spectra, background regions, observation id
893  """
894  # Retrieve observation from container
895  onoff = None
896  bkg_reg = None
897  obs = self.obs()[i]
898 
899  # Retrieve dedicated Off observation if it exists
900  if not self.obs_off().is_empty():
901  obs_off = self.obs_off()[i]
902  # Otherwise use the same as On
903  else:
904  obs_off = self.obs()[i]
905 
906  # Log header
907  self._log_header3(gammalib.NORMAL,'%s observation "%s"' % \
908  (obs.instrument(), obs.id()))
909 
910  # Skip non CTA observations
911  if obs.classname() != 'GCTAObservation':
912  self._log_string(gammalib.NORMAL, ' Skip because not a "GCTAObservation"')
913 
914  # Otherwise calculate On/Off spectra
915  else:
916  # Get background model usage flag and log flag
917  use_model_bkg = self['use_model_bkg'].boolean()
918  if use_model_bkg:
919  msg = ' Use background model.'
920  else:
921  msg = ' Background model not used, assume constant backround rate.'
922  self._log_string(gammalib.NORMAL, msg)
923 
924  # Get offset angle of source
925  pnt_dir = obs.pointing().dir()
926  offset = pnt_dir.dist_deg(self._src_dir)
927 
928  # Skip observation if it is pointed too far from the source
929  if offset >= self['maxoffset'].real():
930  msg = ' Skip because observation is pointed at %.3f deg >= ' \
931  '"maxoffset=%.3f" from source.' \
932  % (offset, self['maxoffset'].real())
933  self._log_string(gammalib.NORMAL, msg)
934  # ... otherwise continue to process
935  else:
936 
937  # Set background regions for this observation
938  bkg_reg = self._set_background_regions(obs,obs_off)
939 
940  # If there are any background regions then create On/Off observation
941  # and append it to the output container
942  if bkg_reg.size() >= 0:
943 
944  # Create On/Off observation
945  onoff = gammalib.GCTAOnOffObservation(obs, obs_off,
946  self._models,
947  self._srcname,
948  self._etruebounds,
949  self._ebounds,
950  self._src_reg,
951  bkg_reg,
952  use_model_bkg)
953 
954  # Set On/Off observation ID
955  onoff.id(obs.id())
956 
957  # Otherwise log observation skipped
958  else:
959  msg = ' Skip because no valid Off regions could be determined'
960  self._log_string(gammalib.NORMAL, msg)
961 
962  # Construct dictionary with results
963  result = {'onoff' : onoff,
964  'bkg_reg' : bkg_reg,
965  'instrument': obs.instrument(),
966  'id' : obs.id()}
967 
968  # Return results
969  return result
970 
971  def _unpack_result(self, outobs, result):
972  """
973  Unpack result from calculation of On/Off regions
974 
975  Parameters
976  ----------
977  outobs : `~gammalib.GObservations`
978  Observation container
979  result : dict
980  On/Off spectra, background regions, observation id
981 
982  Returns
983  -------
984  outobs : `~gammalib.GObservations`
985  Observation container with result appended
986  """
987  # Continue only if result is valid
988  if result['onoff'] != None:
989 
990  # If the results contain an On/Off observation
991  if result['onoff'].classname() == 'GCTAOnOffObservation':
992 
993  # Set statistic according to background model usage
994  obs = self._set_statistic(result['onoff'])
995 
996  # Append observation to observation container
997  outobs.append(obs)
998 
999  # Append background regions
1000  self._bkg_regs.append({'regions': result['bkg_reg'],
1001  'id': result['id']})
1002 
1003  # Return observation container
1004  return outobs
1005 
1006 
1007  # Public methods
1008  def process(self):
1009  """
1010  Process the script
1011  """
1012  # Get parameters
1013  self._get_parameters()
1014 
1015  # Write observation into logger
1016  self._log_observations(gammalib.NORMAL, self.obs(), 'Observation')
1017  if not self.obs_off().is_empty():
1018  self._log_observations(gammalib.NORMAL, self._obs_off, 'Off Observation')
1019 
1020  # Set true energy bins
1021  self._etruebounds = self._etrue_ebounds()
1022 
1023  # Write header
1024  self._log_header1(gammalib.TERSE, 'Spectral binning')
1025 
1026  # Log reconstructed energy bins
1027  for i in range(self._ebounds.size()):
1028  value = '%s - %s' % (str(self._ebounds.emin(i)),
1029  str(self._ebounds.emax(i)))
1030  self._log_value(gammalib.TERSE, 'Bin %d' % (i+1), value)
1031 
1032  # Write header
1033  self._log_header1(gammalib.NORMAL,
1034  'Generation of source and background spectra')
1035 
1036  # Initialise run variables
1037  outobs = gammalib.GObservations()
1038  self._bkg_regs = []
1039  results = []
1040 
1041  # If there is more than one observation and we use multiprocessing
1042  if self._nthreads > 1 and self.obs().size() > 1:
1043 
1044  # Compute observations
1045  args = [(self, '_process_observation', i)
1046  for i in range(self.obs().size())]
1047  poolresults = mputils.process(self._nthreads, mputils.mpfunc, args)
1048 
1049  # Construct results
1050  for i in range(self.obs().size()):
1051  result = poolresults[i][0]
1052  outobs = self._unpack_result(outobs, result)
1053  results.append(result)
1054  self._log_string(gammalib.TERSE, poolresults[i][1]['log'], False)
1055 
1056  # Otherwise, loop through observations and generate pha, arf, rmf files
1057  else:
1058  for i in range(self.obs().size()):
1059  # Process individual observation
1060  result = self._process_observation(i)
1061  outobs = self._unpack_result(outobs, result)
1062  results.append(result)
1063 
1064  # Stack observations
1065  if outobs.size() > 1 and self['stack'].boolean():
1066 
1067  # Write header
1068  self._log_header1(gammalib.NORMAL, 'Stacking %d observations' %
1069  (outobs.size()))
1070 
1071  # Stack observations
1072  stacked_obs = gammalib.GCTAOnOffObservation(outobs)
1073 
1074  # Set statistic according to background model usage
1075  stacked_obs = self._set_statistic(stacked_obs)
1076 
1077  # Put stacked observations in output container
1078  outobs = gammalib.GObservations()
1079  outobs.append(stacked_obs)
1080 
1081  # Create models that allow On/Off fitting
1082  models = self._set_models(results)
1083 
1084  # Set models in output container
1085  outobs.models(models)
1086 
1087  # Set observation container
1088  self.obs(outobs)
1089 
1090  # Return
1091  return
1092 
1093  def save(self):
1094  """
1095  Save data
1096  """
1097  # Write header
1098  self._log_header1(gammalib.TERSE, 'Save data')
1099 
1100  # Get XML output filename, prefix and clobber
1101  outobs = self['outobs'].filename()
1102  outmodel = self['outmodel'].filename()
1103  prefix = self['prefix'].string()
1104  clobber = self['clobber'].boolean()
1105 
1106  # Loop over all observation in container
1107  for obs in self.obs():
1108 
1109  # Set filenames
1110  if self['stack'].boolean():
1111  onname = prefix + '_stacked_pha_on.fits'
1112  offname = prefix + '_stacked_pha_off.fits'
1113  arfname = prefix + '_stacked_arf.fits'
1114  rmfname = prefix + '_stacked_rmf.fits'
1115  elif self.obs().size() > 1:
1116  onname = prefix + '_%s_pha_on.fits' % (obs.id())
1117  offname = prefix + '_%s_pha_off.fits' % (obs.id())
1118  arfname = prefix + '_%s_arf.fits' % (obs.id())
1119  rmfname = prefix + '_%s_rmf.fits' % (obs.id())
1120  else:
1121  onname = prefix + '_pha_on.fits'
1122  offname = prefix + '_pha_off.fits'
1123  arfname = prefix + '_arf.fits'
1124  rmfname = prefix + '_rmf.fits'
1125 
1126  # Set background and response file names in On spectrum
1127  obs.on_spec().backfile(offname)
1128  obs.on_spec().respfile(rmfname)
1129  obs.on_spec().ancrfile(arfname)
1130 
1131  # Save files
1132  obs.on_spec().save(onname, clobber)
1133  obs.off_spec().save(offname, clobber)
1134  obs.arf().save(arfname, clobber)
1135  obs.rmf().save(rmfname, clobber)
1136 
1137  # Stamp files
1138  self._stamp(onname)
1139  self._stamp(offname)
1140  self._stamp(arfname)
1141  self._stamp(rmfname)
1142 
1143  # Log file names
1144  self._log_value(gammalib.NORMAL, 'PHA on file', onname)
1145  self._log_value(gammalib.NORMAL, 'PHA off file', offname)
1146  self._log_value(gammalib.NORMAL, 'ARF file', arfname)
1147  self._log_value(gammalib.NORMAL, 'RMF file', rmfname)
1148 
1149  # Save observation definition XML file
1150  self.obs().save(outobs)
1151 
1152  # Save model definition XML file
1153  self.obs().models().save(outmodel)
1154 
1155  # Log file names
1156  self._log_value(gammalib.NORMAL, 'Obs. definition XML file', outobs.url())
1157  self._log_value(gammalib.NORMAL, 'Model definition XML file', outmodel.url())
1158 
1159  # Save ds9 On region file
1160  regname = prefix + '_on.reg'
1161  self._src_reg.save(regname)
1162  self._log_value(gammalib.NORMAL, 'On region file', regname)
1163 
1164  # Save ds9 Off region files
1165  for bkg_reg in self._bkg_regs:
1166 
1167  # Set filename
1168  if len(self._bkg_regs) > 1:
1169  regname = prefix + '_%s_off.reg' % (bkg_reg['id'])
1170  else:
1171  regname = prefix + '_off.reg'
1172 
1173  # Save ds9 region file
1174  bkg_reg['regions'].save(regname)
1175 
1176  # Log file name
1177  self._log_value(gammalib.NORMAL, 'Off region file', regname)
1178 
1179  # Return
1180  return
1181 
1182  def exclusion_map(self, object=None):
1183  """
1184  Return and optionally set the exclusion regions map
1185 
1186  Parameters
1187  ----------
1188  object : `~gammalib.GSkyRegion` or `~gammalib.GSkyMap` or `~gammalib.GFilename`
1189  Exclusion regions object
1190 
1191  Returns
1192  -------
1193  region : `~gammalib.GSkyRegionMap`
1194  Exclusion regions map
1195  """
1196  # If a regions object is provided then set the regions ...
1197  if object is not None:
1198  self._excl_reg = gammalib.GSkyRegionMap(object)
1199 
1200  # Return
1201  return self._excl_reg
1202 
1203  def obs_off(self, obs=None):
1204  """
1205  Return and optionally set the Off observations
1206 
1207  Parameters
1208  ----------
1209  obs : `~gammalib.GCTAObservations`
1210  Off observations container
1211 
1212  Returns
1213  -------
1214  observation container : `~gammalib.GCTAObservations`
1215  Off observations container
1216  """
1217  # If an observation container is provided then set the Off observations ...
1218  if obs is not None:
1219  self._obs_off = obs
1220 
1221  # Return
1222  return self._obs_off
1223 
1224 
1225 # ======================== #
1226 # Main routine entry point #
1227 # ======================== #
1228 if __name__ == '__main__':
1229 
1230  # Create instance of application
1231  app = csphagen(sys.argv)
1232 
1233  # Execute application
1234  app.execute()