diff --git a/ChangeLog b/ChangeLog index 65bf1cd..aa316c2 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,9 +1,10 @@ -2017-08-08 +2017-10-10 * Version 1.5.0 released ======================== Use t-type parameters for all "tmin" and "tmax" parameters (#1864) + Add keywords to sky maps produced by ctskymap (#2194) 2017-07-28 diff --git a/NEWS b/NEWS index 2940dcf..9422c63 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ New Features and Important Changes in ctools 1.5.0 -8 August 2017 +10 October 2017 Introduction @@ -53,7 +53,9 @@ None ctskymap - CTA sky mapping -------------------------- -None +Header keywords informing about the dates and times covered by the observations +used to generate the skymap, as well as the event energy range, are now appended +as header keywords to the FITS file (#2194). ctmodel - CTA counts cube modelling diff --git a/doc/source/admin/release_history/1.5.rst b/doc/source/admin/release_history/1.5.rst index 9e4ab6b..1f96d96 100644 --- a/doc/source/admin/release_history/1.5.rst +++ b/doc/source/admin/release_history/1.5.rst @@ -22,7 +22,8 @@ Bug fixes Improvements ------------ -* TBW +* [`2194 `_] - + Add keywords to sky maps produced by :ref:`ctskymap` New features 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 --- /dev/null +++ b/doc/source/users/tutorials/howto/comptel/howto_comptel_1.8MeV_tsmap.png 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 --- /dev/null +++ b/doc/source/users/tutorials/howto/comptel/howto_comptel_1.8MeV_tsmap_5vp.png 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..75c8e81 --- /dev/null +++ b/doc/source/users/tutorials/howto/comptel/howto_comptel_dri.rst @@ -0,0 +1,297 @@ +.. _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) + + # 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'))) + + # Load observation definition XML file 'obs.xml' (the file you created above) + obsdef = gammalib.GObservations('obs.xml') + + # Compute DRE, DRG and DRX for the first observation in the XML file + dre.compute_dre(obsdef[0]) + drg.compute_drg(obsdef[0]) + drx.compute_drx(obsdef[0]) + + # 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 diff --git a/pyext/ctobservation.i b/pyext/ctobservation.i index 5791958..67e9966 100644 --- a/pyext/ctobservation.i +++ b/pyext/ctobservation.i @@ -1,7 +1,7 @@ /*************************************************************************** * ctobservation - Base class for observation tools * * ----------------------------------------------------------------------- * - * copyright (C) 2016 by Juergen Knoedlseder * + * copyright (C) 2016-2017 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -89,6 +89,7 @@ public: // Public methods const GObservations& obs(void) const; + void write_ogip_keywords(GFitsHDU* hdu) const; // Make methods private in Python by prepending an underscore %rename(_first_unbinned_observation) first_unbinned_observation; diff --git a/src/ctskymap/ctskymap.cpp b/src/ctskymap/ctskymap.cpp index 33ea26f..6ecaa1a 100644 --- a/src/ctskymap/ctskymap.cpp +++ b/src/ctskymap/ctskymap.cpp @@ -283,23 +283,35 @@ void ctskymap::save(void) GFits fits; // Write sky map into FITS file - m_skymap.write(fits); + GFitsHDU* hdu = m_skymap.write(fits); + + // Write keywords into sky map extension + write_ogip_keywords(hdu); + write_hdu_keywords(hdu); // If background subtraction is requested then write background map // and significance map to FITS file if (m_bkgsubtract != "NONE") { // Write background map into FITS file - m_bkgmap.write(fits); + hdu = m_bkgmap.write(fits); // Set background map extension name - fits[fits.size()-1]->extname("BACKGROUND"); + hdu->extname("BACKGROUND"); + + // Write keywords into background extension + write_ogip_keywords(hdu); + write_hdu_keywords(hdu); // Write significance map into FITS file - m_sigmap.write(fits); + hdu = m_sigmap.write(fits); // Set significance map extension name - fits[fits.size()-1]->extname("SIGNIFICANCE"); + hdu->extname("SIGNIFICANCE"); + + // Write keywords into significance extension + write_ogip_keywords(hdu); + write_hdu_keywords(hdu); } // endif: background subtraction was requested @@ -673,6 +685,31 @@ void ctskymap::map_background_irf(GCTAObservation* obs) /***********************************************************************//** + * @brief Write keywords in FITS HDU + * + * @param[in,out] hdu Pointer to FITS HDU. + * + * Writes keywords in FITS HDU. + ***************************************************************************/ +void ctskymap::write_hdu_keywords(GFitsHDU* hdu) const +{ + // Continue only if pointer is valid + if (hdu != NULL) { + + // Set keywords + hdu->card("BKGSUB", m_bkgsubtract, "Background substraction method"); + hdu->card("E_MIN", m_emin, "[TeV] Lower energy boundary"); + hdu->card("E_MAX", m_emax, "[TeV] Upper energy boundary"); + hdu->card("EUNIT", "TeV", "Units for E_MIN and E_MAX"); + + } // endif: pointer was valid + + // Return + return; +} + + +/***********************************************************************//** * @brief Estimates the background in sky map based on IRF template * * @param[in] lnE Natural logarithm of energy in MeV. diff --git a/src/ctskymap/ctskymap.hpp b/src/ctskymap/ctskymap.hpp index cf16a85..da62e77 100644 --- a/src/ctskymap/ctskymap.hpp +++ b/src/ctskymap/ctskymap.hpp @@ -71,6 +71,7 @@ protected: void map_events(GCTAObservation* obs); void map_background(GCTAObservation* obs); void map_background_irf(GCTAObservation* obs); + void write_hdu_keywords(GFitsHDU* hdu) const; // Background integration kernel class irf_kern : public GFunction { diff --git a/src/support/ctobservation.cpp b/src/support/ctobservation.cpp index 4006a9b..2e03c6d 100644 --- a/src/support/ctobservation.cpp +++ b/src/support/ctobservation.cpp @@ -192,6 +192,53 @@ ctobservation& ctobservation::operator=(const ctobservation& app) = = ==========================================================================*/ +/***********************************************************************//** + * @brief Write OGIP keywords in FITS HDU + * + * @param[in,out] hdu Pointer to FITS HDU. + * + * Writes OGIP keywords in FITS HDU. + ***************************************************************************/ +void ctobservation::write_ogip_keywords(GFitsHDU* hdu) const +{ + // Continue only if pointer is valid + if (hdu != NULL) { + + // Set creator + std::string creator = this->name() + " v" + this->version(); + + // Compute times + std::string utc_obs = m_ogip_tstart.utc(); + std::string utc_end = m_ogip_tstop.utc(); + std::string date_obs = utc_obs.substr(0, 10); + std::string time_obs = utc_obs.substr(11, 8); + std::string date_end = utc_end.substr(0, 10); + std::string time_end = utc_end.substr(11, 8); + + // Compute deadtime correction + double deadc = (m_ogip_ontime > 0.0) ? m_ogip_livetime / m_ogip_ontime : 1.0; + + // Set observation information + hdu->card("CREATOR", creator, "Program which created the file"); + hdu->card("TELESCOP", m_ogip_telescope, "Telescope"); + + // Set observation time information + hdu->card("DATE_OBS", date_obs, "Observation start date"); + hdu->card("TIME_OBS", time_obs, "Observation start time"); + hdu->card("DATE_END", date_end, "Observation end date"); + hdu->card("TIME_END", time_end, "Observation end time"); + hdu->card("TELAPSE", m_ogip_telapse, "[s] Elapsed time"); + hdu->card("ONTIME", m_ogip_ontime, "[s] Total good time including deadtime"); + hdu->card("LIVETIME", m_ogip_livetime, "[s] Total livetime"); + hdu->card("EXPOSURE", m_ogip_exposure, "[s] Exposure time"); + hdu->card("DEADC", deadc, "Deadtime correction factor"); + + } // endif: pointer was valid + + // Return + return; +} + /*========================================================================== = = @@ -208,6 +255,13 @@ void ctobservation::init_members(void) m_obs.clear(); // Initialise private members + m_ogip_telescope.clear(); + m_ogip_tstart.clear(); + m_ogip_tstop.clear(); + m_ogip_telapse = 0.0; + m_ogip_exposure = 0.0; + m_ogip_ontime = 0.0; + m_ogip_livetime = 0.0; m_index_unbinned = 0; // Return @@ -226,6 +280,13 @@ void ctobservation::copy_members(const ctobservation& app) m_obs = app.m_obs; // Copy private members + m_ogip_telescope = app.m_ogip_telescope; + m_ogip_tstart = app.m_ogip_tstart; + m_ogip_tstop = app.m_ogip_tstop; + m_ogip_telapse = app.m_ogip_telapse; + m_ogip_exposure = app.m_ogip_exposure; + m_ogip_ontime = app.m_ogip_ontime; + m_ogip_livetime = app.m_ogip_livetime; m_index_unbinned = app.m_index_unbinned; // Return @@ -258,6 +319,15 @@ const GCTAObservation* ctobservation::first_unbinned_observation(void) const // Initialise index m_index_unbinned = 0; + // Initialise OGIP members + m_ogip_telescope.clear(); + m_ogip_tstart.clear(); + m_ogip_tstop.clear(); + m_ogip_telapse = 0.0; + m_ogip_exposure = 0.0; + m_ogip_ontime = 0.0; + m_ogip_livetime = 0.0; + // Get next unbinned CTA observation const GCTAObservation* obs = next_unbinned_observation(); @@ -323,6 +393,31 @@ const GCTAObservation* ctobservation::next_unbinned_observation(void) const } // endfor: looped over all observation + // Update OGIP members + if (obs != NULL) { + if (m_ogip_telapse == 0.0) { + m_ogip_telescope = obs->instrument(); + m_ogip_tstart = obs->gti().tstart(); + m_ogip_tstop = obs->gti().tstop(); + m_ogip_telapse = obs->gti().telapse(); + m_ogip_exposure = obs->ontime(); + m_ogip_ontime = obs->ontime(); + m_ogip_livetime = obs->livetime(); + } + else { + if (obs->gti().tstart() < m_ogip_tstart) { + m_ogip_tstart = obs->gti().tstart(); + } + if (obs->gti().tstop() > m_ogip_tstop) { + m_ogip_tstop = obs->gti().tstop(); + } + m_ogip_telapse += obs->gti().telapse(); + m_ogip_exposure += obs->ontime(); + m_ogip_ontime += obs->ontime(); + m_ogip_livetime += obs->livetime(); + } + } + // Return next CTA observation return obs; } diff --git a/src/support/ctobservation.hpp b/src/support/ctobservation.hpp index 3e6dc7a..7b0bee0 100644 --- a/src/support/ctobservation.hpp +++ b/src/support/ctobservation.hpp @@ -63,6 +63,7 @@ public: // Methods const GObservations& obs(void) const; + void write_ogip_keywords(GFitsHDU* hdu) const; #ifndef SWIG protected: @@ -82,11 +83,18 @@ protected: void save_events_xml(void); // Protected members - GObservations m_obs; //!< Observation container + GObservations m_obs; //!< Observation container private: // Private members - mutable int m_index_unbinned; //!< Current index of unbinned observation + mutable std::string m_ogip_telescope; //!< Name of telescope + mutable GTime m_ogip_tstart; //!< Start time for OGIP keywords + mutable GTime m_ogip_tstop; //!< Stop time for OGIP keywords + mutable double m_ogip_telapse; //!< Elapsed time + mutable double m_ogip_exposure; //!< Exposure time + mutable double m_ogip_ontime; //!< Ontime for OGIP keywords + mutable double m_ogip_livetime; //!< Livetime for OGIP keywords + mutable int m_index_unbinned; //!< Current index of unbinned observation };