Support pickling of COMPTEL classes (#1938)

Jürgen Knödlseder
1 parent 24f4336c
Showing 51 changed files with 1035 additions and 195 deletions
NEWS
inst/com/include/GCOMD1Response.hpp
inst/com/include/GCOMD2Response.hpp
inst/com/include/GCOMEventBin.hpp
inst/com/include/GCOMIaq.hpp
inst/com/include/GCOMObservation.hpp
inst/com/include/GCOMResponse.hpp
inst/com/include/GCOMTim.hpp
inst/com/pyext/GCOMD1Response.i
inst/com/pyext/GCOMD2Response.i
inst/com/pyext/GCOMDri.i
inst/com/pyext/GCOMEventAtom.i
inst/com/pyext/GCOMEventBin.i
inst/com/pyext/GCOMEventCube.i
inst/com/pyext/GCOMEventList.i
inst/com/pyext/GCOMIaq.i
inst/com/pyext/GCOMInstDir.i
inst/com/pyext/GCOMModelDRBFitting.i
inst/com/pyext/GCOMOad.i
inst/com/pyext/GCOMOads.i
@@ -92,6 +92,16 @@ The following methods have been added:
 - GCTAModelRadialAcceptance::temporal(GModelTemporal*)
 - GCTAModelIrfBackground::spectral(GModelSpectral*)
 - GCTAModelIrfBackground::temporal(GModelTemporal*)
+- GCOMD1Response::write(GFitsBinTable&)
+- GCOMD2Response::write(GFitsBinTable&)
+- GCOMEventBin::index(int&)
+- GCOMEventBin::dir(GCOMInstDir&)
+- GCOMEventBin::energy(GEnergy&)
+- GCOMEventBin::time(GTime&)
+- GCOMEventBin::counts(double&)
+- GCOMEventBin::solidangle(double&)
+- GCOMEventBin::ewidth(GEnergy&)
+- GCOMEventBin::ontime(double&)
 The following methods have been removed:
 - GSkyRegionMap::filename()
 /***************************************************************************
  *           GCOMD1Response.hpp - COMPTEL D1 module response class         *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2017 by Juergen Knoedlseder                              *
+ *  copyright (C) 2017-2018 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -37,6 +37,7 @@
 /* __ Forward declarations _______________________________________________ */
 class GFitsTable;
+class GFitsBinTable;
 /***********************************************************************//**
@@ -65,6 +66,7 @@ public:
     const GCaldb&   caldb(void) const;
     void            load(const std::string& sdaname);
     void            read(const GFitsTable& table);
+    void            write(GFitsBinTable& table);
     double          position(const double& etrue) const;
     double          sigma(const double& etrue) const;
     double          amplitude(const double& etrue) const;
 /***************************************************************************
  *           GCOMD2Response.hpp - COMPTEL D2 module response class         *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2017 by Juergen Knoedlseder                              *
+ *  copyright (C) 2017-2018 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -38,6 +38,7 @@
 /* __ Forward declarations _______________________________________________ */
 class GFitsTable;
+class GFitsBinTable;
 /***********************************************************************//**
@@ -66,6 +67,7 @@ public:
     const GCaldb&   caldb(void) const;
     void            load(const std::string& sdbname);
     void            read(const GFitsTable& table);
+    void            write(GFitsBinTable& table);
     double          position(const double& etrue) const;
     double          sigma(const double& etrue) const;
     double          amplitude(const double& etrue) const;
 /***************************************************************************
  *                GCOMEventBin.hpp - COMPTEL event bin class               *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2012-2016 by Juergen Knoedlseder                         *
+ *  copyright (C) 2012-2018 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -78,6 +78,13 @@ public:
     const double&  solidangle(void) const;
     const GEnergy& ewidth(void) const;
     const double&  ontime(void) const;
+    void           index(const int& index);
+    void           dir(const GCOMInstDir& dir);
+    void           energy(const GEnergy& energy);
+    void           time(const GTime& time);
+    void           solidangle(const double& solidangle);
+    void           ewidth(const GEnergy& ewidth);
+    void           ontime(const double& ontime);
 protected:
     // Protected methods
@@ -111,77 +118,6 @@ std::string GCOMEventBin::classname(void) const
 /***********************************************************************//**
- * @brief Return instrument direction
- *
- * @return Instrument direction.
- *
- * Returns the instrument direction of the event bin.
- ***************************************************************************/
-inline
-const GCOMInstDir& GCOMEventBin::dir(void) const
-{
-    return (*m_dir);
-}
-
-
-/***********************************************************************//**
- * @brief Return energy
- *
- * @return Energy.
- *
- * Returns the energy of the event bin.
- ***************************************************************************/
-inline
-const GEnergy& GCOMEventBin::energy(void) const
-{
-    return (*m_energy);
-}
-
-
-/***********************************************************************//**
- * @brief Return time
- *
- * @return Time.
- *
- * Returns the time of the event bin.
- ***************************************************************************/
-inline
-const GTime& GCOMEventBin::time(void) const
-{
-    return (*m_time);
-}
-
-
-/***********************************************************************//**
- * @brief Return number of counts
- *
- * @return Number of counts.
- *
- * Returns the number of counts in the event bin.
- ***************************************************************************/
-inline
-double GCOMEventBin::counts(void) const
-{
-    return (*m_counts);
-}
-
-
-/***********************************************************************//**
- * @brief Set number of counts
- *
- * @param[in] counts Number of counts.
- *
- * Set the number of counts in the event bin.
- ***************************************************************************/
-inline
-void GCOMEventBin::counts(const double& counts)
-{
-    *m_counts = counts;
-    return;
-}
-
-
-/***********************************************************************//**
  * @brief Return bin index
  *
  * @return Bin index.
@@ -197,44 +133,18 @@ const int&amp; GCOMEventBin::index(void) const
 /***********************************************************************//**
- * @brief Return solid angle
+ * @brief Set bin index
  *
- * @return Solid angle.
- *
- * Returns the solid angle of the event bin.
- ***************************************************************************/
-inline
-const double& GCOMEventBin::solidangle(void) const
-{
-    return (*m_solidangle);
-}
-
-
-/***********************************************************************//**
- * @brief Return energy width
- *
- * @return Energy width.
- *
- * Returns the energy width of the event bin.
- ***************************************************************************/
-inline
-const GEnergy& GCOMEventBin::ewidth(void) const
-{
-    return (*m_ewidth);
-}
-
-
-/***********************************************************************//**
- * @brief Return ontime
- *
- * @return Ontime.
+ * @return Bin index.
  *
- * Returns the ontime of the event bin.
+ * Set the index of the event bin if the bin is part of an event cube. If
+ * the event is not part of an event cube the index should be set to -1.
  ***************************************************************************/
 inline
-const double& GCOMEventBin::ontime(void) const
+void GCOMEventBin::index(const int& index)
 {
-    return (*m_ontime);
+    m_index = index;
+    return;
 }
 #endif /* GCOMEVENTBIN_HPP */
@@ -55,8 +55,8 @@ public:
     // Constructors and destructors
     GCOMIaq(void);
     GCOMIaq(const GCOMIaq& iaq);
-    GCOMIaq(const double&   phigeo_max, const double& phigeo_bin_size,
-            const double&   phibar_max, const double& phibar_bin_size);
+    GCOMIaq(const double& phigeo_max, const double& phigeo_bin_size,
+            const double& phibar_max, const double& phibar_bin_size);
     ~GCOMIaq(void);
     // Operators
 /***************************************************************************
  *            GCOMObservation.hpp - COMPTEL observation class              *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2012-2017 by Juergen Knoedlseder                         *
+ *  copyright (C) 2012-2018 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -66,6 +66,7 @@ class GCOMObservation : public GObservation {
 public:
     // Constructors and destructors
     GCOMObservation(void);
+    explicit GCOMObservation(const GXmlElement& xml);
     GCOMObservation(const GFilename& drename,
                     const GFilename& drbname,
                     const GFilename& drgname,
 /***************************************************************************
  *                 GCOMResponse.hpp - COMPTEL Response class               *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2012-2017 by Juergen Knoedlseder                         *
+ *  copyright (C) 2012-2018 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -42,6 +42,7 @@ class GEnergy;
 class GTime;
 class GObservation;
 class GFitsImage;
+class GFitsImageFloat;
 /***********************************************************************//**
@@ -89,6 +90,7 @@ public:
     const std::string& rspname(void) const;
     void               load(const std::string& rspname);
     void               read(const GFitsImage& hdu);
+    void               write(GFitsImageFloat& image) const;
 private:
     // Private methods
 /***************************************************************************
  *             GCOMTim.hpp - COMPTEL Good Time Intervals class             *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2017 by Juergen Knodlseder                               *
+ *  copyright (C) 2017-2018 by Juergen Knodlseder                          *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -51,6 +51,7 @@ class GCOMTim : public GBase {
 public:
     // Constructors and destructors
     GCOMTim(void);
+    explicit GCOMTim(const GGti& gti);
     GCOMTim(const GCOMTim& tim);
     GCOMTim(const GFilename& filename, const std::string& usage = "",
                                        const std::string& mode  = "");
 /***************************************************************************
  *           GCOMD1Response.i - COMPTEL D1 module response class           *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2017 by Juergen Knoedlseder                              *
+ *  copyright (C) 2017-2018 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -54,6 +54,7 @@ public:
     const GCaldb&   caldb(void) const;
     void            load(const std::string& sdaname);
     void            read(const GFitsTable& hdu);
+    void            write(GFitsBinTable& table);
     double          position(const double& etrue) const;
     double          sigma(const double& etrue) const;
     double          amplitude(const double& etrue) const;
@@ -72,4 +73,14 @@ public:
     GCOMD1Response copy() {
         return (*self);
     }
+%pythoncode {
+    def __getstate__(self):
+        table = gammalib.GFitsBinTable()
+        self.write(table)
+        state = (table,)
+        return state
+    def __setstate__(self, state):
+        self.__init__()
+        self.read(state[0])
+}
 };
 /***************************************************************************
  *           GCOMD2Response.i - COMPTEL D2 module response class           *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2017 by Juergen Knoedlseder                              *
+ *  copyright (C) 2017-2018 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -54,6 +54,7 @@ public:
     const GCaldb&   caldb(void) const;
     void            load(const std::string& sdaname);
     void            read(const GFitsTable& hdu);
+    void            write(GFitsBinTable& table);
     double          position(const double& etrue) const;
     double          sigma(const double& etrue) const;
     double          amplitude(const double& etrue) const;
@@ -76,4 +77,14 @@ public:
     GCOMD2Response copy() {
         return (*self);
     }
+%pythoncode {
+    def __getstate__(self):
+        table = gammalib.GFitsBinTable()
+        self.write(table)
+        state = (table,)
+        return state
+    def __setstate__(self, state):
+        self.__init__()
+        self.read(state[0])
+}
 };
 /***************************************************************************
  *                   GCOMDri.i - COMPTEL Data Space class                  *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2017 by Juergen Knoedlseder                              *
+ *  copyright (C) 2017-2018 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -116,4 +116,15 @@ public:
     GCOMDri copy() {
         return (*self);
     }
+%pythoncode {
+    def __getstate__(self):
+        fits = gammalib.GFits()
+        self.write(fits)
+        state = (fits,)
+        return state
+    def __setstate__(self, state):
+        self.__init__()
+        if not state[0].is_empty():
+            self.read(state[0][0])
+}
 };
 /***************************************************************************
  *                GCOMEventAtom.i - COMPTEL event atom class               *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2017 by Juergen Knoedlseder                              *
+ *  copyright (C) 2017-2018 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -87,4 +87,28 @@ public:
     GCOMEventAtom copy() {
         return (*self);
     }
+%pythoncode {
+    def __getstate__(self):
+        state = (self.dir(), self.energy(), self.time(),
+                 self.e1(), self.e2(), self.phibar(), self.theta(),
+                 self.phi(), self.eha(), self.psd(), self.tof(),
+                 self.modcom(), self.reflag(), self.veto())
+        return state
+    def __setstate__(self, state):
+        self.__init__()
+        self.dir(state[0])
+        self.energy(state[1])
+        self.time(state[2])
+        self.e1(state[3])
+        self.e2(state[4])
+        self.phibar(state[5])
+        self.theta(state[6])
+        self.phi(state[7])
+        self.eha(state[8])
+        self.psd(state[9])
+        self.tof(state[10])
+        self.modcom(state[11])
+        self.reflag(state[12])
+        self.veto(state[13])
+}
 };
 /***************************************************************************
  *                GCOMEventBin.i  -  COMPTEL event bin class               *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2012-2015 by Juergen Knoedlseder                         *
+ *  copyright (C) 2012-2018 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -58,6 +58,13 @@ public:
     const double&  solidangle(void) const;
     const GEnergy& ewidth(void) const;
     const double&  ontime(void) const;
+    void           index(const int& index);
+    void           dir(const GCOMInstDir& dir);
+    void           energy(const GEnergy& energy);
+    void           time(const GTime& time);
+    void           solidangle(const double& solidangle);
+    void           ewidth(const GEnergy& ewidth);
+    void           ontime(const double& ontime);
 };
@@ -68,4 +75,21 @@ public:
     GCOMEventBin copy() {
         return (*self);
     }
+%pythoncode {
+    def __getstate__(self):
+        state = (self.index(), self.dir(), self.energy(), self.time(),
+                 self.counts(), self.solidangle(), self.ewidth(),
+                 self.ontime())
+        return state
+    def __setstate__(self, state):
+        self.__init__()
+        if state[0] == -1:
+            self.dir(state[1])
+            self.energy(state[2])
+            self.time(state[3])
+            self.counts(state[4])
+            self.solidangle(state[5])
+            self.ewidth(state[6])
+            self.ontime(state[7])
+}
 };
 /***************************************************************************
  *          GCOMEventCube.i  -  COMPTEL event bin container class          *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2012-2017 by Juergen Knoedlseder                         *
+ *  copyright (C) 2012-2018 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -66,9 +66,6 @@ public:
  * @brief GCOMEventCube class extension
  ***************************************************************************/
 %extend GCOMEventCube {
-    GCOMEventCube copy() {
-        return (*self);
-    }
     GCOMEventBin* __getitem__(int index) {
         if (index >= 0 && index < self->size()) {
             return (*self)[index];
@@ -85,4 +82,18 @@ public:
             throw GException::out_of_range("__setitem__(int)", index, self->size());
         }
     }
+    GCOMEventCube copy() {
+        return (*self);
+    }
+%pythoncode {
+    def __getstate__(self):
+        fits = gammalib.GFits()
+        self.write(fits)
+        state = (fits,)
+        return state
+    def __setstate__(self, state):
+        self.__init__()
+        if not state[0].is_empty():
+            self.read(state[0])
+}
 };
 /***************************************************************************
  *                GCOMEventList.i - COMPTEL event list class               *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2017 by Juergen Knoedlseder                              *
+ *  copyright (C) 2017-2018 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -26,7 +26,6 @@
 %{
 /* Put headers and other declarations here that are needed for compilation */
 #include "GCOMEventList.hpp"
-#include "GTools.hpp"
 %}
@@ -129,4 +128,15 @@ public:
     GCOMEventList copy() {
         return (*self);
     }
+%pythoncode {
+    def __getstate__(self):
+        fits = gammalib.GFits()
+        self.write(fits)
+        state = (fits,)
+        return state
+    def __setstate__(self, state):
+        self.__init__()
+        if not state[0].is_empty():
+            self.read(state[0])
+}
 };
 /***************************************************************************
  *          GCOMIaq.i - COMPTEL instrument response representation         *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2017 by Juergen Knoedlseder                              *
+ *  copyright (C) 2017-2018 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -40,8 +40,8 @@ public:
     // Constructors and destructors
     GCOMIaq(void);
     GCOMIaq(const GCOMIaq& iaq);
-    GCOMIaq(const double&   phigeo_max, const double& phigeo_bin_size,
-            const double&   phibar_max, const double& phibar_bin_size);
+    GCOMIaq(const double& phigeo_max, const double& phigeo_bin_size,
+            const double& phibar_max, const double& phibar_bin_size);
     ~GCOMIaq(void);
     // Methods
 /***************************************************************************
  *           GCOMInstDir.i - COMPTEL instrument direction class            *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2012-2017 by Juergen Knoedlseder                         *
+ *  copyright (C) 2012-2018 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -62,4 +62,11 @@ public:
     GCOMInstDir copy() {
         return (*self);
     }
+%pythoncode {
+    def __getstate__(self):
+        state = (self.dir(), self.phibar())
+        return state
+    def __setstate__(self, state):
+        self.__init__(state[0], state[1])
+}
 };
 /***************************************************************************
  *         GCOMModelDRBFitting.i - COMPTEL DRB model fitting class         *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2012-2016 by Juergen Knoedlseder                         *
+ *  copyright (C) 2012-2018 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -66,4 +66,16 @@ public:
     GCOMModelDRBFitting copy() {
         return (*self);
     }
+%pythoncode {
+    def __getstate__(self):
+        xml = gammalib.GXmlElement()
+        self.write(xml)
+        state = (xml,)
+        return state
+    def __setstate__(self, state):
+        if state[0][0].elements('node') == 0:
+            self.__init__()
+        else:
+            self.__init__(state[0][0])
+}
 };
@@ -78,4 +78,22 @@ public:
     GCOMOad copy() {
         return (*self);
     }
+%pythoncode {
+    def __getstate__(self):
+        state = (self.tstart(), self.tstop(), self.zaxis(), self.xaxis(),
+                 self.tjd(), self.tics(), self.gcaz(), self.gcel(),
+                 self.georad())
+        return state
+    def __setstate__(self, state):
+        self.__init__()
+        self.tstart(state[0])
+        self.tstop(state[1])
+        self.zaxis(state[2])
+        self.xaxis(state[3])
+        self.tjd(state[4])
+        self.tics(state[5])
+        self.gcaz(state[6])
+        self.gcel(state[7])
+        self.georad(state[8])
+}
 };
@@ -121,4 +121,15 @@ public:
     GCOMOads copy() {
         return (*self);
     }
+%pythoncode {
+    def __getstate__(self):
+        state = tuple([x for x in self]),
+        return state
+    def __setstate__(self, state):
+        self.__init__()
+        size = len(state[0])
+        self.reserve(size)
+        for x in state[0]:
+            self.append(x)
+}
 };
 /***************************************************************************
  *             GCOMObservation.i - COMPTEL observation class               *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2012-2017 by Juergen Knoedlseder                         *
+ *  copyright (C) 2012-2018 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -38,10 +38,14 @@ class GCOMObservation : public GObservation {
 public:
     // Constructors and destructors
     GCOMObservation(void);
+    explicit GCOMObservation(const GXmlElement& xml);
     GCOMObservation(const GFilename& drename,
                     const GFilename& drbname,
                     const GFilename& drgname,
                     const GFilename& drxname);
+    GCOMObservation(const GFilename&              evpname,
+                    const GFilename&              timname,
+                    const std::vector<GFilename>& oadnames);
     GCOMObservation(const GCOMObservation& obs);
     virtual ~GCOMObservation(void);
@@ -92,4 +96,16 @@ public:
     GCOMObservation copy() {
         return (*self);
     }
+%pythoncode {
+    def __getstate__(self):
+        xml = gammalib.GXmlElement()
+        self.write(xml)
+        state = (xml,)
+        return state
+    def __setstate__(self, state):
+        if state[0].elements() > 0:
+            self.__init__(state[0])
+        else:
+            self.__init__()
+}
 };
 /***************************************************************************
  *                  GCOMResponse.i - COMPTEL Response class                *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2012-2017 by Juergen Knoedlseder                         *
+ *  copyright (C) 2012-2018 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -66,6 +66,7 @@ public:
     const std::string& rspname(void) const;
     void               load(const std::string& rspname);
     void               read(const GFitsImage& hdu);
+    void               write(GFitsImageFloat& image) const;
 };
@@ -76,4 +77,14 @@ public:
     GCOMResponse copy() {
         return (*self);
     }
+%pythoncode {
+    def __getstate__(self):
+        image = gammalib.GFitsImageFloat()
+        self.write(image)
+        state = (image,)
+        return state
+    def __setstate__(self, state):
+        self.__init__()
+        self.read(state[0])
+}
 };
 /***************************************************************************
  *               GCOMRoi.i - COMPTEL region of interest class              *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2017 by Juergen Knoedlseder                              *
+ *  copyright (C) 2017-2018 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -69,4 +69,11 @@ public:
     GCOMRoi copy() {
         return (*self);
     }
+%pythoncode {
+    def __getstate__(self):
+        state = (self.centre(), self.radius(), self.phibar_min(), self.phibar_max())
+        return state
+    def __setstate__(self, state):
+        self.__init__(state[0], state[1], state[2], state[3])
+}
 };
 /***************************************************************************
  *              GCOMTim.i - COMPTEL Good Time Intervals class              *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2017 by Juergen Knodlseder                               *
+ *  copyright (C) 2017-2018 by Juergen Knodlseder                          *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -39,6 +39,7 @@ class GCOMTim : public GBase {
 public:
     // Constructors and destructors
     GCOMTim(void);
+    explicit GCOMTim(const GGti& gti);
     GCOMTim(const GCOMTim& tim);
     GCOMTim(const GFilename& filename, const std::string& usage = "",
                                        const std::string& mode  = "");
@@ -67,4 +68,11 @@ public:
     GCOMTim copy() {
         return (*self);
     }
+%pythoncode {
+    def __getstate__(self):
+        state = (self.gti(),)
+        return state
+    def __setstate__(self, state):
+        self.__init__(state[0])
+}
 };
 /***************************************************************************
  *          GCOMD1Response.cpp - COMPTEL D1 module response class          *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2017 by Juergen Knoedlseder                              *
+ *  copyright (C) 2017-2018 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -30,6 +30,10 @@
 #endif
 #include "GMath.hpp"
 #include "GCOMD1Response.hpp"
+#include "GFitsTable.hpp"
+#include "GFitsBinTable.hpp"
+#include "GFitsTableFloatCol.hpp"
+#include "GFitsTableDoubleCol.hpp"
 /* __ Method name definitions ____________________________________________ */
@@ -359,6 +363,63 @@ void GCOMD1Response::read(const GFitsTable&amp; table)
 /***********************************************************************//**
+ * @brief Write COMPTEL D1 module response.
+ *
+ * @param[in] table FITS table.
+ *
+ * Write the COMPTEL D1 module response into a SDA FITS table.
+ ***************************************************************************/
+void GCOMD1Response::write(GFitsBinTable& table)
+{
+    // Set extension name
+    table.extname("SDA");
+
+    // Set number of entries
+    int num = m_energies.size();
+
+    // Allocate columns
+    GFitsTableFloatCol  col_energy("ENERGY", num);
+    GFitsTableDoubleCol col_position("POSITION", num);
+    GFitsTableDoubleCol col_width("WIDTH", num);
+    GFitsTableDoubleCol col_amplitude("AMPLITUDE", num);
+    GFitsTableDoubleCol col_emin("EMIN", num);
+    GFitsTableDoubleCol col_ewidth("EWIDTH", num);
+    GFitsTableDoubleCol col_emax("EMAX", num);
+
+    // Set column units
+    col_energy.unit("MeV");
+    col_position.unit("MeV");
+    col_width.unit("MeV");
+    col_emin.unit("MeV");
+    col_ewidth.unit("MeV");
+    col_emax.unit("MeV");
+
+    // Copy data from vectors into table
+    for (int i = 0; i < num; ++i) {
+        col_energy(i)    = m_energies[i];
+        col_position(i)  = m_positions[i];
+        col_width(i)     = m_sigmas[i];
+        col_amplitude(i) = m_amplitudes[i];
+        col_emin(i)      = m_emins[i];
+        col_ewidth(i)    = m_ewidths[i];
+        col_emax(i)      = m_emaxs[i];
+    }
+
+    // Append columns to table
+    table.append(col_energy);
+    table.append(col_position);
+    table.append(col_width);
+    table.append(col_amplitude);
+    table.append(col_emin);
+    table.append(col_ewidth);
+    table.append(col_emax);
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
  * @brief Print COMPTEL D1 module response information
  *
  * @param[in] chatter Chattiness.
 /***************************************************************************
  *          GCOMD2Response.cpp - COMPTEL D2 module response class          *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2017 by Juergen Knoedlseder                              *
+ *  copyright (C) 2017-2018 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -31,6 +31,10 @@
 #include "GMath.hpp"
 #include "GIntegral.hpp"
 #include "GCOMD2Response.hpp"
+#include "GFitsTable.hpp"
+#include "GFitsBinTable.hpp"
+#include "GFitsTableFloatCol.hpp"
+#include "GFitsTableDoubleCol.hpp"
 /* __ Method name definitions ____________________________________________ */
@@ -352,6 +356,75 @@ void GCOMD2Response::read(const GFitsTable&amp; table)
 /***********************************************************************//**
+ * @brief Write COMPTEL D2 module response.
+ *
+ * @param[in] table FITS table.
+ *
+ * Write the COMPTEL D2 module response into a SDB FITS table.
+ ***************************************************************************/
+void GCOMD2Response::write(GFitsBinTable& table)
+{
+    // Set extension name
+    table.extname("SDB");
+
+    // Set number of entries
+    int num = m_energies.size();
+
+    // Allocate columns
+    GFitsTableFloatCol  col_energy("ENERGY", num);
+    GFitsTableDoubleCol col_position("POSITION", num);
+    GFitsTableDoubleCol col_width("WIDTH", num);
+    GFitsTableDoubleCol col_amplitude("AMPLITUDE", num);
+    GFitsTableDoubleCol col_escape1("ESCAPE1", num);
+    GFitsTableDoubleCol col_escape2("ESCAPE2", num);
+    GFitsTableDoubleCol col_tail("TAIL", num);
+    GFitsTableDoubleCol col_background("BACKGROUND", num);
+    GFitsTableDoubleCol col_emin("EMIN", num);
+    GFitsTableDoubleCol col_ewidth("EWIDTH", num);
+    GFitsTableDoubleCol col_emax("EMAX", num);
+
+    // Set column units
+    col_energy.unit("MeV");
+    col_position.unit("MeV");
+    col_width.unit("MeV");
+    col_emin.unit("MeV");
+    col_ewidth.unit("MeV");
+    col_emax.unit("MeV");
+
+    // Copy data from vectors into table
+    for (int i = 0; i < num; ++i) {
+        col_energy(i)     = m_energies[i];
+        col_position(i)   = m_positions[i];
+        col_width(i)      = m_sigmas[i];
+        col_amplitude(i)  = m_amplitudes[i];
+        col_escape1(i)    = m_escapes1[i];
+        col_escape2(i)    = m_escapes2[i];
+        col_tail(i)       = m_tails[i];
+        col_background(i) = m_backgrounds[i];
+        col_emin(i)       = m_emins[i];
+        col_ewidth(i)     = m_ewidths[i];
+        col_emax(i)       = m_emaxs[i];
+    }
+
+    // Append columns to table
+    table.append(col_energy);
+    table.append(col_position);
+    table.append(col_width);
+    table.append(col_amplitude);
+    table.append(col_escape1);
+    table.append(col_escape2);
+    table.append(col_tail);
+    table.append(col_background);
+    table.append(col_emin);
+    table.append(col_ewidth);
+    table.append(col_emax);
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
  * @brief Print COMPTEL D2 module response information
  *
  * @param[in] chatter Chattiness.
@@ -861,7 +861,9 @@ void GCOMDri::write(GFits&amp; fits, const std::string&amp; extname) const
     GFitsHDU *image = m_dri.write(fits, extname);
     // Write DRI attributes
-    write_attributes(image);
+    if (image != NULL) {
+        write_attributes(image);
+    }
     // Return
     return;
@@ -1124,7 +1126,9 @@ void GCOMDri::read_attributes(const GFitsHDU* hdu)
     GTime tstop  = com_time(hdu->integer("VIEDAY"), hdu->integer("VIETIM"));
     // Set Good Time Intervals
-    m_gti = GGti(tstart, tstop);
+    if (tstop > tstart) {
+        m_gti = GGti(tstart, tstop);
+    }
     // Optionally read Phibar attributes
     if (hdu->has_card("CDELT3") &&
 /***************************************************************************
  *                GCOMEventBin.cpp - COMPTEL event bin class               *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2012-2016 by Juergen Knoedlseder                         *
+ *  copyright (C) 2012-2018 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -34,6 +34,20 @@
 #include "GTools.hpp"
 /* __ Method name definitions ____________________________________________ */
+#define G_DIR_GET                                       "GCOMEventBin::dir()"
+#define G_ENERGY_GET                                 "GCOMEventBin::energy()"
+#define G_TIME_GET                                     "GCOMEventBin::time()"
+#define G_COUNTS_GET                                 "GCOMEventBin::counts()"
+#define G_SOLIDANGLE_GET                         "GCOMEventBin::solidangle()"
+#define G_EWIDTH_GET                                 "GCOMEventBin::ewidth()"
+#define G_ONTIME_GET                                 "GCOMEventBin::ontime()"
+#define G_DIR_SET                           "GCOMEventBin::dir(GCOMInstDir&)"
+#define G_ENERGY_SET                         "GCOMEventBin::energy(GEnergy&)"
+#define G_TIME_SET                               "GCOMEventBin::time(GTime&)"
+#define G_COUNTS_SET                          "GCOMEventBin::counts(double&)"
+#define G_SOLIDANGLE_SET                  "GCOMEventBin::solidangle(double&)"
+#define G_EWIDTH_SET                         "GCOMEventBin::ewidth(GEnergy&)"
+#define G_ONTIME_SET                          "GCOMEventBin::ontime(double&)"
 /* __ Macros _____________________________________________________________ */
@@ -190,6 +204,98 @@ double GCOMEventBin::size(void) const
 /***********************************************************************//**
+ * @brief Return instrument direction of event bin
+ *
+ * @return Instrument direction of event bin
+ *
+ * @exception GException::invalid_value
+ *            Invalid instrument direction pointer.
+ *
+ * Returns reference to the instrument direction of the event bin.
+ ***************************************************************************/
+const GCOMInstDir& GCOMEventBin::dir(void) const
+{
+    // Throw an exception if instrument direction pointer is not valid
+    if (m_dir == NULL) {
+        std::string msg = "No valid instrument direction found in event bin";
+        throw GException::invalid_value(G_DIR_GET, msg);
+    }
+
+    // Return instrument direction
+    return *m_dir;
+}
+
+
+/***********************************************************************//**
+ * @brief Return energy of event bin
+ *
+ * @return Energy of event bin
+ *
+ * @exception GException::invalid_value
+ *            Invalid energy pointer.
+ *
+ * Returns reference to the energy of the event bin.
+ ***************************************************************************/
+const GEnergy& GCOMEventBin::energy(void) const
+{
+    // Throw an exception if energy pointer is not valid
+    if (m_energy == NULL) {
+        std::string msg = "No valid energy found in event bin";
+        throw GException::invalid_value(G_ENERGY_GET, msg);
+    }
+
+    // Return energy
+    return *m_energy;
+}
+
+
+/***********************************************************************//**
+ * @brief Return time of event bin
+ *
+ * @return Time of event bin
+ *
+ * @exception GException::invalid_value
+ *            Invalid time pointer.
+ *
+ * Returns reference to the time of the event bin.
+ ***************************************************************************/
+const GTime& GCOMEventBin::time(void) const
+{
+    // Throw an exception if time pointer is not valid
+    if (m_energy == NULL) {
+        std::string msg = "No valid time found in event bin";
+        throw GException::invalid_value(G_TIME_GET, msg);
+    }
+
+    // Return time
+    return *m_time;
+}
+
+
+/***********************************************************************//**
+ * @brief Return number of counts in event bin
+ *
+ * @return Number of counts in event bin
+ *
+ * @exception GCTAException::invalid_value
+ *            Invalid counts pointer.
+ *
+ * Returns reference to the number of counts in the event bin.
+ ***************************************************************************/
+double GCOMEventBin::counts(void) const
+{
+    // Throw an exception if counts pointer is not valid
+    if (m_counts == NULL) {
+        std::string msg = "No valid counts found in event bin";
+        throw GException::invalid_value(G_COUNTS_GET, msg);
+    }
+
+    // Return counts
+    return *m_counts;
+}
+
+
+/***********************************************************************//**
  * @brief Return error in number of counts
  *
  * @return Error in number of counts in event bin.
@@ -213,6 +319,257 @@ double GCOMEventBin::error(void) const
 /***********************************************************************//**
+ * @brief Return solid angle of event bin
+ *
+ * @return Solid angle of event bin
+ *
+ * @exception GException::invalid_value
+ *            Invalid solid angle pointer.
+ *
+ * Returns reference to the solid angle of the event bin.
+ ***************************************************************************/
+const double& GCOMEventBin::solidangle(void) const
+{
+    // Throw an exception if counts pointer is not valid
+    if (m_solidangle == NULL) {
+        std::string msg = "No valid solid angle found in event bin";
+        throw GException::invalid_value(G_SOLIDANGLE_GET, msg);
+    }
+
+    // Return solid angle
+    return *m_solidangle;
+}
+
+
+/***********************************************************************//**
+ * @brief Return energy width of event bin
+ *
+ * @return Energy width of event bin
+ *
+ * @exception GException::invalid_value
+ *            Invalid energy width pointer.
+ *
+ * Returns reference to the energy width of the event bin.
+ ***************************************************************************/
+const GEnergy& GCOMEventBin::ewidth(void) const
+{
+    // Throw an exception if energy width pointer is not valid
+    if (m_ewidth == NULL) {
+        std::string msg = "No valid energy width found in event bin";
+        throw GException::invalid_value(G_EWIDTH_GET, msg);
+    }
+
+    // Return energy width
+    return *m_ewidth;
+}
+
+
+/***********************************************************************//**
+ * @brief Return ontime of event bin
+ *
+ * @return Ontime of event bin
+ *
+ * @exception GException::invalid_value
+ *            Invalid ontime pointer.
+ *
+ * Returns reference to the ontime of the event bin.
+ ***************************************************************************/
+const double& GCOMEventBin::ontime(void) const
+{
+    // Throw an exception if energy width pointer is not valid
+    if (m_ontime == NULL) {
+        std::string msg = "No valid ontime found in event bin";
+        throw GException::invalid_value(G_ONTIME_GET, msg);
+    }
+
+    // Return ontime
+    return *m_ontime;
+}
+
+
+/***********************************************************************//**
+ * @brief Set instrument direction of event bin
+ *
+ * @param[in] dir Instrument direction of event bin
+ *
+ * @exception GException::invalid_value
+ *            No memory available to hold instrument direction.
+ *
+ * Sets the instrument direction of the event bin.
+ ***************************************************************************/
+void GCOMEventBin::dir(const GCOMInstDir& dir)
+{
+    // Throw an exception if no memory has been allocated
+    if (m_dir == NULL) {
+        std::string msg = "No memory available to hold instrument direction.";
+        throw GException::invalid_value(G_DIR_SET, msg);
+    }
+
+    // Set instrument direction
+    *m_dir = dir;
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Set energy of event bin
+ *
+ * @param[in] energy Energy of event bin
+ *
+ * @exception GException::invalid_value
+ *            No memory available to hold energy.
+ *
+ * Sets the energy of the event bin.
+ ***************************************************************************/
+void GCOMEventBin::energy(const GEnergy& energy)
+{
+    // Throw an exception if no memory has been allocated
+    if (m_energy == NULL) {
+        std::string msg = "No memory available to hold energy.";
+        throw GException::invalid_value(G_ENERGY_SET, msg);
+    }
+
+    // Set energy
+    *m_energy = energy;
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Set time of event bin
+ *
+ * @param[in] time Time of event bin
+ *
+ * @exception GException::invalid_value
+ *            No memory available to hold instrument direction.
+ *
+ * Sets the time of the event bin.
+ ***************************************************************************/
+void GCOMEventBin::time(const GTime& time)
+{
+    // Throw an exception if no memory has been allocated
+    if (m_time == NULL) {
+        std::string msg = "No memory available to hold time.";
+        throw GException::invalid_value(G_TIME_SET, msg);
+    }
+
+    // Set time
+    *m_time = time;
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Set number of counts in event bin
+ *
+ * @param[in] counts Number of counts.
+ *
+ * @exception GException::invalid_value
+ *            No memory available to hold counts.
+ *
+ * Set the number of counts in the event bin.
+ ***************************************************************************/
+void GCOMEventBin::counts(const double& counts)
+{
+    // Throw an exception if counts pointer is not valid
+    if (m_counts == NULL) {
+        std::string msg = "No memory available to hold counts.";
+        throw GException::invalid_value(G_COUNTS_SET, msg);
+    }
+
+    // Set number of counts in event bin
+    *m_counts = counts;
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Set solid angle of event bin
+ *
+ * @param[in] solidangle Solid angle of event bin
+ *
+ * @exception GException::invalid_value
+ *            No memory available to hold solid angle.
+ *
+ * Sets the solid angle of the event bin.
+ ***************************************************************************/
+void GCOMEventBin::solidangle(const double& solidangle)
+{
+    // Throw an exception if no memory has been allocated
+    if (m_solidangle == NULL) {
+        std::string msg = "No memory available to hold solid angle.";
+        throw GException::invalid_value(G_SOLIDANGLE_SET, msg);
+    }
+
+    // Set solid angle
+    *m_solidangle = solidangle;
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Set energy width of event bin
+ *
+ * @param[in] ewidth Energy width of event bin
+ *
+ * @exception GException::invalid_value
+ *            No memory available to hold energy width.
+ *
+ * Sets the energy width of the event bin.
+ ***************************************************************************/
+void GCOMEventBin::ewidth(const GEnergy& ewidth)
+{
+    // Throw an exception if no memory has been allocated
+    if (m_ewidth == NULL) {
+        std::string msg = "No memory available to hold energy width.";
+        throw GException::invalid_value(G_EWIDTH_SET, msg);
+    }
+
+    // Set energy width
+    *m_ewidth = ewidth;
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Set ontime of event bin
+ *
+ * @param[in] ontime Ontime of event bin (sec).
+ *
+ * @exception GException::invalid_value
+ *            No memory available to hold ontime.
+ *
+ * Sets the ontime of the event bin.
+ ***************************************************************************/
+void GCOMEventBin::ontime(const double& ontime)
+{
+    // Throw an exception if no memory has been allocated
+    if (m_ontime == NULL) {
+        std::string msg = "No memory available to hold ontime.";
+        throw GException::invalid_value(G_ONTIME_SET, msg);
+    }
+
+    // Set solid angle
+    *m_ontime = ontime;
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
  * @brief Print event information
  *
  * @param[in] chatter Chattiness (defaults to NORMAL).
@@ -333,7 +333,7 @@ void GCOMEventList::read(const GFits&amp; file)
  *
  * @param[in] file FITS file.
  *
- * Write the LCOMPTELAT event list into FITS file.
+ * Write the COMPTEL event list into FITS file.
  *
  * @todo Implement method.
  ***************************************************************************/
@@ -77,7 +77,7 @@ GCOMModelDRBFitting::GCOMModelDRBFitting(void) : GModelData()
 /***********************************************************************//**
- * @brief Constructor
+ * @brief XML constructor
  *
  * @param[in] xml XML element.
  *
 /***************************************************************************
  *             GCOMObservation.cpp - COMPTEL Observation class             *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2012-2017 by Juergen Knoedlseder                         *
+ *  copyright (C) 2012-2018 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -85,6 +85,27 @@ GCOMObservation::GCOMObservation(void) : GObservation()
 /***********************************************************************//**
+ * @brief XML constructor
+ *
+ * @param[in] xml XML element.
+ *
+ * Constructs a COMPTEL observation from the information that is found in an
+ * XML element.
+ ***************************************************************************/
+GCOMObservation::GCOMObservation(const GXmlElement& xml) : GObservation()
+{
+    // Initialise members
+    init_members();
+
+    // Read XML
+    read(xml);
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
  * @brief Binned observation constructor
  *
  * @param[in] drename Event cube name.
@@ -487,7 +508,7 @@ void GCOMObservation::write(GXmlElement&amp; xml) const
         par->attribute("file", gammalib::xml_file_reduce(xml, m_drxname));
         // Set IAQ parameter
-        par = gammalib::xml_need_par(G_WRITE, xml, "DRX");
+        par = gammalib::xml_need_par(G_WRITE, xml, "IAQ");
         par->attribute("value", m_response.rspname());
     } // endif: observation was binned
 /***************************************************************************
  *                GCOMResponse.cpp - COMPTEL Response class                *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2012-2017 by Juergen Knoedlseder                         *
+ *  copyright (C) 2012-2018 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -41,6 +41,7 @@
 #include "GTime.hpp"
 #include "GObservation.hpp"
 #include "GFitsImage.hpp"
+#include "GFitsImageFloat.hpp"
 #include "GModelSky.hpp"
 #include "GModelSpatialPointSource.hpp"
 #include "GCOMResponse.hpp"
@@ -592,47 +593,107 @@ void GCOMResponse::load(const std::string&amp; rspname)
  ***************************************************************************/
 void GCOMResponse::read(const GFitsImage& image)
 {
-    // Store IAQ dimensions
-    m_phigeo_bins = image.naxes(0);
-    m_phibar_bins = image.naxes(1);
-
-    // Store IAQ axes definitions
-    m_phigeo_ref_value = image.real("CRVAL1");
-    m_phigeo_ref_pixel = image.real("CRPIX1");
-    m_phigeo_bin_size  = image.real("CDELT1");
-    m_phibar_ref_value = image.real("CRVAL2");
-    m_phibar_ref_pixel = image.real("CRPIX2");
-    m_phibar_bin_size  = image.real("CDELT2");
-
-    // Get axes minima (values of first bin)
-    m_phigeo_min = m_phigeo_ref_value + (1.0-m_phigeo_ref_pixel) * m_phigeo_bin_size;
-    m_phibar_min = m_phibar_ref_value + (1.0-m_phibar_ref_pixel) * m_phibar_bin_size;
-
-    // Compute IAQ size. Continue only if size is positive
-    int size = m_phigeo_bins * m_phibar_bins;
-    if (size > 0) {
-
-        // Allocate memory for IAQ
-        m_iaq.assign(size, 0.0);
-
-        // Copy over IAQ values
-        for (int i = 0; i < size; ++i) {
-            m_iaq[i] = image.pixel(i);
+    // Continue only if there are two image axis
+    if (image.naxis() == 2) {
+
+        // Store IAQ dimensions
+        m_phigeo_bins = image.naxes(0);
+        m_phibar_bins = image.naxes(1);
+
+        // Store IAQ axes definitions
+        m_phigeo_ref_value = image.real("CRVAL1");
+        m_phigeo_ref_pixel = image.real("CRPIX1");
+        m_phigeo_bin_size  = image.real("CDELT1");
+        m_phibar_ref_value = image.real("CRVAL2");
+        m_phibar_ref_pixel = image.real("CRPIX2");
+        m_phibar_bin_size  = image.real("CDELT2");
+
+        // Get axes minima (values of first bin)
+        m_phigeo_min = m_phigeo_ref_value + (1.0-m_phigeo_ref_pixel) * m_phigeo_bin_size;
+        m_phibar_min = m_phibar_ref_value + (1.0-m_phibar_ref_pixel) * m_phibar_bin_size;
+
+        // Compute IAQ size. Continue only if size is positive
+        int size = m_phigeo_bins * m_phibar_bins;
+        if (size > 0) {
+
+            // Allocate memory for IAQ
+            m_iaq.assign(size, 0.0);
+
+            // Copy over IAQ values
+            for (int i = 0; i < size; ++i) {
+                m_iaq[i] = image.pixel(i);
+            }
+
+        } // endif: size was positive
+
+        // Convert IAQ matrix from probability per Phigeo bin into a
+        // probability per steradian
+        double omega0 = gammalib::fourpi *
+                        std::sin(0.5 * m_phigeo_bin_size * gammalib::deg2rad);
+        for (int iphigeo = 0; iphigeo < m_phigeo_bins; ++iphigeo) {
+            double phigeo = iphigeo * m_phigeo_bin_size + m_phigeo_min;
+            double omega  = omega0 * std::sin(phigeo * gammalib::deg2rad);
+            for (int iphibar = 0; iphibar < m_phibar_bins; ++iphibar) {
+                m_iaq[iphigeo+iphibar*m_phigeo_bins] /= omega;
+            }
         }
-    } // endif: size was positive
-
-    // Convert IAQ matrix from probability per Phigeo bin into a
-    // probability per steradian
-    double omega0 = gammalib::fourpi *
-                    std::sin(0.5 * m_phigeo_bin_size * gammalib::deg2rad);
-    for (int iphigeo = 0; iphigeo < m_phigeo_bins; ++iphigeo) {
-        double phigeo = iphigeo * m_phigeo_bin_size + m_phigeo_min;
-        double omega  = omega0 * std::sin(phigeo * gammalib::deg2rad);
-        for (int iphibar = 0; iphibar < m_phibar_bins; ++iphibar) {
-            m_iaq[iphigeo+iphibar*m_phigeo_bins] /= omega;
+    } // endif: image had two axes
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Write COMPTEL response into FITS image.
+ *
+ * @param[in] image FITS image.
+ *
+ * Writes the COMPTEL response into an IAQ FITS file.
+ ***************************************************************************/
+void GCOMResponse::write(GFitsImageFloat& image) const
+{
+    // Continue only if response is not empty
+    if (m_phigeo_bins > 0 && m_phibar_bins > 0) {
+
+        // Initialise image
+        image = GFitsImageFloat(m_phigeo_bins, m_phibar_bins);
+
+        // Convert IAQ matrix from probability per steradian into
+        //probability per Phigeo bin
+        double omega0 = gammalib::fourpi *
+                        std::sin(0.5 * m_phigeo_bin_size * gammalib::deg2rad);
+        for (int iphigeo = 0; iphigeo < m_phigeo_bins; ++iphigeo) {
+            double phigeo = iphigeo * m_phigeo_bin_size + m_phigeo_min;
+            double omega  = omega0 * std::sin(phigeo * gammalib::deg2rad);
+            for (int iphibar = 0; iphibar < m_phibar_bins; ++iphibar) {
+                int i          = iphigeo + iphibar*m_phigeo_bins;
+                image(iphigeo, iphibar) = m_iaq[i] * omega;
+            }
         }
-    }
+
+        // Set header keywords
+        image.card("CTYPE1", "Phigeo", "Geometrical scatter angle");
+        image.card("CRVAL1", m_phigeo_ref_value,
+                   "[deg] Geometrical scatter angle for reference bin");
+        image.card("CDELT1", m_phigeo_bin_size,
+                   "[deg] Geometrical scatter angle bin size");
+        image.card("CRPIX1", m_phigeo_ref_pixel,
+                   "Reference bin of geometrical scatter angle");
+        image.card("CTYPE2", "Phibar", "Compton scatter angle");
+        image.card("CRVAL2", m_phibar_ref_value,
+                   "[deg] Compton scatter angle for reference bin");
+        image.card("CDELT2", m_phibar_bin_size, "[deg] Compton scatter angle bin size");
+        image.card("CRPIX2", m_phibar_ref_pixel,
+                   "Reference bin of Compton scatter angle");
+        image.card("BUNIT", "Probability per bin", "Unit of bins");
+        image.card("TELESCOP", "GRO", "Name of telescope");
+        image.card("INSTRUME", "COMPTEL", "Name of instrument");
+        image.card("ORIGIN", "GammaLib", "Origin of FITS file");
+        image.card("OBSERVER", "Unknown", "Observer that created FITS file");
+
+    } // endif: response was not empty
     // Return
     return;
 /***************************************************************************
  *             GCOMTim.cpp - COMPTEL Good Time Intervals class             *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2017 by Juergen Knodlseder                               *
+ *  copyright (C) 2017-2018 by Juergen Knodlseder                          *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -63,6 +63,24 @@ GCOMTim::GCOMTim(void)
 /***********************************************************************//**
+ * @brief Good Time Interval constructor
+ *
+ * @param[in] gti Good Time Intervals.
+ ***************************************************************************/
+GCOMTim::GCOMTim(const GGti& gti)
+{
+    // Initialise class members
+    init_members();
+
+    // Set GTIs
+    m_gti = gti;
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
  * @brief Copy constructor
  *
  * @param[in] tim COMPTEL Good Time Intervals.
 # ==========================================================================
 # This module performs unit tests for the GammaLib COMPTEL module.
 #
-# Copyright (C) 2012-2015 Juergen Knoedlseder
+# Copyright (C) 2012-2018 Juergen Knoedlseder
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -17,7 +17,9 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
 # ==========================================================================
+import os
 import gammalib
+import test_support
 # ====================================== #
@@ -25,12 +27,12 @@ import gammalib
 # ====================================== #
 class Test(gammalib.GPythonTestSuite):
     """
-    Test class for GammaLib COMPTEL module.
+    Test class for GammaLib COMPTEL module
     """
     # Constructor
     def __init__(self):
         """
-        Constructor.
+        Constructor
         """
         # Call base class constructor
         gammalib.GPythonTestSuite.__init__(self)
@@ -38,24 +40,127 @@ class Test(gammalib.GPythonTestSuite):
         # Return
         return
-    # Set test functions
-    def set(self):
+    # Test class pickeling
+    def _test_pickeling(self):
         """
-        Set all test functions.
+        Test class pickeling
         """
-        # Set test name
-        self.name("COMPTEL")
+        # Set CALDB
+        os.environ['CALDB'] = os.environ['TEST_COM_DATA']+'/../caldb'
-        # Append tests
-        self.append(self.test, "COMPTEL dummy test")
+        # Perform pickeling tests of empty classes
+        test_support._pickeling(self, gammalib.GCOMD1Response())
+        test_support._pickeling(self, gammalib.GCOMD2Response())
+        test_support._pickeling(self, gammalib.GCOMDri())
+        test_support._pickeling(self, gammalib.GCOMEventAtom())
+        test_support._pickeling(self, gammalib.GCOMEventBin())
+        test_support._pickeling(self, gammalib.GCOMEventCube())
+        test_support._pickeling(self, gammalib.GCOMEventList())
+        #test_support._pickeling(self, gammalib.GCOMIaq())
+        #test_support._pickeling(self, gammalib.GCOMInstChars())
+        test_support._pickeling(self, gammalib.GCOMInstDir())
+        test_support._pickeling(self, gammalib.GCOMModelDRBFitting())
+        test_support._pickeling(self, gammalib.GCOMOad())
+        test_support._pickeling(self, gammalib.GCOMOads())
+        test_support._pickeling(self, gammalib.GCOMObservation())
+        test_support._pickeling(self, gammalib.GCOMResponse())
+        test_support._pickeling(self, gammalib.GCOMRoi())
+        #test_support._pickeling(self, gammalib.GCOMSelection())
+        #test_support._pickeling(self, gammalib.GCOMStatus())
+        test_support._pickeling(self, gammalib.GCOMTim())
+
+        # Setup test
+        caldb   = gammalib.GCaldb('cgro','comptel')
+        dir     = gammalib.GSkyDir()
+        instdir = gammalib.GCOMInstDir(dir, 10.0)
+        energy  = gammalib.GEnergy(1.0,'MeV')
+        time    = gammalib.GTime(1.0,'sec')
+        map     = gammalib.GSkyMap('CAR','GAL',0.0,0.0,1.0,1.0,10,10,10)
+        dri     = gammalib.GCOMDri(map, 0.0, 2.0, 10)
+        atom    = gammalib.GCOMEventAtom()
+        ebds    = gammalib.GEbounds(1,gammalib.GEnergy(0.75,'MeV'),gammalib.GEnergy(30.0,'MeV'))
+        gti     = gammalib.GGti(gammalib.GTime(1.0,'secs'),gammalib.GTime(3.0,'secs'))
+        dre     = gammalib.GCOMDri(map, 0.0, 2.0, 10)
+        obs     = gammalib.GObservations(os.environ['TEST_COM_DATA']+'/obs.xml')
+        model   = gammalib.GModels(os.environ['TEST_COM_DATA']+'/crab.xml')
+        dre.ebounds(ebds)
+        dre.gti(gti)
+        atom.dir(instdir)
+        atom.energy(energy)
+        atom.time(time)
+        atom.e1(100.0)
+        atom.e2(200.0)
+        atom.phibar(11.1)
+        atom.theta(22.2)
+        atom.phi(33.3)
+        atom.eha(44.4)
+        atom.psd(55.5)
+        atom.tof(66.6)
+        atom.modcom(1)
+        atom.reflag(2)
+        atom.veto(3)
+        list = gammalib.GCOMEventList()
+        list.append(atom)
+        bin = gammalib.GCOMEventBin()
+        bin.dir(instdir)
+        bin.energy(energy)
+        bin.time(time)
+        bin.counts(11.0)
+        bin.solidangle(11.0)
+        bin.ewidth(energy)
+        bin.ontime(111.0)
+        oad = gammalib.GCOMOad()
+        oad.tstart(gammalib.GTime(1.0,'secs'))
+        oad.tstop(gammalib.GTime(3.0,'secs'))
+        oad.zaxis(dir)
+        oad.xaxis(dir)
+        oad.tjd(1000)
+        oad.tics(100)
+        oad.gcaz(10.0)
+        oad.gcel(20.0)
+        oad.georad(30.0)
+        oads = gammalib.GCOMOads()
+        oads.append(oad)
+
+        # Take provision for rounding errors
+        fits = gammalib.GFits()
+        dre.write(fits)
+        dre.read(fits[0])
+
+        # Perform pickeling tests of filled classes
+        test_support._pickeling(self, gammalib.GCOMD1Response(caldb,'DEFAULT'))
+        test_support._pickeling(self, gammalib.GCOMD2Response(caldb,'DEFAULT'))
+        test_support._pickeling(self, gammalib.GCOMDri(dri))
+        test_support._pickeling(self, gammalib.GCOMEventAtom(atom))
+        test_support._pickeling(self, gammalib.GCOMEventBin(bin))
+        test_support._pickeling(self, gammalib.GCOMEventCube(dre))
+        test_support._pickeling(self, gammalib.GCOMEventList(list))
+        #test_support._pickeling(self, gammalib.GCOMIaq())
+        #test_support._pickeling(self, gammalib.GCOMInstChars(caldb,'DEFAULT'))
+        test_support._pickeling(self, gammalib.GCOMInstDir(instdir))
+        test_support._pickeling(self, gammalib.GCOMModelDRBFitting(model['Background']))
+        test_support._pickeling(self, gammalib.GCOMOad(oad))
+        test_support._pickeling(self, gammalib.GCOMOads(oads))
+        test_support._pickeling(self, gammalib.GCOMObservation(obs[0]))
+        test_support._pickeling(self, gammalib.GCOMResponse(caldb, 'UNH(1.0-3.0)MeV'))
+        test_support._pickeling(self, gammalib.GCOMRoi(instdir, 10.0, 1.0, 5.0))
+        #test_support._pickeling(self, gammalib.GCOMSelection())
+        #test_support._pickeling(self, gammalib.GCOMStatus())
+        test_support._pickeling(self, gammalib.GCOMTim(gti))
         # Return
         return
-    # Test function
-    def test(self):
+    # Set test functions
+    def set(self):
         """
-        Test function.
+        Set all test functions
         """
+        # Set test name
+        self.name('COMPTEL')
+
+        # Append tests
+        self.append(self._test_pickeling, 'Test COMPTEL class pickeling')
+
         # Return
         return
@@ -277,7 +277,7 @@ const GTime&amp; GCTAEventBin::time(void) const
  *
  * @return Number of counts in event bin
  *
- * @exception GCTAException::no_member
+ * @exception GCTAException::invalid_value
  *            Invalid counts pointer.
  *
  * Returns reference to the number of counts in the event bin.
@@ -165,6 +165,25 @@ if WITH_INST_CTA
   testdc1cta_DATA   = $(testdc1cta)
 endif
+# Distribute COM test data (optionally)
+if WITH_INST_COM
+  testdatacom       = $(top_srcdir)/inst/com/test/data/obs.xml \
+                      $(top_srcdir)/inst/com/test/data/crab.xml \
+                      $(top_srcdir)/inst/com/test/data/m32171_drx.fits.gz \
+                      $(top_srcdir)/inst/com/test/data/m34997_drg.fits.gz \
+                      $(top_srcdir)/inst/com/test/data/m50439_dre.fits.gz
+  testdatacomdir    = $(pyexecdir)/gammalib/tests/com/data
+  testdatacom_DATA  = $(testdatacom)
+  testcaldbcom      = $(top_srcdir)/inst/com/test/caldb/data/cgro/comptel/caldb.indx
+  testcaldbcomdir   = $(pyexecdir)/gammalib/tests/com/caldb/data/cgro/comptel
+  testcaldbcom_DATA = $(testcaldbcom)
+  testirfcom        = $(top_srcdir)/inst/com/test/caldb/data/cgro/comptel/bcf/m10371_sda.fits.gz \
+                      $(top_srcdir)/inst/com/test/caldb/data/cgro/comptel/bcf/m10564_sdb.fits.gz \
+                      $(top_srcdir)/inst/com/test/caldb/data/cgro/comptel/bcf/u47569_iaq.fits.gz
+  testirfcomdir     = $(pyexecdir)/gammalib/tests/com/caldb/data/cgro/comptel/bcf
+  testirfcom_DATA   = $(testirfcom)
+endif
+
 # If building of Python wrappers was requested then add them to the build
 # source
 if BUILD_WRAPPERS