Remove phase curve parameter for ctphase (#1982)

Jürgen Knödlseder
1 parent 0eeff30e
Showing 6 changed files with 44 additions and 29 deletions
AUTHORS
ChangeLog
NEWS
src/ctphase/ctphase.cpp
src/ctphase/ctphase.par
test/test_ctphase.py
 Johan Bregeon (johan.bregeon@lupm.in2p3.fr)
 Rolf Buehler (rolf.buehler@desy.de)
+Joshua Cardenzana (jcardenzana@irap.omp.eu)
 Jean-Baptiste Cayrou (cayrou@etud.insa-toulouse.fr)
 Christoph Deil (Christoph.Deil@mpi-hd.mpg.de)
 Leonardo Di Venere (leonardo.divenere@ba.infn.it)
@@ -3,6 +3,7 @@
         * Version 1.3.0 released
         ========================
+        Add ctphase tool (#1982)
         Add ctprob tool (#1991)
         ctlike stores covariance matrix in either CSV or FITS format (#1863)
         Use correct ROOT histograms in csroot2caldb
@@ -142,6 +142,13 @@ event list a column for each model component that will contain the probability
 that the event originates from that component (#1991).
+ctphase - Event phase computation
+---------------------------------
+The ctphase tool has been added that compute the phase for each event based on
+a temporal phase curve model. This allows for example to perform a pulsar or
+binary analysis (#1982).
+
+
 cscripts
 --------
 The csebins script has been added that generates an energy boundaries
@@ -463,19 +463,23 @@ void ctphase::get_parameters(void)
     else {
         // Read phase curve user parameters
-        GFilename phasecurve = (*this)["phasecurve"].filename();
-        double    mjd_value  = (*this)["mjd"].real();
-        double    phase      = (*this)["phase"].real();
-        double    f0         = (*this)["f0"].real();
-        double    f1         = (*this)["f1"].real();
-        double    f2         = (*this)["f2"].real();
+        double mjd_value = (*this)["mjd"].real();
+        double phase     = (*this)["phase"].real();
+        double f0        = (*this)["f0"].real();
+        double f1        = (*this)["f1"].real();
+        double f2        = (*this)["f2"].real();
         // Set reference time
         GTime mjd;
         mjd.mjd(mjd_value);
         // Set phase curve
-        m_phase = GModelTemporalPhaseCurve(phasecurve, mjd, phase, f0, f1, f2);
+        m_phase = GModelTemporalPhaseCurve();
+        m_phase.mjd(mjd);
+        m_phase.phase(phase);
+        m_phase.f0(f0);
+        m_phase.f1(f1);
+        m_phase.f2(f2);
     } // endelse: read phase curve information from User parameters
@@ -25,12 +25,11 @@ source,  s, a, &quot;Crab&quot;,,, &quot;Source in model file to use for phase information&quot;
 #
 # Phase information parameters
 #=====================
-phasecurve, f, a, phase.fits,,, "Phase curve FITS file"
-mjd,        r, a, 51544.5,,, "Reference MJD for phase information"
-phase,      r, a, 0.0,0.0,1.0, "Phase at reference MJD"
-f0,         r, a, 1.0,,, "Frequency at reference MJD (Hz)"
-f1,         r, a, 0.1,,, "First freq. derivative at reference MJD (s^-1)"
-f2,         r, a, 0.01,,, "Second freq. derivative at reference MJD (s^-2)"
+mjd,   r, a, 51544.5,,, "Reference MJD for phase information"
+phase, r, a, 0.0,0.0,1.0, "Phase at reference MJD"
+f0,    r, a, 1.0,,, "Frequency at reference MJD (Hz)"
+f1,    r, a, 0.1,,, "First freq. derivative at reference MJD (s^-2)"
+f2,    r, a, 0.01,,, "Second freq. derivative at reference MJD (s^-3)"
 #
 # Standard parameters
@@ -41,7 +41,6 @@ class Test(test):
         # Set test data
         self._model_file     = self._datadir + '/model_temporal_phasecurve.xml'
-        self._phase_file     = self._datadir + '/model_temporal_phasecurve.fits'
         self._invalid_events = self._datadir + '/invalid_event_list.fits'
         # Return
@@ -153,17 +152,16 @@ class Test(test):
         self._check_result_file('ctphase_py1.fits')
         # Now run ctphase tool without a model file
-        phase['inmodel']    = 'NONE'
-        phase['source']     = 'NONE'
-        phase['phasecurve'] = self._phase_file
-        phase['mjd']        = 51544.5
-        phase['phase']      = 0.0
-        phase['f0']         = 1.0
-        phase['f1']         = 0.1
-        phase['f2']         = 0.01
-        phase['outobs']     = 'ctphase_py2.fits'
-        phase['logfile']    = 'ctphase_py2.log'
-        phase['chatter']    = 3
+        phase['inmodel'] = 'NONE'
+        phase['source']  = 'NONE'
+        phase['mjd']     = 51544.5
+        phase['phase']   = 0.0
+        phase['f0']      = 1.0
+        phase['f1']      = 0.1
+        phase['f2']      = 0.01
+        phase['outobs']  = 'ctphase_py2.fits'
+        phase['logfile'] = 'ctphase_py2.log'
+        phase['chatter'] = 3
         # Execute ctphase tool
         phase.logFileOpen()   # Make sure we get a log file
@@ -204,11 +202,16 @@ class Test(test):
         nevents : int, optional
             Expected number of events
         """
-        # Open result file
-        events      = gammalib.GCTAEventList(filename)
+        # Open event list
+        events = gammalib.GCTAEventList(filename)
+
+        # Check event list
+        self.test_value(events.size(), nevents, 'Check number of events')
+
+        # Check phases
         phases      = [events[i].phase() for i in range(events.size())]
-        test_result = int(sum(phases)>0)
-        self.test_value(test_result, 1, 'Check number of events')
+        test_result = int(sum(phases) > 0)
+        self.test_value(test_result, 1, 'Check whether phase information is set')
         # Return
         return