Merge branch '1706-composite-models' into integration

# Conflicts: # AUTHORS # include/GTools.hpp # src/support/GTools.cpp

Merge branch '1706-composite-models' into integration
# Conflicts: # AUTHORS # include/GTools.hpp # src/support/GTools.cpp
Jürgen Knödlseder
2 parents 221a703e 89470207
Showing 26 changed files with 2259 additions and 11 deletions
AUTHORS
include/GModelSpatialComposite.hpp
include/GModelSpectralComposite.hpp
include/GTools.hpp
include/GTypemaps.hpp
include/GammaLib.hpp
include/Makefile.am
inst/cta/include/GCTAResponseCube.hpp
inst/cta/include/GCTAResponseIrf.hpp
inst/cta/src/GCTAResponseCube.cpp
inst/cta/src/GCTAResponseIrf.cpp
pyext/GModelSpatialComposite.i
pyext/GModelSpectralComposite.i
pyext/gammalib/model.i
src/model/GModelSky.cpp
src/model/GModelSpatialComposite.cpp
src/model/GModelSpatialEllipticalDisk.cpp
src/model/GModelSpatialPointSource.cpp
src/model/GModelSpectralComposite.cpp
src/model/Makefile.am
@@ -16,4 +16,5 @@ Michael Mayer (michael.mayer@physik.hu-berlin.de)
 Ellis Owen (Ellis.Owen@mpi-hd.mpg.de)
 Anneli Schulz (anneli.schulz@desy.de)
 Hubert Siejkowski (h.siejkowski@cyfronet.pl)
+Domenico Tiziani (domenico.tiziani@fau.de)
 Alexander Ziegler (alexander.ziegler@fau.de)
+/***************************************************************************
+ *     GModelSpatialPointSource.hpp - Spatial point source model class     *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2016 by Domenico Tiziani                                 *
+ * ----------------------------------------------------------------------- *
+ *                                                                         *
+ *  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 GModelSpatialComposite.hpp
+ * @brief Spatial composite model class interface definition
+ * @author Domenico Tiziani
+ */
+
+#ifndef GMODELSPATIALCOMPOSITE_HPP
+#define GMODELSPATIALCOMPOSITE_HPP
+
+/* __ Includes ___________________________________________________________ */
+#include <string>
+#include "GModelSpatial.hpp"
+#include "GModelPar.hpp"
+#include "GSkyDir.hpp"
+#include "GXmlElement.hpp"
+
+
+/***********************************************************************//**
+ * @class GModelSpatialComposite
+ *
+ * @brief Spatial composite model
+ *
+ * This class implements the spatial model for a composition of spatial
+ * models.
+ ***************************************************************************/
+class GModelSpatialComposite : public GModelSpatial {
+
+public:
+    // Constructors and destructors
+    GModelSpatialComposite(void);
+    GModelSpatialComposite(const bool& dummy, const std::string& type);
+    explicit GModelSpatialComposite(const GXmlElement& xml);
+    GModelSpatialComposite(const GModelSpatialComposite& model);
+    virtual ~GModelSpatialComposite(void);
+
+    // Operators
+    virtual GModelSpatialComposite& operator=(const GModelSpatialComposite& model);
+
+    // Implemented pure virtual methods
+    virtual void                    clear(void);
+    virtual GModelSpatialComposite* clone(void) const;
+    virtual std::string             classname(void) const;
+    virtual std::string             type(void) const;
+    virtual GClassCode              code(void) const;
+    virtual double                  eval(const GPhoton& photon,
+                                         const bool& gradients = false) const;
+    virtual GSkyDir                 mc(const GEnergy& energy,
+                                       const GTime& time,
+                                       GRan& ran) const;
+    virtual double                  mc_norm(const GSkyDir& dir,
+                                            const double&  radius) const;
+    virtual bool                    contains(const GSkyDir& dir,
+                                             const double&  margin = 0.0) const;
+    virtual void                    read(const GXmlElement& xml);
+    virtual void                    write(GXmlElement& xml) const;
+    virtual std::string             print(const GChatter& chatter = NORMAL) const;
+
+    // Other methods
+    int                  components(void) const;
+    void                 append(const GModelSpatial& component,
+                                const std::string&   name = "");
+    const GModelSpatial* component(const int& index) const;
+    const GModelSpatial* component(const std::string& name) const;
+
+
+protected:
+    // Protected methods
+    void init_members(void);
+    void copy_members(const GModelSpatialComposite& model);
+    void free_members(void);
+
+    // Protected members
+    std::string                 m_type;            //!< Model type
+    std::vector<GModelSpatial*> m_components;      //!< Components
+    std::vector<std::string>    m_component_names; //!< Component names
+};
+
+
+/***********************************************************************//**
+ * @brief Return class name
+ *
+ * @return String containing the class name ("GModelSpatialComposite").
+ ***************************************************************************/
+inline
+std::string GModelSpatialComposite::classname(void) const
+{
+    return ("GModelSpatialComposite");
+}
+
+
+/***********************************************************************//**
+ * @brief Return model type
+ *
+ * @return Model type.
+ *
+ * Returns the type of the spatial model.
+ ***************************************************************************/
+inline
+std::string GModelSpatialComposite::type(void) const
+{
+    return (m_type);
+}
+
+
+/***********************************************************************//**
+ * @brief Return class code
+ *
+ * @return GModelSpatialComposite.
+ *
+ * Returns the code GModelSpatialComposite of the class.
+ ***************************************************************************/
+inline
+GClassCode GModelSpatialComposite::code(void) const
+{
+    return GMODEL_SPATIAL_COMPOSITE;
+}
+
+
+/***********************************************************************//**
+ * @brief Return normalization of composite source for Monte Carlo
+ *        simulations
+ *
+ * @param[in] dir Centre of simulation cone.
+ * @param[in] radius Radius of simulation cone (degrees).
+ * @return Normalization.
+ *
+ * Returns unity.
+ ***************************************************************************/
+inline
+double GModelSpatialComposite::mc_norm(const GSkyDir& dir,
+                                       const double&  radius) const
+{
+    return 1.0;
+}
+
+
+/***********************************************************************//**
+ * @brief Return number of model components
+ *
+ * @return Number of model components.
+ ***************************************************************************/
+inline
+int GModelSpatialComposite::components(void) const
+{
+    return m_components.size();
+}
+
+#endif /* GMODELSPATIALCOMPOSITE_HPP */
+/***************************************************************************
+ *         GModelSpectralComposite.hpp - Spectral power law model class    *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2016 by Michael Mayer                               *
+ * ----------------------------------------------------------------------- *
+ *                                                                         *
+ *  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 GModelSpectralComposite.hpp
+ * @brief Composite spectral model class interface definition
+ * @author Michael Mayer
+ */
+
+#ifndef GMODELSPECTRALCOMPOSITE_HPP
+#define GMODELSPECTRALCOMPOSITE_HPP
+
+/* __ Includes ___________________________________________________________ */
+#include <string>
+#include "GModelSpectral.hpp"
+#include "GModelPar.hpp"
+#include "GEnergy.hpp"
+
+/* __ Forward declarations _______________________________________________ */
+class GRan;
+class GTime;
+class GXmlElement;
+
+
+/***********************************************************************//**
+ * @class GModelSpectralComposite
+ *
+ * @brief Composite spectral model class
+ *
+ * This class implements a composite spectrum. The model is defined by the sum
+ * of several individual spectral models. Each spectral model can be added via
+ * the XML interface or using the append() method.
+ *
+ ***************************************************************************/
+class GModelSpectralComposite : public GModelSpectral {
+
+public:
+    // Constructors and destructors
+	GModelSpectralComposite(void);
+    explicit GModelSpectralComposite(const GXmlElement& xml);
+    GModelSpectralComposite(const GModelSpectralComposite& model);
+    virtual ~GModelSpectralComposite(void);
+
+    // Operators
+    virtual GModelSpectralComposite& operator=(const GModelSpectralComposite& model);
+
+    // Implemented pure virtual methods
+    virtual void                     clear(void);
+    virtual GModelSpectralComposite* clone(void) const;
+    virtual std::string              classname(void) const;
+    virtual std::string              type(void) const;
+    virtual double                   eval(const GEnergy& srcEng,
+                                         const GTime&   srcTime = GTime(),
+                                         const bool&    gradients = false) const;
+    virtual double                   flux(const GEnergy& emin,
+                                         const GEnergy& emax) const;
+    virtual double                   eflux(const GEnergy& emin,
+                                         const GEnergy& emax) const;
+    virtual GEnergy                  mc(const GEnergy& emin,
+                                        const GEnergy& emax,
+                                        const GTime&   time,
+                                        GRan&          ran) const;
+    virtual void                     read(const GXmlElement& xml);
+    virtual void                     write(GXmlElement& xml) const;
+    virtual std::string              print(const GChatter& chatter = NORMAL) const;
+
+    // Other methods
+    void                  append(const GModelSpectral& spec, const std::string& name="");
+    int                   components(void) const;
+    const GModelSpectral* component(const int& index) const;
+    const GModelSpectral* component(const std::string& name) const;
+
+protected:
+    // Protected methods
+    void init_members(void);
+    void copy_members(const GModelSpectralComposite& model);
+    void free_members(void);
+    void update_mc_cache(const GEnergy& emin, const GEnergy& emax) const;
+
+    // Protected members
+    std::string                  m_type;        //!< Model type
+    std::vector<GModelSpectral*> m_spectral;    //!< Container of spectral models
+    std::vector<std::string>     m_components;  //!< Names of components
+
+    // MC cache
+    mutable double               m_mc_flux;     //!< Flux cache
+    mutable std::vector<double>  m_mc_probs;    //!< Probailities of individual components
+    mutable GEnergy              m_mc_emin;     //!< Last minimum energy
+    mutable GEnergy              m_mc_emax;     //!< Last maximum energy
+
+};
+
+
+/***********************************************************************//**
+ * @brief Return class name
+ *
+ * @return String containing the class name ("GModelSpectralComposite").
+ ***************************************************************************/
+inline
+std::string GModelSpectralComposite::classname(void) const
+{
+    return ("GModelSpectralComposite");
+}
+
+
+/***********************************************************************//**
+ * @brief Return model type
+ *
+ * @return Model type.
+ *
+ * Returns the type of the spectral composite model.
+ ***************************************************************************/
+inline
+std::string GModelSpectralComposite::type(void) const
+{
+    return (m_type);
+}
+
+
+/***********************************************************************//**
+ * @brief Return number of spectral components
+ *
+ * @return Number of model components.
+ *
+ * Returns the number of spectral components.
+ ***************************************************************************/
+inline
+int GModelSpectralComposite::components(void) const
+{
+    return (m_spectral.size());
+}
+
+
+#endif /* GMODELSPECTRALCOMPOSITE_HPP */
@@ -109,7 +109,7 @@ namespace gammalib {
     bool                     is_notanumber(const double& x);
     bool                     contains(const std::string& str,
                                       const std::string& substring);
-    bool                     contains(const std::vector<std::string>& strings,
+    bool                     contains(const std::vector<std::string> strings,
                                       const std::string& string);
     void                     warning(const std::string& origin,
                                      const std::string& message);
@@ -43,7 +43,8 @@ typedef enum {
     GMODEL_SPATIAL_POINT_SOURCE,
     GMODEL_SPATIAL_RADIAL,
     GMODEL_SPATIAL_ELLIPTICAL,
-    GMODEL_SPATIAL_DIFFUSE
+    GMODEL_SPATIAL_DIFFUSE,
+    GMODEL_SPATIAL_COMPOSITE
 } GClassCode;
  
 /* __ Typemaps ___________________________________________________________ */
@@ -208,6 +208,7 @@
 #include "GModelSpatialDiffuseConst.hpp"
 #include "GModelSpatialDiffuseCube.hpp"
 #include "GModelSpatialDiffuseMap.hpp"
+#include "GModelSpatialComposite.hpp"
 #include "GModelSpectral.hpp"
 #include "GModelSpectralRegistry.hpp"
 #include "GModelSpectralBrokenPlaw.hpp"
@@ -223,6 +224,7 @@
 #include "GModelSpectralPlaw.hpp"
 #include "GModelSpectralPlawPhotonFlux.hpp"
 #include "GModelSpectralPlawEnergyFlux.hpp"
+#include "GModelSpectralComposite.hpp"
 #include "GModelTemporal.hpp"
 #include "GModelTemporalRegistry.hpp"
 #include "GModelTemporalConst.hpp"
@@ -166,6 +166,7 @@ pkginclude_HEADERS = GammaLib.hpp \
                      GModelSpatialDiffuseConst.hpp \
                      GModelSpatialDiffuseCube.hpp \
                      GModelSpatialDiffuseMap.hpp \
+                     GModelSpatialComposite.hpp \
                      GModelSpectral.hpp \
                      GModelSpectralRegistry.hpp \
                      GModelSpectralBrokenPlaw.hpp \
@@ -181,6 +182,7 @@ pkginclude_HEADERS = GammaLib.hpp \
                      GModelSpectralPlaw.hpp \
                      GModelSpectralPlawPhotonFlux.hpp \
                      GModelSpectralPlawEnergyFlux.hpp \
+                     GModelSpectralComposite.hpp \
                      GModelTemporal.hpp \
                      GModelTemporalRegistry.hpp \
                      GModelTemporalConst.hpp
@@ -136,6 +136,9 @@ private:
     double irf_diffuse(const GEvent&       event,
                        const GSource&      source,
                        const GObservation& obs) const;
+    double irf_composite(const GEvent&       event,
+                         const GSource&      source,
+                         const GObservation& obs) const;
  
     // Private data members
     GCTACubeExposure   m_exposure;    //!< Exposure cube
@@ -176,6 +176,9 @@ private:
     double    irf_diffuse(const GEvent&       event,
                           const GSource&      source,
                           const GObservation& obs) const;
+    double    irf_composite(const GEvent&       event,
+                            const GSource&      source,
+                            const GObservation& obs) const;
     double    nroi_ptsrc(const GModelSky&    model,
                          const GEnergy&      srcEng,
                          const GTime&        srcTime,
@@ -200,6 +203,12 @@ private:
                            const GEnergy&      obsEng,
                            const GTime&        obsTime,
                            const GObservation& obs) const;
+    double    nroi_composite(const GModelSky&    model,
+                             const GEnergy&      srcEng,
+                             const GTime&        srcTime,
+                             const GEnergy&      obsEng,
+                             const GTime&        obsTime,
+                             const GObservation& obs) const;
  
     // Private data members
     GCaldb          m_caldb;          //!< Calibration database
@@ -40,6 +40,7 @@
 #include "GModelSpatialRadial.hpp"
 #include "GModelSpatialRadialShell.hpp"
 #include "GModelSpatialElliptical.hpp"
+#include "GModelSpatialComposite.hpp"
 #include "GModelSpatialDiffuse.hpp"
 #include "GPhoton.hpp"
 #include "GSource.hpp"
@@ -417,6 +418,9 @@ double GCTAResponseCube::irf(const GEvent&amp;       event,
         case GMODEL_SPATIAL_DIFFUSE:
             irf = irf_diffuse(event, source, obs);
             break;
+        case GMODEL_SPATIAL_COMPOSITE:
+            irf = irf_composite(event, source, obs);
+            break;
         default:
             break;
     }
@@ -1493,3 +1497,47 @@ double GCTAResponseCube::irf_diffuse(const GEvent&amp;       event,
     // Return IRF value
     return irf;
 }
+
+
+/***********************************************************************//**
+ * @brief Return instrument response to composite source
+ *
+ * @param[in] event Observed event.
+ * @param[in] source Source.
+ * @param[in] obs Observation.
+ * @return Instrument response to composite source.
+ *
+ * Returns the instrument response to a specified composite source.
+ ***************************************************************************/
+double GCTAResponseCube::irf_composite(const GEvent&       event,
+                                       const GSource&      source,
+                                       const GObservation& obs) const
+{
+    // Initialise IRF
+    double irf = 0.0;
+
+    // Get pointer to composite model
+    const GModelSpatialComposite* model =
+        dynamic_cast<const GModelSpatialComposite*>(source.model());
+
+    // Loop over model components
+    for (int i = 0; i < model->components(); ++i) {
+
+        // Get pointer to spatial component
+        GModelSpatial* spat = const_cast<GModelSpatial*>(model->component(i));
+
+        // Create new GSource object
+        GSource src(source.name(), spat, source.energy(), source.time());
+
+        // Compute irf value
+        irf += this->irf(event, source, obs);
+
+    }
+
+    // Divide by number of model components
+    if (model->components() > 0) {
+        irf /= double(model->components());
+    }
+
+    return irf;
+}
@@ -44,6 +44,7 @@
 #include "GModelSpatialRadial.hpp"
 #include "GModelSpatialRadialShell.hpp"
 #include "GModelSpatialElliptical.hpp"
+#include "GModelSpatialComposite.hpp"
 #include "GCTAObservation.hpp"
 #include "GCTAResponseIrf.hpp"
 #include "GCTAResponse_helpers.hpp"
@@ -431,6 +432,9 @@ double GCTAResponseIrf::irf(const GEvent&amp;       event,
         case GMODEL_SPATIAL_DIFFUSE:
             irf = irf_diffuse(event, source, obs);
             break;
+        case GMODEL_SPATIAL_COMPOSITE:
+            irf = irf_composite(event, source, obs);
+            break;
         default:
             break;
     }
@@ -1864,6 +1868,9 @@ double GCTAResponseIrf::nroi(const GModelSky&amp;    model,
         case GMODEL_SPATIAL_DIFFUSE:
             nroi = nroi_diffuse(model, srcEng, srcTime, obsEng, obsTime, obs);
             break;
+        case GMODEL_SPATIAL_COMPOSITE:
+            nroi = nroi_composite(model, srcEng, srcTime, obsEng, obsTime, obs);
+            break;
         default:
             break;
     }
@@ -2846,6 +2853,51 @@ double GCTAResponseIrf::irf_diffuse(const GEvent&amp;       event,
  
  
 /***********************************************************************//**
+ * @brief Return instrument response to composite source
+ *
+ * @param[in] event Observed event.
+ * @param[in] source Source.
+ * @param[in] obs Observation.
+ * @return Instrument response to composite source.
+ *
+ * Returns the instrument response to a specified composite source.
+ ***************************************************************************/
+double GCTAResponseIrf::irf_composite(const GEvent&       event,
+                                      const GSource&      source,
+                                      const GObservation& obs) const
+{
+    // Initialise IRF
+    double irf = 0.0;
+
+    // Get pointer to composite model
+    const GModelSpatialComposite* model =
+        dynamic_cast<const GModelSpatialComposite*>(source.model());
+
+    // Loop over model components
+    for (int i = 0; i < model->components(); ++i) {
+
+        // Get pointer to spatial component
+        GModelSpatial* spat = const_cast<GModelSpatial*>(model->component(i));
+
+        // Create new GSource object
+        GSource src(source.name(), spat, source.energy(), source.time());
+
+        // Compute irf value
+        irf += this->irf(event, src, obs);
+
+    }
+
+    // Divide by number of model components
+    if (model->components() > 0) {
+        irf /= double(model->components());
+    }
+
+    // Return IRF value
+    return irf;
+}
+
+
+/***********************************************************************//**
  * @brief Return spatial integral of point source model
  *
  * @param[in] model Sky Model.
@@ -3432,3 +3484,68 @@ double GCTAResponseIrf::nroi_diffuse(const GModelSky&amp;    model,
     // Return Nroi
     return nroi;
 }
+
+
+/***********************************************************************//**
+ * @brief Return spatial integral of composite source model
+ *
+ * @param[in] model Sky Model.
+ * @param[in] srcEng True photon energy.
+ * @param[in] srcTime True photon arrival time.
+ * @param[in] obsEng Observed event energy.
+ * @param[in] obsTime Observed event arrival time.
+ * @param[in] obs Observation.
+ *
+ * Computes the integral
+ *
+ * \f[
+ *    N_{\rm ROI}(E',t'|E,t) = \int_{\rm ROI} P(p',E',t'|E,t) dp'
+ * \f]
+ *
+ * of
+ *
+ * \f[
+ *    P(p',E',t'|E,t) = \int
+ *                      S(p,E,t) \times R(p',E',t'|p,E,t) \, dp
+ * \f]
+ *
+ * over the Region of Interest (ROI) for a composite source model
+ * \f$S(p,E,t)\f$ and the response function \f$R(p',E',t'|p,E,t)\f$.
+ ***************************************************************************/
+double GCTAResponseIrf::nroi_composite(const GModelSky&    model,
+                                       const GEnergy&      srcEng,
+                                       const GTime&        srcTime,
+                                       const GEnergy&      obsEng,
+                                       const GTime&        obsTime,
+                                       const GObservation& obs) const
+{
+    // Initialise nroi
+    double nroi = 0.0;
+
+    // Get composite model
+    GModelSpatialComposite* comp =
+        dynamic_cast<GModelSpatialComposite*>(model.spatial());
+
+    // Loop over model components
+    for (int i = 0; i < comp->components(); ++i) {
+
+        // Create new sky model
+        GModelSky sky = GModelSky(model);
+
+        // Assign spatial part
+        sky.spatial(comp->component(i));
+
+        // Compute nroi
+        nroi += this->nroi(sky, srcEng, srcTime, obsEng, obsTime, obs);
+
+    }
+
+    // Divide by number of model components
+    if (comp->components() > 0) {
+        nroi /= double(comp->components());
+    }
+
+    // Return nroi
+    return nroi;
+
+}
+/***************************************************************************
+ *        GModelSpatialComposite.i - Spatial point source model class      *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2016 by Domenico Tiziani                                 *
+ * ----------------------------------------------------------------------- *
+ *                                                                         *
+ *  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 GModelSpatialComposite.i
+ * @brief Spatial composite model class Python interface
+ * @author Domenico Tiziani
+ */
+%{
+/* Put headers and other declarations here that are needed for compilation */
+#include "GModelSpatialComposite.hpp"
+%}
+
+
+/***********************************************************************//**
+ * @class GModelSpatialComposite
+ *
+ * @brief Spatial composite model
+ ***************************************************************************/
+class GModelSpatialComposite  : public GModelSpatial {
+
+public:
+    // Constructors and destructors
+    GModelSpatialComposite(void);
+    explicit GModelSpatialComposite(const GXmlElement& xml);
+    GModelSpatialComposite(const GModelSpatialComposite& model);
+    virtual ~GModelSpatialComposite(void);
+
+    // Implemented virtual methods
+    virtual void                    clear(void);
+    virtual GModelSpatialComposite* clone(void) const;
+    virtual std::string             classname(void) const;
+    virtual std::string             type(void) const;
+    virtual double                  eval(const GPhoton& photon,
+                                         const bool& gradients = false) const;
+    virtual GSkyDir                 mc(const GEnergy& energy,
+                                       const GTime& time,
+                                       GRan& ran) const;
+    virtual double                  mc_norm(const GSkyDir& dir,
+                                            const double&  radius) const;
+    virtual bool                    contains(const GSkyDir& dir,
+                                             const double&  margin = 0.0) const;
+    virtual void                    read(const GXmlElement& xml);
+    virtual void                    write(GXmlElement& xml) const;
+
+    // Other methods
+    int                  components(void) const;
+    void                 append(const GModelSpatial& component,
+                                const std::string&   name = "");
+    const GModelSpatial* component(const int& index) const;
+    const GModelSpatial* component(const std::string& name) const;
+};
+
+
+/***********************************************************************//**
+ * @brief GModelSpatialComposite class extension
+ ***************************************************************************/
+%extend GModelSpatialComposite {
+    GModelSpatialComposite copy() {
+        return (*self);
+    }
+};
+/***************************************************************************
+ *          GModelSpectralComposite.i - Spectral power law model class     *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2016 Michael Mayer                                       *
+ * ----------------------------------------------------------------------- *
+ *                                                                         *
+ *  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 GModelSpecComposite spectral model class interface definition
+ * @author Michael Mayer
+ */
+%{
+/* Put headers and other declarations here that are needed for compilation */
+#include "GModelSpectralComposite.hpp"
+%}
+
+
+/***********************************************************************//**
+ * @class GModelSpectralComposite
+ *
+ * @brief Composite spectral model class
+ ***************************************************************************/
+class GModelSpectralComposite : public GModelSpectral {
+public:
+    // Constructors and destructors
+    GModelSpectralComposite(void);
+    explicit GModelSpectralComposite(const GXmlElement& xml);
+    GModelSpectralComposite(const GModelSpectralComposite& model);
+    virtual ~GModelSpectralComposite(void);
+
+    // Implemented pure virtual methods
+    virtual void                clear(void);
+    virtual GModelSpectralComposite* clone(void) const;
+    virtual std::string         classname(void) const;
+    virtual std::string         type(void) const;
+    virtual double              eval(const GEnergy& srcEng,
+                                     const GTime&   srcTime = GTime(),
+                                     const bool&    gradients = false) const;
+    virtual double              flux(const GEnergy& emin,
+                                     const GEnergy& emax) const;
+    virtual double              eflux(const GEnergy& emin,
+                                      const GEnergy& emax) const;
+    virtual GEnergy             mc(const GEnergy& emin,
+                                   const GEnergy& emax,
+                                   const GTime&   time,
+                                   GRan&          ran) const;
+    virtual void                read(const GXmlElement& xml);
+    virtual void                write(GXmlElement& xml) const;
+
+    // Other methods
+    void                append(const GModelSpectral& spec, const std::string& name="");
+    int                 components(void) const;
+    GModelSpectral*     component(const int& index) const;
+    GModelSpectral*     component(const std::string& name) const;
+
+};
+
+
+/***********************************************************************//**
+ * @brief GModelSpectralComposite class extension
+ ***************************************************************************/
+%extend GModelSpectralComposite {
+    GModelSpectralComposite copy() {
+        return (*self);
+    }
+};
@@ -141,6 +141,7 @@
 %include "GModelSpatialDiffuseConst.i"
 %include "GModelSpatialDiffuseCube.i"
 %include "GModelSpatialDiffuseMap.i"
+%include "GModelSpatialComposite.i"
 %include "GModelSpectral.i"
 %include "GModelSpectralRegistry.i"
 %include "GModelSpectralBrokenPlaw.i"
@@ -155,6 +156,7 @@
 %include "GModelSpectralPlaw.i"
 %include "GModelSpectralPlawPhotonFlux.i"
 %include "GModelSpectralPlawEnergyFlux.i"
+%include "GModelSpectralComposite.i"
 %include "GModelTemporal.i"
 %include "GModelTemporalRegistry.i"
 %include "GModelTemporalConst.i"
@@ -48,9 +48,11 @@
 const GModelSky         g_pointsource_seed("PointSource");
 const GModelSky         g_extendedsource_seed("ExtendedSource");
 const GModelSky         g_diffusesource_seed("DiffuseSource");
+const GModelSky         g_compositesource_seed("CompositeSource");
 const GModelRegistry    g_pointsource_registry(&g_pointsource_seed);
 const GModelRegistry    g_extendedsource_registry(&g_extendedsource_seed);
 const GModelRegistry    g_diffusesource_registry(&g_diffusesource_seed);
+const GModelRegistry    g_compositesource_registry(&g_compositesource_seed);
  
 /* __ Method name definitions ____________________________________________ */
 #define G_NPRED           "GModelSky::npred(GEnergy&, GTime&, GObservation&)"
+/***************************************************************************
+ *       GModelSpatialComposite.cpp - Spatial composite model class        *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2016 by Domenico Tiziani                                 *
+ * ----------------------------------------------------------------------- *
+ *                                                                         *
+ *  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 GModelSpatialComposite.cpp
+ * @brief Spatial composite model class implementation
+ * @author Domenico Tiziani
+ */
+
+/* __ Includes ___________________________________________________________ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+#include "GException.hpp"
+#include "GTools.hpp"
+#include "GMath.hpp"
+#include "GModelSpatialComposite.hpp"
+#include "GModelSpatialRegistry.hpp"
+
+/* __ Constants __________________________________________________________ */
+
+
+/* __ Globals ____________________________________________________________ */
+const GModelSpatialComposite g_spatial_comp_seed;
+const GModelSpatialRegistry  g_spatial_comp_registry(&g_spatial_comp_seed);
+
+/* __ Method name definitions ____________________________________________ */
+#define G_READ                   "GModelSpatialComposite::read(GXmlElement&)"
+#define G_WRITE                 "GModelSpatialComposite::write(GXmlElement&)"
+#define G_COMPONENT_INDEX           "GModelSpatialComposite::component(int&)"
+#define G_COMPONENT_NAME    "GModelSpatialComposite::component(std::string&)"
+#define G_APPEND            "GModelSpatialComposite::append(GModelSpatial&, "\
+                                                              "std::string&)"
+
+/* __ Macros _____________________________________________________________ */
+
+/* __ Coding definitions _________________________________________________ */
+
+/* __ Debug definitions __________________________________________________ */
+
+
+/*==========================================================================
+ =                                                                         =
+ =                        Constructors/destructors                         =
+ =                                                                         =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Void constructor
+ *
+ * Constructs empty spatial composite model.
+ ***************************************************************************/
+GModelSpatialComposite::GModelSpatialComposite(void) : GModelSpatial()
+{
+    // Initialise members
+    init_members();
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Model type constructor
+ *
+ * @param[in] dummy Dummy flag.
+ * @param[in] type Model type.
+ *
+ * Constructs empty spatial composite model by specifying a model @p type.
+ ***************************************************************************/
+GModelSpatialComposite::GModelSpatialComposite(const bool&        dummy,
+                                               const std::string& type) :
+                        GModelSpatial()
+{
+    // Initialise members
+    init_members();
+
+    // Set model type
+    m_type = type;
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief XML constructor
+ *
+ * @param[in] xml XML element.
+ *
+ * Construct a spatial composite model by extracting information from an
+ * XML element. See the read() method for more information about the expected
+ * structure of the XML element.
+ ***************************************************************************/
+GModelSpatialComposite::GModelSpatialComposite(const GXmlElement& xml) :
+                        GModelSpatial()
+{
+    // Initialise members
+    init_members();
+    
+    // Read information from XML element
+    read(xml);
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Copy constructor
+ *
+ * @param[in] model Spatial composite model.
+ ***************************************************************************/
+GModelSpatialComposite::GModelSpatialComposite(const GModelSpatialComposite& model) :
+                        GModelSpatial(model)
+{
+    // Initialise members
+    init_members();
+
+    // Copy members
+    copy_members(model);
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Destructor
+ ***************************************************************************/
+GModelSpatialComposite::~GModelSpatialComposite(void)
+{
+    // Free members
+    free_members();
+
+    // Return
+    return;
+}
+
+
+/*==========================================================================
+ =                                                                         =
+ =                               Operators                                 =
+ =                                                                         =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Assignment operator
+ *
+ * @param[in] model Spatial composite model.
+ * @return Spatial composite model.
+ ***************************************************************************/
+GModelSpatialComposite& GModelSpatialComposite::operator=(const GModelSpatialComposite& model)
+{
+    // Execute only if object is not identical
+    if (this != &model) {
+
+        // Copy base class members
+        this->GModelSpatial::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 spatial composite model
+ ***************************************************************************/
+void GModelSpatialComposite::clear(void)
+{
+    // Free class members (base and derived classes, derived class first)
+    free_members();
+    this->GModelSpatial::free_members();
+
+    // Initialise members
+    this->GModelSpatial::init_members();
+    init_members();
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Clone spatial composite model
+ *
+ * @return Pointer to deep copy of spatial composite model.
+ ***************************************************************************/
+GModelSpatialComposite* GModelSpatialComposite::clone(void) const
+{
+    // Clone point source model
+    return new GModelSpatialComposite(*this);
+}
+
+
+/***********************************************************************//**
+ * @brief Evaluate function
+ *
+ * @param[in] photon Incident photon.
+ * @param[in] gradients Compute gradients?
+ * @return Model value.
+ *
+ * Evaluates the spatial composite model by summing all model components and
+ * dividing the sum by the number of components.
+ ***************************************************************************/
+double GModelSpatialComposite::eval(const GPhoton& photon,
+                                    const bool&    gradients) const
+{
+    // Initalise value
+    double value = 0.0;
+    
+    // Sum over all components
+    for (int i = 0; i < m_components.size(); ++i) {
+        value += m_components[i]->eval(photon, gradients);
+    }
+    
+    // Normalise sum by the number of components
+    if (m_components.size() > 0) {
+        value /= double(m_components.size());
+    }
+
+    // Return value
+    return value;
+}
+
+
+/***********************************************************************//**
+ * @brief Returns MC sky direction
+ *
+ * @param[in] energy Photon energy.
+ * @param[in] time Photon arrival time.
+ * @param[in,out] ran Random number generator.
+ * @return Sky direction.
+ *
+ * Draws an arbitrary model component and an arbitrary direction from that 
+ * component.
+ ***************************************************************************/
+GSkyDir GModelSpatialComposite::mc(const GEnergy& energy,
+                                   const GTime&   time,
+                                   GRan&          ran) const
+{
+    // Randomly choose one model component
+    int index = int(ran.uniform() * double(m_components.size()));
+    if (index >= m_components.size()) {
+        index = m_components.size();
+    }
+
+    // Draw random sky direction from the selected component
+    GSkyDir sky_dir = m_components[index]->mc(energy, time, ran);
+    
+    // Return sky direction
+    return (sky_dir);
+}
+
+
+/***********************************************************************//**
+ * @brief Checks where model contains specified sky direction
+ *
+ * @param[in] dir Sky direction.
+ * @param[in] margin Margin to be added to sky direction (degrees)
+ * @return True if the model contains the sky direction.
+ *
+ * Signals whether a sky direction is contained in one of the 
+ * model components.
+ ***************************************************************************/
+bool GModelSpatialComposite::contains(const GSkyDir& dir,
+                                      const double&  margin) const
+{
+    // Initialise containment flag
+    bool containment = false;
+
+    // Loop over all components
+    for (int i = 0; i < m_components.size(); ++i) {
+        if (m_components[i]->contains(dir, margin)) {
+            containment = true;
+            break;
+        }
+    }
+
+    // Return containment
+    return containment;
+}
+
+
+/***********************************************************************//**
+ * @brief Read model from XML element
+ *
+ * @param[in] xml XML element.
+ *
+ * Reads the spatial information from an XML element.
+ ***************************************************************************/
+void GModelSpatialComposite::read(const GXmlElement& xml)
+{   
+    // Get number of spatial components
+    int n_comp = xml.elements("spatialModel");
+    
+    // Loop over spatial elements
+    for (int i = 0; i < n_comp; ++i) {
+        
+        // Get spatial XML element
+        const GXmlElement* spec = xml.element("spatialModel", i);
+        
+        // Initialise a spatial registry object
+        GModelSpatialRegistry registry;
+        
+        // Read spatial model
+        GModelSpatial* ptr = registry.alloc(*spec);
+        
+        // Get component attribute from XML file
+        std::string component_name = spec->attribute("component");
+        
+        // Append spatial component to container
+        append(*ptr, component_name);
+        
+    } // endfor: loop over components
+    
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Write model into XML element
+ *
+ * @param[in] xml XML element into which model information is written.
+ *
+ * @exception GException::model_invalid_spatial
+ *            Existing XML element is not of type 'SkyDirFunction'
+ * @exception GException::model_invalid_parnum
+ *            Invalid number of model parameters found in XML element.
+ * @exception GException::model_invalid_parnames
+ *            Invalid model parameter names found in XML element.
+ *
+ * Write the spatial information into an XML element.
+ ***************************************************************************/
+void GModelSpatialComposite::write(GXmlElement& xml) const
+{
+    // Set model type
+    if (xml.attribute("type") == "") {
+        xml.attribute("type", type());
+    }
+
+    // Verify model type
+    if (xml.attribute("type") != type()) {
+        throw GException::model_invalid_spatial(G_WRITE, xml.attribute("type"),
+              "Spatial model is not of type \""+type()+"\".");
+    }
+
+    // Write all model components
+    for (int i = 0; i < m_components.size(); ++i) {
+
+        // Get spatial component
+        GModelSpatial* component = m_components[i];
+        
+        // Fall through if component is empty
+        if (component == NULL) {
+            continue;
+        }
+
+        // Find XML element with matching name
+        GXmlElement *matching_model = NULL;
+        for (int k = 0; k < xml.elements("spatialModel"); ++k) {
+            if (xml.element("spatialModel", k)->attribute("name") ==
+                m_component_names[i]) {
+                matching_model = xml.element("spatialModel", k);
+                break;
+            }
+        } // endfor: loop over all XML elements
+
+        // Create copy of the spatial model
+        GModelSpatial* spatial = component->clone();
+
+        // Loop over all parameters of component and strip prefix
+        for (int k = 0; k < spatial->size(); ++k) {
+            
+            // Get model parameter
+            GModelPar& par = (*spatial)[k];
+
+            // If name contains colon then strip the prefix and the colon from
+            // the name
+            std::string name  = par.name();
+            int         found = name.find(":");
+            if (found != std::string::npos) {
+                name = name.substr(found+1, std::string::npos);
+                par.name(name);
+            }
+
+        } // endfor looped over all parameters
+
+        // If an XML element exists then use it ...
+        if (matching_model != NULL) {
+            spatial->write(*matching_model);
+        }
+
+        // ... otherwise create new XML element
+        else {
+            GXmlElement element("spatialModel");
+            spatial->write(element);
+            xml.append(element);
+        }
+
+        // Delete clone
+        delete spatial;
+
+    } // endfor: loop over all components
+
+    // TODO: Remove unused xml elements
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Append spatial component
+ *
+ * @param[in] component Spatial model component to append
+ * @param[in] name Name of spatial model (can be empty)
+ * 
+ * Appends a spatial component to the composite model
+ ***************************************************************************/
+void GModelSpatialComposite::append(const GModelSpatial& component,
+                                    const std::string&   name)
+{
+    // Append model container
+    m_components.push_back(component.clone());
+
+    // Get index of latest model
+    int index = m_components.size()-1;
+    
+    // Use model index if component name is empty
+    std::string component_name = !name.empty() ? name
+                                               : gammalib::str(m_components.size());
+    
+    // Check if component name is unique, throw exception if not
+    if (gammalib::contains(m_component_names, component_name)) {
+    	std::string msg = "Attempt to append component \""+component_name+"\" "
+                          "to composite spatial model, but a component with the "
+                          "same name exists already. Each component needs a "
+                          "unique name.";
+        throw GException::invalid_value(G_APPEND, msg);
+    }
+
+    // Add component name
+    m_component_names.push_back(component_name);
+    
+    // Get number of spectral parameters from model
+    int npars = m_components[index]->size();
+    
+    // Loop over model parameters
+    for (int ipar = 0; ipar < npars; ++ipar) {
+        
+        // Get model parameter
+        GModelPar* par = &(m_components[index]->operator[](ipar));
+        
+        // Modify parameter name
+        par->name(component_name+":"+par->name());
+        
+        // Append model parameter with new name to internal container
+        m_pars.push_back(par);
+        
+    } // endfor: loop over model parameters
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Returns pointer to spatial component element
+ *
+ * @param[in] index Index of spatial component [0,...,components()-1].
+ * @return Spatial model.
+ *
+ * @exception GException::out_of_range
+ *            Index is out of range.
+ *
+ * Returns a spatial component to the composite model
+ ***************************************************************************/
+const GModelSpatial* GModelSpatialComposite::component(const int& index) const
+{
+	// Check if index is in validity range
+	if (index < 0 || index >= m_components.size()) {
+		throw GException::out_of_range(G_COMPONENT_INDEX, "Component Index",
+                                       index, m_components.size());
+	}
+
+	// Return spatial component
+	return m_components[index];
+}
+
+
+/***********************************************************************//**
+ * @brief Returns pointer to specific spatial component
+ *
+ * @param[in] name Name of spatial component.
+ * @return Spatial model.
+ *
+ * @exception GException::invalid_argument
+ *            Spatial model not found.
+ *
+ * Returns a spatial component of the composite model
+ ***************************************************************************/
+const GModelSpatial* GModelSpatialComposite::component(const std::string& name) const
+{
+	// Check if model name is found
+	int index = -1;
+	for (int i = 0; i < m_component_names.size(); ++i) {
+        if (m_component_names[i] == name) {
+			index = i;
+			break;
+		}
+	}
+
+	// Check if component name was found
+	if (index == -1) {
+        std::string msg = "Model component \""+name+"\" not found in composite "
+                          "spatial model.";
+		throw GException::invalid_argument(G_COMPONENT_NAME, msg);
+	}
+
+	// Return spatial component
+	return m_components[index];
+}
+
+
+/***********************************************************************//**
+ * @brief Print composite spatial model information
+ *
+ * @param[in] chatter Chattiness.
+ * @return String containing composite spatial model information.
+ ***************************************************************************/
+std::string GModelSpatialComposite::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("=== GModelSpatialComposite ===");
+
+        // Append information
+		result.append("\n"+gammalib::parformat("Number of components"));
+		result.append(gammalib::str(components()));
+		result.append("\n"+gammalib::parformat("Number of parameters"));
+		result.append(gammalib::str(size()));
+
+		// Print parameter information
+		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 GModelSpatialComposite::init_members(void)
+{
+    // Initialise model type
+    m_type = "Composite";
+
+    // Initialise models vector
+    m_components.clear();
+    m_component_names.clear();
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Copy class members
+ *
+ * @param[in] model Spatial composite model.
+ ***************************************************************************/
+void GModelSpatialComposite::copy_members(const GModelSpatialComposite& model)
+{
+    // Copy members
+    m_type            = model.m_type;
+    m_component_names = model.m_component_names;
+
+    // Initialise components
+    m_components.clear();
+    m_pars.clear();
+
+    // Copy components
+    for (int i = 0; i < model.m_components.size(); ++i) {
+
+        // Clone component
+        m_components.push_back(model.m_components[i]->clone());
+
+        // Get number of parameters to append
+        int npars = m_components[i]->size();
+    
+        // Append parameter references
+        for (int ipar = 0; ipar < npars; ++ipar) {
+            GModelPar& par = (*m_components[i])[ipar];
+            m_pars.push_back(&par);
+        }
+
+    } // endfor: looped over all components
+    
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Delete class members
+ ***************************************************************************/
+void GModelSpatialComposite::free_members(void)
+{
+    // Free model components
+    for (int i = 0; i < m_components.size(); ++i) {
+        
+        // Delete component i
+        delete m_components[i];
+
+        // Signal free pointer
+        m_components[i] = NULL;
+
+    }
+    
+    // Return
+    return;
+}
@@ -599,7 +599,7 @@ void GModelSpatialEllipticalDisk::write(GXmlElement&amp; xml) const
 /***********************************************************************//**
  * @brief Print information
  *
- * @param[in] chatter Chattiness (defaults to NORMAL).
+ * @param[in] chatter Chattiness.
  * @return String containing model information.
  ***************************************************************************/
 std::string GModelSpatialEllipticalDisk::print(const GChatter& chatter) const
@@ -501,6 +501,7 @@ void GModelSpatialPointSource::write(GXmlElement&amp; xml) const
 /***********************************************************************//**
  * @brief Print point source information
  *
+ * @param[in] chatter Chattiness.
  * @return String containing point source information.
  ***************************************************************************/
 std::string GModelSpatialPointSource::print(const GChatter& chatter) const
+/***************************************************************************
+ *         GModelSpectralComposite.cpp - Spectral composite model class    *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2016 by Michael Mayer                                    *
+ * ----------------------------------------------------------------------- *
+ *                                                                         *
+ *  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 GModelSpectralComposite.cpp
+ * @brief Composite spectral model class implementation
+ * @author Michael Mayer
+ */
+
+/* __ Includes ___________________________________________________________ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+#include <cmath>
+#include "GException.hpp"
+#include "GTools.hpp"
+#include "GRan.hpp"
+#include "GModelSpectralComposite.hpp"
+#include "GModelSpectralRegistry.hpp"
+
+/* __ Constants __________________________________________________________ */
+
+/* __ Globals ____________________________________________________________ */
+const GModelSpectralComposite     g_spectral_comp_seed;
+const GModelSpectralRegistry g_spectral_comp_registry(&g_spectral_comp_seed);
+
+/* __ Method name definitions ____________________________________________ */
+#define G_MC               "GModelSpectralComposite::mc(GEnergy&, GEnergy&, GTime&, GRan&)"
+#define G_WRITE            "GModelSpectralComposite::write(GXmlElement&)"
+#define G_COMPONENT_INDEX  "GModelSpectralComposite::component(const int&)"
+#define G_COMPONENT_NAME   "GModelSpectralComposite::component(const std::string&)"
+#define G_APPEND           "GModelSpectralComposite::append(const GModelSpectral&, const std::string&)"
+
+/* __ Macros _____________________________________________________________ */
+
+/* __ Coding definitions _________________________________________________ */
+
+/* __ Debug definitions __________________________________________________ */
+
+
+/*==========================================================================
+ =                                                                         =
+ =                        Constructors/destructors                         =
+ =                                                                         =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Void constructor
+ ***************************************************************************/
+GModelSpectralComposite::GModelSpectralComposite(void) : GModelSpectral()
+{
+    // Initialise private members for clean destruction
+    init_members();
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief XML constructor
+ *
+ * @param[in] xml XML element containing position information.
+ *
+ * Constructs a power law spectral model by extracting information from an
+ * XML element. See the read() method for more information about the expected
+ * structure of the XML element.
+ ***************************************************************************/
+GModelSpectralComposite::GModelSpectralComposite(const GXmlElement& xml) : GModelSpectral()
+{
+    // Initialise members
+    init_members();
+
+    // Read information from XML element
+    read(xml);
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Copy constructor
+ *
+ * @param[in] model Spectral power law model.
+ ***************************************************************************/
+GModelSpectralComposite::GModelSpectralComposite(const GModelSpectralComposite& model)
+                                       : GModelSpectral(model)
+{
+    // Initialise members
+    init_members();
+
+    // Copy members
+    copy_members(model);
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Destructor
+ ***************************************************************************/
+GModelSpectralComposite::~GModelSpectralComposite(void)
+{
+    // Free members
+    free_members();
+
+    // Return
+    return;
+}
+
+
+/*==========================================================================
+ =                                                                         =
+ =                                Operators                                =
+ =                                                                         =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Assignment operator
+ *
+ * @param[in] model Spectral power law model.
+ * @return Spectral power law model.
+ ***************************************************************************/
+GModelSpectralComposite& GModelSpectralComposite::operator=(const GModelSpectralComposite& model)
+{
+    // Execute only if object is not identical
+    if (this != &model) {
+
+        // Copy base class members
+        this->GModelSpectral::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 spectral power law model
+ ***************************************************************************/
+void GModelSpectralComposite::clear(void)
+{
+    // Free class members (base and derived classes, derived class first)
+    free_members();
+    this->GModelSpectral::free_members();
+
+    // Initialise members
+    this->GModelSpectral::init_members();
+    init_members();
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Clone spectral power law model
+ *
+ * @return Pointer to deep copy of spectral power law model.
+ ***************************************************************************/
+GModelSpectralComposite* GModelSpectralComposite::clone(void) const
+{
+    // Clone spectral power law model
+    return new GModelSpectralComposite(*this);
+}
+
+
+/***********************************************************************//**
+ * @brief Evaluate function
+ *
+ * @param[in] srcEng True photon energy.
+ * @param[in] srcTime True photon arrival time.
+ * @param[in] gradients Compute gradients?
+ * @return Model value (ph/cm2/s/MeV).
+ *
+ * Evaluates
+ *
+ * \f[
+ *    \sum_{i=0}^{N} {M_{\rm i}}(\rm srcEng, srcTime)
+ * \f]
+ *
+ * where \f${M_{\rm i}}\f$ is the model evaluation of the i-th model
+ *
+ * If the @p gradients flag is true the method will also compute the
+ * partial derivatives of each model component with respect to the parameters using
+ *
+ * @todo The method expects that energy!=0. Otherwise Inf or NaN may result.
+ ***************************************************************************/
+double GModelSpectralComposite::eval(const GEnergy& srcEng,
+                                const GTime&   srcTime,
+                                const bool&    gradients) const
+{
+    // Initialise result
+	double value = 0.0;
+
+	// Loop over model components
+	for (int i = 0; i < m_spectral.size(); ++i) {
+
+		// Add model component
+		value += m_spectral[i]->eval(srcEng, srcTime, gradients);
+
+	} // endfor: loop over model components
+
+    // Compile option: Check for NaN/Inf
+    #if defined(G_NAN_CHECK)
+    if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
+        std::cout << "*** ERROR: GModelSpectralComposite::eval";
+        std::cout << "(srcEng=" << srcEng;
+        std::cout << ", srcTime=" << srcTime << "):";
+        std::cout << " NaN/Inf encountered";
+        std::cout << " (value=" << value;
+        std::cout << ")" << std::endl;
+    }
+    #endif
+
+    // Return
+    return value;
+}
+
+
+/***********************************************************************//**
+ * @brief Returns model photon flux between [emin, emax] (units: ph/cm2/s)
+ *
+ * @param[in] emin Minimum photon energy.
+ * @param[in] emax Maximum photon energy.
+ * @return Photon flux (ph/cm2/s).
+ *
+ * Computes the sum of all photon fluxes of each individual component
+ *
+ ***************************************************************************/
+double GModelSpectralComposite::flux(const GEnergy& emin,
+                                const GEnergy& emax) const
+{
+    // Initialise flux
+    double flux = 0.0;
+
+    // Compute only if integration range is valid
+    if (emin < emax) {
+
+    	// Loop over model components
+    	for (int i = 0; i < m_spectral.size(); ++i) {
+
+    		// Add model component
+    		flux += m_spectral[i]->flux(emin, emax);
+
+    	} // endfor: loop over model components
+
+    } // endif: integration range was valid
+
+    // Return flux
+    return flux;
+}
+
+
+/***********************************************************************//**
+ * @brief Returns model energy flux between [emin, emax] (units: erg/cm2/s)
+ *
+ * @param[in] emin Minimum photon energy.
+ * @param[in] emax Maximum photon energy.
+ * @return Energy flux (erg/cm2/s).
+ *
+ * Computes the sum of all energy fluxes of each individual component
+ *
+ ***************************************************************************/
+double GModelSpectralComposite::eflux(const GEnergy& emin,
+                                 const GEnergy& emax) const
+{
+   // Initialise eflux
+	double eflux = 0.0;
+
+	// Compute only if integration range is valid
+	if (emin < emax) {
+
+		// Loop over model components
+		for (int i = 0; i < m_spectral.size(); ++i) {
+
+			// Add model component
+			eflux += m_spectral[i]->flux(emin, emax);
+
+		} // endfor: loop over model components
+
+	} // endif: integration range was valid
+
+	// Return flux
+	return eflux;
+}
+
+
+/***********************************************************************//**
+ * @brief Returns Monte Carlo energy between [emin, emax]
+ *
+ * @param[in] emin Minimum photon energy.
+ * @param[in] emax Maximum photon energy.
+ * @param[in] time True photon arrival time.
+ * @param[in,out] ran Random number generator.
+ * @return Energy.
+ *
+ * @exception GException::erange_invalid
+ *            Energy range is invalid (emin < emax required).
+ *
+ * Returns Monte Carlo energy by randomly drawing from a power law.
+ ***************************************************************************/
+GEnergy GModelSpectralComposite::mc(const GEnergy& emin,
+                               const GEnergy& emax,
+                               const GTime&   time,
+                               GRan&          ran) const
+{
+    // Throw an exception if energy range is invalid
+    if (emin >= emax) {
+        throw GException::erange_invalid(G_MC, emin.MeV(), emax.MeV(),
+              "Minimum energy < maximum energy required.");
+    }
+
+    // Throw exception if model container is empty
+    if (m_spectral.size() == 0) {
+    	throw GException::runtime_error(G_MC,
+    	              "Composite spectral model is empty. It is required to have at"\
+					  " least one spectral model to run simulations");
+    }
+
+    // Update MC cache
+    update_mc_cache(emin, emax);
+
+    // Get random value
+    double u = ran.uniform();
+
+    // Initialise index of model to use for simulation
+    int index = 0;
+
+    // Loop over spectral components and compute relative probabilites
+    for (int i = 0; i < m_mc_probs.size(); ++i) {
+
+    	if (u <= m_mc_probs[i]) {
+    		index = i;
+    		break;
+    	}
+    }
+
+    // Set energy
+    GEnergy energy = m_spectral[index]->mc(emin, emax, time, ran);
+
+    // Return energy
+    return energy;
+}
+
+
+/***********************************************************************//**
+ * @brief Read model from XML element
+ *
+ * @param[in] xml XML element.
+ *
+ * Reads the spectral information from an XML element.
+ ***************************************************************************/
+void GModelSpectralComposite::read(const GXmlElement& xml)
+{
+	// Get number of spectral components
+	int n_spectrals = xml.elements("spectrum");
+
+	// Loop over spectral elements
+	for(int i = 0; i < n_spectrals; ++i) {
+
+		// Get spectral XML element
+		const GXmlElement* spec = xml.element("spectrum", i);
+
+		// Initialise a spectral registry object
+	    GModelSpectralRegistry registry;
+
+		// Read spectral model
+		GModelSpectral* ptr = registry.alloc(*spec);
+
+		// Get component attribute from XML file
+		std::string component_name = spec->attribute("component");
+
+		// Append spectral component to container
+		append(*ptr, component_name);
+
+	} // endfor: loop over components
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Write model into XML element
+ *
+ * @param[in] xml XML element.
+ *
+ * @exception GException::model_invalid_spectral
+ *            Existing XML element is not of the expected type.
+ *
+ * Writes the spectral information into an XML element.
+ ***************************************************************************/
+void GModelSpectralComposite::write(GXmlElement& xml) const
+{
+    // Set model type
+    if (xml.attribute("type") == "") {
+        xml.attribute("type", type());
+    }
+
+    // Verify model type
+    if (xml.attribute("type") != type()) {
+        throw GException::model_invalid_spectral(G_WRITE, xml.attribute("type"),
+              "Spectral model is not of type \""+type()+"\".");
+    }
+
+    // Loop over model components
+    for(int i = 0; i < m_spectral.size(); i++) {
+
+    	// Write spectral model
+		if (m_spectral[i] != NULL) {
+
+			// Create new spectrum node
+			xml.append(GXmlElement("spectrum"));
+
+			// Get new spectrum node
+			GXmlElement* spec = xml.element("spectrum", xml.elements("spectrum")-1);
+
+			// Create temporary copy of the spectral model
+			// This is a kluge to write out the original parameters.
+			GModelSpectral* cpy = m_spectral[i]->clone();
+
+			// Loop over parameters of model
+			for (int j = 0; j < cpy->size(); ++j) {
+
+				// Get model parameter and name
+				GModelPar& par = (*cpy)[j];
+				std::string parname = par.name();
+
+				// Check if name contains colon
+				if (gammalib::contains(parname, ":")) {
+
+					// Split at the colon
+					std::vector<std::string> splits = gammalib::split(parname, ":");
+
+					// Use second part of the string to recover original parameter name
+					par.name(splits[1]);
+				}
+
+			} // endfor: loop over parameters
+
+			// Write spectral component
+			cpy->write(*spec);
+
+			// Add component name if previously available
+			if (m_components[i] != gammalib::str(i+1)) {
+				spec->attribute("component", m_components[i]);
+			}
+
+			// Remove temporary copy
+			delete cpy;
+
+		} // endif: spectral model was not NULL
+
+    } // endfor: loop over model components
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Print powerlaw information
+ *
+ * @param[in] chatter Chattiness (defaults to NORMAL).
+ * @return String containing model information.
+ ***************************************************************************/
+std::string GModelSpectralComposite::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("=== GModelSpectralComposite ===");
+
+        // Append information
+		result.append("\n"+gammalib::parformat("Number of components"));
+		result.append(gammalib::str(components()));
+		result.append("\n"+gammalib::parformat("Number of parameters"));
+		result.append(gammalib::str(size()));
+
+		// Print parameter information
+		for (int i = 0; i < size(); ++i) {
+			result.append("\n"+m_pars[i]->print(chatter));
+		}
+
+    } // endif: chatter was not silent
+
+    // Return result
+    return result;
+}
+
+
+/***********************************************************************//**
+ * @brief Append spectral component
+ *
+ * @param[in] spec Spectral model component
+ * @param[in] name Name of spectral component (can be empty)
+ *
+ * Appends a spectral component to the composite model
+ ***************************************************************************/
+void GModelSpectralComposite::append(const GModelSpectral& spec, const std::string& name)
+{
+    // Append model container
+    m_spectral.push_back(spec.clone());
+
+    // Get index of latest model
+    int index = m_spectral.size()-1;
+
+    // Use model index if component name is empty
+    std::string component_name = !name.empty() ? name : gammalib::str(m_spectral.size());
+
+    // Check if component name is unique, throw exception if not
+    if (gammalib::contains(m_components, component_name)) {
+    	std::string msg =
+			"Attempt to append component with name \""+component_name+"\" to composite"
+			" spectral model container, but a component with the same name exists already.\n"
+			"Every component in the container needs a unique name. On default the system will"
+			"increment an integer if no component name is provided.";
+		throw GException::invalid_value(G_APPEND, msg);
+    }
+
+	// Add component name (for now simple number)
+	m_components.push_back(component_name);
+
+	// Get number of spectral parameters from model
+	int npars = m_spectral[index]->size();
+
+	// Loop over model parameters
+	for (int ipar = 0; ipar < npars; ++ipar) {
+
+		// Get model parameter
+		GModelPar* par = &(m_spectral[index]->operator[](ipar));
+
+		// Modify parameter name
+		par->name(component_name+":"+par->name());
+
+		// Append model parameter with new name to internal container
+		m_pars.push_back(par);
+
+	} // endfor: loop over model parameters
+}
+
+
+/***********************************************************************//**
+ * @brief Returns spectral component element
+ *
+ * @param[in] index Index of spectral component.
+ * @return Spectral model.
+ *
+ * Returns a spectral component to the composite model
+ ***************************************************************************/
+const GModelSpectral* GModelSpectralComposite::component(const int& index) const
+{
+	// Check if index is in validity range
+	if (index >= m_spectral.size() || index < 0) {
+		throw GException::out_of_range(G_COMPONENT_INDEX, "Component Index", index, m_spectral.size(),"");
+	}
+
+	// Return spectral component
+	return m_spectral[index];
+
+}
+
+
+/***********************************************************************//**
+ * @brief Returns pointer to specific spectral component
+ *
+ * @param[in] name Name of spectral component.
+ * @return Spectral model.
+ *
+ * Returns a spectral component of the composite model
+ ***************************************************************************/
+const GModelSpectral* GModelSpectralComposite::component(const std::string& name) const
+{
+	// Check if model name is found
+	int index = -1;
+	for(int i = 0; i < m_components.size(); ++i) {
+		if (m_components[i] == name) {
+			index = i;
+			break;
+		}
+	}
+
+	// Check if component name was found
+	if (index == -1) {
+		throw GException::model_not_found(G_COMPONENT_NAME, name,"");
+	}
+
+	// Return spectral component
+	return m_spectral[index];
+
+}
+
+/*==========================================================================
+ =                                                                         =
+ =                             Private methods                             =
+ =                                                                         =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Initialise class members
+ ***************************************************************************/
+void GModelSpectralComposite::init_members(void)
+{
+    // Initialise model type
+    m_type = "Composite";
+
+    // Clear spectral models
+    m_spectral.clear();
+    m_components.clear();
+
+    // Clear MC cache
+    m_mc_probs.clear();
+    m_mc_emin.clear();
+    m_mc_emax.clear();
+    m_mc_flux = 0.0;
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Copy class members
+ *
+ * @param[in] model GModelSpectralComposite members which should be copied.
+ ***************************************************************************/
+void GModelSpectralComposite::copy_members(const GModelSpectralComposite& model)
+{
+    // Copy members
+    m_type       = model.m_type;
+    m_components = model.m_components;
+
+    // Copy MC cache
+    m_mc_probs   = model.m_mc_probs;
+    m_mc_emin    = model.m_mc_emin;
+    m_mc_emax    = model.m_mc_emax;
+    m_mc_flux    = model.m_mc_flux;
+
+    // Copy pointer(s) of spectral component
+    m_spectral.clear();
+
+    // Clear parameters and copy the pointers from the clone
+    m_pars.clear();
+    for(int i = 0; i < model.components(); ++i) {
+    	m_spectral.push_back(model.m_spectral[i]->clone());
+
+    	// Retrieve spectral model
+    	GModelSpectral* spec = m_spectral[i];
+
+    	// Loop over parameters and store pointers
+		for (int ipar = 0; ipar < spec->size(); ++ipar) {
+
+			// Get model parameter reference
+			GModelPar& par = spec->operator[](ipar);
+
+			// Append model parameter pointer to internal container
+			m_pars.push_back(&par);
+
+		} // endfor: loop over parameters
+
+    } // endfor: loop over spectral models
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Delete class members
+ ***************************************************************************/
+void GModelSpectralComposite::free_members(void)
+{
+    // Free memory
+	for(int i = 0; i < m_spectral.size(); ++i) {
+
+		// Delete component i
+		if (m_spectral[i] != NULL) delete m_spectral[i];
+
+		// Signal free pointer
+		m_spectral[i] = NULL;
+	}
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Update Monte Carlo pre computation cache
+ *
+ * @param[in] emin Minimum photon energy.
+ * @param[in] emax Maximum photon energy.
+ *
+ * Updates the precomputation cache for Monte Carlo simulations.
+ ***************************************************************************/
+void GModelSpectralComposite::update_mc_cache(const GEnergy& emin,
+                                         const GEnergy& emax) const
+
+{
+	// Update cache if energy range has changed
+	if (emin != m_mc_emin || emax != m_mc_emax) {
+		m_mc_emin = emin;
+		m_mc_emax = emax;
+		m_mc_flux = flux(m_mc_emin, m_mc_emax);
+
+		// Initialise sum and probabilites
+	    double sum = 0.0;
+	    m_mc_probs.clear();
+
+	    // Loop over spectral components and compute relative probabilites
+	    for (int i = 0; i < m_spectral.size(); ++i) {
+
+	    	// Relative probability
+	    	double prob = m_spectral[i]->flux(emin, emax) / m_mc_flux;
+
+	    	// Add probability
+	    	m_mc_probs.push_back(prob + sum);
+
+	    	// Increment sum
+	    	sum += prob;
+
+	    } //endfor: looped over spectral components
+
+	} // endif: emin and emax have changed
+
+    // Return
+    return;
+}
@@ -45,6 +45,7 @@ sources = GModelPar.cpp \
           GModelSpatialDiffuseConst.cpp \
           GModelSpatialDiffuseCube.cpp \
           GModelSpatialDiffuseMap.cpp \
+          GModelSpatialComposite.cpp \
           GModelSpectral.cpp \
           GModelSpectralRegistry.cpp \
           GModelSpectralBrokenPlaw.cpp \
@@ -60,6 +61,7 @@ sources = GModelPar.cpp \
           GModelSpectralPlaw.cpp \
           GModelSpectralPlawPhotonFlux.cpp \
           GModelSpectralPlawEnergyFlux.cpp \
+          GModelSpectralComposite.cpp \
           GModelTemporal.cpp \
           GModelTemporalRegistry.cpp \
           GModelTemporalConst.cpp \
@@ -43,7 +43,6 @@
 #include "GEbounds.hpp"       // will become obsolete
 #include "GObservation.hpp"
 #include "GModelSky.hpp"
-#include "GModelSpatialPointSource.hpp"
  
 /* __ Method name definitions ____________________________________________ */
 #define G_IRF_RADIAL               "GResponse::irf_radial(GEvent&, GSource&,"\
@@ -322,9 +321,8 @@ double GResponse::eval_prob(const GModelSky&amp;    model,
  
         // Set source
         GSource source(model.name(), model.spatial(), srcEng, srcTime);
-        
-        // Get IRF value. This method returns the spatial component of the
-        // source model.
+
+        // Compute IRF value
         double irf = this->irf(event, source, obs);
  
         // If required, apply instrument specific model scaling
@@ -1180,7 +1180,7 @@ bool gammalib::contains(const std::string&amp; str, const std::string&amp; substring)
     bool result = false;
  
     // checks if substring is in str
-    if (str.find(substring) != std::string::npos){
+    if (str.find(substring) != std::string::npos) {
         result = true;
     }
  
@@ -1193,13 +1193,13 @@ bool gammalib::contains(const std::string&amp; str, const std::string&amp; substring)
  * @brief Checks if a string is contained in a vector of strings
  *
  * @param[in] strings Vector of strings you want to search in.
- * @param[in] string String you are looking for in strings.
+ * @param[in] string string you are looking for in strings.
  * @return True if a string is contained a vector.
  *
  * Checks if a string is contained in a vector of strings
  ***************************************************************************/
-bool gammalib::contains(const std::vector<std::string>& strings,
-                        const std::string&              string)
+bool gammalib::contains(const std::vector<std::string> strings,
+                        const std::string&             string)
 {
     // Compute result
     bool result = std::find(strings.begin(), strings.end(), string) !=
+<?xml version="1.0" standalone="no"?>
+<source_library title="source library">
+  <source name="Crab" type="CompositeSource">
+    <spectrum type="PowerLaw">
+      <parameter name="Prefactor"   scale="1e-16" value="5.7"  min="1e-07" max="1000.0" free="1"/>
+      <parameter name="Index"       scale="-1"    value="2.48" min="0.0"   max="+5.0"   free="1"/>
+      <parameter name="PivotEnergy" scale="1e6"   value="0.3"  min="0.01"  max="1000.0" free="0"/>
+    </spectrum>
+    <spatialModel type="Composite">
+      <spatialModel type="PointSource" component="PointSource">
+        <parameter name="RA"    scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
+        <parameter name="DEC"   scale="1.0" value="22.0145" min="-90"  max="90"  free="1"/>
+      </spatialModel>
+      <spatialModel type="RadialGaussian">
+        <parameter name="RA"    scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
+        <parameter name="DEC"   scale="1.0" value="22.0145" min="-90"  max="90"  free="1"/>
+        <parameter name="Sigma" scale="1.0" value="0.20"    min="0.01" max="10"  free="1"/>
+      </spatialModel>
+    </spatialModel>
+  </source>
+</source_library>
+<?xml version="1.0" standalone="no"?>
+<source_library title="source library">
+  <source name="Crab" type="PointSource">
+    <spectrum type="Composite">
+      <spectrum type="PowerLaw" component="SoftComponent">     
+        <parameter name="Prefactor"   scale="1e-17" value="3"  min="1e-07" max="1000.0" free="1"/>
+        <parameter name="Index"       scale="-1"    value="3.5" min="0.0"   max="+5.0"   free="1"/>
+        <parameter name="PivotEnergy" scale="1e6"   value="1"  min="0.01"  max="1000.0" free="0"/>
+      </spectrum>
+      <spectrum type="PowerLaw" component="HardComponent">     
+        <parameter name="Prefactor"   scale="1e-17" value="5"  min="1e-07" max="1000.0" free="1"/>
+        <parameter name="Index"       scale="-1"    value="2.0" min="0.0"   max="+5.0"   free="1"/>
+        <parameter name="PivotEnergy" scale="1e6"   value="1"  min="0.01"  max="1000.0" free="0"/>
+      </spectrum>
+    </spectrum>
+    <spatialModel type="PointSource">
+      <parameter free="0" max="360" min="-360" name="RA" scale="1" value="83.6331" />
+      <parameter free="0" max="90" min="-90" name="DEC" scale="1" value="22.0145" />
+    </spatialModel>
+  </source>
+</source_library>
@@ -67,6 +67,7 @@ void TestGModel::set(void)
     m_xml_point_multiplicative     = datadir + "/model_point_multiplicative.xml";
     m_xml_model_point_nodes        = datadir + "/model_point_nodes.xml";
     m_xml_model_point_filefct      = datadir + "/model_point_filefct.xml";
+    m_xml_model_spectral_composite = datadir + "/model_spectral_composite.xml";
     m_xml_model_diffuse_const      = datadir + "/model_diffuse_const.xml";
     m_xml_model_diffuse_cube       = datadir + "/model_diffuse_cube.xml";
     m_xml_model_diffuse_map        = datadir + "/model_diffuse_map.xml";
@@ -75,6 +76,7 @@ void TestGModel::set(void)
     m_xml_model_radial_shell       = datadir + "/model_radial_shell.xml";
     m_xml_model_elliptical_disk    = datadir + "/model_elliptical_disk.xml";
     m_xml_model_elliptical_gauss   = datadir + "/model_elliptical_gauss.xml";
+    m_xml_model_spatial_composite  = datadir + "/model_spatial_composite.xml";
  
     // Set legacy model definition XML files
     m_xml_legacy_radial_disk       = datadir + "/legacy_radial_disk.xml";
@@ -116,6 +118,8 @@ void TestGModel::set(void)
            "Test GModelSpatialDiffuseCube");
     append(static_cast<pfunction>(&TestGModel::test_diffuse_map),
            "Test GModelSpatialDiffuseMap");
+    append(static_cast<pfunction>(&TestGModel::test_spatial_composite),
+           "Test GModelSpatialComposite");
     append(static_cast<pfunction>(&TestGModel::test_spatial_model),
            "Test spatial model XML I/O");
  
@@ -146,6 +150,8 @@ void TestGModel::set(void)
            "Test GModelSpectralNodes");
     append(static_cast<pfunction>(&TestGModel::test_filefct),
            "Test GModelSpectralFunc");
+    append(static_cast<pfunction>(&TestGModel::test_spectral_composite),
+           "Test GModelSpectralComposite");
     append(static_cast<pfunction>(&TestGModel::test_spectral_model),
            "Test spectral model XML I/O");
  
@@ -1402,6 +1408,49 @@ void TestGModel::test_elliptical_gauss(void)
 }
  
  
+/***********************************************************************//**
+ * @brief Test GModelSpatialComposite class
+ ***************************************************************************/
+void TestGModel::test_spatial_composite(void)
+{
+    // Test void constructor
+    GModelSpatialComposite model1;
+    test_value(model1.type(), "Composite", "Check void model type");
+
+    // Set sky direction
+    GSkyDir dir = GSkyDir();
+    dir.radec_deg(83.6331, 22.01);
+
+    // Test append method
+    model1.append(GModelSpatialPointSource(dir));
+    model1.append(GModelSpatialRadialGauss(dir, 0.2));
+    test_value(model1.components(), 2);
+    test_value(model1.size(), 5);
+
+    // Test XML constructor
+    GXml         xml(m_xml_model_spatial_composite);
+    GXmlElement* element = xml.element(0)->element(0)->element("spatialModel", 0);
+    GModelSpatialComposite model2(*element);
+    test_value(model2.size(), 5);
+    test_value(model2.type(), "Composite", "Check model type");
+      
+    // Test access of individual parameters
+    test_value(model2["2:RA"].value(), 83.6331);
+    test_value(model2["2:DEC"].value(), 22.0145);
+    test_value(model2["2:Sigma"].value(), 0.2);
+    test_value(model2["PointSource:RA"].value(), 83.6331);
+    test_value(model2["2:RA"].value(), 83.6331);
+    model2["2:RA"].value(83.1331);
+    test_value(model2["2:RA"].value(), 83.1331);
+    model2["2:DEC"].value(22.51);
+    test_value(model2["2:DEC"].value(), 22.51);
+    model2["2:Sigma"].value(0.6);
+    test_value(model2["2:Sigma"].value(), 0.6);
+
+    // Exit test
+    return;
+}
+
  
 /***********************************************************************//**
  * @brief Test GModelSpectralConst class
@@ -2194,6 +2243,51 @@ void TestGModel::test_filefct(void)
     return;
 }
  
+/***********************************************************************//**
+ * @brief Test GModelSpectralComposite class
+ ***************************************************************************/
+void TestGModel::test_spectral_composite(void)
+{
+    // Test void constructor
+    GModelSpectralComposite model1;
+    test_value(model1.type(), "Composite", "Check void model type");
+
+    // Test append method
+    model1.append(GModelSpectralPlaw(3e-17,-3.5, GEnergy(1, "TeV")));
+    model1.append(GModelSpectralPlaw(5e-17,-2.0, GEnergy(1, "TeV")));
+    test_value(model1.size(), 6);
+    test_value(model1.components(), 2);
+    test_value(model1.eval(GEnergy(1.0, "TeV")), 8.0e-17);
+
+    // Test XML constructor
+    GXml                      xml(m_xml_model_spectral_composite);
+    GXmlElement*              element = xml.element(0)->element(0)->element("spectrum", 0);
+    GModelSpectralComposite model2(*element);
+    test_value(model2.size(), 6);
+    test_value(model2.components(), 2);
+    test_value(model2.type(), "Composite", "Check model type");
+
+    // Test access of individual parameters
+    test_value(model2["HardComponent:Prefactor"].value(), 5e-17);
+    test_value(model2["HardComponent:Index"].value(), -2.0);
+    test_value(model2["HardComponent:PivotEnergy"].value(), 1.0e6);
+
+    // Test prefactor method
+    model2["SoftComponent:Prefactor"].value(2.3e-16);
+    test_value(model2["SoftComponent:Prefactor"].value(), 2.3e-16);
+
+    // Test index method
+    model2["SoftComponent:Index"].value(-2.6);
+    test_value(model2["SoftComponent:Index"].value(), -2.6);
+
+    // Test pivot method
+    model2["SoftComponent:PivotEnergy"].value(0.5e6);
+    test_value(model2["SoftComponent:PivotEnergy"].value(), 0.5e6);
+
+    // Exit test
+    return;
+}
+
  
 /***********************************************************************//**
  * @brief Test GModelTemporalConst class
@@ -2297,6 +2391,7 @@ void TestGModel::test_spectral_model(void)
     test_xml_model("GModelSpectralLogParabola",    m_xml_model_point_logparabola);
     test_xml_model("GModelSpectralNodes",          m_xml_model_point_nodes);
     test_xml_model("GModelSpectralFunc",           m_xml_model_point_filefct);
+    test_xml_model("GModelSpectralComposite",      m_xml_model_spectral_composite);
  
     // Return
     return;
@@ -58,6 +58,7 @@ public:
     void                test_radial_shell(void);
     void                test_elliptical_disk(void);
     void                test_elliptical_gauss(void);
+    void                test_spatial_composite(void);
     void                test_spatial_model(void);
     void                test_const(void);
     void                test_gauss(void);
@@ -72,6 +73,7 @@ public:
     void                test_multiplicative(void);
     void                test_nodes(void);
     void                test_filefct(void);
+    void                test_spectral_composite(void);
     void                test_spectral_model(void);
     void                test_temp_const(void);
     void                test_model(void);
@@ -115,6 +117,7 @@ private:
     std::string m_xml_point_multiplicative;
     std::string m_xml_model_point_nodes;
     std::string m_xml_model_point_filefct;
+    std::string m_xml_model_spectral_composite;
     std::string m_xml_model_diffuse_const;
     std::string m_xml_model_diffuse_cube;
     std::string m_xml_model_diffuse_map;
@@ -123,6 +126,7 @@ private:
     std::string m_xml_model_radial_shell;
     std::string m_xml_model_elliptical_disk;
     std::string m_xml_model_elliptical_gauss;
+    std::string m_xml_model_spatial_composite;
  
     // Legacy model definition XML files
     std::string m_xml_legacy_radial_disk;