Veronika Vodeb / ctools

Commit 9955078b9ca0983db5c530be452edcee372610d3

Authored by Jürgen Knödlseder
1 parent 3705c8a9

Add keywords to sky maps produced by ctskymap (#2194)

Inline Side-by-side
Showing 8 changed files with 158 additions and 12 deletions
ChangeLog
  View file @9955078
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 @9955078
... ... @@ -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
... ...
pyext/ctobservation.i
  View file @9955078
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 @9955078
... ... @@ -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 @9955078
... ... @@ -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 @9955078
... ... @@ -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 @9955078
... ... @@ -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  
... ...