Enhance unit test of ctprob and ctphase tools

Jürgen Knödlseder
1 parent 3c202bd5
Showing 4 changed files with 109 additions and 15 deletions
src/ctphase/ctphase.cpp
test/test_ctbutterfly.py
test/test_ctphase.py
test/test_ctprob.py
@@ -504,8 +504,8 @@ void ctphase::get_parameters(void)
  * @exception GException::invalid_value
  *            No events found in observation.
  *
- * Append phase information to events from a FITS file by adding a "PHASE"
- * column to the FITS file.
+ * Computes the phase for all events in an observation based on a temporal
+ * phase curve model.
  ***************************************************************************/
 void ctphase::phase_events(GCTAObservation* obs)
 {
@@ -525,20 +525,29 @@ void ctphase::phase_events(GCTAObservation* obs)
     // Continue only if there are events
     if (events->size() > 0) {
  
-        // Print a warning into the log file if there are some precision concerns
+        // Get MJD information of observation
         GCTAEventAtom* event = (*events)[0];
         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 ***");
+        double         f1    = m_phase.f1() * dt;
+        double         f2    = m_phase.f2() * dt * dt;
+
+        // Log information about MJD of the observation
+        log_value(NORMAL, "Reference MJD", m_phase.mjd().mjd());
+        log_value(NORMAL, "First event MJD", event->time().mjd());
+        log_value(NORMAL, "First event offset", gammalib::str(dt)+" s");
+        log_value(NORMAL, "F1 * dt", gammalib::str(f1)+" Hz");
+        log_value(NORMAL, "F2 * dt * dt", gammalib::str(f2)+" Hz");
+
+        // Print a warning into the log file if there are some precision
+        // concerns
+        if ((f1 > 1.0e8) || (f2 > 1.0e8)) {
+            std::string msg = "\nWARNING: Values supplied for reference MJD "
+                              "and/or F1 and/or F2 are large and may\n"
+                              "         result in numerical precision issues. "
+                              "Consider using an ephemeris derived\n"
+                              "         from a date closer to the data being "
+                              "analyzed.\n";
+            log_string(TERSE, msg);
         }
  
         // Compute the phase for all events
@@ -2,7 +2,7 @@
 # ==========================================================================
 # This scripts performs unit tests for the ctbutterfly tool.
 #
-# Copyright (C) 2014-2016 Michal Mayer
+# Copyright (C) 2014-2017 Michal Mayer
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -144,6 +144,7 @@ class Test(test):
         butterfly['srcname'] = 'Crab'
         butterfly['caldb']   = self._caldb
         butterfly['irf']     = self._irf
+        butterfly['method']  = 'ENVELOPE'
         butterfly['emin']    = 0.1
         butterfly['emax']    = 100.0
         butterfly['fit']     = True
@@ -191,6 +191,72 @@ class Test(test):
         # Check result file
         self._check_result_file('ctphase_py3.fits')
  
+        # Now run ctphase tool without a model file and a significant
+        # MJD offset
+        phase['inobs']   = self._events
+        phase['inmodel'] = 'NONE'
+        phase['srcname'] = 'NONE'
+        phase['mjd']     = 51500.0
+        phase['phase']   = 0.0
+        phase['f0']      = 1.0
+        phase['f1']      = 0.1
+        phase['f2']      = 0.0001
+        phase['outobs']  = 'ctphase_py4.fits'
+        phase['logfile'] = 'ctphase_py4.log'
+        phase['chatter'] = 3
+
+        # Execute ctphase tool
+        phase.logFileOpen()   # Make sure we get a log file
+        phase.execute()
+
+        # Check result file
+        self._check_result_file('ctphase_py4.fits')
+
+        # Test invalid source name
+        self.test_try('Test ctphase with invalid source name')
+        try:
+            test = ctools.ctphase()
+            test['inobs']   = self._events
+            test['inmodel'] = self._model_file
+            test['srcname'] = 'Venus'
+            test['outobs']  = 'ctphase_py5.fits'
+            test['logfile'] = 'ctphase_py5.log'
+            test.logFileOpen()
+            test.execute()
+            self.test_try_failure('Exception not thrown')
+        except (ValueError):
+            self.test_try_success()
+
+        # Test invalid model type
+        self.test_try('Test ctphase with invalid model')
+        try:
+            test = ctools.ctphase()
+            test['inobs']   = self._events
+            test['inmodel'] = self._model
+            test['srcname'] = 'Background'
+            test['outobs']  = 'ctphase_py6.fits'
+            test['logfile'] = 'ctphase_py6.log'
+            test.logFileOpen()
+            test.execute()
+            self.test_try_failure('Exception not thrown')
+        except (ValueError):
+            self.test_try_success()
+
+        # Test invalid model without temporal phase curve model
+        self.test_try('Test ctphase with invalid model')
+        try:
+            test = ctools.ctphase()
+            test['inobs']   = self._events
+            test['inmodel'] = self._model
+            test['srcname'] = 'Crab'
+            test['outobs']  = 'ctphase_py7.fits'
+            test['logfile'] = 'ctphase_py7.log'
+            test.logFileOpen()
+            test.execute()
+            self.test_try_failure('Exception not thrown')
+        except (ValueError):
+            self.test_try_success()
+
         # Return
         return
  
@@ -139,6 +139,7 @@ class Test(test):
         prob['caldb']   = self._caldb
         prob['irf']     = self._irf
         prob['edisp']   = False
+        prob['publish'] = True
         prob['logfile'] = 'ctprob_py1.log'
         prob['chatter'] = 2
  
@@ -189,6 +190,23 @@ class Test(test):
         # Check result file
         self._check_result_file('ctprob_py3_crab_events.fits')
  
+        # Copy ctprob tool
+        cpy_prob = prob.copy()
+
+        # Clear phase tool
+        prob.clear()
+
+        # Execute copy
+        cpy_prob['prefix']  = 'ctprob_py4_'
+        cpy_prob['outobs']  = 'ctprob_py4_obs.xml'
+        cpy_prob['logfile'] = 'ctprob_py4.log'
+        cpy_prob['chatter'] = 4
+        cpy_prob.logFileOpen()  # Needed to get a new log file
+        cpy_prob.execute()
+
+        # Check result file
+        self._check_result_file('ctprob_py4_ctprob_py3_crab_events.fits')
+
         # Return
         return