Added warning and updated test cases for ctphase (#1982)

Added a warning in ctphase that is written into the log file when either the first or second derivative term in the phase computation is above 1e8, as at this point precision issues can manifest in the computed phase values. Updated test cases and fixed some bugs that allow the tests of ctphase to pass.

Added warning and updated test cases for ctphase (#1982)
Added a warning in ctphase that is written into the log file when either the first or second derivative term in the phase computation is above 1e8, as at this point precision issues can manifest in the computed phase values. Updated test cases and fixed some bugs that allow the tests of ctphase to pass.
Josh Cardenzana
1 parent a1db5073
Showing 3 changed files with 30 additions and 16 deletions
src/ctphase/ctphase.cpp
src/ctphase/ctphase.par
test/test_ctphase.py
@@ -633,16 +633,30 @@ void ctphase::phase_events(GCTAObservation*   obs,
                              const std::string& evtname,
                              const std::string& gtiname)
 {
-    // Write header into logger
-    log_header3(NORMAL, "Events phasing");
-
     // Get CTA event list pointer
     GCTAEventList* list =
         static_cast<GCTAEventList*>(const_cast<GEvents*>(obs->events()));
+    GCTAEventAtom* event = (*list)[0];
+    
+    // Print a warning into the log file if there are some precision concerns
+    double dt = std::abs(event->time() - m_phase.mjd());
+    if (((m_phase.f1() * dt) > 1.0e8) ||
+        (m_phase.f2() * dt * dt) > 1.0e8) {
+        log_string(NORMAL,"*** WARNING ***");
+        log_string(NORMAL,"   Values supplied for reference MJD and/or f1 and/or f2");
+        log_string(NORMAL,"   are large and may result in numerical precision issues.");
+        log_string(NORMAL,"   Consider using an ephemeris derived from a date closer");
+        log_string(NORMAL,"   to the data being analyzed.");
+        log_value(NORMAL, "     Observation ID", obs->obs_id());
+        log_value(NORMAL, "     Reference MJD", m_phase.mjd().mjd());
+        log_value(NORMAL, "     First event MJD", event->time().mjd());
+        log_string(NORMAL,"*** WARNING ***");
+    }
+    
     for (int evnt=0; evnt<list->size(); evnt++) {
         // Get the next event
-        GCTAEventAtom* event = (*list)[evnt];
+        event = (*list)[evnt];
         // Get the event time
         GTime time = event->time();
@@ -24,11 +24,11 @@ prefix, s, h, &quot;phased_&quot;,,, &quot;Prefix for event lists in observation definition XML
 #
 # Phase information parameters
 #=====================
-mjd,    r, a, 52400,,, "Reference MJD for phase information"
-f0,     r, a, 1,,, "Frequency at reference MJD [Hz]"
-f1,     r, a, -0.2,,, "First freq. derivative at reference MJD [s^-1]"
-f2,     r, a, 0,,, "Second freq. derivative at reference MJD [s^-2]"
+mjd,    r, a, 51544.5,,, "Reference MJD for phase information"
 p0,     r, a, 0.1,0.0,1.0, "Phase at reference MJD"
+f0,     r, a, 1,,, "Frequency at reference MJD [Hz]"
+f1,     r, a, 0.0,,, "First freq. derivative at reference MJD [s^-1]"
+f2,     r, a, 0.0,,, "Second freq. derivative at reference MJD [s^-2]"
 #
 # Standard parameters
@@ -71,7 +71,7 @@ class Test(test):
         # Setup ctphase command
         cmd = ctphase+' inobs="'+self._events+'"'+ \
                        ' outobs="ctphase_cmd1.fits"'+ \
-                       ' p0=0.1 mjd=52400'+ \
+                       ' p0=0.1 mjd=51544.5'+ \
                        ' f0=0.1 f1=0. f2=0.'+ \
                        ' logfile="ctphase_cmd1.log" chatter=1'
@@ -89,7 +89,7 @@ class Test(test):
         # Setup ctphase command
         cmd = ctphase+' inobs="event_file_that_does_not_exist.fits"'+ \
                        ' outobs="ctphase_cmd2.fits"'+ \
-                       ' p0=0.1 mjd=52400'+ \
+                       ' p0=0.1 mjd=51544.5'+ \
                        ' f0=0.1 f1=0. f2=0.'+ \
                        ' logfile="ctphase_cmd2.log" chatter=1'
@@ -137,9 +137,9 @@ class Test(test):
         # Now set ctphase parameters
         phase['inobs']   = self._events
         phase['p0']      = 0.
-        phase['mjd']     = 52400
+        phase['mjd']     = 51544.5
         phase['f0']      = 1.
-        phase['f1']      = -0.2
+        phase['f1']      = 0.0
         phase['f2']      = 0.
         phase['outobs']  = 'ctphase_py1.fits'
         phase['logfile'] = 'ctphase_py1.log'
@@ -158,18 +158,18 @@ class Test(test):
         # Execute copy of ctphase tool again, now with a negative value
         # of initial phase to test for invalid parameter value
+        self.test_try('Test invalid parameter value')
         try:
-            self.test_try('Test invalid parameter value')
-            cpy_phase['p0']      = -0.5
-            cpy_phase['mjd']     = 52400
+            cpy_phase['mjd']     = 51544.5
             cpy_phase['f0']      = 1.
-            cpy_phase['f1']      = -0.2
+            cpy_phase['f1']      = 0.2
             cpy_phase['f2']      = 0.
             cpy_phase['outobs']  = 'ctphase_py2.fits'
             cpy_phase['logfile'] = 'ctphase_py2.log'
             cpy_phase['chatter'] = 3
             cpy_phase['publish'] = True
             cpy_phase.logFileOpen()  # Needed to get a new log file
+            cpy_phase['p0']      = -0.5
             cpy_phase.execute()
             self.test_try_failure('Exception not thrown')
         except ValueError: