Adapt comobsbin script to new phase interface and add parameters for pulsar phase selection

Jürgen Knödlseder
1 parent 1aa478ab
Showing 3 changed files with 87 additions and 55 deletions
modules/comscripts/comobsbin.par
modules/comscripts/comobsbin.py
test/test_comobsbin.py
@@ -23,31 +23,33 @@ outfolder, s, a, $COMDATASTORE,,, &quot;Output folder for files&quot;
 #
 # Script parameters
 #==================
-ebinalg,  s, a, FILE,FILE|LIN|LOG,, "Algorithm for defining energy bins"
-emin,     r, a, 0.75,0.05,50.0, "Minimum energy (MeV)"
-emax,     r, a, 30.0,0.05,50.0, "Maximum energy (MeV)"
-phase,    s, a, NONE,,, "Phase expression in the format phasemin0-phasemax0;phasemin1-phasemax1;..."
-phase0,   t, a, JD 2451943.09,,, "Date of Phase 0 (UTC string, JD, MJD or MET in seconds)"
-period,   r, a, 3.90603,,, "Period (days)"
-enumbins, i, a, 1,1,20, "Number of energy bins"
-ebinfile, f, a, ebounds.fits,,, "Name of the file containing the energy bin definition"
-coordsys, s, a, GAL,CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)"
-proj,     s, a, TAN,AIT|AZP|CAR|GLS|MOL|SFL|SIN|STG|TAN,, "Projection method"
-usepnt,   b, h, yes,,, "Use pointing instead of chi0/psi0 parameters?"
-chi0,     r, a, 83.63,0,360, "Chi of DRI centre in degrees (RA or Galactic longitude)"
-psi0,     r, a, 22.51,-90,90, "Psi of DRI centre in degrees (DEC or Galactic latitude)"
-dchi,     r, h, 1.0,0.1,10.0, "Bin size in Chi direction (deg)"
-dpsi,     r, h, 1.0,0.1,10.0, "Bin size in Psi direction (deg)"
-dphibar,  r, h, 2.0,0.1,10.0, "Bin size in Phibar direction (deg)"
-nchi,     i, h, 80,1,800, "Number of bins in Chi direction"
-npsi,     i, h, 80,1,800, "Number of bins in Psi direction"
-nphibar,  i, h, 25,1,500, "Number of bins in Phibar direction"
-psdmin,   i, h, 0,0,110, "Minimum PSD value"
-psdmax,   i, h, 110,0,110, "Maximum PSD value"
-zetamin,  r, h, 5.0,0.0,10.0, "Minimum Earth horizon - Phibar (zeta) angle (deg)"
-fpmtflag, i, h, 0,0,2, "Handling of D2 modules with failed PMT flag (0: exclude, 1: include, 2: exclude PMT)"
-d1use,    s, h, 1111111,,, "D1 module usage (1: use, 0: don't use)"
-d2use,    s, h, 11111111111111,,, "D2 module usage (1: use, 0: don't use)"
+ebinalg,     s, a, FILE,FILE|LIN|LOG,, "Algorithm for defining energy bins"
+emin,        r, a, 0.75,0.05,50.0, "Minimum energy (MeV)"
+emax,        r, a, 30.0,0.05,50.0, "Maximum energy (MeV)"
+phase,       s, a, NONE,,, "Phase expression in the format phasemin0-phasemax0;phasemin1-phasemax1;..."
+psrname,     s, a, NONE,,, "Pulsar name"
+ephemerides, f, a, NONE,,, "Pulsar ephemerides file"
+phase0,      t, a, JD 2451943.09,,, "Date of Phase 0 (UTC string, JD, MJD or MET in seconds)"
+period,      r, a, 3.90603,,, "Period (days)"
+enumbins,    i, a, 1,1,20, "Number of energy bins"
+ebinfile,    f, a, ebounds.fits,,, "Name of the file containing the energy bin definition"
+coordsys,    s, a, CEL,CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)"
+proj,        s, a, TAN,AIT|AZP|CAR|GLS|MOL|SFL|SIN|STG|TAN,, "Projection method"
+usepnt,      b, h, yes,,, "Use pointing instead of chi0/psi0 parameters?"
+chi0,        r, a, 83.63,0,360, "Chi of DRI centre in degrees (RA or Galactic longitude)"
+psi0,        r, a, 22.51,-90,90, "Psi of DRI centre in degrees (DEC or Galactic latitude)"
+dchi,        r, h, 1.0,0.1,10.0, "Bin size in Chi direction (deg)"
+dpsi,        r, h, 1.0,0.1,10.0, "Bin size in Psi direction (deg)"
+dphibar,     r, h, 2.0,0.1,10.0, "Bin size in Phibar direction (deg)"
+nchi,        i, h, 80,1,800, "Number of bins in Chi direction"
+npsi,        i, h, 80,1,800, "Number of bins in Psi direction"
+nphibar,     i, h, 25,1,500, "Number of bins in Phibar direction"
+psdmin,      i, h, 0,0,110, "Minimum PSD value"
+psdmax,      i, h, 110,0,110, "Maximum PSD value"
+zetamin,     r, h, 5.0,0.0,10.0, "Minimum Earth horizon - Phibar (zeta) angle (deg)"
+fpmtflag,    i, h, 0,0,2, "Handling of D2 modules with failed PMT flag (0: exclude, 1: include, 2: exclude PMT)"
+d1use,       s, h, 1111111,,, "D1 module usage (1: use, 0: don't use)"
+d2use,       s, h, 11111111111111,,, "D2 module usage (1: use, 0: don't use)"
  
 #
 # Standard parameters
@@ -2,7 +2,7 @@
 # ==========================================================================
 # Bin COMPTEL observations
 #
-# Copyright (C) 2019-2021 Juergen Knoedlseder
+# Copyright (C) 2019-2022 Juergen Knoedlseder
 #
 # 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
@@ -62,8 +62,11 @@ class comobsbin(ctools.csobservation):
         self._get_ebounds()
         self._get_phases()
         if not self._phases.is_empty():
-            self['phase0'].time()
-            self['period'].real()
+            if gammalib.toupper(self['psrname'].string()) != 'NONE':
+                self['ephemerides'].filename()
+            else:
+                self['phase0'].time()
+                self['period'].real()
         self['coordsys'].string()
         self['proj'].string()
         if not self['usepnt'].boolean():
@@ -80,7 +83,7 @@ class comobsbin(ctools.csobservation):
         self['psdmax'].integer()
         self['zetamin'].real()
         self['fpmtflag'].integer()
-        
+
         # Get D1 and D2 module usage strings
         d1use = self['d1use'].string()
         d2use = self['d2use'].string()
@@ -183,39 +186,62 @@ class comobsbin(ctools.csobservation):
         # Clear selection set
         self._select.clear()
  
-        # If phases are set then set the phase curve and phases of the
-        # selection set
+        # If phases are set then handle phase selection
         if not self._phases.is_empty():
  
-            # Raise error if period is not positive
-            if self['period'].real() <= 0.0:
-                msg = 'Period %f is not positive. Please specify a positive period.'
-                raise RuntimeError(msg)
+            # If pulsar string is valid then handle pulsar
+            psrname = self['psrname'].string()
+            if gammalib.toupper(psrname) != 'NONE':
+
+                # Log pulsar phase selection
+                self._log_string(gammalib.NORMAL, 'Phase selection for pulsar "%s"' %
+                                 psrname)
+
+                # Read ephemerides
+                pulsar = gammalib.GPulsar(self['ephemerides'].filename(), psrname)
+
+                # Set pulsar ephemerides
+                self._select.pulsar(pulsar)
+
+                # Set pulsar phases
+                self._select.pulsar_phases(self._phases)
  
-            # Get frequency of variability in units of Hz
-            f0 = 1.0 / (self['period'].real() * gammalib.sec_in_day)
+                # Set phase suffix
+                suffix = '_pulphase'
+
+                # Set DRE, DRG and DRX suffix
+                phases           = self['phase'].string()
+                self._dre_suffix = suffix+phases.replace(';', '_')
+                self._drg_suffix = suffix
+                self._drx_suffix = suffix
+
+            # ... otherwise handle orbital phase selection
+            else:
+
+                # Raise error if period is not positive
+                period = self['period'].real()
+                if period <= 0.0:
+                    msg = 'Period %f is not positive. Please specify a positive period.'
+                    raise RuntimeError(msg)
  
-            # Set temporal phase curve
-            phase_curve = gammalib.GModelTemporalPhaseCurve()
+                # Log orbital phase selection
+                self._log_string(gammalib.NORMAL, 'Orbital selection for period %f days' %
+                                 period)
  
-            # Set phase curve elements
-            phase_curve.mjd(self['phase0'].time()) # Reference time
-            phase_curve.phase(0.0)                 # Phase @ reference time
-            phase_curve.f0(f0)                     # Frequency @ reference time
-            phase_curve.f1(0.0)                    # No frequency derivative
-            phase_curve.f2(0.0)                    # No 2nd frequency derivative
+                # Set orbital period
+                self._select.orbital_period(period, self['phase0'].time())
  
-            # Set phase curve
-            self._select.phase_curve(phase_curve)
+                # Set orbital phases
+                self._select.orbital_phases(self._phases)
  
-            # Set phases
-            self._select.phases(self._phases)
+                # Set phase suffix
+                suffix = '_orbphase'
  
-            # Set DRE, DRG and DRX suffix
-            phases           = self['phase'].string()
-            self._dre_suffix = '_phases'+phases.replace(';', '_')
-            self._drg_suffix = '_phases'+phases.replace(';', '_')
-            self._drx_suffix = '_phases'+phases.replace(';', '_')
+                # Set DRE, DRG and DRX suffix
+                phases           = self['phase'].string()
+                self._dre_suffix = suffix+phases.replace(';', '_')
+                self._drg_suffix = suffix+phases.replace(';', '_')
+                self._drx_suffix = suffix+phases.replace(';', '_')
  
         # If PSD interval differs from standard interval then set the interval
         # and append flag to suffix
@@ -2,7 +2,7 @@
 # ==========================================================================
 # This scripts performs unit tests for the comobsbin script.
 #
-# Copyright (C) 2021 Juergen Knoedlseder
+# Copyright (C) 2021-2022 Juergen Knoedlseder
 #
 # 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
@@ -81,6 +81,7 @@ class Test(test):
         # Setup command
         cmd = comobsbin+' inobs="'+self._obs+'" ebinalg="FILE" '+ \
                         'ebinfile="'+self._ebounds+'" phase="NONE" '+ \
+                        'coordsys="GAL" proj="TAN" '+ \
                         'outfolder="comobsbin_cmd1" outobs="comobsbin_cmd1.xml" '+ \
                         'logfile="comobsbin_cmd1.log" chatter=1'
  
@@ -94,6 +95,7 @@ class Test(test):
         # Setup command
         cmd = comobsbin+' inobs="file_that_does_not_exist.xml" ebinalg="FILE" '+ \
                         'ebinfile="'+self._ebounds+'" phase="NONE" '+ \
+                        'coordsys="GAL" proj="TAN" '+ \
                         'outfolder="comobsbin_cmd2" outobs="comobsbin_cmd2.xml" '+ \
                         'logfile="comobsbin_cmd2.log" chatter=1'
  
@@ -118,6 +120,8 @@ class Test(test):
         bin['ebinalg']   = 'FILE'
         bin['ebinfile']  = self._ebounds
         bin['phase']     = 'NONE'
+        bin['coordsys']  = 'GAL'
+        bin['proj']      = 'TAN'
         bin['outfolder'] = 'comobsbin_py1'
         bin['outobs']    = 'comobsbin_py1.xml'
         bin['logfile']   = 'comobsbin_py1.log'