Revise handling of MET time references.

Remove ctools::time_reference constant and handle the time reference for MET either by using the time reference from the input observation, or by adding a hidden "mjdref" parameter. The following tools transparently use the time reference of the input observation: * ctselect (apply time reference of each input observation) * ctbin (apply time reference of first input observation) * csobsselect (apply time reference of each input observation) A "mjdref" parameter was added the following tools: * ctobssim * cslightcrv * cspull * csviscube * csobsdef

Revise handling of MET time references.
Remove ctools::time_reference constant and handle the time reference for MET either by using the time reference from the input observation, or by adding a hidden "mjdref" parameter. The following tools transparently use the time reference of the input observation: * ctselect (apply time reference of each input observation) * ctbin (apply time reference of first input observation) * csobsselect (apply time reference of each input observation) A "mjdref" parameter was added the following tools: * ctobssim * cslightcrv * cspull * csviscube * csobsdef
Jürgen Knödlseder
1 parent 2624e245
Showing 28 changed files with 170 additions and 103 deletions
ChangeLog
NEWS
cscripts/cslightcrv.par
cscripts/cslightcrv.py
cscripts/csobsdef.par
cscripts/csobsdef.py
cscripts/csobsinfo.py
cscripts/csobsselect.py
cscripts/cspull.par
cscripts/csviscube.par
cscripts/csviscube.py
cscripts/obsutils.py
doc/source/admin/release_history/1.6.rst
doc/source/users/reference_manual/cslightcrv.rst
doc/source/users/reference_manual/csobsdef.rst
doc/source/users/reference_manual/csobsselect.rst
doc/source/users/reference_manual/cspull.rst
doc/source/users/reference_manual/csviscube.rst
doc/source/users/reference_manual/ctobssim.rst
doc/source/users/reference_manual/ctselect.rst
-2018-06-12
+2018-06-13
         * Version 1.6.0 released
         ========================
@@ -12,6 +12,8 @@
         Allow specifying of input sky map in ctskymap (#2311)
         Add ras() and decs() methods to csobsinfo (#2313)
         Accelerate IRF and RING background methods in ctskymap (#2309)
+        Define MET with respect to the observation time reference or "mjdref"
+          parameter (#2054)
 2018-04-26
 New Features and Important Changes in ctools 1.6.0
-12 June 2018
+13 June 2018
 Introduction
@@ -26,13 +26,25 @@ None
 Base classes
 ------------
+The ctool::time_reference constant was removed since it is not required
+anymore. The time reference is now transparently handled by the tools, and
+if there is a choice, a "mjdref" parameter was added to the tools (#2054).
+
+The ctool::create_cta_obs() method now queries the optional "mjdref" parameter
+so that a MJD reference time can be specified by the tool. This functionality
+is now used by ctobssim (#2054).
+
+The ctools::get_gti() method now takes GTimeReference as argument to enable
+tools to retrieve Mission Elapsed Time (MET) in the time reference of the
+observation (#2054).
+
 The ctools::is_valid_filename() method was removed since the validity check
 of the filename is now done by the GApplicationPar class.
 ctbin - CTA event binning
 -------------------------
-None
+Copy the time reference from the first input observation (#2054).
 ctbkgcube - CTA background cube generation
@@ -82,7 +94,8 @@ None
 ctobssim - CTA observation simulation
 -------------------------------------
-None
+Add "mjdref" parameter to allow specification of the MJD reference time
+for the simulation (#2054).
 ctphase - Event phase computation
@@ -102,7 +115,8 @@ None
 ctselect - CTA event selection
 ------------------------------
-None
+Interpret MET specified by "tmin" and "tmax" parameters in the time reference
+of the input observations (#2054). 
 ctskymap - CTA sky mapping
@@ -172,7 +186,8 @@ None
 cslightcrv - Generate light curve
 ---------------------------------
-None
+Add hidden "mjdref" parameter for the specification of the MET reference
+(#2054).
 csmodelinfo - Return model information
@@ -202,7 +217,8 @@ None
 csobsdef - Generate observation definition XML file
 ---------------------------------------------------
-None
+Add hidden "mjdref" parameter for the specification of the MET reference
+(#2054).
 csobsinfo - Return observation definition information
@@ -212,7 +228,8 @@ Add ras() and decs() methods (#2313)
 csobsselect - Select observations
 ---------------------------------
-None
+Use time reference of input observations if time selection is given in MET
+(#2054).
 csphagen - Generate On/Off observations
@@ -233,7 +250,8 @@ None
 cspull - Generate pull distributions
 ------------------------------------
-None
+Add hidden "mjdref" parameter for the specification of the MET reference
+(#2054).
 csresmap - Generate residual map
@@ -284,7 +302,8 @@ None
 csviscube - Generate visibility cube
 ------------------------------------
-None
+Add hidden "mjdref" parameter for the specification of the MET reference
+(#2054).
 csworkflow - Execute workflow
@@ -309,7 +328,7 @@ None
 obsutils - Observation utilities
 --------------------------------
-None
+Add "mjdref" parameter to set_obs() function (#2054).
 Examples
@@ -31,6 +31,7 @@ outfile,     f, a, lightcurve.fits,,, &quot;Output light curve file&quot;
 tbinalg,   s, a, GTI,FILE|LIN|GTI,, "Algorithm for defining time bins"
 tmin,      t, a, MJD 51544.5,,,"Lightcurve start time (UTC string, JD, MJD or MET in seconds)"
 tmax,      t, a, MJD 51544.6,,,"Lightcurve stop time (UTC string, JD, MJD or MET in seconds)"
+mjdref,    r, h, 51544.5,,, "Reference MJD for MET"
 tbins,     i, a, 5,1,10000, "Number of time bins"
 tbinfile,  f, a, NONE,,, "File defining the time binning"
 method,    s, a, 3D, 3D|ONOFF,, "Analysis method"
@@ -137,8 +137,11 @@ class cslightcrv(ctools.csobservation):
         gti : `~gammalib.GGti`
             Light curve bins in form of Good Time Intervals
         """
+        # Get MET time reference
+        tref = gammalib.GTimeReference(self['mjdref'].real(),'s','TT','LOCAL')
+
         # Initialise Good Time Intervals
-        gti = gammalib.GGti()
+        gti = gammalib.GGti(tref)
         # Get algorithm to use for defining the time intervals. Possible
         # values are "FILE", "LIN" and "GTI". This is enforced at
@@ -173,8 +176,8 @@ class cslightcrv(ctools.csobservation):
         elif algorithm == 'LIN':
             # Get start and stop time and number of time bins
-            time_min = self['tmin'].time(ctools.time_reference)
-            time_max = self['tmax'].time(ctools.time_reference)
+            time_min = self['tmin'].time(tref)
+            time_max = self['tmax'].time(tref)
             nbins    = self['tbins'].integer()
             # Compute time step in seconds and setup the GTIs
@@ -385,8 +388,8 @@ class cslightcrv(ctools.csobservation):
         select = ctools.ctselect(self.obs())
         select['emin'] = self['emin'].real()
         select['emax'] = self['emax'].real()
-        select['tmin'] = tmin.convert(ctools.time_reference)
-        select['tmax'] = tmax.convert(ctools.time_reference)
+        select['tmin'] = tmin
+        select['tmax'] = tmax
         select['rad']  = 'UNDEFINED'
         select['ra']   = 'UNDEFINED'
         select['dec']  = 'UNDEFINED'
@@ -444,8 +447,8 @@ class cslightcrv(ctools.csobservation):
             select = ctools.ctselect(self.obs())
             select['emin'] = self['emin'].real()
             select['emax'] = self['emax'].real()
-            select['tmin'] = tmin.convert(ctools.time_reference)
-            select['tmax'] = tmax.convert(ctools.time_reference)
+            select['tmin'] = tmin
+            select['tmax'] = tmax
             select['rad']  = 'UNDEFINED'
             select['ra']   = 'UNDEFINED'
             select['dec']  = 'UNDEFINED'
@@ -30,6 +30,7 @@ emin,     r, h, UNDEF,,, &quot;Lower energy limit (TeV)&quot;
 emax,     r, h, UNDEF,,, "Upper energy limit (TeV)"
 rad,      r, h, UNDEF,,, "Radius of ROI (degrees)"
 deadc,    r, h, 0.98,,, "Deadtime correction factor"
+mjdref,   r, h, 51544.5,,, "Reference MJD for observations"
 #
 # Standard parameters
@@ -272,11 +272,12 @@ class csobsdef(ctools.cscript):
                 duration = float(self._pntdef[row, header.index('duration')])
             else:
                 duration = self['duration'].real()
+            tref       = gammalib.GTimeReference(self['mjdref'].real(),'s')
             tmin       = self._tmin
             tmax       = self._tmin + duration
-            gti        = gammalib.GGti(ctools.time_reference)
-            tstart     = gammalib.GTime(tmin, ctools.time_reference)
-            tstop      = gammalib.GTime(tmax, ctools.time_reference)
+            gti        = gammalib.GGti(tref)
+            tstart     = gammalib.GTime(tmin, tref)
+            tstop      = gammalib.GTime(tmax, tref)
             self._tmin = tmax
             gti.append(tstart, tstop)
             obs.ontime(gti.ontime())
@@ -289,10 +289,11 @@ class csobsinfo(ctools.csobservation):
         self._log_value(gammalib.NORMAL, 'Maximum energy', max_value)
         # Log time range
+        ref = obs.gti().reference()
         mjd = '%.3f - %.3f' % (self._gti.tstart().mjd(),  self._gti.tstop().mjd())
         utc = '%s - %s'     % (self._gti.tstart().utc(),  self._gti.tstop().utc())
-        met = '%.3f - %.3f' % (self._gti.tstart().convert(ctools.time_reference),
-                               self._gti.tstop().convert(ctools.time_reference))
+        met = '%.3f - %.3f' % (self._gti.tstart().convert(ref),
+                               self._gti.tstop().convert(ref))
         self._log_header3(gammalib.NORMAL, 'Time range')
         self._log_value(gammalib.NORMAL, 'MJD (days)', mjd)
         self._log_value(gammalib.NORMAL, 'UTC',  utc)
@@ -68,8 +68,10 @@ class csobsselect(ctools.csobservation):
             self['width'].real()
             self['height'].real()
-        # Query time interval
-        self._get_gti()
+        # Query time interval. The GammaLib time reference is specified as a
+        # dummy argument since the relevant intervals will be queried later
+        # using the time reference for each observation
+        self._get_gti(gammalib.GTimeReference())
         # Query ahead output model filename
         if self._read_ahead():
@@ -113,7 +115,7 @@ class csobsselect(ctools.csobservation):
         # Return
         return
-    def _select_observation(self, obs, gti):
+    def _select_observation(self, obs):
         """
         Select observation
@@ -121,8 +123,6 @@ class csobsselect(ctools.csobservation):
         ----------
         obs : `~gammalib.GObservation`
             Observation
-        gti : `~gammalib.GGti`
-            Good time intervals
         Returns
         -------
@@ -140,6 +140,10 @@ class csobsselect(ctools.csobservation):
             tstart  = obs_gti.tstart()
             tstop   = obs_gti.tstop()
+            # Get User start and stop time in the time reference of the CTA
+            # observations
+            gti = self._get_gti(obs_gti.reference())
+
             # If the observation is not within the time selection then skip it
             if gti.size() > 0 and not gti.contains(tstart) and \
                                   not gti.contains(tstop):
@@ -294,9 +298,6 @@ class csobsselect(ctools.csobservation):
         # Get parameters
         self._get_parameters()
-        # Get Good Time Interval (empty if no time selection is applied)
-        gti = self._get_gti()
-
         # Initialise empty observation container for selected observations
         selected_obs = gammalib.GObservations()
@@ -310,7 +311,7 @@ class csobsselect(ctools.csobservation):
         for obs in self.obs():
             # If observation is selected then append observation
-            if self._select_observation(obs, gti):
+            if self._select_observation(obs):
                 selected_obs.append(obs)
         # Copy selected observations into observation
@@ -40,6 +40,7 @@ emax,      r, a, 100.0,,, &quot;Upper energy limit (TeV)&quot;
 enumbins,  i, a, 0,,, "Number of energy bins (0=unbinned)"
 tmin,      t, a, 2020-01-01T00:00:00,,, "Start time (UTC string, JD, MJD or MET in seconds)"
 tmax,      t, a, 2020-01-01T00:30:00,,, "Stop time (UTC string, JD, MJD or MET in seconds)"
+mjdref,    r, h, 51544.5,,, "Reference MJD for MET"
 npix,      i, a, 200,,, "Number of pixels for binned"
 binsz,     r, a, 0.05,,, "Pixel size for binned (degrees/pixel)"
 statistic, s, h, DEFAULT,DEFAULT|CSTAT|WSTAT|CHI2,, "Optimization statistic"
@@ -18,6 +18,7 @@
 #===================
 tmin,       t, a, 2020-01-01T00:00:00,,, "Start time of observation (UTC string, JD, MJD or MET in seconds)"
 tmax,       t, a, 2021-01-01T00:00:00,,, "Stop time of observation (UTC string, JD, MJD or MET in seconds)"
+mjdref,     r, h, 51544.5,,, "Reference MJD for MET"
 geolon,     r, a,  79.4041,,, "Geographic longitude of array (deg)"
 geolat,     r, a, -24.6272,,, "Geographic latitude of array (deg)"
 sunzenith,  r, h, 105.0,,, "Minimum Sun zenith angle for observations (deg)"
@@ -64,6 +64,7 @@ class csviscube(ctools.cscript):
         Get all parameters
         """
         # Query parameters
+        self['mjdref'].real()
         self['tmin'].time()
         self['tmax'].time()
         self['geolon'].real()
@@ -259,9 +260,12 @@ class csviscube(ctools.cscript):
         # Get array geographic longitude
         #geolon = self['geolon'].real()
+        # Get MET time reference
+        tref = gammalib.GTimeReference(self['mjdref'].real(),'s','TT','LOCAL')
+
         # Get time interval
-        tmin = self['tmin'].time(ctools.time_reference)
-        tmax = self['tmax'].time(ctools.time_reference)
+        tmin = self['tmin'].time(tref)
+        tmax = self['tmax'].time(tref)
         # Initialise hour angle list
         hour_angles = []
@@ -224,7 +224,7 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0,
 # ======================= #
 def set_obs(pntdir, tstart=0.0, duration=1800.0, deadc=0.98, \
             emin=0.1, emax=100.0, rad=5.0, \
-            irf='South_50h', caldb='prod2', obsid='000000'):
+            irf='South_50h', caldb='prod2', obsid='000000', mjdref=51544.5):
     """
     Set a single CTA observation
@@ -283,9 +283,10 @@ def set_obs(pntdir, tstart=0.0, duration=1800.0, deadc=0.98, \
     roi.radius(rad)
     # Set GTI
-    gti = gammalib.GGti(ctools.time_reference)
-    gti.append(gammalib.GTime(tstart, ctools.time_reference),
-               gammalib.GTime(tstart+duration, ctools.time_reference))
+    ref = gammalib.GTimeReference(mjdref,'s','TT','LOCAL')
+    gti = gammalib.GGti(ref)
+    gti.append(gammalib.GTime(tstart, ref),
+               gammalib.GTime(tstart+duration, ref))
     # Set energy boundaries
     ebounds = gammalib.GEbounds(gammalib.GEnergy(emin, 'TeV'),
@@ -36,6 +36,9 @@ Improvements
   Allow specifying of input skymap to :ref:`ctskymap` tool
 * [`2309 <https://cta-redmine.irap.omp.eu/issues/2309>`_] -
   Accelerate `IRF` and `RING` background methods in :ref:`ctskymap`
+* [`2045 <https://cta-redmine.irap.omp.eu/issues/2045>`_] -
+  Define MET with respect to the observation time reference or the ``mjdref``
+  user parameter (#2054)
 New features
@@ -48,3 +51,5 @@ Clean up
 --------
 * Remove ``ctools::is_valid_filename()`` method
+* [`2045 <https://cta-redmine.irap.omp.eu/issues/2045>`_] -
+  Remove ``ctools::time_reference`` constant
@@ -77,6 +77,9 @@ General parameters
 ``tmax [time]``
     Lightcurve stop time (UTC string, JD, MJD or MET in seconds).
+``(mjdref = 51544.5) [real]``
+    Reference Modified Julian Day (MJD) for Mission Elapsed Time (MET).
+
 ``tbins [integer]``
     Number of time bins.
@@ -82,6 +82,9 @@ General parameters
 ``(deadc = 0.95) [real]``
     Average deadtime correction factor.
+``(mjdref = 51544.5) [real]``
+    Reference Modified Julian Day (MJD) for observations.
+
 Standard parameters
 -------------------
@@ -52,12 +52,16 @@ General parameters
     Height of selection box (degrees).
 ``tmin [time]``
-    Start time (UTC string, JD, MJD or time in seconds).
+    Start time for observation selection (UTC string, JD, MJD or MET in seconds).
+    Start times given in MET seconds are counted with respect to the time
+    reference of the input observation(s).
     If ``INDEF``, ``NONE``, ``UNDEF`` or ``UNDEFINED`` is passed as value, no time
     selection will be performed.
 ``tmax [time]``
-    Stop time (UTC string, JD, MJD or time in seconds).
+    Stop time for observation selection (UTC string, JD, MJD or MET in seconds).
+    Stop times given in MET seconds are counted with respect to the time
+    reference of the input observation(s).
     If ``INDEF``, ``NONE``, ``UNDEF`` or ``UNDEFINED`` is passed as value, no time
     selection will be performed.
@@ -106,6 +106,9 @@ General parameters
 ``tmax [time]``
     Stop time (UTC string, JD, MJD or MET in seconds).
+``(mjdref = 51544.5) [real]``
+    Reference Modified Julian Day (MJD) for MET.
+
 ``npix [integer]``
     Number of pixels for binned analysis.
@@ -30,6 +30,9 @@ General parameters
 ``tmax [time]``
     Stop time of observation (UTC string, JD, MJD or MET in seconds).
+``(mjdref = 51544.5) [real]``
+    Reference Modified Julian Day (MJD) for MET.
+
 ``geolon [real]``
     Geographic longitude of the IACT array in degrees.
@@ -89,6 +89,10 @@ General parameters
 ``tmax [time]``
     Stop time (UTC string, JD, MJD or MET in seconds).
+``(mjdref = 51544.5) [real]``
+    Reference Modified Julian Day (MJD) for simulated events. The times in
+    seconds for each event are counted from this reference time on.
+
 ``emin [real]``
     Lower energy limit of simulated events (in TeV).
@@ -66,12 +66,16 @@ General parameters
     selection will be performed.
 ``tmin [time]``
-    Start time (UTC string, JD, MJD or time in seconds).
+    Start time for event selection (UTC string, JD, MJD or MET in seconds).
+    Start times given in MET seconds are counted with respect to the time
+    reference of the input observation(s).
     If ``INDEF``, ``NONE``, ``UNDEF`` or ``UNDEFINED`` is passed as value, no time
     selection will be performed.
 ``tmax [time]``
-    Stop time (UTC string, JD, MJD or time in seconds).
+    Stop time for event selection (UTC string, JD, MJD or MET in seconds).
+    Stop times given in MET seconds are counted with respect to the time
+    reference of the input observation(s).
     If ``INDEF``, ``NONE``, ``UNDEF`` or ``UNDEFINED`` is passed as value, no time
     selection will be performed.
@@ -57,10 +57,6 @@ public:
 %include "argcargv.i"
 /* __ Constants __________________________________________________________ */
-namespace ctools {
-    const GTimeReference time_reference =
-          GTimeReference(G_CTA_MJDREF, "s", "TT", "LOCAL");
-}
 /***********************************************************************//**
@@ -142,7 +138,7 @@ public:
     // Protected methods that extract user parameters
     GCTARoi           get_roi(const GCTAPointing& pnt = GCTAPointing());
     GEbounds          get_ebounds(void);
-    GGti              get_gti(void);
+    GGti              get_gti(const GTimeReference& ref);
     GCTAPointing      get_pointing(void);
     GSkyDir           get_skydir(void);
     std::string       set_outfile_name(const std::string& filename);
@@ -385,9 +385,6 @@ void ctbin::init_members(void)
     // Initialise cache members
     m_dirs.clear();
-    // Set CTA reference time
-    m_gti.reference(ctools::time_reference);
-
     // Return
     return;
 }
@@ -586,6 +583,12 @@ void ctbin::fill_cube(GCTAObservation* obs)
     // Update time information
     #pragma omp critical(ctbin_fill_cube)
     {
+        // If GTI is empty then set its time reference from the observation.
+        // From then on we keep that time reference
+        if (m_gti.is_empty()) {
+            m_gti.reference(events->gti().reference());
+        }
+
         // Append GTIs
         m_gti.extend(events->gti());
@@ -784,9 +784,10 @@ void ctobssim::simulate_source(GCTAObservation* obs,
         // Loop over all Good Time Intervals
         for (int it = 0; it < events->gti().size(); ++it) {
-            // Extract time interval
-            GTime tmin = events->gti().tstart(it);
-            GTime tmax = events->gti().tstop(it);
+            // Extract time interval and time reference
+            GTime          tmin = events->gti().tstart(it);
+            GTime          tmax = events->gti().tstop(it);
+            GTimeReference tref = events->gti().reference();
             // Log time interval. Increment indentation if there are
             // several Good Time Intervals.
@@ -796,9 +797,9 @@ void ctobssim::simulate_source(GCTAObservation* obs,
                     wrklog->indent(indent);
                 }
                 *wrklog << gammalib::parformat("Time interval", indent);
-                *wrklog << tmin.convert(ctools::time_reference);
+                *wrklog << tmin.convert(tref);
                 *wrklog << " - ";
-                *wrklog << tmax.convert(ctools::time_reference);
+                *wrklog << tmax.convert(tref);
                 *wrklog << " s" << std::endl;
             }
@@ -1063,8 +1064,8 @@ void ctobssim::simulate_interval(GCTAObservation*       obs,
                     wrklog->indent(indent);
                 }
                 *wrklog << gammalib::parformat("Time slice", indent);
-                *wrklog << tstart.convert(ctools::time_reference) << " - ";
-                *wrklog << tstop.convert(ctools::time_reference) << " s";
+                *wrklog << tstart.convert(obs->gti().reference()) << " - ";
+                *wrklog << tstop.convert(obs->gti().reference()) << " s";
                 if (model->name().length() > 0) {
                     *wrklog << " [" << model->name() << "]";
                 }
@@ -34,6 +34,7 @@ dec,     r, a, 22.51,-90,90, &quot;Dec of pointing (degrees)&quot;
 rad,     r, a, 5.0,0,180, "Radius of FOV (degrees)"
 tmin,    t, a, 2020-01-01T00:00:00,,, "Start time (UTC string, JD, MJD or MET in seconds)"
 tmax,    t, a, 2020-01-01T00:30:00,,, "Stop time (UTC string, JD, MJD or MET in seconds)"
+mjdref,  r, h, 51544.5,,, "Reference MJD for simulated events"
 emin,    r, a, 0.1,,, "Lower energy limit (TeV)"
 emax,    r, a, 100.0,,, "Upper energy limit (TeV)"
 deadc,   r, h, 0.98,0,1, "Average deadtime correction factor"
@@ -465,7 +465,6 @@ void ctselect::init_members(void)
 {
     // Initialise parameters
     m_outobs.clear();
-    m_gti.clear();
     m_emin   = 0.0;
     m_emax   = 0.0;
     m_expr.clear();
@@ -476,8 +475,6 @@ void ctselect::init_members(void)
     m_infiles.clear();
     m_evtname.clear();
     m_gtiname.clear();
-    m_timemin.clear();
-    m_timemax.clear();
     m_phases.clear();
     m_select_energy = false;
     m_select_phase  = false;
@@ -496,7 +493,6 @@ void ctselect::copy_members(const ctselect&amp; app)
 {
     // Copy parameters
     m_outobs   = app.m_outobs;
-    m_gti      = app.m_gti;
     m_emin     = app.m_emin;
     m_emax     = app.m_emax;
     m_expr     = app.m_expr;
@@ -507,8 +503,6 @@ void ctselect::copy_members(const ctselect&amp; app)
     m_infiles       = app.m_infiles;
     m_evtname       = app.m_evtname;
     m_gtiname       = app.m_gtiname;
-    m_timemin       = app.m_timemin;
-    m_timemax       = app.m_timemax;
     m_phases        = app.m_phases;
     m_select_energy = app.m_select_energy;
     m_select_phase  = app.m_select_phase;
@@ -555,8 +549,10 @@ void ctselect::get_parameters(void)
     GCTAPointing pnt;
     GCTARoi      roi = get_roi(pnt);
-    // Get Good Time Intervals (empty if no selection is done)
-    m_gti = get_gti();
+    // Query Good Time Intervals. The GammaLib time reference is specified
+    // as a dummy argument since the relevant intervals will be queried later
+    // using the time reference for each observation
+    get_gti(GTimeReference());
     // Get energy selection parameters
     if ((*this)["emin"].is_valid() && (*this)["emax"].is_valid()) {
@@ -653,11 +649,6 @@ void ctselect::get_parameters(void)
         (*this)["prefix"].query();
     }
-    // Set time interval with input times given in CTA reference
-    // time (in seconds)
-    m_timemin = m_gti.tstart();
-    m_timemax = m_gti.tstop();
-
     // Write parameters into logger
     log_parameters(TERSE);
@@ -685,7 +676,7 @@ void ctselect::get_parameters(void)
  * observation is replaced by selected event list read from the FITS file.
  *
  * Good Time Intervals of the observation will be limited to the time
- * interval [m_timemin, m_timemax].
+ * interval specified by the User parameters tmin and tmax.
  ***************************************************************************/
 void ctselect::select_events(GCTAObservation*   obs,
                              const std::string& filename,
@@ -706,6 +697,12 @@ void ctselect::select_events(GCTAObservation*   obs,
     GCTAEventList* list =
         static_cast<GCTAEventList*>(const_cast<GEvents*>(obs->events()));
+    // Get User GTI with MET in the time reference that is specified for the
+    // observation
+    GGti  user_gti = get_gti(list->gti().reference());
+    GTime timemin  = user_gti.tstart();
+    GTime timemax  = user_gti.tstop();
+
     // Get existing Roi and energy bounds for possible later use
     // (will be empty if unavailable)
     GCTARoi  old_roi     = list->roi();
@@ -736,12 +733,12 @@ void ctselect::select_events(GCTAObservation*   obs,
     // Set time selection interval. We make sure here that the time selection
     // interval cannot be wider than the GTIs covering the data. This is done
     // using GGti's reduce() method.
-    if (!m_gti.is_empty()) {
+    if (!user_gti.is_empty()) {
         // Reduce GTIs to specified time interval. The complicated cast is
         // necessary here because the gti() method is declared const, so
         // we're not officially allowed to modify the GTIs.
-        ((GGti*)(&list->gti()))->reduce(m_timemin, m_timemax);
+        ((GGti*)(&list->gti()))->reduce(timemin, timemax);
     } // endif: time selection was required
@@ -749,7 +746,7 @@ void ctselect::select_events(GCTAObservation*   obs,
     GGti gti = list->gti();
     // Make time selection
-    if (!m_gti.is_empty()) {
+    if (!user_gti.is_empty()) {
         // If there are no GTIs then request removal of all events
         if (list->gti().is_empty()) {
 /***************************************************************************
  *                      ctselect - Data selection tool                     *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2010-2017 by Juergen Knoedlseder                         *
+ *  copyright (C) 2010-2018 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -81,7 +81,7 @@ protected:
     // User parameters
     std::string m_outobs;     //!< Output event list or XML file
-    GGti        m_gti;        //!< Time selection
+    //GGti        m_gti;        //!< Time selection
     double      m_emin;       //!< Lower energy
     double      m_emax;       //!< Upper energy
     std::string m_expr;       //!< Selection expression
@@ -92,8 +92,8 @@ protected:
     std::vector<std::string> m_infiles;       //!< Input event filenames
     std::vector<std::string> m_evtname;       //!< Event extension names
     std::vector<std::string> m_gtiname;       //!< GTI extension names
-    GTime                    m_timemin;       //!< Earliest time
-    GTime                    m_timemax;       //!< Latest time
+    //GTime                    m_timemin;       //!< Earliest time
+    //GTime                    m_timemax;       //!< Latest time
     GPhases                  m_phases;        //!< Phase intervals
     bool                     m_select_energy; //!< Perform energy selection
     bool                     m_select_phase;  //!< Perform phase selection
@@ -751,14 +751,18 @@ GCTAEventCube ctool::create_cube(const GObservations&amp; obs)
  *
  * The following parameters are read
  *
- *      ra:    Right Ascension of pointing and RoI centre (deg)
- *      dec:   Declination of pointing and RoI centre (deg)
- *      rad:   Radius of RoI (deg)
- *      deadc: Deadtime correction factor
- *      tmin:  Start time
- *      tmax:  Stop time
- *      emin:  Minimum energy (TeV)
- *      emax:  Maximum energy (TeV)
+ *      ra:     Right Ascension of pointing and RoI centre (deg)
+ *      dec:    Declination of pointing and RoI centre (deg)
+ *      rad:    Radius of RoI (deg)
+ *      deadc:  Deadtime correction factor
+ *      tmin:   Start time
+ *      tmax:   Stop time
+ *      emin:   Minimum energy (TeV)
+ *      emax:   Maximum energy (TeV)
+ *      mjdref: Time reference (optional)
+ *
+ * If the time reference parametere "mjdref" is not available, the CTA time
+ * reference will be assumed.
  ***************************************************************************/
 GCTAObservation ctool::create_cta_obs(void)
 {
@@ -775,8 +779,13 @@ GCTAObservation ctool::create_cta_obs(void)
     // Set RoI
     GCTARoi roi = get_roi(pnt);
+    // Set time reference
+    GTimeReference ref = (has_par("mjdref"))
+                       ? GTimeReference((*this)["mjdref"].real(), "s", "TT", "LOCAL")
+                       : GTimeReference(G_CTA_MJDREF, "s", "TT", "LOCAL");
+
     // Set GTI
-    GGti gti = get_gti();
+    GGti gti = get_gti(ref);
     // Set energy boundaries
     GEbounds ebounds = get_ebounds();
@@ -1057,6 +1066,7 @@ GEbounds ctool::get_ebounds(void)
 /***********************************************************************//**
  * @brief Return Good Time Intervals from User parameter
  *
+ * @param[in] ref Time reference.
  * @return Good Time Intervals.
  *
  * Returns Good Time Intervals from the User parameters `tmin` and `tmax`,
@@ -1064,17 +1074,18 @@ GEbounds ctool::get_ebounds(void)
  * If one of `tmin` or `tmax` is not valid, an empty Good Time Interval
  * object is returned. If `tmin` is not valid, `tmax` is not queried.
  ***************************************************************************/
-GGti ctool::get_gti(void)
+GGti ctool::get_gti(const GTimeReference& ref)
 {
-    // Initialise Good Time Intervals with CTA reference time
-    GGti gti(ctools::time_reference);
+    // Initialise Good Time Intervals with reference time
+    GGti gti(ref);
     // Continue only if "tmin" has a valid time value
     if ((*this)["tmin"].is_valid() && (*this)["tmax"].is_valid()) {
-        // Get minimum and maximum time for CTA time reference
-        GTime tmin = (*this)["tmin"].time(ctools::time_reference);
-        GTime tmax = (*this)["tmax"].time(ctools::time_reference);
+        // Get minimum and maximum time using the time reference for MET
+        // times
+        GTime tmin = (*this)["tmin"].time(ref);
+        GTime tmax = (*this)["tmax"].time(ref);
         // Append GTI
         gti.append(tmin, tmax);
@@ -1630,11 +1641,8 @@ GObservations ctool::get_observations(const bool&amp; get_response)
         // and load file into observation
         if (GFilename(filename).is_fits()) {
-            // Allocate empty CTA observation
-            GCTAObservation cta_obs;
-
-            // Load data
-            cta_obs.load(filename);
+            // Load CTA observation
+            GCTAObservation cta_obs(filename);
             // Get response if required
             if (get_response) {
@@ -34,10 +34,6 @@
 /* __Definitions _________________________________________________________ */
 /* __ Constants __________________________________________________________ */
-namespace ctools {
-    const GTimeReference time_reference =
-          GTimeReference(G_CTA_MJDREF, "s", "TT", "LOCAL");
-}
 /***********************************************************************//**
@@ -102,7 +98,7 @@ protected:
     // Protected methods that extract user parameters
     GCTARoi           get_roi(const GCTAPointing& pnt = GCTAPointing());
     GEbounds          get_ebounds(void);
-    GGti              get_gti(void);
+    GGti              get_gti(const GTimeReference& ref);
     GCTAPointing      get_pointing(void);
     GSkyDir           get_skydir(void);
     std::string       set_outfile_name(const std::string& filename);