25 from cscripts
import mputils
26 from cscripts
import obsutils
32 class csscs(ctools.csobservation):
34 Spectral component separation script
45 List of IRAF command line parameter strings of the form
49 self._init_csobservation(self.__class__.__name__, ctools.__version__, argv)
65 Extend ctools.csobservation __getstate__ method
73 state = {
'base': ctools.csobservation.__getstate__(self),
86 Extend ctools.csobservation __setstate__ method
94 ctools.csobservation.__setstate__(self, state[
'base'])
98 self.
_fits = state[
'fits']
108 Get On/Off analysis parameters.
110 This is done here rather than in ctools.is_onoff because the exclusion
111 map needs to be handled differently as an automatic parameter but
112 queried only if not already set and we need to verify is not empty.
113 Also many parameters are already queried for all methods
115 TODO: verify if this can be made more uniform with other scripts
118 if (self.
_excl_reg_map is not None)
and (self._excl_reg_map.map().npix() > 0):
121 elif self[
'inexclusion'].is_valid():
125 msg =
'csscs in On/Off mode requires input exclusion region.'
126 raise RuntimeError(msg)
129 self[
'enumbins'].integer()
130 if self[
'bkgmethod'].string() ==
'REFLECTED':
131 self[
'bkgregmin'].integer()
132 self[
'bkgregskip'].integer()
133 self[
'maxoffset'].real()
134 self[
'etruemin'].real()
135 self[
'etruemax'].real()
136 self[
'etruebins'].integer()
143 Get parameters from parfile
146 if self.obs().is_empty():
147 self.obs(self._get_observations())
150 self._set_obs_statistic(gammalib.toupper(self[
'statistic'].string()))
157 for obs
in self.obs():
158 if obs.classname() ==
'GCTAObservation':
159 if obs.eventtype() ==
'CountsCube':
163 elif obs.classname() ==
'GCTAOnOffObservation':
165 n_cta = n_unbinned + n_binned + n_onoff
166 n_other = self.obs().size() - n_cta
169 if n_onoff > 0
or n_other > 0:
170 msg =
'On/Off or non-CTA observations found. ' \
171 'csscs does not support this type of observations.'
172 raise RuntimeError(msg)
175 elif n_unbinned > 0
and n_binned > 0:
176 msg =
'Mix of unbinned and binned CTA observations ' \
177 'found in observation container. csscs does not ' \
179 raise RuntimeError(msg)
182 if self.obs().models().size() == 0:
183 self.obs().models(self[
'inmodel'].filename())
188 self[
'coordsys'].string()
189 self[
'proj'].string()
190 self[
'nxpix'].integer()
191 self[
'nypix'].integer()
203 self.
_method = self[
'method'].string()
212 srcnames = self[
'srcnames'].string()
217 for s, name
in enumerate(self.
_srcnames):
223 src_in_model = src_in_model
and self.obs().models().contains(name)
224 if src_in_model ==
False:
225 msg =
'Not all target sources are present in input model.'
226 raise RuntimeError(msg)
235 for model
in self.obs().models():
236 if model.classname() ==
'GModelSky':
242 msg =
'Background gamma-ray sources found in the model. ' \
243 'On/Off analysis does not support this feature.'
244 raise RuntimeError(msg)
251 self[
'edisp'].boolean()
252 self[
'calc_ulim'].boolean()
253 self[
'calc_ts'].boolean()
254 self[
'fix_bkg'].boolean()
255 self[
'fix_srcs'].boolean()
258 if self._read_ahead():
259 self[
'outfile'].query()
262 self._log_parameters(gammalib.TERSE)
272 Create a gammalib.GSkyMap object to store the output
276 skymap : `~gammalib.GSkyMap`
280 skymap = gammalib.GSkyMap(self[
'proj'].string(), self[
'coordsys'].string(),
281 self[
'xref'].real(), self[
'yref'].real(),
282 -self[
'binsz'].real(), self[
'binsz'].real(),
283 self[
'nxpix'].integer(), self[
'nypix'].integer(),
291 Adjust model parameters
294 self._log_header1(gammalib.TERSE,
'Adjust model parameters')
297 for model
in self.obs().models():
303 self._log_header3(gammalib.EXPLICIT, model.name())
312 gammalib.GModelSky(model)
314 if model.spatial() ==
'DiffuseIsotropic':
319 msg =
' Setting spatial model to diffuse isotropic'
320 self._log_string(gammalib.EXPLICIT, msg)
321 model.spatial(gammalib.GModelSpatialDiffuseConst())
327 normpar = model.spectral()[0]
329 if par.name() == normpar.name():
333 self._log_string(gammalib.EXPLICIT,
' Fixing "' + par.name() +
'"')
341 if normpar.is_fixed():
342 self._log_string(gammalib.EXPLICIT,
' Freeing "' + normpar.name() +
'"')
346 if self[
'calc_ts'].boolean():
350 elif self[
'fix_bkg'].boolean()
and not model.classname() ==
'GModelSky':
352 if normpar.is_free():
353 self._log_string(gammalib.EXPLICIT,
' Fixing "' + normpar.name() +
'"')
358 elif self[
'fix_srcs'].boolean()
and model.classname() ==
'GModelSky':
360 if normpar.is_free():
361 self._log_string(gammalib.EXPLICIT,
' Fixing "' + normpar.name() +
'"')
369 Mask cube observation
373 obs : `gammalib.GCTAObservation`
374 Input observation of type cube
376 Right Ascension (deg)
384 new_obs : `gammalib.GCTAObservations`
385 Output observations with masked cubes
388 cubemask = ctools.ctcubemask(obs)
390 cubemask[
'dec'] = dec
391 cubemask[
'rad'] = rad
392 cubemask[
'regfile'] =
'NONE'
393 cubemask[
'emin'] = self[
'emin'].real()
394 cubemask[
'emax'] = self[
'emax'].real()
398 if self._logVerbose()
and self._logDebug():
399 cubemask[
'debug'] =
True
405 new_obs = cubemask.obs().copy()
416 obs : `gammalib.GCTAObservations`
417 Input observations of type event list
419 Right Ascension (deg)
427 new_obs : `gammalib.GCTAObservations`
428 Output observations with masked event lists
431 select = ctools.ctselect(obs)
432 select[
'usepnt'] =
False
436 select[
'emin'] = self[
'emin'].real()
437 select[
'emax'] = self[
'emax'].real()
438 select[
'tmin'] =
'INDEF'
439 select[
'tmax'] =
'INDEF'
443 if self._logVerbose()
and self._logDebug():
444 select[
'debug'] =
True
450 new_obs = select.obs().copy()
457 Create observations restricted to circular ROI
462 Right Ascension (deg)
470 new_obs : `~gammalib.GObservations`
471 Observations in circular ROI
478 new_obs = self.
_mask_cube(self.obs(), ra, dec, rad)
480 new_obs = obsutils.get_onoff_obs(self, self.obs(), nthreads=1,
489 Performs analysis over the region of interest corresponding to a
490 single pixel of output map
507 pixdir = outmap.inx2dir(inx)
509 dec = pixdir.dec_deg()
512 msg =
'Spatial bin (%d) centred on (R.A.,Dec) = (%f,%f) deg' % (inx, ra, dec)
513 self._log_header2(gammalib.TERSE, msg)
518 result[name] = {
'flux': 0.0,
524 self._log_header3(gammalib.NORMAL,
'Masking observations')
528 if masked_obs.size() > 0:
531 self._log_header3(gammalib.NORMAL,
'Performing fit in region')
532 like = ctools.ctlike(masked_obs)
533 like[
'edisp'] = self[
'edisp'].boolean()
538 if self._logVerbose()
and self._logDebug():
545 self._log_string(gammalib.EXPLICIT, str(like.obs().models()))
548 if like.obs().logL() != 0.0:
552 centre = gammalib.GSkyDir()
553 centre.radec_deg(ra, dec)
554 roi = gammalib.GSkyRegionCircle(centre, self[
'rad'].real())
557 emin = gammalib.GEnergy(self[
'emin'].real(),
'TeV')
558 emax = gammalib.GEnergy(self[
'emax'].real(),
'TeV')
564 source = like.obs().models()[name]
568 flux = source.spectral().
flux(emin, emax)
572 corr_factor = source.spatial().
flux(roi)
573 corr_factor /= gammalib.twopi * (1.0 - math.cos(math.radians(self[
'rad'].real())))
577 result[name][
'flux'] = flux
581 normpar = source.spectral()[0]
586 if normpar.value() > 0.:
587 flux_error = flux * normpar.error() / normpar.value()
590 result[name][
'flux_error'] = flux_error
593 if self[
'calc_ts'].boolean():
594 result[name][
'TS'] = source.ts()
597 if self[
'calc_ulim'].boolean():
600 self._log_header3(gammalib.NORMAL,
601 'Computing upper limit for source ' + name)
604 ulimit = ctools.ctulimit(like.obs())
605 ulimit[
'srcname'] = name
606 ulimit[
'emin'] = self[
'emin'].real()
607 ulimit[
'emax'] = self[
'emax'].real()
611 if self._logVerbose()
and self._logDebug():
612 ulimit[
'debug'] =
True
617 ulimit_value = ulimit.flux_ulimit()
619 ulimit_value *= corr_factor
620 result[name][
'ulimit'] = ulimit_value
622 self._log_string(gammalib.NORMAL,
'Upper limit '
623 'calculation failed.')
627 value =
'Likelihood is zero. Bin is skipped.'
628 self._log_value(gammalib.TERSE,
'(R.A.,Dec) = (%f,%f) deg' % (ra, dec), value)
632 value =
'Size of masked observations is zero. Bin is skipped.'
633 self._log_value(gammalib.TERSE,
'(R.A.,Dec) = (%f,%f) deg' % (ra, dec), value)
644 hdu : `~gammalib.GFitsHDU`
649 hdu.card(
'BUNIT',
'ph/cm2/s/sr',
'Photon flux')
652 hdu.card(
'E_MIN', self[
'emin'].real(),
'[TeV] Lower energy boundary')
653 hdu.card(
'E_MAX', self[
'emax'].real(),
'[TeV] Upper energy boundary')
654 hdu.card(
'EUNIT',
'TeV',
'Units for E_MIN and E_MAX')
657 hdu.card(
'ROISZ', self[
'rad'].real(),
658 '[deg] Size of ROI used for component separation')
661 hdu.card(
'METHOD', self.
_method,
'Analysis method')
668 hdu.card(
'EXCLMAP', self._excl_reg_name.url(),
'Exclusion map name')
671 hdu.card(
'BKGMODEL', self[
'use_model_bkg'].boolean(),
'Use background model?')
674 hdu.card(
'STAT', self[
'statistic'].string(),
'Fit statistic')
681 Fill FITS object to store the results
690 fits : `~gammalib.GFits`
694 fits = gammalib.GFits()
697 fits.append(gammalib.GFitsImageDouble())
700 for s, name
in enumerate(self.
_srcnames):
711 for inx
in range(fmap.npix()):
712 fmap[inx, 0] = results[inx][name][
'flux']
713 errmap[inx, 0] = results[inx][name][
'flux_error']
714 tsmap[inx, 0] = results[inx][name][
'TS']
715 ulmap[inx, 0] = results[inx][name][
'ulimit']
718 Name = gammalib.toupper(name)
721 fhdu = fmap.write(fits, Name +
' FLUX')
723 errhdu = errmap.write(fits, Name +
' FLUX ERROR')
727 if self[
'calc_ts'].boolean():
728 tshdu = tsmap.write(fits, Name +
' TS')
732 if self[
'calc_ulim'].boolean():
733 ulhdu = ulmap.write(fits, Name +
' FLUX UPPER LIMIT')
741 Fetches skymap from FITS container
752 skymap : `~gammalib.GSkyMap`
764 if self.
_fits !=
None:
770 Name = gammalib.toupper(name)
773 hduname = Name +
' ' + quantity
777 hdu = self.
_fits[hduname]
779 skymap = gammalib.GSkyMap(hdu)
782 msg =
'HDU "' + hduname +
'" not found in FITS container.'
783 raise RuntimeError(msg)
787 print(
'ERROR: Source "' + name +
'" not found in list of target sources.')
788 raise NameError(name)
802 self._log_observations(gammalib.NORMAL, self.obs(),
'Input observation')
814 args = [(self,
'_pixel_analysis', i)
815 for i
in range(npix)]
816 poolresults = mputils.process(self.
_nthreads, mputils.mpfunc, args)
820 for i
in range(npix):
821 results.append(poolresults[i][0])
822 self._log_string(gammalib.TERSE, poolresults[i][1][
'log'],
False)
827 for i
in range(npix):
834 self._stamp(self.
_fits)
841 Save maps to Fits flie
844 self._log_header1(gammalib.TERSE,
'Save maps')
847 if self.
_fits !=
None:
849 outfile = self[
'outfile'].filename()
852 self._log_value(gammalib.NORMAL,
'Output file', outfile.url())
855 self._fits.saveto(outfile, self[
'clobber'].boolean())
862 Return and optionally set the exclusion regions map
866 object : `~gammalib.GSkyRegion` or `~gammalib.GSkyMap` or `~gammalib.GFilename`
867 Exclusion regions object
871 region : `~gammalib.GSkyRegionMap`
872 Exclusion regions map
875 if object
is not None:
883 Return copy of the fits container
887 fits : `~gammalib.GFits`
888 FITS file containing the maps
891 return self._fits.copy()
904 skymap : `~gammalib.GSkyMap`
915 Return flux error skymap
924 skymap : `~gammalib.GSkyMap`
944 skymap : `~gammalib.GSkyMap`
948 if self[
'calc_ts'].boolean():
953 msg =
'TS computation not requested. Change user parameter ' \
954 '"calc_ts" to True to obtain TS.'
955 raise RuntimeError(msg)
962 Return flux upper limit skymap
971 skymap : `~gammalib.GSkyMap`
972 Flux upper limit skymap
975 if self[
'calc_ulim'].boolean():
976 skymap = self.
_get_skymap(name,
'FLUX UPPER LIMIT')
980 msg =
'Upper limit computation not requested. Change user parameter ' \
981 '"calc_ulimit" to True to obtain upper limit.'
982 raise RuntimeError(msg)
991 if __name__ ==
'__main__':
def _get_onoff_parameters