diff --git a/cscripts/csspec.py b/cscripts/csspec.py index 3ed8f0b..de0deb6 100755 --- a/cscripts/csspec.py +++ b/cscripts/csspec.py @@ -2,7 +2,7 @@ # ========================================================================== # Generates a spectrum. # -# Copyright (C) 2014-2017 Michael Mayer +# Copyright (C) 2014-2018 Michael Mayer # # 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 @@ -100,26 +100,29 @@ class csspec(ctools.csobservation): n_other = self.obs().size() - n_cta # If spectrum method is not "NODES" then set spectrum method and - # script mode according to type of CTA observations + # script mode according to type of observations if self._method != 'NODES': - if n_unbinned == 0 and n_binned != 0 and n_onoff == 0: - self._binned_mode = True - self._method = 'SLICE' - elif n_unbinned == 0 and n_binned == 0 and n_onoff != 0: - self._onoff_mode = True - self._method = 'SLICE' - elif n_unbinned == 0 and n_binned != 0 and n_onoff != 0: - msg = 'Mix of binned and On/Off CTA observations found in ' \ - 'observation container. csscript does not support this ' \ - 'mix.' - raise RuntimeError(msg) - elif n_unbinned != 0 and (n_binned != 0 or n_onoff != 0): - msg = 'Mix of unbinned and binned or On/Off CTA observations ' \ - 'found in observation container. csscript does not support ' \ - 'this mix.' - raise RuntimeError(msg) - elif n_unbinned != 0: - self._method = 'SLICE' + if n_other > 0: + self._method = 'NODES' + else: + if n_unbinned == 0 and n_binned != 0 and n_onoff == 0: + self._binned_mode = True + self._method = 'SLICE' + elif n_unbinned == 0 and n_binned == 0 and n_onoff != 0: + self._onoff_mode = True + self._method = 'SLICE' + elif n_unbinned == 0 and n_binned != 0 and n_onoff != 0: + msg = 'Mix of binned and On/Off CTA observations found ' \ + 'in observation container. csscript does not support ' \ + 'this mix.' + raise RuntimeError(msg) + elif n_unbinned != 0 and (n_binned != 0 or n_onoff != 0): + msg = 'Mix of unbinned and binned or On/Off CTA observations ' \ + 'found in observation container. csscript does not ' \ + 'support this mix.' + raise RuntimeError(msg) + elif n_unbinned != 0: + self._method = 'SLICE' # Set ebounds self._set_ebounds() @@ -148,10 +151,9 @@ class csspec(ctools.csobservation): # If there are a mix of CTA and non-CTA observations and the method # is 'SLICE' then log a warning that non-CTA observations will be # ignored + warning = False if n_cta > 0 and n_other > 0 and self._method == 'SLICE': - self._log_string(gammalib.TERSE, '\nWARNING: Only CTA observation ' - 'can be handled with the "SLICE" method, all ' - 'non-CTA observation will be ignored.') + warning = True # If there are only non-CTA observations and the method is 'SLICE' # then stop now @@ -167,6 +169,12 @@ class csspec(ctools.csobservation): # Log selected spectrum method self._log_value(gammalib.TERSE, 'Selected spectrum method', self._method) + # Signal warning + if warning: + self._log_string(gammalib.TERSE, ' WARNING: Only CTA observation ' + 'can be handled with the "SLICE" method, all ' + 'non-CTA observation will be ignored.') + # Return return diff --git a/doc/source/users/tutorials/howto/howto_combine_instruments.png b/doc/source/users/tutorials/howto/howto_combine_instruments.png new file mode 100644 index 0000000..cdd7ead --- /dev/null +++ b/doc/source/users/tutorials/howto/howto_combine_instruments.png diff --git a/doc/source/users/tutorials/howto/howto_combine_instruments.rst b/doc/source/users/tutorials/howto/howto_combine_instruments.rst new file mode 100644 index 0000000..6b33c28 --- /dev/null +++ b/doc/source/users/tutorials/howto/howto_combine_instruments.rst @@ -0,0 +1,210 @@ +.. _howto_combine_instrument: + +How to jointly analyse data from different instruments? +------------------------------------------------------- + + .. admonition:: What you will learn + + You will learn how to **jointly analyse data from different instruments**. + +In this tutorial you will learn how to jointly analyse data observed with +Fermi-LAT and simulated data for CTA for the Vela pulsar and pulsar wind +nebula. + +Prepare Fermi-LAT data +^^^^^^^^^^^^^^^^^^^^^^ + +The preparation of Fermi-LAT data will follow the one +`described in this section `_. + + +Prepare CTA data +^^^^^^^^^^^^^^^^ + +For CTA the simulations of the Galactic Plane Scan from the first Data Challenge +will be used, but if you do not have these data at hand you can create your own +simulated CTA dataset by +`following the description here `_. + +First, bin the simulated data into a counts cube + +.. code-block:: bash + + $ ctbin + Input event list or observation definition XML file [events.fits] obs.xml + First coordinate of image center in degrees (RA or galactic l) (0-360) [83.63] 128.838 + Second coordinate of image center in degrees (DEC or galactic b) (-90-90) [22.51] -45.178 + Projection method (AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN) [CAR] + Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL] + Image scale (in degrees/pixel) [0.02] + Size of the X axis in pixels [200] + Size of the Y axis in pixels [200] + Algorithm for defining energy bins (FILE|LIN|LOG) [LOG] + Lower energy limit (TeV) [0.1] 0.03 + Upper energy limit (TeV) [100.0] 150.0 + Number of energy bins (1-200) [20] 40 + Output counts cube file [cntcube.fits] + +Then compute the response cubes as follows + +.. code-block:: bash + + $ ctexpcube + Input event list or observation definition XML file [NONE] obs.xml + Input counts cube file to extract exposure cube definition [NONE] cntcube.fits + Output exposure cube file [expcube.fits] + +.. code-block:: bash + + $ ctpsfcube + Input event list or observation definition XML file [NONE] obs.xml + Input counts cube file to extract PSF cube definition [NONE] + First coordinate of image center in degrees (RA or galactic l) (0-360) [83.63] 128.838 + Second coordinate of image center in degrees (DEC or galactic b) (-90-90) [22.51] -45.178 + Projection method (AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN) [CAR] + Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL] + Image scale (in degrees/pixel) [1.0] + Size of the X axis in pixels [10] + Size of the Y axis in pixels [10] + Algorithm for defining energy bins (FILE|LIN|LOG) [LOG] + Lower energy limit (TeV) [0.1] 0.03 + Upper energy limit (TeV) [100.0] 150.0 + Number of energy bins (1-200) [20] 40 + Output PSF cube file [psfcube.fits] + +.. code-block:: bash + + $ ctbkgcube + Input event list or observation definition XML file [NONE] obs.xml + Input counts cube file to extract background cube definition [NONE] cntcube.fits + Input model definition XML file [NONE] $CTOOLS/share/models/bkg_irf.xml + Output background cube file [bkgcube.fits] + Output model definition XML file [NONE] bkgcube.xml + + +Combining the observations +^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Now you have all the data and hand. You have to create an +:ref:`observation definition file ` +to combine the data for an analysis: + +.. code-block:: xml + + + + + + + + + + + + + + + + + +The +:ref:`observation definition file ` +contains two observations of the Vela pulsar, a first done with CTA and a +second done with Fermi-LAT. The ``instrument`` attribute distinguishes between +both instruments. + + +Generate a spectral energy distribution +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Before being able to generate a spectral energy distribution (SED) you have +to define a +:ref:`model definition file ` +that will be used to model the events for both observations. +The model that will be used in this analysis is shown below: + +.. code-block:: xml + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +The model contains a point source located at the position of the Vela pulsar +with a power law spectrum. In addition, it contains two ``DiffuseSource`` +components that are only applied for Fermi-LAT observations, which is indicated +by their ``instrument="LAT"`` attribute. Both components model the diffuse +background that prevails at GeV energies. Finally, the model contains a +``CTACubeBackground`` component that applies to CTA. You may have noticed that +the point source is the only component that has no ``instrument`` attribute, +meaning that this component applies to both instruments. + +Now you are ready to generate the spectral energy distribution for the +combined data set. You do this using the :ref:`csspec` script as follows + +.. code-block:: bash + + $ csspec + Input event list, counts cube, or observation definition XML file [events.fits] obs.xml + Input model definition XML file [$CTOOLS/share/models/crab.xml] models.xml + Source name [Crab] Vela + Spectrum generation method (SLICE|NODES|AUTO) [AUTO] + Algorithm for defining energy bins (FILE|LIN|LOG) [LOG] + Start value for first energy bin in TeV [0.1] 0.0001 + Stop value for last energy bin in TeV [100.0] 150.0 + Number of energy bins (1-200) [20] + Output spectrum file [spectrum.fits] + +This will generate a logarithmically spaced spectrum composed of 20 energy bins +comprised within 100 MeV and 150 TeV. The :ref:`csspec` tool is run in the ``AUTO`` +mode, which for different instruments corresponds to the ``NODES`` method. +The resulting spectral energy distribution is shown below: + +.. figure:: howto_combine_instruments.png + :width: 600px + :align: center + + *Vela spectrum derived using csspec from real Fermi-LAT and simulated CTA data* + +.. note:: + The figure was created using the ``show_spectrum.py`` script that is + located in the ctools example folder. The example script requires the + `matplotlib `_ Python module for display. + You may reproduce the plot by typing + + .. code-block:: bash + + $ $CTOOLS/share/examples/python/show_spectrum.py spectrum.fits diff --git a/doc/source/users/tutorials/howto/index.rst b/doc/source/users/tutorials/howto/index.rst index cbf32f7..fa0780d 100644 --- a/doc/source/users/tutorials/howto/index.rst +++ b/doc/source/users/tutorials/howto/index.rst @@ -46,3 +46,4 @@ Analysis of data from other instruments How to analyse Fermi-LAT data with ctools? How to analyse COMPTEL data with ctools? How to analyse H.E.S.S. data with ctools? + howto_combine_instruments