From aa443cc0124bc21707d6c897fb281f7a643cb6ed Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Sat, 9 Sep 2017 10:18:23 +0200 Subject: [PATCH] Add COMPTEL DRI tutorial (#2187) --- doc/source/users/tutorials/howto/comptel/howto_comptel_1.8MeV_tsmap.png | Bin 0 -> 89738 bytes doc/source/users/tutorials/howto/comptel/howto_comptel_1.8MeV_tsmap_5vp.png | Bin 0 -> 60539 bytes doc/source/users/tutorials/howto/comptel/howto_comptel_dri.rst | 300 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ doc/source/users/tutorials/howto/comptel/index.rst | 1 + 4 files changed, 301 insertions(+), 0 deletions(-) create mode 100644 doc/source/users/tutorials/howto/comptel/howto_comptel_1.8MeV_tsmap.png create mode 100644 doc/source/users/tutorials/howto/comptel/howto_comptel_1.8MeV_tsmap_5vp.png create mode 100644 doc/source/users/tutorials/howto/comptel/howto_comptel_dri.rst diff --git a/doc/source/users/tutorials/howto/comptel/howto_comptel_1.8MeV_tsmap.png b/doc/source/users/tutorials/howto/comptel/howto_comptel_1.8MeV_tsmap.png new file mode 100644 index 0000000..94528bf Binary files /dev/null and b/doc/source/users/tutorials/howto/comptel/howto_comptel_1.8MeV_tsmap.png differ diff --git a/doc/source/users/tutorials/howto/comptel/howto_comptel_1.8MeV_tsmap_5vp.png b/doc/source/users/tutorials/howto/comptel/howto_comptel_1.8MeV_tsmap_5vp.png new file mode 100644 index 0000000..6412900 Binary files /dev/null and b/doc/source/users/tutorials/howto/comptel/howto_comptel_1.8MeV_tsmap_5vp.png differ diff --git a/doc/source/users/tutorials/howto/comptel/howto_comptel_dri.rst b/doc/source/users/tutorials/howto/comptel/howto_comptel_dri.rst new file mode 100644 index 0000000..5972d75 --- /dev/null +++ b/doc/source/users/tutorials/howto/comptel/howto_comptel_dri.rst @@ -0,0 +1,300 @@ +.. _howto_comptel_dri: + +Generating your own event selection +----------------------------------- + + .. admonition:: What you will learn + + You will learn how to use the GammaLib classes to generate your own + data cubes and response functions for a COMPTEL analysis. + +The data in the +`HEASARC archive `_ +comprises counts cubes for the standard COMPTEL energy ranges (0.75-1 MeV, +1-3 MeV, 3-10 MeV and 10-30 MeV) per viewing period. In case that you want +to analyse a different energy range, or that you want to perform a different +event selection, you need to create ``DRE``, ``DRG``, ``DRB`` and ``DRX`` +data and the corresponding ``IAQ`` response function yourself. You can do +this with GammaLib. + +Let's assume you want to make a 1.8 MeV line map of the Galactic centre. +You can use for this the data of the COMPTEL viewing period 16 in phase 1 of +the mission. Go to the +`HEASARC archive `_ +and download the ``EVP``, ``TIM`` and ``OAD`` datasets. Once this is done, +create the following +:ref:`observation definition file `: + +.. code-block:: xml + + + + + + + + + + + + + + + + + + + + + + + + +.. + + .. note:: + + There exist two formats for the + :ref:`observation definition file ` + for COMPTEL data, one that describes the locations of the COMPTEL + data space cubes and response functions, and one that desribes + the location of an ``EVP`` event list, a ``TIM`` Good Time + Interval file, and an arbitrary number of ``OAD`` Orbit Aspect + Data files. In general there exists one ``OAD`` per day, while + there is one ``EVP`` event file and one ``TIM`` Good Time Interval + file per viewing period. + + In the example above, the environment variable ``$COMDATA`` was set + to the root directory of the COMPTEL archive. + +Now you are ready to generate the ``DRE``, ``DRG``, ``DRB``, ``DRX`` and +``IAQ`` datasets that you need for a ctools analysis of COMPTEL data. +This will take a few minutes. Start +`Python `_ +in your terminal and type the following: + +.. code-block:: python + + # You need GammaLib, of course ;) + import gammalib + + # Create sky map for the definition of the DRE and DRG data cubes + cube = gammalib.GSkyMap('TAN', 'GAL', 0.0, 20.0, 1.0, 1.0, 80, 80) + + # Create sky map for the definition of the DRX exposure map + expo = gammalib.GSkyMap('CAR', 'GAL', 0.0, 0.0, 1.0, 1.0, 360, 180) + + # Load observation definition XML file 'obs.xml' (the file you created above) + obsdef = gammalib.GObservations('obs.xml') + + # Extract COMPTEL observation from that file (it's the first and unique observation) + obs = obsdef[0] + + # Allocate DRE, DRG and DRX. DRE and DRG have a 3rd Phibar dimension, starting at 0.0 deg, + # with a binsize of 2.0 deg, and 25 Phibar layers. + dre = gammalib.GCOMDri(cube, 0.0, 2.0, 25) + drg = gammalib.GCOMDri(cube, 0.0, 2.0, 25) + drx = gammalib.GCOMDri(expo) + + # Set DRE energy range to 1.7-1.9 MeV, DRG and DRX are energy independent + dre.ebounds(gammalib.GEbounds(gammalib.GEnergy(1.7, 'MeV'), + gammalib.GEnergy(1.9, 'MeV'))) + + # Compute DRE, DRG and DRX (here is where the real magic happens) + dre.compute_dre(obs.events(), obs.oads(), obs.tim()) + drg.compute_drg(obs.oads(), obs.tim()) + drx.compute_drx(obs.oads(), obs.tim()) + + # Save DRE, DRG and DRX ('True' indicates to overwrite any existing file) + dre.save('dre.fits', True) + drg.save('drg.fits', True) + drx.save('drx.fits', True) + + # Load DRE and DRG to generate a DRB background model cube + dre = gammalib.GCOMDri('dre.fits') + drb = gammalib.GCOMDri('drg.fits') + + # Normalise DRB on the Phibar distribution of the DRE cube + npix = dre.nchi() * dre.npsi() + for k in range(dre.nphibar()): + sum_dre = 0.0 + sum_drb = 0.0 + for i in range(npix): + index = i + k * npix + sum_dre += dre[index] + sum_drb += drb[index] + if sum_drb > 0: + for i in range(npix): + index = i + k * npix + drb[index] *= sum_dre / sum_drb + + # Save DRB + drb.save('drb.fits', True) + + # Initialise IAQ that will hold the response function + iaq = gammalib.GCOMIaq(55.0, 1.0, 50.0, 2.0) + + # Compute IAQ for a line energy of 1.809 MeV and an energy band of 1.7-1.9 MeV + iaq.set(gammalib.GEnergy(1.809, 'MeV'), + gammalib.GEbounds(gammalib.GEnergy(1.7, 'MeV'), + gammalib.GEnergy(1.9, 'MeV'))) + + # Save IAQ and you are done + iaq.save('iaq.fits', True) + +Now you have everything at hand to perform a COMPTEL maximum likelihood +analysis. For that purpoe you need to gather all the datasets that you +just created in a new +:ref:`observation definition file ` +that should look as follows: + +.. code-block:: xml + + + + + + + + + + + +.. + + .. warning:: + + Be aware that the attribute for the ``IAQ`` parameter is ``value`` and + not ``file`` since the ``IAQ`` parameter is not necessarily a file + but can be also a response name of the calibration database. + +Before doing a model fit you need a model. Create the following +:ref:`model definition file `: + +.. code-block:: xml + + + + + + + + + + + + + + + + + + + + + + ... + + + + + + + +This file contains a single point source at the position of the Galactic +Centre. The spectral model is a simple constant normalisation that will +return the gamma-ray line flux in units of +:math:`{\rm photons}\,{\rm cm}^{-2}\,{\rm s}^{-1}`. +For the background we do a Phibar fitting of the ``DRB`` cube. Since the +first four layers of the ``DRE`` cube are empty we start the nodes at the +fifth layer which corresponds to a Phibar value of 9 degrees. There are +subsequent nodes spaced by 2 degrees (not shown) up to a Phibar value of +49 degrees. + +Now it's time for model fitting. You can produce for example a Test Statistic +map of the region around the Galactic centre as follows: + +.. code-block:: bash + + $ cttsmap + Input event list, counts cube or observation definition XML file [events.fits] obs_dri.xml + Test source name [Crab] GC + Input model definition XML file [$CTOOLS/share/models/crab.xml] models.xml + First coordinate of image center in degrees (RA or galactic l) (0-360) [83.63] 0.0 + Second coordinate of image center in degrees (DEC or galactic b) (-90-90) [22.01] 0.0 + Projection method (AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN) [CAR] + Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL] GAL + Image scale (in degrees/pixel) [0.02] 1.0 + Size of the X axis in pixels [200] 50 + Size of the Y axis in pixels [200] 30 + Output Test Statistic map file [tsmap.fits] tsmap.fits + +Below is the resulting Test Statistic map that shows 1.8 MeV emission +following the Galactic plane and peaking near the Galactic centre. + +.. figure:: howto_comptel_1.8MeV_tsmap.png + :width: 500px + :align: center + + *1.8 MeV Test Statistic map for COMPTEL viewing period 16 of the Galactic Centre* + +It is left to you as an exercise to extend this example. To get a better +statistics you may for example combine observations. For that purpose you +can add the ``DRI`` datasets of multiple observations into single ``DRE``, +``DRG``, ``DRB`` and ``DRX`` datasets, +that's the way how the COMPTEL analysis was done. You can however also leave +the ``DRI`` datasets separate and combine their description in the +:ref:`observation definition file `. +To combine for example the data of viewing periods 5, 7.5, 13, 16 and 27 +the +:ref:`observation definition file ` +should look like this: + +.. code-block:: xml + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +The corresponding Test Statistic map is shown below. + +.. figure:: howto_comptel_1.8MeV_tsmap_5vp.png + :width: 500px + :align: center + + *1.8 MeV Test Statistic map for COMPTEL viewing periods 5, 7.5, 13, 16 and 27 of the Galactic Centre* diff --git a/doc/source/users/tutorials/howto/comptel/index.rst b/doc/source/users/tutorials/howto/comptel/index.rst index 74d5820..c9400d9 100644 --- a/doc/source/users/tutorials/howto/comptel/index.rst +++ b/doc/source/users/tutorials/howto/comptel/index.rst @@ -22,3 +22,4 @@ In this tutorial you will learn how to do that with a simple analysis. howto_comptel_spectrum howto_comptel_tsmap howto_comptel_ulimit + howto_comptel_dri -- libgit2 0.21.4