Nathan Kelley-Hoskins / ctools

Commit 3b6811f5d8019592324d3d2092eee4671bce5cd0

Authored by Jürgen Knödlseder
2 parents 5e03ca0d 9955078b

Merge branch '2194-add-exposure-keywords' into integration

Inline Side-by-side
Showing 12 changed files with 456 additions and 12 deletions
ChangeLog
  View file @3b6811f
1   -2017-08-08
  1 +2017-10-10
2 2  
3 3 * Version 1.5.0 released
4 4 ========================
5 5  
6 6 Use t-type parameters for all "tmin" and "tmax" parameters (#1864)
  7 + Add keywords to sky maps produced by ctskymap (#2194)
7 8  
8 9  
9 10 2017-07-28
... ...
1 1 New Features and Important Changes in ctools 1.5.0
2 2  
3   -8 August 2017
  3 +10 October 2017
4 4  
5 5  
6 6 Introduction
... ... @@ -53,7 +53,9 @@ None
53 53  
54 54 ctskymap - CTA sky mapping
55 55 --------------------------
56   -None
  56 +Header keywords informing about the dates and times covered by the observations
  57 +used to generate the skymap, as well as the event energy range, are now appended
  58 +as header keywords to the FITS file (#2194).
57 59  
58 60  
59 61 ctmodel - CTA counts cube modelling
... ...
doc/source/admin/release_history/1.5.rst
  View file @3b6811f
... ... @@ -22,7 +22,8 @@ Bug fixes
22 22 Improvements
23 23 ------------
24 24  
25   -* TBW
  25 +* [`2194 <https://cta-redmine.irap.omp.eu/issues/2194>`_] -
  26 + Add keywords to sky maps produced by :ref:`ctskymap`
26 27  
27 28  
28 29 New features
... ...
doc/source/users/tutorials/howto/comptel/howto_comptel_1.8MeV_tsmap.png 0 → 100644
View file @3b6811f

87.6 KB
doc/source/users/tutorials/howto/comptel/howto_comptel_1.8MeV_tsmap_5vp.png 0 → 100644
View file @3b6811f

59.1 KB
doc/source/users/tutorials/howto/comptel/howto_comptel_dri.rst 0 → 100644
  View file @3b6811f
  1 +.. _howto_comptel_dri:
  2 +
  3 +Generating your own event selection
  4 +-----------------------------------
  5 +
  6 + .. admonition:: What you will learn
  7 +
  8 + You will learn how to use the GammaLib classes to generate your own
  9 + data cubes and response functions for a COMPTEL analysis.
  10 +
  11 +The data in the
  12 +`HEASARC archive <https://heasarc.gsfc.nasa.gov/xamin/>`_
  13 +comprises counts cubes for the standard COMPTEL energy ranges (0.75-1 MeV,
  14 +1-3 MeV, 3-10 MeV and 10-30 MeV) per viewing period. In case that you want
  15 +to analyse a different energy range, or that you want to perform a different
  16 +event selection, you need to create ``DRE``, ``DRG``, ``DRB`` and ``DRX``
  17 +data and the corresponding ``IAQ`` response function yourself. You can do
  18 +this with GammaLib.
  19 +
  20 +Let's assume you want to make a 1.8 MeV line map of the Galactic centre.
  21 +You can use for this the data of the COMPTEL viewing period 16 in phase 1 of
  22 +the mission. Go to the
  23 +`HEASARC archive <https://heasarc.gsfc.nasa.gov/xamin/>`_
  24 +and download the ``EVP``, ``TIM`` and ``OAD`` datasets. Once this is done,
  25 +create the following
  26 +:ref:`observation definition file <glossary_obsdef>`:
  27 +
  28 +.. code-block:: xml
  29 +
  30 + <?xml version="1.0" standalone="no"?>
  31 + <observation_list title="observation library">
  32 + <observation name="GC" id="000160" instrument="COM">
  33 + <parameter name="EVP" file="$COMDATA/phase01/vp0016_0/m28846_evp.fits"/>
  34 + <parameter name="TIM" file="$COMDATA/phase01/vp0016_0/m11765_tim.fits"/>
  35 + <parameter name="OAD" file="$COMDATA/phase01/vp0016_0/m20745_oad.fits"/>
  36 + <parameter name="OAD" file="$COMDATA/phase01/vp0016_0/m20750_oad.fits"/>
  37 + <parameter name="OAD" file="$COMDATA/phase01/vp0016_0/m20753_oad.fits"/>
  38 + <parameter name="OAD" file="$COMDATA/phase01/vp0016_0/m20759_oad.fits"/>
  39 + <parameter name="OAD" file="$COMDATA/phase01/vp0016_0/m20760_oad.fits"/>
  40 + <parameter name="OAD" file="$COMDATA/phase01/vp0016_0/m20761_oad.fits"/>
  41 + <parameter name="OAD" file="$COMDATA/phase01/vp0016_0/m20762_oad.fits"/>
  42 + <parameter name="OAD" file="$COMDATA/phase01/vp0016_0/m20764_oad.fits"/>
  43 + <parameter name="OAD" file="$COMDATA/phase01/vp0016_0/m20770_oad.fits"/>
  44 + <parameter name="OAD" file="$COMDATA/phase01/vp0016_0/m20771_oad.fits"/>
  45 + <parameter name="OAD" file="$COMDATA/phase01/vp0016_0/m20772_oad.fits"/>
  46 + <parameter name="OAD" file="$COMDATA/phase01/vp0016_0/m20785_oad.fits"/>
  47 + <parameter name="OAD" file="$COMDATA/phase01/vp0016_0/m20786_oad.fits"/>
  48 + <parameter name="OAD" file="$COMDATA/phase01/vp0016_0/m20787_oad.fits"/>
  49 + <parameter name="OAD" file="$COMDATA/phase01/vp0016_0/m20788_oad.fits"/>
  50 + <parameter name="OAD" file="$COMDATA/phase01/vp0016_0/m20789_oad.fits"/>
  51 + </observation>
  52 + </observation_list>
  53 +..
  54 +
  55 + .. note::
  56 +
  57 + There exist two formats for the
  58 + :ref:`observation definition file <glossary_obsdef>`
  59 + for COMPTEL data, one that describes the locations of the COMPTEL
  60 + data space cubes and response functions, and one that desribes
  61 + the location of an ``EVP`` event list, a ``TIM`` Good Time
  62 + Interval file, and an arbitrary number of ``OAD`` Orbit Aspect
  63 + Data files. In general there exists one ``OAD`` per day, while
  64 + there is one ``EVP`` event file and one ``TIM`` Good Time Interval
  65 + file per viewing period.
  66 +
  67 + In the example above, the environment variable ``$COMDATA`` was set
  68 + to the root directory of the COMPTEL archive.
  69 +
  70 +Now you are ready to generate the ``DRE``, ``DRG``, ``DRB``, ``DRX`` and
  71 +``IAQ`` datasets that you need for a ctools analysis of COMPTEL data.
  72 +This will take a few minutes. Start
  73 +`Python <https://www.python.org>`_
  74 +in your terminal and type the following:
  75 +
  76 +.. code-block:: python
  77 +
  78 + # You need GammaLib, of course ;)
  79 + import gammalib
  80 +
  81 + # Create sky map for the definition of the DRE and DRG data cubes
  82 + cube = gammalib.GSkyMap('TAN', 'GAL', 0.0, 20.0, 1.0, 1.0, 80, 80)
  83 +
  84 + # Create sky map for the definition of the DRX exposure map
  85 + expo = gammalib.GSkyMap('CAR', 'GAL', 0.0, 0.0, 1.0, 1.0, 360, 180)
  86 +
  87 + # Allocate DRE, DRG and DRX. DRE and DRG have a 3rd Phibar dimension, starting at 0.0 deg,
  88 + # with a binsize of 2.0 deg, and 25 Phibar layers.
  89 + dre = gammalib.GCOMDri(cube, 0.0, 2.0, 25)
  90 + drg = gammalib.GCOMDri(cube, 0.0, 2.0, 25)
  91 + drx = gammalib.GCOMDri(expo)
  92 +
  93 + # Set DRE energy range to 1.7-1.9 MeV, DRG and DRX are energy independent
  94 + dre.ebounds(gammalib.GEbounds(gammalib.GEnergy(1.7, 'MeV'),
  95 + gammalib.GEnergy(1.9, 'MeV')))
  96 +
  97 + # Load observation definition XML file 'obs.xml' (the file you created above)
  98 + obsdef = gammalib.GObservations('obs.xml')
  99 +
  100 + # Compute DRE, DRG and DRX for the first observation in the XML file
  101 + dre.compute_dre(obsdef[0])
  102 + drg.compute_drg(obsdef[0])
  103 + drx.compute_drx(obsdef[0])
  104 +
  105 + # Save DRE, DRG and DRX ('True' indicates to overwrite any existing file)
  106 + dre.save('dre.fits', True)
  107 + drg.save('drg.fits', True)
  108 + drx.save('drx.fits', True)
  109 +
  110 + # Load DRE and DRG to generate a DRB background model cube
  111 + dre = gammalib.GCOMDri('dre.fits')
  112 + drb = gammalib.GCOMDri('drg.fits')
  113 +
  114 + # Normalise DRB on the Phibar distribution of the DRE cube
  115 + npix = dre.nchi() * dre.npsi()
  116 + for k in range(dre.nphibar()):
  117 + sum_dre = 0.0
  118 + sum_drb = 0.0
  119 + for i in range(npix):
  120 + index = i + k * npix
  121 + sum_dre += dre[index]
  122 + sum_drb += drb[index]
  123 + if sum_drb > 0:
  124 + for i in range(npix):
  125 + index = i + k * npix
  126 + drb[index] *= sum_dre / sum_drb
  127 +
  128 + # Save DRB
  129 + drb.save('drb.fits', True)
  130 +
  131 + # Initialise IAQ that will hold the response function
  132 + iaq = gammalib.GCOMIaq(55.0, 1.0, 50.0, 2.0)
  133 +
  134 + # Compute IAQ for a line energy of 1.809 MeV and an energy band of 1.7-1.9 MeV
  135 + iaq.set(gammalib.GEnergy(1.809, 'MeV'),
  136 + gammalib.GEbounds(gammalib.GEnergy(1.7, 'MeV'),
  137 + gammalib.GEnergy(1.9, 'MeV')))
  138 +
  139 + # Save IAQ and you are done
  140 + iaq.save('iaq.fits', True)
  141 +
  142 +Now you have everything at hand to perform a COMPTEL maximum likelihood
  143 +analysis. For that purpoe you need to gather all the datasets that you
  144 +just created in a new
  145 +:ref:`observation definition file <glossary_obsdef>`
  146 +that should look as follows:
  147 +
  148 +.. code-block:: xml
  149 +
  150 + <?xml version="1.0" standalone="no"?>
  151 + <observation_list title="observation library">
  152 + <observation name="GC" id="000160" instrument="COM">
  153 + <parameter name="DRE" file="dre.fits"/>
  154 + <parameter name="DRB" file="drb.fits"/>
  155 + <parameter name="DRG" file="drg.fits"/>
  156 + <parameter name="DRX" file="drx.fits"/>
  157 + <parameter name="IAQ" value="iaq.fits"/>
  158 + </observation>
  159 + </observation_list>
  160 +..
  161 +
  162 + .. warning::
  163 +
  164 + Be aware that the attribute for the ``IAQ`` parameter is ``value`` and
  165 + not ``file`` since the ``IAQ`` parameter is not necessarily a file
  166 + but can be also a response name of the calibration database.
  167 +
  168 +Before doing a model fit you need a model. Create the following
  169 +:ref:`model definition file <glossary_moddef>`:
  170 +
  171 +.. code-block:: xml
  172 +
  173 + <?xml version="1.0" encoding="UTF-8" standalone="no"?>
  174 + <source_library title="source library">
  175 + <source name="GC" type="PointSource" tscalc="1">
  176 + <spectrum type="Constant">
  177 + <parameter name="Normalization" scale="1.0e-5" value="1.0" min="-100.0" max="100.0" free="1"/>
  178 + </spectrum>
  179 + <spatialModel type="PointSource">
  180 + <parameter name="RA" scale="1.0" value="266.40" min="-360" max="360" free="0"/>
  181 + <parameter name="DEC" scale="1.0" value="-28.94" min="-90" max="90" free="0"/>
  182 + </spatialModel>
  183 + </source>
  184 + <source name="Background" type="DRBFitting" instrument="COM">
  185 + <node>
  186 + <parameter name="Phibar" value="9" scale="1" min="0" max="50" free="0" />
  187 + <parameter name="Normalization" value="1.0" scale="1" min="0" max="1000" free="1" />
  188 + </node>
  189 + <node>
  190 + <parameter name="Phibar" value="11" scale="1" min="0" max="50" free="0" />
  191 + <parameter name="Normalization" value="1.0" scale="1" min="0" max="1000" free="1" />
  192 + </node>
  193 + ...
  194 + <node>
  195 + <parameter name="Phibar" value="49" scale="1" min="0" max="50" free="0" />
  196 + <parameter name="Normalization" value="1.0" scale="1" min="0" max="1000" free="1" />
  197 + </node>
  198 + </source>
  199 + </source_library>
  200 +
  201 +This file contains a single point source at the position of the Galactic
  202 +Centre. The spectral model is a simple constant normalisation that will
  203 +return the gamma-ray line flux in units of
  204 +:math:`{\rm photons}\,{\rm cm}^{-2}\,{\rm s}^{-1}`.
  205 +For the background we do a Phibar fitting of the ``DRB`` cube. Since the
  206 +first four layers of the ``DRE`` cube are empty we start the nodes at the
  207 +fifth layer which corresponds to a Phibar value of 9 degrees. There are
  208 +subsequent nodes spaced by 2 degrees (not shown) up to a Phibar value of
  209 +49 degrees.
  210 +
  211 +Now it's time for model fitting. You can produce for example a Test Statistic
  212 +map of the region around the Galactic centre as follows:
  213 +
  214 +.. code-block:: bash
  215 +
  216 + $ cttsmap
  217 + Input event list, counts cube or observation definition XML file [events.fits] obs_dri.xml
  218 + Test source name [Crab] GC
  219 + Input model definition XML file [$CTOOLS/share/models/crab.xml] models.xml
  220 + First coordinate of image center in degrees (RA or galactic l) (0-360) [83.63] 0.0
  221 + Second coordinate of image center in degrees (DEC or galactic b) (-90-90) [22.01] 0.0
  222 + Projection method (AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN) [CAR]
  223 + Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL] GAL
  224 + Image scale (in degrees/pixel) [0.02] 1.0
  225 + Size of the X axis in pixels [200] 50
  226 + Size of the Y axis in pixels [200] 30
  227 + Output Test Statistic map file [tsmap.fits] tsmap.fits
  228 +
  229 +Below is the resulting Test Statistic map that shows 1.8 MeV emission
  230 +following the Galactic plane and peaking near the Galactic centre.
  231 +
  232 +.. figure:: howto_comptel_1.8MeV_tsmap.png
  233 + :width: 500px
  234 + :align: center
  235 +
  236 + *1.8 MeV Test Statistic map for COMPTEL viewing period 16 of the Galactic Centre*
  237 +
  238 +It is left to you as an exercise to extend this example. To get a better
  239 +statistics you may for example combine observations. For that purpose you
  240 +can add the ``DRI`` datasets of multiple observations into single ``DRE``,
  241 +``DRG``, ``DRB`` and ``DRX`` datasets,
  242 +that's the way how the COMPTEL analysis was done. You can however also leave
  243 +the ``DRI`` datasets separate and combine their description in the
  244 +:ref:`observation definition file <glossary_obsdef>`.
  245 +To combine for example the data of viewing periods 5, 7.5, 13, 16 and 27
  246 +the
  247 +:ref:`observation definition file <glossary_obsdef>`
  248 +should look like this:
  249 +
  250 +.. code-block:: xml
  251 +
  252 + <?xml version="1.0" standalone="no"?>
  253 + <observation_list title="observation library">
  254 + <observation name="GC" id="000050" instrument="COM">
  255 + <parameter name="DRE" file="dre_0050.fits"/>
  256 + <parameter name="DRB" file="drb_0050.fits"/>
  257 + <parameter name="DRG" file="drg_0050.fits"/>
  258 + <parameter name="DRX" file="drx_0050.fits"/>
  259 + <parameter name="IAQ" value="iaq.fits"/>
  260 + </observation>
  261 + <observation name="GC" id="000075" instrument="COM">
  262 + <parameter name="DRE" file="dre_0075.fits"/>
  263 + <parameter name="DRB" file="drb_0075.fits"/>
  264 + <parameter name="DRG" file="drg_0075.fits"/>
  265 + <parameter name="DRX" file="drx_0075.fits"/>
  266 + <parameter name="IAQ" value="iaq.fits"/>
  267 + </observation>
  268 + <observation name="GC" id="000130" instrument="COM">
  269 + <parameter name="DRE" file="dre_0130.fits"/>
  270 + <parameter name="DRB" file="drb_0130.fits"/>
  271 + <parameter name="DRG" file="drg_0130.fits"/>
  272 + <parameter name="DRX" file="drx_0130.fits"/>
  273 + <parameter name="IAQ" value="iaq.fits"/>
  274 + </observation>
  275 + <observation name="GC" id="000160" instrument="COM">
  276 + <parameter name="DRE" file="dre_0160.fits"/>
  277 + <parameter name="DRB" file="drb_0160.fits"/>
  278 + <parameter name="DRG" file="drg_0160.fits"/>
  279 + <parameter name="DRX" file="drx_0160.fits"/>
  280 + <parameter name="IAQ" value="iaq.fits"/>
  281 + </observation>
  282 + <observation name="GC" id="000270" instrument="COM">
  283 + <parameter name="DRE" file="dre_0270.fits"/>
  284 + <parameter name="DRB" file="drb_0270.fits"/>
  285 + <parameter name="DRG" file="drg_0270.fits"/>
  286 + <parameter name="DRX" file="drx_0270.fits"/>
  287 + <parameter name="IAQ" value="iaq.fits"/>
  288 + </observation>
  289 + </observation_list>
  290 +
  291 +The corresponding Test Statistic map is shown below.
  292 +
  293 +.. figure:: howto_comptel_1.8MeV_tsmap_5vp.png
  294 + :width: 500px
  295 + :align: center
  296 +
  297 + *1.8 MeV Test Statistic map for COMPTEL viewing periods 5, 7.5, 13, 16 and 27 of the Galactic Centre*
... ...
doc/source/users/tutorials/howto/comptel/index.rst
  View file @3b6811f
... ... @@ -22,3 +22,4 @@ In this tutorial you will learn how to do that with a simple analysis.
22 22 howto_comptel_spectrum
23 23 howto_comptel_tsmap
24 24 howto_comptel_ulimit
  25 + howto_comptel_dri
... ...
pyext/ctobservation.i
  View file @3b6811f
1 1 /***************************************************************************
2 2 * ctobservation - Base class for observation tools *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2016-2017 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -89,6 +89,7 @@ public:
89 89  
90 90 // Public methods
91 91 const GObservations& obs(void) const;
  92 + void write_ogip_keywords(GFitsHDU* hdu) const;
92 93  
93 94 // Make methods private in Python by prepending an underscore
94 95 %rename(_first_unbinned_observation) first_unbinned_observation;
... ...
src/ctskymap/ctskymap.cpp
  View file @3b6811f
... ... @@ -283,23 +283,35 @@ void ctskymap::save(void)
283 283 GFits fits;
284 284  
285 285 // Write sky map into FITS file
286   - m_skymap.write(fits);
  286 + GFitsHDU* hdu = m_skymap.write(fits);
  287 +
  288 + // Write keywords into sky map extension
  289 + write_ogip_keywords(hdu);
  290 + write_hdu_keywords(hdu);
287 291  
288 292 // If background subtraction is requested then write background map
289 293 // and significance map to FITS file
290 294 if (m_bkgsubtract != "NONE") {
291 295  
292 296 // Write background map into FITS file
293   - m_bkgmap.write(fits);
  297 + hdu = m_bkgmap.write(fits);
294 298  
295 299 // Set background map extension name
296   - fits[fits.size()-1]->extname("BACKGROUND");
  300 + hdu->extname("BACKGROUND");
  301 +
  302 + // Write keywords into background extension
  303 + write_ogip_keywords(hdu);
  304 + write_hdu_keywords(hdu);
297 305  
298 306 // Write significance map into FITS file
299   - m_sigmap.write(fits);
  307 + hdu = m_sigmap.write(fits);
300 308  
301 309 // Set significance map extension name
302   - fits[fits.size()-1]->extname("SIGNIFICANCE");
  310 + hdu->extname("SIGNIFICANCE");
  311 +
  312 + // Write keywords into significance extension
  313 + write_ogip_keywords(hdu);
  314 + write_hdu_keywords(hdu);
303 315  
304 316 } // endif: background subtraction was requested
305 317  
... ... @@ -673,6 +685,31 @@ void ctskymap::map_background_irf(GCTAObservation* obs)
673 685  
674 686  
675 687 /***********************************************************************//**
  688 + * @brief Write keywords in FITS HDU
  689 + *
  690 + * @param[in,out] hdu Pointer to FITS HDU.
  691 + *
  692 + * Writes keywords in FITS HDU.
  693 + ***************************************************************************/
  694 +void ctskymap::write_hdu_keywords(GFitsHDU* hdu) const
  695 +{
  696 + // Continue only if pointer is valid
  697 + if (hdu != NULL) {
  698 +
  699 + // Set keywords
  700 + hdu->card("BKGSUB", m_bkgsubtract, "Background substraction method");
  701 + hdu->card("E_MIN", m_emin, "[TeV] Lower energy boundary");
  702 + hdu->card("E_MAX", m_emax, "[TeV] Upper energy boundary");
  703 + hdu->card("EUNIT", "TeV", "Units for E_MIN and E_MAX");
  704 +
  705 + } // endif: pointer was valid
  706 +
  707 + // Return
  708 + return;
  709 +}
  710 +
  711 +
  712 +/***********************************************************************//**
676 713 * @brief Estimates the background in sky map based on IRF template
677 714 *
678 715 * @param[in] lnE Natural logarithm of energy in MeV.
... ...
src/ctskymap/ctskymap.hpp
  View file @3b6811f
... ... @@ -71,6 +71,7 @@ protected:
71 71 void map_events(GCTAObservation* obs);
72 72 void map_background(GCTAObservation* obs);
73 73 void map_background_irf(GCTAObservation* obs);
  74 + void write_hdu_keywords(GFitsHDU* hdu) const;
74 75  
75 76 // Background integration kernel
76 77 class irf_kern : public GFunction {
... ...
src/support/ctobservation.cpp
  View file @3b6811f
... ... @@ -192,6 +192,53 @@ ctobservation&amp; ctobservation::operator=(const ctobservation&amp; app)
192 192 = =
193 193 ==========================================================================*/
194 194  
  195 +/***********************************************************************//**
  196 + * @brief Write OGIP keywords in FITS HDU
  197 + *
  198 + * @param[in,out] hdu Pointer to FITS HDU.
  199 + *
  200 + * Writes OGIP keywords in FITS HDU.
  201 + ***************************************************************************/
  202 +void ctobservation::write_ogip_keywords(GFitsHDU* hdu) const
  203 +{
  204 + // Continue only if pointer is valid
  205 + if (hdu != NULL) {
  206 +
  207 + // Set creator
  208 + std::string creator = this->name() + " v" + this->version();
  209 +
  210 + // Compute times
  211 + std::string utc_obs = m_ogip_tstart.utc();
  212 + std::string utc_end = m_ogip_tstop.utc();
  213 + std::string date_obs = utc_obs.substr(0, 10);
  214 + std::string time_obs = utc_obs.substr(11, 8);
  215 + std::string date_end = utc_end.substr(0, 10);
  216 + std::string time_end = utc_end.substr(11, 8);
  217 +
  218 + // Compute deadtime correction
  219 + double deadc = (m_ogip_ontime > 0.0) ? m_ogip_livetime / m_ogip_ontime : 1.0;
  220 +
  221 + // Set observation information
  222 + hdu->card("CREATOR", creator, "Program which created the file");
  223 + hdu->card("TELESCOP", m_ogip_telescope, "Telescope");
  224 +
  225 + // Set observation time information
  226 + hdu->card("DATE_OBS", date_obs, "Observation start date");
  227 + hdu->card("TIME_OBS", time_obs, "Observation start time");
  228 + hdu->card("DATE_END", date_end, "Observation end date");
  229 + hdu->card("TIME_END", time_end, "Observation end time");
  230 + hdu->card("TELAPSE", m_ogip_telapse, "[s] Elapsed time");
  231 + hdu->card("ONTIME", m_ogip_ontime, "[s] Total good time including deadtime");
  232 + hdu->card("LIVETIME", m_ogip_livetime, "[s] Total livetime");
  233 + hdu->card("EXPOSURE", m_ogip_exposure, "[s] Exposure time");
  234 + hdu->card("DEADC", deadc, "Deadtime correction factor");
  235 +
  236 + } // endif: pointer was valid
  237 +
  238 + // Return
  239 + return;
  240 +}
  241 +
195 242  
196 243 /*==========================================================================
197 244 = =
... ... @@ -208,6 +255,13 @@ void ctobservation::init_members(void)
208 255 m_obs.clear();
209 256  
210 257 // Initialise private members
  258 + m_ogip_telescope.clear();
  259 + m_ogip_tstart.clear();
  260 + m_ogip_tstop.clear();
  261 + m_ogip_telapse = 0.0;
  262 + m_ogip_exposure = 0.0;
  263 + m_ogip_ontime = 0.0;
  264 + m_ogip_livetime = 0.0;
211 265 m_index_unbinned = 0;
212 266  
213 267 // Return
... ... @@ -226,6 +280,13 @@ void ctobservation::copy_members(const ctobservation&amp; app)
226 280 m_obs = app.m_obs;
227 281  
228 282 // Copy private members
  283 + m_ogip_telescope = app.m_ogip_telescope;
  284 + m_ogip_tstart = app.m_ogip_tstart;
  285 + m_ogip_tstop = app.m_ogip_tstop;
  286 + m_ogip_telapse = app.m_ogip_telapse;
  287 + m_ogip_exposure = app.m_ogip_exposure;
  288 + m_ogip_ontime = app.m_ogip_ontime;
  289 + m_ogip_livetime = app.m_ogip_livetime;
229 290 m_index_unbinned = app.m_index_unbinned;
230 291  
231 292 // Return
... ... @@ -258,6 +319,15 @@ const GCTAObservation* ctobservation::first_unbinned_observation(void) const
258 319 // Initialise index
259 320 m_index_unbinned = 0;
260 321  
  322 + // Initialise OGIP members
  323 + m_ogip_telescope.clear();
  324 + m_ogip_tstart.clear();
  325 + m_ogip_tstop.clear();
  326 + m_ogip_telapse = 0.0;
  327 + m_ogip_exposure = 0.0;
  328 + m_ogip_ontime = 0.0;
  329 + m_ogip_livetime = 0.0;
  330 +
261 331 // Get next unbinned CTA observation
262 332 const GCTAObservation* obs = next_unbinned_observation();
263 333  
... ... @@ -323,6 +393,31 @@ const GCTAObservation* ctobservation::next_unbinned_observation(void) const
323 393  
324 394 } // endfor: looped over all observation
325 395  
  396 + // Update OGIP members
  397 + if (obs != NULL) {
  398 + if (m_ogip_telapse == 0.0) {
  399 + m_ogip_telescope = obs->instrument();
  400 + m_ogip_tstart = obs->gti().tstart();
  401 + m_ogip_tstop = obs->gti().tstop();
  402 + m_ogip_telapse = obs->gti().telapse();
  403 + m_ogip_exposure = obs->ontime();
  404 + m_ogip_ontime = obs->ontime();
  405 + m_ogip_livetime = obs->livetime();
  406 + }
  407 + else {
  408 + if (obs->gti().tstart() < m_ogip_tstart) {
  409 + m_ogip_tstart = obs->gti().tstart();
  410 + }
  411 + if (obs->gti().tstop() > m_ogip_tstop) {
  412 + m_ogip_tstop = obs->gti().tstop();
  413 + }
  414 + m_ogip_telapse += obs->gti().telapse();
  415 + m_ogip_exposure += obs->ontime();
  416 + m_ogip_ontime += obs->ontime();
  417 + m_ogip_livetime += obs->livetime();
  418 + }
  419 + }
  420 +
326 421 // Return next CTA observation
327 422 return obs;
328 423 }
... ...
src/support/ctobservation.hpp
  View file @3b6811f
... ... @@ -63,6 +63,7 @@ public:
63 63  
64 64 // Methods
65 65 const GObservations& obs(void) const;
  66 + void write_ogip_keywords(GFitsHDU* hdu) const;
66 67  
67 68 #ifndef SWIG
68 69 protected:
... ... @@ -82,11 +83,18 @@ protected:
82 83 void save_events_xml(void);
83 84  
84 85 // Protected members
85   - GObservations m_obs; //!< Observation container
  86 + GObservations m_obs; //!< Observation container
86 87  
87 88 private:
88 89 // Private members
89   - mutable int m_index_unbinned; //!< Current index of unbinned observation
  90 + mutable std::string m_ogip_telescope; //!< Name of telescope
  91 + mutable GTime m_ogip_tstart; //!< Start time for OGIP keywords
  92 + mutable GTime m_ogip_tstop; //!< Stop time for OGIP keywords
  93 + mutable double m_ogip_telapse; //!< Elapsed time
  94 + mutable double m_ogip_exposure; //!< Exposure time
  95 + mutable double m_ogip_ontime; //!< Ontime for OGIP keywords
  96 + mutable double m_ogip_livetime; //!< Livetime for OGIP keywords
  97 + mutable int m_index_unbinned; //!< Current index of unbinned observation
90 98 };
91 99  
92 100  
... ...