Add GCOMModelDRM class (#3910)

Jürgen Knödlseder
1 parent 784f3185
Showing 11 changed files with 820 additions and 8 deletions
ChangeLog
NEWS
README.md
doc/source/admin/release_history/2.0.rst
inst/com/Makefile.am
inst/com/include/GCOMModelDRM.hpp
inst/com/pyext/GCOMModelDRM.i
inst/com/pyext/com.i
inst/com/src/GCOMModelDRBFitting.cpp
inst/com/src/GCOMModelDRM.cpp
inst/com/test/test_COM.py
-2021-11-25
+2021-11-28
         * Version 2.0.0 released
         ========================
+        Add GCOMModelDRM class (#3910)
         Add angle(GVector&,GVector&) friend function (#3908)
         Make extension name check case insensitive in GFits::extno()
         Add numerical gradient step size member to GObservation
 New Features and Important Changes in GammaLib 2.0.0
-25 November 2021
+28 November 2021
 1. Introduction
@@ -22,6 +22,7 @@ The following classes have been added:
 - GModelSpectralBins
 - GCTAModelSkyCube
 - GCTABackground2D
+- GCOMModelDRM
 The following classes have been removed:
 - GCTAExceptionHandler
@@ -378,7 +379,9 @@ spatial model gradients (#3203).
 22. COMPTEL interface
 ---------------------
-Add phase selection to GCOMSelection and GCOMDri (#3564).
+Added GCOMModelDRM class to enable fitting of DRM data cubes to the data (#3910).
+
+Added phase selection to GCOMSelection and GCOMDri (#3564).
 Store DRB, DRG and DRX as instances of GCOMDri instead of GSkyMap in GCOMObservation
 (#3560).
 GammaLib information
 ====================
-* Version: 2.0.0.dev (25 November 2021)
+* Version: 2.0.0.dev (28 November 2021)
 [![Build Status](https://cta-jenkins.irap.omp.eu/buildStatus/icon?job=gammalib-integrate-os)](https://cta-jenkins.irap.omp.eu/job/gammalib-integrate-os/)
 .. _2.0:
-GammaLib 2.0 release (28 May 2021)
+GammaLib 2.0 release (28 November 2021)
 ==================================
 Overview
@@ -11,6 +11,7 @@ GammaLib 2.0 is a major release that adds significant functionality.
 In particular, this release provides:
 * vectorised response computation
+* full support of all spatial models for COMPTEL
 * ...
@@ -35,6 +36,10 @@ Bug fixes
 Improvements
 ------------
+* [`3908 <https://cta-redmine.irap.omp.eu/issues/3908>`_] -
+  Add ``angle(GVector&,GVector&)`` friend function
+* Make extension name check case insensitive in ``GFits::extno()``
+* Add numerical gradient step size member to ``GObservation``
 * [`3662 <https://cta-redmine.irap.omp.eu/issues/3662>`_] -
   Add ``stamp()`` methods to ``GApplication``
 * [`3630 <https://cta-redmine.irap.omp.eu/issues/3630>`_] -
@@ -83,6 +88,8 @@ Improvements
 New features
 ------------
+* [`3910 <https://cta-redmine.irap.omp.eu/issues/3910>`_] -
+  Add ``GCOMModelDRM`` class
 * [`3606 <https://cta-redmine.irap.omp.eu/issues/3606>`_] -
   Add ``GModelSpectralBins`` class
 * [`3527 <https://cta-redmine.irap.omp.eu/issues/3527>`_] -
@@ -55,7 +55,8 @@ sources = src/com_helpers_response_vector.cpp \
           src/GCOMInstChars.cpp \
           src/GCOMIaq.cpp \
           src/GCOMObservation.cpp \
-          src/GCOMModelDRBFitting.cpp
+          src/GCOMModelDRBFitting.cpp \
+          src/GCOMModelDRM.cpp
 # Define headers to be installed
 pkginclude_HEADERS = include/GCOMTools.hpp \
@@ -78,6 +79,7 @@ pkginclude_HEADERS = include/GCOMTools.hpp \
                      include/GCOMIaq.hpp \
                      include/GCOMObservation.hpp \
                      include/GCOMModelDRBFitting.hpp \
+                     include/GCOMModelDRM.hpp \
                      include/GCOMLib.hpp
 # Calibration files to distribute
+/***************************************************************************
+ *                GCOMModelDRM.hpp - COMPTEL DRM model class               *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2021 by 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   *
+ *  the Free Software Foundation, either version 3 of the License, or      *
+ *  (at your option) any later version.                                    *
+ *                                                                         *
+ *  This program is distributed in the hope that it will be useful,        *
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of         *
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          *
+ *  GNU General Public License for more details.                           *
+ *                                                                         *
+ *  You should have received a copy of the GNU General Public License      *
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.  *
+ *                                                                         *
+ ***************************************************************************/
+/**
+ * @file GCOMModelDRM.hpp
+ * @brief COMPTEL DRM model class interface definition
+ * @author Juergen Knoedlseder
+ */
+
+#ifndef GCOMMODELDRM_HPP
+#define GCOMMODELDRM_HPP
+
+/* __ Includes ___________________________________________________________ */
+#include <string>
+#include "GModelData.hpp"
+#include "GModelPar.hpp"
+#include "GEvent.hpp"
+#include "GObservation.hpp"
+#include "GNodeArray.hpp"
+#include "GXmlElement.hpp"
+#include "GCOMEventCube.hpp"
+
+
+/***********************************************************************//**
+ * @class GCOMModelDRM
+ *
+ * @brief COMPTEL DRM model class
+ *
+ * This class implements a COMPTEL data space model. The model has a
+ * normalization parameters.
+ ***************************************************************************/
+class GCOMModelDRM : public GModelData {
+
+public:
+    // Constructors and destructors
+    GCOMModelDRM(void);
+    explicit GCOMModelDRM(const GXmlElement& xml);
+    GCOMModelDRM(const GCOMModelDRM& model);
+    virtual ~GCOMModelDRM(void);
+
+    // Operators
+    virtual GCOMModelDRM& operator=(const GCOMModelDRM& model);
+
+    // Implemented pure virtual methods
+    virtual void           clear(void);
+    virtual GCOMModelDRM*  clone(void) const;
+    virtual std::string    classname(void) const;
+    virtual std::string    type(void) const;
+    virtual bool           is_constant(void) const;
+    virtual double         eval(const GEvent& event,
+                                const GObservation& obs,
+                                const bool& gradients = false) const;
+    virtual double         npred(const GEnergy& obsEng, const GTime& obsTime,
+                                 const GObservation& obs) const;
+    virtual GCOMEventCube* mc(const GObservation& obs, GRan& ran) const;
+    virtual void           read(const GXmlElement& xml);
+    virtual void           write(GXmlElement& xml) const;
+    virtual std::string    print(const GChatter& chatter = NORMAL) const;
+
+protected:
+    // Protected methods
+    void init_members(void);
+    void copy_members(const GCOMModelDRM& model);
+    void free_members(void);
+
+    // Proteced data members
+    GModelPar m_norm;     //!< Normalisation
+    GCOMDri   m_drm;      //!< Model
+    GFilename m_filename; //!< Name of DRM FITS file
+};
+
+
+/***********************************************************************//**
+ * @brief Return class name
+ *
+ * @return String containing the class name ("GCOMModelDRM").
+ ***************************************************************************/
+inline
+std::string GCOMModelDRM::classname(void) const
+{
+    return ("GCOMModelDRM");
+}
+
+
+/***********************************************************************//**
+ * @brief Return model type
+ *
+ * @return Model type.
+ *
+ * Returns the type of the model. The type for a DRM model is "DRM".
+ ***************************************************************************/
+inline
+std::string GCOMModelDRM::type(void) const
+{
+    return ("DRM");
+}
+
+
+/***********************************************************************//**
+ * @brief Signals if model is temporally constant
+ *
+ * @return True if model is temporally constant, false otherwise.
+ *
+ * Signals if the model is temporally constant. By definition, a DRM model
+ * is always temporally constant.
+ ***************************************************************************/
+inline
+bool GCOMModelDRM::is_constant(void) const
+{
+    return (true);
+}
+
+#endif /* GCOMMODELDRM_HPP */
+/***************************************************************************
+ *                GCOMModelDRM.i - COMPTEL DRM model class                 *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2021 by 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   *
+ *  the Free Software Foundation, either version 3 of the License, or      *
+ *  (at your option) any later version.                                    *
+ *                                                                         *
+ *  This program is distributed in the hope that it will be useful,        *
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of         *
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          *
+ *  GNU General Public License for more details.                           *
+ *                                                                         *
+ *  You should have received a copy of the GNU General Public License      *
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.  *
+ *                                                                         *
+ ***************************************************************************/
+/**
+ * @file GCOMModelDRM.hpp
+ * @brief COMPTEL DRM model class interface definition
+ * @author Juergen Knoedlseder
+ */
+%{
+/* Put headers and other declarations here that are needed for compilation */
+#include "GCOMModelDRM.hpp"
+%}
+
+
+/***********************************************************************//**
+ * @class GCOMModelDRM
+ *
+ * @brief COMPTEL DRM model class
+ ***************************************************************************/
+class GCOMModelDRM : public GModelData {
+public:
+    // Constructors and destructors
+    GCOMModelDRM(void);
+    explicit GCOMModelDRM(const GXmlElement& xml);
+    GCOMModelDRM(const GCOMModelDRM& model);
+    virtual ~GCOMModelDRM(void);
+
+    // Implemented pure virtual methods
+    virtual void           clear(void);
+    virtual GCOMModelDRM*  clone(void) const;
+    virtual std::string    classname(void) const;
+    virtual std::string    type(void) const;
+    virtual bool           is_constant(void) const;
+    virtual double         eval(const GEvent& event,
+                                const GObservation& obs,
+                                const bool& gradients = false) const;
+    virtual double         npred(const GEnergy& obsEng, const GTime& obsTime,
+                                 const GObservation& obs) const;
+    virtual GCOMEventCube* mc(const GObservation& obs, GRan& ran) const;
+    virtual void           read(const GXmlElement& xml);
+    virtual void           write(GXmlElement& xml) const;
+};
+
+
+/***********************************************************************//**
+ * @brief GCOMModelDRM class extension
+ ***************************************************************************/
+%extend GCOMModelDRM {
+    GCOMModelDRM copy() {
+        return (*self);
+    }
+%pythoncode {
+    def __getstate__(self):
+        xml = gammalib.GXmlElement()
+        self.write(xml)
+        state = {'xml': xml}
+        return state
+    def __setstate__(self, state):
+        if state['xml'].elements('parameter') == 0:
+            self.__init__()
+        else:
+            self.__init__(state['xml'])
+}
+};
@@ -114,6 +114,7 @@
 %include "GCOMInstChar.i"
 %include "GCOMIaq.i"
 %include "GCOMModelDRBFitting.i"
+%include "GCOMModelDRM.i"
 %include "GCOMTim.i"
 %include "GCOMOad.i"
 %include "GCOMOads.i"
@@ -44,6 +44,8 @@ const GModelRegistry      g_com_drb_fitting_registry(&amp;g_com_drb_fitting_seed);
 /* __ Method name definitions ____________________________________________ */
 #define G_EVAL     "GCOMModelDRBFitting::eval(GEvent&, GObservation&, bool&)"
+#define G_NRED  "GCOMModelDRBFitting::npred(GEnergy&, GTime&, GObservation&)"
+#define G_MC                  "GCOMModelDRBFitting::mc(GObservation&, GRan&)"
 #define G_READ                      "GCOMModelDRBFitting::read(GXmlElement&)"
 #define G_WRITE                    "GCOMModelDRBFitting::write(GXmlElement&)"
@@ -52,7 +54,6 @@ const GModelRegistry      g_com_drb_fitting_registry(&amp;g_com_drb_fitting_seed);
 /* __ Coding definitions _________________________________________________ */
 /* __ Debug definitions __________________________________________________ */
-//#define G_DUMP_MC                                  //!< Dump MC information
 /*==========================================================================
@@ -374,6 +375,10 @@ double GCOMModelDRBFitting::eval(const GEvent&amp;       event,
  * @param[in] obsEng Measured event energy.
  * @param[in] obsTime Measured event time.
  * @param[in] obs Observation.
+ * @return Spatially integrated data model.
+ *
+ * @exception GException::feature_not_implemented
+ *            Method is not implemented.
  *
  * Spatially integrates the data model for a given measured event energy and
  * event time.
@@ -387,6 +392,10 @@ double GCOMModelDRBFitting::npred(const GEnergy&amp;      obsEng,
     // Initialise result
     double npred = 0.0;
+    // Method is not implemented
+    std::string msg = "Spatial integration of data model is not implemented.";
+    throw GException::feature_not_implemented(G_NRED, msg);
+
     // Return
     return npred;
 }
@@ -397,6 +406,10 @@ double GCOMModelDRBFitting::npred(const GEnergy&amp;      obsEng,
  *
  * @param[in] obs Observation.
  * @param[in] ran Random number generator.
+ * @return COMPTEL event cube.
+ *
+ * @exception GException::feature_not_implemented
+ *            Method is not implemented.
  *
  * Draws a sample of events from the COMPTEL DRB fitting model using a Monte
  * Carlo simulation.
@@ -409,6 +422,10 @@ GCOMEventCube* GCOMModelDRBFitting::mc(const GObservation&amp; obs,
     // Initialise new event cube
     GCOMEventCube* cube = new GCOMEventCube;
+    // Method is not implemented
+    std::string msg = "Monte Carlo simulation of data model is not implemented.";
+    throw GException::feature_not_implemented(G_MC, msg);
+
     // Return
     return cube;
 }
+/***************************************************************************
+ *                GCOMModelDRM.cpp - COMPTEL DRM model class               *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2021 by 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   *
+ *  the Free Software Foundation, either version 3 of the License, or      *
+ *  (at your option) any later version.                                    *
+ *                                                                         *
+ *  This program is distributed in the hope that it will be useful,        *
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of         *
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          *
+ *  GNU General Public License for more details.                           *
+ *                                                                         *
+ *  You should have received a copy of the GNU General Public License      *
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.  *
+ *                                                                         *
+ ***************************************************************************/
+/**
+ * @file GCOMModelDRM.cpp
+ * @brief COMPTEL DRM model class implementation
+ * @author Juergen Knoedlseder
+ */
+
+/* __ Includes ___________________________________________________________ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+#include <typeinfo>
+#include "GException.hpp"
+#include "GTools.hpp"
+#include "GModelRegistry.hpp"
+#include "GCOMModelDRM.hpp"
+#include "GCOMObservation.hpp"
+#include "GCOMEventBin.hpp"
+
+/* __ Constants __________________________________________________________ */
+
+/* __ Globals ____________________________________________________________ */
+const GCOMModelDRM   g_com_drm_seed;
+const GModelRegistry g_com_drm_registry(&g_com_drm_seed);
+
+/* __ Method name definitions ____________________________________________ */
+#define G_EVAL            "GCOMModelDRM::eval(GEvent&, GObservation&, bool&)"
+#define G_NRED         "GCOMModelDRM::npred(GEnergy&, GTime&, GObservation&)"
+#define G_MC                         "GCOMModelDRM::mc(GObservation&, GRan&)"
+#define G_READ                             "GCOMModelDRM::read(GXmlElement&)"
+#define G_WRITE                           "GCOMModelDRM::write(GXmlElement&)"
+
+/* __ Macros _____________________________________________________________ */
+
+/* __ Coding definitions _________________________________________________ */
+
+/* __ Debug definitions __________________________________________________ */
+
+
+/*==========================================================================
+ =                                                                         =
+ =                        Constructors/destructors                         =
+ =                                                                         =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Void constructor
+ *
+ * Constructs an empty COMPTEL DRM model.
+ ***************************************************************************/
+GCOMModelDRM::GCOMModelDRM(void) : GModelData()
+{
+    // Initialise members
+    init_members();
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief XML constructor
+ *
+ * @param[in] xml XML element.
+ *
+ * Constructs a COMPTEL DRM model from the information that is found in an
+ * XML element. Please refer to the method GCOMModelDRM::read to learn more
+ * about the information that is expected in the XML element.
+ ***************************************************************************/
+GCOMModelDRM::GCOMModelDRM(const GXmlElement& xml) : GModelData(xml)
+{
+    // Initialise members
+    init_members();
+
+    // Read XML
+    read(xml);
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Copy constructor
+ *
+ * @param[in] model COMPTEL DRM model.
+ *
+ * Constructs a COMPTEL DRM model by copying information from an existing
+ * model. Note that the copy is a deep copy, so the original object can be
+ * destroyed after the copy without any loss of information.
+ ***************************************************************************/
+GCOMModelDRM::GCOMModelDRM(const GCOMModelDRM& model) : GModelData(model)
+{
+    // Initialise private members for clean destruction
+    init_members();
+
+    // Copy members
+    copy_members(model);
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Destructor
+ *
+ * Destroys a COMPTEL DRM model.
+ ***************************************************************************/
+GCOMModelDRM::~GCOMModelDRM(void)
+{
+    // Free members
+    free_members();
+
+    // Return
+    return;
+}
+
+
+/*==========================================================================
+ =                                                                         =
+ =                               Operators                                 =
+ =                                                                         =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Assignment operator
+ *
+ * @param[in] model COMPTEL DRM model.
+ * @return COMPTEL DRM model
+ *
+ * Assigns the information from a COMPTEL DRM model to the actual object.
+ * Note that a deep copy of the information is performed, so the original
+ * object can be destroyed after the assignment without any loss of
+ * information.
+ ***************************************************************************/
+GCOMModelDRM& GCOMModelDRM::operator=(const GCOMModelDRM& model)
+{
+    // Execute only if object is not identical
+    if (this != &model) {
+
+        // Copy base class members
+        this->GModelData::operator=(model);
+
+        // Free members
+        free_members();
+
+        // Initialise members
+        init_members();
+
+        // Copy members
+        copy_members(model);
+
+    } // endif: object was not identical
+
+    // Return
+    return *this;
+}
+
+
+/*==========================================================================
+ =                                                                         =
+ =                            Public methods                               =
+ =                                                                         =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Clear instance
+ *
+ * Resets the object to a clean initial state. All information that resided
+ * in the object will be lost.
+ ***************************************************************************/
+void GCOMModelDRM::clear(void)
+{
+    // Free class members (base and derived classes, derived class first)
+    free_members();
+    this->GModelData::free_members();
+    this->GModel::free_members();
+
+    // Initialise members
+    this->GModel::init_members();
+    this->GModelData::init_members();
+    init_members();
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Clone instance
+ *
+ * @return Pointer to deep copy of COMPTEL DRM model.
+ *
+ * Clone a COMPTEL DRM model. Cloning performs a deep copy of the
+ * information, so the original object can be destroyed after cloning without
+ * any loss of information.
+ ***************************************************************************/
+GCOMModelDRM* GCOMModelDRM::clone(void) const
+{
+    return new GCOMModelDRM(*this);
+}
+
+
+/***********************************************************************//**
+ * @brief Evaluate function
+ *
+ * @param[in] event Observed event.
+ * @param[in] obs Observation.
+ * @param[in] gradients Compute gradients?
+ * @return Background model value.
+ *
+ * @exception GException::invalid_argument
+ *            Observation is not a COMPTEL observation.
+ *            Event is not a COMPTEL event bin.
+ *
+ * Evaluates the COMPTEL DRM model.
+ ***************************************************************************/
+double GCOMModelDRM::eval(const GEvent&       event,
+                          const GObservation& obs,
+                          const bool&         gradients) const
+{
+    // Extract COMPTEL observation
+    const GCOMObservation* observation = dynamic_cast<const GCOMObservation*>(&obs);
+    if (observation == NULL) {
+        std::string cls = std::string(typeid(&obs).name());
+        std::string msg = "Observation of type \""+cls+"\" is not a COMPTEL "
+                          "observations. Please specify a COMPTEL observation "
+                          "as argument.";
+        throw GException::invalid_argument(G_EVAL, msg);
+    }
+
+    // Extract COMPTEL event bin
+    const GCOMEventBin* bin = dynamic_cast<const GCOMEventBin*>(&event);
+    if (bin == NULL) {
+        std::string cls = std::string(typeid(&event).name());
+        std::string msg = "Event of type \""+cls+"\" is  not a COMPTEL event. "
+                          "Please specify a COMPTEL event as argument.";
+        throw GException::invalid_argument(G_EVAL, msg);
+    }
+
+    // Initialise value
+    double value = 0.0;
+
+    // Optionally initialise gradients
+    if (gradients) {
+        m_norm.factor_gradient(0.0);
+    }
+
+    // Get bin index
+    int index = bin->index();
+
+    // Get bin size.
+    double size = bin->size();
+    
+    // Continue only if bin size is positive
+    if (size > 0.0) {
+
+        // Get DRM model value
+        value = m_drm.map().pixels()[index] / size;
+
+        // Optionally compute gradients
+        if (gradients) {
+
+            // Compute partial derivative
+            double grad = (m_norm.is_free()) ? value * m_norm.scale() : 0.0;
+
+            // Set gradient
+            m_norm.factor_gradient(grad);
+
+        } // endif: gradient computation was requested
+
+        // Compute model value
+        value *= m_norm.value();
+
+        // Compile option: Check for NaN/Inf
+        #if defined(G_NAN_CHECK)
+        if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
+            std::cout << "*** ERROR: GCOMModelDRM::eval";
+            std::cout << "(index=" << index << "):";
+            std::cout << " NaN/Inf encountered";
+            std::cout << " (value=" << value;
+            std::cout << ")" << std::endl;
+        }
+        #endif
+
+    } // endif: binsize was positive
+
+    // Return
+    return value;
+}
+
+
+/***********************************************************************//**
+ * @brief Return spatially integrated DRM model
+ *
+ * @param[in] obsEng Measured event energy.
+ * @param[in] obsTime Measured event time.
+ * @param[in] obs Observation.
+ * @return Spatially integrated DRM model.
+ *
+ * @exception GException::feature_not_implemented
+ *            Method is not implemented.
+ *
+ * Spatially integrates the DRM model for a given measured event energy and
+ * event time.
+ *
+ * @todo Implement method.
+ ***************************************************************************/
+double GCOMModelDRM::npred(const GEnergy&      obsEng,
+                           const GTime&        obsTime,
+                           const GObservation& obs) const
+{
+    // Initialise result
+    double npred = 0.0;
+
+    // Method is not implemented
+    std::string msg = "Spatial integration of DRM model is not implemented.";
+    throw GException::feature_not_implemented(G_NRED, msg);
+
+    // Return
+    return npred;
+}
+
+
+/***********************************************************************//**
+ * @brief Return simulated list of events
+ *
+ * @param[in] obs Observation.
+ * @param[in] ran Random number generator.
+ * @return COMPTEL event cube.
+ *
+ * @exception GException::feature_not_implemented
+ *            Method is not implemented.
+ *
+ * Draws a sample of events from the COMPTEL DRM model using a Monte Carlo
+ * simulation.
+ *
+ * @todo Implement method.
+ ***************************************************************************/
+GCOMEventCube* GCOMModelDRM::mc(const GObservation& obs, GRan& ran) const
+{
+    // Initialise new event cube
+    GCOMEventCube* cube = new GCOMEventCube;
+
+    // Method is not implemented
+    std::string msg = "Monte Carlo simulation of DRM model is not implemented.";
+    throw GException::feature_not_implemented(G_MC, msg);
+
+    // Return
+    return cube;
+}
+
+
+/***********************************************************************//**
+ * @brief Read model from XML element
+ *
+ * @param[in] xml XML element.
+ *
+ * Read the COMPTEL DRM model from an XML element.
+ *
+ * The model has the following XML file syntax:
+ *
+ *     <source name="Model" type="DRM" file="drm.fits" instrument="COM">
+ *       <parameter name="Normalization" ../>
+ *     </source>
+ ***************************************************************************/
+void GCOMModelDRM::read(const GXmlElement& xml)
+{
+    // Verify number of model parameters
+    gammalib::xml_check_parnum(G_READ, xml, 1);
+
+    // Get parameter pointers
+    const GXmlElement* norm = gammalib::xml_get_par(G_READ, xml, m_norm.name());
+
+    // Read parameters
+    m_norm.read(*norm);
+
+    // Read filename
+    m_filename = gammalib::xml_file_expand(xml, xml.attribute("file"));
+
+    // Load DRM
+    m_drm.load(m_filename);
+
+    // Read model attributes
+    read_attributes(xml);
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Write model into XML element
+ *
+ * @param[in] xml XML element.
+ *
+ * Write the COMPTEL DRM model into an XML element.
+ *
+ * The model has the following XML file syntax:
+ *
+ *     <source name="Model" type="DRM" file="drm.fits" instrument="COM">
+ *       <parameter name="Normalization" ../>
+ *     </source>
+ ***************************************************************************/
+void GCOMModelDRM::write(GXmlElement& xml) const
+{
+    // Initialise pointer on source
+    GXmlElement* src = NULL;
+
+    // Search corresponding source
+    int n = xml.elements("source");
+    for (int k = 0; k < n; ++k) {
+        GXmlElement* element = xml.element("source", k);
+        if (element->attribute("name") == name()) {
+            src = element;
+            break;
+        }
+    }
+
+    // If no source with corresponding name was found then append one.
+    // Set also the type and the instrument.
+    if (src == NULL) {
+        src = xml.append("source");
+        src->attribute("name", name());
+        src->attribute("type", type());
+        if (instruments().length() > 0) {
+            src->attribute("instrument", instruments());
+        }
+    }
+
+    // Verify model type
+    gammalib::xml_check_type(G_WRITE, *src, type());
+
+    // Get XML parameters
+    GXmlElement* norm = gammalib::xml_need_par(G_WRITE, *src, m_norm.name());
+
+    // Write parameters
+    m_norm.write(*norm);
+
+    // Set DRM file name
+    src->attribute("file", gammalib::xml_file_reduce(*src, m_filename));
+
+    // Write model attributes
+    write_attributes(*src);
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Print model information
+ *
+ * @param[in] chatter Chattiness.
+ * @return String containing model information.
+ ***************************************************************************/
+std::string GCOMModelDRM::print(const GChatter& chatter) const
+{
+    // Initialise result string
+    std::string result;
+
+    // Continue only if chatter is not silent
+    if (chatter != SILENT) {
+
+        // Append header
+        result.append("=== GCOMModelDRM ===");
+
+        // Append attributes
+        result.append("\n"+print_attributes());
+
+        // Append parameters
+        result.append("\n"+gammalib::parformat("DRM file"));
+        result.append(m_filename);
+        result.append("\n"+gammalib::parformat("Number of parameters"));
+        result.append(gammalib::str(size()));
+        for (int i = 0; i < size(); ++i) {
+            result.append("\n"+m_pars[i]->print(chatter));
+        }
+
+    } // endif: chatter was not silent
+
+    // Return result
+    return result;
+}
+
+
+/*==========================================================================
+ =                                                                         =
+ =                            Private methods                              =
+ =                                                                         =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Initialise class members
+ ***************************************************************************/
+void GCOMModelDRM::init_members(void)
+{
+    // Initialise Value
+    m_norm.clear();
+    m_norm.name("Normalization");
+    m_norm.value(1.0);
+    m_norm.scale(1.0);
+    m_norm.range(0.0, 1.0e6);
+    m_norm.gradient(0.0);
+    m_norm.free();
+    m_norm.has_grad(true);
+
+    // Set parameter pointer(s)
+    m_pars.clear();
+    m_pars.push_back(&m_norm);
+
+    // Initialise other members
+    m_drm.clear();
+    m_filename.clear();
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Copy class members
+ *
+ * @param[in] model Model.
+ ***************************************************************************/
+void GCOMModelDRM::copy_members(const GCOMModelDRM& model)
+{
+    // Copy members
+    m_norm     = model.m_norm;
+    m_drm      = model.m_drm;
+    m_filename = model.m_filename;
+
+    // Set parameter pointer(s)
+    m_pars.clear();
+    m_pars.push_back(&m_norm);
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Delete class members
+ ***************************************************************************/
+void GCOMModelDRM::free_members(void)
+{
+    // Return
+    return;
+}
 # ==========================================================================
 # This module performs unit tests for the GammaLib COMPTEL module.
 #
-# Copyright (C) 2012-2018 Juergen Knoedlseder
+# Copyright (C) 2012-2021 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
@@ -60,6 +60,7 @@ class Test(gammalib.GPythonTestSuite):
         #test_support.pickeling(self, gammalib.GCOMInstChars())
         test_support.pickeling(self, gammalib.GCOMInstDir())
         test_support.pickeling(self, gammalib.GCOMModelDRBFitting())
+        test_support.pickeling(self, gammalib.GCOMModelDRM())
         test_support.pickeling(self, gammalib.GCOMOad())
         test_support.pickeling(self, gammalib.GCOMOads())
         test_support.pickeling(self, gammalib.GCOMObservation())