diff --git a/src/ctphase/ctphase.cpp b/src/ctphase/ctphase.cpp index 34e1bd2..4d5da3b 100644 --- a/src/ctphase/ctphase.cpp +++ b/src/ctphase/ctphase.cpp @@ -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(const_cast(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; evntsize(); evnt++) { // Get the next event - GCTAEventAtom* event = (*list)[evnt]; + event = (*list)[evnt]; // Get the event time GTime time = event->time(); diff --git a/src/ctphase/ctphase.par b/src/ctphase/ctphase.par index c3bfd80..bf62251 100644 --- a/src/ctphase/ctphase.par +++ b/src/ctphase/ctphase.par @@ -24,11 +24,11 @@ prefix, s, h, "phased_",,, "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 diff --git a/test/test_ctphase.py b/test/test_ctphase.py index d02646c..2412ff4 100644 --- a/test/test_ctphase.py +++ b/test/test_ctphase.py @@ -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: