#! /usr/bin/env python # ========================================================================== # Computes the PHA spectra for source/background and ARF/RMF files using the # reflected region method # # 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 # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program. If not, see . # # ========================================================================== import gammalib import ctools import math import sys # =============== # # csfindobs class # # =============== # class csphagen(ctools.csobservation): """ Generate PHA, ARF and RMF files for classical IACT spectral analysis """ # Constructor def __init__(self, *argv): """ Constructor """ # Initialise application by calling the appropriate class constructor self._init_csobservation(self.__class__.__name__, ctools.__version__, argv) # Initialise other variables self._ebounds = gammalib.GEbounds() self._src_dir = gammalib.GSkyDir() self._src_reg = gammalib.GSkyRegions() self._bkg_regs = [] self._excl_reg = None self._has_exclusion = False self._srcshape = '' self._rad = 0.0 # Return return # Private methods def _query_src_direction(self): """ Set up the source direction parameter. Query relevant parameters. """ self._src_dir = gammalib.GSkyDir() coordsys = self['coordsys'].string() if coordsys == 'CEL': ra = self['ra'].real() dec = self['dec'].real() self._src_dir.radec_deg(ra, dec) elif coordsys == 'GAL': 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)) # 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.classname() == 'GCTAObservation': 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._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()) self._has_exclusion = True else: self._has_exclusion = False # Query remaining parametersStacking self['stack'].boolean() # Query ahead output parameters if (self._read_ahead()): self['outobs'].string() self['prefix'].string() # If there are no observations in container then get them from the # parameter file 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) # Return return def _reflected_regions(self, obs): """ Calculate list of reflected regions for a single observation (pointing) Parameters ---------- obs : `~gammalib.GCTAObservation()` CTA observation Returns ------- regions : `~gammalib.GSkyRegions` List of reflected regions """ # Initialise list of reflected regions regions = gammalib.GSkyRegions() # Get offset angle of source pnt_dir = obs.pointing().dir() offset = pnt_dir.dist_deg(self._src_dir) # If ... if offset <= self._rad or offset >= self['maxoffset'].real(): self._log_string(gammalib.EXPLICIT, 'Observation %s pointed at %.3f ' 'deg from source' % ((obs.id(), offset))) # ... otherwise else: posang = pnt_dir.posang_deg(self._src_dir) if self._srcshape == 'CIRCLE': # angular separation of reflected regions wrt camera center # and number alpha = 1.05 * 2.0 * self._rad / offset # 1.05 ensures background regions do not overlap due to # numerical precision issues N = int(2.0 * math.pi / alpha) if N < self['bkgregmin'].integer() + 3: self._log_string(gammalib.EXPLICIT, 'Observation %s: ' 'insufficient regions for background ' 'estimation' % (obs.id())) # Otherwise loop over position angle to create reflected # regions else: alpha = 360.0 / N for s in range(2, N - 1): dphi = s * alpha ctr_dir = pnt_dir.clone() ctr_dir.rotate_deg(posang + dphi, offset) region = gammalib.GSkyRegionCircle(ctr_dir, self._rad) if self._has_exclusion: if self._excl_reg.overlaps(region): self._log_string(gammalib.VERBOSE, 'Observation ' '%s: reflected region overlaps ' 'with exclusion region' % (obs.id())) else: regions.append(region) else: regions.append(region) # Return reflected regions return regions def _set_models(self, obs): """ Set models in observation container The method replaces all "CTA", "HESS", "VERITAS", and "MAGIC" background models by "CTAOnOff" background models. Parameters ---------- obs : `~gammalib.GObservations()` Observation container Returns ------- obs : `~gammalib.GObservations()` Observation container """ # Initialise model container models = gammalib.GModels() # Loop over all models and replace CTA/HESS/VERITAS/MAGIC background # models by CTAOnOff background model for model in self.obs().models(): if 'GCTA' in model.classname(): if 'CTA' in model.instruments() or \ 'HESS' in model.instruments() or \ 'VERITAS' in model.instruments() or \ 'MAGIC' in model.instruments(): model.instruments('CTAOnOff') models.append(model) # Append model to observation container obs.models(models) # Return observation container return obs def _etrue_ebounds(self): """ Set true energy boundaries Returns ------- ebounds : `~gammalib.GEbounds()` True energy boundaries """ # Determine minimum and maximum energies emin = self._ebounds.emin() * 0.5 emax = self._ebounds.emax() * 1.2 if emin.TeV() < self['etruemin'].real(): emin = gammalib.GEnergy(self['etruemin'].real(), 'TeV') if emax.TeV() > self['etruemax'].real(): emax = gammalib.GEnergy(self['etruemax'].real(), 'TeV') # Determine number of energy bins n_decades = (emax.log10TeV() - emin.log10TeV()) n_bins = int(n_decades * float(self['etruebins'].integer())) + 1 # Set energy boundaries ebounds = gammalib.GEbounds(n_bins, emin, emax) # Write header self._log_header1(gammalib.TERSE, 'True energy binning') # Log true energy bins for i in range(ebounds.size()): value = '%s - %s' % (str(ebounds.emin(i)), str(ebounds.emax(i))) self._log_value(gammalib.TERSE, 'Bin %d' % (i+1), value) # 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): """ Run the script """ # Switch screen logging on in debug mode if self._logDebug(): self._log.cout(True) # Get parameters self._get_parameters() # Write observation into logger self._log_observations(gammalib.NORMAL, self.obs(), 'Observation') # Set true energy bins etrue_ebounds = self._etrue_ebounds() # Write header self._log_header1(gammalib.TERSE, 'Spectral binning') # Log reconstructed energy bins for i in range(self._ebounds.size()): value = '%s - %s' % (str(self._ebounds.emin(i)), str(self._ebounds.emax(i))) self._log_value(gammalib.TERSE, 'Bin %d' % (i+1), value) # Write header self._log_header1(gammalib.NORMAL, 'Generation of source and background spectra') # Initialise run variables outobs = gammalib.GObservations() self._bkg_regs = [] # Loop through observations and generate pha, arf, rmf files for obs in self.obs(): # Skip non CTA observations #if obs.instrument() != 'CTA': if obs.classname() != 'GCTAObservation': 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())) # Set background regions for this observation bkg_reg = self._set_background_regions(obs) # 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, etrue_ebounds, self._ebounds, self._src_reg, bkg_reg) onoff.id(obs.id()) outobs.append(onoff) self._bkg_regs.append({'regions': bkg_reg, 'id': obs.id()}) else: self._log_string(gammalib.NORMAL, 'Observation %s not included ' 'in spectra generation' % (obs.id())) # Stack observations if outobs.size() > 1 and self['stack'].boolean() == True: # Write header self._log_header1(gammalib.NORMAL, 'Stacking %d observations' % (outobs.size())) # Stack observations stacked_obs = gammalib.GCTAOnOffObservation(outobs) # Put stacked observations in output container outobs = gammalib.GObservations() outobs.append(stacked_obs) # Set models in output container outobs = self._set_models(outobs) # Set observation container self.obs(outobs) # Return return def save(self): """ Save data """ # Write header self._log_header1(gammalib.TERSE, 'Save data') # Get XML output filename, prefix and clobber outobs = self['outobs'].string() prefix = self['prefix'].string() clobber = self['clobber'].boolean() # Loop over all observation in container for obs in self.obs(): # Set filenames if self['stack'].boolean(): onname = prefix + '_stacked_pha_on.fits' offname = prefix + '_stacked_pha_off.fits' arfname = prefix + '_stacked_arf.fits' rmfname = prefix + '_stacked_rmf.fits' elif self.obs().size() > 1: onname = prefix + '_%s_pha_on.fits' % (obs.id()) offname = prefix + '_%s_pha_off.fits' % (obs.id()) arfname = prefix + '_%s_arf.fits' % (obs.id()) rmfname = prefix + '_%s_rmf.fits' % (obs.id()) else: onname = prefix + '_pha_on.fits' offname = prefix + '_pha_off.fits' arfname = prefix + '_arf.fits' rmfname = prefix + '_rmf.fits' # Save files obs.on_spec().save(onname, clobber) obs.off_spec().save(offname, clobber) obs.arf().save(arfname, clobber) obs.rmf().save(rmfname, clobber) # Log file names self._log_value(gammalib.NORMAL, 'PHA on file', onname) self._log_value(gammalib.NORMAL, 'PHA off file', offname) self._log_value(gammalib.NORMAL, 'ARF file', arfname) self._log_value(gammalib.NORMAL, 'RMF file', arfname) # Save observation definition XML file self.obs().save(outobs) # Log file name self._log_value(gammalib.NORMAL, 'Obs. definition XML file', outobs) # Save ds9 On region file regname = prefix + '_on.reg' self._src_reg.save(regname) self._log_value(gammalib.NORMAL, 'On region file', regname) # Save ds9 Off region files for bkg_reg in self._bkg_regs: # Set filename if len(self._bkg_regs) > 1: regname = prefix + '_%s_off.reg' % (bkg_reg['id']) else: regname = prefix + '_off.reg' # Save ds9 region file bkg_reg['regions'].save(regname) # Log file name self._log_value(gammalib.NORMAL, 'Off region file', regname) # Return return # ======================== # # Main routine entry point # # ======================== # if __name__ == '__main__': # Create instance of application app = csphagen(sys.argv) # Execute application app.execute()