Add handling of D2 modules with failed PMT flag to comobsbin and support DRBPhib…

…arBins in comobsmodel script

Add handling of D2 modules with failed PMT flag to comobsbin and support DRBPhib…
…arBins in comobsmodel script
Jürgen Knödlseder
1 parent 79f16c8e
Showing 9 changed files with 116 additions and 34 deletions
ChangeLog
NEWS
README.md
doc/source/users/reference_manual/comobsbin.rst
doc/source/users/reference_manual/comobsmodel.rst
modules/comscripts/comobsbin.par
modules/comscripts/comobsbin.py
modules/comscripts/comobsmodel.par
modules/comscripts/comobsmodel.py
-2021-11-29
+2021-12-01
         * Version 2.0.0 released
         ========================
-        Add phi_first and phi_last parameters to comlixfit and comlixmap (#3912)
         Add csadd2caldb script
         Add restart of bisection in case of bad convergence in ctulimit
         Add COMPTEL analysis scripts (#3738)
 New Features and Important Changes in ctools 2.0.0
-29 November 2021
+1 December 2021
 Introduction
@@ -399,17 +399,11 @@ regions (#3248).
 comlixfit - COMPTEL maximum likelihood fitting using SRCLIX algorithm
 ---------------------------------------------------------------------
-Added phi_first and phi_last parameters to allow selection of the Phibar layers
-during the maximum likelihood fit (#3912).
-
 Added script (#3738).
 comlixmap - Generation of COMPTEL Test Statistic map using SRCLIX algorithm
 ---------------------------------------------------------------------------
-Added phi_first and phi_last parameters to allow selection of the Phibar layers
-during the maximum likelihood fit (#3912).
-
 Added script (#3738).
 ctools information
 ==================
-* Version:             2.0.0.dev (29 November 2021)
+* Version:             2.0.0.dev (1 December 2021)
 * GammaLib dependency: 2.0.0.dev
 [![Build Status](https://cta-jenkins.irap.omp.eu/buildStatus/icon?job=ctools-integrate-os)](https://cta-jenkins.irap.omp.eu/job/ctools-integrate-os/)
@@ -66,6 +66,18 @@ General parameters
 ``(nphibar = 25) [integer]``
     Number of bins in Phibar direction.
+``(psdmin = 0) [integer]``
+    Minimum PSD value.
+
+``(psdmax = 110) [integer]``
+    Maximum PSD value.
+
+``(zetamin = 5.0) [real]``
+    Minimum Earth horizon - Phibar (zeta) angle (deg).
+
+``(fpmtflag = 0) [integer]``
+    Handling of D2 modules with failed PMT flag (0: exclude, 1: include, 2: exclude PMT).
+
 Standard parameters
 -------------------
@@ -41,9 +41,12 @@ General parameters
 ``iso <NONE|MAP|CUBE> [string]``
     Isotropic component.
-``modeltype <NODES|BINS> [string]``
+``diffusetype <NODES|BINS> [string]``
     Diffuse model type.
+``bkgtype <NODES|BINS> [string]``
+    Background model type.
+
 ``(bremsmap = $COMDATA/../skymaps/galprop/map_comptel_bremsstrahlung.fits) [file]``
     Bremsstrahlung map file name.
@@ -39,7 +39,8 @@ npsi,     i, h, 80,1,500, &quot;Number of bins in Psi direction&quot;
 nphibar,  i, h, 25,1,100, "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"
+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)"
 #
 # Standard parameters
@@ -74,6 +74,7 @@ class comobsbin(ctools.csobservation):
         self['psdmin'].integer()
         self['psdmax'].integer()
         self['zetamin'].real()
+        self['fpmtflag'].integer()
         # Query ahead output model filename
         if self._read_ahead():
@@ -191,13 +192,14 @@ class comobsbin(ctools.csobservation):
             # Set phases
             self._select.phases(self._phases)
-            # Set suffix
+            # Set DRE, DRG and DRX suffix
             phases           = self['phase'].string()
-            self._drx_suffix = '_phases'+phases.replace(';', '_')
-            self._drg_suffix = '_phases'+phases.replace(';', '_')
             self._dre_suffix = '_phases'+phases.replace(';', '_')
+            self._drg_suffix = '_phases'+phases.replace(';', '_')
+            self._drx_suffix = '_phases'+phases.replace(';', '_')
         # If PSD interval differs from standard interval then set the interval
+        # and append flag to suffix
         if self['psdmin'].integer() != 0 or self['psdmax'].integer() != 110:
             # Set PSD interval
@@ -208,10 +210,22 @@ class comobsbin(ctools.csobservation):
             self._dre_suffix += '_psd%d-%d' % (self['psdmin'].integer(),
                                                self['psdmax'].integer())
-        # If zeta angle differs from standard value then append zeta angle to suffix
+        # If zeta angle differs from standard value then append zeta angle to DRE
+        # and DRG suffix
         if self['zetamin'].real() != 5.0:
-            self._drg_suffix += '_zeta%.1f' % (self['zetamin'].real())
             self._dre_suffix += '_zeta%.1f' % (self['zetamin'].real())
+            self._drg_suffix += '_zeta%.1f' % (self['zetamin'].real())
+
+        # If handling of D2 modules with failed PMT flag differs from standard then
+        # set the flag and append flag to suffix
+        if self['fpmtflag'].integer() != 0:
+
+            # Set handling of D2 modules with failed PMT flag
+            self._select.fpmtflag(self['fpmtflag'].integer())
+
+            # Set DRE and DRG suffix
+            self._dre_suffix += '_fpmt%1d' % (self['fpmtflag'].integer())
+            self._drg_suffix += '_fpmt%1d' % (self['fpmtflag'].integer())
         # Log input observations
         self._log_string(gammalib.NORMAL, str(self._select))
@@ -22,18 +22,19 @@ outmodel, f, a, models.xml,,, &quot;Output model definition file&quot;
 #
 # Script parameters
 #==================
-ra,        r, a, NONE,,, "Right Ascension of point source (deg)"
-dec,       r, a, NONE,,, "Declinaison of point source (deg)"
-srcname,   s, a, NONE,,, "Name of point source"
-brems,     s, a, NONE,NONE|MAP|CUBE,, "Bremsstrahlung component"
-ic,        s, a, NONE,NONE|MAP|CUBE,, "Inverse Compton component"
-iso,       s, a, NONE,NONE|CONST|CONSTFIX,, "Isotropic component"
-modeltype, s, h, BINS,NODES|BINS,, "Diffuse model type"
-ebinfile,  f, h, NONE,,, "Name of the file containing the diffuse model energy bin definition"
-bremsmap,  f, h, $COMDATA/../skymaps/galprop/map_comptel_bremsstrahlung.fits,,, "Bremsstrahlungs map"
-bremscube, f, h, $COMDATA/../skymaps/galprop/bremss_mapcube_54_77Xvarh7S.fits,,, "Bremsstrahlungs cube"
-icmap,     f, h, $COMDATA/../skymaps/galprop/map_comptel_ic.fits,,, "Inverse Compton map"
-iccube,    f, h, $COMDATA/../skymaps/galprop/ics_isotropic_mapcube_54_77Xvarh7S.fits,,, "Inverse Compton cube"
+ra,          r, a, NONE,,, "Right Ascension of point source (deg)"
+dec,         r, a, NONE,,, "Declinaison of point source (deg)"
+srcname,     s, a, NONE,,, "Name of point source"
+brems,       s, a, NONE,NONE|MAP|CUBE,, "Bremsstrahlung component"
+ic,          s, a, NONE,NONE|MAP|CUBE,, "Inverse Compton component"
+iso,         s, a, NONE,NONE|CONST|CONSTFIX,, "Isotropic component"
+diffusetype, s, h, BINS,NODES|BINS,, "Diffuse model type"
+bkgtype,     s, h, BINS,NODES|BINS,, "Background model type"
+ebinfile,    f, h, NONE,,, "Name of the file containing the diffuse model energy bin definition"
+bremsmap,    f, h, $COMDATA/../skymaps/galprop/map_comptel_bremsstrahlung.fits,,, "Bremsstrahlungs map"
+bremscube,   f, h, $COMDATA/../skymaps/galprop/bremss_mapcube_54_77Xvarh7S.fits,,, "Bremsstrahlungs cube"
+icmap,       f, h, $COMDATA/../skymaps/galprop/map_comptel_ic.fits,,, "Inverse Compton map"
+iccube,      f, h, $COMDATA/../skymaps/galprop/ics_isotropic_mapcube_54_77Xvarh7S.fits,,, "Inverse Compton cube"
 #
 # Standard parameters
@@ -64,7 +64,8 @@ class comobsmodel(ctools.csobservation):
         self['brems'].string()
         self['ic'].string()
         self['iso'].string()
-        self['modeltype'].string()
+        self['diffusetype'].string()
+        self['bkgtype'].string()
         if self['ebinfile'].is_valid():
             self['ebinfile'].filename()
         self['bremsmap'].filename()
@@ -111,9 +112,9 @@ class comobsmodel(ctools.csobservation):
         # Return energy boundaries
         return ebounds
-    def _generate_background_model(self, obs):
+    def _generate_background_nodes_model(self, obs):
         """
-        Generate background model component for observation
+        Generate nodes background model component for observation
         Parameters
         ----------
@@ -164,6 +165,60 @@ class comobsmodel(ctools.csobservation):
         # Return model
         return model
+    def _generate_background_bins_model(self, obs):
+        """
+        Generate bins background model component for observation
+
+        Parameters
+        ----------
+        obs : `~gammalib.GCOMObservation`
+            COMPTEL observation
+
+        Returns
+        -------
+        model : `~gammalib.GCOMModelDRBPhibarBins`
+            Background model
+        """
+        # Set model attributes
+        name = 'Background_%s' % obs.id()
+        id   = obs.id()
+
+        # Extract attributes from observation
+        drb     = obs.drb()
+        npix    = drb.nchi() * drb.npsi()
+        nphibar = drb.nphibar()
+        dphibar = 2.0         # Fixed so far as we have no method to recover
+
+        # Create XML model list object
+        xml = gammalib.GXml()
+        lib = xml.append('source_library title="source library"')
+        src = lib.append('source name="%s" type="DRBPhibarBins" instrument="COM" id="%s"' % (name, id))
+
+        # Loop over all phibar layers
+        for k in range(nphibar):
+
+            # Compute DRB content for phibar layer
+            sum_drb = 0.0
+            for i in range(npix):
+                index    = i + k * npix
+                sum_drb += drb[index]
+
+            # Free or fix parameter
+            if sum_drb > 0.0:
+                free = '1'
+            else:
+                free = '0'
+
+            # Add bin
+            src.append('parameter name="Normalization" value="1.00" '
+                       'scale="1" min="0" max="1000" free="%s"' % free)
+
+        # Create background model from XML source element
+        model = gammalib.GCOMModelDRBPhibarBins(src)
+
+        # Return model
+        return model
+
     def _generate_source_model(self):
         """
         Generate source model
@@ -231,7 +286,7 @@ class comobsmodel(ctools.csobservation):
             ebds = ebounds
         # Re-sample model for 'NODES' model type
-        if self['modeltype'].string() == 'NODES':
+        if self['diffusetype'].string() == 'NODES':
             # Setup energies vector for bin mean energies
             energies = gammalib.GEnergies()
@@ -485,7 +540,10 @@ class comobsmodel(ctools.csobservation):
             self._log_string(gammalib.NORMAL, str(obs))
             # Generate model
-            model = self._generate_background_model(obs)
+            if self['bkgtype'].string() == 'NODES':
+                model = self._generate_background_nodes_model(obs)
+            else:
+                model = self._generate_background_bins_model(obs)
             # Log model
             self._log_string(gammalib.EXPLICIT, str(model))